**DOI:**10.1128/AAC.01712-10

## ABSTRACT

Pharmacological mechanism-based modeling was refined and used to develop an *in silico* model of antimalarial drug treatment validated against clinical and field data. We used this approach to investigate key features of antimalarial drug action and effectiveness, with emphasis on the current generation of artemisinin combination therapies. We made the following conclusions. (i) The development of artemisinin tolerance and resistance will, unless checked, have an immediate, large impact on the protection afforded to its partner drug and on the likely clinical efficacy of artemisinin combination therapies. (ii) Long follow-up periods are required in clinical trials to detect all drug failures; the follow-up periods of 28 days recommended by the World Health Organization are likely to miss at least 50% of drug failures, and we confirmed recent suggestions that 63 days would be a more appropriate follow-up period. (iii) Day 7 serum drug concentrations are a significant risk factor of failure, although, paradoxically, receiver operating characteristic curve analysis revealed that their predictive power is relatively poor. (iv) The pharmacokinetic properties of the partner drugs in artemisinin-containing combination therapies are the most important determinants of treatment outcome, particularly the maximum killing rate. We discuss the assumptions made in such modeling approaches and how similar approaches may be refined in future work.

## INTRODUCTION

Malaria remains an important public health problem in countries where malaria is endemic, and prompt treatment with effective antimalarials is crucial to patient survival. Correctly dosing patients with malaria is a constant problem, and the standard dosage regimens are usually based on information from adults with uncomplicated malaria. This neglects the differences seen in patients' age and weight, the use of concomitant medication, genetic variation in pharmacokinetics (PK), especially in the groups most at risk (pregnant women, children, and patients with severe malaria), and natural variation in parasite drug susceptibility. Furthermore, dosing regimens are determined during efficacy trials in which exact doses are calculated according to weight and taken under supervision, with strict adherence to dose timings. In real life, effectiveness is more important than efficacy. Drug effectiveness in the field is invariably less than the efficacy seen during trials for several reasons. (i) Drug regimens consisting of more than one dose are vulnerable to the effects of patient compliance. Some treatments, such as artemether-lumefantrine, require up to 6 doses per course, so compliance may be much worse than single-dosage therapies such as sulfadoxine-pyrimethamine (SP). (ii) Drug trials tend to be small and may not capture the full natural variation in PK and pharmacodynamics (PD) within a population. Interindividual variation in PK/PD may result in patients with subtherapeutic drug concentrations. (iii) Health systems in countries where malaria is endemic are not equipped to treat patients according to their weight, which poses a particular problem in dosing children. A child's dose is determined according to an age band, and while this is more practical, it inevitably leads to a large proportion of dosages either above or below those recommended (61).

Clinical trials are expensive and relatively small and cannot ethically measure effects of factors such as poor compliance and underdosing. Subsequent field studies on effectiveness may infer the impact of factors such as compliance (52), but these can be done only after the drugs are deployed. Neither type of study can determine how robust the regimens are in the face of small changes in parasite drug sensitivity or how vulnerable the regimens are to the evolution of drug resistance. In this paper, we investigate whether we can develop pharmacological models of antimalarial drug treatment that are sufficiently compelling that these questions can be addressed usefully *in silico*. The consensus method of modeling the effect of drug treatment on an infection is to track the change in the total number of parasites in the body over time following treatment. A review of this “mechanism-based” PK/PD modeling of antimicrobial drugs by Czock and Keller (12) lists 28 papers that have used (with modifications) this model, including its application to malaria, by Austin et al. (3), to identify the optimal pattern of drug administration to clear an infection; by Hoshen et al. (30, 31), to determine whether treatment outcomes are improved with a split chloroquine dose; and by Simpson et al. (55), to compare the development of resistance after the *de novo* use of 2 of the most widely used antimalarial doses of mefloquine (MQ). We investigated four drugs, chloroquine (CQ), lumefantrine (LF), MQ, and piperaquine (PQ), given as monotherapies and as components of artemisinin-containing combination therapies (ACTs) with artesunate (AS), artemether (AR), or dihydroartemisinin (DHA) as appropriate. Chloroquine is no longer officially deployed through the formal health sector for treatment of Plasmodium falciparum (with the interesting exception of Guinea-Bissau [27, 35, 63]). It is included here for historical comparison and because it results in higher failure rates (see below) than the ACTs, allowing us to compare the different dynamics of drug failure. Following convention, CQ, LF, MQ, and PQ are collectively known as the partner drugs within these ACTs. We demonstrate that this modeling approach generates results consistent with field data and try to identify the most important factors that determine drug treatment success and failure.

## MATERIALS AND METHODS

Basic model.The rate of change in the total number of parasites present in a single patient over time depends on the parasite multiplication rate discounted by the proportions killed by antimalarial drugs and cleared by host defenses such as immunity. The methodology is based on the following standard differential equation:
*P* is the number of parasites in the infection, *t* is the time after treatment (days in these simulations), *a* is the parasite growth rate (per day), *f*(*C*) represents the drug killing rate for parasites, which depends on the drug concentration *C*, and *f*(*I*) represents the host's background immunity to the infection. The drug killing function *f*(*C*) is given by the following equation:
*C* is the drug concentration (mg/liter), *V* is the maximal parasite-killing rate constant (per day), *K* is the concentration at which 50% of the maximal killing rate occurs (mg/liter), and *n* is the slope of the dose-response curve. The drug concentration decays over time, as follows:
*k* is the terminal elimination rate constant (obtained from the drug half-life) and *C*_{0} is the drug concentration at time zero, i.e., immediately after treatment (drug absorption and conversion are assumed to be instantaneous in this methodology; see later discussion). The immunity function *f*(*I*) is assumed to be time independent (i.e., immunity is not acquired over the course of treatment) and is a constant value depending on age and transmission settings (see part 1 of the supplemental material).

