Semimechanistic Pharmacokinetic and Pharmacodynamic Modeling of Piperaquine in a Volunteer Infection Study with Plasmodium falciparum Blood-Stage Malaria

Dihydroartemisinin-piperaquine is a recommended first-line artemisinin combination therapy for Plasmodium falciparum malaria. Piperaquine is also under consideration for other antimalarial combination therapies.

pharmacokinetic properties of piperaquine were best described by a three-compartment disposition model, characterizing the triphasic disposition of the drug (difference in objective function value [DOFV] = -41.3; 2 degrees of freedom [df], compared with two-compartment disposition model). An additional disposition compartment for piperaquine did not improve the model fit (DOFV = 0). The absorption process was best described by two-transit compartments, resulting in a substantial improvement compared with a traditional first-order absorption model (DOFV = -230; 0 df). Interoccasion variability (IOV) was evaluated on absorption parameters, i.e., relative bioavailability (F) and mean absorption transit time (MTT), which improved the model fit significantly (DOFV = -37.8; 2 df). An additive error model described the data accurately, and a combined additive and proportional error model did not improve the model fit (DOFV = 0.05). Implementation of body weight, using an allometric function, improved the model (DOFV = -4.44; 0 df). No additional significant covariates were found during the stepwise covariate search. The model evaluations indicated satisfactory results, with no obvious trends in the goodness of fit plots (see Fig. S1A to D in the supplemental material). Similarly, the visual predictive check (Fig. S1E) demonstrated a good predictive performance of the model. The numerical predictive check (n = 2,000) showed that 3.2% (95% confidence interval [CI], 1.5% to 10.1%) of piperaquine observations were below the simulated 90% prediction interval and 5.1% (95% CI, 1.3% to 9.9%) were above. The bootstrap results demonstrated acceptable robustness of the piperaquine population pharmacokinetic model. The final pharmacokinetic parameter estimates of piperaquine, with precision and shrinkage are summarized in Table 2.
P. falciparum parasite growth model. A recently published pooled analysis of parasite growth data from malaria volunteer infection studies, using a larger data set including results from this study, reported a parasite life cycle of 38.8 h (17). Thus, the parasite life cycle was fixed to 38.8 h in all growth models evaluated in the current study. A total of 8.82% of the observed parasite density values were below the lower limit of detection (LOD), and these data were successfully handled using the M3 method (15). Three parasite growth models were evaluated, as described below.
Log-linear growth model and sine-wave growth model. The parasite growth rate (k G ) estimated from a log-linear growth model was 0.066 h 21 , equivalent to a parasite multiplication rate per life cycle (PMR LC ; parasite multiplication rate, equivalent to the number of new merozoites released from 1 schizont) of 12.9-fold (95% CI, 11.5 to 14.3). Adding interindividual variability on parasite life cycle improved the model fit (DOFV = -6.29; 1 df). Additional interindividual variability on k G did not improve the model fit  Table S1 and S2, respectively, in the supplemental material. The goodness of fit of the growth phase data using log-linear and sine-wave growth models showed an acceptable overall goodness of fit with no obvious systematic model misspecifications (see Fig. S2 in the supplemental material).
Semimechanistic growth model. A semimechanistic growth model mimicking the natural P. falciparum malaria parasite life cycle was developed. The rate constant of parasite maturation (k MAT ) and the rate constant of schizont rupture (k RUP ) were fixed to an arbitrary high value of 2 in order to ensure that all parasites matured from being young rings to mature parasites and that schizont rupture took place at 38.8 h (resulting in parasite multiplication). The fraction of parasite sequestration (F SQ ) was not identifiable with the data available. A sensitivity analysis of F SQ was performed by fixing this parameter, ranging from 40% to 90% with a 10% increase each time. No alteration was observed in the OFV or in model goodness-of-fit diagnostics, suggesting identifiability issues of F SQ . However, the observed growth phase data showed a recurring 10-fold   drop in total circulating parasites (P CIR ), suggesting that the fraction of parasite sequestration should be approximately 90%. Thus, the F SQ was fixed to 90% according to these observations. The onset of sequestration (T SQ ) was estimated at 29.1 h (95% CI, 27.7 to 29.7 h). Parasite growth rate was estimated as 0.0710 h 21 (95% CI, 0.0682 to 0.0771 h 21 ), equivalent to a PMR LC of 15.7-fold (95% CI, 14.1 to 19.9). Adding interindividual variability on the parasite life cycle and PMR LC improved the model fit (DOFV = -48.6; 2 df).
The developed semimechanistic growth model described the observed data well and the goodness of fit of the semimechanistic growth model demonstrated a better model fit than the log-linear growth model (DOFV = 233.2, DBIC = 222.9) and also in terms of describing the net decrease in total circulating parasite number, a consequence of parasite sequestration (Fig. S2). The semimechanistic growth model did not result in a better model fit, in terms of DOFV and DBIC, compared with the sine-wave growth model (DOFV = 36.5, DBIC = 36.5), but demonstrated a similar goodness of fit. However, the semimechanistic growth model gave some advantages by describing the maturation of parasites, sequestration of mature parasites, synchronicity of infections (i.e., synchronous parasite maturation and multiplication resulting in periodic bursts of red blood cells and the release of young parasites), and multiplication of parasites, as described in natural infections with P. falciparum. These advantages provided additional flexibility for investigating drug effects on specific stages of the parasite life cycle. Thus, the semimechanistic growth model was carried forward for further investigation. The final parameter estimates with precision and shrinkage of model parameters are presented in Table 3.
In vivo parasiticidal effect of piperaquine. The pharmacokinetic-pharmacodynamic model was developed to investigate the in vivo parasiticidal effect of piperaquine. The schematic of the final pharmacokinetic-pharmacodynamic model is shown in Fig. 2  p . c Based on nonparametric bootstrap diagnostics (n = 1,000). Parameter precision is presented as relative standard deviation (% RSE), calculated as 100 Â standard deviation/ mean value. 7.33 ng/ml), resulting in a median parasite reduction ratio per life cycle (PRR LC ) of 2.68 Â 10 2 , 2.89 Â 10 2 , and 3.02Â 10 2 when piperaquine was given as a single dose of 480, 640, and 960 mg piperaquine phosphate, respectively. The median minimum inhibitory concentration (MIC) of piperaquine derived from the 24 healthy volunteers was 2.87 ng/ml (95% CI, 1.87 to 18.3 ng/ml). The goodness-of-fit diagnostics of the final pharmacokinetic-pharmacodynamic model are presented in Fig. S3 in the supplemental material, simulation-base diagnostics (i.e., visual predictive checks) are presented in Fig. 3, and individual plots of the final pharmacokinetic-pharmacodynamic model are presented in Fig. S4 in the supplemental material. Parameter estimates from the final pharmacokinetic-pharmacodynamic model are presented in Table 3.
Simulations of clinical scenarios. In order to simulate treatment scenarios of drugresistant infections, dihydroartemisinin was used as a representative artemisinin derivative. The parasiticidal effect of dihydroartemisinin was adjusted based on the observed parasite clearance half-life from the Tracking Resistance to Artemisinin Collaboration (TRAC) study data (18). The adjusted dihydroartemisinin E max values were 0.477 h 21 and 0.216 h 21 for sensitive and resistant infections, respectively (see Fig. S5 in the supplemental material).
Various degrees of piperaquine resistance were represented by doubling (10.9 ng/ ml), tripling (16.3 ng/ml), and quadrupling (21.7 ng/ml) the estimated EC 50 of sensitive parasites (5.43 ng/ml). In a symptomatic infection (initial total circulating parasites of 10 10 ), these simulations predicted a similar probability of treatment failure in piperaquine resistance alone compared with dihydroartemisinin resistance alone (2.58% versus 1.81%). The probability of treatment failure increased with multidrug-resistant infections, resulting in 23.6% treatment failures in the presence of dihydroartemisinin resistance and a high-degree of piperaquine resistance (EC 50 of 21.7 ng/ml). A similar trend of treatment failure was predicted in asymptomatic infections (initial total circulating parasites of 10 6 ). A summary of the treatment failure probability of each drug-

