An improved pharmacometric model for artesunate treatment of falciparum malaria

Sompob Saralamba (  sompob@tropmedres.ac ) Mahidol Oxford Tropical Medicine Research Unit, Mahidol University Julie A Simpson Centre for Epidemiology and Biostatistics, Melbourne School of Population and Global Health, The University of Melbourne https://orcid.org/0000-0002-2660-2013 Noppon Choosri Center of Data Analytics and Knowledge Synthesis for Healthcare, Chiang Mai University Lisa White Big Data Institute, Li Ka Shing Centre for Health Information and Discovery, Nu eld Department of Medicine, University of Oxford https://orcid.org/0000-0002-6523-185X Wirichada Pan-Ngum ) Department of Tropical Hygiene, Faculty of Tropical Medicine, Mahidol University https://orcid.org/0000-0002-9839-5359 Arjen M Dondorp Mahidol Oxford Tropical Medicine Research Unit, Faculty of Tropical Medicine, Mahidol University Nicholas J White (  nickwdt@tropmedres.ac ) Mahidol Oxford Tropical Medicine Research Unit, Faculty of Tropical Medicine, Mahidol University


Introduction
The artemisinin derivatives are the cornerstone of current antimalarial therapies. These drugs are highly effective in killing both circulating and sequestered malaria parasites but are eliminated very rapidly from the body. The elimination half-life of the common biologically active metabolite, dihydroartemisinin (DHA), is usually reported as less than one hour. Despite this, in most regimens, the artemisinin derivatives are given only once daily. Simple pharmacometric models, in which parasite killing is a direct function of plasma concentration, predict that with greater exposure to the artemisinins there should be increased parasite killing. This was one of the motivations behind the development of more slowly eliminated synthetic peroxide antimalarial drugs, and was proposed as a solution to the reduced parasite killing associated with artemisinin resistance. However, administering artemisinins more than once daily does not accelerate parasite clearance or improve therapeutic responses [1][2][3] . It has been suggested that parasite clearance lags behind killing of malaria parasites, and so damaged or dead parasites accumulate relative to live parasites in the circulation 4 . From this hypothesis it was claimed that parasite clearance rate is not a suitable indicator of artemisinin resistance, and that the failure of split dosing to accelerate parasite clearance could still obscure a bene t in cure rates 4 . This was addressed in a metaanalysis of monotherapy trials which showed that split dosing was not associated with higher cure rates 5,6 . Thus, the results of clinical trials do not support this hypothesis, and current pharmacometric models do not explain the observed therapeutic responses satisfactorily.
Taken together these observations are not compatible with simple concentration-dependent malaria parasite killing by artemisinins and current mathematical models do not explain satisfactorily the in-vivo malaria parasite density dynamics following treatment with artemisinin derivatives. Data suggest that potentially reversible parasite injury and delayed parasite stress responses after exposure to the drug might contribute to more complex pharmacodynamics [7][8][9][10][11][12][13][14][15] . One explanation for the failure of artemisinin monotherapies to achieve 100% cure rates, even in artemisinin sensitive infections, is the persistence of temporarily drug insensitive or "dormant" parasites 8, 11,15,16 .
In the current study we evaluate the hypothesis that, following exposure to artemisinin drugs, parasites can be damaged and rendered temporarily refractory to further injury, and that a fraction of these parasites could recover afterwards. We assessed whether applying this hypothesis in a mathematical model could describe satisfactorily parasite clearance dynamics after treatment with artemisinin derivatives.

Results
The measured DHA drug concentrations of patients given 2 mg/kg of oral artesunate daily were rst analysed to estimate individual speci c pharmacokinetic parameters. Supplementary Table S1 and S2 summarises the estimated pharmacokinetic parameters for DHA for the 39 patients. A sequential pharmacokinetic-pharmacodynamic modelling approach was then performed where DHA concentrations were simulated for each individual patient at the corresponding times of the parasitaemia measurements.