Combining equations 1, 2, and 3 and integrating allows us to predict the number of parasites at any given time point *t* after treatment, as follows:
*P*_{0} is the initial number of parasites present in the body at time zero, when treatment time commences. It is assumed that if and when *P*_{t} falls below 1, an infection has been cleared. Note that this methodology tracks the total number of parasites in the body, whereas clinical observations are given as percentages of infected red blood cells (RBC); simple arithmetic using the number of RBC appropriate for patient age and weight allows conversion between the scales.

Extension of basic methodology for malaria.Antimalarial regimens typically contain several doses spread over 3 days (Table 1). The concentration at time zero (*C*_{0}) is dependent on the existing drug concentration in the blood, augmented by any new drug dosage administered, as follows:
*C*′ (*C*′ = 0 if it is the first dosage), *D* is the drug dosage (mg) given, *V*_{d} is the volume of distribution (liters/kg), and *W* is the weight of the patient (kg). The new dosage is converted to the concentration in the blood, assuming instantaneous absorption and distribution; the dosage can be reduced to reflect a drug bioavailability of less than 100%.

Equations 4 and 5 are applied to each separate time interval in the drug regimen. For example, in a 3-day regimen of CQ (Table 1), there would be 3 intervals, i.e., day 1, day 2, and day 3 onwards; for artemether-lumefantrine, there would be 5 half-day periods, followed by dosage 6 onwards.

Antimalarials are now invariably deployed as combination therapies to improve therapeutic efficacy and to delay the development of drug resistance, although we do recognize that many other drugs, including monotherapies, are still available through informal health services. We also note that the use of quinine, particularly via continuous intravenous infusion in the case of severe malaria, is still administered frequently. However, this model focuses on the treatment of uncomplicated malaria with orally administered drugs. It is assumed that two drugs act independently of each other in the combination (see Discussion), and we simply expand equation 1 as follows:
*K*, *C*_{0}, *n*, *k*, and *V* indicate whether the parameter value refers to drug 1 or 2. This equation can be used to find the change in parasite number between any two given time points; in this model, it was used alongside equation 3 to update the parasite load and drug concentration daily. For clarity, *K* is referred to here as the 50% inhibitory concentration (IC_{50}).

Acquired immunity was incorporated using the data of Pongtavornpinyo et al. (46), who identified three possible measures of acquired immunity for use in modeling, based on parasite biomass, the incidence of severe malaria, and the proportion of symptomatic infections. For the purposes of this model, the equations were standardized to return a minimum value of 0 (when age = 0 and the entomological inoculation rate [ EIR] = 0) and a maximum value of 1 and then scaled by a factor φ to produce the *f*(*I*) used in equations 1 and 6 (see part 2 of the supplemental material for details). We use immunity only briefly here (to simulate clinical trials) but include it as a proof of principle.

The dosing regimens of the drugs investigated and estimates of the PK/PD parameters employed to model treatment are given in Table 1 and shown in part 1 of the supplemental material. This model has been implemented in R (version 2.9.2) (48), although earlier versions were run in Excel and Maple 12. All packages generated the same results, and the results presented here were generated in R. The model was run in half-day time steps, using equation 7, for the first 7 days to allow for multiple-dose regimens and in 1-day time steps thereafter to speed up simulations. It is possible to use equation 4 to find the treatment outcome algebraically after the final dose (30), but this approach of single daytime steps is more explicit, allows easy calculation of factors such as parasite clearance time and period of chemoprophylaxis, and allows easy incorporation of factors such as stochastic variation in predicted parasite numbers. The discrete time step and algebraic analysis do, of course, give the same result.

Model validation and analysis.The first step in model validation is to check that the outputs match observations made in the field. The chief criteria are that the regimens give a reliable cure rate (except in the case of short-term artemisinin monotherapies), that parasite clearance times (PCT) are plausible, and that times until new infections are noted (period of chemoprophylaxis [PoC]) are reasonable. The first term is self-explanatory: PCT is the time taken for the infection to fall below the limit of microscopic detection, which we assumed to be 10^{8}. The PoC describes the length of time that a drug treatment suppresses the appearance of new infections; this is often reported in clinical trials. We assumed that drugs do not affect parasites during the liver phase, as is believed to be the case for the drugs considered here (13). The earliest parasites able to cause a patent reinfection emerge from the liver at precisely the time when the drug concentration from previous treatments has declined sufficiently that parasite reduction turns into an increase, i.e., when *dP/dt* becomes positive. Successful reinfection can occur before this point (i.e., when *dP/dt* is <0), but the number of parasites will initially fall, and hence reinfection will occur fastest if parasites emerge on the day that *dP/dt* first becomes positive. The time when *dP/dt* becomes greater than zero is found by evaluating equation 4 to obtain *P*_{t} at the end of each daily time-step and finding when *P*_{t} is >*P*_{t−1}. Only the partner drug was considered in calculating PoC because we assumed that the artemisinin in an ACT had disappeared by the time the first drug became permissive for parasite growth. We assumed that 10,000 merozoites emerge from the liver and become detectable when their numbers exceed 10^{8}; the time taken for this to occur can be estimated by solving equation 4 with respect to *t*, where *P*_{0} = 10^{5}, *P*_{t} = 10^{8}, and *C*_{0} is the drug concentration on the day that *dP/dt* becomes >0. The PoC is then found by adding together the time taken for *dP/dt* to become positive and the time required for a new infection to become patent. To obtain a rough estimate of when new infections may reasonably be observed, we increased the PD parameters and parasite growth rate by 1 standard deviation (assuming a coefficient of variation [CV] of 0.3 [see below]) to represent parasites that are slightly above average in their ability to resist the drug. The PoC will vary from person to person, depending on variations in human PK, parasite drug sensitivity, and growth rates. We later allowed for variation in all PK and PD parameters, resulting in the PoC being different for each person; in this case, the 5th centile value of PoC is reported as being indicative of when new infections may plausibly first start to be observed in a high-transmission setting.

