**DOI:**10.1128/AAC.01463-12

## ABSTRACT

Murine models are used to study erythrocytic stages of malaria infection, because parasite morphology and development are comparable to those in human malaria infections. Mechanism-based pharmacokinetic-pharmacodynamic (PK-PD) models for antimalarials are scarce, despite their potential to optimize antimalarial combination therapy. The aim of this study was to develop a mechanism-based growth model (MBGM) for Plasmodium berghei and then characterize the parasiticidal effect of dihydroartemisinin (DHA) in murine malaria (MBGM-PK-PD). Stage-specific (ring, early trophozoite, late trophozoite, and schizont) parasite density data from Swiss mice inoculated with Plasmodium berghei were used for model development in S-ADAPT. A single dose of intraperitoneal DHA (10 to 100 mg/kg) or vehicle was administered 56 h postinoculation. The MBGM explicitly reflected all four erythrocytic stages of the 24-hour P. berghei life cycle. Merozoite invasion of erythrocytes was described by a first-order process that declined with increasing parasitemia. An efflux pathway with subsequent return was additionally required to describe the schizont data, thus representing parasite sequestration or trapping in the microvasculature, with a return to circulation. A 1-compartment model with zero-order absorption described the PK of DHA, with an estimated clearance and distribution volume of 1.95 liters h^{−1} and 0.851 liter, respectively. Parasite killing was described by a turnover model, with DHA inhibiting the production of physiological intermediates (IC_{50}, 1.46 ng/ml). Overall, the MBGM-PK-PD described the rise in parasitemia, the nadir following DHA dosing, and subsequent parasite resurgence. This novel model is a promising tool for studying malaria infections, identifying the stage specificity of antimalarials, and providing insight into antimalarial treatment strategies.

## INTRODUCTION

Malaria remains a major health concern in countries where the disease is endemic, and optimized dosing strategies are therefore required for effective treatment (1, 2). The artemisinins are currently first-line malaria treatment, due to their more rapid killing of parasites than any other available antimalarial drug (3, 4). A further benefit is that the artemisinins are active against parasites resistant to traditional antimalarials, although emerging resistance has been reported recently (5–7). However, their short elimination half-lives *in vivo* make them inadequate for use as monotherapy, due to recrudescence of infections after a standard 3-day dosing regimen (8). For this reason, artemisinin combination therapy (ACT) with longer-acting partner drugs has been recommended for the treatment of malaria (1). Nevertheless, ACT dosing regimens have been challenging to develop and have at best been empirically derived, primarily due to the lack of pharmacokinetic (PK) and pharmacodynamic (PD) data to guide dose optimization (9).

Murine models are useful tools for investigating malaria infection, because parasite morphology and development in mice are comparable to those in human malaria infections (10). The limitations of murine models are well documented and have repercussions for pathogenesis (including cerebral malaria), immunity, and vaccine testing (11, 12). Nevertheless, murine models have a well-established role in drug discovery and chemotherapy research programs and may support the development of future mathematical models for patients (11, 12). Plasmodium berghei belongs to a group of Plasmodium species that infect rodents and is extensively used in preclinical drug testing (13). For translational PK-PD, the main difference is the time to complete asexual schizogony, which takes 18 to 24 h for P. berghei in rodents and 48 to 72 h for Plasmodium falciparum in humans (13). The key advantage of using P. berghei in murine models is that all erythrocytic stages of the parasite life cycle are observable in thin blood films and are easily differentiated by light microscopy (10). Thus, murine malaria studies can generate robust PK-PD data, as has been reported for dihydroartemisinin (DHA) (14, 15), chloroquine (16), and piperaquine (17), which could be used to optimize dosing. However, traditional animal studies do not allow the evaluation of all doses and are unable to dissect immune elimination of parasite or the contributions of individual antimalarials in combination therapy. Thus, the sole use of animal models to inform ACT dosing has limitations, and mechanism-based mathematical models offer significant advantages in the development of antimalarial dosage regimens.

