Causes of Variation in Malaria Infection Dynamics: Insights from Theory and Data
Abstract
Parasite strategies for exploiting host resources are key determinants of disease severity (i.e., virulence) and infectiousness (i.e., transmission between hosts). By iterating the development of theory and empirical tests, we investigated whether variation in parasite traits across two genetically distinct clones of the rodent malaria parasite, Plasmodium chabaudi, explains differences in within-host infection dynamics and virulence. First, we experimentally tested key predictions of our earlier modeling work. As predicted, the more virulent genotype produced more progeny parasites per infected cell (burst size), but in contrast to predictions, invasion rates of red blood cells (RBCs) did not differ between the genotypes studied. Second, we further developed theory by confronting our earlier model with these new data, testing a new set of models that incorporate more biological realism, and developing novel theoretical tools for identifying differences between parasite genotypes. Overall, we found robust evidence that differences in burst sizes contribute to variation in dynamics and that differential interactions between parasites and host immune responses also play a role. In contrast to previous work, our model predicts that RBC age structure is not important for explaining dynamics. Integrating theory and empirical tests is a potentially powerful way of progressing understanding of disease biology.
Introduction
In both natural and laboratory host-parasite systems there is widespread variation in within-host infection dynamics (i.e., patterns of parasite densities, availability of host resources, and immune responses over the course of an infection) and disease outcomes (e.g., Holmes and Burch 2000; Paterson and Viney 2003; Reece et al. 2009; Morrison et al. 2010), yet the ecological and evolutionary causes and consequences of this variation are poorly understood. Since within-host interactions between parasite genotypes and between hosts and parasite genotypes (G × G interactions) are major drivers of disease evolution, a better understanding of variation within parasite species is important. A defining feature of parasites is that they exploit host resources for their own survival and transmission, with obvious implications for host health, so parasite traits involved in host exploitation could underlie observed variation in virulence (i.e., reduction in host fitness as a result of infection). Disentangling the roles of parasite traits, host traits, and their interactions is a major challenge in the evolutionary biology of infectious diseases with important implications for disease severity, patterns of transmission and control (e.g., Mackinnon et al. 2000; Nesse and Stearns 2008; Read and Huijben 2009; Pollitt et al. 2011b).
Rodent malaria parasites are a useful system with which to investigate the existence and implications of genetic variation in parasite traits. Experiments with a number of genetically distinct clones of the rodent malaria Plasmodium chabaudi reveal substantially different dynamics as well as different morbidity and mortality outcomes (e.g., Mackinnon and Read 1999; Bell et al. 2006; Grech et al. 2006; Pollitt et al. 2011a). Malaria parasites exploit host red blood cells (RBCs) during the disease-causing stage of their life cycle. After invading a host RBC, an asexual blood stage malaria parasite replicates, eventually bursting the RBC and releasing multiple progeny parasites (merozoites), each with the capacity to invade another RBC and replicate. Like malaria parasites of humans, different species of rodent malaria parasites are known to differ in the number of merozoites produced per infected cell (burst size) and in the ages of RBCs they infect (Garnham 1966; Paul et al. 2003; Reece et al. 2009). Differences in these “host exploitation” traits are predicted to underlie variation in infection dynamics and virulence across genotypes of P. chabaudi (Antia et al. 2008; Mideo et al. 2008b). Experiments suggest that dynamics of infections with multiple genotypes may be partially resource driven (Råberg et al. 2006; Barclay et al. 2008; Pollitt et al. 2011a); therefore, differences in host exploitation traits could also explain the outcome of competitive interactions between parasite genotypes (e.g., de Roode et al. 2005; Bell et al. 2006; Barclay et al. 2008; Pollitt et al. 2011a). Developing an understanding of the consequences of variation in host exploitation traits is extremely challenging, in part due to the dynamical nature of infections: multiple causal and confounding factors are likely to be varying simultaneously. A powerful way of generating this understanding is by combining theoretical approaches that aim to tease apart mechanisms by formalizing biological hypotheses with experimental approaches that unambiguously test these hypotheses (e.g., Haydon et al. 2003; Antia et al. 2008; Mideo et al. 2008a; Kochin et al. 2010).
As a laboratory disease model, P. chabaudi offers an excellent opportunity for iterating model building and experimentation and progress has been made on the first stages of this process (i.e., fitting models to data; Haydon et al 2003; Antia et al. 2008; Mideo et al. 2008b; Kochin et al. 2010; Miller et al. 2010). We previously developed mathematical models of the within-host dynamics of P. chabaudi infections and fitted these models, using a maximum likelihood method, to data from hosts that were experimentally depleted of their CD4+ T cells (Mideo et al. 2008b). CD4+ T cells play a major role in fighting malaria infections (Langhorne et al. 1990; Good and Doolan 1999; Pombo et al. 2002; Urban et al. 2005; Stephens and Langhorne 2006), so depleting this immune component removes some of the complexity of the system, reducing the number of factors that can regulate infection dynamics. Our specific aim in fitting models to these data was to determine whether red blood cell limitation was sufficient for explaining the acute phase of infection dynamics in this CD4+-depleted setting and whether interactions between parasites and host resources differed across P. chabaudi genotypes (Mideo et al. 2008b). The genotypes we studied (AS and DK) differ in their relative virulence, as defined by weight loss; AS is more virulent in single infections and capable of competitively suppressing the less virulent DK in mixed infections (Barclay et al. 2008). Our results (Mideo et al. 2008b) suggested that early within-host infection dynamics could be explained by resource limitation alone and that there was important structuring of the resource pool by age, that is, parasite interactions with RBCs differed across ages of RBC (see also Cromer et al. 2006, 2009; Antia et al. 2008).
We predicted that the different patterns of infection dynamics observed across the genotypes could be explained by parasite variation in traits involved in the exploitation of these resources. In particular, the more virulent genotype AS was predicted to have higher burst sizes in all ages of RBCs and produce more progeny parasites when infecting young RBCs (reticulocytes) than when infecting fully mature RBCs, and vice versa for the less virulent genotype DK. Further, while both genotypes were predicted to have a “preference” for invading mature RBCs—with invasion rates of these cells predicted to be an order of magnitude higher than for reticulocytes—the more virulent genotype was predicted to invade mature RBCs at a faster rate than the less virulent genotype. Using an independent data set for validation, we also showed that these differences between genotypes and their interactions with RBCs of different ages were sufficient to generate patterns that qualitatively match with the early dynamics of mixed genotype infections (Mideo et al. 2008b).
Here, we test the predictions of our previous modeling work by repeating the experimental procedures used to generate the original data modeled (i.e., single-genotype infections in CD4+ T cell–depleted mice; Barclay et al. 2008), this time using novel empirical methods to measure invasion rates and burst sizes. We find that some of our original predictions are supported while others are not. We therefore investigate if the model of Mideo et al. (2008b) can generate the within-host infection dynamics observed in the new data set, given our empirical estimates of some of the model parameters, or if the dynamics are better captured by a model with a more complex mathematical description of the underlying biology.
Methods
Experiment: Estimating Invasion Rates and Burst Sizes
Our original model (Mideo et al. 2008b) was fitted to data from Barclay et al. (2008). We repeated the same experimental procedures as in that study, including experimental depletion of CD4+ T cells. From these new infections, we estimated RBC invasion rates and generated a new data set to inform the next iteration of models by tracking the dynamics of RBC and parasite densities over the course of individual infections. Although our model (Mideo et al. 2008b) generated specific predictions for parasites in different ages of RBC, it is not possible to directly measure burst sizes for infected RBCs of different known ages. To overcome this, we added additional treatment groups to the experiment in which we used phenylhydrazine (PHZ)—which chemically lyses RBCs, thus inducing the production of new RBCs—to manipulate host RBC age structure before infection. This enabled us to compare average burst sizes in hosts with the increased proportions of reticulocytes induced by PHZ treatment with average burst sizes in PHZ-control mice in which reticulocytes are relatively scarce.
Parasites and Hosts
We initiated experimental infections with one of two genetically distinct Plasmodium chabaudi genotypes denoted “AS” and “DK” (World Health Organization Registry of Standard Malaria Parasites, University of Edinburgh), originally isolated from thicket rats Thamnomys rutilans, in the Central African Republic (Beale et al. 1978). All mice were inbred female C57BL/6JolaHsd mice aged 6–8 weeks old (Harlan-Olac) and received an intraperitoneal inoculation of 106 parasitized red blood cells in a 0.1-mL dose consisting of 47.5% Ringer’s (27 mM KCl, 27 mM CaCl2, 0.15 M NaCl), 50% heat-inactivated calf serum, and 2.5% heparin (200 U/mL). Mice were housed randomly in groups of 5 at 20°C with a 12L∶12D photperiod and provided food (41B, Harlan-Teklad) and water with 0.05% PABA (to enhance parasite growth) ad lib.
Manipulating Immunity
As described in Barclay et al. (2008), mice received intraperitoneal injections of an anti-CD4+ antibody (GK1.5) dissolved in phosphate-buffered saline (PBS), according to the following regimen: 500 μg 5 days before infections, 250 μg 4 days before infections, 250 μg 1 day before infection, and 250 μg weekly after infection. To confirm the effects of the depleting treatment, we also set up immune-intact control infections where mice (5 replicates for each parasite genotype) received injections of a nondepleting antibody of the same isotype (IgG I4131, Sigma) following the same regimen. Using a fluorescence-activated cell sorter (as desaccribed in Barclay et al. 2008), we verified that CD4+ depletion was successful. The proportion of circulating CD4+ T cells was reduced by 2 orders of magnitude in the depleted versus immune-intact control infections, and as expected, parasite densities were higher in depleted mice (see fig. A1 ).
Experimental Setup and Sampling
We ran a fully cross-factored experiment with two treatments in which we manipulated host RBC age structure and parasite genotype. A total of 20 mice were given CD4+ T cell–depleting antibodies, and groups of 5 of these mice were treated with one of the following regimes: (a) 106 AS parasites and PBS (PHZ control, henceforth referred to as “control”), (b) 106 DK parasites and PBS, (c) 106 AS parasites and PBS+PHZ, (d) 106 DK parasites and PBS+PHZ.
Treatment groups a and b represent a repeat of the experiment in Barclay et al. (2008). Infections from these groups were used to estimate invasion rates for comparison with model predictions (Mideo et al. 2008b), and to generate a new data set with which to test and develop models. Treatments c and d were necessary for testing predictions about RBC-age-specific burst sizes. We took daily measurements of individual mass, and blood samples from the tail from all individual hosts until 18 days postinfection to prepare Giemsa-stained thin blood smears, to estimate RBC density (using flow cytometry; Beckman Coulter), and to measure parasite density with genotype-specific real-time quantitative PCR (as described in Barclay et al. 2008). Additional measurements from infections in treatment groups a and b were required for estimating invasion rates. On the relevant days postinfection, we used thin blood smears under ×1,000 magnification to calculate (i) the proportion of RBCs that were infected by estimating the total number of RBCs in at least five microscope fields, counting at minimum 100 RBCs, and (ii) the proportions of infected and uninfected reticulocytes and mature RBCs by counting at least 40 of each type of cell, or up to 20 microscopic fields. All measurements were converted to densities by multiplying proportions by the measured RBC densities. One CD4+-depleted mouse, infected with AS parasites, had very unusual dynamics, with zero parasite density until day 8 postinfection. We exclude this mouse from all analyses and model fitting.
Measuring RBC Invasion Rate
Expressions for the invasion rates of reticulocytes, βR, and mature RBCs (normocytes), βN, can be obtained by rearranging equations (A14) and (A15) of Mideo et al. (2008b) to give
and 
are the densities of infected reticulocytes and mature RBCs on day i, respectively. A hat indicates a density immediately before the invasion process begins and the absence of a hat indicates a density after the invasion process is finished, when all merozoites are assumed to have invaded an RBC or died. Our experimental procedures provided measures of 
, 
, 
, 
, and 
. We estimate invasion rates on day 6 postinfection, that is, 
, from each individual infection, to be consistent with our measurements of burst sizes (see next section).
To estimate the number of merozoites that are present on day 6 to start the invasion cycle, 
, we multiply the density of infected RBCs from day 5 by the average burst size calculated on day 5. Plasmodium chabaudi has a 24-h cycle, meaning that merozoites released from RBCs that were infected on day 5 give rise to the infected RBCs observed on day 6. We assume that no RBC death takes place before or during a given round of parasite invasion, so that the total density of RBCs we observe in the blood smears, after invasion is complete, is equal to the density of RBCs that the parasites were exposed to at the start of the invasion cycle. Since we are sampling about 8 h into the parasites’ synchronous 24-h cycle (O’Donnell et al. 2011), there is no RBC loss due to parasites rupturing cells. Thus, 
and 
are the densities of infected reticulocytes and mature RBCs on day 6, respectively; 
is estimated as the sum of the densities of uninfected and infected reticulocytes on day 6 and, similarly, 
the sum of the densities of uninfected and infected mature RBCs on day 6. Finally, since the exact value of the merozoite life span, 1/μ, is not known, we rearrange equations 1 and 2 and estimate the invasion rate multiplied by the merozoite life span, that is, βR/μ and βN/μ (as in Hetzel and Anderson 1996). Any conclusions about differences in invasion rates thus implicitly assume that merozoite life span is equal across parasite genotypes.
Measuring Burst Size
It is not possible to sample mature parasites (schizonts) in vivo, because most asexual parasites leave the circulation in late stages of development to sequester in host tissues and reach maturity. Therefore to measure burst size without destructive sampling we collected parasites from tail snips at trophozoite stage, before sequestration, and cultured them in vitro to maturity when the number of progeny merozoites inside individual schizonts can be counted (Reece and Thompson 2008). Specifically, we cultured 10 μL parasitized blood from each mouse in 1 mL complete culture medium (RPMI medium with NaHCO3, Hepes, and l-glutamine (Invitrogen), containing 25% heat-inactivated fetal calf serum (Gibco), at pH 7.25) in the presence of 10% O2, 5% CO2, 85% N2, at 30°C. We sampled cultures at 18.5, 20, and 21 h to enable us to capture the optimal time for assaying burst size (i.e., to sample when the frequency of mature schizonts is high but before they rupture and release merozoites that rapidly reinvade). We cultured parasites from all infections on both days 5 and 6 postinfection to ensure we captured the optimal percent parasitemia for schizont culture (<10%). Cultures from day 5 were more successful than day 6, when several treatment groups exceeded 10% parasitemia, and 20 h proved to be best for capturing the highest proportion of mature schizonts.
Since schizonts fill entire RBCs and obscure staining patterns, it is not possible to reliably determine whether a schizont is in a reticulocyte or a mature RBC. To overcome this, we manipulated the RBC age-structure of a subset of hosts with PHZ. Treated mice received an intraperitoneal injection of 40 mg/kg of phenylhydrazine hydrochloride dissolved in PBS, 2 days before infection. This dose of PHZ causes a rapid decrease in RBC densities and a subsequent influx of reticulocytes as the system “re-equilibrates” (Savill et al. 2009). Control mice received equal-volume injections of PBS alone. RBC densities in our experiment show a marked decrease in response to PHZ treatment (see fig. A1). We assessed the effects of PHZ on RBC age structure with nonparametric Wilcoxon rank sum tests. Mice treated with PHZ had significantly higher proportions of bloodstream reticulocytes than control mice over the relevant days postinfection (days 5 to 7: 
in all cases; 
). Further, the proportion of infected RBCs that were reticulocytes was also significantly higher in PHZ-treated mice across these days (
in all cases; 
). These effects are plotted in figure A2. Average burst sizes for individual mice were estimated by counting the number of merozoites in at least 25 mature schizonts from thin blood smears of cultured blood. Burst sizes for one PHZ-treated mouse that was infected with DK could not be estimated since blood smears contained no schizonts, and an absence of parasites in the blood was confirmed by PCR.
Statistical Analyses
All statistical analyses were carried out using R, version 2.7.1 (R Foundation for Statistical Computing,
Mathematical Models
In addition to empirically testing the predictions of Mideo et al. (2008b), our experiment provides a new data set (RBC and parasite dynamics from replicate AS and DK infections, in CD4+ T cell–depleted mice) for challenging our previous model and informing the next iteration of model building. We use this data set in two ways: first, we test the qualitative fit of the Mideo et al. (2008b) model, given our experimental estimates of some of the model parameters; second, we test whether that model provides a quantitatively good fit in comparison to more complex models.
Confronting the Model of Mideo et al. (2008b) with New Data
The model of within-host malaria infection dynamics we developed in Mideo et al. (2008b) tracks the changes in densities of parasites and red blood cells of different ages. The model does not include host immune responses, and therefore infection dynamics are solely regulated by resource (i.e., RBC) availability. In turn, resource availability is controlled by host traits, such as the rate of RBC replenishment, and parasite traits, such as burst sizes and invasion rates. Using our experimental estimates of genotype-specific burst sizes and invasion rates we asked, could the model of Mideo et al. (2008b) generate dynamics that qualitatively match the average dynamics observed in our new experiment? All other model parameters, for example, those governing RBC production, were set to the published genotype-specific best estimates (table C3 in Mideo et al. 2008b), and we chose initial parasite and RBC densities so as to best capture early dynamics. Comparing the dynamics generated from this model with the average dynamics for each genotype is only a rough assessment of the model’s suitability—we expected the qualitative patterns of infection to be robust but the quantitative match to be imperfect since we ignored individual variation in infections which is expected to be important (Molineaux and Dietz 1999; Mideo et al. 2008b; Miller et al. 2010).
The Next Iteration of Model Development and Comparisons
Given a poor qualitative fit of the model of Mideo et al. (2008b) to the average dynamics in our new experiment (see below), we next asked whether that model could provide a quantitatively good fit to individual data sets, when allowing for variation in model parameters across individual infections, or if we needed a more complex description of the underlying biology of the system. A main component of P. chabaudi infections that is not included in the model of Mideo et al. (2008b) is host immunity (the focus of other models, e.g., Haydon et al. 2003; Kochin et al. 2010; Miller et al. 2010). Recent work has shown that models require the inclusion of up to three different immune responses—independently targeting merozoites, infected RBCs, and uninfected RBCs—to capture the dynamics of malaria infections in mice with intact immune systems (Miller et al. 2010). While CD4+ depletion has a clear influence on dynamics (Barclay et al. 2008), there could be an important role for innate immunity (Stevenson and Riley 2004) and other CD4+ T-cell independent mechanisms in regulating parasite densities in our experimental infections. For example, parasite densities may be controlled early in infections by macrophages through phagocytosis of infected RBCs and merozoites (e.g., Su et al. 2002) and platelets through direct killing of parasites within infected RBCs (McMorran et al. 2009). We therefore incorporated the immune responses of the model of Miller et al. (2010) into our original RBC age-structured model (Mideo et al. 2008b) and refer to this as the “hybrid” model. A brief derivation of this model is given in appendix A. Determining how best to describe the multifactorial complexity of host immune responses is challenging, because there are very few time-series data on the effector cells and molecules to fit the models to, there is a poor quantitative understanding of what activates these cells, at what rates, and how their densities translate to actual parasite killing. Instead of modeling the mechanistic details of the immune system per se, we model immune responses phenomenologically as clearance rate functions over the course of infection (as in Miller et al. 2010). Each function is zero, except during “windows” of immune activity defined by four parameters: day postinfection when immune clearance begins, maximum clearance rate, rise time to maximum rate, and total duration of immune activity window (see fig. 1 for a schematic example). In this way, we allow empirical data to determine activation time of immune responses rather than having to impose arbitrary assumptions on the mechanisms of activation. Our experimental data shows two clear peaks of parasites that could represent different antigenic variants. Since each of these could potentially have a unique, specific immune response, we allow for two windows of immune activity.
Schematic of immune clearance rates of parasitized red blood cells (RBCs). Fitted model parameters are labeled for one of the two “windows” of immune activity. These include day postinfection when immune clearance begins (sp), maximum clearance rate (cp), rise time to maximum rate (rp), and total duration of immune activity window (lp). We allow for two windows of immune activity since experimental infection data show two clear peaks of parasites that could each potentially have a unique, specific immune response. For further details, see appendix A.
We compared the fit of this hybrid model with that of the Mideo et al. (2008b) model that excludes immunity. We also study reduced versions of this hybrid model by systematically removing regulatory components, thus independently testing the role of RBC age structure and each of the three immune responses. All models were fit to our newly collected parasite and RBC density dynamics from individual mice using the adaptive, population-based Markov chain Monte Carlo (MCMC) method developed in Savill et al. (2009) and Miller et al. (2010). The MCMC method rarely becomes trapped in local maxima (common in other model-fitting methods), is excellent at sampling from multiple modes of the posterior distribution that can arise in nonlinear dynamical systems and automatically determines the full variance-covariance structure of the parameters. This theoretical approach also allows us to incorporate prior beliefs about model parameters from our previous model-fitting studies and naturally updates these beliefs with the new empirical estimates presented here. In addition, we estimate the genotype-level (hyper)parameters, allowing us to infer genotype-level effects in the parameters (see app. A). Rather than comparing single estimates from fitting to average dynamics (Antia et al. 2008; Kochin et al. 2010), or means or medians of estimated parameter distributions from individually fitted data sets (Mideo et al. 2008b; Miller et al. 2010), this is a more conservative approach for determining whether the model predicts differences in parasite traits (burst sizes and invasion rates) between genotypes or host traits (RBC production) between infections with different parasite genotypes. It also bypasses the need for an arbitrary description of significant differences (Mideo et al. 2008b). Models were compared using Bayes factors, as described in detail in the supplementary material of Miller et al. (2010).
Results
Experiment
Estimating Invasion Rates
Analysis of invasion rates shows a significant interaction between genotype and RBC age (
, 
; see table B1; fig. 2A), with the more virulent genotype AS having higher invasion rates of mature RBCs than reticulocytes and vice versa for the less virulent genotype DK. But, the lack of clear variation in invasion rates across ages of RBCs and genotypes contradicts our predictions that invasion rates of mature RBCs would be higher for the more virulent genotype and that invasion rates of reticulocytes would be lower than for mature RBCs for both genotypes (fig. 2A). Most striking is this last discrepancy; while all measured rates were higher than predicted, reticulocyte invasion rates were higher by an order of magnitude.
Predictions from the model of Mideo et al. (2008b; left) and experimental data (right) for invasion rates (A) and burst sizes (B). A, left, The invasion rate of reticulocytes (Retic) is predicted to be similar for both parasite genotypes and much lower than for mature red blood cells (RBCs; Normo). The more virulent genotype (AS, triangles) is predicted to have a higher invasion rate of mature RBCs than the less virulent genotype (DK, circles). A, right, We find little support for our predictions. While invasion rates of both genotypes are similar in reticulocytes, they are an order of magnitude higher than expected, and higher invasion rates in mature RBCs across genotypes were not observed. Means of log10-transformed invasion rate estimates × merozoite life span (assumed to be 30 min) are presented from four or five individual infections. Values along the Y-axis have units of cells per microliter. B, left, Burst sizes are predicted to be higher for AS (more virulent) than for DK (less virulent). The model predicts that AS has a higher burst size in reticulocytes than in mature RBCs, so the average burst size in AS infections should increase when the proportion of infected RBCs that are reticulocytes is experimentally increased in a subset of hosts, while the opposite is predicted for DK. B, right, As predicted, AS had a higher burst size in control infections, but increasing the proportion of reticulocytes did not result in the predicted changes. Means are calculated from average burst sizes in four or five individual infections from each treatment. Error bars = ±1 SEM.
Estimating Burst Sizes
Analysis of burst sizes across both control and PHZ-treated hosts shows that, in support of one of the predictions from Mideo et al. (2008b), burst sizes are significantly higher for the more virulent genotype AS (
, 
; fig. 2B). A second prediction was that the interaction between genotype and PHZ would be significant: we predicted that PHZ would result in an increase of the burst size of the more virulent genotype AS and a decrease in the burst size of the less virulent genotype DK. However, our results show that the interaction between genotype and PHZ was not significant (
, 
; fig. 2B). Instead, a marginally significant main effect of PHZ suggests that across genotypes, PHZ treatment resulted in an increase in average burst sizes (
, 
). While the effect of PHZ is marginal, we include it in the minimal model since it is a main effect and the minimal model explains 10% more variation than the model that excludes PHZ (adjusted 
vs. 0.36). Averages for each treatment are given in table 1, third column. All of the covariates fitted in the maximal model were nonsignificant (see table B2).
Estimating the RBC age-specific burst sizes for individual genotypes
| Genotype, treatment | Average burst size | Relative proportion of infected reticulocytes | Estimated reticulocyte burst size | Estimated mature RBC burst size |
|---|---|---|---|---|
| More virulent, AS | 36.79 | 5.38 | ||
| Control | 5.72 | .011 | ||
| PHZ | 6.69 | .042 | ||
| Less virulent, DK | 19.39 | 4.81 | ||
| Control | 5.04 | .015 | ||
| PHZ | 5.32 | .035 |
Our predictions for the effects of PHZ were based on two assumptions: first, as predicted by the model of Mideo et al. (2008b), burst sizes differ for reticulocytes and mature RBCs, and second, PHZ increases the relative abundance of circulating reticulocytes and the relative proportion of infected RBCs that are reticulocytes. Therefore, average burst sizes should change in a predictable way under PHZ treatment. Note that the predictions presented in fig. 2B for the parasite genotype-specific effects of PHZ on burst size are qualitative only. Useful quantitative predictions were precluded since we did not know a priori what the relative proportions of infected reticulocytes and infected mature RBCs would be, nor by how much PHZ would change this. (While there are good estimates for the effects of PHZ in mice [Savill et al. 2009], there are no relevant measurements of the combined effects of PHZ and malaria infection for our host-parasite system.) Given that we observed an increase in burst sizes in response to PHZ treatment, we can a posteriori ask whether our original assumptions were valid, that is, whether this result could be driven entirely by the PHZ-induced change in the proportion of infected reticulocytes and a higher burst size in these younger RBCs. Assuming this was the case, we can estimate what the RBC age-specific burst sizes would have to be if it was the sole cause of the change in the average burst sizes observed across the PHZ and control treatments (table 1, third column), given the associated change in the average proportion of all infected RBCs that were reticulocytes (fig. A2). These estimates of RBC-age specific burst sizes for each parasite genotype (∼37 for AS and 20 for DK; table 1) are higher than the parameter estimates from the Mideo et al. (2008b) model and given that the maximum burst size observed in our experiment was 12 (see fig. B1 for histograms of burst size measurements), it is extremely unlikely that the change in burst size observed in the PHZ treatment can be attributed solely to higher burst sizes in reticulocytes. More plausible is that there is plasticity in burst sizes and that PHZ, for some reason, results in higher burst sizes in both RBC age classes for both genotypes.
Mathematical Models
Confronting the Model of Mideo et al. (2008b) with New Data
Given that some of our experimental measurements do not support the predictions from our previous model (Mideo et al. 2008b), we first ask whether that model can still capture infection dynamics when we set the appropriate model parameters to their experimentally measured values. We find that, regardless of how we deal with the burst sizes (i.e., assume reticulocyte and mature RBC burst sizes are the same and equal to the measured average, or use the estimated age-specific burst sizes from table 1), the dynamics do not qualitatively match the data (cf. fig. 3A, 3B). Secondary peaks of parasites occur far too frequently, and this is likely to be explained by the high reticulocyte invasion rates we estimated experimentally. The model output suggests that after the first peak of parasites, when RBC densities are low and starting to be replenished by more reticulocytes, these new RBCs are being infected too quickly compared to real infections. The model comes closer to generating qualitatively reasonable dynamics when we use the estimates of invasion rates from Mideo et al. (2008b), where reticulocyte invasion rates are an order of magnitude lower (fig. 3C). Assuming our experimental measurements are accurate (an issue we address in the discussion), this suggests that (i) the remaining parameters in the model of Mideo et al. (2008b) need to change to accommodate the experimentally measured invasion rates, or (ii) the model framework itself is incomplete. We explore these options by fitting a number of model variations to our newly collected infection data as well as to the data set we fitted in Mideo et al. (2008b).
Experimental and predicted infection dynamics of a relatively more virulent genotype (AS, top row) and less virulent genotype (DK, bottom row). A, Average densities ± 1 SEM for red blood cells (RBCs; solid lines, filled circles) and parasites (dashed lines, open circles). B, Infection dynamics as predicted by the model of Mideo et al. (2008b), using experimentally estimated invasion rates and average burst sizes (in control mice, col. 3, table 1). All other parameters are set to the genotype-specific values published in Mideo et al. (2008b). Secondary parasite peaks are predicted to occur with greater frequency and magnitude than is seen in the data. Even more extreme and less congruent results occur when burst sizes are set to the age-specific values listed in the last two columns of table 1, since these high reticulocyte burst sizes further enhance parasite growth when RBC densities start to rebound and the relative density of reticulocytes increases (not shown). C, Infection dynamics as predicted when invasion rates are set to the published best estimates in Mideo et al. (2008b) rather than experimentally estimated values. It this case, the reticulocyte invasion rate is an order of magnitude lower than for mature RBCs, which slows and dampens secondary parasite peaks.
The Next Iteration of Model Development and Comparisons
We use residual analysis to assess model fit to data (app. B). Outlying residuals more than about 3 standard deviations from 0 and serially correlated residuals suggest poor fit. Mean standardized residuals outside of their 95% predictive intervals suggest model bias and again poor fit. It can be seen in figure B2 that the hybrid model (with immunity and RBC age structure) provides a better fit than the model of Mideo et al. (2008b; with RBC age structure only). The likelihood ratio (Bayes factor) of the hybrid model to the model of Mideo et al. (2008b) is 
(quoted as with 1 standard error; table 2). Thus, there is overwhelming evidence that the data are better explained by a model that includes immune responses and that these responses are important regulators of malaria infections in CD4+-depleted mice. By removing each immune component and refitting the model, we can determine the evidence for each component and for each parasite genotype (table 2). The Bayes factors reported in table 2 show that there is overwhelming support for merozoite clearance in AS-infected mice but no support in DK-infected mice; there is some support for parasitized RBC clearance in AS- and DK-infected mice; and there is some support for unparasitized RBC clearance in AS-infected mice but none in DK-infected mice. Finally, we tested whether RBC age structure was necessary to explain the data, that is, if there was evidence for RBC age-specific invasion rates and burst sizes. The likelihood ratio of the hybrid model (with age structure) to the model without age structure is 
(table 2). Thus, we find no evidence that RBC age structure is required to explain the data. Overall, the most likely model is therefore one that does not contain RBC age structure, does allow for higher death rates of multiply parasitized RBCs (table 2), and does contain immune responses that independently target merozoites, parasitized RBCs, and unparasitized RBCs, although not all of these responses are necessary to explain the dynamics of every individual mouse. (Despite this, we refer to the model with all forms of immunity as the “most likely” since it is the minimal model required to explain all of the data, given a particular choice of parameters for each individual data set. In other words, setting key immune parameters to zero will generate infection dynamics that match some, but not all, of the individual data sets.)
Model comparisons using Bayes factors (likelihood ratios of two competing models) and biological interpretations
| Genotype | BF ± 1 SE | Interpretation of BF |
|---|---|---|
| A. Immune response required? Comparing hybrid model with Mideo et al. (2008b) model: | ||
| AS + DK | (5.5 ± 3.0) × 1052 | Immune responses required to explain data |
| B. Age-structure required? Comparing hybrid model with Miller et al. (2010) model: | ||
| AS + DK | .44 ± .08 | Age structure not required to explain data |
| C. Removal of extra death rate for multiply parasitized RBCs: | ||
| AS + DK | (1.1 ± .2) × 1016 | Increased death rate required to explain data |
| D. Removal of immune-mediated merozoite clearance from model: | ||
| AS | (1.1 ± 1.1) × 1016 | Decisive evidence of clearance of AS merozoites |
| DK | 1.0 ± .1 | No evidence of clearance of DK merozoites |
| E. Removal of immune-mediated parasitized RBC clearance from model: | ||
| AS | 3.6 ± .4 | Some evidence of parasitized RBC clearance |
| DK | 2.8 ± .7 | Some evidence of parasitized RBC clearance |
| F. Removal of immune-mediated unparasitized RBC clearance from model: | ||
| AS | 4.9 ± .5 | Some evidence of unparasitized RBC clearance |
| DK | .40 ± .05 | No evidence of unparasitized RBC clearance |
The fits of the most likely model to individual data sets are shown in figure 4. This model generates new inferences about differences between genotypes that can be compared to those in Mideo et al. (2008b). First, our original model predicted that the burst size of the more virulent genotype, AS, is higher than the less virulent genotype, DK, and this prediction has proven robust to changes in the model structure. AS-infected RBCs have a median burst size of 5.9 merozoites (5.4–6.5, 95% credible interval [CI]; fig. 5, left) and DK parasites have a median burst size of 5.1 merozoites (4.7–5.6, 95% CI). The likelihood ratio of the model allowing for genotype specific burst sizes to the model with genotype nonspecific burst sizes is 
. Thus, genotype specific burst sizes are about 10 times more probable than genotype nonspecific burst sizes. Second, the likelihood ratio of the model allowing genotype specific invasion rates to the model with genotype nonspecific invasion rates is 
. Thus, in contrast to our previous findings (Mideo et al. 2008b), the current model provides no support for the prediction that the more virulent genotype has a higher invasion rate of mature RBCs (fig. 5, right). These two inferences are in line with what we have measured experimentally (fig. 2, right panels). Third, in Mideo et al. (2008b), as in previous work (Haydon et al. 2003), we found evidence that RBC production differed across hosts infected with different parasite genotypes, in particular that hosts infected with DK increased RBC production more in response to anemia. Here we find little evidence for a difference: hosts infected with AS parasites have a median upregulation rate of 0.19 h−1 (0.091–0.31, 95% CI; fig. 5, middle) and hosts infected with DK parasites have a median upregulation rate of 0.26 h−1 (0.16–0.37, 95% CI). The likelihood ratio of the model with genotype specific upregulation rates to the model with genotype nonspecific upregulation rates is 
. This discrepancy from our original predictions about RBC production in Mideo et al. (2008b) is largely a reflection of our stricter statistical criteria for assessing differences between genotypes (hyperparameters) compared to our previous work where we used an arbitrary criterion, that is, a difference of 10% or greater in the median predicted parameter between genotypes—like the difference in upregulation rates just reported—was considered evidence of a significant effect.
The next iteration of model fitting. Fits of reduced hybrid model (the “most likely” model) to newly collected parasite and red blood cell (RBC) densities. Crosses are data, gray regions correspond to 95% posterior predictive intervals of the predicted dynamics, and solid lines give the best-fit solutions for each individual mouse.
Parasite genotype-specific model inferences. Plots show marginal distributions of parasite genotype-level (hyper)parameters. Histograms represent 4,000 best-fit estimates from refitting the model to all data. By comparing the fit of the most likely model that allows for these genotype-specific parameters, to the model that includes only a single parameter for both genotypes, the support for genotype-specific differences can be determined. We find that there is important variation in burst sizes across genotypes, but not in the upregulation of host red blood cell production or in invasion rates.
Finally, some important genotypic differences in the interactions between parasites and host immunity emerge. As outlined above, immune responses are generally not necessary to explain the dynamics of infections with the less virulent genotype, DK, while they are of key importance in infections with the more virulent genotype, AS (table 2). There is considerable individual variation in the importance of the various immune responses in infections (fig. B3–B5), but overall there is strong evidence for an effect of merozoite clearance on the dynamics of AS infections and a weak effect of parasitized and unparasitized RBC clearance (table 2). These same differences emerge when we refit the data that was used in Mideo et al. (2008b), that is, the data from Barclay et al. (2008). Specifically, the data from the more virulent genotype is about 1021 times more likely given the model with immunity, while immune responses are not necessary to explain the data for the less virulent genotype (table B3).
Discussion
We have experimentally tested the predictions made by our earlier modeling work (Mideo et al. 2008b) about what underlies variation in the within-host dynamics of rodent malaria infections. We found support for some predictions, for example, burst size was higher for the more virulent genotype, but we failed to find support for others. Assuming our experimental estimates of parasite traits were accurate, we showed that, given these parasite trait estimates as parameters, the original model (Mideo et al. 2008b) was unable to reproduce the average dynamics observed in our new experiment. We used our new infection dynamic data to inform another iteration of model development and selection, once again allowing parameters that govern parasite traits to vary and finding the best estimates of these. We find the most support for a model where infection dynamics are regulated by RBC availability and multiple forms of immune responses, but not by RBC age structure (i.e., the model of Miller et al. 2010).
Despite a substantially more complex model, we again find that burst sizes play a role in underlying observed variation across the rodent malaria genotypes studied, with the more virulent genotype predicted to have a higher burst size. However, our new modeling results also suggest that differences in the interactions between parasites and host immune responses are of key importance, despite the experimental depletion of CD4+ T cells. Immune responses in these CD4+-depleted mice are necessary to explain the dynamics of the more virulent genotype, which suffers a noticeable drop in parasite densities around day 10 postinfection before recrudescing, but not the less virulent genotype, for which changes in parasite densities after peak are much more subtle. Studying infections with a different pair of parasite genotypes in a different host immune background, Miller et al. (2010) also point to a role for immune responses in generating variation in infection dynamics. Specifically, they found that differences in dynamics were explained by initial inoculation dose, which they hypothesized could be the result of variation in the ability of parasites to evade fast-acting innate immunity after inoculation.
In addition to affecting dynamics, the role differential immune responses play in explaining variation in virulence across the two parasite genotypes studied here remains unclear. The more virulent genotype, AS, leads to “sicker” hosts (e.g., greater weight loss; data not shown), yet differences in parasite and RBC dynamics in infected hosts seem less straightforward. For instance, the less virulent genotype, DK, maintains high parasite densities for longer. However immune-mediated virulence, or immunopathology, is increasingly being acknowledged as an important phenomenon in malaria infections and more broadly (e.g., Stevenson and Riley 2004; Graham et al. 2005; Lamb et al. 2006). Previous experimental work on Plasmodium chabaudi has demonstrated that infections with more virulent parasite genotypes result in a greater induction of innate immune responses, even after controlling for differences in parasite density (Long et al. 2006, 2008). This suggests that immune-mediated virulence explains some of the variation in damage done to hosts by different genotypes (Graham et al. 2005; Long and Graham 2011). Ideally, a next step would be to test our new model inferences by quantifying immune killing in infections with different parasite genotypes, but determining what to measure and when is nontrivial. Consequently, immunological data remains a relatively untapped source of information for testing and refining models.
Our experimental and theoretical results describe how parasite traits and host factors that regulate infection dynamics differ across infections with two malaria parasite genotypes (e.g., the more virulent genotype produces more progeny parasites in infected cells and experiences greater regulation by immune responses). This raises the question of whether parasite phenotypes are the cause or the consequence of differences in host factors. In other words, does the higher burst size of the more virulent genotype result in these infections being more immunogenic, while the less virulent genotype has a more “prudent” strategy of flying under the radar of immune responses? Or might a functional relationship work in the other direction, with higher burst size being a strategy for dealing with a stronger immune response, that is, by providing safety in numbers from immune attack (Mideo 2009; O’Donnell et al. 2011)? It is interesting to note that when alone in a host, the less virulent genotype can establish a relatively successful acute infection, while in mixed infections this genotype is outcompeted by the more virulent one (Barclay et al. 2008). To what extent this is the result of resource- or immune-mediated competition is unclear. In CD4+-depleted hosts, competition was more severe than in hosts with normal immune systems (Barclay et al. 2008), which could be because either resource competition is exacerbated (since parasite densities are higher) or innate immune responses are stronger (since antigen load is higher). In either case, the higher burst size of the more virulent genotype would seem to be advantageous in multiple infections. Virulent genotypes are better competitors (de Roode et al. 2005; Bell et al. 2006), and it is an interesting possibility that differences in burst sizes could explain why. Whether burst sizes correlate with competitive ability and virulence across a greater range of genotypes, and in hosts with intact immune systems, is an open question.
In addition to proposing some novel parasite biology, our results highlight the next major challenge for integrating theory and data for studying within-host infection dynamics. In Mideo et al. (2008b), we predicted that reticulocytes provided a better habitat for malaria parasites, resulting in higher burst sizes in these younger cells. Biologically, this seemed plausible since these cells have residual intracellular components that may allow more efficient remodeling of the RBC, which is required for parasite growth (Haldar and Mohandas 2007). Given the results in table 1 and the new modeling results presented here, we do not find evidence to support this prediction. Rather, our experimental results suggest that burst sizes are variable (plastic) and that burst sizes increase across all ages of RBC in response to PHZ or some correlated effect. In addition to altering RBC age structure, PHZ changes RBC densities and elicits cytokine production, potentially altering environmental constraints on parasite replication. Plasticity in malaria parasite traits in response to changes in the within-host environment is well documented (Dixon et al. 2008; Reece et al. 2009; Pollitt et al. 2011a). This presents at least two potential issues. First, should parasite traits be modeled as constant parameters, as we have done for burst sizes and invasion rates, rather than as functions that change over the course of an infection? Second, how reliable are our experimental estimates of parasite traits given that they are from only a single snapshot of the course of infection? On the theory side, adding parasite plasticity into dynamical models is technically feasible, though conceptually difficult. A key challenge will be to mathematically define parasite traits in ways that do not just add complexity to models arbitrarily but rather are motivated by known biology. Currently, our models are in essence generating predictions about “average” traits over the course of infections. On the experimental side, measuring many parasite traits over the entire course of individual infections is a major technical and logistical challenge.
Our estimates of burst sizes and invasion rates are dependent on data derived from in vitro cultures, since capturing sufficient numbers of mature schizonts would require destructive sampling and prevent us from simultaneously tracking infection dynamics. The reliability of our estimates of parasite traits is unknown, but this is currently the best available method for obtaining the data required, and is standard practice for studies of human malaria parasites (Trager and Jensen 1997). However, even if our experimental estimates of parasite traits are inaccurate, the inferences we draw from our new modeling work should be largely unaffected. These estimates were used to further constrain the priors on model parameters, but the parameters were still fitted (rather than set).
Even models that turn out to be inaccurate have the potential to generate new questions and guide further work (e.g., Epstein 2008). Indeed, the model of Mideo et al. (2008b) raised the question of whether parasite exploitation traits explain variation in infection dynamics and provided specific hypotheses that could be empirically tested. Ultimately, if one is interested in generating new understanding then model refinement is a key step in this process. The “best” model presented here not only reliably reproduces the observed dynamics, but it also generates inferences about genotype-specific burst sizes and invasion rates that fit with what we have measured experimentally (and quantitative estimates of burst sizes match up remarkably well). Our modeling approach has also been improved in a number of ways (described in “Methods” and apps. A and B ). For example, here we fit models to both RBC and parasite data rather than just to RBC data (as in Mideo et al. 2008b), so we have a stricter test of the models. This is due to the fact that in Mideo et al. (2008b) we had data from a single experiment that we were trying to explain. Therefore, we fit models to only half of that data set in order to “test” the model with the other half. Now, we are developing an approach where models and their predictions are tested by further experimental iterations and so are best when informed by as much data as possible. Even a visual inspection reveals that while the original fits in Mideo et al. (2008b) are surprisingly good considering the model was not fit to parasite data, the new model fitted to this same data is much better at capturing both RBC and parasite dynamics (fig. B7).
Mathematical models provide a tool for formalizing hypotheses and identifying plausible explanations, but of course “the ability of a model to describe the data does not make it correct” (Kochin et al. 2010, p. 1). Clearly models need to be tested, but experiments that can dissect the dynamic interactions between hosts and parasites are extremely challenging. Finding simple perturbations that can help to distinguish between competing models is a key goal. For example, a number of recent modeling studies have suggested that preferential invasion of old RBCs by P. chabaudi parasites has a significant effect on infection dynamics (Antia et al. 2008; Mideo et al. 2008b), while our results call into question the role of RBC age structure. To resolve this discrepancy, experimentally manipulating RBC age structure in a way that is not also likely to impact immunity (e.g., using the hormone erythropoietin; Reece et al. 2005) could provide useful data. Finally, there remains much scope for integrating other genetically controlled parasite traits into mathematical models, for example, production of transmissible parasite forms (gametocytes), in order to gain a better understanding of variation in infection dynamics. Our previous work showed that conversion to gametocytes was not an important regulator of infection dynamics, but this was modeled as a static parasite trait (Mideo et al. 2008b), and experiments have demonstrated considerable variation in this trait across parasite genotypes and within infections (Wargo et al. 2007; Reece et al. 2008; Pollitt et al. 2011a). Rodent malaria parasites are probably unique in the availability of relevant data, techniques, and protocols, and while major challenges remain, we have shown that integration between empirical and modeling approaches can be achieved. Through this process, a clearer picture of what regulates infection dynamics, and what generates variation thereof, in this system is emerging.
We thank S. Babayan, K. Fairlie-Clarke, A. Graham, and L. Pollitt for discussion, S. Brown for reagents and assistance with fluorescence-activated cell sorter analysis, and two anonymous reviewers for helpful comments. This work was funded by the Wellcome Trust (N.J.S. was supported by Wellcome Trust Project Grant 082601/Z/07/Z), the Natural Sciences and Engineering Research Council of Canada, and the Canada Research Chairs Program.
Appendix A. Supplementary Methods
Effects of Experimental Manipulation
Some effects of our experimental manipulation of host immunity and RBC age structure can be seen in figure A1. Mean parasite densities for each treatment over the initial peak of infection are shown in figure A1, top row. As expected, parasite densities are higher in CD4+ T-depleted mice than in immune-intact control mice (Barclay et al. 2008). Mean RBC densities for each treatment over the initial peak of infection are shown in figure A1, bottom row. PHZ treatment results in a marked decrease in RBC densities. Figure A2 illustrates that PHZ-treated mice had significantly higher proportions of bloodstream reticulocytes and that a significantly higher proportion of infected RBCs were reticulocytes, across the relevant days postinfection.
Model Derivation
Basic Structure
Incorporating the age structure of the Mideo et al. (2008b) model into the model of Miller et al. (2010) yields the following basic model structure. On the (i+1)th day postinfection, just after all infected red blood cells (RBCs) burst, the densities of merozoites, 
, reticulocytes on their jth day in the bloodstream, 
, and mature RBCs, Ni+1 are given by
is the total RBC density τ days before i. The parameter τ allows RBC production in response to anemia to be time-lagged, since RBC precursors take time to develop in the bone marrow. The parameters 
and 
describe the immune clearance rates of parasitized and unparasitized RBCs as a function of time i. These parameters are described in more detail below. The parameters λR and λN define the average number of merozoites (that survive clearance in the bloodstream) per reticulocyte and mature RBC, respectively. To find expressions for these parameters, we note that an individual merozoite has a probability of infecting a reticulocyte given by
describes the immune clearance rate of merozoites as a function of time, i. Multiplying these probabilities by the initial density of merozoites and dividing by the respective densities of RBCs gives the average number of merozoites per reticulocyte:
The model of Mideo et al. (2008b) discussed in the main text is recovered by setting all immune clearance parameters (
, 
, and 
) to 0. This is slightly different than the published model in Mideo et al. (2008b) since it tracks infections in time steps of thirds of a day rather than whole days and multiply infected RBCs are tracked separately. This allows the model to make predictions about the densities of unparasitized, singly parasitized, and multiply parasitized cells that should more accurately correspond to what is measured experimentally. However, all underlying biological hypotheses described by the model remain unchanged.
Immune Clearance
Three independent functions describe the immune clearance rates of merozoites, parasitized RBCs, and unparasitized RBCs over the course of infection (denoted 
, 
, and 
). We assume that there are at maximum two “windows” of immune activity, and each window is described by four parameters (so each function is fully described by eight parameters). A schematic “clearance rate function” for parasitized RBCs is given in figure 1 of the main text. The parameter sx defines the day postinfection when immune clearance begins, cx represents the maximum clearance rate, rx represents the time it takes to reach that maximum rate (or “rise time”), and lx represents the duration of immune clearance (note that 
). The subscript x denotes the cells being cleared (either m for merozoites, p for parasitized RBCs, or u for unparasitized RBCs). For each immune clearance function, the two windows of immune activity will be described by a different set of parameters.
Model Fitting
The fitted model parameters and prior distributions are given in table A1. The prior distributions were either taken from the literature or based on our experimental measurements. The prior distributions for burst size, RBC upregulation rate, and invasion rate are specified with hyperparameters (table A2); hyperparameters essentially specify the distribution from which the individual-level parameters are randomly drawn. We do this for these parameters because we want to examine if they are parasite genotype-specific (see below).
As well as fitting the hybrid model (Mideo et al. 2008b; Miller et al. 2010) with age structure and immune responses, we fit a model without age structure (Miller et al. 2010) by setting reticulocyte and normocyte burst sizes and invasion rates equal (
, 
) and a model without immune responses (Mideo et al. 2008b) by setting immune-mediated clearance rates to 0 (
). We also fit models without age structure and with each of the three immune components removed by setting their respective clearance rates to 0.
Miller et al. (2010) showed that there is a nonidentifiability between RBC invasion rate β and natural death rate of merozoites μ. This meant that it was impossible to obtain separate estimates for them and all that could be obtained was their ratio. In Mideo et al. (2008b), however, we did not fit μ but assumed it had a value of 48 d−1; thus, an estimate for β was obtainable. In order to test for genotype-specific differences in invasion rates we must therefore fix μ and assume that it is non-genotype-specific. Miller et al. (2010) estimated μ/β to be of the order 106, taking 
as in Mideo et al. (2008b) gives β−1 to be of the order 104. The inverse of β is estimated rather than β because we have no prior knowledge about its upper bound.
In order to determine if there are genotype-specific differences in parameters we estimate the hyperparameters of the two strains in the model without age structure. For example, the prior distribution of burst size for AS-infected mice is 
and for DK-infected mice 
. Terms 
, 
, 
are hyperparameters that we estimate. If there are genotype specific differences in parameter estimates we expect to see a difference in the posterior distributions of 
and 
, that is, the mean burst sizes of mice infected with AS and DK parasites, respectively. The variance hyperparameter 
, is assumed to be non-strain-specific as it is of little interest here. We can also suppose that the burst sizes for all mice are randomly drawn from the same population, that is, non-strain-specific burst sizes, which means we take the prior as 
over all mice. We must also specify hyperpriors on the hyperparameters; these are given in table A2 and are chosen to be conjugate to the priors.
By estimating genotype-specific and non-genotype-specific hyperparameters we can calculate the (marginalized) likelihood ratio of the data supposing genotype-specific hyperparameters to the data supposing non-genotype-specific hyperparameters. Without any prior preference for these two suppositions their prior odds are 1∶1. Thus, by Bayes’s theorem, our posterior odds for these suppositions are equal to their likelihood ratio. If 
is the likelihood of mouse m’s data given model M, then our odds on genotype-specific hyperparameters is given by the ratio
Infection dynamics of different treatments. Treatments are grouped by parasite genotype: A, More virulent parasite genotype AS; B, less virulent parasite genotype DK. Top row shows mean parasite density, and bottom row shows mean RBC densities over time. Error bars, ±1 SEM. Colors indicate what treatments in addition to infection the hosts received: black, none (immune-intact control); blue, CD4+ depletion (CD4−); red, PHZ and CD4+ depletion. Over the course of the initial peak, parasite densities were higher in CD4+-depleted hosts as compared with immune-intact hosts. PHZ treatment (2 days before infection, i.e., day −2) resulted in a decrease in red blood cell densities. This effect coincided with the initial stages of infection.
Effect of phenylhydrazine (PHZ) on red blood cell age structure in CD4+ depleted mice. The mean proportion of all red blood cells (RBCs) that are reticulocytes (A) and infected RBCs that are reticulocytes (B) over 3 consecutive days in untreated mice and those that received PHZ treatment. Error bars, ±1 SEM. At the early stages of infection, mice treated with PHZ had a significantly higher proportion of reticulocytes, as expected, and subsequently, a significantly higher proportion of all infected cells were also reticulocytes.
Model parameters and prior distributions of hybrid model
| Parameter | Description | Value or prior | Source |
|---|---|---|---|
| ωR, ωN | Burst sizes in reticulocytes and normocytes | N(νω, σω2) | See table A2 |
| βR, βN | Invasion rate of reticulocytes, normocytes ([cells/μL]−1 s−1) | Not fitted | |
| ρ | log10(βR/βN) | N(0, .32) | Hetzel and Anderson 1996; Antia et al. 2008; Mideo et al. 2008b |
| μ | Natural death rate of merozoites (day−1) | 48 | Garnham 1966; Mideo et al. 2008b |
| β−1 | Inverse of invasion rate (cells s/μL) | Exp(l) | See table A2 |
| θ | Rate of upregulation of erythropoeisis (day−1) | N(νθ, σθ2) | See table A2 |
| τ | Time lag in erythropoeisis (day) | U(0, 6) | Chang et al. 2004; Mideo et al. 2008b |
| d | Natural death rate of RBCs (day−1) | .025 | van Putten 1958; Bannerman 1983 |
| dm | Increased death rate of multiply-parasitized RBCs (day−1) | Exp(1) | Miller et al. 2010 |
| sm, sp, su | Start day of immunity targeting merozoites, parasitized RBCs, unparasitized RBCs | U(0, 21) | Miller et al. 2010 |
| rm, rp, ru | Rise time of immunity targeting merozoites, parasitized RBCs, unparasitized RBCs | U(0, 21) | Miller et al. 2010 |
| cm | Maximum level of immunity targeting merozoites (cells/s) | Exp(108) | Miller et al. 2010 |
| cp, cu | Maximum clearance rate of immunity targeting parasitized RBCs, unparasitized RBCs (day−1) | Exp(1) | Miller et al. 2010 |
| lm, lp, lu | Duration of immunity targeting merozoites, parasitized RBCs, unparasitized RBCs | U(0, 21) | Miller et al. 2010 |
| P0 | Initial parasite density (parasites/μL) |  |
Miller et al. 2010 |
| N0 | Initial RBC density (RBCs/μL) | NT(6.5 × 106, 1012) |
Hyperparameters and their hyperprior distributions.
| Hyperparameter | Description | Hyperprior | Source |
|---|---|---|---|
| νω | Mean burst size | N(6, .52) | Garnham 1966; Carter and Walliker 1975; Carter and Diggs 1977; Mideo et al. 2008b |
| σω2 | Variance of burst size | InvGam(2, 12) | Garnham 1966; Carter and Walliker 1975; Carter and Diggs 1977; Mideo et al. 2008b |
| νβ | Mean inverse invasion rate | InvGam(1, 104) | Mideo et al. 2008b |
| νθ | Mean red blood cell upregulation rate | N(0.4, 0.22) | Haydon et al. 2003; Mideo et al. 2008b |
| σθ2 | Variance of red blood cell upregulation rate | InvGam(2, 0.22) | Haydon et al. 2003; Mideo et al. 2008b |
Appendix B. Supplementary Results
Statistical Analyses of Experimental Data
Tables B1 and B2 present statistics from linear mixed effects models of invasion rate variation and linear models of burst size variation. Histograms showing the burst size observations across treatment groups are shown in figure B1.
Assessing Model “Goodness of Fit”
We compare the fit of our original model (Mideo et al. 2008b) with that of the hybrid model, which includes multiple forms of immunity, by plotting the overlaid standardized residuals for parasite and RBC densities for each model (fig. B2). When the mean of the standardized residuals lies outside the 95% (Bonferroni corrected for multiple tests) predictive interval (i.e., the red line falls outside of the dashed lines), this suggests that the model is over- or underestimating the data, and is suggestive of a poor fit. While there are clearly some time points for which neither model captures the data well, overall the hybrid model is fitting better (fig. B2, right ). The early time points for which the RBC densities are fitted poorly is not particularly surprising since there is a lot of unexplained variation in the data at these very early days postinfection.
Model Inferences
The most likely model includes immune responses that independently target merozoites, parasitized RBCs, and unparasitized RBCs, although not every response is necessary to explain the dynamics of every individual mouse. This individual variation is evident in the posterior predictive intervals for the different immune responses, depicted below (figs. B3–B5). Overall, immune responses targeting merozoites and parasitized RBCs are more important for explaining the dynamics of infections with the more virulent clone AS than with the less virulent clone DK. The marginal posterior distributions for all other parameters are given in figure B6.
Refitting Original Data
Given that the model we presented in Mideo et al. (2008b) was assessed as providing a “good fit” to the data used in that study, we refit that data (originally from Barclay et al. 2008) to the Mideo et al. (2008b) model and the Miller et al. (2010) model using the Bayesian framework of this study. As with the new data we explore in the main text, the immune response is necessary to explain the data from infections with the more virulent genotype, but not the avirulent genotype (table B3).
In figure B6, we compare the best fits from Mideo et al. (2008b), where models were fitted only to parasite data, with the new hybrid model fitted to this same data.
Histogram of burst size observations across treatments. Burst sizes were estimated by counting the number of merozoites in at least 25 mature schizonts for each individual infected mouse. The distributions here are plotted from pooling this data according to the genotype of the infecting parasites (DK or AS) and phenylhydrazine treatment (PHZ or none).
Standardized residuals of the fit of the age-structured model of Mideo et al. (2008b; left) and the fit to the hybrid model (right) to newly collected data. Top, Parasite density; bottom, red blood cell density. Each cross represents the standardized residual on a particular day for an individual mouse. The red line joins the means of the standardized residuals for each day, and the dashed lines denote the 95% (Bonferroni corrected for multiple tests) predictive intervals of the mean standardized residual assuming the model is true. The Y-axis is scaled to units of standard deviations.
Posterior predictive interval (PPI) of immune responses targeting parasitized red blood cells (RBCs). Solid lines give best-fit function describing clearance rate. Light gray regions correspond to 95% PPI; dark gray regions correspond to 50% PPI.
Posterior predictive interval (PPI) of immune responses targeting unparasitized red blood cells (RBCs; i.e., bystander death). Solid lines give best-fit function describing clearance rate. Light gray regions correspond to 95% PPI; dark gray regions correspond to 50% PPI.
Posterior predictive interval (PPI) of immune responses targeting merozoites. Solid lines give best-fit function describing total clearance. Light gray regions correspond to 95% PPI; dark gray regions correspond to 50% PPI.
Marginal distributions of fitted parameters for the most likely reduced hybrid model. Each row corresponds to an individual mouse. White panels are for individuals infected with the more virulent genotype AS; gray panels are for individuals infected with the less virulent genotype DK. Dashed lines indicate the prior distributions on each parameter. Units are given in table A1.
Comparison of model fits. A, Original best fit of red blood cell (RBC) age-structured model with no immunity, fitted to data from Barclay et al. (2008). Redrawn from Mideo et al. (2008b). B, Fit of hybrid model (including RBC age-structure and immune responses) to the same data set.
Analysis of red blood cell (RBC) invasion rates in CD4+-depleted mice
| LRT (χ2 ) | P | |
|---|---|---|
| Minimal model: | ||
| RBC age | NA | |
| Genotype | NA | |
| Genotype∶RBC age |  |
.004 |
| Nonsignificant terms deleted from maximal model: | ||
| Mass of mouse |  |
.538 |
Analysis of burst sizes in CD4+-depleted mice
| F | P | |
|---|---|---|
| Minimal model: | ||
| Genotype | F1, 15 = 11.021 | .005 |
| Phenylhydrazine | F1, 15 = 4.067 | .062 |
| Nonsignificant terms deleted from maximal model: | ||
| Mass of mouse | F1, 14 = .893 | .361 |
| Parasite density | F1, 13 = .181 | .677 |
| Uninfected red blood cell density | F1, 12 = .317 | .584 |
| Genotype∶phenylhydrazine | F1, 11 = .975 | .345 |
Comparing the fit of original model from Mideo et al. (2008b) with the Miller et al. (2010) model of the main text using Bayes factors (BFs)
| Genotype | BF | Interpretation of BF |
|---|---|---|
| AS | 1.1 × 1040 | AS data are overwhelming more likely under the Miller et al. (2010) model |
| DK | 4.5 × 10−6 | DK data are overwhelmingly more likely under the Mideo et al. (2008b) model, suggesting that the Miller et al. (2010) model overfits the data |
| All | 4.9 × 1036 | Conclusions are the same as those for the newly collected data: immune responses are important for explaining dynamics of more virulent (AS) but not for the less virulent (DK) genotype |
Literature Cited
Antia, R., A. Yates, and J. C. de Roode. 2008. The dynamics of acute malaria infections. I. Effect of the parasite’s red blood cell preference. Proceedings of the Royal Society B: Biological Sciences 275:1449–1458. Bannerman, R. M. 1983. Hematology. Pages 293–312 in J. D. Small, H. L. Foster, and J. G. Fox, eds. The mouse in biomedical research. Academic Press, New York. Barclay, V. C., L. Råberg, B. H. K. Chan, S. Brown, D. Gray, and A. F. Read. 2008. CD4 + T cells do not mediate within-host competition between genetically diverse malaria parasites. Proceedings of the Royal Society B: Biological Sciences 275:1171–1179. Beale, G. H., R. Carter, and D. Walliker. 1978. Genetics. Pages 213–245 in R. Killick-Kendrick and W. Peters, eds. Rodent malaria. Academic Press, London. Bell, A. S., J. C. de Roode, D. Sim, and A. F. Read. 2006. Within-host competition in genetically diverse malaria infections: parasite virulence and competitive success. Evolution 60:1358–1371. Carter, R., and C. L. Diggs. 1977. Plasmodia of rodents. Pages 359–461 in J. M. Kreier, ed. Parasitic protozoa. Vol. 3. Academic Press, New York. Carter, R., and D. Walliker. 1975. New observations on malaria parasites of rodents of Central African Republic: Plasmodium vinckei petteri subsp. nov. and Plasmodium chabaudi Landau, 1965. Annals of Tropical Medicine and Parasitology 69:187–196. Chang, K. H., M. Tam, and M. M. Stevenson. 2004. Modulation of the course and outcome of blood-stage malaria by erythropoietin-induced reticulocytosis. Journal of Infectious Diseases 189:735–743. Cromer, D., K. J. Evans, L. Schofield, and M. P. Davenport. 2006. Preferential invasion of reticulocytes during late-stage Plasmodium berghei infection accounts for reduced circulating reticulocyte levels. International Journal for Parasitology 36:1389–1397. Cromer, D., J. Stark, and M. P. Davenport. 2009. Low red cell production may protect against severe anemia during a malaria infection: insights from modeling. Journal of Theoretical Biology 257:533–542. de Roode, J. C., R. Pansini, S. J. Cheesman, M. E. H. Helinski, S. Huijben, A. R. Wargo, A. S. Bell, et al. 2005. Virulence and competitive ability in genetically diverse malaria infections. Proceedings of the National Academy of Sciences of the USA 102:7624–7628. Dixon, M. W. A., J. Thompson, D. L. Gardiner, and K. R. Trenholme. 2008. Sex in Plasmodium: a sign of commitment. Trends in Parasitology 24:168–175. Epstein, J. M. 2008. Why model? Journal of Artificial Societies and Social Simulation 11:12. Garnham, P. C. C. 1966. Malaria parasites and other haemosporidia. Blackwell Scientific, Oxford. Good, M. F., and D. L. Doolan. 1999. Immune effector mechanisms in malaria. Current Opinion in Immunology 11:412–419. Graham, A. L., J. E. Allen, and A. F. Read. 2005. Evolutionary causes and consequences of immunopathology. Annual Review of Ecology, Evolution and Systematics 36:373–397. Grech, K., K. Watt, and A. F. Read. 2006. Host-parasite interactions for virulence and resistance in a malaria model system. Journal of Evolutionary Biology 19:1620–1630. Haldar, K., and N. Mohandas. 2007. Erythrocyte remodeling by malaria parasites. Current Opinion in Hematology 14:203–209. Haydon, D. T., L. Matthews, R. Timms, and N. Colegrave. 2003. Top-down or bottom-up regulation of intra-host blood-stage malaria: do malaria parasites most resemble the dynamics of prey or predator? Proceedings of the Royal Society B: Biological Sciences 270:289–298. Hetzel, C., and R. M. Anderson. 1996. The within-host cellular dynamics of bloodstage malaria: theoretical and experimental studies. Parasitology 113:25–38. Holmes, E. C., and S. S. Burch. 2000. The causes and consequences of genetic variation in dengue virus. Trends in Microbiology 8:74–77. Kochin, B. F., A. J. Yates, J. C. de Roode, and R. Antia. 2010. On the control of acute rodent malaria infections by innate immunity. PLoS One 5:e10444. Lamb, T. J., D. E. Brown, A. J. Potocnik, and J. Langhorne. 2006. Insights into the immunopathogenesis of malaria using mouse models. Expert Reviews in Molecular Medicine 8:1–22. Langhorne, J., B. Simon-Haarhaus, and S. J. Meding. 1990. The role of CD4 + T cells in the protective immune response to Plasmodium chabaudi in vivo. Immunology Letters 25:101–108. Long, G. H., and A. L. Graham. 2011. Consequences of immunopathology for pathogen virulence evolution and public health: malaria as a case study. Evolutionary Applications 4:278–291. Long, G. H., B. H. K. Chan, J. E. Allen, A. F. Read, and A. L. Graham. 2006. Parasite genetic diversity does not influence TNF-mediated effects on the virulence of primary rodent malaria infections. Parasitology 133:673–684. ———. 2008. Experimental manipulation of immune-mediated disease and its fitness costs for rodent malaria parasites. BMC Evolutionary Biology 8:128. Mackinnon, M. J., D. M. Gunawardena, J. Rajakaruna, S. Weerasingha, K. N. Mendis, and R. Carter. 2000. Quantifying genetic and nongenetic contributions to malarial infection in a Sri Lankan population. Proceedings of the National Academy of Sciences of the USA 97:12661–12666. Mackinnon, M. J., and A. F. Read. 1999. Genetic relationships between parasite virulence and transmission in the rodent malaria Plasmodium chabaudi. Evolution 53:689–703. McMorran, B. J., V. M. Marshall, C. de Graaf, K. E. Drysdale, M. Shabbar, G. K. Smyth, J. E. Corbin, et al. 2009. Platelets kill intraerythrocytic malarial parasites and mediate survival to infection. Science 323:797–800. Mideo, N. 2009. Parasite adaptations to within-host competition. Trends in Parasitology 25:261–268. Mideo, N., T. Day, and A. F. Read. 2008 a. Modelling malaria pathogenesis. Cellular Microbiology 10:1947–1955. Mideo, N., V. C. Barclay, B. H. K. Chan, N. J. Savill, A. F. Read, and T. Day. 2008 b. Understanding and predicting strain-specific patterns of pathogenesis in the rodent malaria Plasmodium chabaudi. American Naturalist 172:E214–E238. Miller, M. R., L. Råberg, A. F. Read, and N. J. Savill. 2010. Quantitative analysis of immune response and erythropoiesis during rodent malarial infection. PLoS Computational Biology 6:e1000946. Molineaux, L., and K. Dietz. 1999. Review of intra-host models of malaria. Parasitologia 41:221–231. Morrison, L. J., S. McLellan, L. Sweeney, C. N. Chan, A. MacLeod, A. Tait, and C. M. R. Turner. 2010. Role for parasite genetic diversity in differential host responses to Trypanosoma brucei infection. Infection and Immunity 78:1096–1108. Nesse, R. M., and S. C. Stearns. 2008. The great opportunity: evolutionary applications in medicine. Evolutionary Applications 1:28–48. O’Donnell, A., P. Schneider, H. G. McWatters, and S. E. Reece. 2011. Fitness costs of disrupting circadian rhythms in malaria parasites. Proceedings of the Royal Society B: Biological Sciences 278:2429–2436. Paterson, S., and M. E. Viney. 2003. Functional consequences of genetic diversity in Strongyloides ratti infections. Proceedings of the Royal Society B: Biological Sciences 270:1023–1032. Paul, R. E. L., F. Ariey, and V. Robert. 2003. The evolutionary ecology of Plasmodium. Ecology Letters 6:866–880. Pollitt, L. C., N. Mideo, D. R. Drew, P. Schneider, N. Colegrave, and S. E. Reece. 2011 a. Competition and the evolution of reproductive restraint in malaria parasites. American Naturalist 177:358–367. Pollitt, L. C., P. MacGregor, K. Matthews, and S. E. Reece. 2011 b. Malaria and trypanosome transmission: different parasites, same rules? Trends in Parasitology 27:197–203. Pombo, D. J., G. Lawrence, C. Hirunpetcharat, C. Rzepczyk, M. Bryden, N. Cloonan, K. Anderson, et al. 2002. Immunity to malaria after administration of ultra-low doses of red cells infected with Plasmodium falciparum. Lancet 360:610–617. Råberg, L., J. C. de Roode, A. S. Bell, P. Stamou, D. Gray, and A. F. Read. 2006. The role of immune-mediated apparent competition in genetically diverse malaria infections. American Naturalist 168:41–53. Read, A. F., and S. Huijben. 2009. Evolutionary biology and the avoidance of antimicrobial resistance. Evolutionary Applications 2:40–51. Reece, S. E., and J. Thompson. 2008. Transformation of the rodent malaria parasite Plasmodium chabaudi and generation of a stable fluorescent line PcGFP(CON). Malaria Journal 7:183. Reece, S. E., A. B. Duncan, S. A. West, and A. F. Read. 2005. Host cell preferences and variable transmission strategies in malaria parasites. Proceedings of the Royal Society B: Biological Sciences 272:511–517. Reece, S. E., D. R. Drew, and A. Gardner. 2008. Sex ratio adjustment and kin discrimination in malaria parasites. Nature 453:609–614. Reece, S. E., R. S. Ramiro, and D. H. Nussey. 2009. Plastic parasites: sophisticated strategies for survival and reproduction? Evolutionary Applications 2:11–23. Savill, N. J., W. Chadwick, and S. E. Reece. 2009. Quantitative analysis of mechanisms that govern red blood cell age structure and dynamics during anaemia. PLoS Computational Biology 5: e1000416. Stephens, R., and J. Langhorne. 2006. Priming of CD4 + T cells and development of CD4 + T cell memory; lessons for malaria. Parasite Immunology 28:25–30. Stevenson, M. M., and E. M. Riley. 2004. Innate immunity to malaria. Nature Reviews Immunology 4:169–180. Su, Z., A. Fortin, P. Gros, and M. M. Stevenson. 2002. Opsonin-independent phagocytosis: an effector mechanism against acute blood-stage Plasmodium chabaudi AS infection. Journal of Infectious Diseases 186:1321–1329. Trager, W., and J. B. Jensen. 1997. Continuous culture of Plasmodium falciparum: its impact on malaria research. International Journal for Parasitology 27:989–1006. Urban, B. C., R. Ing, and A. M. Stevenson. 2005. Early interactions between blood-stage Plasmodium parasites and the immune system. Current Topics in Microbiology 297:25–70. van Putten, L. M. 1958. The life span of red cells in the rat and the mouse as determined by labelling with Dfp32 in vivo. Blood 13:789–794. Wargo, A. R., J. C. de Roode, S. Huijben, D. R. Drew, and A. F. Read. 2007. Transmission stage investment of malaria parasites in response to in-host competition. Proceedings of the Royal Society B: Biological Sciences 274:2759–2768.
Associate Editor: Pejman Rohani
Editor: Judith L. Bronstein
