A Sub-Microscopic Gametocyte Reservoir Can Sustain Malaria Transmission

Background Novel diagnostic tools, including PCR and high field gradient magnetic fractionation (HFGMF), have improved detection of asexual Plasmodium falciparum parasites and especially infectious gametocytes in human blood. These techniques indicate a significant number of people carry gametocyte densities that fall below the conventional threshold of detection achieved by standard light microscopy (LM). Methodology/Principal Findings To determine how low-level gametocytemia may affect transmission in present large-scale efforts for P. falciparum control in endemic areas, we developed a refinement of the classical Ross-Macdonald model of malaria transmission by introducing multiple infective compartments to model the potential impact of highly prevalent, low gametocytaemic reservoirs in the population. Models were calibrated using field-based data and several numerical experiments were conducted to assess the effect of high and low gametocytemia on P. falciparum transmission and control. Special consideration was given to the impact of long-lasting insecticide-treated bed nets (LLIN), presently considered the most efficient way to prevent transmission, and particularly LLIN coverage similar to goals targeted by the Roll Back Malaria and Global Fund malaria control campaigns. Our analyses indicate that models which include only moderate-to-high gametocytemia (detectable by LM) predict finite eradication times after LLIN introduction. Models that include a low gametocytemia reservoir (requiring PCR or HFGMF detection) predict much more stable, persistent transmission. Our modeled outcomes result in significantly different estimates for the level and duration of control needed to achieve malaria elimination if submicroscopic gametocytes are included. Conclusions/Significance It will be very important to complement current methods of surveillance with enhanced diagnostic techniques to detect asexual parasites and gametocytes to more accurately plan, monitor and guide malaria control programs aimed at eliminating malaria.


