All animal experiments were performed according to European regulations concerning FELASA category B and GV-SOLAS standard guidelines. Animal experiments were approved by German authorities (Regierungspräsidium Karlsruhe, Germany), § 8 Abs. 1 Tierschutzgesetz (TierSchG) under the license G-260/12 and were performed according to National and European regulations. For all experiments, female C57BL/6 mice (6- to 8-week-old) were purchased from Janvier laboratories, France. All mice were kept under specified pathogen-free (SPF) conditions within the animal facility at Heidelberg University (IBF).

In the first set of experiments, C57BL/6 mice were intravenously infected with 10⁶ infected red blood cells (iRBC) taken from mice infected either with wild-type PbGFP Luc_con (P. berghei line 676m1c11) (WT), a GFP-luciferase transgenic derivative of P. berghei ANKA (Franke-Fayard et al., 2005), or the mutant PbmaLS_05 (–) parasites (KO) generated in the wild-type PbGFP Luc_con strain (Fernandes et al., submitted manuscript). An additional group of age-matched mice was left uninfected and treated as naïve controls. Daily blood samples of 10 μl were taken from all mice from the day of infection until day 4 post infection (p.i.). The total red blood cell count and reticulocyte percentage were measured using a Coulter counter and FACS analysis of CD71 (CD71-PE, eBioscience, Clone R17217) labeled reticulocytes, respectively. Parasitemia was determined by FACS analysis of GFP positive infected red blood cells. A sketch of the experimental protocol is shown in Figure 1A. Mice were sacrificed at day 5 p.i., when mice infected with WT parasites showed first symptoms of ECM.

A second set of mice were pretreated with two doses of phenylhydrazine (PHZ, 40 mg/kg), on two consecutive days prior to infection with 10⁶ iRBC using the same groups of mice as before. Again, daily blood samples of 10 μl were taken from each mouse and analyzed up to day 5 p.i. before sacrificing the mice on day 6 p.i.

To describe the blood stage-infection dynamics of the murine malaria parasite accounting for RBC age, we used a mathematical model for erythropoiesis as described before (Mackey, 1997). The age-structured model follows the population density of RBCs of age τ at time t based on a system of coupled ordinary differential equations that breaks the age ranges of RBC into n = τ_RBC/h different compartments with h being the compartment size and τ_RBC the maximal lifespan of RBCs. The concentration of RBCs within each compartment is denoted by x_i(t), i = 1,… n. New RBCs are constantly produced by the bone marrow that enter the population of RBCs after a delay T, with the actual influx at each time point determined by a Hill-function dependent on the maximal production rate of RBCs in the bone marrow, F₀, and the concentration of RBCs at time t-T, X(t-T). Mathematically, the model is then described by the following system of ordinary differential equations:

Hereby, the parameter θ describes the concentration of RBC where the production rate is half of the maximum and k the Hill-coefficient (Mackey, 1997). In addition, we also assumed that in each compartment x_i RBCs are lost by an age-independent loss-rate 1/τ_RBC to have at least 85% of RBC lost until their assumed maximal lifespan τ_RBC. Equations (1–3) represent a mean-field approximation of the originally developed system relying on partial differential equations, thereby transforming assumed fixed, constant lifespans of RBC into gamma-distributed lifetimes (Mackey, 1997; Antia et al., 2008).

This basic model for erythropoiesis is then extended to account for malaria blood-stage infection as done previously (McQueen and McKenzie, 2004; Antia et al., 2008; Figure 1B). Uninfected RBCs can get infected by free merozoites, z, at a rate β(τ), which is dependent on the age-preference of the infecting parasite strain. Each infected RBC releases a number of merozoites, m, by bursting after having reached an infection maturation time, t_m. In addition, free merozoites are assumed to have an average lifetime of 1/d_m. As for uninfected RBC, the concentration of infected cells, Y(t), is broken down into g = t_m/h different age compartments, y_i(t), i = 1,…,g leading to a system of coupled ordinary differential equations with a gamma-distributed maturation time with mean t_m. The basic model for erythropoiesis (Equations 1–3) is then extended to:

Hereby, z(t) describes the concentration of merozoites at time t and RF the so called reticulocyte factor, i.e., the fold-change in infectivity of the parasite for reticulocytes compared to the general infection rate assumed for normocytes, β₀ (see Cromer et al., 2006). The parameter τ_Reti defines the maturation time of reticulocytes into normocytes.