Piperaquine pharmacokineƟc model
Parasite dynamic model FIG 2 The schematic of final pharmacokinetic-pharmacodynamic model describing the semimechanistic model of P. falciparum malaria parasites and the final pharmacokinetic model of piperaquine. In the piperaquine pharmacokinetic model (left), F represents relative bioavailability, k tr represents transit rate constant, CL/F represents apparent oral clearance, V C /F represents apparent central volume of distribution (PK sampling compartment), Q/F represents intercompartmental clearance from central compartment to peripheral compartment, and V P /F represents apparent peripheral volume of distribution. In the parasite dynamic model (right), circulating parasites (P1 1 P2) represent the observed parasitemia, k MAT represents first-order rate constant of parasite maturation, k SQ represents first-order rate constant of parasite sequestration, and k RUP represents first-order rate constant of schizont rupture. The killing effect of piperaquine (EFF) was described by an E max function; E max represents the maximum parasite killing rate of piperaquine, C P represents piperaquine plasma concentration, and EC 50 represents plasma concentration of piperaquine associated with half of maximum parasite killing rate.
resistant scenario is presented in Table 4. Additionally, simulations to inform on suitable pharmacokinetic-pharmacodynamic characteristics for candidate partner drugs used in triple-combination therapy demonstrated that the addition of a hypothetical drug with a PRR LC of $10 2 and total therapeutic duration of $2 weeks had a .99% probability of successful treatment. The summary of treatment failure probability with simulated hypothetical partner drug activities against multidrug-resistant infections is presented in Table 5.