Pharmacodynamic model
The proposed pharmacodynamic model was tted to the parasite clearance pro les from patients treated with artesunate monotherapy, using a Bayesian hierarchical model 17 Table 2.
The clinical study 18 was conducted in two areas, one where artemisinin resistant parasites were prevalent (Pailin, western Cambodia) and the other where parasites were still sensitive to artemisinin (Wang Pha, western border of Thailand). The population means of estimated recovery and death rates of each parasite asexual stage showed that parasite recovery rates were slightly different between stages and sites. Recovery rates from Pailin, western Cambodia were slightly lower than the recovery rates from Wang Pha, western Thailand. As expected the estimated death rates of ring stage parasites in Pailin were substantially (over ten fold) lower than those of Wang Pha ( Table 2). The estimated death rates for trophozoite and schizont stage parasites were not different between the locations. This is supported by concomitant standard 48 hour in-vitro susceptibility studies which showed little difference in artesunate susceptibility between the sites 18 . The box-and-whisker plots of recovery and death rates for ring, trophozoite and schizont stage parasites compared between the clinical sites are shown in Fig. 3. The population means of the recovery lag-times for damaged parasites to recover to normal state were 2.42 (0.04-7.00) hours in Wang Pha and 2.01 (0.04-7.32) hours in Pailin.

Split dose simulation
For the split-dose simulation, the DHA concentration pro les were incorporated in the model every 12 hours for 7 days (i.e. a dosing regimen of 2 mg/kg oral artesunate every 12 hours for 7 days) and parasite clearance was determined from the corresponding simulated parasitaemia pro les (Fig. 4). Figure 5 compares the parasite clearance pro les of the daily and split dose (every 12 hours) artesunate dosing regimens. These distributions of parasite clearance show that, within the same clinical site, the median and spread of parasite clearance following daily and 12 hourly drug administration are not signi cantly different ( Fig. 6 and Table S3).

Discussion
Mathematical models describing the pharmacometric properties of the artemisinin antimalarials have assumed a simple and direct relationship between drug exposure and malaria parasite killing. Because these drugs are eliminated very rapidly these models naturally predicted that more sustained exposure, created by slowing drug clearance, administering by constant infusion or by frequent dosing, would enhance parasite clearance, and thus improve therapeutic responses 4,19,20 . This was proposed as a solution to the challenge posed by artemisinin resistant P. falciparum infections. However, clinical studies did not con rm these predictions 3 and there was no signi cant advantage in terms of parasite clearance or cure rate from giving the artemisinin derivatives more than once daily despite their rapid elimination.
The simple models thus appear to be wrong.
Artemisinin resistance in P. falciparum is characterised by a reduced parasite clearance rate in vivo 18 .
This re ects reduced susceptibility of the circulating ring stage parasites to the drug 21 . Artemisinin resistance in eld isolates is associated with mutations in the propeller region of the Pfkelch gene on chromosome 13 (K13). This causal association has been con rmed in transfection studies, although the contribution of other genes (often described collectively as the genetic background) is substantial 22 .
Artemisinin resistance is thought to involve altered parasite cellular responses rather than receptor or transporter alterations, but the exact mechanism remains unclear 23 . As the more mature stages of K13 mutant P. falciparum isolates remain sensitive to the artemisinin derivatives the drugs do remain e cacious in clinical practice, but they kill less parasites per asexual cycle and so the therapeutic responses are diminished. Longer exposures from longer courses over three or four asexual cycles improve therapeutic responses in artemisinin resistant infections, but giving the drugs more frequently than once daily does not 1-3 .
In the present study, a new mathematical model expanding the earlier simpler models (20,21) was evaluated. This model was constructed to account for the dose-response relationship and could explain the lack of additional effect by more frequent administration. The proposed model incorporated the hypothesis that, in the presence of artemisinins, some parasites are injured, stop growing and become temporarily insensitive to the drug. These parasites can, however, recover and return to their normal, drugsensitive state. Whether this unresponsive ("injury") state causing a temporary arrest in development is the same or a similar process as the well described "dormancy" phenomenon has not been speci ed, and both processes can explain recrudescences following standard treatments is not speci ed 9 . The model was tted to parasite count data obtained from patients in western Cambodia (where artemisinin resistance was prevalent) and western Thailand (before the main emergence of artemisinin resistance there) 18 . The model captured satisfactorily the dynamics of P. falciparum parasites in patients receiving artesunate monotherapy every 24 hours for 7 days. As reported previously 19 the major difference between parasites from Pailin, which were artemisinin resistant, and those from Wang Pha, which were sensitive, was the more than ten fold estimated reduction in ring stage killing. Trophozoite stage killing was similar and schizont killing was reduced by about 45%. Recovery rates were slightly but not signi cantly different between the two sites. This suggests, if this hypothesis is correct, that recovery rates from parasite injury are not greatly affected by the mechanisms involved in artemisinin resistance.
The lower predicted death rates of the asexual ring stage parasites in the artemisinin resistant infections 19 is supported by extensive experimental investigations 7-12, 21,24,25 and the strong correlation between speci c ring stage in-vitro susceptibility evaluations, K13 mutations and slow parasite clearance. The unresponsive ("injury") state postulated in this model suggested that the parasites were in a damaged state between 2 and 50 hours prior to dying. The estimated parasite recovery rate from this model was about 1 in 740-850 unresponsive parasites per 24 hours. The estimate of the recovery lag-time was 2-2.4 hours suggesting that the parasites became damaged or dormant shortly after being exposed to the drug.
This model was developed to account for the failure of frequent dosing regimens to accelerate parasite clearance or enhance cure rates in artemisinin containing antimalarial drug regimens. It did this satisfactorily but there are several limitations to this modelling exercise. Many of the pharmacodynamic parameter values in this model have not been measured directly so the system is unidenti able. It has also simpli ed the complex relationship between parasite stage of development, and time and intensity of drug exposure, and assumed homogeneous parasite stage distributions and multiplication and elimination kinetics. Although the model can reproduce the observed data, this does not mean that it has explained the underlying biology (i.e. the model may not be correct). Nevertheless, it does represent a relatively simple hypothesis that is consistent with observations and is testable. But there may well be other hypotheses whose models would t equally well with the data. We propose that, as a minimum, such models should be capable of reproducing both delayed clearance and the unchanged clearance rates under frequent dosing in order to be consistent with observed data.
In conclusion, a new within-host pharmacometric model is proposed, which supports the hypothesis that parasites enter a temporary drug refractory injury state after contact with artemisinin antimalarials, which is followed by delayed death or reactivation. The model tted the observed sequential parasitaemia data from patients with artemisinin resistant and sensitive P. falciparum infections and con rmed the known dose-response relationship for reduced ring stage activity in artemisinin resistant infections.