Mechanism-based models describe the time course of drug effects and disease, thereby maximizing the information gained from experimental data (18–20). However, mechanism-based models that account for multiple stages of the parasite life cycle are lacking (21), despite the potential for optimizing ACT strategies (9, 22–25). Murine models with P. berghei are particularly valuable in this regard, since the effects of antimalarials on specific stages of the parasite life cycle can be well resolved. Mechanism-based models of murine malaria must therefore simultaneously consider the parasite life cycle, the delayed onset of drug effect in relation to its plasma PK, and the stage specificity of antimalarials (26, 27). Once developed, these models can be used to simulate effective dosing strategies to guide the choice of antimalarial monotherapy and combination therapy in patients (9).

The purpose of this study was to develop a mechanism-based growth model (MBGM) for P. berghei in murine malaria and incorporate the parasite-killing effect of DHA (MBGM-PK-PD). We discuss the utility of the model for other antimalarials and the potential for dose optimization of ACT in humans.

## MATERIALS AND METHODS

Source of data.Parasite density data were obtained from a previously reported study that investigated P. berghei growth and the subsequent effect of DHA in malaria-infected mice (15). Briefly, male Swiss mice (6 to 7 weeks old) were infected with an inoculum of 10^{7} P. berghei parasitized erythrocytes by intraperitoneal (i.p.) injection. The mice then received a single i.p. dose of either 0 (control), 10, 30, or 100 mg kg^{−1} DHA (dissolved in 60:40 DMSO-polysorbate 80) at 56 h following inoculation (parasitemia, 2 to 5%). Peripheral blood smears for the determination of parasitemia were prepared every 12 h prior to DHA dosing, every 4 h for the first 52 h after drug treatment, and then twice daily until euthanasia (15). Parasitemia was determined, and organisms were classified as one of four erythrocytic stages (ring, early trophozoite, late trophozoite, and schizont [13]) by an experienced microscopist. To facilitate the modeling, a linear correlation was used to convert the percentage of infected erythrocytes to parasite density (the number of parasitized erythrocytes per μl of whole blood) (14). The limit of detection was 1,000 parasites/μl.

In an independent study, plasma concentration-time data were obtained following i.p. administration of DHA (100 mg kg^{−1}) to mice (28). Blood was collected by cardiac puncture at 2, 5, 10, 15, 30, 45, 60, 90, and 120 min after DHA dosing, and plasma was analyzed using a validated high-performance liquid chromatography (HPLC) method (29). The limit of quantitation (LOQ) was 90 μg/liter, and assay precision was 9%.

Mathematical modeling.Initially, a mechanism-based growth model (MBGM) for P. berghei was developed independently of that for DHA PK. Logarithmically transformed parasite density-versus-time data were analyzed by nonlinear mixed-effects modeling in S-ADAPT (version 1.57) (30), using the SADAPT-TRAN facilitator tool (31). Parasite density data from control and treatment mice were then combined and linked to the DHA PK model, to simultaneously describe parasite growth and drug killing (MBGM-PK-PD). Data below the LOQ were handled using an established likelihood-based (M3) method (32).

Model selection was based on visual inspection of diagnostic scatter plots, the objective function value (OBJ; reported as −1 × log likelihood in S-ADAPT), and the biological plausibility of parameter estimates. Statistical comparison of nested models was performed using a χ^{2} test, in which a decrease in the OBJ of 1.92 units (α = 0.05) was considered significant.

The final models were evaluated by performing a visual predictive check (VPC) (33). For this, 1,000 data sets were simulated from the final parameter estimates using the original data. The median and the 25th and 75th percentiles of simulated predictions were then computed and plotted against observed values. For the PK model only, a nonparametric bootstrap method was used to assess the uncertainty of all parameter estimates in the model (34). The 2.5th, 50th, and 97.5th percentiles were calculated from the empirical posterior distribution of 1,000 bootstrap replicates. However, this bootstrapping method was not performed for the MBGM-PK-PD, because of the long estimation time (1 to 2 days) associated with a single run of the full model. All reported parameter estimates and VPCs are from the simultaneous modeling of parasite growth and DHA PK-PD.

Mechanism-based growth model (MBGM).A model reflecting all four erythrocytic stages (ring, early trophozoite, late trophozoite, and schizont) of the parasite life cycle was initially developed to describe P. berghei growth. The structure of the final growth model is illustrated in Fig. 1.

