Is a reproduction number of one a threshold for Plasmodium falciparum malaria elimination?
© The Author(s) 2016
Received: 1 May 2016
Accepted: 10 July 2016
Published: 26 July 2016
The basic reproduction number (R 0) is an important summary of the dynamics of an infectious disease. It is a threshold parameter: an infection can only invade a population if R 0 is greater than 1. However, a number of studies using simple models have suggested that for malaria, it is in theory possible for infection to persist indefinitely even if an intervention has reduced R 0 below 1. Such behaviour is known as a bistable equilibrium. Using two published mathematical models which have both been fitted to detailed, age-stratified data on multiple outcomes, the article investigates whether these more complex models behave in such a way, and hence whether a bistable equilibrium might be a real feature of Plasmodium falciparum malaria in Africa.
With the best-fitting parameter values, neither model has a bistable state, because immunity reduces onwards infectiousness. The results imply that there is a threshold such that if interventions can reduce transmission so that R 0 is below 1 for long enough, then malaria will be locally eliminated.
This means that calculations of the reduction in R 0 that interventions can achieve (the effect size) have a useful and straightforward interpretation, whereas if the theoretical possibility of a bistable equilibrium were the real behaviour, then such effect size calculations would not have a clear interpretation.
KeywordsMalaria Mathematical model Elimination Reproduction number
The basic reproduction number R 0 of an infectious disease is defined as the number of secondary cases produced by a typical infected case in an otherwise susceptible population. It is important as a threshold quantity: if R 0 is less than 1 then it is not possible for an infection to invade a susceptible population. The same threshold in many cases also applies in the other direction: starting from an endemic state, if R 0 is reduced below 1 and kept there then the infection will die out. However a number of studies have shown that for malaria, in theory it is possible that there are stable endemic states with R 0 below 1 that can persist indefinitely. This phenomenon is known as a bistable equilibrium since there can be a set of model parameters for which the disease-free equilibrium and an endemic equilibrium are both stable.
Águas et al.  formulated a model of malaria transmission with immunity, where the effect of immunity is to reduce the probability of clinical malaria. They used the model to show that if immune individuals have a longer duration of infection, and hence more overall onwards infectiousness than non-immunes, then there are parameter regions with a bistable equilibrium. Keegan and Dushoff  further explored the dynamics of a similar model, and derived a quantity based on the ratio of reproduction numbers with and without immunity that is a threshold for whether there is a bistable equilibrium. In a review of countries which had eliminated malaria, Chikaya et al. found that once elimination is reached, a resurgence to high endemic levels is rare, whereas it is common in countries that pushed malaria to low levels without eliminating it [3, 4]. Smith et al.  investigated this stability of elimination, looking at the evidence for several possible explanations. One explanation was a bistable equilibrium: it could be that in some endemic settings R 0 was below 1, perhaps having been reduced previously by changes such as improved housing and increased levels of treatment. Then if interventions further reduced transmission for some years, malaria could be eliminated and the population would lose their immunity. When the interventions are removed, the reproduction number will return to the same value it had when there was endemic transmission, but that value is below 1 and so if malaria is imported, it will die out.
In these models of malaria transmission with a bistable equilibrium, partial immunity to malaria has the effect of increasing an individual’s onward transmission over the course of an infection. The most plausible mechanism for this is that immunity reduces the probability of clinical symptoms. With no symptoms and hence no treatment, the infection can persist for longer, whereas effective treatment clears the infection. However, immunity to malaria not only reduces the probability of symptoms, but also lowers parasite densities and hence reduces onwards infectiousness per unit time, and so it is not clear a priori whether overall, immunity increases or reduces onwards transmission. This article investigates whether models which capture immunity in a more realistic way show bistability, using two mathematical models which both incorporate immunity at several stages, and which have both been fitted to a wide range of types of data on Plasmodium falciparum malaria in Africa, but which have quite different internal structures.
Two models of human infection and immunity are considered. The first is the model of Griffin et al. [6, 7], which is referred to here as the Griffin model. It is implemented as both a deterministic compartmental model that was fitted to data and a corresponding individual-based model for flexibility in looking at the impact of interventions: the former is used for the results here. The second model is the OpenMalaria model described in [8–11] and summarized with updated fitted parameter values in .
Λ M is the force of infection experienced by mosquitoes, and S M , E M and I M are the numbers (per person) which are susceptible, latently infected and infectious, respectively.
Reproduction number formulae
α is the biting rate on humans and γ = e−μτ ⁄μ is the expected duration of infectiousness of an infected mosquito. b is the probability of infection if bitten by an infectious mosquito. C is the expected overall human infectiousness to mosquitoes if infected, integrating over the human infectious period. For example if there is a human infectious period of average length d with a constant probability c of infecting each susceptible mosquito that bites during this time, then C = dc. The vectorial capacity is V = mγα 2, and R 0 = VbC.
For all the results in this paper, the reproduction number may include treatment for symptomatic malaria, which changes the onward infectiousness C by clearing some infections. As treatment could be considered as an intervention, the term R C for “reproduction number under control” could be used in these cases, but for simplicity the term R 0 has been used throughout. It is necessary to include treatment in the model as it is treatment which in some models may lead to bistable behaviour.
Deriving R0 from the endemic equilibrium
It is assumed that the quantity which varies to give different values for R 0 and the EIR is the mosquito density per person m. For a given EIR ε, the human equilibrium states can be found. For the Griffin model, which is a compartmental model stratified by age and biting rates, this can be done analytically, as described in Additional file 1. For the individual-based OpenMalaria model it is done by running the simulation model with a forced, fixed EIR over one human lifetime with a large enough population size that stochastic variation is negligible.
Note that Λ M is found from the endemic equilibrium, but E H (α 2 bC) is for people with no immunity. ε here is the mean EIR experienced by the population, not the EIR for adults.
Details of R0 for the Griffin model
δ a and δ h are corrections due to the variation in biting rates by age and other heterogeneity which both increase R 0.
Taking 1/η as 21 years, to approximate the cross-sectional age distribution in sub-Saharan Africa, with ρ = 0.85 and a 0 = 8 years, gives δ a = 1.1.
So this amount of heterogeneity in biting increases R 0 by more than fivefold for a given value of m.
Each r is the recovery rate from the corresponding model state, and each c is the infectiousness of a state, with subscripts T, D, A and U referring to the states shown in Fig. 1. With no immunity, c A = c D .
Equations (4) and (5) ignore the fact that during an infection, individuals age and so the rate at which they are bitten by mosquitoes changes. Correctly accounting for ageing during an infection changes R 0 by around 1 %, but greatly complicates the formulae when there are multiple infection states since it removes the separation between time since infection and age. The formulae that correctly account for ageing are given in Additional file 1, and were used for the results.
The force of infection on mosquitoes at the endemic equilibrium is found by calculating the equilibrium proportions in each state, stratified by age and exposure to mosquitoes, and summing their infectiousness to mosquitoes, weighted by the relative biting rates.
Checking stability in the Griffin model
Once the equilibrium states were found in the Griffin model for a given EIR ε 0, the equilibrium was checked for stability by randomly perturbing the model states by a small amount and then numerically solving the model differential equations to see whether the system returned to the staring equilibrium or to a different equilibrium. At equilibrium, all the human and mosquito model states were multiplied by 1 + z, where z is a random draw from a normal distribution with mean zero and standard deviation 0.01 (different for each age group, heterogeneity level and state). The human model states were renormalized to have the correct age distribution. The differential equations were then run for 2000 years, and the EIR implied by the final model states was calculated, denoted by ε 1: such a long time is needed because near the boundary between stable and unstable regions, the model may be extremely slow to either move away from or return to the original equilibrium. For each parameter set and starting EIR this was repeated five times. If the maximum value over the five trials of |log(ε 1/ε 0)| was less than 0.003, then the equilibrium was considered stable, otherwise it was unstable. This criterion was chosen by trial and error as it discriminates well between stable and unstable regions.
Details of R0 for the OpenMalaria model
In the OpenMalaria model, as well as variation in the biting rate by age and in some model variants between people, C also varies by age: there is maternal immunity which reduces the parasite density and which is taken to be independent of transmission intensity, and so for the model as implemented this is present even in populations with no prior exposure; the incidence of severe malaria and its case fatality also vary by age, even with no immunity, and the resulting mortality reduces the expected infectiousness; and for R 0 with treatment, there is additional variation in C because the probability of clinical malaria and hence treatment varies by age due to an age-dependent parasite density threshold determining the probability of symptoms (the pyrogenic threshold).
With no immunity, all infectious bites result in infection in this model when only considering a single infectious bite per person, which is the case when calculating R 0, and so b = 1. (The probability of infection per infectious bite is assumed to decrease if people receive many infectious bites per unit time). C was found for a range of age groups by running the simulation model with no transmission, infecting all individuals at time 0, and following them over time to calculate their infectivity to mosquitoes at each time-step (but still not allowing any onward transmission), storing the output by date of birth. Variation in the biting rate by age and due to other factors were included by outputting these quantities for each individual in the population, and then these were combined to find E H (α 2 bC).
For more complex models with multiple types of immunity, there is no simple formula that will tell us whether there is a bistable region. However it is still possible to see how a model behaves by plotting R 0 against the EIR as in Fig. 2. Two models of Plasmodium falciparum malaria are considered which include several kinds of immunity and which have both been fitted to detailed, age-stratified data on multiple epidemiological outcomes. The first is the model described of Griffin et al. . This model is referred to as the Griffin model. The other model is the individual-based OpenMalaria model developed at the Swiss Tropical and Public Health Institute [8–12], whose source code has been made available for download .
For both of these models, there is a unique equilibrium state corresponding to any given EIR. It is assumed that the quantity which varies to give different values for the EIR and R 0 is the mosquito density per person m. There is always a single value of m and R 0 for any given EIR, but there may be more than one EIR that implies the same value of m and R 0 (while keeping the model parameters other than m fixed). If R 0 is a non-monotonic function of the EIR then there is a bistable region; whereas if R 0 is an increasing function of the EIR and is above 1 for any positive EIR then there is no bistable region.
The results here use the same mosquito model combined with both human models, and so the calculated values of R 0 for the OpenMalaria model will differ a little from those implied by the mosquito model developed by Chitnis et al. . However this will not affect the qualitative pattern of how R 0 varies with the EIR, which is determined by human immunity. All results are for non-seasonal settings, and R 0 is calculated with the same probability of treatment for clinical malaria as is used to find the endemic equilibrium states.
R 0 monotonically increases with EIR;
R 0 increases, decreases and then increases, but is always above 1 for any positive EIR;
R 0 increases, decreases and then increases, and goes below 1 for some positive EIRs.
A small number (fewer than 0.1 %) of parameter sets had additional small wiggles in these curves. The posterior probability of each pattern can be calculated by counting how many parameter sets show each type of behaviour. If the proportion of clinical cases treated is up to 80 %, then the posterior probability that R 0 monotonically increases with EIR is over 99 % (Fig. 4d). The posterior probability of pattern 3, where R 0 is below 1 for some EIRs, is 1, 4 or 13 % if 90, 95 or 100 % of cases are effectively treated, respectively. The stability of each equilibrium state was tested as described in the “Methods” section. For all three patterns, regions where R 0 was increasing as a function of EIR were stable and those where R 0 was decreasing were unstable. Both patterns 2 and 3 have a bistable region, i.e. a range of values of R 0 with two stable states. For R 0 > 1, the bistable region consists of two stable endemic states, whereas for R 0 < 1 it consists of a stable endemic state plus the infection-free equilibrium.
For all parameter sets and all treatment levels, in this model R 0 increases above 1 for small EIRs, i.e. there is a forward bifurcation, even for parameter sets that do have a bistable region. So a local stability analysis at small EIRs around the region R 0 = 1 would not have revealed any bistable behaviour. This is in contrast to the models illustrated in Fig. 2, in which R 0 either increases monotonically (a forward bifurcation with no bistable region) or initially decreases below 1 at low EIRs before then increasing (a backward bifurcation and with a bistable region), and so a stability analysis at small EIRs would detect the bistability.
Both the Griffin and OpenMalaria transmission models also show a decrease in infectiousness with age (Fig. 6c, d). Only the base OpenMalaria model is plotted in this figure, but the other model variants show the same pattern. The Griffin model has an empirical function to capture immunity that reduces the detectability of infections, and hence reduces the recorded parasite prevalence. Infectiousness is modelled as a transformation of the probability of detection, with the model fitted to the feeding studies plotted in Fig. 6a as well as to data on parasite prevalence by age from a much wider range of settings. Infectiousness in the OpenMalaria model is modelled differently, and is based on different data: this model explicitly tracks parasite densities as they change during an infection, along with the reduction in these densities due to immunity. Infectiousness to mosquitoes as a function of asexual parasite densities was fitted to data from malaria-therapy studies, in which non-immune patients were infected with malaria as a treatment for syphilis, and also took part in mosquito feeding studies .
Previous work has shown that for malaria, in theory it is possible for R 0 to be reduced below 1 but for infection to persist indefinitely, if there is a bistable equilibrium. The results presented here suggest that P. falciparum does not behave in this way, and the most likely behaviour is that R 0 increases monotonically as the EIR increases, and is above 1 for all positive EIRs. This implies that the value R 0 = 1 is a threshold for malaria elimination, i.e. starting from an endemic state, if interventions can reduce R 0 below 1 and keep it there for long enough, then the infection will be locally eliminated unless there is re-importation.
The uncertainty in these conclusions can be assessed for the Griffin model using the posterior distribution of human model parameters. For treatment levels above 80 %, the results become less certain, but there is only a substantial posterior probability (above 5 %) of the existence of endemic states with R 0 below 1 if more than 95 % of clinical cases are effectively treated, and even with 100 % treatment a monotonically increasing R 0 is more likely. For the OpenMalaria model the uncertainty was assessed in a different way, as there are 14 different model variants. All of these behaved in the same way, with R 0 an increasing function of the EIR for all treatment levels, and so there were no endemic states with R 0 below 1.
The results in this paper include the same probability of treatment for symptomatic malaria in the calculation of R 0 as is used to calculate the endemic equilibrium states. The reason for this is that treatment for symptomatic malaria is conceptually different to other interventions: one would try to treat people’s symptoms regardless of the transmission setting, but would take the setting into account when considering other interventions. For example when moving towards elimination, if malaria has been locally eliminated in a certain area, policy makers may consider scaling back some interventions and switching to a more reactive policy. But they would not scale back treatment for symptomatic malaria, and so the reproduction number with treatment is the relevant metric.
The methods used here rely on finding the model equilibrium state. Hence the results are a simplification as they do not account for the dynamic changes over time in immunity that would occur if transmission was rapidly reduced . Depending on how immunity alters someone’s onward infectiousness, it is possible that there would be a period after transmission has been reduced when the population still has immunity from the pre-intervention period and that elimination would be possible during this time, even if the reproduction number with the intervention is greater than 1. As the results are for equilibrium states, they also do not include seasonality in transmission.
Of the previously published models which suggested there could be a bistable equilibrium, the model of Águas et al. was fitted to data on the age patterns of severe malaria, while the models of Keegan and Dushoff and Smith et al. were not fitted to data but explored the model behaviour over a range of parameter values [1, 2, 5]. None of the models were fitted to data on infectiousness to mosquitoes, or used any data on parasite prevalence or parasite densities as these vary by age and transmission intensity. All three models showed regions of parameter space with and other regions without a bistable equilibrium. These models all have a single type of immunity. This is modelled in different ways but in all cases, the average duration of infection with immunity is increased compared to non-immunes. Keegan and Dushoff also allow the relative infectiousness of immunes and non-immunes to vary when exploring the model behaviour. They derive a threshold that determines whether or not there is a bistable equilibrium: this occurs if immunes are more infectious over the course of their infection than non-immunes by a sufficiently large margin. The model of Águas et al. behaves in a similar way. It is difficult to determine directly from currently available data whether the assumptions that lead to bistable behaviour are satisfied.
The results in the present paper are based on two models which include immunity at multiple points, and have been fitted to a wide range of different types of data. In particular, as well as having immunity which reduces the probability of symptoms, both models allow for a lower infectiousness per unit time in infected people with some immunity compared to infected people with no immunity. The model structures for onward infectiousness to mosquitoes and how this varies with past exposure and for immunity to clinical malaria differ between the two models: these are the main aspects that determine whether there is a bistable region. The OpenMalaria model explicitly tracks individuals’ parasite densities, which are decreased by immunity. Symptomatic malaria occurs with a certain probability that depends on the parasite density, and the probability of infecting a mosquito also depends on the recent parasite density, and so both of these probabilities are reduced because immunity reduces the parasite density. The Griffin model instead uses separate empirical functions for each kind of immunity. The data the models were fitted to are also different, in particular for onward infectiousness to mosquitoes. These differences make the findings more robust than if the results were only supported by a single model.
For understanding whether there is a bistable equilibrium, some information about immunity when there is low transmission intensity is needed. The datasets that the two models were fitted to were from a range of transmission intensities, although there were more sites with moderate to high transmission. The Griffin model was fitted to data on the age patterns of clinical malaria and of parasite prevalence in sites with parasite prevalence in 2–10 years olds ranging from below 1 to 90 %. Furthermore, people acquire immunity over time as their cumulative exposure increases, and so since there was age-stratified data on multiple end-points, these age patterns carry some additional information about how immunity changes with exposure. Most importantly, older children and adults are less infectious to mosquitoes than young children, averaging over those with detectable infections and those without. This suggests that in terms of the overall contribution to onwards infectiousness, the effect of immunity in reducing parasite densities and hence reducing onward infectiousness when infected, outweighs the effect of a reduction in the probability of symptoms and hence of the infection being cleared by treatment.
These results are important for using mathematical models to assess the contribution that interventions could make to elimination. The relative change in R 0 due to interventions is known as the effect size [21, 22], and is often straightforward to calculate. The theoretical studies which suggested there was bistable behaviour and stable equilibrium states with R 0 below 1 cast doubt on the usefulness of such effect size calculations because if that prediction was correct then R 0 = 1 would not be a threshold for elimination.
The author declare that he has no competing interests.
Availability of data and material
The work was funded by a fellowship from the UK Medical Research Council, Grant No. G1002284. http://www.mrc.ac.uk/index.htm. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Águas R, White LJ, Snow RW, Gomes MGM. Prospects for malaria eradication in sub-Saharan Africa. PLoS ONE. 2008;3:e1767.View ArticlePubMedPubMed CentralGoogle Scholar
- Keegan LT, Dushoff J. Population-level effects of clinical immunity to malaria. BMC Infect Dis. 2013;13:428.View ArticlePubMedPubMed CentralGoogle Scholar
- Chiyaka C, Tatem A, Cohen J, Gething P, Johnston G, Gosling R, et al. The stability of malaria elimination. Science. 2013;339:909–10.View ArticlePubMedGoogle Scholar
- Cohen JM, Smith DL, Cotter C, Ward A, Yamey G, Sabot OJ, et al. Malaria resurgence: a systematic review and assessment of its causes. Malar J. 2012;11:122.View ArticlePubMedPubMed CentralGoogle Scholar
- Smith DL, Cohen JM, Chiyaka C, Johnston G, Gething PW, Gosling R, et al. A sticky situation: the unexpected stability of malaria elimination. Philos Trans R Soc Lond B Biol Sci. 2013;368:20120145.View ArticlePubMedPubMed CentralGoogle Scholar
- Griffin JT, Ferguson NM, Ghani AC. Estimates of the changing age-burden of Plasmodium falciparum malaria disease in sub-Saharan Africa. Nat Commun. 2014;5:3136.View ArticlePubMedPubMed CentralGoogle Scholar
- Griffin JT, Hollingsworth TD, Okell LC, Churcher TS, White M, Hinsley W, et al. Reducing Plasmodium falciparum malaria transmission in Africa: a model-based evaluation of intervention strategies. PLoS Med. 2010;7:e1000324.View ArticlePubMedPubMed CentralGoogle Scholar
- Smith T, Maire N, Dietz K, Killeen GF, Vounatsou P, Molineaux L, et al. Relationship between the entomologic inoculation rate and the force of infection for Plasmodium falciparum malaria. Am J Trop Med Hyg. 2006;75(2 suppl):11–8.PubMedGoogle Scholar
- Maire N, Smith T, Ross A, Owusu-Agyei S, Dietz K, Molineaux L. A model for natural immunity to asexual blood stages of Plasmodium falciparum malaria in endemic areas. Am J Trop Med Hyg. 2006;75(2 suppl):19–31.PubMedGoogle Scholar
- Ross A, Killeen G, Smith T. Relationships between host infectivity to mosquitoes and asexual parasite density in Plasmodium falciparum. Am J Trop Med Hyg. 2006;75(2 suppl):32–7.PubMedGoogle Scholar
- Smith T, Ross A, Maire N, Rogier C, Trape J-F, Molineaux L. An epidemiologic model of the incidence of acute illness in Plasmodium falciparum malaria. Am J Trop Med Hyg. 2006;75(2 suppl):56–62.PubMedGoogle Scholar
- Smith T, Ross A, Maire N, Chitnis N, Studer A, Hardy D, et al. Ensemble modeling of the likely public health impact of a pre-erythrocytic malaria vaccine. PLoS Med. 2012;9:e1001157.View ArticlePubMedPubMed CentralGoogle Scholar
- Diekmann O, Heesterbeek H, Britton T. Mathematical tools for understanding infectious disease dynamics. Princeton: Princeton University Press; 2012.View ArticleGoogle Scholar
- OpenMalaria. 2014 [cited 2014 3rd November]. https://code.google.com/p/openmalaria/.
- Chitnis N, Schapira A, Smith T, Steketee R. Comparing the effectiveness of malaria vector-control interventions through a mathematical model. Am J Trop Med Hyg. 2010;83:230–40.View ArticlePubMedPubMed CentralGoogle Scholar
- Boudin C, Olivier M, Molez J-F, Chiron J-P, Ambroise-Thomas P. High human malarial infectivity to laboratory-bred Anopheles gambiae in a village in Burkina Faso. Am J Trop Med Hyg. 1993;48:700–6.PubMedGoogle Scholar
- Bonnet S, Gouagna LC, Paul RE, Safeukui I, Meunier JY, Boudin C. Estimation of malaria transmission from humans to mosquitoes in two neighbouring villages in south Cameroon: evaluation and comparison of several indices. Trans R Soc Trop Med Hyg. 2003;97:53–9.View ArticlePubMedGoogle Scholar
- Githeko AK, Brandling-Bennett AD, Beier M, Atieli F, Owaga M, Collins FH. The reservoir of Plasmodium falciparum malaria in a holoendemic area of western Kenya. Trans R Soc Trop Med Hyg. 1992;86:355–8.View ArticlePubMedGoogle Scholar
- Churcher TS, Bousema T, Walker M, Drakeley C, Schneider P, Ouédraogo AL, et al. Predicting mosquito infection from Plasmodium falciparum gametocyte density and estimating the reservoir of infection. Elife. 2013;2:e00626.View ArticlePubMedPubMed CentralGoogle Scholar
- Smith T, Schapira A. Reproduction numbers in malaria and their implications. Trends Parasitol. 2012;28:3–8.View ArticlePubMedGoogle Scholar
- Smith DL, Hay SI, Noor AM, Snow RW. Predicting changing malaria risk after expanded insecticide-treated net coverage in Africa. Trends Parasitol. 2009;25:511–6.View ArticlePubMedPubMed CentralGoogle Scholar
- Griffin JT. The interaction between seasonality and pulsed interventions against malaria in their effects on the reproduction number. PLoS Comput Biol. 2015;11:e1004057.View ArticlePubMedPubMed CentralGoogle Scholar
- R Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2013.Google Scholar