Parasite clearance data and pharmacokinetic data
We re-evaluated data used to generate an earlier simpler mathematical model in which there was a direct concentration-dependent malaria parasite killing (18). The data used in this evaluation were serial parasite count data (asexual parasites) and drug concentration data (dihydroartemisinin (DHA)) from 20 patients from Pailin, western Cambodia and 19 patients from Wang Pha, western Thailand. These patients were treated with artesunate monotherapy as part of a clinical trial to investigate artemisinin resistance conducted during 2007 and 2008 by Dondorp et al. 18 . Patients received 2 mg/kg oral artesunate every 24 h for 7 days. Parasite counts were determined by microscopy at 0, 4, 8 and 12 h, and then every 6 h until two consecutive negative slides were recorded. DHA concentrations were measured at 0, 0.25, 0.5, 1, 1.5, 2, 3, 4, 5, 6, 8 and 12 h following the rst dose.

Mathematical model
The base model is that of White et al. describing the within-patient dynamics of P. falciparum parasites 19,25,26 . The age distribution of the parasites at any time before treatment was assumed to be unimodal and Gaussian. The asexual life cycle was assumed to be 48 hours. This distribution shifts to the right as parasites become older. The infection multiplication factor is the population average number of merozoites which successfully infect new red blood cells (RBCs) per one infected RBC. During the expansion phase of the infection this is a positive number, limited by the average number of merozoites per schizont. Time in the model was discretised and xed at 1 hour intervals. Asexual P. falciparum parasites were therefore divided into 48 individual age-groups, from 1 to 48 hours. As the parasites become older by 1 hour the number of parasites moves to the next age-group, until at age 47-48 hours at schizogony, after which they become ring stage parasites with aged 0-1 hour at a density de ned by the multiplication factor.

