Quantitative Analysis of Immune Response and Erythropoiesis during Rodent Malarial Infection

Malarial infection is associated with complex immune and erythropoietic responses in the host. A quantitative understanding of these processes is essential to help inform malaria therapy and for the design of effective vaccines. In this study, we use a statistical model-fitting approach to investigate the immune and erythropoietic responses in Plasmodium chabaudi infections of mice. Three mouse phenotypes (wildtype, T-cell-deficient nude mice, and nude mice reconstituted with T-cells taken from wildtype mice) were infected with one of two parasite clones (AS or AJ). Under a Bayesian framework, we use an adaptive population-based Markov chain Monte Carlo method and fit a set of dynamical models to observed data on parasite and red blood cell (RBC) densities. Model fits are compared using Bayes' factors and parameter estimates obtained. We consider three independent immune mechanisms: clearance of parasitised RBCs (pRBC), clearance of unparasitised RBCs (uRBC), and clearance of parasites that burst from RBCs (merozoites). Our results suggest that the immune response of wildtype mice is associated with less destruction of uRBCs, compared to the immune response of nude mice. There is a greater degree of synchronisation between pRBC and uRBC clearance than between either mechanism and merozoite clearance. In all three mouse phenotypes, control of the peak of parasite density is associated with pRBC clearance. In wildtype mice and AS-infected nude mice, control of the peak is also associated with uRBC clearance. Our results suggest that uRBC clearance, rather than RBC infection, is the major determinant of RBC dynamics from approximately day 12 post-innoculation. During the first 2–3 weeks of blood-stage infection, immune-mediated clearance of pRBCs and uRBCs appears to have a much stronger effect than immune-mediated merozoite clearance. Upregulation of erythropoiesis is dependent on mouse phenotype and is greater in wildtype and reconstitited mice. Our study highlights the informative power of statistically rigorous model-fitting techniques in elucidating biological systems.


Measurement errors in RBC and parasite densities
We fit to the RBC densities using the model: We fit to the measured parasite densities using the model: We estimated the variances in the measurement errors as follows.RBC densities were measured using flow cytometry [1].We have previously estimated the standard deviation of the error, σ N , of this method to be 0.61 × 10 6 (95% CI: 0.44 × 10 6 , 0.84 × 10 6 ) cells/µl [2].Parasite densities were measured using quantitative PCR [1].In another experiment using the same protocol as here (A.Bell, personal communication) two blood samples were simultaneously taken from infected mice at various time points (Figure S1).Assuming that parasite density is log-normally distributed with standard deviation σ p , the difference in two simultaneously taken samples is log-normally distributed with mean 0 and standard deviation √ 2σ p .We therefore estimate σ p to be 0.15.The errors in the samples are consistent with normality (Andersen-Darling statistic: A 2 = 1.08, p > 0.1).Taking a prior on σ N as N(0.6 × 10 6 , 0.1) and a prior on σ p as N(0.15, 0.1) we re-estimated their values by fitting to the data in this paper.The new estimates are σ N = 5.0 × 10 5 ± 1.3 × 10 5 cells/µl and σ p = 0.130 ± 0.002.

Priors
Burst size, ω, is estimated to be 6-8 merozoites per pRBC [3].We thus take its prior as N T (7, 0.5).The rate of upregulation of erythropoiesis, θ, is estimated to be about 0.2 day −1 in [4] (this previous study uses different parasite strains) and 0.4 day −1 in [5].We thus take its prior to be N T (0.28, 0.1).The delay prior to upregulation, δ, is estimated to be about 2.5 days (95% CI: 0-5 days) after anaemia is chemically induced in mice [2].In malaria infections, anaemia begins from about day 5.We thus take a wide prior on δ of U(0, 12).The erythropoietic time lag, τ , is estimated to be between 1-3 days [5,6] but with significant variation between mice.We take a prior of U(0, 6) to cover this range.The normal RBC density, K, can be approximated from the first few days of the experiment using all mice.The mean density is 1.06 × 10 7 cells/µl and the standard deviation is 0.10 × 10 7 cells/µl.We thus take a prior of N T (1.06 × 10 7 , 0.1 × 10 7 ).Each mouse was infected with 10 5 parasites.Mice have roughly 3 ml of blood, therefore the initial parasite density was about 33 parasites/µl, or approximately 1.5 on a log-scale.Given uncertainty about how many parasites enter the circulation from the peritonea into which they are innoculated, we take a broad prior on log(P 0 ) of N T (1.5, 0.5).The death rate of multiply-parasitised RBCs, d 2 , is unknown.We therefore assume a prior on α of U(0, 1).The natural loss rate of merozoites µ is unknown.However, μ must be of the same order of magnitude as the RBC density, N , to have an effect on the dynamics (see Equation 6in the main text).Given that the maximum of N is K, which is of order 10 7 , we assume a broad prior on μ of U(0, 2 × 10 7 ).
There are no estimates in the literature for the immune response parameters, and we therefore assumed broad priors.For merozoite clearance to have an effect on the dynamics, c m must be of the order of K (see Equation 6in the main text).We therefore take a prior of N T (10 7 , 10 7 ).For pRBC and uRBC clearance to have an effect on the dynamics, c p and c u must be on the order of λ, the average number of parasites per RBC (see Equations 16 and 17 in the main text).Near to the peak of parasite density, this is of the order 1.We therefore assume priors on c p and c m of N T (1, 1).The start times, rise times and durations of the immune responses could be anywhere between the start and end of the experiment.We therefore take their priors as U(0, 19).Here T N is the number of data points for RBC density, N data,j is the measured RBC density at time point j, and N model,j (Φ) is the model solution of RBC density at time point j.Similarly, T p is the number of data points for parasite density, p data,j is the measured parasite density at time point j, and p model,j (Φ) is the model solution of parasite density at time point j.

