A Dynamic Stress Model Explains the Delayed Drug Effect in Artemisinin Treatment of Plasmodium falciparum

Construction of the dynamic stress model.In order to identify the key features that motivate the development of our model, we first review the in vitro experimental procedure (see reference 7 for details). Cultures containing equal quantities of tightly age-synchronized P. falciparum parasites (3D7 laboratory strain; over 80% of parasites synchronized within a 1-h age window) were treated with a specified dose of dihydroartemisinin (DHA) for a duration of 1, 2, 4, or 6 h before washing (to remove all drug). To quantify the effect of drug, a viability assay was performed. Viable parasites were defined as those able to reproduce and enter the next cycle of replication (thus excluding dead and dormant populations), assessed by measuring the parasitemia (P) in the trophozoite stage in the following life cycle 48 h later. In order to calculate the viability, parasitemia was also measured for two special cases: the control case (P_control), where no drug was applied, and the background case (P_background) with supermaximal DHA concentration (>10× the 50% lethal dose of 3 days, in nanomolars), applied for more than 48 h. Viability (V), a unitless ratio, then was given by subtracting the unviable population, V=P−PbackgroundPcontrol−Pbackground(2) To study stage-specific drug effects, Klonis et al. tested four different parasite ages (by using different age-synchronized groups): 2 h postinfection (h p.i.; early ring stage), 7.5 h p.i. (mid-ring stage), 24 h p.i. (early trophozoite stage), and 34 h p.i. (late trophozoite stage). Two examples of viability data are given in Fig. 1 and show the viability for different durations of drug exposure (1 h, 2 h, 4 h, and 6 h) with an initial DHA concentration of approximately 39 nM or 300 nM. Note that DHA concentration also decays in vitro with a half-life of approximately 8 h (which is much longer than that which occurs in vivo; see Discussion). Experiments were performed in technical replicates for each combination of initial DHA concentration and drug exposure duration.

• Open in new tab
• Download powerpoint

FIG 1

Representative experimental data showing how the fraction of viable parasite (i.e., viability) changes with the duration of drug exposure for two different initial DHA concentrations (39 nM [left] and 300 nM [right]) and four different parasite life stages. For one parasite stage and one drug exposure time in each panel, duplicate viability measurements are provided, which means in total 40 data points are measured in each panel (note that viability for zero exposure time is always equal to 1 due to normalizing parasitemia to itself). Insets indicate the in vitro decay of DHA concentration. Empty circles display the raw viability data, and the curves pass through the arithmetic means of the paired data points. (Data are sourced from reference 7.)

We take as our fundamental conceptualization of antimalarial action that the drug kills or otherwise prevents parasites within infected red blood cells (iRBC) from being able to produce viable merozoites (which would go on to invade and infect other RBC at the end of the first life cycle). We model the number of iRBC of age a (i.e., RBC that have been infected with parasites for a hours) surviving drug exposure [N(a,t)] by the first-order partial differential equation ∂N(a,t)∂t+∂N(a,t)∂a=−kN(a,t)(3) where k is a drug-induced parasite killing rate and may depend on other factors, such as drug concentration, parasite life stage, or even drug exposure duration (which will be explicitly indicated once we formally introduce those dependencies later). We have the boundary condition N(0,t) = rN(48,t), where r is the parasite multiplication factor, indicating the average number of newly infected RBC generated from merozoites released from a single iRBC at the end of the preceding life cycle.

The in vitro experiments use tightly age-synchronized parasites, allowing for further simplification to an ordinary differential equation system, which is sufficient for determining the time dependency in the maximum killing rate (k_max) and the half-maximal killing concentration (K_c). We track only the number of newly infected RBC generated from parasites first exposed to drug at age ā [denoted by N̄(t)]: dN¯(t)dt=−kN¯(t)(4)