The initial analysis simply involved varying each individual parameter value in turn (while keeping other values constant) to find the size of change that resulted in treatment failure. This gave an initial indication of the relative importance of each parameter in determining treatment success or failure. In reality, all parameter values vary simultaneously. This was incorporated by adding variation to seven different model parameters: the five PK/PD parameters listed in Table 1, the maximal parasite growth rate *a*, and the number of parasites present at time of treatment *P*_{0}. With the exception of the initial parasite number, all variables were assumed to be distributed normally, with a CV of 30%; the mean values are given in Table 1. Values must be positive, so the program checked if this was true; if values were less than zero, the model generated another random number in the same way until a positive number was chosen. The number of parasites present at the start of treatment was chosen randomly from a uniform distribution of between 10^{10} and 10^{12}. Populations of 10,000 patients were simulated for the following drug regimens: CQ, LF, MQ, and PQ monotherapies and CQ-AS, MQ-AS, LF-AS, LF-AR, and PQ-DHA ACTs. The results were analyzed by logistic regression (LR) with treatment failure as the outcome, using the PASW Statistics 18 package. The independent variables included the initial parasitemia (*P*_{0}), parasite growth rate (*a*), drug elimination rate constant (*k*), slope factor (*n*), IC_{50}, volume of distribution (*Vd*), and maximal parasite-killing rate (*V*). LR analyzed the effect of a one-unit change in the value of each parameter, but the scales of the parameter values were so diverse that it made comparisons impossible (for example, a 1-unit change in IC_{50} is proportionally much larger than a 1-unit increase in volume of distribution) (Table 1). The parameter values were therefore converted to *z* scores because the odds ratio (OR) statistic on this scale indicates the increased or decreased risk of failure associated with a 1-standard-deviation increase in the parameter value, allowing easy comparison between variables. The changes in log likelihood (LL) and Wald statistics measured the reduction in model fit when individual parameters were omitted: the larger the value, the worse the fit was, and hence the more important the parameter. Thus, Wald and LL statistics measured the relative importance of each parameter's variation in drug failure.

There is considerable interest in how the evolution of resistance to the artemisinin component may compromise the long-term effectiveness of ACTs. We investigated two aspects of this process. First, we measured how an increasing IC_{50} of the artemisinin component would reduce the protection afforded by its partner drug. This was achieved by increasing artemisinin IC_{50} values above their default values and measuring how much of a change in the partner drug IC_{50} would be required for drug failure to occur. Second, we measured how increasing the IC_{50}s of artemisinins would increase ACT failure rates if resistance to the partner drug were already present. As before, we allowed a 30% CV in the default parameter values and measured failure rates while increasing the IC_{50} of the artemisinin component (and still allowing a CV of 0.3 around this new artemisinin IC_{50} value). Both of these measures (change in partner drug IC_{50} and drug failure rates) were standardized to allow direct comparison between different ACTs by constructing a standardized “protection index” so that a value of 1 indicates the basal value at the original artemisinin IC_{50} value and a value of 0 means that the artemisinin component is useless and that the measure has become “maximal,” i.e., identical to that of the partner drug monotherapy; for example, if *f* is the failure rate, then its standardized protection index is 1 − [(*f* − basal)/(maximal − basal)].

“Clinical trial” simulations were run to compare model outputs to field data. The data recorded were those typically measured in trials and available to investigators, i.e., outcome (success/failure), number of parasites present at time of treatment, age of patient, place of residence, and day 7 serum level. Place of residence and age of subjects were included to demonstrate the effects of immunity (see part 2 of the supplemental material) with a scaling factor of the immune function, Φ, set to 0.8: village A had an EIR of 10 and village B an EIR of 100. The clinical trial simulations included 400 individuals, 200 from village A and 200 from village B (aged 6 months to 15 years), and all were treated with artesunate-mefloquine, using the default parameters in Table 1, except that the mean IC_{50} of mefloquine was increased to 0.8 to allow a proportion of failures more typical of clinical trials where a drug is failing. These parameters were varied with a CV of 0.3 to allow treatment failures (see above; there was no point in analyzing a clinical trial where everyone was cured). The results were analyzed using logistic regression with treatment outcome as the dependent variable; independent variables were age (≤5 years or >5 years), location, low day 7 drug serum levels (defined as values below the 15th centile), and initial parasitemia.

Low drug serum levels on day 7 have been shown to be a significant risk factor for treatment failure (see later discussion), so we (arbitrarily) defined low levels as being below the 15th centile and determined the risk associated with low day 7 serum levels by the odds ratio and population attributable risk percentage (PAR%), a measure of the percentage of treatment failures that could be avoided if adequate drug levels were achieved throughout the population. The diagnostic accuracy of the day 7 serum concentration was further explored using PASW Statistics 18.0 to generate receiver operating characteristic (ROC) curves. The correlation between drug serum concentrations on days 3, 5, 7, and 10 and the areas under the drug concentration curve (AUC) for days 0 to 25, 0 to 50, 0 to 100, and 0 to infinity was investigated. The AUC between any two time points was found using equation 8, while AUC_{0-∞} was found using equation 9:

## RESULTS

Use of the default parameter values without variation resulted in CQ, LF, MQ, and PQ monotherapies reliably clearing malaria infections when the drug dosage schedules simulated were those recommended by the WHO (73) (Table 2), assuming that initial parasite numbers were below 10^{12}. The use of artemisinins (AS, AR, and DHA) as monotherapies over 3 days failed to clear infections (Table 2), although they did cure as monotherapies if they were given over 7 days (data not shown). When the initial parasite number was 10^{10}, the parasite clearance times for LF, MQ, and PQ began at 3 days and increased by approximately 1 day with each 10-fold increase in initial parasite number. The parasite clearance time for CQ began at 4 days when initial parasitemia was set to 10^{10} and increased to 6 and 7 days with increasing initial parasite numbers. For all partner drugs, the addition of either AS or DHA reduced PCT by 1 to 2 days, while the addition of AR to LF reduced the PCT by up to 3 days. Estimates of PoC with the nonvaried default parameters (Table 2) were always larger than expected (see later discussion).

