Population Pharmacokinetics of Intravenous Artesunate: A Pooled Analysis of Individual Data From Patients With Severe Malaria

There are ~660,000 deaths from severe malaria each year. Intravenous artesunate (i.v. ARS) is the first-line treatment in adults and children. To optimize the dosing regimen of i.v. ARS, the largest pooled population pharmacokinetic study to date of the active metabolite dihydroartemisinin (DHA) was performed. The pooled dataset consisted of 71 adults and 195 children with severe malaria, with a mixture of sparse and rich sampling within the first 12 h after drug administration. A one-compartment model described the population pharmacokinetics of DHA adequately. Body weight had the greatest impact on DHA pharmacokinetics, resulting in lower DHA exposure for smaller children (6–10 kg) than adults. Post hoc estimates of DHA exposure were not significantly associated with parasitological outcomes. Comparable DHA exposure in smaller children and adults after i.v. ARS was achieved under a dose modification for intramuscular ARS proposed in a separate analysis of children.

combined additive and proportional) was determined for each dataset. Before the PK data from the six studies were pooled each PK dataset, as well as PK datasets from each period of the cross-over trials, was analysed separately to determine whether the same structural model and error model is appropriate for all datasets.
Structural model selection. A NLME model with either a one-compartment or twocompartment IV bolus structural model and additive error on the natural log-scale were fitted to each dataset. Initially both NLME models were fitted by first order conditional estimation (FOCE) with interaction and assuming no below quantification limit (BQL) samples. In order to fit the NLME with a two-compartment structural model to each dataset, no covariance between the individual parameters was estimated and the variance of the following individual parameters was set to zero (i.e. assumed not to vary between patients) for the following studies:  variance of the inter-compartmental clearance was fixed at zero for all studies; Based on the visual predictive checks (VPCs), it was decided that the simpler onecompartment structural model captures the central tendency of the observations adequately and would be appropriate for all datasets.
The fit of the IV bolus one compartment model (which assumes instantaneous absorption) was then compared to a one compartment model with zero order absorption for each dataset.
In order to fit the one compartment model with zero order absorption, absorption duration was assumed not to vary between patients and only the population average absorption parameter could be estimated. For all studies except Kremsner et al. 2012 2 , either the standard error of the population pharmacokinetic (PK) parameter estimates could not be calculated if the population average absorption parameter was included in the model or the standard error of the population mean absorption duration parameter was very small (indicating lack of convergence of NONMEM's covariance step). The estimate of the absorption duration was also very sensitive to the initial value. The objective function value was similar for both models or tended to be greater for the zero absorption model and the population mean clearance and volume of distribution estimates were similar. Due to these fitting issues and the similarity of the population PK parameter estimates from the zero absorption model to the IV bolus model it was decided not to pursue to the zero absorption model.

Error model selection.
Next NLME models with a one-compartment structural model and either an additive error on the natural log-scale or combined error model were fitted to each PK dataset. The additive error on the natural log-scale model was assumed appropriate for all datasets based on the VPCs.
Allometric scaling and BQL samples. Allometric scaling of the population mean PK parameters to be for a patient of median weight (15kg) further improved model fit according to the objective function value and VPCs. To examine if modelling below quantification limit (BQL) samples influences parameter estimation the NLME model with one-compartment structural model, exponential error model and body weight as an allometric function on the population mean clearance (CL) and volume of distribution (V) was fitted using the M3 method. The M3 method maximises the likelihood for the data above the limit of quantification (LOQ) and treats BQL concentrations as censored. The number of BQL concentrations for each study is given in Supplementary Table 4, along with the assay method and LOQ for the assay. The population mean PK parameter estimates derived by FOCE with interaction without any assumed BQL samples were similar to those derived by the M3 method for BQL samples and consequently, the VPCs were comparable to those produced using the parameter estimates derived from the simplified estimation procedure, which assumes no BQL samples.
Meta-analytic approach. The meta-analytic approach involved comparing the population mean CL and V estimates for each study and period of the cross-over trials using a forest plot (Supplementary Figure 2). The population mean CL and V estimates plotted in Supplementary Figure 2 were derived by fitting a NLME model with one-compartment structural model, exponential error and body weight as an allometric function on the population mean CL and V parameters using the M3 method for BQL concentrations.

Supplementary Information 3
A 3-level Bayesian hierarchical model was fitted to the pooled pharmacokinetic (PK) data, which includes an extra level to account for between-study variability, in addition to the nonlinear mixed-effects (NLME) model which characterizes the within-patient processes of distribution and elimination, and the degree to which these processes vary between patients.
The fit of this model was compared to the fit of a 2-level model including study as a fixed effect on both population mean clearance (CL) and volume of distribution (V) parameters, and a 2-level model that did not account for between-study variability. For each model, the residual variability was allowed to vary with study. The mathematical forms of the three hierarchical models fitted to the pooled PK data are given in Supplementary Table 5.
Bayesian estimation: In Bayesian statistics parameter estimates and interval estimates for each NLME model parameter are derived from the posterior distribution, which is the probability distribution of the NLME model parameters (population PK, individual PK and betweensubject variability (BSV) parameters) given the data. if the model fits the data well then the observed percentiles should be contained in the 95% credible intervals for the corresponding predicted percentile). To generate replicated datasets from the posterior predictive distribution for the 3-level model, the 1000 parameter values sampled from the posterior distribution for the population mean clearance (CL) and volume of distribution (V) (averaged across studies) and between-study variability ((along with each patient's covariate information and sampling times) were used to simulate population PK parameters for each study (CL k and V k ), individual PK parameters (CL ki and V ki ) and concentrations from level 1 of the 3-level hierarchical model (see Supplementary Table 5 for   model definitions and corresponding table footnotes for parameter definitions).

Supplementary Information 4
Inclusion of covariates: The mathematical form of the individual clearance (CL i ) and volume (V i ) functions are given in equation (1) and (2), respectively.
In equations (1) and (2) Stepwise covariate selection: For each subgroup analysis, potential covariates were investigated by adding covariates to their respective base models using a stepwise forward addition (p-value of 0.05) and backward elimination (p-value of 0.01) approach. Continuous covariates were included in equations (1) and (2) as linear functions, whose form is given in The assumption of linearity was examined using plots of the post-hoc estimates versus covariate values and also with the likelihood ratio test.
In equations (3),  i represents either the individual CL or V parameter for patient i;  pop,i represents either the population mean clearance or volume of distribution after scaling by body weight and age (see equations (1) and (2)); if COV i is a continuous covariate centred at its median (medianCOV), then  , COV is the fractional change in  with each unit change from the median covariate value; if COV i is cateogrical, then  , COV is the fractional change in  from the reference category to the non-reference categories; ,i is the IIV for . Covariates found to be statistically significantly associated with the population mean PK parameters in Subgroup 1, were added to the base model for the stepwise covariate selection procedure in the subsequent subgroup analyses. The fit of the final models was examined using VPCs.

Supplementary Tables
Supplementary    Table S5: Mathematical form of hierarchical models fitted to the pooled PK data to investigate whether modeling between-study differences in the population mean PK parameters influenced the predictive properties of the model.