Introduction
Recent regional roll back of malaria from Zanzibar, Eritrea, Ethiopia and Rwanda has been achieved by up-scaled and improved malaria control efforts through distribution of long lasting insecticide treated bed nets (LLIN) and effective treatment with artemisinin combination therapy (ACT) [1][2][3]. These successes plant new hope for malaria elimination given necessary levels of support [4]. In the progress report of the UN Secretary-General's Special Envoy for Malaria from 2009, it was announced that the continent of Africa now had surpassed average LLIN coverage of 40% with numerous African nations exceeding 60% [5]. In the malaria-endemic Southwest Pacific region, LLIN coverage rates of 70% have been reported from Papua New Guinea (PNG) and the Solomon Islands [6]. The importance of continued surveillance of progress towards the goal of near zero preventable deaths from malaria by 2015 has been stressed [5]. However for successful economic and strategic planning of largescale malaria control, it is essential to be able to predict how long an intensified malaria control effort will have to be sustained before the threat of resurgent malaria transmission is removed.
Mathematical approaches, mainly based on fundamental work by Ross [7] and Macdonald [8], have been used to model malaria transmission and to predict the impact of malaria control during the malaria elimination campaigns in the 20th century [9]. These models utilize a population-based methodology whereby the relevant human population is divided into several subpopulations representing different states of infection. The subpopulations are usually called Susceptible, Exposed/Latent, Infective and Immune and they are connected by appropriate transitions [10][11][12][13][14][15][16]. Purely susceptible and immune states are rarely attained and therefore can be combined into a single 'partially immune -partially susceptible' subpopulation. The most convenient way to estimate the sizes of the subpopulations is to rely on observations of asexual parasite and gametocyte forms in peripheral blood. Since gametocytes are the sole parasite form possessing the ability to transmit malaria through the mosquito vector, only hosts harbouring these forms are part of an infective subpopulation. However, the fraction of individuals with detectable asexual parasites and gametocytes is dependent on the method of diagnosis. It has been estimated that approximately 50% of all infections with asexual parasites are not detected by standard light microscopic (LM) examination of thick blood films while for gametocytes, the percentage of false negative observations can be as high as 90% [17]. In a Kenyan study, for example, Bousema et al. reported that gametocyte prevalence observed by LM and quantitative nucleic acid sequence based amplification (QT-NASBA) were 22.3% and 91.1%, respectively [18], while a gametocyte prevalence of 7.6% by LM contrasted with 47% by magnetic deposition microscopy in a PNG survey [19]. Furthermore, recent studies have provided evidence that individuals who are not gametocyte positive by LM can still be infective. Gametocyte densities detectable by LM are usually reported to be .5-10 mL 21 but molecular and magnetic fractionation methods detect gametocytes at densities as low as 0.1 mL 21 [18,20].
The differences between gametocyte rates detected using advanced methods and LM as well as improved data on the relationship between submicroscopic gametocyte densities and mosquito infection may have substantial consequences for our understanding of malaria epidemiology. Whereas low-level gametocytemia corresponds to reduced likelihood of mosquito infection, the contribution that microscopically non-identifiable gametocyte carriers make to malaria transmission is still considerable [21][22][23]. Mathematical models which use gametocyte rates determined by LM can considerably underestimate the size of the infective compartment and thus the required scale and duration of malaria control needed to achieve elimination. There is a risk of premature termination of control measures followed by resurgence of disease.
To examine the effect of a subpopulation harboring submicroscopic gametocytes on predicted malaria transmission and predicted impact of malaria control two mathematical models were compared and contrasted. The first model (M1) does not include any observations of submicroscopic asexual parasites or gametocytes. The second model (M2) includes submicroscopic data. The approach presented here is substantially different from the classical Ross-Macdonald setup as it attempts to account for superinfection and for the differences in susceptibility and infectivity of the subpopulations. To estimate model parameters, available data from a low and a high transmission environment [24][25][26][27], as well as individual cases involved in malaria therapy (MT) studies [19,22,23,[28][29][30][31][32][33][34][35][36][37], were used. Several hypothetical control experiments simulating malaria control based on the distribution of LLIN were conducted using the calibrated models.
The present study shows that mathematical models based on field data collected by LM alone are likely to deliver overly optimistic predictions. Planning malaria control and monitoring the progress of interventions aimed at eliminating malaria would therefore benefit significantly from wide application of improved diagnosis of asexual parasites and especially gametocytes.

Model description
All symbols used in the models presented in this study are explained in Table 1.
The model, including only data collected by LM, from now on referred to as M1, consists of three subpopulations for the human host (X-uninfected, Yinfected with only asexual parasites, Zinfected with asexual parasite and gametocyte forms). The model which includes observations on sub-microscopic asexual parasite and gametocyte rates, subsequently referred to as M2, has an additional infective subpopulation (U -infected with asexual parasites and submicroscopic gametocyte forms) representing submicroscopic gametocyte carriers with a lower probability to infect mosquitoes (b U ). Both models are schematically shown in Figure 1. The differential equation systems for the human part of the models are given by equation systems (1) for M1 and (2) In equation systems (1) and (2) the following notations are used: lmosquito-to-human force of infection, 1/clatency period (duration until stage V gametocyte appearance), r,rrecovery rates of LM detectable and sub-microscopic gametocytemia, c Z ,c U -population fraction of Y that move to Z and U respectively. In both, M1 and M2, the part for the mosquito vector consists of three mosquito subpopulations (x-uninfected,y-latent,z-infective) with the usual transitions between them [12,38]. The delayed differential equation system (3) describes the mosquito part of the model.
In equation system (3), the following additional notations are used: L -human-to-mosquito force of infection, m -mosquito mortality rate, t -latency period and s -mosquito survival through the extrinsic incubation period. This mosquito part (Equations (3)) is often replaced by its quasiequilibrated infective prevalence S known as the sporozoite rate and given by equation (4).
The relationship between human biting rate (h), entomologic inoculation rate EIR and sporozoite rate S is given by equation (5).

EIR~hS ð5Þ
The mosquito-to-human force of infection (l) and the human-tomosquito force of infection (L) then become functions of human infectivity (B) and can be expressed as shown in equations (6) and (7) respectively.
L~vB ð7Þ In equations (6) and (7) a is the probability that an infectious mosquito bite results in human infection; v is the mean number of bites per mosquito per day and B is human infectivity as defined in the next section.

Graded human to mosquito infectivity
It has been shown that the probability of submicroscopic gametocyte carriers infecting mosquitoes is considerably reduced [22,39]. We attempted to account for this graded human infectivity by introducing of a second infective population (U) into M2 representing humans with submicroscopic gametocytemia. Therefore, whereas only one infective population (Z) exists in M1, there are two (Z and U) in M2 representing microscopically detectable and submicroscopic gametocyte carriers, respectively.
The probabilities of mosquito infection for these populations were termed b Z and b U , respectively, with b U ,b Z . Thus, mean infectivity B in M1 and M2 is given by equations (8) and (9), respectively.

Superinfection
To account for superinfection, the gametocyte clearance rates r and r in equations (1) and (2) were defined, as shown in equations (10).
, with w(i)~i In equations (10), r 0 and r 0 are the natural recovery rates without possibility of superinfection infection, a is the relative susceptibility to superinfection and w(i) is the superinfection function as described by Dietz [12]. Different levels of susceptibility to superinfection were assigned to each of the human subpopulations due to gradual loss of acquired immunity after an infection. The uninfected population (X) is most likely to be infected and its relative susceptibility compared to the other subpopulations is equal to the value of 1. Therefore, the full l is exerted on X. The subpopulation infected only with asexual parasites (Y) is assumed to be immune to superinfection and therefore has a susceptibility value of 0. The infective populations (Z and U) are susceptible to superinfection, but their recent infection provides them with greater protection. The relative susceptibility a for Z and U therefore lies between 0 and 1. The mosquito-to-human force of infection (l) exerted on these populations then becomes al. We assume that a is dependent on transmission intensity. In a low transmission setting a is high since the population will have developed only a small degree of immunity and it is thus more likely that superinfection is successfully established after an infective bite. In a high transmission setting a is lower since the degree of immunity will be higher. A table summarizing susceptibility, infectivity and clearance for both models can be found in Appendix S1.

Long lasting insecticide treated bed net control
For LLIN-based malaria control, the entire population was divided into a fraction 'q' benefiting from LLIN and a fraction '1-q' not benefiting from LLIN. LLIN create a barrier between human and mosquito and therefore reduce human biting rate (h) for the proportion q of the population using LLIN. Along with h, mosquito force of infection (l) is also is decreased, which will in turn affect clearance rates r and r. However the reduction in h is not 100%. It has been shown that use of LLIN is linked to regional, socio-economic and cultural circumstances, and to the perceived nuisance of mosquito bites and threat of malaria and other diseases [6,[40][41][42]. Therefore the parameter 'LLIN efficiency' (e), was introduced, which determines reduction of the original h for the LLIN possessing part (q) of the population as shown in equation (11).
The factor 'e' can be understood as a measure of compliance with LLIN usage and incorporates the fact that not all distributed LLIN will be used to their maximum capacity. The parameters 'q' and 'e' can be varied between 0 and 1 to mimic different malaria control scenarios. The human-to-mosquito force of infection is then represented by equation (12) where B is weighted for the LLIN parameters q and e.

L~v(eqBz(1{q)B) ð12Þ
Here, B is again defined as b Z Z in M1 and by equation (9) in M2.

Basic reproductive number
The basic reproductive number (R 0 ) sets apart model predictions of sustained versus disrupted malaria transmission. If R 0 ,1, a model predicts malaria eradication, otherwise sustained transmission. In the present study R 0 was analytically derived and computed for all model runs. LLIN usage was taken into account. Only the results by equations (13) for M1 and M2, respectively, are given (see Appendix S1 for details on the derivation of the equations).
The calculations in the present study were conducted using Wolfram Mathematica 7. The source codes can be downloaded by clicking on the following links http://www.cwru.edu/artsci/ math/gurarie/Malaria/High%20Transmission.nb and http:// www.cwru.edu/artsci/math/gurarie/Malaria/Low%20Transmission.nb for the two scenarios, respectively.

Model calibration
The equilibrium equations for the two models along with available epidemiological data and known transmission parameters were used to estimate unknown parameters specific to a malaria-endemic environment. This section only deals with the most important aspects of model calibration. Further information and equations can be found in Appendix S1.
Some model parameters were fixed at values or within ranges published in the literature or using estimates from malaria therapy data. Table 2 shows all values which were used as fixed inputs. . X -uninfected population with highest susceptibility to infection exposed to mosquito-to-human force of infection l. Y -population with only asexual stages present, not susceptible to superinfection and developing gametocytes at a rate c. A fraction of Y (c Z ) develops microscopically detectable gametocytes. Additionally in M2, a fraction c U develops sub-microscopic gametocytes. A certain fraction of Y never develops detectable gametocytes (1-c U -c Z ). Z and U are the populations with microscopically detectable and submicroscopic gametocytes respectively. They have a limited susceptibility to superinfection and are exposed only to a fraction of l, namely a (0,a,1). This limited susceptibility is incorporated into the models by using delayed clearance rates r and r for clearance of microscopically detectable and sub-microscopic gametocytes respectively. Thus, clearance rates r and r are functions of l and a as well as the original clearance rates without superinfection r 0 and r 0 as shown in equations (10). doi:10.1371/journal.pone.0020805.g001 Several studies have correlated probability of mosquito infection (b Z and b U ) with gametocyte density. We used data from these studies to estimate the infectivity of hosts with different gametocyte density ranges. The infectivity for subpopulations Z and U was  [22,23,32]. The gametocyte density corresponding to b U was ,10 mL 21 , which is the limit of conventional LM gametocyte detection [43].
Studies which report on microscopic and submicroscopic parasite (PR M and PR S ) and gametocyte rates (GR M and GR S ) and entomological data such as the entomological inoculation rate (EIR) are rare. We chose a low transmission setting where data were available from studies by Mwerinde (2005) and Shekalaghe et al. (2007) [24,27] which report on all of the above parameters for an area of low transmission in Moshi, Tanzania, and a high transmission setting in Burkina Faso studied by Paganotti et al. (2004 and [25,26]. Table 3 summarizes all values which were used from these two scenarios. To derive the sizes of the subpopulations (X,Y,Z,U) at endemic equilibrium from the data reported in [24][25][26][27], we made two assumptions: (1) all microscopically detectable parasite and gametocyte carriers are also detectable by the respective submicroscopic methods (i.e. there are no false negatives when microscopy is regarded the gold standard) and (2) all gametocyte carriers also have detectable asexual parasites and therefore contribute to the parasite rate. Furthermore M1 uses exclusively data collected by microscopy (microscopic parasite rate (PR M ) and microscopic gametocyte rate (GR M )) and M2 uses data collected by microscopy and more sensitive methods (submicroscopic parasite rate (PR S ) and submicroscopic gametocyte rate (GR S )). This resulted in equations (14) and (15)  The population fractions c Z and c U were used for calibration using the equilibrium equations of M1 and M2 (the equilibrium equations and the entire calibration scheme can be found in Appendix S1).
Little is known about the relative susceptibility (a). The higher the value of a, the more likely is superinfection in the respective population. We assume that a is high in a low transmission area and low in a high transmission area. Thus we chose a = 0.8-0.9 for model calibration with the data from [24,27] and a = 0.01-0.05 for model calibration with data from [25,26].
We did not aim to perform a best fit through more than two data sets as this was beyond the scope of the present work and very few studies report on all necessary input parameters. Furthermore the models presented here are simplistic and relatively inflexible. However, fixed and calibrated parameters should lie within reasonable boundaries.
The model calibration was run 1000 times with starting values chosen from within suitable ranges (also shown in Tables 3 and 4)  Table 3. Model input values from the low transmission scenario described by Shekhalage (2007) and Mwerinde (2005), and the high transmission scenario described by Paganotti (2004 and [24][25][26][27]. to obtain the mean and a minimum to maximum range of predictions. Parasite Rate (PR) in our plots corresponds to either Y+Z in M1 or Y+Z+U in M2. The cut-off PR that was set to indicate malaria elimination was 10 24 (,1 case of malaria in 10,000 people).

Results
The basic reproductive number (R 0 ), which is essential for sustained infection, showed marked differences between M1 and M2 in both high and low transmission settings. The predicted R 0 at endemic equilibrium (without LLIN control) was also considerably higher for M2 in both the low (R 0 for M1: 1.001-1.025 vs. R 0 for M2: 1.47-1.64) and high (R 0 for M1: 3.8-10.4 vs. R 0 for M2: 20.6-57.00) transmission settings. Because of the different stability levels of the endemic equilibria in the two models, similar malaria control efforts will have different, modeldependent, predicted impacts. While a defined control attempt could disrupt transmission for M1, bringing R 0 to values ,1 (unsustainable transmission), M2 could still predict persistent transmission at another equilibrium state with R 0 .1. Figure 2 shows specific examples for cases of malaria control with LLIN in which M1 predicts malaria eradication while M2 predicts persistent transmission. Figure 2A shows the scenario of low transmission with LLIN coverage of 50% and usage of 40%. These values resemble LLIN coverage and usage achieved across Africa and other endemic areas in recent years [5,6]. In this case M1 predicts disrupted malaria transmission and malaria eradication in 1.9 years (R 0 is reduced to 0.68-0.69 in M1), whereas M2 predicts persistent transmission at a lower level (R 0 is reduced to 1.00-1.11 in M2). In the high transmission setting these values of LLIN coverage and usage do not result in predicted malaria eradication by either M1 or M2. Therefore, a hypothetical control scenario set in the high transmission environment is presented in Figure 2B, in which LLIN coverage is 92% and LLIN usage is 88%. In this case M1 predicts malaria eradication from the high transmission area in approx. 6.8 years (R 0 0.35 to 0.97 is reduced to in M1), whereas M2 only predicts a sustained 30% reduction in malaria prevalence (R 0 is reduced to 1.9-5.3 in M2).
Whereas Figure 2 shows parasite rates for specific malaria control scenarios, Figure 3 shows R 0 for all possible combinations of LLIN coverage (0,q,1) and LLIN usage (0,e,1) calculated using equations (13). It can be seen that much higher levels of LLIN coverage and usage are required as predicted by M2 compared to M1 to reduce R 0 to values below 1.

Discussion
The present study emphasizes the important role that host reservoirs with low-level gametocyte and asexual parasite density can play in sustaining malaria transmission despite LLIN-based malaria control. For any given control scenario, the model which includes a large pool of hosts with sub-microscopic gametocyte and asexual parasite densities predicts considerably higher transmission and therefore a need for greater coverage and adherence to LLIN usage. When LLIN coverage is set to levels at or around values currently achieved with the up-scaled efforts enabled by The Global Fund to Fight AIDS, Tuberculosis and Malaria (GFATM) [5], only M1 in a low transmission setting predicts malaria eradication. M2 does not predict interrupted transmission even in the low transmission setting (Figure 2A). For Table 4. Parameters used in M1 and M2, which were either calibrated or derived from [25][26][27].

Value Low Transmission Scenario
High Transmission Scenario Subpopulations X,Y,Z and U were derived from the parasite rates (PR) and gametocyte rates (GR) shown in Table 3 using equations (14) and (15). Mosquito-to-human force of infection (l) was calculated using equation (S15) Clearance rates r and r are prolonged by superinfection as shown in equations (8). Fractions c Z and c U were used for calibration using equations (S16) and (S17) which can be found in Appendix S1. *although a result of the calibration process (equation S16), we constrain the probability that an infective bite results in human infection (a) to lie between 0.05 and 0.5, which we assume to be a reasonable range for a highly uncertain parameter such as a. high transmission environments, the models presented here predict the need for much higher LLIN coverage and usage. Furthermore our predictions emphasize the importance of implementing sensitive diagnostic approaches that are used for monitoring and evaluating malaria control strategies and the progress of malaria control campaigns. Although previous models acknowledge submicroscopic gametocytemia [11,44], our current approach is, to the best of our knowledge, the first to evaluate the impact of these infections on malaria transmission and control. In our models, LLIN is the only malaria control strategy. In real-life field situations, efficient malaria control utilizes LLIN as part of a range of co-ordinated strategies. The next phase of worldwide intensified malaria control aims to build on the success of LLIN distribution in 2008 by deploying rapid diagnostic tests and effective antimalarial therapies as widely as possible in order to reduce preventable deaths from malaria to near zero by 2015 [5]. It should be noted that efficient treatment of symptomatic malaria cases, while undoubtedly an urgent priority, will only partially reduce malaria transmission because asymptomatic carriers will continue to transmit the disease'.
Adult populations in holo-or hyperendemic regions have high levels of acquired immunity to the severe manifestations of the infection. Although parasitemic, these people do not usually feel ill and thus see no reason to visit a health center or seek treatment, but they are likely to harbor infectious gametocytes that will sustain malaria transmission regardless of a successful local reduction in symptomatic malaria. Asymptomatic malaria infection is normally characterized by very low asexual parasite densities and the absence of gametocytes, especially in adults [45][46][47]. A recent review by Okell et al. [17] revealed that standard LM failed to detect asexual malaria parasites in 49.2% of all cases and 91.3% of gametocytemic cases. Furthermore, Bonnet et al. estimate the contribution of hosts with undetectable gametocytemia to malaria transmission at 23.7% [21,22] while Schneider et al. determined an even higher value at 54.8% [21,22]. Sustained active case detection using rapid diagnostic tests and LM will, therefore, be inadequate because the sensitivity of these methods is too low. It is important to maximize detection of cases infectious to mosquitoes for a prolonged period, presumably decades, to achieve true malaria eradication [48].
If a fraction of malaria cases go undetected and the intensified efforts of malaria control are reduced prematurely, malaria will increase such as it did in India [49], Sri Lanka [50], Zanzibar [50] and Madagascar [51] after the halt of the eradication programs in the last century. Resurgent malaria can have devastating effects on populations with levels of immunity that have waned following partially effective control. For example, the incidence of malaria in the highlands of Madagascar fell dramatically during the control program in the 20 th century but more than 40,000 fatalities were attributed to malaria during the 5 years after the control program was stopped [51].  Figure 2. Note that R 0 = 1 sets apart 'sustained transmission' from 'malaria eradication'. It can clearly be seen that in Panels A and C (M1) R 0 is well below this discriminating value of 1. In Panels B and D (M2) the values for R 0 are clearly below the value of 1. Therefore, while M1 predicts malaria eradication, M2 still predict persistent transmission. doi:10.1371/journal.pone.0020805.g003 It will be very important to complement current methods of surveillance with enhanced diagnostic techniques to more accurately plan, monitor progress and guide malaria control which is aimed at eliminating malaria. Recently developed methods to detect gametocytes at low densities such as RT-PCR, QT-NASBA and HFGMF can reduce the threshold of gametocyte detection to cover a wider epidemiologically relevant range. Application of these methods and their potential refinement to maximize time and cost effectiveness appears a vital component of successful malaria elimination.
It is particularly important to detect gametocytes because they are more resilient to treatment. While there is evidence of a gametocytocidal effect of artemisinin drugs, gametocyte clearance times are still very long in ACT-treated patients with a mean of 49 h for asexual stages and 220 h (.9 days) for gametocytes in a study of 559 patients conducted in Thailand in from 1998 to 2006 [52]. Another recent study detected submicroscopic gametocytes up to 48 days after ACT treatment [53]. Primaquine, the most potent gametocytocidal drug, can provoke hemolysis in glucose-6phosphate dehydrogenase (G-6-PD) deficient patients and is thus not used routinely in the many tropical countries in which G-6-PD deficiency has been observed [54,55]. This means that some patients can be infective for weeks following clearance of asexual parasites even after ACT treatment. ACT or other treatments were not included in the current analyses, in part because of factors such as drug-specific post-treatment gametocyte viability [56,57] and a variable if brief period of carriage of gametocytes without asexual forms. We believe that the problem of submicroscopic gametocytemia sustaining malaria transmission is mostly relevant in asymptomatic malaria carriers. Since they will rarely seek treatment, ACT/other therapy will not have a great impact on this part of the population. The exception may be mass drug administration which is not likely to play a big role in malaria control efforts, especially with the fear of emerging artemisinin resistance [58].
The aim of the present study was to highlight the importance of sensitive gametocyte and asexual parasite detection. For this we have used data from two very different transmission settings to calibrate the models. Other calibration procedures that constrain the range of admissible parameters could also be used. The models presented here have several limitations. Ross-Macdonald methodology is a very simplistic approach, since the real-world dynamics of malaria infection depend on many more variables (such as the age distribution of the population, seasonality of transmission and protective host factors including red cell polymorphisms) that cannot be encompassed by such models. Nevertheless, our adaptation of the most commonly used form of mathematical models for malaria can be viewed as a necessary refinement in their further development and an important step in identifying practical strategies for parasite detection necessary to enhance control and elimination efforts.

Supporting Information
Appendix S1 The appendix contains more detailed information on equilibrium equations, model calibration, graded susceptibility and infectivity and the analytical derivation of the basic reproductive numbers.