The percent change that can be tolerated in a given drug parameter before treatment fails is shown in Table 3, assuming an initial parasitemia of 10^{10}, to crudely illustrate how sensitive a drug regimen is to the natural variation between individuals (recall that the partner drugs are able to clear parasites as monotherapy, so it was pointless to vary individual parameter values of the artemisinin component). Studying the parameters without variation suggested that chloroquine and piperaquine regimens were most sensitive to change, while mefloquine and lumefantrine regimens were the most robust. All regimens appeared to be particularly sensitive to changes in the parasite-killing rate constant (*V*); for example, a reduction of only 40% resulted in treatment failure for CQ monotherapy. The slope of the concentration-effect curve had little or no effect on drug treatment outcome.

More realistically, allowing a 30% CV in parameter values reduced overall true cure rates to 85% (CQ), 90% (LF), 96% (MQ), and 90% (PQ) for monotherapies (Table 2). The addition of artemisinins increased these to 88% (AS-CQ), 93% (AS-LF), 95% (AR-LF), 97% (AS-MQ), and 91% (DHA-PQ). The apparent cure rates at 28 and 63 days were recorded to investigate the effects of differing periods of follow-up in clinical trails (Table 2). The cure rates were overestimated at both time points, but the apparent cure on day 63 gave a closer approximation of the true cure rate than the day 28 estimates. The apparent CQ cure rate on day 63 was 6% higher than the true cure rates; all other day 63 cure rates were within 2% of the true cure rate (Table 2). The parasite clearance times found in the presence of variation were all reduced with the addition of an artemisinin to the partner drug (Table 2). The PoC estimates were normally distributed (see Fig. S3.3 in the supplemental material), with the 5th centile values presented in Table 2.

Allowing a CV of 30% in PK/PD parameters allowed the most important parameters to be identified and ranked using the Wald statistic for logistic regression (Table 4; note that overall cure rates are given in Table 2) and then compared to the results obtained by varying each parameter individually (Table 3). All model variables for monotherapies (except the MQ initial parasite number) were found to be significant factors in determining treatment outcome. For CQ, LF, MQ, and PQ, the Wald and log likelihood statistics identified the maximal parasite-killing rate (*V*) to be the most important variable and the slope factor (*n*) to be the least important variable. The volume of distribution (*V*_{d}), IC_{50}, and slope factor (*n*) of the artemisinin component in combinations (with the exception of the AS IC_{50} in the AS-LF combination) were not found to be significant determining factors in treatment success/failure. Additionally, in using the AS-MQ combination, both the initial parasite number and the artesunate maximal parasite-killing rate constant (*V*) were not significant. The most important variable in all combination therapies was the rate of parasite killing (*V*) of the partner drug, closely followed in all cases by the rate of parasite growth (*a*). Where significant, the number of parasites present when treatment began was found to be the least influential factor in combination therapies.

The impact of increasing artemisinin IC_{50} values is shown in Fig. 1, in terms of reduced protection afforded to the partner drug (Fig. 1A) and in terms of the likely increased failure rate of the ACT (Fig. 1B). Both results suggest that increasing resistance to artemisinin causes a rapid decline in protection.

The “clinical trial” simulations examined the importance of four variables typically measured in the field, namely, initial parasite number, location, patient age, and day 7 serum level. Results for artesunate-mefloquine are shown in Table 5. People from a high-transmission village were more immune and, as expected, were less likely to fail treatment; interestingly, this effect was always associated with very small confidence intervals. Young age was associated with increased risk of failure, although this was not significant in our simulation. Several different field studies have reported day 7 serum concentration as a predictor of treatment failures (8, 18, 66); our results confirmed that individuals with low day 7 serum levels have a 3-fold higher risk of treatment failure than those with normal levels.

Individuals were categorized as having either normal or low day 7 serum levels by using the 15th centile value, and its effect on treatment success/failure is reported in Table 6. These 15th centile values were found to be <0.056 (mg/liter) for CQ, <1.647 for LF, <0.637 for MQ, and <0.232 for PQ. The increased risk of treatment failure in those with low day 7 serum levels was largest with DHA-PQ (OR = 2.62) and smallest with AS-MQ (OR = 1.24) (note that we disregarded the OR values in Table 5 because these were obtained after inflating the MQ IC_{50} value [see Materials and Methods]). The PAR% results suggest that increasing drug levels sufficiently to achieve normal day 7 serum levels in all individuals could reduced failure rates by as much as 17% with DHA-PQ but by only 3% with AS-MQ. The area under the ROC curve (AUROC) found low day 7 serum levels to be a poor predictor of treatment outcome, with AUROC values ranging from 0.538 for AS-MQ to 0.646 for PQ (an AUROC value of 0.5 indicates no predictive power, and an AUROC value of 1 indicates perfect sensitivity and specificity). These poor predictive values occurred despite a strong correlation between AUROC and AUC (see Table S3.3 in the supplemental material). The estimates of OR were not noticeably affected by our choice of the 15th centile as the definition of low serum levels: the same ORs were obtained when the cutoff was defined as the 10th, 20th, or 30th centile (see part 3 of the supplemental material).

## DISCUSSION

All models need to make simplifications to make them tractable. We therefore explicitly identify and discuss the key simplifications and assumptions that underlie these analyses before moving on to discuss how well the model predictions fit with clinical data and what the model results imply for the current generations of antimalarial drugs.

The mechanism-based modeling approach (12) assumes that drugs are absorbed instantaneously (across the gut wall for most antimalarials) and converted instantaneously, if necessary, to their active metabolite. This is probably reasonable for CQ, LF, MQ, and PQ because their long half-life dominates the time taken for relatively rapid absorption and conversion, but this approach is less satisfactory for artemisinins, where the half-life is so short that time lags in absorption and conversion may arguably become important. It would be possible to remodel this component using a standard PK compartment model, as discussed below.