In order to compare parasite strains with possible different values for the infection rate, β₀, and the reticulocyte factor RF, we calculated an average infectivity β, which is defined as the infection rate of a single merozoite when placed into the erythropoietic system at the initiation of infection. In a naïve mouse, on average 5.8% of the RBC are reticulocytes, thus the average infectivity is calculated by β = β₀(0.058RF + 0.942).

Besides the reticulocyte factor, RF, the reticulocyte preference, RP, is calculated based on the ratio between the percentage of infected reticulocytes (relative to all reticulocytes) and the percentage of infected normocytes (relative to all normocytes). Thus, if R and I_R define the concentration of reticulocytes and infected reticulocytes, respectively, and N and I_N the corresponding concentrations for normocytes, the reticulocyte preference is calculated by RP = (I_R/R)/(I_N/N). In contrast to the reticulocyte factor, the reticulocyte preference can be directly calculated from experimental measurements.

Treatment with Phenylhydrazine (PHZ) is used for experimental induction of anemia in animal models to study hemolytic anemia or anemia caused by destruction or removal of RBCs from the bloodstream (Berger, 2007). Previous studies developed mathematical models to determine and quantify the effect of PHZ on the RBC age distribution and altered erythropoiesis (Savill et al., 2009). However, these models were inadequate to describe our experimental data suggesting that they incompletely addressed the effects of PHZ. To this end, we tested several different known hypotheses for the effect of PHZ on erythropoiesis (Jain and Hochstein, 1980; Berger, 2007; Savill et al., 2009; Moreau et al., 2012) by fitting them to the data of the PHZ-control group (see Supplementary Material Text S3). The models best explaining the experimental data included the following drug effects: (i) Treatment by PHZ leads to instantaneous lyses of a fraction ρ(τ) of RBCs at the time of treatment, t_p. Hereby, the effect of lysis depends on the age of the RBC, τ, with normocytes being more strongly affected than reticulocytes (Jain and Hochstein, 1980). (ii) An additional influx of reticulocytes from extra medullary sites is considered at a constant rate N_p with a time-delay T_p after the initiation of treatment to account for stress-induced erythropoiesis. Under severe anemia, such as that induced by PHZ-treatment, extra-medullary sites of erythropoiesis such as the spleen and liver are observed to show an increased contribution of RBCs to circulation (Spivak et al., 1973; Ploemacher et al., 1977; Kim, 2010). Thus, under PHZ-treatment, Equations (1, 2) describing RBC turnover are changed as follows:

Hereby, ρ₀ defines the fraction of normocytes lysed by PHZ and γ represents the relative comparison of this fraction for reticulocytes. In addition, I(t = T_p) defines the Indicator function, i.e., with I(t = T_p) = 1 if t = T_p and 0 otherwise. A sketch of the effects of PHZ treatment on erythropoiesis is shown in Figure 4A. A detailed derivation of the model can be found in the Supplementary Material. During infection, we assume that malaria induced changes to RBC production affects both sources of novel reticulocytes, i.e., bone marrow and extra medullary sites alike.

The mathematical models described above were implemented and analyzed using the R language of statistical computing (R Development Core Team, 2017). As indicated, the age of uninfected and infected RBC was compartmentalized leading to a tractable system of coupled ordinary differential equations with gamma-distributed lifetimes and maturation times for RBC and infected cells, respectively (Antia et al., 2008). In the following we used a compartment size of 4 h.

The differential equations were solved using the deSolve package and models were fitted to the experimental data using the optim-fitting routine in R. In cases where a strong correlation between parameters hindered convergence of fitting algorithms, a parameter sweep was performed to find combinations of parameters that fit the data. Proportion data (parasitemia levels and proportion of reticulocytes) were logit- transformed to allow for normally distributed residuals. Model performance was assessed based on simultaneous fitting for all obtained measurements including RBC concentration, reticulocyte proportion and, where applicable, parasitemia. Blood stage infection dynamics of parasites were determined in a stepwise approach: Parameters describing erythropoiesis were fixed to the indicated values obtained from the naïve control group before analyzing infection dynamics (Table 1). Therefore, measurements for the infection groups, i.e., reticulocyte proportion and RBC count, were scaled relative to the naïve group data when estimating parasite infectivity. To evaluate model performance, the average residual sum of squares (aRSS) was used which is the residual sum of squares divided by the number of data points.