DISCUSSION
In the present study, a semimechanistic growth model describing P. falciparum parasite dynamics was successfully developed. The integration of piperaquine pharmacokinetics and parasite dynamics using data from a study conducted using the IBSM model allowed the estimation of the in vivo pharmacodynamic parameters of piperaquine.
The PMR LC estimated from the log-linear growth model (12.9; 95% CI, 11.5 to 14.3), the sine-wave growth model (12.6; 95% CI, 11.1 to 14.0), and the semimechanistic growth model (15.7; 95% CI, 14.1 to 19.9) were all similar to those reported in a pooled analysis of growth data (16.4; 95% CI, 15.1 to 17.8) (17) and were also similar to the approximately 10-fold multiplication rate reported in natural infections (19,20). During parasite growth model development, the log-linear and cosine-wave models estimated F SUR as less than 1% with high relative standard errors of 57% and 41%, respectively. The semimechanistic growth model estimated F SUR as 5.4% with a relative standard error of 19%. Pooling data from several studies conducted at QIMR Berghofer Medical Research Institute (QIMR), a separate pharmacokinetic/pharmacodynamic model estimated the fraction of survival for inoculated parasites as 5%. In order to avoid bias when comparing parasite growth models and the difficulties in estimating this parameter from the available data, F SUR was fixed to 5% based on our semimechanistic growth model and the unpublished model from QIMR.
In the semimechanistic growth model, the mean starting time point of parasite sequestration was estimated at 29.1 h. Sequestration is not likely to start at the same time in all individuals, but the estimated interindividual variability in this parameter was large with low precision. This result is most likely due to the limited number of measurements of parasitemia across the parasite life cycle. Thus, interindividual variability was retained only in the length of the parasite life cycle and parasite multiplication rate in the final semimechanistic growth model.
The parasiticidal effect of piperaquine on P. falciparum predominantly affects latestage parasites (trophozoites and schizonts) (18,21,22). Therefore, the parasiticidal effect of piperaquine was assumed to be against only the late-stage parasites in the semimechanistic growth model, both in peripheral circulation and in the sequestered vascular compartment. The goodness of fits and individual fits demonstrated that the developed model described the observed data adequately both in recrudescent and nonrecrudescent individuals. Simulation-based diagnostics of the final pharmacokinetic-pharmacodynamic model demonstrated a good predictive performance of cured individuals, but the simulated prediction intervals of individuals with recrudesce showed relatively large variability. This finding could be explained by the small proportion of recrudescent data available in this study.
The estimated in vivo EC 50 of 5.43 ng/ml (95% CI, 1.68 to 7.33 ng/ml) was similar to the mean IC 50 values of piperaquine against the 3D7 parasite strain, which was reported in two in vitro susceptibility studies (1.69 to 7.34 ng/ml) (23,24). However, this estimated in vivo EC 50 was lower than the in vitro IC 50 value reported for sensitive infections in the field. The in vitro IC 50 of piperaquine reported in 2011 in three sites in Cambodia were 10.7 ng/ml (interquartile range [IQR], 7.34 to 15.5 ng/ml) in Ratanakiri, 10.3 ng/ml (IQR, 8.09 to 14.0 ng/ml) in Preah Vihear, and 10.5 ng/ml (IQR, 6.37 to 18.2 ng/ml) in Pursat (5). The lower estimated EC 50 in this study than IC 50 values reported in the field might be partially explained by the different genetic background of the parasites, i.e., the 3D7 strain used in the current study was originally isolated from Africa and is sensitive to chloroquine and several antimalarial drugs, including piperaquine (23-25). The median PRR LC derived from the final pharmacodynamic model was 2.68 Â 10 2 to 3.02 Â 10 2 at standard therapeutic doses. These values are similar to the 48-h parasite reduction ratio (PRR 48 ) of 10 2 to 10 5 reported in the literature (26,27), which correspond to a PRR LC of approximately 4.14 Â 10 1 to 1.10 Â 10 4 when corrected for a shorter life cycle length. However, the reduction ratio reported in the present study is in the lower end of what has been reported in the literature. This difference can be explained by several factors such as the synchronicity of the infection, potential differences between the laboratory strain used in this model and clinical isolates (e.g., parasite life cycle, parasite multiplication rate, and parasite susceptibility to piperaquine), and the difference in background immunity in malaria-naive volunteers versus patients from regions of malaria endemicity with pre-exposure to malaria. Simulations of multidrug-resistant scenarios predicted a probability of treatment failure that was somewhat lower than that reported in the clinical settings in Cambodia, especially in the scenario with an extremely resistant infection (i.e., 4-fold increase in EC 50 as reported in Pursat). The rates of clinical recrudescence reported from this study in 2013 were 2%, 16%, and 46% in Ratanakiri, Preah Vihear, and Pursat, respectively. The in vitro IC 50 of piperaquine in these 3 study sites was increased by approximately 2-, 3-, and 4-fold, respectively, in 2013 compared with the values reported in 2011 (4,5). Assuming the same increase in EC 50 in the present study resulted in simulated failure rates of 8.1%, 15.2%, and 23.6%. The difference in the predicted proportion of treatment failures versus observed clinical failure rates could possibly be a result of a difference in parasite dynamics as discussed above and/or the pharmacokinetic properties of piperaquine, which might differ between healthy volunteers and patients. Differences in treatment adherence could also partly explain the difference. Moreover, the assumption of the model with respect to the dihydroartemisinin parasite killing effect (E max ), based on the parasite clearance half-life values associated with artesunate in sensitive and resistant infections from a previous study (28), might not be a perfect representation of dihydroartemisinin resistance occurring in the field. The simulations in the current study enabled a prediction of the probability of treatment failure when an additional partner drug was added to the conventional dihydroartemisinin-piperaquine regimen. This information could help with the selection of new combination therapies and optimization of dosing regimens. In all simulations of treatment failures, an additive drug effect was assumed for drugs included in the treatment. We believe this is the most conservative approach, considering the limited information available on drug synergisms/antagonisms of antimalarial drugs. However, applying a different drug interaction could yield different results. A previous study (29) showed that a model incorporating a different magnitude of interactions between antimalarial drugs, using dihydroartemisinin-piperaquine plus mefloquine as an example, predicted a reduced probability of treatment cure when the level of antagonism between piperaquine and mefloquine was high. Some limitations of the semimechanistic model developed in the current study include characterization of the fraction of sequestered parasite (F SQ ) and the relatively high uncertainty in estimation of the time to sequestration. These issues could be overcome by more frequent measurements of parasite densities during the growth phase. Furthermore, stage-specific parasite measurements could also enhance the robustness and reliability of the model. The semimechanistic growth model resulted in improved model fit compared with the log-linear model but did not result in a better model fits, in terms of DOFV and DBIC, compared with the sine-wave growth model. Nevertheless, we believe that a semimechanistic growth model, based on observed biological processes in the parasite life cycle, is preferable to an empirical description of the data and a more useful tool for translational simulations, especially since the developed semimechanistic model can be modified to include antimalarial drugs with different mechanisms of action (i.e., drug effect can be incorporated at different stages in the parasite life cycle depending on the mechanism of action). Another limitation is the potential differences between the P. falciparum strain used in the IBSM model (3D7) and those in patients with clinical malaria. If the P. falciparum strain used in the IBSM model is more drug susceptible than the strains in clinical infections, this would lead to an overestimation of the drug-mediated killing of parasites. However, the pharmacokinetic-pharmacodynamic model structure, incorporating a semimechanistic growth model, allows the flexibility to evaluate drug effects to a specific stage of the parasite life cycle. Furthermore, the model explained the processes described in natural P. falciparum malaria infections, including maturation of parasites, sequestration of mature parasites, synchronicity of infections, and multiplication of parasites. The implementation of this model structure to growth phase data from a large pool of malaria volunteer infection studies could further confirm the robustness of the model and hopefully allow for all model parameters to be estimated.
In conclusion, an in vivo semimechanistic model of parasite growth and clearance was developed in participants inoculated with P. falciparum malaria, and model parameters (E max , EC 50 , and PRR LC ) associated with piperaquine pharmacokinetic-pharmacodynamic effects were estimated. This semimechanistic parasite model provides important insights and could be an important tool in the development of novel triple-combination therapies and for dose optimization of piperaquine and other antimalarial drugs.