Evolutionary Monte Carlo method
We use the Evolutionary Monte Carlo (EMC) method (without crossover) to sample posterior densities [7,8].Two significant advantages of this method are its excellent characterisation of multi-modal posterior densities and reliable estimates of marginalised likelihoods of the data, and thus Bayes' factors [9,10].
In the EMC method there are B tempered densities with inverse temperatures t 1 , . . ., t B , which form a ladder with the ordering t 1 > . . .> t B .A challenge of EMC is to find the optimal form of this ladder and the number of densities B such that the minimum state-exchange rate between neighbouring densities remains high (above 50% [7]).Below this value, marginalised likelihoods become inaccurate and imprecise.For this data set we have found that B = 40 and a temperature ladder of the form t i = t (i−1) A 2 for i = 3, . . ., B − 2, where A = ln(ln(t B−1 )/ ln(t 2 ))/ ln(B − 2), t 1 = 1, t 2 = 0.9, t B−1 = 10 −4 and t B = 0 gives satisfactory results.
We block-update all model parameters using the Metropolis algorithm.The transition distributions of the algorithm are multivariate normal with mean 0 and covariance matrices S i , for i = 1, . . ., B. The covariance matrices are estimated from the empirical covariance matrices of the Markov chains during an adaptive McMC phase [11].The adaptive algorithm works as follows: 1: Set off-diagonal elements of S i to 0; set diagonal elements to 1/500th the scale of the prior densities (the densities at the initial parameter guesses).
2: Every 100 iterations calculate the empirical S i matrices and their Cholesky decompositions.Tune their scales (by factors of 0.5 and 2) to keep the chains' acceptance rates above 20% and below 40% respectively [12,13].Exit when 500 unique samples have been collected and the acceptance rates have stabilised between 20% and 40%, i.e., have not moved out of this range over 200 iterations.
Four Markov chains are run simultaneously and convergence to the posterior is assessed using the Gelman-Rubin statistic on each parameter [14].We use 2,000,000 iterations of a non-adaptive McMC phase; the first 50% are discarded as burn-in and the the rest thinned to 1 in 500, giving 2,000 samples on which to base inferences.
Models are assessed by three methods: calculation of residual deviance, visual inspection of fits to data, and plotting overlaid standardised residuals.Visual inspection of fits highlights gross inadequacies within and across replicates, but is not particularly sensitive.Standardised residuals provide a more sensitive assessment of fit.Poor fits exhibit outliers and serial correlation.By overlaying plots, systematic biases at single time points are also revealed (see below).
The EMC algorithm, the models and the various diagnostics and statistics are implemented in the C programming language on the Linux operating system, and complied using the Intel C/C++ compiler.Diagnostic visualisation is coded in the Python programming language coupled to the GRACE graphics plotting package.We use the OpenMP Application Program Interface.

Assessment of model adequacy
We assess a model's ability to predict the data using marginal standardised residuals.For data point j the marginal standardised residual is: where x is either N for RBC density or p for parasite density.
For the true model, standardised residuals are approximately normally distributed with mean 0 and standard deviation 1.Although they are not strictly independent and identically distributed, we can treat them as such for the purposes of model assessment, which by its nature is subjective.Poor fits to data therefore exhibit outliers (larger than 3 standard deviations from 0) and serially correlated residuals.
We can also look for systematic biases in fits across mice because their RBC and parasite dynamics are sufficiently similar over time.We do this by overlaying the standardised residual plots for each mouse and plotting the mean of the standardised residuals at each time point.For the true model the mean of n standardised residuals is approximately normally distributed, with mean 0 and standard deviation 1/ √ n.We can therefore examine, for each time point, if the mean of the standardised residuals of n mice lies within the 95% prediction interval (−1.96/ √ n, 1.96/ √ n).The interval is Bonferonni corrected to adjust for multiple time points.Mean standardised residuals that lie above (below) the interval suggest systematic underestimation (overestimation) of the data across all mice.

Model comparison
Model fits are compared using Bayes' factors, which impose a penalty for additional parameters [15].Each model, Π j , is compared to the baseline, Π 0 .Bayes' factors are calculated by evaluating the marginal likelihood of each model, Pr(data|model).The Bayes' factor for model Π 0 over model Π j is: (Bayes' factor) j = Pr(data|Π 0 ) Pr(data|Π j ) (S.18)This is more conveniently expressed in deciBans (tenths of a power of 10): (dB) j = 10(log Pr(data|Π 0 ) − log Pr(data|Π j )) (S.19) Positive deciBan values imply that Π 0 explains the data better than Π j .Negative deciBan values imply that Π j explains the data better than Π 0 [16].We adopt the scale of interpretation given by Jeffreys [15], which is reproduced in the main text.
Because we use EMC to estimate marginal likelihoods, there is inevitably sampling error in their estimates.This feeds through into estimates of deciBans.We estimated the error in deciBans by running 100 simulations of the models and calculating the standard deviation of the marginal likelihoods for each mouse.

Assessment of Markov chain convergence
Figure S2 shows the Gelman-Rubin statistics [14] for each parameter sorted by mouse (top panel) and by parameter (bottom panel).Almost all values are close to 1, suggesting good convergence of the four Markov chains to the posterior.

ForN 2 +
each mouse, the log-likelihood of a model with parameter vector Φ given the data is, up to a constant of proportionality ln L(Φ) data,j − N model,j (Φ) σ N Tp j=1 log p data,j − log p model,j (Φ)