The 95%-confidence intervals, as well as identifiability of parameter estimates were assessed by profile likelihood analysis (Raue et al., 2009). For the measured data, we report mean and standard error.

To determine the dynamics of erythropoiesis in our experimental system, we fitted a mathematical model describing RBC production and subsequent aging (see Equations 1–3 in Materials and Methods; Mary et al., 1980; Mackey, 1997) to the observed progression of RBC concentration in uninfected mice that were sampled daily for 10 μl of blood (see Materials and Methods and Figures 1A–C). In general, bleeding leads to a decrease in the RBC concentration triggering the production of novel RBCs in the bone marrow that will enter the blood circulation after a time delay T. Thereby, the magnitude of the feedback depends on the severity of the anemia, i.e., the larger the loss of blood the larger the subsequent RBC production, which is accounted for in our model by a Hill-type function (Mackey, 1997). Assuming a maximal lifespan for RBC of τ_RBC = 40 days (Bannerman, 1983) and a Hill-coefficient of k = 7.6 (Mackey, 1997), we estimated a maximal RBC production rate in the bone marrow of F₀ = 5.95 × 10⁴ cells μl⁻¹ h⁻¹ [4.02, 6.82] with half of the maximal production rate reached at a RBC concentration of θ = 6.65 × 10⁶ cells μl⁻¹[5.28, 6.84], which is approximately 95% of the RBC concentration at steady state. Newly produced red blood cells are estimated to appear in the circulation after a lag-time of T = 2 days, testing different possible lag-times including T = 0, 1, 2, and 2.5 days. All our estimates are in agreement to parameters that have been determined previously for erythropoiesis in mice (Mary et al., 1980; Mackey, 1997; Figure 1D and Table 1).

As we were especially interested in the dynamics of reticulocytes, i.e., immature red blood cells, we compared model predictions for the proportion of different RBC age classes to the measured proportion of reticulocytes in order to determine the time these cells spend in the blood. We found that a maturation time for reticulocytes into normocytes in the blood of τ_Reti = 36 h best described our measured proportion of reticulocytes (Figure 1E), which is in agreement to previous calculations determining a maturation time for reticulocytes between 1 and 3 days (Ganzoni et al., 1969; Gronowicz et al., 1984; Wiczling and Krzyzanski, 2008). Thus, for the following analyses we assume that after appearance in the blood, a reticulocyte will take on average 1.5 days to develop into a normocyte.

In order to compare the blood-stage infection dynamics of the two Plasmodium berghei strains investigated, mice were either infected with PbANKA (WT) or PbmaLS_05 (-) (KO) infected red blood cells and sampled daily for 10 μl of blood. For both strains, we observe a substantial loss in the proportion of reticulocytes around day 3 post infection (p.i.) coinciding with an increase in parasitemia (Figure 1C). At day 4 p.i., when mice infected with WT show first signs of ECM, the parasitemia was approximately twice as high as the one measured for mice infected with the KO (0.63 ±0.05% WT compared to 0.29 ±0.03% KO) (Figure 1C).

To determine systematic differences in the infection dynamics between the two parasite strains, we extended our mathematical model describing erythropoiesis to include malaria blood-stage infection dynamics (see Equations 4–9). Hereby, RBCs get infected by merozoites at an infection rate β and infected RBC (iRBC) will release new merozoites m after a certain maturation time t_m (see Figure 1B and Materials and Methods for a detailed explanation of the extended model). Assuming the average lifespan of a merozoite of 1/d_m = 30 min (Garnham, 1966), a maturation time of an iRBC of t_m = 24 h (Cox, 1988; De Roode, 2004) and that an infected RBC releases on average m = 9 merozoites after bursting (Cox, 1988; De Roode, 2004; Reilly et al., 2007) we find that the observed loss in the proportion of reticulocytes around day 3 p.i. cannot be explained by the increased parasitemia when using the standard parameterization for erythropoiesis (Table 1). This observation is independent of the assumed infectivity of the parasite strain (Supplementary Figure S1) and is also the case if we assume that the infectivity for reticulocytes is substantially higher than for normocytes. This indicates that the reason for the observed decrease in reticulocyte proportion is not mainly due to reticulocytes being parasitized.