Within a human, drug bioavailability and the extent of absorption are important contributors to the variability of drug outcomes. Lumefantrine oral bioavailability is particularly variable and highly dependent on food intake; it is consequently poor in acute malaria cases but improves markedly with recovery (18). By assuming that all patients modeled had uncomplicated malaria and followed the dosing recommendations, we ignored complications caused by bioavailability absorption. We also ignored any toxic effects resulting from high concentrations and any impact they might have on PK/PD parameters.

We investigated two modifications of the basic methodology. We first examined the effects of assuming density-dependent growth of P. falciparum (12) (see part 1 of the supplemental material for details) and then the effects of adding stochastic variation to *P*_{t} at the end of each day by choosing the value from a Poisson distribution. Neither modification had a significant effect on the outcome, and both were subsequently removed; all results presented use a constant growth rate *a*, as detailed in part 1 of the supplemental material.

The effects of acquired immunity have largely been ignored and introduced only briefly during the simulation of clinical trials. The development and action of human immunity against malaria are vast topics of current research, and as far as we know, there is no consensus mathematical description of its acquisition. Hence, we used proxies to illustrate its basic effect and await more sophisticated descriptions to incorporate its effect (see part 2 of the supplemental material). Its omission can be justified because it is important that drugs act effectively in all humans, irrespective of immune status. Most malaria mortality is in nonimmune or poorly immune African children, so the results presented here are most appropriate to this group.

The assumptions described above concern the structure of the model. The interpretation of the results requires a secondary assumption, namely, that we have properly calibrated the PK/PD values of each drug (see part 1 of the supplemental material). Differences in drug assay sensitivity and reproducibility, combined with natural PK variation between people, resulted in a range of PK estimates. We incorporated natural variation in the model parameters to reflect interindividual variation by simply assuming that each parameter was normally distributed, with a CV of 0.3. Ninety-five percent of parameter values will lie within 2 standard deviations of the mean, i.e., from 0.4 to 1.6, assuming a standard mean of 1, a roughly 4-fold range of interindividual variability. More sophisticated calibrations could be made; for example, lumefantrine bioavailability depends critically on its ingestion with fatty food so may vary widely depending on diet, drug IC_{50} values typically vary 10-fold (40), and so on. At this stage, we were less concerned with exact calibration of the model than with its construction and evaluation, and we leave it to readers to calibrate the equations as they see fit. We do assert that PK/PD parameter values are known to vary widely, so our approach of setting the CV to 0.3 seems a reasonable first approach for investigating the general properties and robustness of the drug regimens.

Despite these assumptions made during model construction and calibration, it is gratifying that the results presented above closely match observations made in the field, particularly as described in the following. (i) CQ, LF, MQ, and PQ monotherapies gave reliable cure rates when modeling infections of drug-sensitive parasites. (ii) Calibrating the CQ model with IC_{50} values typical of resistant parasites gave high treatment failure rates. Doubling the CQ dosage (as was done in Guinea-Bissau [27, 35, 63]) restored the efficacy of CQ against these “resistant” parasites (data not shown). (iii) Artemisinin monotherapies over 3 days are unable to clear infection, but extending the treatment time to 7 days does reliably clear parasites (1). (iv) Adding an artemisinin to a failing drug did not completely prevent treatment failure, but it did reduce the failure rate by around 50% (Table 2), in line with field observations (1). The addition of an artemisinin derivative also appeared to increase the robustness of the partner drug in the face of interindividual variability (Table 3). (v) Estimates of parasite clearance times for ACTs were in the region of 2 to 4 days (Table 2), which closely matches field estimates for LF-AR (19, 33, 34, 49, 57, 64), CQ-AS (43), DHA-PQ (23, 32, 33, 42), and MQ-AS (19, 49). (vi) The predictions for PoC were similarly consistent with field data when allowing variation in PK/PD parameters (Table 2). The model predicted a PoC of 32 days for LF, which is very close to the reported range in fully sensitive parasites of 24 to 30 days (56, 58), although we do note that PoC is reduced as resistance spreads (26). The PoC for CQ was longer than expected, at 52 days, but this may be due to the long terminal elimination half-life values reported for CQ (21, 47), which may not properly reflect elimination rates at higher, physiologically active concentrations (16). Since MQ and PQ are not yet regularly deployed in areas of high transmission, we were unable to find published estimates of PoC for either drug, but we did find papers comparing reinfection rates between two different drug regimens used in the same setting. Grande et al. (23) found that the incidence of reinfections in their study site in Mali was higher for the DHA-PQ group than for the AS-MQ group and tentatively attributed this to the shorter posttreatment prophylactic effect seen with PQ than with MQ. Our model predicts a similar relationship between the prophylactic effects of PQ and MQ, with reinfections occurring up to 19 days earlier following treatment with PQ than following treatment with MQ, despite the longer half-life of PQ. This is possible because the PoC, and hence the time to reinfection, depends not just on the drug half-life but also on its dosage (increasing dosage increases the protection time) and on parasite sensitivity to the drug (56, 58). Our model prediction of a longer PoC after MQ treatments than that with LF implies that reinfections will occur sooner following LF treatment. This is in agreement with the findings of Sagara et al. (49), who found reinfections in Peru to occur more frequently following AR-LF treatment than following AS-MQ treatment. Finally, a study comparing the safety, efficacy, and tolerability of AR-LF and DHA-PQ in Zambian children (41) observed more new infections during the follow-up of patients receiving AR-LF than during that of patients receiving DHA-PQ, which is consistent with our model predictions.

In summary, the model provides a good qualitative and quantitative fit to clinical observations, so with the caveats noted above, it seems reasonable to ask what the modeling implies for current drug regimens.