Antimalarial pharmacodynamics
Antimalarial drugs kill malaria parasites, but the pharmacodynamic effect depends both on the plasma concentration of drug and the susceptibility of the predominant stage of parasite development. When the antimalarial drug (in this case an artemisinin derivative) is present, a fraction of parasites at each susceptible stage will be damaged and some of these die (and subsequently "pitted" from the infected RBC by the spleen). It is hypothesised that a fraction of the remainder of the damaged parasites can recover potentially, but their development is arrested, during which period the parasites are refractory to further injury by the antimalarial drug. In the presence of the drug, the fraction of the damaged parasites ( f(t)) at time t is calculated by: where E m is the maximum extent (%) that the parasites can be damaged by the drug, c(t)is the concentration of DHA at time t hours, Ec 50 is the concentration that gives 50% of the damage effect, and γis the slope constant. Damaged parasites that recover after a lag time return to their pre-treatment stage, with normal aging, multiplication and pre-treatment drug susceptibility. Parasites die at a rate of ν (/hour) and surviving parasites will produce new merozoites. The fraction of damaged parasites that recover to their pre-treatment normal state at time t, κ(t), is calculated from where r m is the maximum fraction of the damaged parasites that can recover, f(t) is the fraction of the damaged parasites at time t, and lag is the recovery lag time. These recovery and death fractions are assumed to depend on the stage of development of the parasites. The susceptibility of malaria parasites to all antimalarial drugs is strongly dependent on their stage of development. For simplicity we divide the parasites into three developmental stages; rings (0-26 hours), trophozoites (27-38 hours) and schizonts (39-48 hours). Laboratory studies suggest that young ring stages are most likely to enter a state of dormancy, particularly following exposure to artemisinin derivatives 8, 10,11,15 . A diagram of the model is shown in Fig. 1. The plasma concentration of DHA, c(t), following each oral dose was modelled as follows: The tting was implemented in Stan using the Hamiltonian Monte Carlo (HMC) method 17 . The likelihood function was derived assuming that log10 transformed of each observed parasite density measurement ( Φ i ) at time i was sampled from a normal distribution, of which the mean was the model output (M i ) for that observed time i and its standard deviation was ϱ. That is, the likelihood can be written as: The rst 1,000 parameter values sampled for each chain were discarded as burn-in and the subsequent 1,000 samples from each chain (n = 3) were used to estimate the posterior distributions and also the posterior predictive checks. The convergence of the chains of each parameter was assessed using trace plots and the R statistic 31 . The prior distribution of each model parameter was the uniform distribution with the boundaries as shown in Table 1. The Stan and C# code of the model can be downloaded from https://github.com/slphyx/DamagedParasites.

Split-dose simulation
To predict the effect of administering artesunate to patients every 12 hours for 7 days, the models were re-run using parameter values sampled from the posterior distributions estimated from tting the model to parasite clearance pro les following artesunate monotherapy once daily, except that drug concentration pro les were changed from every 24 hours to every 12 hours for 7 days. For each hypothetical patient, parasite clearance times were calculated from the parasite versus time pro les which were sampled from the tted hourly pro les of each patient. Here the parasite clearance time was de ned as the time from rst dose of artesunate to the time that the simulated circulating parasite count was below the detection limit (L), which differed between the hypothetical patients and was estimated from: L ≈ l × 80,000 × weight 5 where l is the lowest parasite count number (parasites/µL), and weight is the recorded body weight of the patient (kg) 19 . All patients were assumed to have a total blood volume 80 mL/kg.
To test whether the model could reproduce the trial results from Das et al. 3 , which showed no improvement in parasite clearance time whether patients were given 2 mg/kg artesunate every 12 or 24 hours for 7 days, clearance times predicted by simulating both dose regimens were compared visually.   Figure 1 A diagram of the proposed model, where f(t) is the fraction of sensitive parasites at time t that can be damaged, κ is the recovery rate of each asexual stage, and ν is the death rate of each asexual stage.  Box-and-whisker plots generated from the posterior distributions of estimated population means of the asexual stage speci c death rates (top) and recovery rates (bottom) for each study site (purple = Pailin, western Cambodia; blue = Wang Pha, western Thailand).   Comparison of distribution of parasite clearance times (hours) derived from simulated parasitaemia pro les for a single dose of artesunate every 24 hours (thick lines) versus every 12 hours (dashed lines) for 7 days from Pailin, western Cambodian (Purple) and Wang Pha, western Thailand (Blue).

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. Supplementary.pdf