Changes in the numbers of parasitized erythrocytes (*A*_{ip_par}) after i.p. inoculation and the associated initial condition (IC) were calculated as follows:
*k*_{absp} is the first-order rate constant for absorption of parasitized erythrocytes and *F*_{total} is the sum of bioavailability for rings (*F*_{rng}), early trophozoites (*F*_{etrp}), and late trophozoites (*F*_{ltrp}). The number of schizonts in the initial inoculum (dose_{par} = 10^{7} parasites) was assumed to be negligible, since this parasite stage is a small component of the inoculum (<5%) and is less likely to survive the centrifugation and passaging process due to erythrocyte fragility (35, 36).

The rates of merozoite (*A*_{mer}) formation from schizonts (*A*_{sch}), their subsequent infection of erythrocytes to form rings, and loss of noninfecting merozoites are calculated as follows:
*k*_{sch} is the first-order rate constant for schizont transfer. Infection of erythrocytes (quantified as observed rings) was described by a first-order rate constant (*k*_{inf}), which declined with increasing parasite density (merozoites in the model). In this term, *A*_{mer50} is the amount of merozoite at half-maximal infection rate and infhill is the Hill coefficient. The elimination of noninfecting merozoites is represented by the first-order rate constant *k*_{minj}.

Changes in the numbers of rings (*A*_{rng}), early trophozoites (*A*_{etrp}), late trophozoites (*A*_{ltrp}), and schizonts (*A*_{sch}) through one life cycle is described by equations 3 to 6:
where *k*_{rng}, *k*_{etrp}, *k*_{ltrp}, and *k*_{sch} are the first-order rate constants for parasite maturation through the ring, early trophozoite, late trophozoite, and schizont stages, respectively. First-order rate constants for elimination by the host immune system of injured late trophozoites (*k*_{ltinj}) and schizonts (*k*_{sinj}) were also included. The last two terms in equation 6 describe saturable efflux of schizonts with first-order return to the central circulation. In the absence of this pathway, schizont amounts were drastically underpredicted at late times (>100 h) after parasite inoculation. The inclusion of these two terms resolved this underprediction, where *V*_{max_eff_sch} is the maximal rate of schizont efflux, *A*_{sch50} is the schizont amount at half-maximal efflux, and *k*_{trans} is the first-order rate constant for return from 3 transit compartments. Schizont transfer through the 3 transit compartments (tr1 to tr3) is calculated as follows:

All first-order rate constants (*k*) were parameterized and estimated as mean times (1/*k*), to facilitate the interpretation of parameter estimates. The parasite amount at each time was scaled to the apparent volume into which it distributes (*V*_{rbc}), to obtain a prediction in parasite density units (parasites/μl) (for example, the predicted concentration of rings equals *A*_{rng}/*V*_{rbc}).

DHA PK and PD model.One- and two-compartment models were tested to describe the disposition of DHA in mouse plasma. Drug input from the i.p. dosing compartment was modeled as zero- or first-order absorption. The final PK model was then linked to the above structural growth model (MBGM) to simultaneously model parasite growth and DHA killing (MBGM-PK-PD).

Drug killing of parasites was described by a turnover model, to account for the delay between the DHA dose and parasite death. In this model, DHA inhibited the production of a hypothetical physiological intermediate. Changes in the amount of the physiological intermediate (*A*_{phys}) were modeled using a turnover compartment, arbitrarily initialized at a steady-state value of 1. In this system, the effect of DHA is:

where *k*_{in} is the zero-order rate constant for *A*_{phys} production, IC_{50} is the DHA concentration at half-maximal *A*_{phys} inhibition, and *k*_{out} is the first-order rate constant for *A*_{phys} loss. At steady state, the rate of production of the physiological intermediate equals that of its loss, resulting in a unity ratio of rate constants (*k*_{in}/k_{out} = 1). Thus, the apparent lack in *A*_{phys} is

Equation 11 was then incorporated into a “kill” function that stimulated the DHA elimination of parasite (*A*_{par}) as follows:
*k*_{DHA} is the first-order rate constant for parasite elimination by DHA and kill_{50} is the amount of physiological intermediates at half-maximal killing. For each of the four erythrocytic stages, *A*_{par} corresponded to the respective amount of parasite in that compartment. The kill function is incorporated in large parentheses and ensures zero killing of parasite in the absence of DHA, where *A*_{phys} is arbitrarily defined at a steady state of 1 (lack_{phys} = 0). A schematic representation of the DHA PK-PD model is provided in Fig. 2.

Between-subject variability and residual error models.The between-subject variability (BSV) was assumed to follow a lognormal distribution and was estimated for all parameters. For convenient interpretation, BSV is reported as the square root of the estimated variance, as this approximates the apparent coefficient of variation (CV) of a normal distribution on log scale.

Residual unexplained variability (RUV) for logarithmically transformed parasite density data was estimated using exponential random error. For the DHA PK model, exponential and additive components were fixed to assay precision (9%) and LOQ (90 μg/liter), respectively. This allowed calculation of BSV, as only a single observation was obtained from each mouse due to sampling by terminal bleeds (28).

## RESULTS

P. berghei growth model (MBGM).Parasite density-time data from seven untreated malaria-infected mice were used to develop the growth model for P. berghei. The final model explicitly incorporated all four maturation stages of the parasite life cycle, as well as release of merozoites and their subsequent infection of erythrocytes. Absorption of parasitized erythrocytes was rapid and first order (*T*_{absp} = 0.0522 h; BSV, 22%), with the bioavailability of rings (*F*_{rng} = 0.568) far exceeding that for early (*F*_{etrp} = 0.101) and late (*F*_{ltrp} = 0.0289) trophozoites. The apparent volume into which the parasite distributes (*V*_{rbc}) was estimated at 1180 μl (BSV, 18%), which approximates the red blood cell volume in mice.

The estimated mean time to complete one life cycle of asexual schizogony for P. berghei was 24.6 h. Merozoite infection of erythrocytes to rings was described by a first-order process (*T*_{inf} = 0.0446 h; BSV 12%) that declined with increasing parasitemia. This decline function significantly improved the model (ΔOBJ, −31.8) for merozoite invasion of erythrocytes and the corresponding formation of rings. Immune elimination of merozoites (*T*_{minj} = 5.99 h) was faster than that for schizonts (*T*_{sinj} = 115 h), which in turn was more rapid than the estimate for late trophozoite loss (*T*_{ltinj} = 2,390 h). The inclusion of similar first-order elimination pathways for rings and early trophozoites did not significantly improve the structural growth model (ΔOBJ, −0.725), and these pathways were therefore not included in the final model. Saturable schizont efflux with first-order return through 3 transit compartments was included to resolve the drastic underprediction of schizont amounts at late times (>100 h) after inoculation. Statistical inferiority was obtained when this route was excluded (ΔOBJ, −108) from the MBGM, or when schizont efflux was first-order (ΔOBJ, −65.4), suggesting that inclusion of these processes significantly improved the final model (Fig. 1). This pathway represents parasite sequestration or trapping in the microvasculature, with subsequent return to the circulation. The incorporation of this efflux return pathway to other stages of the P. berghei life cycle did not improve the model, and it was therefore not included.

DHA PK model.Plasma concentration-time data from 65 healthy and 33 malaria-infected mice were used to independently develop the PK model for DHA. Drug absorption was fitted by a zero-order process, with an estimated duration of 0.224 h (BSV 48%). The disposition of DHA in plasma was best described by a one-compartment model, which included estimation of BSV on all parameters. Addition of a second compartment did not significantly improve the disposition model (ΔOBJ, −0.128), and other structural models were therefore not tested. In malaria-infected mice, the estimated clearance of DHA was 1.95 liters h^{−1}, which was approximately 1.3 times faster than that in healthy mice (CL_{healthy} = 1.48 liters h^{−1}). Similarly, the estimated volume of distribution (*V*_{malaria} = 0.851 liter) was 1.6 times larger in malaria-infected mice than in healthy controls (*V*_{healthy} = 0.535 liter). Table 1 provides a summary of parameter estimates from the final PK model.

DHA PK-PD model (MBGM-PK-PD).Parasite density data from 7 untreated and 20 treated mice were combined to simultaneously model P. berghei growth and DHA killing. The concentration-effect relationship used for driving DHA elimination of parasite was defined by linking the developed structural growth model (MBGM) to that for drug PK in malaria-infected mice. A delayed effects model described DHA PK-PD and was statistically superior (ΔOBJ, −148) to a corresponding model in which drug killing was immediate. The delay in onset of parasite killing in relation to maximal drug concentration was best described using a turnover model, with DHA inhibiting the production of hypothetical physiological intermediates (IC_{50}, 1.46 ng/ml; BSV, 37%). This turnover was incorporated into a kill function that stimulated parasite elimination (equation 12 and Fig. 2), and was necessary for all stages of the P. berghei erythrocytic cycle. While stage-specific differences in the rates of parasite killing were evaluated, a single first-order rate constant (mean time for parasite elimination by DHA [*T*_{DHA}] = 5.66 h; BSV, 17%) adequately described the elimination of all erythrocytic maturation stages. Goodness-of-fit plots from the simultaneous modeling of P. berghei growth and DHA killing (MBGM-PK-PD) suggested an acceptable performance of the final model (see Fig. S1 to S4 in the supplemental material), which simultaneously fitted four dependent variables. The parameter estimates for the final model are provided in Table 2.

Model evaluation.The VPCs for DHA in plasma for healthy and malaria-infected mice are presented in Fig. 3. Representative VPCs at the intermediate dose of 30 mg kg^{−1} DHA are also provided for the P. berghei growth and parasite killing model (Fig. 4 and 5). These plots illustrate that 20 to 25% of the data lie outside the predicted 25th and 75th percentiles, with the majority of observations being evenly distributed around the median. In addition, all of the prediction percentiles broadly match their corresponding observed percentiles. These VPCs therefore demonstrated a reasonable predictive performance of the models, with prediction percentiles describing the central tendency and variability in the data. For all PK parameters from the bootstrap analysis, the median values were similar to that in the final model (Table 1).

## DISCUSSION

In the present study, we developed a novel, mechanism-based population model that describes the rise in parasitemia following murine inoculation with P. berghei, the nadir after DHA dosing, and subsequent parasite resurgence. This model explicitly incorporated four distinct erythrocytic stages of P. berghei, merozoite release from schizonts and subsequent erythrocyte invasion, as well as the PD effect of parasite killing by DHA. While the parasiticidal effect of DHA was demonstrated for all four erythrocytic stages, our model offers the advantage of testing antimalarial stage specificity on the parasite life cycle.

The structural growth model for P. berghei (MBGM) in murine malaria demonstrated biological and physiological plausibility for a number of reasons. First, parasites and/or parasitized erythrocytes are rapidly absorbed after i.p. inoculation (*T*_{absp} = 0.0522 h), possibly via the lymphatic system (37–39). Furthermore, the dominance in bioavailability of rings (*F*_{rng} = 0.568) is expected, since these are most likely to survive the inoculation process relative to early or late trophozoites (*F* < 0.11). The latter are more susceptible to physical and chemical breakage, due to their larger size and maturation in the parasite life cycle (14). In addition, the apparent red blood cell volume into which the parasite distributes (*V*_{rbc} = 1180 μl) is consistent with the theoretical erythrocyte volume in mice (40). However, a possible limitation is that this *V*_{rbc} estimate does not include changes in red blood cell homeostasis, which may occur with malaria pathophysiology. Second, the high rate of merozoite infection (*T*_{inf} = 0.0446 h) and elimination of noninfecting merozoites is likely, given the known mechanism of invasion of erythrocytes (41). A decline in the rate of infection with increasing merozoite is also appropriate, since the pool of uninfected erythrocytes will be depleted with increasing parasite density (14, 26). Third, a faster elimination of schizonts (*T*_{sinj} = 115 h) is expected relative to late trophozoites (*T*_{ltinj} = 2,390 h), as schizonts are more readily detected by the host reticuloendothelial system (13, 26). In contrast, immune elimination of rings and early trophozoites was not supported by the model for the present data set in mice, since these are considered pliable when passing through the spleen (15, 42). Finally, the sequestration or trapping of parasites in the microvasculature is a known phenomenon and is a potential mechanism by which the parasite avoids immune detection (43, 44).

Inclusion of schizont efflux into the model was a relatively minor quantitative component, because schizonts comprised an average of 5 to 7% of the total parasite density. Moreover, while sequestration of P. berghei is considered negligible in Swiss mice, parasite trapping could explain our finding that schizont transfer through transit compartments improved the growth model (15, 43, 45). However, in studies where sequestration is a more significant factor, schizont efflux or a similar strategy may be an essential component of the model. Furthermore, in the DHA PD model, this schizont efflux pathway could potentially include parasite dormancy, which is induced by the artemisinin drugs (46–48).

Population PK analysis demonstrated that a 1-compartment model adequately described the disposition of DHA in mice. The estimated clearance (1.95 liters h^{−1}) and volume of distribution (*V* = 0.851 liter) in malaria-infected mice were consistent with the original study (28), hence the use of population PK parameters to develop the MBGM-PK-PD model. The time course of DHA killing of parasite was described by a turnover model, with DHA inhibiting the production of hypothetical physiological intermediates. This concentration-effect relationship is indirect and is a well-established tool for describing the delay between drug PK and PD (18). Although there are several proposed mechanisms of action for artemisinin drugs, the DHA turnover model utilized in this model encompasses all probable theories associated with its antimalarial activity (49–51). However, the IC_{50} (1.46 ng/ml) was difficult to estimate in this study, because the concentrations of DHA achieved after dosing greatly exceeded its independently determined IC_{50} (1 to 4 nM) *in vitro* (52, 53). Given its low *in vitro* IC_{50} estimate and the high plasma protein binding (90% [54]), further experimentation at lower doses is required to elucidate the IC_{50} of DHA *in vivo*. Nonetheless, the kill function for DHA was required for all stages of parasite elimination, consistent with its broad stage specificity (52).

Evaluation of the final mechanism-based model by VPC (Fig. 4 and 5), as well as the goodness-of-fit plots (see Fig. S1 to S4 in the supplemental material), demonstrated a reasonable predictive performance of the model. However, some model misspecification is apparent at early times (24 h postinoculation), with an underprediction of measured parasite density data. This limitation is likely related to the lack of data to accurately identify the kinetics of parasite absorption after inoculation, which, combined with the asynchronicity of infection, may explain the higher estimated multiplication factor (MF = 22), compared to reports of 6 to 20 (55) and 12 to 18 merozoites per schizont (13) for P. berghei. Nonetheless, while low parasite counts are difficult to obtain because they are close to the limit of detection, the estimated bioavailability and first-order absorption of parasite appear to be plausible. Further model misspecification is apparent at 60 to 65 h and 80 h (nadir) for rings and early trophozoites (Fig. 5; also, see Fig. S1 and S2 in the supplemental material). A probable explanation is that at all administered doses (10 to 30 mg/kg), the concentrations of DHA achieved exceeded the concentration required for near-maximal parasite killing. As a consequence, the estimation of IC_{50} and *k*_{DHA} was challenging in this study and could be better carried out with the administration of DHA at lower (and potentially more clinically relevant) doses. These experiments would also assist with specifying the schizont PD model, which appears to overpredict parasite killing where the data are at the limit of parasite density detection (log_{10} of 3.0).

Our mechanism-based model for P. berghei growth and DHA PD is potentially applicable to the development of ACT dosing guidelines in malaria chemotherapy. Its structure is adaptable to the quantitative evaluation of other parasite strains and antimalarials for both preclinical and human studies. However, for clinical studies, the counts of circulating parasites must be considered because, unlike P. berghei, not all of the erythrocytic stages of P. falciparum can be measured. Other differences in malaria infections between mice and humans also may need to be considered for translation of the present model to humans. Nevertheless, the proposed model could be a valuable tool in the development of antimalarial treatment strategies, because there is scope to simulate and design effective single-agent, multiple-dose or combination therapy studies. A further benefit is the capacity of the model to dissect the contributions of individual antimalarials in combination therapy, as well as parasite elimination by the host immune system. Without a modeling approach, it is essentially impossible to obtain this information from traditional preclinical studies or from human clinical trials. Therefore, this model can be a cost-effective tool to expedite the translation of antimalarial drug development and to rationally design effective ACT regimens.

## ACKNOWLEDGMENTS

B.R.M. is the recipient of a National Health and Medical Research Council (Australia) Early Career Fellowship; CJ Martin Overseas Biomedical Fellowship (APP1036951). J.B.B. is the recipient of a DECRA fellowship (DE120103084) from the Australian Research Council.

## FOOTNOTES

- Received 16 July 2012.
- Returned for modification 16 August 2012.
- Accepted 3 November 2012.
- Accepted manuscript posted online 12 November 2012.
Supplemental material for this article may be found at http://dx.doi.org/10.1128/AAC.01463-12.

- Copyright © 2013, American Society for Microbiology. All Rights Reserved.