It is known that malarial-induced anemia causes erythropoietic suppression, starting during the early stages of infection (Villeval et al., 1990; Sexton et al., 2004; Thawani et al., 2014). By analyzing the expression levels of previously studied genes (Sexton et al., 2004), we found that the fold change in the expression of the genes most strongly associated with erythropoiesis, i.e., α-globin, β-globin major and β-1-globin, can be described by a logistic-loss function given by

where λ defines the loss-rate of gene-expression, i.e., the loss of RBC production and t₀ the time point at which half of the maximal gene expression is reached. We estimate λ = 2.22 d⁻¹ (95%-CI [1.31, 3.05]) and t₀ = 3.70 d [3.28, 4.23] (Figure 2A, Table 1). This parameterization is then used to account for malaria-induced modulation of RBC production during the analyses of blood stage infection dynamics in WT and KO infected mice.

To analyze the infectivity of both parasite strains, we fitted our extended mathematical model (Equations 4–9 with Equation 14) to the experimental data on RBC count, reticulocyte proportion and parasitemia. Additionally accounting for a modulation of RBC production due to infection (i.e., replacing F₀ by F₀F(t) with λ = 2.22 d⁻¹ and t₀ = 3.70 d in Equation 3) improves model predictions, especially regarding the substantial loss in the proportion of reticulocytes starting 3 days p.i. (compare Figure 2B and Supplementary Figure S1A).

By estimating the infectivity for each parasite strain characterized by the rate of infection (β₀) and the reticulocyte factor (RF), our analysis indicates that the WT parasites have a higher preference for infecting reticulocytes than normocytes (Figure 2B and Table 2). During this early infection phase, we estimate a more than 22-fold higher infectivity for reticulocytes than for normocytes i.e., RF > 22 (Table 2). In contrast, a similar preference for reticulocytes could not be found explicitly for the KO parasite. Here, a model assuming equal infectivities for reticulocytes and normocytes, i.e., RF = 1, performs equally well as a model that assumes a reticulocyte preference (AIC 40.7 vs. AIC 42.7). However, our time courses are too short to clearly identify such a reticulocyte preference for both parasite strains. As a high infection rate β₀ can be compensated by a small value of RF and vice versa, several combinations of β₀ and RF can explain the observed dynamics (Supplementary Figure S2).

To compare the infectivity of WT and KO parasites, we calculated an average infectivity β based on the estimates of β₀ and RF, which is defined as the infection rate of a single merozoite when placed into the erythropoietic system at the start of infection (see Materials and Methods for a detailed calculation). We find that KO parasites have a reduced average infectivity compared to WT parasites leading to less productive infections (β = 1.13 × 10⁻⁷ mz⁻¹μl⁻¹ h⁻¹ [1.08, 1.16] for WT vs. β = 0.78 × 10⁻⁷ mz⁻¹μl⁻¹ h⁻¹ [0.74, 0.83] for KO, numbers in brackets represent 95%-confidence intervals; Table 2). This reduced average infectivity can explain the slower increase in the parasitemia observed for the KO strain (Figure 2B).

If we assume that the infection rate β₀ does not differ between the two parasite strains, we find a consistently lower reticulocyte factor for the KO compared to the WT (Figure 2C and Supplementary Figure S2). Thus, our analysis indicates that KO parasites might have a particularly impaired ability to productively infect reticulocytes in comparison to the WT during the early erythrocytic stage of infection.

