1 Institut Pasteur de Dakar, Dakar, Senegal.
2 Institut National de la Santé et de la Recherche Médicale (INSERM), Unité 444, Faculté de Médecine Saint-Antoine, Université Pierre et Marie Curie, Paris, France.
3 Institut de Médecine TropicaleService de Santé des Armées, Marseille, France.
4 Institut de Recherche pour le Développement (formerly ORSTOM), Dakar, Senegal.
5 Institut National de l'Environnement Industriel et des Risques, Verneuil en Halatte, France.
![]() |
ABSTRACT |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
disease transmission; immunity; malaria; models; statistical; Plasmodium falciparum
Abbreviations: CI, confidence interval; EIR, entomologic inoculation rate
![]() |
INTRODUCTION |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
However, malaria is a complex disease. Plasmodium falciparum infection confers only labile and partial immunity. Acquisition of such immunity only reduces the incidence of clinical malaria attacks without preventing infection (47
); an infected subject may acquire a new infection before recovering from a previous one (a phenomenon known as "superinfection") (8
); and immunity is acquired slowly and is a function of exposure to infecting mosquitoes (9
). To further complicate the problem, exposure to infectious mosquito bites is difficult to measure (10
). Practically, exposure cannot be assessed for each subject, and all individuals are usually considered identically exposed; this constitutes a potential confounding factor for acquisition of immunity. The presence of parasites in individuals gives little information on exposure, because parasitemia can last a long time in the absence of treatment (11
). Finally, in some geographic areas, exposure shows marked seasonality and can be very different from one year to the next (4
).
An epidemiologic model accounting for these complexities would be a useful tool with which to describe the dynamics of malaria and assess the epidemiologic status of exposed populations (12), answering the needs of both public health planners and malaria researchers. Among the relevant models reported in the literature (8
, 12
25
), few have been statistically calibrated (i.e., formally fitted to data) because of a lack of extensive longitudinal data and adequate statistical techniques. In two cases (16
, 21
), three or four of the model parameters were estimated through minimization of
2 or G2 (akin to entropy) criteria of fit. Simulation results from these models have never been presented with confidence intervals allowing an assessment of their reliability. Newly developed Bayesian numerical techniques (e.g., Markov chain Monte Carlo methods) (26
28
) offer the possibility of an extensive statistical treatment of complex epidemiologic models, including inference about model parameter values, calculation of confidence intervals for model predictions, model-checking, and hypothesis-testing. In this paper, we apply these techniques to the malaria model described by Struchiner et al. (23
). This model embodies a series of coherent, plausible research hypotheses about the natural history of P. falciparum malaria and offers a synthesis of previous work carried out by MacDonald (8
), Dietz et al. (16
), and Nedelman (21
). We use Markov chain Monte Carlo simulations to fit the model to longitudinal data gathered from the village of Ndiop in a meso-endemic area of Senegal. Model simulations of the underlying infection dynamics are presented.
![]() |
MATERIALS AND METHODS |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
Human data.
Only the 176 villagers (of a total of 396) who were continuously present in the village during the study period were included. Exclusions can be considered random with respect to exposure (in particular, it is very unlikely that people left the village because of mosquito bites or to obtain treatment elsewhere), and no bias should have been introduced by the removal of the traveling villagers. On the other hand, traveling villagers had unmeasured exposure to mosquitoes during their outings, and they could have introduced bias if included. A local medical care unit was created to support the study after agreement was reached with the population of the village and the public health authorities. Informed consent was obtained individually from the participants or from their parents (for children); approval was obtained from the Ministère du Plan et de la Coopération and from the Ministère de la Santé Publique. The local medical care unit was provided with basic equipment for malaria diagnosis. A team of four physicians, two technicians, and four fieldworkers stationed in the village was in charge of the contacts with the community and clinical follow-up. They were all trained in the clinical and laboratory diagnosis of malarial infection. Active surveillance consisted of: 1) a daily visit to each villager at home, to record body temperature and clinical symptoms of any type that had occurred during the previous 24 hours; and 2) collection of thick blood smears once per week from July to October of 1993 and once per month from November 1993 to August 1994. Passive surveillance included collection of thick blood smears from subjects reporting to the health unit on any occasion. Special attention was paid to the use of antimalarial drugs in the community, and the population was asked not to use any such drugs without prescription. The only antimalarial drug used was quinine (Quinimax; Sanofi-Labaz, Paris, France), administered at a dose of 25 mg/kg/day for 7 days to the following types of patients: children under 10 years of age with fever (temperature higher than 38°C) and high parasitemia (over 30 trophozoites per 100 leukocytes); pregnant women with clinical symptoms suggestive of malaria; persons with fever and very high parasitemia (over 200 trophozoites per 100 leukocytes); and persons with severe malaria symptoms (coma, etc.).
A total of 5,736 thick blood smears were collected from the 176 villagers studied. To study the natural evolution of parasitemia without interfering with antimalarial treatments, we excluded from the analysis thick blood smears collected within 15 days after the beginning of treatment. As a result, only 5,000 thick blood smears were considered in our analysis. All smears were double-read, once in the field by the technicians and once by an expert microscopist at the Institut de Recherche pour le Développement laboratory in Dakar, whose reading was definitive. Slides were stained using 4 percent Giemsa stain, and up to 200 microscopic oil-immersion fields were examined at a magnification of 100. Asexual-stage parasite (trophozoite) densities were reported as parasite count for 100 leukocytes (detection limit: 0.01 trophozoites for 100 leukocytes). The number, D(t), of trophozoite-positive thick blood smears on any given day for which subjects were seen, given the total number of thick blood smears examined, M(t), constituted a single data point. The entire data set is denoted D below. The observed prevalences presented in figure 2 were obtained by dividing D(t) by M(t) at each time t. Table 1 gives the distribution of subjects, positive thick blood smears, and total thick blood smears examined, by age and sex.
|
|
Dynamic model
The deterministic compartmental model developed by Struchiner et al. (23) was used to describe the natural history of the infection in humans. A brief summary is given here, since the model has been described elsewhere in full detail (23
, 32
, 33
). The model equations and the definition of each parameter are given in the Appendix.
In the model, the human population is divided into four epidemiologic classes or compartments (figure 1): nonimmune subjects and immune negative subjects, in proportions X1(t) and X3(t), respectively; and nonimmune subjects and immune positive subjects, in proportions Y2(t) and Y3(t), respectively. These proportions are time-varying, and the model can predict their full time course. Nonimmune positive subjects infectious to mosquitoes, in proportion Y1(t), are the subset of the nonimmune positives showing sexual-stage (gametocyte) parasitemia. It is assumed that immunity does not totally protect against infection but reduces the probability of becoming infected. Individuals are deemed positive if they show trophozoite parasitemia (as assessed by thick blood smear examination). For our study, the birth and death rate, , was equal to zero.
|
An effective inoculation by an infectious mosquito produces one brood of parasites within the human host. Different mosquito bites can each inject a brood into an individual, and those broods are cleared independently of each other. We did not model the mosquito part of the parasite cycle: The EIR defined in Struchiner et al.'s model (23) was replaced by our measurements, linearly interpolated.
Assuming all subjects, at t = 0, to be nonimmune and negative (X1(0) = 1) is not realistic for Ndiop. However, neither initial conditions nor the actual "initial time" are known a priori. Preliminary simulations showed that, when starting with realistic conditions, equilibrium is reached in approximately 2.5 years. In that period of time, the system has also practically lost memory of its initial state. Consequently, the initial time was chosen to be December 22, 1990. To acknowledge uncertainty in their values, the initial proportions of individuals in the compartments and initial average numbers of broods in individuals were considered as parameters to estimate. To respect the constraint of the summation to 1 of the initial proportions, we used the reparameterization given in the Appendix (equations 1719). No data on EIR were available from December 22, 1990, to July 1, 1993, before the beginning of parasitologic monitoring. To supply realistic weekly EIR values for input to the model (equations 9, 10, 14, and 16 in the Appendix) during that period, we used the average of the values recorded in Ndiop from July 1993 to December 1996 for each week.
The model differential equations are nonlinear and include delays. They were integrated numerically, using the "Lsodes" algorithm provided by MCSim software, version 4.2 (34). For given parameter values, integration of the equations shown in the Appendix gives the time courses of its variables (X1(t), X3(t), Y1(t), Y2(t), Y3(t), z1(t), z2(t), and z3(t)) over any period of time starting on December 22, 1990. The sum,
, of Y2(t) and Y3(t) is a model-computed estimate of the instantaneous prevalence of P. falciparum parasitemia in humans.
Statistical computations
A Bayesian approach was used to calibrate the model using the counts of P. falciparum trophozoite-positive thick blood smears in Ndiop (the data, D). Basically, each model parameter (see the list in table 2 and in the Appendix), i, was considered as a random variable and assigned an independent prior distribution, p(
i). Those distributions were updated together to yield a joint posterior distribution, p(
|D), such that model-computed time courses of prevalence (using parameter values drawn from that joint posterior) would be compatible with the data. According to Bayes' rule, p(|
D) is proportional to the product of the prior distributions by the data likelihood, p(D|
), under the model (35
, 36
).
|
To define the data likelihood, the observed number, D(t), of trophozoite-positive thick blood smears at time t was assumed to be binomially distributed with parameters , the model-predicted prevalence (0 <
< 1), and M(t), the total number of thick blood smears counted at t. The joint posterior is therefore of the form
![]() |
Unfortunately, because the dynamic model is nonlinear, there is no known analytical form for p(|D). It is impossible to describe it and report inference about the parameters in a direct way. It remains possible to summarize that distribution by drawing random sets of parameter values using Metropolis sampling (26
). This iterative procedure belongs to a class of Markov chain Monte Carlo techniques which has recently received much interest (27
, 37
39
). Briefly, the algorithm was as follows: At the start of a sampling "chain," all parameters were assigned values sampled from the priors. For any following iteration of the sampler, each component,
i, of the parameter vector was eventually updated by drawing a "proposed" new value,
i', out of a Gaussian "proposal" distribution centered on
i. Values of the joint posterior density at
i and
i' were then computed using the displayed equation (this required running the differential model to obtain all needed values for
). The two obtained density values were labeled
and
'. If
'/
exceeded 1, the new value
i' was accepted and replaced
i; otherwise,
i' was accepted only with probability
'/
. In the case of rejection of
i', the value
i was kept. After updating (eventually) all model parameters sequentially (the updating order does not matter in the long run), we recorded their values, therefore completing one iteration of the chain. Iterations were performed until the chain had reached equilibrium, i.e., until all parameters had approximately converged in distribution to p(
|D). The standard deviation of the Gaussian proposal distribution was adjusted periodically to yield an acceptance rate of 25 percent (40
). The convergence of several Markov chains to p(
|D) was assessed using Gelman and Rubin's
diagnostic (41
). The parameter sets recorded after equilibrium was reached were used to form histograms or compute summary statistics of the posterior distributions for estimands of interest (e.g., marginal parameter distributions, combinations of parameters, or model predictions). Obtaining the posterior distribution of model predictions required running the malaria model once for each parameter set recorded. All of the above computations were performed using MCSim software, version 4.2 (34
).
RESULTS
Model fit
Convergence of three independent Markov chain Monte Carlo chains was reached after approximately 40,000 iterations () diagnostic at 1.07 on average, ranging from 1 to 1.2). Fifteen thousand parameter sets were sampled, by keeping one out of every six iterations from an additional 30,000 of each chain (each iteration kept yielded a parameter set). All simulations and inferences presented below were made using this final sample from p(
|D).
A good fit to the data was obtained, while maintaining scientifically plausible parameter values. Figure 2 shows the daily trophozoite parasitemia prevalences (D(t)/M(t)) together with the corresponding model predictions, , made with the parameter set having the highest posterior density in the final sample. This model prediction is the "best" of all, but it is also quite representative of the set. Ten other model-predicted time courses of prevalence, obtained using parameter vectors randomly drawn from p(
|D), are presented. The differing model trajectories for all of these curves reflect uncertainty in model predictions. However, they all have similar behavior. The posterior 95 percent confidence interval for predictions is also displayed. The rise to the peak (from
20 percent prevalence to 7080 percent) is somewhat jagged and is mostly driven by the random biting rate of mosquitoes. The subsequent decrease is smoother and is driven by the gradual recovery of the infected subjects in the dry season. Prevalence returns to approximately 20 percent at the end of the dry season. Running a standard smoothing curve through the data would be purely descriptive and would provide no insight into the underlying dynamics. Our goal was not so much to "fit" the data as to extract from them information about the model parameters.
Posterior parameter distributions
The joint posterior distribution of all parameters can be viewed in several dimensions, but for simplicity only the marginal distribution of each parameter is described here. Table 3 summarizes these distributions on the basis of the final parameter sample. For all biologic parameters, the posterior location is noticeably different from the corresponding prior mean. Posterior standard deviations are much lower than specified a priori (compare tables 2 and 3), because important information about those parameters has been extracted from the data.
|
The recovery rate from infectiousness to mosquitoes among nonimmune positive hosts (1) is high but poorly identified. It corresponds to a median half-life of 6 days (95 percent CI: 0.5, 125). The window of infectivity is therefore quite small (as expected from the relative brevity of the prevalence peak during the year).
As indicated by the median value of 2, 213 days (1/0.0047) are necessary for a human host to acquire immunity to P. falciparum. This estimation is quite precise (coefficient of variation of
10 percent) and is much lower than the a priori value. Simulations indicate that for an inhabitant of Ndiop exposed to seasonal meso-endemic transmission, the numbers of infectious and noninfectious parasite broods present at any time in nonimmune hosts (z1 and z2) average approximately 10 (data not shown). Under such conditions (see equation 13 in the Appendix), the actual immunity acquisition delay (A1) is equal to 214 days on average, which is close to its minimum value. Immunity appears more quickly than expected a priori, but it still takes at least half a year to be in effect.
The interval of time, , until an immune host loses immunity in the absence of exposure to infectious mosquito bites is approximately 230 days (95 percent CI: 180, 290). The incubation period, N1, lasts 22 days, on average (95 percent CI: 15, 28).
Few infectious bites (f = 3 percent) are able to boost immunity, and this boosting does not seem to be very important a posteriori. The proportion of potentially infectious bites actually resulting in infection is much higher in nonimmune hosts (b1 = 40 percent; 95 percent CI: 26, 91) than in immune subjects (b2 = 2.5 percent; 95 percent CI: 1.5, 3.5). Immunity, although it is progressively acquired, seems to efficiently protect against infection.
The marginal posterior distributions of the sampled initial state variables or their reparameterizations (X3(0), FY2(0), FY3(0), z1(0), z2(0), and z3(0)) are very close to the corresponding priors. This shows the insensitivity of the model to those parameters: Their values do not appreciably affect the results.
Model predictions of the underlying dynamics
After fitting, the model can be used to simulate various scenarios to better understand the dynamics of P. falciparum infection in our study population. Figure 3 shows predictions of the time course of malaria immunity status during the year of our study and the following wet season (July 1993 to December 1994), together with the measured EIR. At the end of the dry season, the population composition is as follows: 63 percent (95 percent CI: 35, 85) nonimmune negative (X1), 12 percent (95 percent CI: 10, 15) nonimmune positive (Y2), 1 percent (95 percent CI: 0.7, 2) immune positive (Y3), and 24 percent (95 percent CI: 5, 50) immune negative (X3). These proportions vary from year to year as mosquito biting fluctuates in timing and intensity. Nonetheless, some behaviors appear stable. The proportion of nonimmune negative individuals falls very quickly, and practically to zero, as soon as mosquito biting increases. Infected subjects initially transfer to a nonimmune positive status, whose proportion reaches a peak at approximately the same time as the biting rate (with a short delay imposed by incubation). Nonimmune positive subjects then transfer mostly to the immune status. However, the fraction of immune positive individuals does not increase much, and individuals quickly eliminate parasites to become immune negative. Near the end of the dry season, immune negative individuals lose immunity and the fraction of nonimmune negative persons increases quickly. Overall, most subjects move through the four states as the year progresses.
|
The incidence of malarial infection is difficult to measure, since it requires identification of new infections in individuals who are potentially already parasitemic. The model can easily provide an estimate for various incidence rates. Figure 4 presents the model-reconstructed instantaneous incidence rates (number of new cases per day) of trophozoite parasitemia contributed by nonimmune and immune individuals during the 1993 wet season in Ndiop (incidence was null during the dry season). These rates correspond to the products 1(t) x X1 x 396 and
2(t) x X3 x 396, respectively (396 being the total size of the population of Ndiop). The time-weighted average incidences, over the year, contributed by nonimmune and immune individuals are approximately the same0.52 (95 percent CI: 0.40, 0.63) and 0.47 (95 percent CI: 0.30, 0.67) cases per day, respectively. However, the incidence time profiles for these two subpopulations do differ. Early incident cases are contributed mostly by nonimmune individuals and late cases by immune subjects. This is explained by the progressive decline of the nonimmune negative population as the wet season progresses (figure 3). For the whole population, the time-weighted average is estimated at 0.99 (95 percent CI: 0.87, 1.2) cases per day.
|
|
The data were obtained through intensive follow-up. A potential bias in subject recruitment arose as individuals were offered the opportunity to present themselves for clinical examination. P. falciparum-infected individuals with presenting symptoms may have been overrepresented. However, blood samples were analyzed at all self-motivated visits, regardless of whether the visits were related to malaria, and only 16 percent of the samples analyzed were obtained during such visits. We also verified, through examination of residuals after model calibration, that the prevalence of P. falciparum parasitemia in self-motivated consultations was not higher than that in systematic screenings. The impact of a potential bias in self-motivated consultations should therefore have been small or nonexistent. We only analyzed data on the prevalence of P. falciparum trophozoite parasitemia, but it would be interesting to extend the model to also consider data on the number of clinical malaria attacks.
The model developed by Struchiner et al. (23) offers a reasonable, albeit simplified, description of malaria physiopathology. We did not include the original description by Struchiner et al. of the parasite cycle in mosquitoes. That was not needed, since our data included EIR throughout the year. However, that rate was assumed to be precisely measured and identical for all subjects. This assumption was needed because the full treatment of "error in variables" problems is difficult in the context of large and computationally intensive models. Another important set of modeling assumptions concerns immunity. In the model, an infected person does not necessarily acquire immunity after one inoculation, and immunity can be lost with time. Although the hypothesis of a definitive acquisition of immunity to the different antigenic strains of parasite seen during an individual's lifetime (18
) is not explicitly considered, the model does assume that the acquisition of immunity is a function of the number of coinfecting strains.
Our analysis was performed by pooling data on all subjects (but still preserving the longitudinal aspect of the data at the population level). This could be improved by taking into account the age structure of the populationfor example, through a hierarchical statistical model (44). This could shed light on age-related differences in susceptibility to P. falciparum infection. At this occasion, it might be possible to take into account the fact that the feeding behavior of mosquitoes is affected by a number of host- or environment-related factors (45
48
).
According to the model, under conditions similar to those in Ndiop, the fraction of susceptible subjects is the highest at the very end of the dry season, when mosquitoes start biting again. This makes sense given what is known of the natural history of malaria. The advantage of using a calibrated model, assuming it is correct or sufficiently robust, is that it offers a quantitative estimate of this fraction and of the associated uncertainty. Use of such information in vaccination trial design can help researchers assess and improve statistical power. Power calculations show that the effective size of a trial is proportional to the fraction of susceptible subjects (e.g., a study with 10,000 person-days and 50 percent susceptible subjects in each group has the same power as a study with 5,000 person-days and 100 percent susceptible subjects). Location-specific EIRs could be used for input to the model to assess the best time of the year for a study in areas other than Ndiop.
A dynamic perspective on malaria, as embodied in an epidemiologic model able to disentangle time-varying exposures, superinfections, and complex immunity acquisition processes, is essential for a proper analysis of malaria field study data. Too many pitfalls of confounding and bias, difficult to avoid, await standard data analyses. The model analyzed here is by no means complete or perfect, but it offers a reasonable basis for extension and improvement. Several research teams worldwide are currently attempting to improve malaria models. These efforts would benefit from the statistical techniques presented here. Calibrated models can be powerful predictive tools for experimental design and exploration of public health measures.
![]() |
APPENDIX |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
![]() | (1) |
![]() | (2) |
![]() | (3) |
![]() | (4) |
![]() | (5) |
![]() | (6) |
![]() | (7) |
![]() | (8) |
![]() | (9) |
![]() | (10) |
![]() | (11) |
![]() | (12) |
![]() | (13) |
![]() | (14) |
![]() |
![]() | (15) |
![]() | (16) |
![]() | (17) |
![]() | (18) |
![]() | (19) |
The symbols used in the above equations and in the text are defined below (in alphabetical order).
A1: rate at which immunity to Plasmodium falciparum infection is acquired by a human host.
b1: proportion of bites by infectious mosquitoes on negative nonimmune hosts actually resulting in infection.
b2: proportion of bites by infectious mosquitoes on negative immune hosts actually resulting in infection.
f: boosting factor (i.e., proportion of bites by infectious mosquitoes on immune hosts resulting in boosted immunity).
FY2(0): deconstraining parameter for Y2(0).
FY3(0): deconstraining parameter for Y3(0).
he(·): entomologic inoculation rate (EIR): number of P. falciparum infectious bites per human per day.
N1: P. falciparum parasitemia incubation period in humans (in days).
r1: rate constant of elimination of a brood of parasites by nonimmune positive hosts (in days-1).
r2: rate constant of elimination of a brood of parasites by immune positive hosts (in days-1).
R1: recovery rate for nonimmune positive individuals (in days-1).
R2: recovery rate for immune positive individuals (in days-1).
X1: proportion of nonimmune negative (i.e., naive) individuals in the population.
X1(0): initial value (at time 0, i.e., December 22, 1990) of X1.
X3: proportion of immune negative individuals.
X3(0): initial value of X3.
Y1: proportion of nonimmune positive individuals potentially infectious for mosquitoes.
Y2: proportion of nonimmune positive individuals.
Y2(0): initial value of Y2.
Y3: proportion of immune positive individuals.
Y3(0): initial value of Y3.
z1: average number of infectious broods of the parasite per nonimmune positive human host.
z1(0): initial value of z1.
z2: average number of noninfectious broods of the parasite per nonimmune positive human host.
z2z(0): initial value of z2.
z3: average number of noninfectious broods of the parasite per immune positive human host.
z3(0): initial value of z3.
1: recovery rate from infectiousness to mosquitoes among nonimmune positive hosts (in days-1).
2: maximum rate at which immunity to P. falciparum infection can be acquired by a human host (in days-1).
: death and birth rate in the human population (in days-1).
3: daily fraction of immune negative subjects losing immunity.
1: infection rate for nonimmune negative subjects (probability per day of such a subjectÕs becoming infected).
2: infection rate for immune negative subjects (probability per day of such a subject's becoming infected).
: time delay needed for an immune host to lose immunity in the absence of exposure to infection (in days).
![]() |
ACKNOWLEDGMENTS |
---|
The authors acknowledge the administrative and technical support of the Institut Pasteur and Institut de Recherche pour le Développement in Dakar, Senegal, particularly the help of E.-H. Ba. The authors are grateful to the villagers of Ndiop and to the local study team for their collaboration. They also thank Drs. K. Dietz, J. Eisenberg, D. Fontenille, O. Garraud, S. Gupta, and L. Lochouarn for their helpful comments.
![]() |
NOTES |
---|
![]() |
REFERENCES |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|