All drugs achieved reliable cures under the default PK/PD parameters, with the addition of artemisinins rapidly reducing patient parasitemia (PCT) (Table 2) and hence helping in the resolution of symptoms, as expected. Monotherapy drug failures occurred at rates between 4% (MQ) and 10% (LF and PQ) when variation in PK/PD was introduced, but these are within the 10% limit suggested by the WHO to indicate a need to change the drug regimen. The failure rate with “CQ-sensitive” parasites was relatively high, at 15%. We believe this arose because good-quality parameter estimates are readily available only from fairly recent studies, where some degree of CQ resistance had already evolved. We are comfortable with this failure rate because we wanted to investigate at least one situation where resistance had already arisen naturally (rather than by our manipulation of PD parameters) and to compare this with patterns noted for highly effective drugs such as the ACTs. Adding artemisinins improved cure rates but did not eliminate all drug failures, again in line with field observations (1).

Efficacy is measured in the field by using clinical trials where continued patient follow-up is operationally challenging and problematic in most areas of endemicity. How this may affect trial results is addressed in Table 2, which reports apparent cure rates at days 28 and 63 (the proportion of patients with no detectable parasites, which may include patients with subpatent infections that will later recrudesce), as well as true cure rates. It is generally accepted that 14 days is the minimum follow-up period, although current WHO guidelines mandate at least 28 days (59, 74), and there have been suggestions that the long half-lives of partner drugs in ACTs make it necessary to follow patients for up to 63 days (59). The magnitude of errors likely to be caused by short follow-up periods was not quantified by these authors, and field data are difficult to assess: in principle, drug failures can be distinguished from new infections during the follow-up period by use of molecular typing (75), but in practice this process is hampered by poor detection of minor parasite clones present at the time of treatment (50) which, combined with technical limitations, generates ambiguous results (e.g., see reference 25). The results presented in Table 2 add quantification to the argument about the appropriate length of patient follow-up. There was little difference between apparent cure rates at days 28 and 63 for chloroquine and mefloquine (although, interestingly, many failures occurred after day 63 for CQ), but the results for LF and PQ suggest that following patients for only 28 days would detect fewer than 50% of drug failures and that, consequently, a 63-day follow-up would be extremely beneficial in obtaining good estimates of clinical efficacy.

The robustness of drug regimens faced with changes in individual parameter values is shown in Table 3 and suggests that both CQ and PQ regimens are most sensitive to changes in PK/PD parameters, followed by MQ and, finally, LF; typically, LF could tolerate variations up to 10-fold higher than those with CQ or PQ, implying that it is a much more robust and forgiving drug regimen. It also implies that a mutation(s) would have to encode very large effects to produce an LF drug-resistant phenotype. The addition of artemisinins to these monotherapies increases their power to tolerate variation in PK/PD parameters, helping to protect against the rise and spread of resistance. The size of the protective effect could be large (for example, adding AR to LF more than doubled the LF IC_{50} value at which drug failure occurred [Table 3]) but, in general, was relatively modest in allowing increases on the order of 20 to 30% over the monotherapy value.

Tables 3 and 4 can be used to identify parameters having the largest effects on treatment outcome. The maximal parasite-killing rate (*V*) of the partner drug was consistently ranked as the most important parameter. Unfortunately, it is probably the parameter that is least well estimated (see part 1 of the supplemental material), so better estimates of its magnitude and variance would be valuable data for investigating antimalarial drug action. The parasite intrinsic growth rate constant (*a*) was almost as important as *V* for most treatment regimens modeled. This supports recent suggestions that faster-growing, more virulent parasites may be better able to survive drug treatment (51); interestingly, our results suggest that it is virulence caused by an increased parasite growth rate rather than virulence attributable to high parasitemia that is likely to affect drug sensitivity (see the later discussion about the role of the parasite number present at time of treatment). The relatively low ranks of artemisinin PK/PD parameters in logistic regression analysis (Table 4) are consistent with their stated role as providing protection against resistance to the partner and increasing the speed of resolution of symptoms (reducing PCT) rather than being the primary determinants of treatment outcome.

The analyses shown in Tables 3 and 4 consistently show that the slope of the concentration-effect curve has little effect on treatment outcome. In plotting the concentration-effect curve, the point of inflection is the point on a curve at which the curvature changes sign; this occurs where drug levels have decayed to the IC_{50} value. At drug levels higher than the IC_{50}, the drug effect increases with increasing values of *n*, and at levels lower than the IC_{50}, the drug effect decreases with increasing values of *n*. The results suggest that the two effects cancel out so that the slope has little effect on treatment outcome (we assume that this could be shown algebraically from the above equations but are forced to leave this to researchers more mathematically gifted than ourselves). This observation is important: an influential paper by Shen et al. (53) on HIV treatment argued that drugs with a high value for the slope parameter would be much more effective. Our results suggest that this cannot be extrapolated to antimalarials. The underlying reason is presumably because antiretroviral drugs are taken daily and maintained at high concentrations (when high values of *n* are beneficial), so the penalty paid by a high *n* at low drug concentrations is never incurred. It should be noted that higher values of *n* will increase kill rates at higher concentrations and may therefore be important in rapidly clearing parasites following treatment and hence rapidly alleviating symptoms. We are aware of unpublished suggestions that changes in the value of *n* have been associated with increased resistance, and our analysis suggests that these changes may be an indirect consequence of structural alteration changing the IC_{50} and/or *V* rather than being driven directly by selection on *n*.

An unexpected result from the simulations was the small apparent effect of initial parasitemia on treatment outcome. It was consistently ranked as one of the least important factors influencing outcome (Table 4) and had no significant effect in the clinical trial simulation (Table 5). In contrast, most real clinical trials identify high parasitemia as a strong risk factor for failure. Intuitively, we might also expect an infection of 10^{11} parasites to be about 10 times more difficult to eradicate than an infection of 10^{10} parasites. First note that the OR associated with initial parasitemia is generally around 1.1 to 1.3 (Table 4) and highly significant, so it does have an effect; it is simply not one of the more important parameters. Its nonsignificance in Table 5 can be explained by the simulated clinical trial having only 400 subjects, thereby lacking statistical power to detect the effect. These results suggests that patients generating high drug concentrations and/or infected with more sensitive parasites can reliably clear infections irrespective of initial parasitemia, whereas patients with low drug concentrations and/or more resistant parasites are unable to clear infections with any more than 10^{10} parasites. This argument implies that there is only a small number of patients for whom drug concentrations and parasite drug sensitivity are sufficiently well balanced that initial parasitemia becomes the decisive factor determining outcome. This raises the interesting possibility that the importance of initial parasitemia in clinical trials may not be a direct effect but may be due to confounding with another factor(s). Initial parasitemia is such a good indicator of immunity that we used it as a surrogate for immune status in our simulations (see part 2 of the supplemental material, noting that we allowed immune status and parasite number to be independent variables in the simulations), so it may be that the significance of initial parasitemia is due to its inverse correlation with host immunity. The other plausible confounding factor is that high parasitemia at treatment is often associated with clinical symptoms which may affect drug absorption and metabolism (e.g., see references 36, 42, and 68 to 70).