To elicit possible differences in reticulocyte preferences between the two parasite strains we pre-treated mice with the drug Phenylhydrazine (PHZ) before infecting them with either WT or KO parasites (Figure 3A). PHZ artificially induces anemia in mice causing peroxidation of RBC lipids leading to hemolysis and a change in RBC age distributions (Savill et al., 2009). In uninfected mice that were pre-treated with two doses of 40 mg/kg of PHZ on two consecutive days, we observe a substantial loss in the concentration of red blood cells to roughly ~1/3 of the concentration under homeostatic conditions 2 days after the last treatment with PHZ (2.5 × 10⁶ cells/μl vs. 7.6 × 10⁶ cells/μl, mean values; Figure 3B). There was a corresponding increase in the proportion of reticulocytes to up to 50% of the total RBC count at 5–6 days after the last treatment with PHZ (Figure 3B). Changes in RBC count and reticulocyte proportion of WT or KO infected mice that were pre-treated with PHZ are visible on day 5 p.i. with RBC counts reaching 4.0 ± 0.32 and 3.6 ± 0.15 × 10⁶ cells/μl for WT and KO, respectively, compared to 6.0 ± 0.29 × 10⁶ cells/μl in uninfected animals (Figure 3C). In addition, the proportion of reticulocytes in infected animals is substantially reduced compared to naïve mice; with KO infected mice still having ~3-fold higher levels than WT infected mice [42.6 ± 2.6% (naïve), 4.8 ± 1.2% (WT), 15.6 ± 1.0% (KO); Figures 3B,C]. While parasitemia levels are comparable between both infection groups (22.5 ± 1.2% vs. 21.0 ± 2.0%), the percentage of infected reticulocytes is slightly higher for WT compared to KO (24.3 ± 4.6% vs. 16.3 ± 0.8%; Figure 3C). Given these measurements, the average reticulocyte preference RP, calculated by the proportion of infected reticulocytes among reticulocytes divided by the proportion of infected normocytes among normocytes, is determined by RP_WT = 1.46 and RP_KO = 0.76, respectively. In accordance with our previous results (Figure 2C), these observations suggest that deletion of PbmaLS_05 has a potential effect on the parasite's ability to productively infect reticulocytes.

To determine if the calculated infection characteristics for WT and KO during normal erythropoietic conditions also apply after PHZ treatment, we extended our previous model to account for drug-induced changes to erythropoiesis. The exact mechanisms by which PHZ induces hemolysis and changes in the RBC age distribution have not been determined so far. Several hypotheses including faster aging of RBCs or direct lysis have been suggested and corresponding mathematical models have been proposed (Savill et al., 2009). However, these models fail to fit our experimental data, partly because they are limited to a particular PHZ treatment protocol (Savill et al., 2009). Therefore, we performed a rigorous analysis, testing several different assumptions for the effect of PHZ on erythropoiesis and their ability to explain the observed changes in total RBC count and reticulocyte dynamics in our data (see Materials and Methods and Supplementary Material Text S3 for a detailed description of the different models tested).

Our data indicate an increasing influx of reticulocytes, as well as a decreasing net-loss in normocytes after the last PHZ-treatment (Figures 3D,E). Thereby, the increased production of reticulocytes cannot solely be explained by the anemia-induced production from the bone marrow. We found that the best models explaining the effect of PHZ treatment on erythropoiesis assume (i) instantaneous hemolysis with ~35–50% of the RBC being lysed upon PHZ administration, and (ii) stress-induced erythropoiesis with an additional production of reticulocytes from different sources than the bone marrow (Figure 3F). Thereby, this additional production starts around 4.5 days after the last PHZ-treatment has been given (Table 3 and Text S3). In addition, our analysis indicates that PHZ leads to an increased death rate of RBC, reducing the average lifetime of RBC from τ_RBC~40 days to τ_RBC~8 days (see Figure 3E and Table 3). Besides a constant death rate, a linear decreasing death rate, as indicated by our calculation of the observed net-loss in normocytes (Figure 3E), could also be possible as it shows similar explanatory power for the data (Table 3). By incorporating these effects within our model, we are able to provide a modeling framework that describes PHZ-induced changes on erythropoiesis in our experimental system (Figures 4A,B).

We then simulated the pre-treatment of mice with PHZ and subsequent infection using different assumptions for parasite infectivity, β₀, and reticulocyte preference, RF, and predicted the expected levels of parasitemia and reticulocyte proportion on day 5 post infection (Figure 4C). For the KO strain, relevant parameter combinations as determined previously (Table 2) lead to reticulocyte proportions (~13%) comparable to the ones observed in the experimental data, but result in parasitemia levels of less than 1%. In contrast, combinations of RF and β₀ within the determined ranges for the WT parasite predict reticulocyte proportions that are twice as high as seen in the data (Figure 3C), and parasitemia levels that are only one-tenth of the observed level. However, neglecting previous knowledge and directly estimating RF and β₀ based on the observed parasitemia and reticulocyte proportion under PHZ treatment, both groups expect that nearly all reticulocytes are infected (80–100%), which does not agree with our data (Figure 3C). These findings indicate that there could be disease-induced changes to PHZ treatment effects that cannot be explained by a simple combination of separately determined processes of blood-stage infection kinetics, erythropoiesis and PHZ dynamics.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.