As mentioned in the Introduction, the parasite killing rate, k, is empirically modeled by a Hill function of drug concentration, C(t), k(C(t))=kmaxC(t)γKcγ+C(t)γ(5) where k_max is the maximum killing rate, K_c indicates the drug concentration at which half-maximal killing (k_max/2) is achieved and γ is the Hill coefficient. To capture the time-dependent features of the in vitro data, we generalized the model by allowing k_max and K_c to be dependent on the duration of drug exposure. We considered the time variation to be a function of an auxiliary modulatory variable S(t), which we refer to as a general cell stress. During drug exposure, parasites develop a stress response, the extent of which determines the killing effect (and thus the concentration-killing rate function; equation 5). The stress, S(t), is normalized to vary between 0 and 1 (inclusive). We consider S to increase in the presence of drug above some (very small) threshold level, C*, but decrease once drug concentration, C, is below C*. For the increase phase, we apply a simple first-order differential equation: dSdt=λ(1−S)(6) where λ is a rate constant which sets the time scale for stress development. In the absence of additional experimental data, S(t) is assumed to immediately reset to zero once drug concentration falls below C*. While this is sufficient to capture all available in vitro data, we anticipate that further experimental research will allow us to more closely tie empirical determinations of the mechanisms of stress and its accumulation to our modulatory variable, S (with consequential changes to equation 6).

We then assume that k_max and K_c depend on the development of the stress response, i.e., they are functions of S. This construction allows for k_max and K_c to vary with different drug exposure times, mediated through the stress response. k_max is evidently positively correlated with S and K_c is negatively correlated, indicating that as stress accumulates, the ability of the drug (at a given concentration) to kill parasites increases. In the absence of detailed experimental data, we assume these relationships are linear: kmax=αS(7) and Kc=β1(1−S)+β2(8) where α, β₁, and β₂ are parameters to be determined. α can be viewed as the maximum achievable killing rate (i.e., k = α), achieved only after long exposure to a high drug concentration. β₁ affects how sensitive the killing rate, k, is to a change in the stress variable, S, and β₂ represents the half-maximal activation concentration once stress is at a maximum (i.e., when S = 1). Note that to avoid issues of parameter identifiability as λ→∞, we do not allow for an additional constant term in the expression for k_max. Similarly, and as explored further in Discussion, without direct experimental evidence for a particular underlying mechanism for stress, we have assumed that a single stress response is able to influence both k_max and K_c.

Under this simple formulation of the model for stress accumulation and the linear relationship between stress and killing, we can solve equation 6 to obtain S(t)=1−e−λt(9) and thus kmax=α(1−e−λt)(10) and Kc=β1e−λt+β2(11) Of note, when λ→∞ (i.e., the modulatory variable reaches its steady state instantaneously), k_max = α and K_c = β₂ such that our model reduces to a traditional PK-PD model with a fixed relationship between drug concentration and killing rate.

For finite λ, k_max and K_c become functions of S and so duration of exposure. In particular, for low λ, our model displays a slow development of the stress response and thus is capable of capturing a delayed reduction in viability, indicating its suitability for the in vitro data shown in Fig. 2 and detailed in references 7 and 9. For simplicity, in this paper we refer to the traditional killing rate model with constant (although perhaps stage-dependent) k_max and K_c as the stationary model and refer to our generalized model as a dynamic stress model.

• Open in new tab
• Download powerpoint

FIG 2

Results of fitting the model to viability data (early ring stage). The initially applied DHA concentration is indicated for each panel. Empty circles (appearing in duplicate) are the repeated measures of viability by (initial) drug concentration and exposure duration. For one DHA concentration and one drug exposure time in each panel, duplicate viability measurements are provided, which means in total 10 data points are measured in each panel (note that viability for zero exposure time is always equal to 1 due to normalizing parasitemia to itself). Black curves show the predicted mean viability measurements from the model with fixed γ parameter. Red dashed lines are the 95% confidence intervals (CI) for the predicted mean viability measurements (derived using simulation-estimation of 500 concentration-effect profiles and parametric bootstrap CIs), and blue-shaded regions are 95% prediction intervals (PI; derived 2.5th and 97.5th percentiles of 500 simulated concentration-effect profiles) for a new viability measurement if it were generated under the same experimental conditions (i.e., drug concentration and pulse duration).

The dynamic stress model contains five parameters (λ, γ, α, β₁, and β₂) to be determined by fitting to available data. These five parameters are assumed to be stage specific, and thus the model is fitted separately to the viability data for each parasite stage. Note that by incorporating the phenomenological model of stress through the modulatory variable S, we aim to develop a better understanding of how the killing rate evolves in the presence of drug. While we are as yet unable to explore the underlying mechanisms governing the development of both stress and drug action (e.g., the changes at the cellular or even molecular levels), we discuss possible biological interpretations further in Discussion.