The “clinical trial” simulations allowed us to assess whether the models produce results consistent with field data obtained in clinical trials. It also allowed us to interrogate the simulated field data to assess the usefulness of output measurements in typical clinical trial analyses because, unlike in real trials, we had access to all parameter values that determined outcome. Using four variables commonly measured in the field, we found both the effects of acquired immunity and the concentration of drug in the patient's serum on day 7 to be significant factors affecting the likelihood of treatment failure. A patient's age and the transmission intensity at time of treatment were both used as indications of the effects of immunity on treatment outcome: transmission intensity was invariably associated with treatment outcome, and the effect of age was also sometimes significant. We ran numerous simulations of these clinical trials, and although the basic patterns were consistent, there were large differences in parameter estimates between simulations (data not shown), even though they were based on random variables taken from the same parameter distributions, and this effect may become more pronounced in real clinical trials, which typically have fewer patients and lower failure rates. There have been suggestions that day 7 serum drug levels be collected as a routine part of antimalarial drug effectiveness trials (71). Again, the detailed data produced by our simulations allow this suggestion to be quantified and examined in more detail. The simulations confirmed that low drug levels (defined here as those below the 15th centile) are associated with increased odds of failing treatment (Table 6) but that, somewhat counterintuitively, their predictive ability, measured as sensitivity, specificity, and area under the ROC curve, was generally poor. The reason becomes obvious by considering the case of DHA-PQ. The failure rate for patients with low drug levels was 16.99%, and that for patients with normal drug levels was 7.22%, giving an OR of 2.62. However, many people with normal levels failed treatment, while many people with low levels were treated successfully. Hence, sensitivity and specificity were both fairly poor, at 29% and 86%, respectively. A practical use of day 7 serum levels is for drug effectiveness surveillance: if low levels are associated with increased risk of failure, then consideration should be given to increasing the drug dosage. Table 6 presents the population attributable risk percentage associated with low drug levels; in principle, this predicts the reduction in failure rates that would occur if drug levels were increased so that no one was subsequently exposed to “low” levels. However, this is a crude method and underestimates the true effect. Increasing drug levels would not just prevent people from receiving low levels but also would increase every patient's drug levels, thereby increasing cure rates among patients receiving “normal” levels (and the possibility that some patients may develop high drug levels, with potential risks of adverse events). It is impossible to quantify this effect from analysis of failure rates in clinical trials and is another instance where quantitative modeling of the type described here can contribute to predicting the benefits of increasing the drug dosage rate.

The ROC curve is a graphical representation of the trade-offs between sensitivity and specificity, and its area under the curve (AUROC) quantifies the predictive power of the variable, with a value of 1 indicating perfect predictive capability and a value of 0.5 indicating no predictive ability. The AUROC values given in Table 6 and the ROC curves presented in Fig. S3.2 in the supplemental material suggest the day 7 serum to be a moderate predictor, at best, of treatment outcome. Our analysis allowed both human PK variation and differential levels of drug sensitivity in parasites. We conjecture that the area under the ROC curve may even provide a clue as to why treatments may be failing. If they are failing due to underdosing, then the AUROC should be high (low drug level is a good predictor of failure), while if parasite drug resistance is the main factor causing failure, then drug levels should be less important determinants of treatment failure and the AUROC should be correspondingly low (as in our simulations). Future work on the simulations will explicitly investigate how much information ROC curves may provide on the etiology of drug failure. In summary, we therefore recommend that clinical trials report not just the odds ratio of failure associated with low drug levels but also an ROC curve analysis.

Increasing tolerance and possible resistance to artemisinins have recently been observed (14, 44, 72), leading to intense speculation about how this will affect the overall effectiveness of ACTs (e.g., see references 6 and 17). An obvious question is how the protection afforded by artemisinins to their partner drugs changes as resistance to artemisinins evolves: is there likely to be a safety margin associated with artemisinins whereby large increases in the IC_{50} can be tolerated before the protective value falls, is a linear fall in protection likely, or as a worst-case scenario, is a rapid decline in protective effect likely to occur? The results shown in Fig. 1 suggest the latter scenario, supporting assertions that measures are urgently required to prevent the evolution and spread of artemisinin resistance (7, 15).

Two issues not directly addressed above are those of synergism and cross-resistance between drugs in a combination (note that the two factors are distinct). Modeling the effect of combination therapies on parasite numbers was achieved by assuming that the two drugs act independently (equation 7). This is likely to be the case for most combinations with ACTs, but other combinations, such as SP and atovaquone-proguanil, may well show synergy. Unfortunately, there is no universally accepted approach for determining synergism and antagonism, and the topic is fraught with controversy and confusion. Greco et al. (24) list no fewer than 13 different methods to determine synergism, and in a comprehensive review, Chou (10) says, “it is hard to find any other field in biomedical science that has more controversy and confusion than drug combinations”; he then cites Goldin and Mantel (22) as giving seven different definitions for synergism, with none of them supporting the others. There appears to be no concise mathematical way of describing synergisms. The best way to incorporate synergism is likely to be the empirical approach taken for SP by Gatton et al. (20), who simply used isobolograms to predict the kill rate for any given concentration combination of the constituent drugs. In summary, incorporating synergy into these models is likely to be problematic, both philosophically and practically, so we do not attempt it here. We do note that this precludes use of this methodology to investigate combinations such as SP and atovaquone-proguanil. Cross-resistance may occur between drugs in a combination even though they act independently; for example, there are concerns that parasites show cross-resistance to mefloquine and artemisinins (2, 40). This effect can readily be included by using covariance terms between the IC_{50} (and/or maximum kill rate) values for each drug. This was not addressed here in the interest of simplicity and to avoid adding additional covariance terms, but it is important to recognize that the assumption that PD parameters are independent can easily be relaxed.