MATERIALS AND METHODS
Study design. Pharmacokinetic and pharmacodynamic data were collected from 24 healthy volunteers who participated in a previously described phase-Ib single-center clinical trial (26). In brief, the study was conducted at the contract research organization Q-Pharm (Brisbane, Australia). Malaria-naive subjects aged 18 to 50 years old who met all of the inclusion and none of the exclusion criteria were eligible to enroll in the study. The study protocol was approved by the QIMR Berghofer Human Research Ethics Committee and was registered in the Australian and New Zealand Clinical Trials Registry (ANZCTRN12613000565741). Eligible participants were intravenously inoculated with approximately 1,800 viable P. falciparum-infected human erythrocytes (chloroquine-susceptible 3D7 strain) on day 0. A single dose of piperaquine was given to participants in four different cohorts. The doses of piperaquine phosphate were 960 mg (cohort 1), 640 mg (cohort 2), and 480 mg (cohorts 3a and 3b). Piperaquine was administered on day 7 (cohort 1 and 3b; n = 11) or day 8 (cohort 2 and 3a; n = 13) after inoculation. The parasite density in inoculated participants was quantified using a quantitative PCR (qPCR) that targets the P. falciparum 18S rRNA gene (30). Treatment for malaria recrudescence consisted of artemether-lumefantrine in cohorts 1 and 2 and a second dose of piperaquine (960 mg) in cohorts 3a and 3b. At the end of the study, artemether-lumefantrine was given to all participants as a curative malaria treatment. Plasma concentrations of piperaquine were measured using liquid chromatography and mass spectrometry. The plasma samples were collected before piperaquine administration; at 0.5, 1, 2, 3,4,6,8,12,24,48,72,96, and 144 hours after piperaquine administration; at 8, 11, and 14 days after treatment; and at the end of study (day 28 for cohorts 1 and 2, day 37 for cohort 3a, and day 35 for cohort 3b). For cohorts 3a and 3b, an additional sample was taken 18 days after piperaquine administration. The range of the assay was 0.5 to 1000 mg/liter for piperaquine. The coefficient of variation across four different concentration levels was 1.2% to 4.4% (n = 10) for the intraassay and 2.5% to 6.6% (n = 10) for the interassay comparisons. The accuracy of the assay was 97% to 104% (n = 10).
Pharmacometric analysis. The population pharmacokinetic and pharmacodynamic analysis was performed using nonlinear mixed-effects modeling in NONMEM, version 7.4 (Icon Development Solution, Ellicott City, MD). RStudio version 1.2.1335 (31), Xpose version 4.0, Pirana version 2.9.4 (32), and Pearl-speaks-NONMEM (PsN) version 4.7.0 (33) were used for model diagnostics and visualization of results. Piperaquine base plasma concentrations were transformed into their natural logarithms prior to pharmacokinetic model development. None of the samples had piperaquine concentrations measured below the lower limit of quantification. The first-order conditional estimation method with interaction (FOCE-I) was used throughout the pharmacokinetic model building process. The average parasite density from qPCR measurements (parasites/ml) was transformed into the total circulating parasites (P CIR ; parasites) prior to pharmacodynamic model development (equation 2). The calculated total circulating parasite densities were transformed into their natural logarithms prior to pharmacodynamic model development. qPCR measurements below LOD (1 parasite/ml) were treated as categorical data and were modeled simultaneously with the reported continuous parasite density data above LOD, using the M3 method (34). Pharmacodynamic parameters were estimated using the FOCE-I and the Laplacian method (35). The difference in objective function value (DOFV) was used as a statistical criterion for discrimination between nested models. The difference in Bayesian information criterion (DBIC) was used when comparing nonnested models (36).
The descriptive performance of the model was assessed by goodness-of-fit diagnostics, and the predictive performance of the model was evaluated by simulation-based diagnostics. Eta and epsilon shrinkages were used to evaluate the reliability of the individual estimates and the ability to detect model misspecification in the goodness-of-fit diagnostics (37). The predictive performance of the final model was illustrated by visual and numerical predictive checks (n = 2,000). The 5th, 50th, and 95th percentile of the observed concentrations was overlaid with the 95% CI of each simulated percentile to detect model bias. Model robustness and nonparametric CI were evaluated using a bootstrap methodology (n = 1,000).
Population pharmacokinetic model of piperaquine. Pharmacokinetic parameters were assumed to be log-normally distributed, and interindividual variability was therefore implemented with an exponential function. Interoccasion variability, also implemented with an exponential function, was investigated to reflect the random variability between dosing occasions. An additive error model and a combined additive and proportional error model, both on the logarithmic scale, were evaluated. To evaluate the effect of body size on the pharmacokinetic properties of piperaquine, body weight was implemented as an allometric function on all clearance and volume of distribution parameters (equation 1).
where u i denotes individual clearance or individual volume of distribution parameter, u denotes the typical value (population mean) of clearance or volume parameters, BW i denotes individual body weight, BW median denotes median body weight of the participants, and n was set to be equal to 0.75 for clearance parameters and 1 for volume parameters. Additional covariate relationships including age, sex, and race were examined using a stepwise forward inclusion (P , 0.05, DOFV = -3.84), followed by stepwise backward elimination (P . 0.001, DOFV= -10.83) procedure. P. falciparum parasite growth model. The initial total number of circulating parasites (P CIR ) was fixed to 1,800 parasites (equivalent to the approximate number of viable inoculated parasites) for all investigated parasite growth models, and the fraction of parasite survival (F SUR ) after inoculation was fixed to 5%, based on results from a pharmacokinetic/pharmacodynamic model using pooled induced blood-stage malaria data at QIMR (unpublished work). The parasite life cycle was fixed to 38.8 h (17). The average parasite density from qPCR measurements (parasites/ml) at each time point was transformed to the P CIR based on individual body weight multiplied by the average blood volume which was assume to be 80 ml/kg (equation 2) (38). P CIR ¼ average parasite density ðparasites=mlÞ Â body weight ðkgÞ Â 80 ðml=kgÞ The calculated number of total circulating parasites was transformed into natural logarithms for parasite dynamic model development. Initially, only the growth-phase data were used to develop the parasite growth model (Fig. S6). Three different types of models were evaluated, namely, log-linear growth model, sine-wave growth model, and a semimechanistic growth model. The parasite dynamic model that best described the observed data was carried forward to evaluate parasite dynamics after piperaquine administration. The details of each parasite growth model are described below.
Log-linear growth model and sine-wave growth model. The log-linear growth model used to describe parasite growth data was implemented using a differential equation to explain the change of parasite density with time (equation 3).
where P CIR denotes the number of total circulating parasites (parasites) and k G denotes parasite growth rate (h 21 ). The sine-wave model used in the previously published pooled analysis of parasite growth data from volunteer infection studies (17) was implemented to describe parasite growth data in the current study (equation 4).
where, ln(P CIR ) denotes the natural logarithm of total circulating parasites (parasites), a denotes y-intercept (i.e., the P CIR at time zero), time denotes time after parasite inoculation (h), C denotes sine-wave amplitude, T PC denotes duration of the parasite life cycle (fixed to 38.8 h), and k denotes sine-wave phase shift. The overall parasite multiplication rate per life cycle (PMR LC ) from the log-linear growth and sinewave growth model was calculated using the estimated parasite growth rate and the duration of the parasite life cycle (equation 5).
Semimechanistic growth model. The semimechanistic model was developed based on prior knowledge of the P. falciparum life cycle. Time windows corresponding to P. falciparum parasite stages were based on the microscopic observations previously reported (39). These time windows were corrected for a shorter life cycle length used in the current study (38.8 h). The proposed model consisted of three parasite compartments (Fig. 4). The first parasite compartment (P1) represents the small rings that are circulating in the peripheral blood. The second parasite compartment (P2) represents the large rings, trophozoites, and schizonts that are also circulating in the blood. These two parasite compartments represent the total circulating parasites in the peripheral blood, which can be measured by qPCR and microscopy. However, the total parasite biomass also includes sequestered parasites, which attach to endothelial cells. The third parasite compartment (P3) represents the sequestered parasites. Thus, combining the parasite number in all parasite compartments (P1, P2, and P3) yields a total parasite biomass in the body.
Three square-wave functions were used to regulate the movement of the parasites between each compartment at specific time periods (REG 1 , REG 2 , and REG 3 ). These square-wave functions were set to change from 0 to 1 (i.e., "on-and-off") during certain time periods to move parasites between the different parasite compartments. These square-wave functions were described by the following equations (equation 6 to 8).
REG ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where S 1 denotes the first sine-wave function, S 2 denotes the second sine-wave function, SH 1 denotes time when the function remains in state 0 (h), SH 2 denotes the peak shift time from state 0 to 1, T PC denotes duration of the parasite life cycle (fixed to 38.8 h), and REG denotes the square-wave function regulating parasite movement. SH 1 and SH 2 values used to generate each square-wave function (REG 1 to REG 3 ) for each parasite compartment were based on time windows corresponding to P. falciparum parasite stages (39). Details on how these functions were implemented in the model are shown in the supplemental material (NONMEM code).
In the semimechanistic growth model, the parasite number in P1 was initialized with 1,800 parasites according to the approximate number of viable inoculated parasites. The parasites in P1 mature over time and started moving to P2 at 9.7 to 19.4 h (REG 1 ), and at 19.4 h, the entire parasite population in P1 moved to P2. The parasites in P2 continued to mature to schizonts, and either stayed in P2 or cytoadhered to the vascular endothelium and moved to P3 (sequestration). The sequestration of parasites was regulated (REG 2 ) to occur during the latter half of the parasite life cycle at 19.4 to 38.8 h since sequestration begins at the large ring stage (40,41). At the end of the parasite life cycle (38.8 h), the matured schizonts in P2 and P3 ruptured and released new merozoites back to compartment P1 (REG 3 ). The PMR LC associated with the rupture of schizonts was estimated. This semimechanistic growth model was described by the following differential equations system (equation 9 to 11).
where k MAT denotes the first-order rate constant for small rings to become mature parasites; k SQ denotes the first-order rate constant for parasite sequestration; k RUP denotes the first-order rate constant of schizont rupture; and REG 1 , REG 2 , and REG 3 denote the square-wave functions regulating the time and duration of parasite movement in each parasite compartment. Details of the square-wave functions used to regulate parasite movement are presented in the Fig. 4. In order to describe the parasite sequestration in a quantitative manner, the parasite sequestration rate was parameterized as a fraction of sequestered parasites (equation 12).
where F SQ denotes the fraction of sequestered parasite (%) and T SQ denotes the onset of parasite sequestration (restricted to be between 19.4 and 38.8 h).
In vivo parasiticidal effect of piperaquine. The model that best described the parasite growth dynamics was carried forward and linked to the pharmacokinetic model of piperaquine. The final parameter estimates from the best performing model was used to impute individual growth time profiles. All total circulating parasite data, including the total circulating parasites after piperaquine administration, were used to estimate the parameters associated with the drug-dependent parasite elimination (Fig. S6). The parasiticidal effect of piperaquine was implemented as an E max function (equation 13).
where EFF denotes the parasiticidal effect of piperaquine (h 21 ), E max denotes the maximum parasite killing rate of piperaquine (h 21 ), C P denotes piperaquine plasma concentration (ng/ml), EC 50 denotes the plasma concentration of piperaquine (ng/ml) associated with half of maximum parasite killing rate, and g denotes the hill factor. Piperaquine was assumed to have an effect on the later stages of blood-stage parasites (P2 and P3 compartment) based on previously published information (21,22). The implementation of the drug effect was described by the following differential equations (equations 14 and 15). Cure was assumed to be achieved when the total parasite biomass (P1 1 P2 1 P3) was less than 1 parasite, triggering the PMR LC to be 1 (i.e., resulting in no parasite growth).
The simulation-based diagnostics (i.e., visual predictive checks) of the final pharmacokinetic-pharmacodynamic model was based on 1,000 simulations using the individual parameter estimates from the final pharmacokinetic model and the final parameter estimates from the pharmacokinetic-pharmacodynamic model. Frequent dummy time points were added to the data set for simulations (every 5 h, from 0 h to 576 h). The visual predictive checks were stratified on the day of piperaquine treatment and the predicted treatment outcome (curative versus recrudescent infection). The observed data were overlaid with the 90% prediction interval from 1,000 simulations to evaluate the predictive performance of the model.
Clinical scenario simulations. The final pharmacokinetic-pharmacodynamic model describing the dynamic parasite growth and piperaquine in vivo parasiticidal effects was used to perform population-based simulations of clinical scenarios in NONMEM. Simulations were performed to predict the probability of treatment failure at different levels of drug resistance. The total number of circulating parasites was initialized with 10 6 or 10 10 parasites in order to simulate asymptomatic and symptomatic infections, respectively. The following six different clinical scenarios with resistant infections were simulated to predict the probability of treatment failure: (1) absence of artemisinin resistance and absence of piperaquine resistance, (2) presence of artemisinin resistance and absence of piperaquine resistance, (3) absence of artemisinin resistance and presence of piperaquine resistance (2 Â EC 50 of piperaquine), and (4) to (6) presence of artemisinin resistance and presence of piperaquine resistance (2 Â EC 50 , 3 Â EC 50 , and 4 Â EC 50 of piperaquine).
Dihydroartemisinin was used as a representative of the artemisinin derivative, and pharmacokinetic parameters of dihydroartemisinin were taken from a previously published population pharmacokinetic analysis (42). The parasiticidal effect of dihydroartemisinin was adjusted based on the observed parasite clearance half-life from the Tracking Resistance to Artemisinin Collaboration (TRAC) study data (28). The reported mean parasite clearance half-life for sensitive (2.5 h) and resistant (6.2 h) infections from the TRAC study was used to calculate the parasite clearance slope and generate parasite clearance profiles, starting at an initial total circulating parasite density of 10 11 parasites (Fig. S5). Simulations were performed in Berkeley Madonna (43) using the developed semimechanistic growth model to adjust the E max values of dihydroartemisinin to match the parasite clearance profiles generated from the parasite clearance half-life reported in the TRAC study. The derived adjusted E max values, resulting in equivalent residual total circulating parasites at 72 h as seen in the TRAC study, were used throughout the simulations.
The parasite killing effect of dihydroartemisinin was implemented in all parasite compartments (P1, P2, and P3) using an E max function (equation 13) because dihydroartemisinin has an effect on almost all stages of the parasite life cycle (44). Treatment failure was defined as a predicted total parasite biomass of .1 parasites at 30 days after the first dose of piperaquine was given. Additionally, simulations were conducted to predict the probability of treatment failure when adding an additional hypothetical drug to the conventional dihydroartemisinin-piperaquine regimen. The effect of the hypothetical drug was implemented in the same manner as the parasite killing effect of piperaquine (added on P2 and P3 compartments). Different efficacy and duration of action of the hypothetical drug were investigated, including drugs demonstrating a PRR LC of 10 1 , 10 2 , and 10 3 and a duration of action of 1 week, 2 weeks, 3 weeks, 4 weeks, and 5 weeks. The following differential equations described the drug effects for the simulations (equation 16 to 18) (18) where EFF 1 denotes the parasiticidal effect of dihydroartemisinin, EFF 2 denotes the parasiticidal effect of piperaquine, and EFF 3 denotes the parasiticidal effect of the third hypothetical drug. Each simulation scenario consisted of 4,800 simulated patients, including 48 individuals with various times of first dose with a 1-h difference among each simulated patient from 0 to 48 h (100 simulations). The variation mimics a real-life scenario where patients present to the clinic and receive treatment at different stages of the parasite life cycle. The standard 3-day dose of dihydroartemisinin-piperaquine (120/ 960 mg) for a patient weighing 60 kg was used in all simulations. A schematic illustration of the structural model used for simulations is presented in the supplemental material (Fig. S7).

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, PDF file, 2.5 MB.