For modeling the in vitro experiments of tightly age-synchronized parasites, we use equation 4. The solution to equation 4, subject to an initial condition N̄(0) = N̄₀, is N¯(t)=N0¯e−∫0tk(C(τ),S(τ))dτ(12) where we have explicitly presented the killing rate, k, as a function of DHA concentration, C, and stress, S, both of which are functions of time (following initiation of drug exposure). DHA concentration, C(t), as a function of time due to in vitro decay, is given by C(t)=C0e−ln(2)t1/2t(13) where C₀ is the initial dug concentration and the in vitro half-life of DHA (t_1/2) was measured to be about 8 h (21).

For a drug pulse with a duration of t_d hours in a given stage of the parasite life cycle, the total number of iRBC, N_d, at the time of data collection (during the trophozoite stage in the next life cycle) is given by Nd=rN0¯e−∫0tdkdτ+N0¯(1−e−∫0tdkdτ)(14) where k's dependence on C(τ) and S(τ) is now implicit. The first term represents the number of iRBC with live parasites (having expanded by factor r, the parasite multiplication factor), while the second term represents the number of nonviable parasites. For the control case (no drug), the number of parasites, N_c, is given by equation 14 with k = 0 (N_c = rN̅0̅). For the background case (all parasites killed due to supermaximal exposure), the number of parasites, N_b, is given by equation 14 with k→∞ (N_b = N̅0̅). Substituting these two expressions back into equation 14, we have Nd=Nce−∫0tdkdτ+Nb(1−e−∫0tdkdτ)(15) which can be rearranged to give e−∫0tdkdτ=Nd−NbNc−Nb(16) The right-hand side is precisely the parasite viability (V) as defined in the in vitro experiments (7, 9). Therefore, we have V(C0,td)=e−∫0tdk(C(τ),S(τ))dτ(17) where the dependence of V on the initial DHA concentration, C₀, is established through equation 13. Furthermore, we identify ϕ=∫0tdk(C(τ),S(τ))dτ as the cumulative stress effect. We have identified the parasite viability as a function of the cumulative stress effect arising from a dynamic mechanistic model, allowing us to fit our model to available viability data, and then use the estimated parameters to perform detailed PK-PD simulations in an in vivo context.

Fitting the model to viability data. Figure 2 shows the fitting result for the early ring stage. Results for the other stages are provided in Fig. S1 to S3 in the supplemental material. Parameter estimates and confidence intervals (CI) are given in Table 1. The model captures the data very well, in particular the dependence of viability on drug exposure time, which is the key advance we require. Note, as discussed in the introduction, the standard model (in which there is no stress accumulation) is incapable of explaining the experimental observations for short exposure times. The standard model always predicts a faster decrease in viability for a higher drug concentration, contradicting the data (Fig. 1). This is further confirmed by a statistical comparison of the dynamic stress model and the standard model both in terms of fit visualization (Fig. S4 to S7) and the Akaike information criterion (AIC) (Tables S2 and S3). In particular, the improvement of goodness of fit is most evident for the mid-ring stage (Fig. S5). This finding comes as no surprise; the estimates for λ for the mid-ring stage are small relative to the λ→∞ value, which recovers the standard (no stress accumulation) model. This analysis demonstrates that the introduction of a time-dependent stress response, modulating both the maximal killing rate and the half-maximal drug concentration, provides a substantially improved explanation for the in vitro experimental data.

We also compare the model fits on a logarithmic scale (Fig. S8 to S11). The dynamic stress model outperforms the standard model overall, particularly so for low concentrations and short exposure times, which are the primary point of interest. However, for some cases of high DHA concentration and long exposure times (e.g., for 4-h and 6-h exposures at drug concentrations larger than 100 nM), we identify that the dynamic stress model subtly overestimates the killing effect for early ring and early and late trophozoite stages (Fig. S8, S10, and S11). Of course, the log scale emphasizes the subtle effect while obscuring the significantly improved fit for shorter exposure durations and lower drug concentrations (where the interesting and new dynamic effects are evident). Furthermore, the limit of detection in the in vitro experiment remains uncertain, and while we conservatively assumed a detection limit of 0.005 for data fitting, it may be as high as 0.05. Refitting with a higher detection limit (of 0.05) provides nearly indistinguishable best fits to the data with no meaningful change in parameter estimates (data not shown) but does provide even stronger statistical support for the dynamic stress model over the standard model given the relatively diminished contribution to the likelihood from low-viability data (more of which now falls below the detection threshold).

