Abstract
Lassa fever (Lf) is a viral haemorrhagic disease endemic to West Africa and is caused by the Lassa mammarenavirus. The rodent Mastomys natalensis serves as the primary reservoir and its ecology and behaviour have been linked to the distinct spatial and temporal patterns in the incidence of Lf. Nigeria has experienced an unprecedented epidemic that lasted from January until April of 2018, which has been followed by subsequent epidemics of Lf in the same period every year since. While previous research has modelled the case seasonality within Nigeria, this did not capture the seasonal variation in the reproduction of the zoonotic reservoir and its effect on case numbers. To this end, we introduce an approximate Bayesian computation scheme to fit our model to the case data from 2018–2020 supplied by the NCDC. In this study we used a periodically forced seasonal nonautonomous system of ordinary differential equations as a vector model to demonstrate that the population dynamics of the rodent reservoir may be responsible for the spikes in the number of observed cases in humans. The results show that in December through to March, spillover from the zoonotic reservoir drastically increases and spreads the virus to the people of Nigeria. Therefore to effectively combat Lf, attention and efforts should be concentrated during this period.
Author summary Lassa fever is a viral disease prevalent in West Africa, with Mastomys natalensis serving as the primary reservoir. In Nigeria, annual outbreaks occur from December to March. Using a novel model and data from 2018-2020, we demonstrate that the population dynamics of the reservoir contribute to spikes in human cases. Specifically, spillover transmission increases drastically during this period, highlighting the need for concentrated efforts and interventions. Understanding the seasonal dynamics of the reservoir is crucial for effective Lassa fever control and prevention strategies in Nigeria.
Introduction
Lassa fever (Lf) is a viral zoonotic disease, caused by the Lassa mammarenavirus (LASV), that is endemic to West African countries such as Nigeria, Sierra Leone and Guinea [1, 2]. Lf has a natural reservoir in the rodent Mastomys natalensis in which the virus persists and crossover events to humans occur [3, 4, 5]. The disease, which was first described after two nurses contracted the disease in a hospital in Jos, Nigeria, in 1969, has since been identified as a significant risk to health in West Africa with 300, 000–500, 000 cases per year resulting in approximately 5, 000 deaths annually [3, 6]. Those infected with LASF typically experience acute symptoms of headaches, sore throat, muscle pain, vomiting and diarrhoea, and in severe cases bleeding from the mouth, nose, vagina or gastrointestinal tract [1]. The risk that Lf poses to public health will only increase without widespread intervention and a viable vaccine as growth in inter-border traffic and international travel increases the likelihood of introducing the virus to other regions within and outside of the African continent [7, 8].
In recent years, Nigeria has experienced epidemics with peak incidence occurring between December and April in the years since 2017, just after the rainy season ends in October. Both the number of cases and the exposure rate of Lf increases in certain periods of the year and have been correlated with rainfall patterns [9, 10, 11]. This may be because the reproduction of Mastomys is greatest just after the rainy season, which results in an increase in the size of the rat population and in the spread of LASV from infected to susceptible rats[12, 13, 14, 15]. The ecological dynamics of M. natalensis are relevant to Lf in humans because the majority of infections (80%) are suspected to be spillover events as opposed to human-to-human transmissions [16]. The importance of rat contact with humans is illustrated by data reported by Tobin et al. (2015) for Edo state, Nigeria, which recorded 32.4%(385/1189) of the confirmed national cases of Lf in 2020; 96.1% houses had found the multimammate rat within them in the past 6 months and 58.2% of the resident were seropositive (i.e. tested positive for Lf-specific antibodies) [17]. Hence the relationship between the presence and behaviour of M. natalensis and the prevalence of Lf is critical to understanding and predicting future outbreaks.
Mathematical modelling studies of Lf are rare compared with other diseases, despite the inclusion of Lf in the World Health Organization’s Blueprint list of diseases to be prioritized for research and development [18]. Published models for recent outbreaks in Nigeria can incorporate the population dynamics of the disease reservoir and highlight areas and periods of the year at high risk of transmission [19, 20, 21]. Akhmetzhanov et al. (2019) used suspected case data and a rodent model to inform a time-dependent exposure rate of Lf to susceptible people in Nigeria, but did not include explicit modelling of human infection dynamics [9]. To date, the mechanistic models describing Lf epidemics in Nigeria lacked the focus on time-dependent parameters relating to the rodent population dynamics to explain the relationship between those dynamics and the seasonality of outbreaks.
In this study, we developed an epidemiological model to describe the temporal dynamics of Lf within Nigeria in both human hosts and rodent vectors incorporating seasonal variations in rodent population dynamics. The model was then fitted to confirmed case data from the Nigerian Centre for Disease Control (NCDC) at the national scale in Nigeria using Approximate Bayesian Computation (ABC). This approach allowed us to demonstrate how the annual fluctuations of the rodent population translated into the seasonal outbreaks of Lf cases in human hosts.
We believe that providing an epidemiological model of the dynamics of Lf in humans and the rodent reservoir within Nigeria as a whole that fits to confirmed cases supports the hypothesis that the fluctuations of the rodent population strongly influence the cases observed. In order to inform future public health efforts, it is important to better understand the role of the disease reservoir in the spatio-temporal profile of infections.
Materials and methods
The weekly situation reports on Lf produced by the NCDC provided a stream of publicly available data that were discussed in regular direct communications with NCDC representatives. The data used in this study was the set of dated Lf cases collected by the NCDC’s surveillance network for Lf between 7th January 2018 until 12th July 2020, which covers the 2018, 2019 and 2020 epidemics [22]. Lf cases were categorised as either suspected, confirmed and probable. Confirmed cases were those which had a positive result for IgM antibody, PCR or virus isolation. Suspected cases were any individual experiencing symptoms such as fever, sore throat, vomiting, diarrhoea. Additionally they also met one of the following criteria: if they had a history of contact with either 1) excreta or urine of rodents, 2) with a probable or confirmed Lf case recently, or 3) any person with inexplicable bleeding/haemorrhagia. Due to the uncertainty that would result from use of unconfirmed data we used only confirmed cases.
Model
In order to describe the confirmed Lf cases in Nigeria, we constructed a vector-host model in which the dynamics of the rodent M. natalensis are modelled explicitly in addition to the dynamics of human transmission and disease progression (Fig 1 and equations 1–3).
Blue solid arrows denote recruitment. Black solid arrows denote progression of the disease. Red dashed arrows denote disease transmission. Purple solid arrows denote mortalities. Parameters are detailed in full in table 1 where λh and λr are defined in equations 2 (i) and (ii) respectively, and B(t) is defined in equation 1.
The human population was split into susceptible, exposed, asymptomatic or infected, and recovered (Sh, Eh, Ah, Ih, Rh). The infection pathway for humans was described in the model as follows: Susceptible people acquire infections from the pool of all infectious individuals. This is at a rate of λh. Contact with an infectious host, whether it is a human or vector rodent, transmits the disease to the susceptible individual and then they become exposed. Exposed humans become infectious after a period of 1/ν days, and thus exposed individuals will leave the compartment at a rate of ν. Those who were exposed to the virus will then become Asymptomatic with a probability of p and Infected with a probability of 1 − p. Infected humans are assumed to experience more severe symptoms and be recorded as cases in data; they also experience an infection induced mortality rate . Infected individuals are assumed to be detectable. Infected individuals will recover with an average recovery period of 1/γh days.
The parameters of interest were inferred using algorithm 1 in section 2.4.
All humans are assumed to have an average life expectancy of 1/μh days regardless of their compartment. That is, the number of individuals in each compartment decrease at a rate of μh per capita. Humans are assumed to have a constant birth rate, Bh, and newborns are assumed to be fully susceptible to Lf.
The vector population of rodents was split into susceptible, infected, and recovered (Sr, Ir, Rr). Rats are born susceptible at a time-dependent per capita rate of B(t). The rat infections follow a more simple pathway than that found in humans: Susceptible rats come into contact with infected rats, whereby the disease is transmitted to the susceptible individual. This occurs at a rate of λr. The infected rats do not experience increased mortality and thus the rats have a constant mortality rate of μr. Infectious individuals will recover with an average recovery period of 1/γr.
M. natalensis has repeatedly been observed to have seasonal breeding habits over different areas of Africa [12, 13, 14]. Therefore the recruitment rate per capita of the rats B(t) was chosen as
where k is the magnitude of the function; s is a shape parameter denoting how long the period of low reproduction rates lasts for, a smaller s meaning a close to constant recruitment rate over the year whereas a larger s would equate to a long low period then a sharper change to a high period; and ϕ is the point in the year where the reproduction of the rats is at its minimum. k and s are of positive value, and ϕ is between 0 and 1. Once parameters s and the natural mortality rate of the rats, μr, have been fixed, k may be scaled to keep the reservoir population constant, year-on-year [23].
A population whose yearly dynamics are similar is a better representation of a species that is endemic to the environment and although fluctuations will happen, the data to accurately model the population does not exist.
The rate at which susceptibles become infected, otherwise known as the force of infection, is denoted by λh and λr for humans and rats respectively. In humans, this is defined as a linear combination of the possible contact routes with LASV carriers. To represent the difference in transmissibility of Lf in rats and humans, the contact rates βrh and βhh are separated. Explicitly, these parameters are the rate of successful transmissions from rats to humans per infected rat and humans to humans per infected human respectively. In addition, we have assumed a density dependent contact rate between susceptible humans and the pool of infectious individuals, which means that contacts will occur at an invariable rate irrespective of the size of the human population.
With the transmission rates between compartments defined, the susceptible compartments experience a force of infection from the infectious agents. The total force of infection per susceptible is therefore the sum linear combination of the number of infectious agents that can infect that susceptible host multiplied by the appropriate transmission rate. Therefore the force of infection per individual human and rat, respectively denoted λh and λr, are as follows:
Basic reproduction numbers
The basic reproduction ratio, R0, is the expected number of secondary infections caused by a single infectious agent in an otherwise susceptible population. Over time, the conditions which the infectious agent inhabit change and thus we also calculate the effective reproduction rate, which shows the expected number of secondary infections from one infectious agent in the population (which may have other infectious agents) at a specified time t.
Based on the next generation method for deriving the reproductive ratio, R0, from Diekmann et al (1990) and the particular method used in van den Driessche and Watmough’s work (2002) we produce the effective reproductions rates between species in equation set 4 with derivations in S1 appendix (Section 1.2) [24, 25].
Parameter Selection
In table 1 we show the choices for the model parameters that are described in this subsection.
Outbreak events within Nigeria have a seasonal pattern, thus each epidemic should not be analysed in isolation as this would neglect the effect of seasonal dynamics. Therefore the data fitting takes place over a longer time period and the population dynamics of the people of Nigeria should be taken into account. The natural death rate of humans μh is estimated to be day−1 since 54 years was the average life expectancy to 2 significant figures for Nigerians given by the World Bank for 2018. The human birth rate Bh was approximated as 1.2 × 10−4day−1. Nigeria has experienced close to exponential growth rate in recent years and if it is assumed that the population growth of Nigeria has been stable over this time period then we may assume that
, where Bh and μh are constant. Therefore
. We then obtain
. The growth of Nigeria from 2015 to 2019 was 181.1 million to 201.0 million to 4 significant figures which gives the growth rate of Nigeria to be 7.14 × 10−5 per capita per day hence giving Bh = μh + 7.14 × 10−5 [26].
The incubation period was assumed to be 14 days. There is a wide range for the incubation period and is reported to be around 2 days to 3 weeks. For simplification this assumption was made to be approximately 14 days and thus ν = 1/14. The rate of recovery for humans γh was 0.1 day−1. Ranges for the time to recover are broad, between 2 and 21 days, so similarly to ν a value of 0.1 was assumed. We assume that the probability of becoming asymptomatic p = 0.8 [1].
The infection induced mortality rate μhI = 0.195γ = 0.0195 day−1. The proportion of those that have died, retrospectively, during the outbreak period considered is 196/1006. Therefore the case fatality rate (CFR) is 19.5% to 3 significant figures [22]. To give a rate of mortality the CFR is multiplied by the recovery rate to approximate the rate that individuals die from Lf before they recover[27].
The initial number of infected humans Ih(0) is 2 since there were two recorded confirmed cases in the week commencing 01/01/2018. Therefore Ah(0) = 8 to maintain the ratio between asymptomatic and symptomatic infected persons to 1:4. Eh(0) is 5 times that of the number of cases reported the week after the data being used starts.
This is done to represent that those who show symptoms were likely exposed the week before and that 20% of the exposed will go on to show symptoms. The initial total number of people living in Nigeria was assumed to be 2 × 108. The initial number of susceptible Sh(0) is therefore the remaining number of humans, 2 × 108.
The mortality rate of rats is assumed as μr = 1/500 per day since this value has been previously used as the baseline value in previous works [20, 28]. γr = 1/90 per day [29]. For the population size of the zoonotic reservoir we tested four initial values for size, Nr(0) ∈ {1, 106, 5 × 106, 2.5 × 107}. The number of initially susceptible rats was . This is necessarily bounded by Nr(0) above and below by 0 since betarr may be sampled such that
, the total initial rat population, was taken from an array of values to test for the impact on the system,
. The later value was chosen as it is the endemic equilibrium of an SIR system with a recruitment/birth function equal to μr. The remaining rats were assumed to be recovered.
The remaining parameters, ϕ, s, βrr, βrh and βhh will be estimated within the fitting scheme in the following section.
Fitting and Data
Bayesian estimation techniques involve a suite of statistical inference methods based on the idea that after specifying a prior assumption upon the parameters being investigated, the prior is updated with the introduction of more information from the observed data. Following Bayes’ theorem, the posterior distribution of parameters is obtained by combining the prior beliefs (prior distribution) with the evidence of the data which usually comes in the form of the likelihood function [30].
In the absence of a likelihood function, which may arise because the function is intractable or computationally expensive, Approximate Bayesian Computation (ABC) is a robust method that can be used. ABC can be summarised as a family of techniques where parameters are sampled, in various different ways that are dependent on the specific scheme, and then accepted or rejected if the simulated data given the parameters are sufficiently close to the observed data [31]. ABC schemes are becoming increasing popular due to their relative ease of use.
In this paper we use a modified Approximate Bayesian Computation Sequential Monte Carlo scheme (ABC SMC) to fit our parameters (see Algorithm 1 for psuedocode description of scheme used) since the scheme has been shown to be reliable and converges faster than some of the more primitive schemes [32]. The ABC SMC schema iterates a population of parameter particles over T generations with decreasing tolerances, , allowed between the data, y*,and the data simulated from the model with the particle θ, yθ. This converges to the desired approximate posterior distribution as the distributions of the parameters are sequentially improved upon [33]. ABC SMC fits a model M with unknown model parameters θ to data. The standard algorithm requires one to specify a decreasing sequence of thresholds ε1 ≥ ε2 ≥ … ≥ εT for the T generations. When starting, t = 1, parameters are sampled from prior distributions, π(θ). For each subsequent generation t = 2, …, T parameters will be sampled from a perturbation kernel,
, based on a sampled particle accepted in generation t − 1. The model is then simulated with the sampled parameter particle and using a chosen distance metric to compare the data and the simulated data, the parameters are accepted if the error calculated is smaller than the given tolerance for that generation, i.e. d(y∗, yθ) ≤ εt.
When implementing the SMC algorithm, instead of manually defining the sequence ε1, ε2, …, εT, which is often done by manually calibrating the tolerances after some initial test runs, we instead initialise our algorithm with an integer, K, and a proportion, Q. If the desired number of parameter particles is N then we initialise by sampling KN particles from the prior distributions and reject all the N particles with the smallest errors. This serves as our first generation of sampling. For subsequent generations we have set a desired quantile, Q, where the particle that is the Qth quantile has its error between the data and model set as the tolerance for the next generation. That is, if the set of parameter particles for generation g − 1 is and is ordered with respect to d(y*, yθ) then for generation g the tolerance
. Furthermore the perturbation kernel that we use is a multivariate Gaussian with variance equal to twice that of the co-variance between the previous generation’s particles:
,where i is sampled using the weightings generated.
Prior Distributions
For the parameters that we are unable to calculate, we fit with the scheme from section 2.4 specifically described as Algorithm 1. In this section we explain the choices for the prior distributions from our prior knowledge on the model and disease epidemiology. ϕ is sampled from a uniform distribution between 0 and 1 because this does not give any preference to any period of the year over another. s is sampled from a Log-normal distribution, Lognormal(μ, σ2), where μ = 3 and σ= 1 as this allows the parameter samples to be varied and comparable to the range used in Peel et al (2014) [23]. βrr ∼ Lognormal (μ, σ2). Since the zoonotic system is assumed to be in endemic equilibrium, is assumed to be greater than 1. We therefore set μ so that the mode of the distribution, m1 = exp μ − σwhere σ= 1, would equate to
with a prior flexible enough to sample through parameter space for a variety of values. Therefore μ = log(10(γr + μr)) + σ2 = −1.03 to 3 s.f. As our system is not at the point of introduction of Lf, rather it is close to endemic equilibrium, the rate at which new
This was used to fit the model in section 2.1. Instead of using an arbitrary sequence of tolerances, the tolerances are calculated from the errors produced in the previous generation. For the first generation, the algorithm runs a multiple, K, of the number of desired particles, N, and then accepts the best N. infections are acquired by humans from rats will be proportional to Nr and so to factor this in we set the Log-normal prior distribution for βrh to have a mode m2 to be such that where A = 2, 500. Thus the prior distribution has μ = log(10(γr + μr)) + σ2 + log(A/Nr(0)) = −7.77 to 3 s.f. Since human-to-human transmissions are unlikely outside of nosocomial settings, we set σ= 0.5 and βhh ∼ Lognormal(μ, σ2) so that the mode m3 when equated to βhh and inserted into
would give
. Therefore μ = log(ψ) + σ2 = −2.35 to 3 s.f. where
As previously stated, the data used here are taken directly from our communications with the NCDC. The data points used are laboratory confirmed cases aggregated by week from the 12th of December 2018 until the 4th of April 2020.
Implementation
The fitting algorithm and all models were implemented in MATLAB, with the models using the ODE45 solver for numerical integration of the ODEs.
The algorithm ran for T = 15 generations with each generation consisting of 2500 parameter particles. For the first generation K = 10 multiples of 2500 particles were ran and then the 1/K quantile with the least error was accepted. Subsequent tolerances for each generation were set with the quantile Q = 1/6, where errors from the previous generation would generate the next tolerance value.
Results
We applied an ABC fitting scheme (Algorithm 1) to the model detailed in Materials and methods (Fig 1 and eq. 3) to the confirmed Lassa fever cases in Nigeria. Notably, our fitted model successfully replicated the observed seasonal trends in the data, as demonstrated in Figure 2. The posterior distributions of parameters and vector dynamics (Fig 3 and 4, respectively) reveal the seasonal nature of the epidemics and how the population dynamics of the primary reservoir, M. natalensis, affect the number of observed cases in Nigeria.
The simulated cases compared with the observed data. In orange is the entire range of values Ih takes in the final generation at each time point; the median value in blue. Confirmed case data for Nigeria are in black. The model replicates the sharp increase in case incidences occurring at the start of the year for 3 years.
Figure 3 top left the shape parameter of the rodent recruitment function, s. Figure 3 top right the rodent-to-rodent transmission rate βrr. Figure 3 mid left the human-to-human transmission rate βhh. Figure 3 mid right the rodent-to-human transmission rate βrh. Figure 3 bottom left the date of minimum rat reproduction ϕ.
The figure showcases the evolution of M. natalensis compartments throughout the observed period, simulated using the parameters derived from the final generation of accepted values. The median value is represented by the dashed line, while the colored area illustrates the range. Susceptible rats are depicted in red, infected rats in green, and recovered rats in blue. Notably, the recruitment of susceptible rats progressively rises, providing impetus for the growth of infected rats, reaching its peak in late December. Consequently, this surge in infected rats leads to spillover infections in humans.
Model Evaluation
The model fit resulted in a final generation whose simulations of the number of symptom-presenting humans, Ih, can be seen in Fig 2. The entire range of values that Ih takes in the final generation at each time point is in light orange and the median value in blue. Overlaid in black is the confirmed case data for Nigeria. This shows that the model can replicate the year-on-year trend and that the seasonal epidemics in Nigeria can be explained by the vector population dynamics.
The number of infected rats increases drastically in December when the pool of available susceptibles grows, thereby increasing the spillover rate to humans (Fig 4). As the number of rats in contact with humans decreases over the year due to natural mortality and recovery the number of spillover events decreases. This process starts again just before the next epidemic and continues cyclically.
Marginal Posterior Parameter distributions
The recruitment rate of the vector M. natalensis is determined by the shape parameter, s fig 3 top left, and the date of the minimum rate, ϕ fig 3 bottom left, is focused and seasonal. Since the posterior of s has increased its median substantially (6.08 × 102), this results in a large ratio between the reproduction rate’s lowest and highest value and thus the rat population experiences an influx of susceptible rats at ϕ+6 months in early December. The distribution for ϕ is concentrated around a median of the 7th of June with a 90% credible interval of the 4th – 11th of June (see S1 appendix table 2 for unconverted values). This is observed in the rodent dynamics in fig 4 where the number of susceptible rats increases rapidly.
The time-varying reproduction rate for rat-to-rat transmission crosses the threshold of in early December (fig 5 (a)) and causes the number of infectious rats to increase similarly. This then spills over to humans and causes the spike in incidence data that is observed in Nigeria between January and March.
Rrr(t) and when the threshold for Rrr(t) ≥ 1 is met. In Fig 5 (a) Rrr(t), median dashed-line and range in coloured block, exhibits a sharp increase towards the end of the year, foreshadowing the subsequent outbreaks in the following months. To maintain clarity, the data is limited to the years 2019 and 2020, as no complete earlier records are available. Figure 5 (b) showcases a box diagram illustrating the time of year when Rrr(t) exceeds the threshold of 1, denoting high transmission. The bottom panel of Fig 5 (b) captures the onset of the high transmission period, while the upper panel displays its conclusion. The intermediate phase witnesses a rapid shift in reservoir dynamics, leading to an escalation in the number of infected vectors
The transmission parameters, fig 3 top right and middle, are such that human infections are predominantly the result of a spillover event from the zoonotic reservoir. Despite rat-to-human transmission, βrh, low estimation as seen in S1 appendix table 2 with the expected number of secondary human infections per infected rat in a completely susceptible population, , having a median of 4.70 × 10−3, the number of infected rats during December through to May is sufficient to cause a significant number of infections in humans. Moreover, the proportion of humans infected by rats is estimated to be 97.42% – 99.72% (90% credible interval) with a median of 98.51%.
Biological Implications
In our study, we explicitly modelled the vector dynamics allowing us to investigate biological implications. The seasonal recruitment of the natural reservoir as seen in Fig 4 is a crucial aspect of the disease dynamics and serves as a key driver of the transmission cycle. We can infer from these results that M. natalensis recruitment, which may be an combination of birth and migration, influence the occurrence of spillover events to humans in Nigeria, with rat-to-rat infection peaking just before the epidemics observed in Nigeria (Figs 5 (a) and 5 (b)). This emphasizes the importance of understanding the ecological and reproductive dynamics of the reservoir species, as it directly impacts the risk of disease transmission to human populations.
Sensitivity Analysis
We conducted a sensitivity analysis (see S1 appendix 3.1) to assess the impact of varying the assumed population size of the zoonotic reservoir on the predicted number of reported cases. We conducted the test with populations of sizes 106, 5 × 106 and 2.5 × 107 where each was tested with 1/2, 1/20 and 1/100 of the population as being initially infected. We also tested with the original assumption of . We used the unaltered SMC algorithm with the tolerances set to be those calculated in the original run as seen in S1 appendix Table 1.
The results showed that the size of the zoonotic reservoir can be balanced by the transmission rate βrh and it changes with the inverse change of the population (S1 appendix Fig 1 and 2). I.e. 5 times more total rats meant that the posterior distribution of βrh was approximately a fifth of the original result as seen in 3 middle right.
The initial number of infected rats changed the beginning of the simulated case data but the fitted model was still able to capture the peaks of the epidemics in 2019 and 2020. The different initialisations of the case scenarios showed that the model was not flexible enough to fit to both the later epidemics and the first, smaller epidemic. This prevented case scenarios with 1/2 of their initial rat population being infected from completing more than 5 generations and scenarios with Ir(0) = Nr(0)/100 from completing more than 9 of the 11 generations.
We also completed a fit where we let Nr(0) = 1 and included the number of initial infected rats as a parameter to be fitted. Thus the zoonotic reservoir was represented in proportions of the population. The resulting fit (S1 appendix Fig 3 and 4) was still able to capture the seasonality of the data for Nigeria and the initial number of rats was estimated to be a median of 9.80% of the population (6.79–12.97%, 90% CI).
Discussion
Nigeria has experienced substantial Lf epidemics in recent years. The drivers of these epidemics have not, however, been previously explored in detail using a mechanistic model. The role of the rat reservoir in the seasonality of Lf cases is therefore not well understood. To that end, we developed and fitted a mechanistic model to national data of Lf epidemics in Nigeria from 2018 to 2020. This included previously over-looked rodent population dynamics for which the possibility of seasonal reproduction was investigated. We found that the model qualitatively replicated the weekly confirmed case data, showing that the seasonal peak in cases can be attributed to the population and epidemiological dynamics of the zoonotic rodent reservoir.
By fitting our model to case data from the outbreaks in Nigeria, we inferred that the recruitment rate of Mastomys rats in contact with at-risk Nigerians was highly seasonal, which lead to rapid and substantial shifts in the proportion of the reservoir that were susceptible to and infected with Lf. We found that there was a pulse increase of new susceptible rats in December of each year that quickly become infected by LASV-carrying rodents surviving the previous epidemic. These infected rats then spread the infection to human inhabitants in shared rural environments. This finding corroborates the results of Coetzee (1965) and Coetzee (1975) that showed a marked increase in Mastomys breeding during the rainy season [14, 15]. This drove the increase in observed cases in humans with over 95% transmission events being spillover events from the zoonotic reservoir. This proportion was larger than that seen in other studies, such as Iacono et al (2015), which may reflect the omission of within hospital interactions and the increased risk of transmission in hospitals without appropriate infection control precautions [16].
In our sensitivity analysis, we adjusted the size of the vector population and the initial proportion of infected rodents. The total size of the vector population could be factored into the rat-to-human transmission parameter by a proportionate rescaling (see S1 appendix Fig 2(a)-(c) and 4 mid right) and the initial proportion of infected rats altered the simulated number of cases in a predictable manner. Altering initial values for infected rats predictably changed the trajectory of the simulated number of observed cases for the first year (see S1 appendix Fig 1 and 3) however simulation of data for subsequent years was unaffected.
Our estimates for the basic reproduction number—the expected number of secondary infections from a primary infection in a fully susceptible population—of Lf in rats were well above that which one would consider to be an extremely infective disease, e.g. measles being an order of magnitude lower [34]. This corresponded to the sharp increase in rat reproduction rates in December, when Nigeria typically experiences its dry season. This is in conflict with other studies that have demonstrated a strong correlation between Mastomys reproduction and rainfall with a delay of 1–3 months after the rainy season [12]. It is important to note that the model’s representation of rat dynamics, particularly the recruitment function, may not capture the true complexity of the population dynamics. Using a recruitment function instead of a true birth function allows for a simplified life cycle of Mastomys by omitting a nesting/juvenile stage. Therefore dynamics that occur before recruitment and contact with humans may be missed, and could explain why the basic reproduction number for rats is higher than expected because potential susceptible individuals are excluded from transmission.
Incorporation of fine-scale space or further refinement of the rat population dynamics into the model—such as including time-dependent rat-to-rat transmission rates which may represent a hypothetical change in behaviour and proximity to humans—could improve both the model’s realism and enhance our understanding of rat dynamics during different climate seasons. For example, these developments may explicitly consider increased rat populations near homes in the dry season and higher LASV prevalence in Mastomys during the rainy season [11, 35]. Additionally, a metapopulation model would better capture spatial variations between different administrative areas, improving the representation of disease dynamics. To enhance comprehensiveness, future studies should incorporate hospitalization and treatment options specific to human cases, offering insights into healthcare worker risks and transmission reduction strategies. These enhancements would advance understanding of the disease and improve applicability of the model, thereby helping guide better strategies for Lf control and prevention in Nigeria.
Conclusion
Our model captured the dynamics of weekly confirmed Lf cases over multiple epidemics in Nigeria. Our approach demonstrated that the population of M. natalensis experienced annual Lf outbreaks due to seasonal recruitment rates resulting in an influx of new susceptible rats to fuel the spread of the disease. It is not yet clear what proportion of this recruitment is due to seasonal reproduction or migration to rural homes in the dry season. The high number of infected vectors causes a spillover of infection, resulting in annual epidemics between late December and early April. There are only relatively low levels of human-to-human transmission, which supports the notion that the zoonotic reservoir of M. natalensis is the primary driver in the epidemiology of Lf. Therefore, we conjecture that the single most effective measure of controlling the epidemics would be to reduce human contact with the zoonotic reservoir, either by increased food security and hygiene, or with more effective trapping and culling of rats.
Data Availability
NCDC situation reports for lassa fever are readily available at https://reliefweb.int/country/nga. Codebase and specific time-series data used can be found in the public github repository at https://github.com/JimmyMckendrick1/LassaFever_SeasonalityNigeria
https://reliefweb.int/country/nga
https://github.com/JimmyMckendrick1/LassaFever_SeasonalityNigeria
Supporting information
S1 Appendix. Appendix containing mathematical analysis and sensitivity analysis.
Acknowledgments
We thank the Nigeria Center for Disease Control for support and data supply.