Incorporation of the artemisinin component was the most simplified part of the simulation. The mechanism-based modeling assumes instant absorption and conversion of artemether and artesunate into the active form DHA. As noted above, this can be incorporated using a compartmental model where the compartments are the gut, artesunate/artemether left unconverted in the serum, and DHA in the serum after conversion. The different rates of these processes may explain why the PK/PD parameter estimates differed for the three forms of artemisinins (a fuller discussion of the time course and conversion of artemisinins can be found in the work of Giao and de Vries [21]). We decided against incorporating an explicit compartment model at this stage to minimize model complexity and because calibration of the transfer rates between compartments would be problematic and probably contentious. Instead, we chose to calibrate models separately for each artemisinin variant (Table 1) and to present the results for the partner drugs with alternative artemisinin variants to demonstrate that the results are robust. A second limitation of the model for artemisinins was that it ignored the possibility that parasites enter a drug-induced dormancy stage in which they are unaffected by the drug, as has been suggested to occur for artemisinin (29, 37, 62) and, more recently, for atovaquone (45). Intuitively, it seems unlikely that this will affect the results: artemisinins in 3-day regimens do not clear all parasites, so a small residue persisting in a dormant stage may be negligible, especially as they are likely to recover on a time scale where they are likely to encounter high residual levels of the partner drug. Furthermore, ignoring dormancy leads to overestimating the impact of the artemisinin component in clearing infection, and the results shown above suggest that even an overestimated impact is secondary to the role played by the partner drugs (Table 4). We note that our approach of daily updating the parasite load and drug level was explicitly designed to make the calculations highly flexible, and in principle, we could incorporate the effects of dormancy by augmenting the parasites present each day with those predicted to be exiting the dormancy stage. In summary, it would be possible to make the extension into compartmental models and dormancy, but we leave this to future work.

A second way in which this mechanism-based modeling approach could be developed usefully for malaria is to make the simulations more specific to drugs and their human subjects. Drug absorption, distribution, metabolism, and excretion can differ substantially in young children or infants, pregnant women, patients with severe disease, and those with HIV/AIDS coinfection (4). For ethical reasons, most clinical trials are performed on nonpregnant adults with uncomplicated malaria and no comorbidities, so it seems likely that, at least in the first instance, the impact of these factors on treatment outcome may be addressed initially by a well-constructed and calibrated PK/PD model. Similarly, we incorporated variation in PK/PD by assuming a CV of 30% across all parameters. Some parameters are likely to be much more variable (for example, Mu et al. [40] reported IC_{50}s over a 100-fold range), while others (possibly *V*?) may be much less variable. Incorporating parameter values and associated distributions specific to individual drugs would give us confidence to extend the results into the quantitative domain useful for policy makers. One notable example is the use of fixed treatment dosages based on age or height bands, where there will be considerable variation in drug concentrations within groups as a consequence of variation in body weight: dosages have to reliably cure all people within the band (while avoiding toxic concentrations), and it is not immediately obvious how to identify the appropriate dosage for each band or, consequently, how many bands would be required.

We do not suggest that pharmacological modeling of antimalarial treatment will ever replace the gold standard of clinical trials, but it does appear capable of generating results that are entirely consistent with field observations. It has key advantages of speed, the ability to generate large data sets, the ability to rapidly compare different scenarios, and freedom from ethical restrictions on investigating factors such as poor compliance. The value of modeling is that it can take arguments that are predominantly verbal into a more explicit, quantitative domain; our objectives here were to improve our understanding of how antimalarial drugs act in general, what conclusions can be drawn about, for example, the impact of increasing levels of artemisinin resistance, and what data collected in the field may reveal about ACTs. It seemed reasonable to apply mechanism-based PK/PD modeling to this problem. Future work will develop the methodology in more specific directions. In particular, we will incorporate the absorption, conversion, and distribution phases of the drugs. These are important for artemisinins, where absorption lag times may be significant compared to their half-lives and where conversion of artesunate and artemether to DHA may be relatively rapid. The absorption and distribution phases also determine peak serum concentrations of the partner drugs, which may be important determinants of potential toxicity; this becomes important in designing fixed-dose regimens based on age, weight, or height bands because dosages per kg may vary widely within a band. We also extended the mathematics to investigate drugs given as infusions (such as intravenous quinine) but have not presented this methodology in the interest of brevity. Finally, it would be informative to include the possibility that parasites may enter dormant stages where they are unaffected by drugs such as the artemisinins (29, 37, 62) and atovaquone (45). Meanwhile, we conclude that initial analyses of antimalarial PK/PD models are encouraging and qualitatively improve our understanding of how antimalarial PK/PD factors combine to determine treatment outcomes, and we await future developments with interest.

## ACKNOWLEDGMENTS

We thank Tiago Antao, Melissa Penny, and Steve Ward for help and advice in developing and calibrating the calculations. We also thank the three reviewers for their helpful comments and suggestions.

This work was supported by the Bill and Melinda Gates Foundation (grant 37999.01), the Swiss Tropical and Public Health Institute, and the Liverpool School of Tropical Medicine.

## FOOTNOTES

- Received 9 December 2010.
- Returned for modification 4 March 2011.
- Accepted 13 April 2011.
- Accepted manuscript posted online 2 May 2011.
↵† Supplemental material for this article may be found at http://dx.doi.org/10.1128/AAC.01712-10.

- Copyright © 2011, American Society for Microbiology