Drug concentration-killing rate curves and stage dependency.The overall impact of parasite killing is determined primarily by the drug concentration-killing rate curve, which we now consider to be a function of exposure time, generalizing the usual modeling assumption that the killing rate is an instantaneous function of drug concentration. Figure 3 shows the modeled evolution of the concentration-killing rate curve for the four different life stages of the parasites used in Fig. 1. Except for the early ring stage, for which the curve reaches its steady state very quickly, the delayed process of approaching the steady-state killing rate curve for the other three stages is biologically significant, particularly for the mid-ring stage where a very strong delay is observed.

• Open in new tab
• Download powerpoint

FIG 3

Model results showing the evolution of the drug concentration-killing rate curve with drug exposure duration for different stages. The time after drug exposure, t, is indicated. Note that the y axis scale differs for different stages.

Figure 4 shows the estimates and 95% confidence intervals for the four model parameters λ, α, β₁, and β₂ by stage. The delay in the evolution of the drug concentration-killing rate curve is determined primarily by the parameter λ (Fig. 4A and Table 1), which reflects the accumulation rate for stress (equation 6). For the early ring stage, λ = 6.25 h⁻¹ and thus S(1 h) > 0.99, indicating that early rings rapidly succumb to drug exposure. In contrast, the rate of accumulation of stress for mid-ring-stage parasites is much lower (λ = 0.37 h⁻¹), and it would take over 12 h of continued exposure to drug for S to exceed 0.99. Early and late trophozoite stages display similar characteristics in terms of the rate of accumulation of stress (Fig. 3). The temporal effect can also be interpreted by introducing the half-life (Table 1) of the unstressed state, given by ln(2)/λ. Accordingly, λ can be seen as the rate at which the unstressed state is lost. We see that the stress in mid-ring-stage parasites takes approximately 2 h to reach 50%, while stress accumulates for other stages more rapidly, taking less than half an hour to reach the same level.

• Open in new tab
• Download powerpoint

FIG 4

Dependence of model parameters on parasite life stage. The parameter estimates are provided in Table 1. Error bars show the model-based 95% CI for each parameter. (A) The rate of stress development (λ) is large for early rings, indicating minimal time dependence in killing for this life stage. In contrast for mid-ring-stage parasites, the rate is very small, indicating a substantial accumulation effect. (B) The maximum killing rate is linear in α, indicating that the maximum killing rate is higher in trophozoites than rings. (C and D) Parameters (β₁ and β₂) model the relationship between stress (S) and the half-maximal drug concentration (equation 8).

The maximum killing rate, α (Fig. 4B), shows a significant reduction for (early and mid-) ring-stage parasites than for (early and late) trophozoites, suggesting that young parasites are more resistant to DHA than mature parasites once (or even when) the killing effect has reached a steady state. The parameters related to the concentration required to achieve the half-maximal killing rate, β₁ and β₂, also exhibit stage specificity (Fig. 4C and D). In particular, the stationary half-maximal killing concentration (Fig. 4D, β₂) shows that ring-stage parasites exhibit a higher sensitivity to drug at steady state than trophozoites. However, it must be remembered that, particularly for mid-ring-stage parasites, the progress toward that steady state (governed by λ) is slow, and the net effect of the dynamics of drug-induced killing is best understood through Fig. 3.

Incorporating the delayed drug effect into PK-PD modeling.Having established the applicability of our model to the in vitro data, we now consider the potential implications for the in vivo application of artemisinin-based medication. We do so by incorporating the time dependency on the killing rate into the general PK-PD framework (equation 3), where we allow for realistic drug PK and a general age structure for the parasites.

Because the stress accumulation effect is most pronounced for the mid-ring stage of 3D7 parasites (the subject of our study), we begin by considering mid-ring-stage parasites treated with a single dose of artesunate (2 mg/kg of body weight). The plasma DHA concentration, C(t), displays biphasic behavior (6): dCdt={Cmaxtm,0≤t<tm−ln(2)t1/2C,t≥tm(18) where C_max is the maximum achievable concentration and t_m indicates the time at which that maximum concentration is achieved. Note that the half-life, t_1/2, refers to the in vivo half-life of DHA, which is much smaller than that measured in vitro due to altered physiological conditions (Fig. 2 and Table S2 in reference 21), i.e., C_max of 2,820 nM, t_m of 1 h, and t_1/2 of 0.9 h per references 2 and 21. The simulated PK data are shown in Fig. 5A, upper.

• Open in new tab
• Download powerpoint

FIG 5

Incorporation of the time-dependent killing rate into the PK-PD model. We study the mid-ring stage for illustrative purposes. Parameter values are taken from Table 1. (A) Simulated in vivo DHA concentration profile (upper; the in vivo half-life is approximately 0.9 h [2, 21]), the kinetics for the modulatory stress variable, S (middle; black curve, λ = 0.37), and the transient killing rate S (lower; black curve, λ = 0.37) induced by the drug pulse. The middle and lower panels also show how S and k evolve if λ is higher (blue) or lower (red). (B) Killing rate surface as a function of DHA concentration, C, and the stress, S, and the projection of the trajectory of the effective killing rate (i.e., a projection of the curves in panel A [lower]) onto the surface. (C) Area under the killing rate curve, an indication of the total amount of killing achievable over the course of the drug pulse.

The middle and lower panels of Fig. 5A shows the time series for stress, S, and killing rate, k, as a result of the changing DHA concentration. The black curves are generated using λ = 0.37 h⁻¹, the best-fit estimate from fitting the model to the in vitro data for the mid-ring stage (Table 1). Decreasing λ will delay the increase of S and in turn lead to a slower and shortened killing rate profile (Fig. 5A, red curves), while increasing λ will do the opposite (Fig. 5A, blue curves). Thus, we consider a lowering of λ as a potential manifestation of (stage-specific) artemisinin resistance. Conversely, λ = 1 implies that parasites accumulate stress rapidly and are rendered nonviable at the stationary rate [given C(t)] soon after initiation of the drug pulse.

Under the model, the killing rate, k, is now a function of two variables (as clearly shown by equation 17), the drug concentration, C(t), and the modulatory stress, S(t). We can represent this graphically by displaying the killing rate as a trajectory on a surface in (C,S) space (Fig. 5B). In this representation, we can clearly see the effect of λ on the evolution of the killing rate: the trajectory corresponding to a smaller λ has less time (controlled by S in the model) to climb up the killing rate surface, even when the achievable DHA concentration remains the same.

We can also consider the net cumulative effect of the drug pulse. As described earlier, ϕ=∫0tdk(C(τ),S(τ))dτ represents the cumulative stress effect, which can be used to indicate drug efficacy. In the in vivo simulation, ϕ is simply the area under the effective killing rate curve (i.e., the area under the curve in the lower panel of Fig. 5A). For mid-ring-stage parasites with λ = 0.37, the cumulative stress effect (Fig. 5C, black bar) corresponds to a reduction in viability of approximately 99.75% over the drug pulse. Further numerical exploration indicates that a roughly 3-fold increase or decrease in λ leads to a significant difference in the cumulative stress effect (Fig. 5C, blue and red bars) and, in turn, a difference of a few orders of magnitude in viability.

In summary, the results presented in Fig. 5 indicate that the temporal drug effect significantly affects the in vivo parasite killing and thus should be considered in model-based prediction of clinical treatment. Furthermore, the visualization of the killing rate trajectory on the (C,S) plane surface suggests a clear evolutionary strategy for the parasite to escape drug pressure, particularly given the short elimination half-life of artemisinin and its derivatives. An ability to outlast the short drug pulse provides an effective means of escape quite distinct from any changes in susceptibility, as are typically considered by a change in the maximal killing rate or drug concentration required to achieve half-maximal killing.

To fully explore the consequences of accumulation effects on the pharmacodynamics of antimalarial treatment, we simulated the time course of total viable parasite count under a standard AS7 dosing regimen (i.e., a dose of 2 mg/kg artesunate every 24 h for 7 days) for both the 3D7 strain and a hypothetical strain which exhibits a slower rate of stress development during the mid-ring stage. Some key simulation details are provided in the legend to Fig. 6. We initiated the simulation with 10¹² parasites per patient with a normally distributed age distribution, with a mean of 10 h p.i. and standard deviation of 2 h p.i. (Fig. 6, inset). For the laboratory 3D7 strain, the model predicts that effective parasite clearance is achieved immediately following the third dose of artesunate (at 48 h in the model) (Fig. 6, green curve). In contrast, for the hypothetical strain which exhibits a slower development of the stress response during the mid-ring stage (i.e., λ is reduced for this, but no other, stage), we observe a clear and substantial delay in parasite clearance. In detail, the red curve in Fig. 6 shows the parasitemia curve for a resistant strain that has a λ value of 0.1 h⁻¹ for the mid-ring stage (with all other parameters [across all stages] unchanged). This simulation has an a priori rationale given previous studies that indicate that field isolates from Pailin (western Cambodia) display a reduced sensitivity to artemisinin-based therapies during the ring stage of infection (6, 8, 9). We note that while this simple simulation does not incorporate the process of splenic clearance or the immune response (which are also important for in vivo parasite clearance; see Discussion), its behavior is consistent with clinical observations of a 1.5 to 2 times longer time to clearance for resistant strains compared to sensitive strains (6).

• Open in new tab
• Download powerpoint

FIG 6

Simulation of parasite killing under a standard treatment of 2 mg/kg artesunate every 24 h. The inset shows the age distribution of a total of 10¹² parasites (per patient) at the start of treatment [∼N(10, 2²)]. The parasite multiplication factor is assumed to be 10 (6, 13), which means that 10 new parasites are produced once a parasite reaches 48 h p.i. (i.e., r = 10 in the model). The PK profile is a series of repeated DHA concentration profiles every 24 h (i.e., repeated simulations of DHA concentration profile shown in Fig. 5A, upper). The black triangles indicate when the doses are given. The green curve corresponding to the laboratory 3D7 strain is generated using the parameters in Table 1, while the red curve is generated using the same set of parameters except for reducing λ for mid-ring stage to be 0.1 h⁻¹ to simulate a more resistant strain. With limited information, we simply divide the 48-h life cycle into early ring stage (0 to 6 h p.i.), mid-ring stage (6 to 26 h p.i.), early trophozoite (26 to 34 h p.i.), and late trophozoite (34 to 48 h p.i.) in the simulation. The modulatory variable S is assumed to follow equation 6 only when DHA concentration, C, is ≥0.1 nM (i.e., C* = 0.1 nM), and S is immediately reset to zero when the DHA concentration drops below 0.1 nM.

In addition, we also compare the in vivo simulation result generated by using the standard killing rate model to the results shown in Fig. 6 (i.e., the results generated by using the dynamic stress model) in Fig. S12. The two models produce qualitatively similar results. Parasites are killed in a stepwise manner associated with the time of drug administration. However, the stress model predicts a much faster parasite clearance than the standard model. This is due to the different ways in which the dynamic stress and traditional models explain the in vitro data. Without an accumulation effect, the standard model does its best to explain the in vitro data by having a lower maximal killing rate (of course, this leads to poor fits compared to the dynamic stress model, so little confidence should be placed in the extrapolated in vivo parasite-time curve). In turn, the lower estimates for k_max compared to α (and higher estimates for K_c compared to β₂) lead to slower parasite killing in the in vivo context. We note that, like for all PK-PD studies that attempt to translate from the in vitro to in vivo context, care must be taken as the in vivo environment may result in quantitative changes to the values of key model parameters (13). Accordingly, the key finding from our in vivo simulation is a qualitative, not quantitative, one. We have demonstrated how a change in the rate of accumulation of stress may influence the clearance time and how ignoring stress may provide vastly different predictions for the clearance time. Conducting such a quantitative translation exercise to the in vivo context using the techniques of Zaloumis et al. (13) will be the subject of future work.