A mechanistic model quantifies artemisinin-induced parasite growth   retardation in blood-stage Plasmodium falciparum infection

Pengxing Cao; Nectarios Klonis; Sophie Zaloumis; David S. Khoury,; Deborah Cromer; Miles P. Davenport; Leann Tilley; Julie A. Simpson; James M.; McCaw

arXiv:1701.05302·q-bio.PE·October 24, 2017

A mechanistic model quantifies artemisinin-induced parasite growth retardation in blood-stage Plasmodium falciparum infection

Pengxing Cao, Nectarios Klonis, Sophie Zaloumis, David S. Khoury,, Deborah Cromer, Miles P. Davenport, Leann Tilley, Julie A. Simpson, James M., McCaw

PDF

TL;DR

This paper presents a mechanistic model that quantifies how artemisinin delays parasite growth in blood-stage Plasmodium falciparum, influencing parasite clearance and informing optimal dosing strategies.

Contribution

The study introduces a novel mechanistic model linking artemisinin concentration to parasite growth retardation, enhancing understanding of drug effects on malaria parasites.

Findings

01

Model accurately reproduces in vitro fluorescence data.

02

Growth retardation depends on drug concentration.

03

Growth delay significantly impacts in vivo parasite dynamics.

Abstract

Falciparum malaria is a major parasitic disease causing widespread morbidity and mortality globally. Artemisinin derivatives---the most effective and widely-used antimalarials that have helped reduce the burden of malaria by 60% in some areas over the past decade---have recently been found to induce growth retardation of blood-stage Plasmodium falciparum when applied at clinically relevant concentrations. To date, no model has been designed to quantify the growth retardation effect and to predict the influence of this property on in vivo parasite killing. Here we introduce a mechanistic model of parasite growth from the ring to trophozoite stage of the parasite's life cycle, and by modelling the level of staining with an RNA-binding dye, we demonstrate that the model is able to reproduce fluorescence distribution data from in vitro experiments using the laboratory 3D7 strain. We…

Tables8

Table 1. Table 1: A list of model parameters The unit h p.i. is an abbreviation of hour post-infection. We assume in the model that the SYTO-61 fluorescence intensity is unitless.

Model parameter	Description	Unit
$μ$	mean of parasite age distribution at the time of data collection	h p.i.
$σ$	standard deviation of age distribution	h p.i.
$A_{r}$	the age before which the ring-to-trophozoite transition cannot occur	h p.i.
$λ$	rate of conversion from ring to trophozoite	$h^{- 1}$
$F_{0}$	a parameter in the age-to-signal mapping (Eq. 2)	–
$r_{1}$	the rate at which the age-dependent intensity of SYTO-61 staining increases in ring stage	${[h p . i .]}^{- 1}$
$r_{2}$	the rate at which the age-dependent intensity of SYTO-61 staining increases in trophozoite stage	${[h p . i .]}^{- 1}$

Table 2. Table 2: A list of new parameters in the probability density function of the reduced fluorescence intensity distribution. The new parameters are some combinations of the original model parameters listed in Table 1 and are primarily used for data fitting. The choice of parameter constraints are provided in the main text.

New parameter	Relationship to the model parameters in Table 1	Constraint
$K_{1}$	$\ln (F_{0}) + r_{1} μ$	[6.2, 10.8]
$K_{2}$	$σ r_{1}$	[0.0039, 1]
$K_{3}$	$\ln (F_{0}) + r_{1} A_{r}$	[2, 10.8]
$K_{4}$	$λ / r_{1}$	[0, 300]
$K_{5}$	$r_{2} / r_{1}$	(1, 10]

Table 3. Table 3: *

ART conc. (nM)	Parameter estimate (95% CI)
ART conc. (nM)	$K_{1}$	$K_{2}$	Other parameters
0	9.1501 (9.1232, 9.1771)	0.3287 (0.3068, 0.3506)	$K_{3} =$ 8.6861 (8.6715, 8.7006); $K_{4} =$ 16.6432 (11.5799, 21.7065); $K_{5} =$ 1.6048 (1.5330, 1.6767).
0	9.2022 (9.1734, 9.2310)	0.3192 (0.2978, 0.3406)
1.2	9.1868 (9.1591, 9.2145)	0.3102 (0.2896, 0.3308)
2.4	9.1843 (9.1560, 9.2126)	0.3293 (0.3072, 0.3514)
4.9	9.1825 (9.1549, 9.2100)	0.3168 (0.2957, 0.3378)
9.8	9.1510 (9.1235, 9.1785)	0.3422 (0.3197, 0.3647)
19.5	9.1355 (9.1085, 9.1625)	0.3469 (0.3237, 0.3702)
39.1	9.0883 (9.0615, 9.1150)	0.3598 (0.3355, 0.3840)
78.1	9.0377 (9.0118, 9.0637)	0.3703 (0.3451, 0.3954)
156.3	9.0183 (8.9930, 9.0435)	0.3613 (0.3372, 0.3855)
312.5	8.9805 (8.9576, 9.0034)	0.3572 (0.3346, 0.3798)
625	8.9693 (8.9483, 8.9903)	0.3445 (0.3232, 0.3659)
1250	8.9400 (8.9178, 8.9622)	0.3561 (0.3338, 0.3785)
2500	8.9094 (8.8891, 8.9297)	0.3500 (0.3291, 0.3708)
5000	8.8730 (8.8529, 8.8930)	0.3557 (0.3352, 0.3763)
10000	8.8582 (8.8398, 8.8766)	0.3340 (0.3151, 0.3529)

Table 4. Table 4: *

DHA conc. (nM)	Parameter estimate (95% CI)
DHA conc. (nM)	$K_{1}$	$K_{2}$	Other parameters
0	9.1011 (9.0727, 9.1294)	0.2816 (0.2599, 0.3033)	$K_{3} =$ 8.6243 (8.6084, 8.6402); $K_{4} =$ 12.3188 (9.1899, 15.4477); $K_{5} =$ 1.7590 (1.6677, 1.8502).
0	9.1343 (9.1044, 9.1642)	0.2884 (0.2660, 0.3107)
0.12	9.1425 (9.1122, 9.1729)	0.2908 (0.2682, 0.3134)
0.24	9.1340 (9.1049, 9.1631)	0.2886 (0.2670, 0.3103)
0.5	9.1305 (9.1006, 9.1605)	0.2872 (0.2648, 0.3096)
1	9.1151 (9.0855, 9.1446)	0.2960 (0.2730, 0.3190)
2	9.1054 (9.0762, 9.1346)	0.2981 (0.2750, 0.3212)
3.9	9.0752 (9.0463, 9.1041)	0.3069 (0.2830, 0.3308)
7.8	9.0556 (9.0273, 9.0839)	0.3075 (0.2832, 0.3319)
15.6	9.0204 (8.9925, 9.0483)	0.3227 (0.2977, 0.3477)
31.3	8.9084 (8.8837, 8.9331)	0.3408 (0.3159, 0.3656)
62.5	8.8085 (8.7886, 8.8283)	0.3136 (0.2934, 0.3339)
125	8.6865 (8.6649, 8.7080)	0.3328 (0.3131, 0.3525)
250	8.6622 (8.6427, 8.6817)	0.3248 (0.3071, 0.3425)
500	8.6079 (8.5892, 8.6266)	0.3063 (0.2908, 0.3218)
1000	8.7061 (8.6912, 8.7209)	0.2640 (0.2491, 0.2789)

Table 5. Table 5: *

Parameter	Estimate	95% CI	True value	Corresponding true model parameter
$K_{1}$	8.7572	(8.7499, 8.7644)	8.7557	$F_{0} = 20$ , $r_{1} = 0.24$ , $μ = 24$
$K_{2}$	0.1494	(0.1423, 0.1564)	0.1440	$σ = 0.6$ , $r_{1} = 0.24$
$K_{3}$	7.9702	(7.8422, 8.0982)	8.0357	$F_{0} = 20$ , $r_{1} = 0.24$ , $A_{r} = 21$
$K_{3}$	7.9939	(7.8759, 8.1119)	8.0357	$F_{0} = 20$ , $r_{1} = 0.24$ , $A_{r} = 21$
$K_{3}$	7.9848	(7.8660, 8.1035)	8.0357	$F_{0} = 20$ , $r_{1} = 0.24$ , $A_{r} = 21$
$K_{3}$	8.0013	(7.8823, 8.1202)	8.0357	$F_{0} = 20$ , $r_{1} = 0.24$ , $A_{r} = 21$
$K_{4}$	1.0746	(0.8548, 1.2944)	1.2500	$r_{1} = 0.24$ , $λ = 0.3$
$K_{4}$	2.3377	(1.8877, 2.7877)	2.5000	$r_{1} = 0.24$ , $λ = 0.6$
$K_{4}$	3.3622	(2.6702, 4.0542)	3.7500	$r_{1} = 0.24$ , $λ = 0.9$
$K_{4}$	4.8077	(3.6120, 6.0031)	5.0000	$r_{1} = 0.24$ , $λ = 1.2$
$K_{5}$	2.2561	(2.0861, 2.4260)	2.3333	$r_{1} = 0.24$ , $r_{2} = 0.56$

Table 6. Table 6: *

Parameter	Estimate	95% CI	True value	Corresponding true model parameter
$K_{1}$	8.7572	(8.7538, 8.7606)	8.7557	$F_{0} = 20$ , $r_{1} = 0.24$ , $μ = 24$
$K_{2}$	0.1452	(0.1416, 0.1489)	0.1440	$σ = 0.6$ , $r_{1} = 0.24$
$K_{3}$	6.6771	(6.5006, 6.8537)	6.5957	$F_{0} = 20$ , $r_{1} = 0.24$ , $A_{r} = 15$
$K_{3}$	7.3726	(7.2536, 7.4915)	7.3157	$F_{0} = 20$ , $r_{1} = 0.24$ , $A_{r} = 18$
$K_{3}$	8.0640	(7.9957, 8.1323)	8.0357	$F_{0} = 20$ , $r_{1} = 0.24$ , $A_{r} = 21$
$K_{3}$	8.7345	(7.9957, 8.7901)	8.7557	$F_{0} = 20$ , $r_{1} = 0.24$ , $A_{r} = 24$
$K_{4}$	2.3022	(2.0296, 2.5748)	2.0833	$r_{1} = 0.24$ , $λ = 0.5$
$K_{4}$	2.1793	(1.9382, 2.4203)	2.0833	$r_{1} = 0.24$ , $λ = 0.5$
$K_{4}$	2.1407	(1.8693, 2.4122)	2.0833	$r_{1} = 0.24$ , $λ = 0.5$
$K_{4}$	1.8049	(0.8657, 2.7440)	2.0833	$r_{1} = 0.24$ , $λ = 0.5$
$K_{5}$	2.3748	(2.2617, 2.4878)	2.3333	$r_{1} = 0.24$ , $r_{2} = 0.56$

Table 7. Table 7: *

Parameter	Estimate	95% CI	True value	Corresponding true model parameter
$K_{1}$	8.7549	(8.7509, 8.7590)	8.7557	$F_{0} = 20$ , $r_{1} = 0.24$ , $μ = 24$
$K_{2}$	0.1474	(0.1442, 0.1506)	0.1440	$σ = 0.6$ , $r_{1} = 0.24$
$K_{3}$	7.2199	(7.0792, 7.3606)	7.3157	$F_{0} = 20$ , $r_{1} = 0.24$ , $A_{r} = 18$
$K_{3}$	7.7435	(7.6476, 7.8395)	7.7957	$F_{0} = 20$ , $r_{1} = 0.24$ , $A_{r} = 20$
$K_{3}$	8.2598	(8.1946, 8.3250)	8.2757	$F_{0} = 20$ , $r_{1} = 0.24$ , $A_{r} = 22$
$K_{3}$	8.8085	(8.7722, 8.8448)	8.7557	$F_{0} = 20$ , $r_{1} = 0.24$ , $A_{r} = 24$
$K_{4}$	4.2809	(3.5036, 5.0583)	5.0000	$r_{1} = 0.24$ , $λ = 1.2$
$K_{4}$	3.9442	(3.2968, 4.5917)	4.1667	$r_{1} = 0.24$ , $λ = 1$
$K_{4}$	3.0947	(2.4884, 3.7009)	3.3333	$r_{1} = 0.24$ , $λ = 0.8$
$K_{4}$	5.0574	(2.2743, 7.8405)	2.5000	$r_{1} = 0.24$ , $λ = 0.6$
$K_{5}$	2.2694	(2.1635, 2.3754)	2.3333	$r_{1} = 0.24$ , $r_{2} = 0.56$

Table 8. Table 8: *

Parameter	Estimate	95% CI	True value	Corresponding true model parameter
$K_{1}$	8.7520	(8.7482, 8.7558)	8.7557	$F_{0} = 20$ , $r_{1} = 0.24$ , $μ = 24$
$K_{2}$	0.1458	(0.1425, 0.1491)	0.1440	$σ = 0.6$ , $r_{1} = 0.24$
$K_{3}$	8.7630	(8.7399, 8.7862)	8.7557	$F_{0} = 20$ , $r_{1} = 0.24$ , $A_{r} = 24$
$K_{3}$	8.2444	(8.1862, 8.3026)	8.2757	$F_{0} = 20$ , $r_{1} = 0.24$ , $A_{r} = 22$
$K_{3}$	7.7834	(7.6986, 7.8681)	7.7957	$F_{0} = 20$ , $r_{1} = 0.24$ , $A_{r} = 20$
$K_{3}$	7.2919	(7.1706, 7.4133)	7.3157	$F_{0} = 20$ , $r_{1} = 0.24$ , $A_{r} = 18$
$K_{4}$	5.6203	(3.7485, 7.4921)	5.0000	$r_{1} = 0.24$ , $λ = 1.2$
$K_{4}$	3.8541	(3.1739, 4.5343)	4.1667	$r_{1} = 0.24$ , $λ = 1$
$K_{4}$	3.3537	(2.9171, 3.7903)	3.3333	$r_{1} = 0.24$ , $λ = 0.8$
$K_{4}$	2.4964	(2.2040, 2.7889)	2.5000	$r_{1} = 0.24$ , $λ = 0.6$
$K_{5}$	2.3139	(2.2106, 2.4172)	2.3333	$r_{1} = 0.24$ , $r_{2} = 0.56$

Equations53

f_{c} = f - (1 - V) f_{b},

f_{c} = f - (1 - V) f_{b},

F(a)=\left\{\begin{array}[]{@{}c@{\quad}l@{}}F_{0}e^{r_{1}a},&a<A_{c}\\ F_{0}e^{r_{1}A_{c}+r_{2}(a-A_{c})},&a\geq A_{c}\\ \end{array}\right.

F(a)=\left\{\begin{array}[]{@{}c@{\quad}l@{}}F_{0}e^{r_{1}a},&a<A_{c}\\ F_{0}e^{r_{1}A_{c}+r_{2}(a-A_{c})},&a\geq A_{c}\\ \end{array}\right.

P r (F) = \frac{1}{σ 2 π} e^{- \frac{( a _{1} - μ ) ^{2}}{2 σ ^{2}}} (\frac{1}{r _{1} F}),

P r (F) = \frac{1}{σ 2 π} e^{- \frac{( a _{1} - μ ) ^{2}}{2 σ ^{2}}} (\frac{1}{r _{1} F}),

P r [F, from rings with age a = a_{2}^{'}] = \frac{1}{σ 2 π} e^{- \frac{( a _{2}^{'} - μ ) ^{2}}{2 σ ^{2}}} [e^{- λ (a_{2}^{'} - A_{r})}] (\frac{1}{r _{1} F}) .

P r [F, from rings with age a = a_{2}^{'}] = \frac{1}{σ 2 π} e^{- \frac{( a _{2}^{'} - μ ) ^{2}}{2 σ ^{2}}} [e^{- λ (a_{2}^{'} - A_{r})}] (\frac{1}{r _{1} F}) .

P r [F, from trophozoites] = \int_{a_{2}}^{a_{2}^{'}} P r (F ∣ a) d a,

P r [F, from trophozoites] = \int_{a_{2}}^{a_{2}^{'}} P r (F ∣ a) d a,

\int_{F}^{F + Δ F} P r (F ∣ a) d F = \frac{1}{σ 2 π} e^{- \frac{( a - μ ) ^{2}}{2 σ ^{2}}} [e^{- λ (A_{c} (F + Δ F) - A_{r})} - e^{- λ (A_{c} (F) - A_{r})}],

\int_{F}^{F + Δ F} P r (F ∣ a) d F = \frac{1}{σ 2 π} e^{- \frac{( a - μ ) ^{2}}{2 σ ^{2}}} [e^{- λ (A_{c} (F + Δ F) - A_{r})} - e^{- λ (A_{c} (F) - A_{r})}],

P r (F ∣ a) = \frac{1}{σ 2 π} e^{- \frac{( a - μ ) ^{2}}{2 σ ^{2}}} [λ e^{- λ (A_{c} - A_{r})}] \frac{1}{( r _{2} - r _{1} ) F},

P r (F ∣ a) = \frac{1}{σ 2 π} e^{- \frac{( a - μ ) ^{2}}{2 σ ^{2}}} [λ e^{- λ (A_{c} - A_{r})}] \frac{1}{( r _{2} - r _{1} ) F},

Pr(F)=\left\{\begin{array}[]{@{}c@{\quad}l@{}}\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(a_{1}-\mu)^{2}}{2\sigma^{2}}}\left(\frac{1}{r_{1}F}\right),&F\leq F_{r}\\ \frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(a_{2}^{\prime}-\mu)^{2}}{2\sigma^{2}}}\left[e^{-\lambda(a_{2}^{\prime}-A_{r})}\right]\left(\frac{1}{r_{1}F}\right)+\int_{a_{2}}^{a_{2}^{\prime}}Pr(F|a)da,&F>F_{r}\\ \end{array}\right.

Pr(F)=\left\{\begin{array}[]{@{}c@{\quad}l@{}}\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(a_{1}-\mu)^{2}}{2\sigma^{2}}}\left(\frac{1}{r_{1}F}\right),&F\leq F_{r}\\ \frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(a_{2}^{\prime}-\mu)^{2}}{2\sigma^{2}}}\left[e^{-\lambda(a_{2}^{\prime}-A_{r})}\right]\left(\frac{1}{r_{1}F}\right)+\int_{a_{2}}^{a_{2}^{\prime}}Pr(F|a)da,&F>F_{r}\\ \end{array}\right.

K_{1} = ln (F_{0}) + r_{1} μ,

K_{1} = ln (F_{0}) + r_{1} μ,

K_{2} = σ r_{1} .

K_{2} = σ r_{1} .

P r (F) = \frac{1}{( σ r _{1} ) 2 π} e^{- \frac{[ ln ( F ) - ( ln ( F _{0} ) + r _{1} μ ) ] ^{2}}{2 ( σ r _{1} ) ^{2}}} (\frac{1}{F}) = \frac{1}{K _{2} 2 π} e^{- \frac{[ ln ( F ) - K _{1} ] ^{2}}{2 K _{2}^{2}}} (\frac{1}{F}),

P r (F) = \frac{1}{( σ r _{1} ) 2 π} e^{- \frac{[ ln ( F ) - ( ln ( F _{0} ) + r _{1} μ ) ] ^{2}}{2 ( σ r _{1} ) ^{2}}} (\frac{1}{F}) = \frac{1}{K _{2} 2 π} e^{- \frac{[ ln ( F ) - K _{1} ] ^{2}}{2 K _{2}^{2}}} (\frac{1}{F}),

K_{3} = ln (F_{0}) + r_{1} A_{r},

K_{3} = ln (F_{0}) + r_{1} A_{r},

K_{4} = \frac{λ}{r _{1}},

K_{4} = \frac{λ}{r _{1}},

\frac{1}{σ 2 π} e^{- \frac{( a _{2}^{'} - μ ) ^{2}}{2 σ ^{2}}} [e^{- λ (a_{2}^{'} - A_{r})}] (\frac{1}{r _{1} F}) = \frac{1}{K _{2} 2 π} e^{- \frac{[ ln ( F ) - K _{1} ] ^{2}}{2 K _{2}^{2}}} [e^{- K_{4} (ln (F) - K_{3})}] (\frac{1}{F}) .

\frac{1}{σ 2 π} e^{- \frac{( a _{2}^{'} - μ ) ^{2}}{2 σ ^{2}}} [e^{- λ (a_{2}^{'} - A_{r})}] (\frac{1}{r _{1} F}) = \frac{1}{K _{2} 2 π} e^{- \frac{[ ln ( F ) - K _{1} ] ^{2}}{2 K _{2}^{2}}} [e^{- K_{4} (ln (F) - K_{3})}] (\frac{1}{F}) .

K_{5} = \frac{r _{2}}{r _{1}} .

K_{5} = \frac{r _{2}}{r _{1}} .

P r (F ∣ a) = \frac{1}{K _{2} 2 π} e^{- \frac{( r _{1} a + ln ( F _{0} ) - K _{1} ) ^{2}}{2 K _{2}^{2}}} (e^{- Φ}) \frac{K _{4} r _{1}}{( K _{5} - 1 ) F},

P r (F ∣ a) = \frac{1}{K _{2} 2 π} e^{- \frac{( r _{1} a + ln ( F _{0} ) - K _{1} ) ^{2}}{2 K _{2}^{2}}} (e^{- Φ}) \frac{K _{4} r _{1}}{( K _{5} - 1 ) F},

Φ = \frac{K _{4} [ K _{5} ( r _{1} a + ln ( F _{0} )) - ln ( F ) + K _{3} - K _{3} K _{5} ]}{K _{5} - 1} .

Φ = \frac{K _{4} [ K _{5} ( r _{1} a + ln ( F _{0} )) - ln ( F ) + K _{3} - K _{3} K _{5} ]}{K _{5} - 1} .

\int_{a_{2}}^{a_{2}^{'}} P r (F ∣ a) d a = \int_{b_{2}}^{b_{2}^{'}} P r (F ∣ b) d b = \int_{b_{2}}^{b_{2}^{'}} \frac{1}{K _{2} 2 π} e^{- \frac{[ b - K _{1} ] ^{2}}{2 K _{2}^{2}}} (e^{- Φ}) \frac{K _{4}}{( K _{5} - 1 ) F} d b,

\int_{a_{2}}^{a_{2}^{'}} P r (F ∣ a) d a = \int_{b_{2}}^{b_{2}^{'}} P r (F ∣ b) d b = \int_{b_{2}}^{b_{2}^{'}} \frac{1}{K _{2} 2 π} e^{- \frac{[ b - K _{1} ] ^{2}}{2 K _{2}^{2}}} (e^{- Φ}) \frac{K _{4}}{( K _{5} - 1 ) F} d b,

Φ = \frac{K _{4} [ K _{5} b - ln ( F ) + K _{3} - K _{3} K _{5} ]}{K _{5} - 1},

Φ = \frac{K _{4} [ K _{5} b - ln ( F ) + K _{3} - K _{3} K _{5} ]}{K _{5} - 1},

b_{2} = \frac{ln ( F ) - K _{3} + K _{3} K _{5}}{K _{5}} and b_{2}^{'} = ln (F) .

b_{2} = \frac{ln ( F ) - K _{3} + K _{3} K _{5}}{K _{5}} and b_{2}^{'} = ln (F) .

Pr(F)=\left\{\begin{array}[]{@{}c@{\quad}l@{}}\frac{1}{K_{2}\sqrt{2\pi}}e^{-\frac{[{\rm ln}(F)-K_{1}]^{2}}{2K_{2}^{2}}}\left(\frac{1}{F}\right),&F\leq e^{K_{3}}\\ \frac{1}{K_{2}\sqrt{2\pi}}e^{-\frac{({\rm ln}(F)-K_{1})^{2}}{2K_{2}^{2}}}\left[e^{-K_{4}({\rm ln}(F)-K_{3})}\right]\left(\frac{1}{F}\right)+\int_{b_{2}}^{b_{2}^{\prime}}Pr(F|b)db,&F>e^{K_{3}}\\ \end{array}\right.

Pr(F)=\left\{\begin{array}[]{@{}c@{\quad}l@{}}\frac{1}{K_{2}\sqrt{2\pi}}e^{-\frac{[{\rm ln}(F)-K_{1}]^{2}}{2K_{2}^{2}}}\left(\frac{1}{F}\right),&F\leq e^{K_{3}}\\ \frac{1}{K_{2}\sqrt{2\pi}}e^{-\frac{({\rm ln}(F)-K_{1})^{2}}{2K_{2}^{2}}}\left[e^{-K_{4}({\rm ln}(F)-K_{3})}\right]\left(\frac{1}{F}\right)+\int_{b_{2}}^{b_{2}^{\prime}}Pr(F|b)db,&F>e^{K_{3}}\\ \end{array}\right.

\frac{d R}{d t}

\frac{d R}{d t}

\frac{d T}{d t}

\frac{d a}{d t}

\lambda=\left\{\begin{array}[]{@{}c@{\quad}l@{}}0,&a<A_{r}\\ \frac{22}{34+C},&a\geq A_{r}\\ \end{array}\right.

\lambda=\left\{\begin{array}[]{@{}c@{\quad}l@{}}0,&a<A_{r}\\ \frac{22}{34+C},&a\geq A_{r}\\ \end{array}\right.

k_{r}

k_{r}

k_{t}

\frac{dC}{dt}=\left\{\begin{array}[]{@{}c@{\quad}l@{}}\frac{C_{max}}{t_{m}},&t<t_{m}\\ -\frac{ln(2)}{t_{1/2}}C,&t\geq t_{m}\\ \end{array}\right.

\frac{dC}{dt}=\left\{\begin{array}[]{@{}c@{\quad}l@{}}\frac{C_{max}}{t_{m}},&t<t_{m}\\ -\frac{ln(2)}{t_{1/2}}C,&t\geq t_{m}\\ \end{array}\right.

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

A mechanistic model quantifies artemisinin-induced parasite growth retardation in blood-stage Plasmodium falciparum infection

Pengxing Cao

School of Mathematics and Statistics, The University of Melbourne, Melbourne, Australia.

Nectarios Klonis

Department of Biochemistry and Molecular Biology, Bio21 Molecular Science and Biotechnology Institute, University of Melbourne, Melbourne, Australia.

Sophie Zaloumis

Centre for Epidemiology and Biostatistics, Melbourne School of Population and Global Health, The University of Melbourne, Melbourne, Australia.

David S. Khoury

Infection Analytics Program, Kirby Institute, UNSW Australia, Kensington, New South Wales, Australia.

Deborah Cromer

Infection Analytics Program, Kirby Institute, UNSW Australia, Kensington, New South Wales, Australia.

Miles P. Davenport

Infection Analytics Program, Kirby Institute, UNSW Australia, Kensington, New South Wales, Australia.

Leann Tilley

Department of Biochemistry and Molecular Biology, Bio21 Molecular Science and Biotechnology Institute, University of Melbourne, Melbourne, Australia.

Julie A. Simpson

Centre for Epidemiology and Biostatistics, Melbourne School of Population and Global Health, The University of Melbourne, Melbourne, Australia.

James M. McCaw Correspondence: [email protected] School of Mathematics and Statistics, The University of Melbourne, Melbourne, Australia.

Centre for Epidemiology and Biostatistics, Melbourne School of Population and Global Health, The University of Melbourne, Melbourne, Australia.

Modelling and Simulation, Infection and Immunity Theme, Murdoch Childrens Research Institute, The Royal Children’s Hospital, Parkville, Victoria, Australia.

Abstract

Falciparum malaria is a major parasitic disease causing widespread morbidity and mortality globally. Artemisinin derivatives—the most effective and widely-used antimalarials that have helped reduce the burden of malaria by 60% in some areas over the past decade—have recently been found to induce growth retardation of blood-stage Plasmodium falciparum when applied at clinically relevant concentrations. To date, no model has been designed to quantify the growth retardation effect and to predict the influence of this property on in vivo parasite killing. Here we introduce a mechanistic model of parasite growth from the ring to trophozoite stage of the parasite’s life cycle, and by modelling the level of staining with an RNA-binding dye, we demonstrate that the model is able to reproduce fluorescence distribution data from in vitro experiments using the laboratory 3D7 strain. We quantify the dependence of growth retardation on drug concentration and demonstrate the model’s utility as a platform to propose experimentally-testable mechanisms of growth retardation. Furthermore we illustrate that a drug-induced delay in growth may significantly influence in vivo parasite dynamics, demonstrating the importance of considering growth retardation in the design of optimal artemisinin-based dosing regimens.

Introduction

Plasmodium falciparum malaria is a major parasitic disease which causes severe morbidity and mortality in approximately half a million people annually [1]. Artemisinin (ART) and its derivatives (e.g. artesunate, dihydroartemisinin and artemether), used in combination with partner drugs, provide front-line protection, and have been responsible for dramatic reductions in disease burden over the past few decades [1]. Despite their clinical and public health effectiveness, the emergence of ART resistance and lack of alternative treatments places current control programs at risk [2, 3, 4, 5]. Development of a comprehensive understanding of ART’s mechanism of action and associated effects on infected red blood cells (iRBCs) is therefore critical for development of optimised ART-based treatment regimens and maintenance of control program impact [6].

Recent in vitro experiments [7, 8, 9], combined with advances in pharmacokinetic–pharmacodynamic (PK–PD) modelling [10], have established a platform to probe the parasite’s temporal response to antimalarial drugs. The key experimental advance underlying these in vitro studies was the application of short drug pulses, which enabled fine-scale measurement of the killing effect of drug [10]. A normal life cycle of an iRBC for P. falciparum is approximately 48 hours and is classified based on morphological appearance into three main stages: the ring stage (approximately 0–26 hours post infection (h p.i.)), trophozoite stage (approximately 27–38 h p.i.) and schizont stage (approximately 39–48 h p.i.). Upon rupture at approximately 48 h p.i., iRBCs release merozoites, 8–12 of which successfully invade susceptible RBCs to initiate a new round of infection [11, 12, 13]. Dogovski et al. demonstrated that a short pulse of ART (or dihydroartemisinin (DHA)) can induce growth retardation, prolonging the 48 hour life cycle [8]. Importantly, they found that growth retardation did not stop parasite growth entirely and was thus considered to be distinct from parasite dormancy which “freezes” parasites for days to weeks [14, 15].

Experimental identification of drug-induced growth retardation raises two questions: 1) By how much is the life cycle of a parasite prolonged in response to a short drug exposure pulse?; and 2) How influential is growth retardation when considering in vivo parasite killing using PK–PD models? These two questions are important because we expect any drug-mediated variation in the duration of one or more life stages to impact on the efficacy of the drug (given the well-established finding that drug can exert stage-specific killing effect to parasites [16, 7, 17, 10].) The alteration in drug efficacy may be significant if the prolonged stage(s) covered by the drug pulse exhibit very distinct killing effects. Quantification of growth retardation and assessment of its potential effect on in vivo parasite killing is therefore important in guiding further experimental investigations into drug activity and strategies for optimising ART-based combination therapies.

To the best of our knowledge, no model has yet been designed to quantify the growth retardation effect and predict its influence on in vivo parasite killing. Here we construct a mechanistic model of parasite growth to explain drug-induced growth retardation. We model the relationship between parasite age and the fluorescence intensity of a parasite’s RNA/DNA-binding dye. Through application to fluorescence data from in vitro experiments using the laboratory 3D7 strain by Dogovski et al.[8], we provide the first quantification of growth retardation. We assess its dependence on drug (ART and DHA) concentration. We also suggest specific alternative hypotheses for the mechanism of drug-induced growth retardation and demonstrate that those mechanisms can be reliably identified in future experimental studies. Finally, by incorporating growth retardation into a PK–PD modelling framework, we simulate in vivo parasite killing upon exposure to a single dose of artesunate and show how growth retardation may manifest as a phenotypical indication of drug-resistance.

Materials and Methods

Experiment and data

We first summarise the in vitro experiment in which parasite growth retardation was identified. We provide sufficient information for the purposes of model development and evaluation. For full details on the experimental implementation we refer the reader to the original publication [8].

Fig. 1 presents the experimental process. A culture containing tightly synchronised rings (3D7 strain; over 80% of the population within a one-hour age window) with an average age of 6 h p.i. was equally divided into a number of small cultures, each of which was treated with a different concentration of drug (ART or DHA) for 4 hours. Two cultures unexposed to drug acted as the control. The cultures were stained with SYTO-61 and examined by flow cytometry [18]. SYTO-61 is a nucleic acid stain which stains both DNA and RNA and allows the distinction of infected RBC from uninfected RBC, as well as distinguishing between parasite-infected RBC of different ages due to an increase in nucleic acid content as the parasite ages. SYTO-61 signals from all cultures were collected simultaneously at approximately 72 hours post drug administration (indicated in Fig. 1), corresponding to the period in which parasites transition from the ring to trophozoite stage (during the second cycle). As trophozoites express significantly more nucleic acid, this period of transition exhibited a bimodal SYTO-61 fluorescence distribution. Quantitative analysis of this bimodal distribution is the key to our approach as it enables us to identify the relative populations of parasites in the two successive life stages [8]. As indicated in Fig. 1, if the administration of drug slows parasite growth then the first life cycle would be prolonged and the second life cycle would start later. In consequence, we would observe an increase in the ring population (i.e. the mode with the lower fluorescence) and a corresponding decrease in the trophozoite population (i.e. the mode with the higher fluorescence) in the SYTO-61 fluorescence histogram. Note that the underlying processes of iRBC rupture and merozoite release and re-infection were not observable in the experiment.

The SYTO-61 fluorescence histograms often contain a small population of dead or dormant parasites (due to drug activity or otherwise) that are not involved in the second life cycle. To account for this non-viable population, a background experiment in which supermaximal drug concentration ( $>10\times$ the 50% Lethal dose (3 days)) was applied for over 48 hours was also performed [8]. The high drug concentration and long exposure time guaranteed that all parasites became non-viable. Hence, denoting the SYTO-61 fluorescence frequency histogram under the background condition by $f_{b}$ and the SYTO-61 fluorescence frequency histogram under a 4h drug pulse by $f$ , the corrected SYTO-61 fluorescence frequency histogram, $f_{c}$ , is given by

[TABLE]

where $V$ represents the viability (the fraction of parasites entering the second life cycle; see [7, 8] for details).

The corrected histogram data is shown in Fig. 2 (for various ART concentrations) and Supplementary Fig. S1 (for various DHA concentrations). Experiments were performed in technical replicates for each drug concentration. The histograms in Fig. 2 present all available SYTO-61 fluorescence intensity data. Also note that for each histogram, the samples with fluorescence intensity less than 3000 (indicated by the vertical dashed lines) were considered to include fluorescence signals from uninfected RBCs and were thus excluded in the analysis. Each histogram was generated by distributing log-transformed SYTO-61 fluorescence intensity samples (with magnitude $>$ 3000) into 40 equally spaced bins. Raw SYTO-61 fluorescence intensity data is provided in Dataset S1.

The model

Since the in vitro experiment measures the SYTO-61 fluorescence data in the second life cycle, our model is designed to reproduce the underlying process of parasite growth over the period of the ring-to-trophozoite transition in the second life cycle (see Fig. 1).

Over the period of the ring-to-trophozoite transition, the growth of individual parasites is modelled by two sequential stages — an “immature” ring stage during which the ring-to-trophozoite transition cannot occur (due to incomplete cellular development) followed by a “mature” ring stage where the ring-to-trophozoite transition is possible — inspired by the classic model for the mammalian cell cycle [19, 20]. We introduce $A_{r}$ , the ready-for-change age (which is assumed to be the same for all parasites). Parasites of age $a<A_{r}$ are, by definition, rings. For parasites of age $a>A_{r}$ , the waiting time before entering the trophozoite stage follows a Poisson distribution with transition rate $\lambda$ (denoted as $Pois(\lambda)$ ).

We assume that the age distribution of viable parasites (i.e. the parasites able to asexually reproduce in the second life cycle) at the time of data collection is Gaussian ( $\sim N(\mu,\sigma^{2})$ ). This is reasonable given that the parasites are tightly synchronised and drug is unlikely to differentially kill parasites over such a tight age window (although a spreading in age distribution over time was evident [8].)

The bimodal fluorescence intensity distributions (Fig. 2) suggest that both the SYTO-61 fluorescence intensity and the rate of increase in SYTO-61 fluorescence intensity is higher in trophozoites than in rings. We model this property using a piecewise function mapping from parasite age ( $a$ ) to SYTO-61 fluorescence intensity ( $F$ ) as follows:

[TABLE]

where $A_{c}$ indicates the age of transition (from ring to trophozoite) for a parasite, noting that each individual parasite’s $A_{c}$ is sampled from Pois( $\lambda$ ). Clearly $A_{c}\geq A_{r}$ for all parasites. For each parasite, prior to the transition (i.e. $a<A_{c}$ ), the age-dependent intensity of SYTO-61 staining increases at a rate $r_{1}$ . Following the transition the rate increases to $r_{2}$ . We note that although $F_{0}$ gives the fluorescence intensity at the start of the second life cycle (i.e. age 0 h p.i.), the model is not designed to study the early intraerythrocytic life-stages of the parasites.

All model parameters are listed in Table 1. We note two important limiting cases of the model. When $A_{r}=0$ , the model allows rings to transition to trophozoites immediately (at 0 h p.i.), with transition rate $\lambda$ . On the other hand, when $\lambda\to\infty$ , all rings reaching age $A_{r}$ immediately undergo the transition to the trophozoite stage. Therefore, our model provides a very general framework in which to study the statistics of the ring-to-trophozoite transition and the mechanism of growth retardation.

In preparation for its application to the experimental data, Fig. 3 presents a single stochastic realisation (performed using the Gillespie algorithm [21]) of the model. We initiated the simulation with a population of $4000$ ring-stage parasites with a mean age of 10 h p.i. (and standard deviation of 0.8 h p.i.) in the second life-cycle. We assumed that $A_{r}=21\ \rm h\ p.i.$ , $\lambda=0.7\ \rm h\ p.i.^{-1}$ , $F_{0}=20$ (unitless), $r_{1}=0.24\ \rm h\ p.i.^{-1}$ and $r_{2}=0.56\ \rm h\ p.i.^{-1}$ . Note that we chose these parameters simply to illustrate the model’s behaviour. While the simulated fluorescence histograms may look similar to the experimental data, the parameters do not necessarily represent the true values. The simulation was implemented in MATLAB (version R2014b; The MathWorks, Natick, MA) and the code is provided in the Supporting Information. The left panels in Fig. 3 show the continuous ageing of the parasite population and accumulation of trophozoites within the simulated population. The right panels show how the transition to a population dominated by trophozoites (with their corresponding increased SYTO-61 staining) leads to a distinct shift in the SYTO-61 fluorescence intensity histogram. The bimodality evident in the observed data (Fig. 2) is also evident in the simulated data.

Derivation of the probability density function for the SYTO-61 fluorescence distribution

Our model for the transition from ring to trophozoite stage can be represented graphically by mapping parasite age to SYTO-61 fluorescence intensity (Fig. 4). The intensity signal for a single ring-stage parasite will increase with age along the line with slope $r_{1}$ . On becoming a trophozoite (at a certain age $A_{c}>A_{r}$ ) it will then follow one of the green lines with slope $r_{2}$ . It follows that a fluorescence intensity $F_{1}$ can only originate from rings with age $a_{1}$ . However, a fluorescence intensity $F_{2}$ can result from trophozoites with ages between $a_{2}$ and $a_{2}^{\prime}$ (inclusive) or a ring with age $a_{2}^{\prime}$ . These considerations suggest that we derive the probability density function for the total fluorescence signal (from a population of parasites) in a piecewise way with a critical point of separation, $F_{r}$ , mapped from the ready-for-change age $A_{r}$ .

Case 1: $F\leq F_{r}$ where $F_{r}=F_{0}e^{r_{1}A_{r}}$ . Since $a_{1}\sim N(\mu,\sigma^{2})$ and $a_{1}={\rm ln}(F/F_{0})/r_{1}$ , we have the probability density function for $F$ ,

[TABLE]

where $1/(r_{1}F)$ is the derivative of $a_{1}$ with respect to $F$ due to use of the logarithmic transform on the fluorescence data.

Case 2: $F>F_{r}$ where $F_{r}=F_{0}e^{r_{1}A_{r}}$ . The probability density function for $F$ is a sum of two parts. One is the probability of rings with age $a_{2}^{\prime}$ and the other is the probability of trophozoites with a range of ages from $a_{2}$ to $a_{2}^{\prime}$ . The former is given by

[TABLE]

The first term in this expression comes from the normal distribution of the population. The second term follows from the modelled Poisson distribution, $Pois(\lambda)$ , for the probability that a parasite of age $a_{2}^{\prime}$ remains in ring form. The last term is the derivative of $a_{2}^{\prime}$ with respect to $F$ due to the logarithmic transform $a_{2}^{\prime}={\rm ln}(F/F_{0})/r_{1}$ .

The contribution to the probability density from trophozoites is given by

[TABLE]

where $Pr(F|a)$ is the probability density function of $F$ from trophozoites with age $a$ . To derive $Pr(F|a)$ , we first give the cumulative probability

[TABLE]

where $A_{c}(F)=[r_{2}a-{\rm ln}(F/F_{0})]/(r_{2}-r_{1})$ and $A_{c}(F+\Delta F)=[r_{2}a-{\rm ln}((F+\Delta F)/F_{0})]/(r_{2}-r_{1})$ . The difference of two exponentials indicates the probability of trophozoites with conversion age between $A_{c}(F)$ and $A_{c}(F+\Delta F)$ . By multiplying $1/\Delta F$ on both sides of Eq. 6 and taking the limit $\Delta F\to 0$ , we obtain

[TABLE]

where $A_{c}=[r_{2}a-{\rm ln}(F/F_{0})]/(r_{2}-r_{1})$ . We notice that the second term $\lambda e^{-\lambda(A_{c}-A_{r})}$ is exactly the probability density function of the Poisson distribution, and the last term is the derivative of $A_{c}$ with respect to $F$ (except for a negative sign cancelled during simplification).

Hence, the probability density function for $F$ is given by

[TABLE]

where $Pr(F|a)$ is given by Eq. 7, $F_{r}=F_{0}e^{r_{1}A_{r}}$ , $a_{1}={\rm ln}(F/F_{0})/r_{1}$ , $a_{2}=[{\rm ln}(F/F_{0})+(r_{2}-r_{1})A_{r}]/r_{2}$ , $a_{2}^{\prime}={\rm ln}(F/F_{0})/r_{1}$ and $A_{c}=[r_{2}a-{\rm ln}(F/F_{0})]/(r_{2}-r_{1})$ . Although $a_{1}$ and $a_{2}^{\prime}$ have the same expression, we retain both for clarity.

Parameter identifiability

To apply the analytical expression for the SYTO-61 fluorescence distribution (Eq. 8) to data, we must determine if the equation is well-posed, i.e. is there sufficient information within the data to uniquely identify the model’s parameters. Eq. 8 contains seven parameters, but as we will now show, this equation can be reduced to an expression containing just five parameters. This expression cannot be further reduced and so is sufficient for parameter estimation (Table 2). The relationship between these five parameters and the seven “biological” parameters—with which we are fundamentally concerned—is derived.

We begin with the case $F<F_{r}$ and identify two new parameters,

[TABLE]

$Pr(F)$ can be rearranged to be a two-parameter distribution

[TABLE]

where $K_{1}$ and $K_{2}$ represent the mean and standard deviation respectively. This is consistent with the fact that a log-normal distribution is uniquely determined by two parameters.

Secondly, for the case of $F>F_{r}$ , the first part is also a log-normal distribution except for a Poisson probability density function. Thus we introduce

[TABLE]

such that the first part is given by

[TABLE]

Thirdly, for the integral part, we introduce

[TABLE]

Therefore, Eq. 7 becomes

[TABLE]

where

[TABLE]

Then taking the transform $b=r_{1}a+{\rm ln}(F)$ , the integral part in Eq. 8 becomes

[TABLE]

where

[TABLE]

and the lower and upper bounds are given by

[TABLE]

At last, we find $F_{r}=e^{K_{3}}$ which finalises the process of expressing the fluorescence distribution using five parameters, each of which should be identifiable by application of the model to the fluorescence intensity histogram data. The final form of the model is given by

[TABLE]

where the integral of $Pr(F|b)$ is given by Eq. 17. A summary of the five new parameters is provided in Table 2.

Method for data fitting

The simplified 5-parameter probability density function (Eq. 18) was fitted to the observed SYTO-61 fluorescence data for the 3D7 strain shown in Fig. 2 and Fig. S1. While drug induces growth retardation of parasites in the first life cycle (from the time of drug exposure to the end of that cycle), the second life cycle is drug-free and assumed to progress normally. Accordingly, when modelling parasites in the second life cycle, the model parameter that may depend upon drug concentration is the mean parasite age $\mu$ (and thus $K_{1}$ in the 5-parameter model). Moreover, since the parasite age distribution may change over time in the experiment, the spread of parasite age distribution $\sigma$ (and in turn $K_{2}$ ) is also allowed to vary. To simultaneously fit the model to all the data shown in Fig. 2 (or Fig. S1), we therefore require 35 parameters (16 for $K_{1}$ , 16 for $K_{2}$ , and one for each of $K_{3}$ , $K_{4}$ and $K_{5}$ ).

To optimise the data fitting procedure, we provide plausible ranges for each parameter (summarised in Table 1). Based on the observed fluorescence range of approximately 500–50000 in the observed data, we have $e^{K_{1}}\in[500,50000]$ which gives $K_{1}\in[6.2,10.8]$ . We estimate $\sigma\in[0.39,2]$ The lower bound is based on the the estimated standard deviation of the initial parasite age distribution where approximately 80% of parasites were aged within a one-hour window centred at the mean age. The upper bound is chosen to allow for a much wider age distribution. $r_{1}$ is assumed to be in $[0.01,0.5]$ where the lower bound avoids a very small denominator in $K_{4}$ and $K_{5}$ . These give $K_{2}\in[0.0039,1]$ . We have $K_{3}\in[2,10.8]$ based on a similarity to $K_{1}$ but with a much smaller lower bound in order to capture possible small values of $A_{r}$ . $K_{4}\in[0,300]$ corresponds to $\lambda\in[0,3]$ where the upper bound $\lambda=3$ represents a very fast ring-to-trophozoite transition (over 95% of rings become trophozoites within one hour of ageing to the ready-for-change age $A_{r}$ ). We have $K_{5}\in(1,10]$ where the lower bound is based on the assumption that $r_{2}>r_{1}$ and the upper bound is arbitrarily chosen.

A global search using a combination of Latin Hypercube Sampling (LHS) and non-linear least squares optimisation was performed. A minimum of 4000 parameter sets (each of which contains 35 parameter values, 16 identical values for $K_{1}$ , 16 identical values for $K_{2}$ , and one value for each of $K_{3}$ , $K_{4}$ and $K_{5}$ ) were randomly generated using MATLAB’s built-in LHS function lhsdesign (plus a few additional sets proposed based on our prior knowledge) and were used as initial points for multi-start optimisation. The precision for the initial values was up to the first decimal place for $K_{2}$ and integer values were sampled for all other parameters. For each of these initial points (minimum 4000) generated using the LHS sampling, MATLAB’s built-in functions lsqcurvefit and nlparci with default settings were used to minimise the sum of squared residuals subject to the above constraints and produce the corresponding best-fit parameter values and associated 95% confidence intervals (CIs). The globally optimal solution and associated best-fit parameter values (which are presented in the Results) were given by the parameter set producing the smallest sum of squared residuals. The trapezoidal integration method with a division of 1000 equal sub-intervals was used to approximate the integral over $[b_{2},b_{2}^{\prime}]$ in the probability density function.

Results

Quantifying the dependence of drug-induced parasite growth retardation on drug concentration

Fig. 5 presents the results of fitting the 5-parameter model to the SYTO-61 fluorescence data for ART (the best-fit parameter values are provided in Supplementary Table S1). The model correctly captures the ART-dependent change of shape in SYTO-61 fluorescence histograms — increasing ART concentration leads to an increase in the ring population (i.e. the mode with the lower fluorescence) at the expense of a decreased trophozoite population (i.e. the mode with the higher fluorescence).

To quantify the extent of ART-induced parasite growth retardation, we define $T_{\Delta}=\mu-A_{r}$ , the time period (in the second life cycle) from the ready-for-change age $A_{r}$ to the mean population age at the time of data collection $\mu$ . Recapping Fig. 1, since the delay in rupture time of the first cycle directly affects the mean parasite age at the time of data collection $\mu$ , but not the ready-for-change age $A_{r}$ , the change in $T_{\Delta}$ therefore provides a direct indication of how much the rupture time is delayed in the first life cycle due to drug application. A smaller $T_{\Delta}$ indicates a longer delay in the rupture time (in the first cycle), i.e. parasites are younger when measured. Conversely, a larger $T_{\Delta}$ indicates a shorter delay, i.e. parasites are relatively older at the time of data collection. Given that $r_{1}T_{\Delta}=K_{1}-K_{3}$ and that the age-dependent SYTO-61 staining rate $r_{1}$ is independent of drug concentration, we can use $K_{1}-K_{3}$ to quantify the dependence of growth retardation on drug concentration. A smaller value for $K_{1}-K_{3}$ indicates a longer delay in the rupture time (in the first cycle). Fig. 6A shows that, for an ART concentration of less than 5 nM (applied as a 4h pulse), little if any growth retardation is evident. However, when ART concentration increases above 5 nM, $K_{1}-K_{3}$ (and so $T_{\Delta}$ ) decreases approximately linearly, implying an approximately linear positive correlation between the delay to rupture time (for the first cycle) and drug concentration.

For clinically relevant DHA, the model performs similarly well (Fig. 7). Similar to ART, we observe a region where parasite growth is almost unaffected (e.g. DHA less than 1 nM; see Fig. 8A) followed by a region where $K_{1}-K_{3}$ decreases monotonically with increasing DHA concentration. However, unlike for ART, the relationship is non-linear. For high drug concentration, $K_{1}-K_{3}$ ceases to decrease (although this might be inaccurate due to limited numbers of samples caused by low viabilities), indicating some sort of saturation effect.

Although we cannot determine $T_{\Delta}$ from $K_{1}-K_{3}$ due to the unknown parameter $r_{1}$ (the rate of age-dependent SYTO-61 staining in ring stage), we can provide some useful bounds based on plausible values for $r_{1}$ . For example, if $r_{1}$ is 0.1 $\rm[h\ p.i.]^{-1}$ then an increase in ART concentration from 0 nM to approximately 1000 nM leads to a decrease of $K_{1}-K_{3}$ from approximately 0.5 to 0.26 (Fig. 6A), implying a delay in rupture time of approximately $T_{\Delta}=(0.5-0.26)/0.1=2.4$ hours. Similarly, for the same $r_{1}$ value but an increase in DHA concentration from 0 nM to 100 nM (see Fig. 8A), the estimated delay in rupture time is approximately $((0.5-0.1)/0.1=)$ 4 hours, which is longer than that induced by 1000nM ART. Clearly, if $r_{1}$ is smaller (or larger), the actual rupture time will be delayed (advanced) accordingly. Based on a reasonable range for $r_{1}$ of [0.05–0.25] $\rm h\ p.i.^{-1}$ , the drug-induced delay in rupture time (for the first life cycle) should be no more than 10 hours.

The parasite age distribution is relatively unaffected by the drug pulse

The parameter $K_{2}$ ( $=\sigma r_{1}$ ) directly indicates how the standard deviation (or spread) of the parasite age distribution at the time of data collection ( $\sigma$ ) depends on drug concentration. As shown in Figs. 6B and 8B, there is a weak positive correlation for relatively low drug concentrations. For relatively high drug concentrations, although $K_{2}$ drops substantially for DHA (Fig. 8B, which might be strongly affected by small sample size), no substantial decrease is observed for ART (Fig. 6B). These results suggest that, although varying drug concentration may affect the spread of the parasite age distribution, the effect is usually very limited.

New experiments to test alternative mechanisms for drug-induced growth retardation

Having established that our model can be used to quantify the relationship between applied drug concentration and growth retardation for the in vitro data, we now explore the utility of the model for identifying possible alternative mechanisms responsible for growth retardation. Focusing on the ring-to-trophozoite transition (which is the core part of the model), one hypothesis is that growth retardation may be a result of either a decreased transition rate $\lambda$ or a postponed ready-for-change age $A_{r}$ for parasites in the first life cycle (and thus delayed rupture time and initiation of the second life cycle). To test the hypothesis, we would require fluorescence intensity measurements for viable parasites to be performed during the first life cycle, particularly around the period of the ring-to-trophozoite transition. Such an experimental setup is illustrated in Fig. 9. We note that the in vitro data used thus far, in which fluorescence measurements were taken in the second life cycle, were primarily designed to assess viability (which is, of course, only measurable through study of the second life cycle.) and so are not suitable for testing this hypothesis. We have therefore taken an in silico approach in anticipation of future experimental studies.

First, suppose that increasing ART concentration decreases the ring-to-trophozoite transition rate $\lambda$ (while leaving the ready-for-change age $A_{r}$ unchanged). Fig. 10 shows simulated data and the model fit for four different values of the transition rate $\lambda$ between 0.3 and 1.2. Simulated data is in black (generated by assuming 4000 parasites were present) and the best-fit model in red (dashed lines). The true (model-based) distribution from which the simulated data were obtained is shown in blue. Fig. 11 presents estimates for $K_{3}$ and $K_{4}$ , which determine $A_{r}$ and $\lambda$ respectively (Table 2). The model fitting procedure correctly identifies that $A_{r}$ , the ready-for-change age, was unchanged across the simulations, while $\lambda$ , the ring-to-trophozoite transition rate, was increased. We note the presence of some minor bias in the estimate for $K_{3}$ for this particular simulated dataset. Supplementary Fig. S2 presents summary statistics for 100 replications of the simulation-reestimation procedure and establishes that unbiased estimates were obtained. The relative bias for all 11 parameters ( $K_{1}$ , $K_{2}$ , $4\times K_{3}$ , $4\times K_{4}$ , $K_{5}$ ) was very small, the largest being 1.03%).

In addition to varying $\lambda$ only, we also examined the other three logical possibilities: 1) $\lambda$ is fixed and $A_{r}$ is varied (see Figs. S3 and S4; Supplementary Table S4); 2) both $\lambda$ and $A_{r}$ are varied but in opposite directions, i.e. one increases while the other decreases (see Figs. S5 and S6; Supplementary Table S5); 3) both $\lambda$ and $A_{r}$ are varied but in the same direction, i.e. both increase or decrease (see Figs. S7 and S8; Supplementary Table S6). In all situations, the results consistently show that the possible variations in the key parameters $\lambda$ and $A_{r}$ are identified for a given set of simulated fluorescence intensity histogram data. Hence, with future availability of fluorescence data collected during the first cycle, we anticipate that our model will reliably distinguish between the alternative hypothesised mechanisms driving growth retardation.

Predicting the effect of growth retardation on in vivo parasite killing

Although a short pulse (4 hours) of ART/DHA may prolong the parasite’s life cycle by a relatively short time (e.g. less than 10 hours as estimated above), it could have a significant effect on in vivo parasite killing. For example, if the delay occurs during the ring-to-trophozoite transition (as hypothesised above), then the recent identification that ring-stage parasites exhibit a much lower sensitivity to ART (and DHA) than those in the trophozoite stage [16, 7, 17] directly suggests that delayed ageing may lead to effective drug escape. Accordingly, here we perform a preliminary investigation of the possible impact of growth retardation on in vivo parasite killing using a simple PK–PD model. The model contains just two compartments representing populations of rings and trophozoites, with a transition rate from rings to trophozoites given by the rate $\lambda$ . Rings are killed by ART/DHA at a slower rate than trophozoites, as suggested by previous findings [16, 7, 17]. Model details are provided in the Supporting Information.

Simulation results are shown in Fig. 12. A population of ring-stage parasites were initially distributed with a mean age of 10 h p.i. (note, this is simply the time since the simulated population of RBCs were infected, not the time of clinical exposure for the simulated host) and a standard deviation of 2 h at the start of simulation. With a single dose of artesunate (2mg/kg) applied at 10 hours (see the middle panel of Fig. 12), the rate of parasite death for a scenario in which there is no drug-induced growth retardation (i.e. $\lambda$ is independent of drug concentration) is significantly higher than that for a scenario in which there is drug-induced growth retardation (i.e. $\lambda$ is a function of drug concentration). Growth retardation delays the transition from rings to trophozoites and in turn prevents those parasites from being killed as the drug’s short half-life (approximately 0.9 h) allows them to avoid exposure during the trophozoite stage. This result suggests that, if drug induces a slower transition from ring to trophozoite stage, growth retardation has a substantial adverse effect on efficient parasite killing and may thus be considered as a potential mechanism for decreasing ART sensitivity.

Discussion

In this paper, we have studied ART-induced parasite growth retardation using a mechanistic model that considers the ring-to-trophozoite transition to be a two-stage process and exploits the differing rates of nucleic acid production (measured through SYTO-61 fluorescence intensity) in those two life stages. By fitting the model to fluorescence histogram data, we have been able to identify the dependence of growth retardation on applied drug concentration (Figs. 6A and 8A) and obtain reliable estimates for how much the parasite’s life span is prolonged due to exposure to drug. Our primary findings from this analysis were that: 1) drug-induced parasite growth retardation exhibits a threshold-like behaviour such that growth retardation is evident only when drug concentration is sufficiently large (for example, $>10\ \rm nM$ ART for a 4 hour drug pulse); 2) the parasite life cycle is prolonged by no more than 10 hours due to application of ART/DHA; and 3) the parasite age distribution at the time of data collection in the second life cycle is relatively unaffected by the application of short drug pulses with different concentrations.

Furthermore, we used the model to propose a hypothetical mechanism of ART-induced parasite growth retardation. We considered growth retardation to be due to either a decreased transition rate $\lambda$ or a postponed ready-for-change age $A_{r}$ , or a combination of the two. We have shown that, if fluorescence intensity data were collected in the first life cycle, that the model is able to accurately identify the dependencies of the transition rate $\lambda$ and the ready-for-change age $A_{r}$ on drug concentration. With the availability of new experimental data (e.g. fluorescence data from the first life cycle), we will be able to incorporate growth retardation into our recently developed dynamic-stress model of antimalarial action [10], providing a comprehensive platform for the study and optimisation of ART-based therapies.

ART-induced growth retardation has two competing effects. Firstly, it slows the overall rate of growth by extending the life cycle. However, by delaying the transition to the trophozoite stage, the parasite avoids being exposed to drug at a highly sensitive stage [7, 10], providing a net benefit (to the parasite) due to the delay (as shown in Fig. 12). Our preliminary in vivo simulations suggest that the overall effect is likely to be strongly beneficial (for the parasite): a short delay in the time of ring-to-trophozoite transition will have a minimal impact on the overall growth rate (and so any implications for onwards transmission are likely minimal), yet is extremely effective in avoiding the short ART/DHA pulse and so killing. Therefore, understanding where in the parasite life cycle drug-induced growth retardation acts is important in the context of combating emergent drug resistance through optimisation of dosing regimens, as demonstrated by our in vivo simulations (Fig. 12). Furthermore, these results suggest that developing new longer-lived ART derivatives may overcome resistance to ART and DHA, as a direct consequence of their longer half lives.

Returning to the model we have introduced to analyse the SYTO-61 fluorescence data, it provides a major advance on current methods used to study mixed histograms from fluorescence assays. Current practice is to simply partition the bimodal signal into “low” and “high” components at a chosen threshold intensity. While suitable for well separated peaks, this method is inadequate for determining mixtures of rings and trophozoites as shown in Fig. 2. Our model clearly dissects the contributions of different parasite populations to the fluorescence distribution (Eq. 18) where the integral part represents the contribution from trophozoites while the sum of all other parts represents the contribution from rings. It can therefore be used to estimate the fraction of parasite subpopulations based on fluorescence intensity distribution in a far more reliable and rigorous way. Moreover, given that fluorescence dye staining is a standard assay in experimental biology, we anticipate that our novel methods will be broadly applicable to other problems, for example in the application to the adoptive transfer data utilised in [22].

Funding information

The work was supported by the National Health and Medical Research Centre of Australia (NHMRC) through Project Grants 1100394 and 1060357, the Centre for Research Excellence ViCBiostat (1035261) and the Centre for Research Excellence PRISM2 (1078068). James M. McCaw was supported by an Australian Research Council (ARC) Future Fellowship. Julie A. Simpson was supported by a NHMRC Senior Research Fellowship. Leann Tilley was supported by an ARC Professorial Fellowship.

Competing interests

We have no competing interests.

Supporting Information

The supporting information contains the following:

MATLAB code for the stochastic simulation shown in Fig. 3 in the main text 2. 2.

The model used to simulate in vivo parasite clearance 3. 3.

MATLAB code for solving the model of in vivo parasite clearance 4. 4.

Supplementary figures (S1–S8) 5. 5.

Supplementary tables (S1–S6)

MATLAB code for the stochastic simulation shown in Fig. 3 in the main text

⬇

1clear

2clc

3tic

4

5Ar = 21; % ready-for-change age

6total_Num = 4e+3; % total number of viable paraistes

7

8init_temp = 0.8*randn(1,total_Num)+10;

9init_age = init_temp(init_temp>=0); % initial age distribution

10age = sort(init_age);

11stage = (age>=Ar)+1;

12% stage=1: rings before age Ar; stage=2: rings after Ar; stage=3: troph

13index = 1:length(age);

14

15rand1 = rand(1,length(stage)).*(stage == 2);

16Xi = log(1./rand1);

17tracking = zeros(1,length(rand1)); % tracking variable

18

19lambda0 = 0.7; % ring-to-troph transition rate

20lambda = (stage == 2).*lambda0;

21age_to_change = tracking; % the actual age when a ring change to troph

22

23r1 = 0.24;

24r2 = 0.56;

25F0 = 20;

26

27SYTO61 = F0exp(r1age);

28

29age0=age;

30stage0=stage;

31SYTO610=SYTO61;

32

33dt = 0.01;

34

35t = 0:dt:20;

36

37time_to_record = 11:15; % record the results of mean age=21,22,23,24,25

38

39indext=1:length(time_to_record);

40

41age_all=ones(length(time_to_record),1)*age;

42stage_all=ones(length(time_to_record),1)*stage;

43SYTO61_all=ones(length(time_to_record),1)*SYTO61;

44

45wb=waitbar(0,’please wait’);

46

47for i = 2:length(t)

48

49 % update age

50 age = age+dt;

51

52 % update the status for those whose ring-to-troph transition occurs

53 tracking = tracking + lambda*dt;

54 trs = (tracking-Xi>=0); % transition to troph occurs

55 stage = stage+trs;

56 age_to_change = age_to_change+trs.*age;

57 tracking = tracking.*(1-trs); % update tracking variable

58 rand1 = rand1.*(1-trs);

59

60 % update the status of the rest parasites

61 temp = age.*(stage == 1);

62 indtemp=index(temp>=Ar);

63 rand1(indtemp) = rand(1,length(indtemp));

64 Xi = log(1./rand1);

65 stage(indtemp) = 2;

66 lambda = (stage == 2).*lambda0;

67

68 % SYTO61 fluo cumulation

69 SYTO61 = (stage<3).(F0exp(r1*age))+…

70 (stage==3).(F0exp(r1*age_to_change).exp(r2(age-age_to_change)));

71

72 % recording results

73 indt=indext(time_to_record==t(i));

74 if ~isempty(indt)

75 age_all(indt,:) = age;

76 stage_all(indt,:) = stage;

77 SYTO61_all(indt,:) = SYTO61;

78 end

79

80 waitbar(i/length(t))

81end

82

83close(wb)

84

85toc

After running the above code, the following code is used to generate the result of Fig. 3 in the main text.

⬇

1for kk=1:5;

2

3 bin_edge = 0:0.2:48;

4

5 ydata_age=histcounts(age_all(kk,:),bin_edge);

6 temp = age_all(kk,:);

7 temp1 = stage_all(kk,:);

8 troph_age = temp(temp1==3);

9 ydata_trogh_age=histcounts(troph_age,bin_edge);

10

11

12 % for plot only

13 xnodes = zeros(1,2*length(bin_edge));

14 xnodes(1:2:end) = bin_edge;

15 xnodes(2:2:end) = bin_edge;

16

17 ynodes = zeros(1,2*length(bin_edge));

18 ynodes(2:2:end-1) = ydata_age;

19 ynodes(3:2:end-1) = ydata_age;

20

21 ynodes1 = zeros(1,2*length(bin_edge));

22 ynodes1(2:2:end-1) = ydata_trogh_age;

23 ynodes1(3:2:end-1) = ydata_trogh_age;

24

25 figure(1)

26 subplot(5,2,1+2*(kk-1))

27 plot(xnodes,ynodes1,’-’,’linewidth’,2,’color’,[0 1 0])

28 hold on

29 plot(xnodes,ynodes,’-’,’linewidth’,2,’color’,[0 0 0])

30

31 % xlabel(’parasite age (h p.i.)’)

32 % ylabel(’No. of events’)

33 box off

34

35 set(gca,’xlim’,[16 30],’ylim’,[0 500],’fontsize’,14)

36 set(gca,’TickDir’,’out’,’LineWidth’,2)

37 set(gca,’ticklength’,[0.02 0.01])

38 text(20,200,[’mean parasite age = ’,num2str(time_to_record(kk)+10),’ h p.i.’],’fontsize’,14)

39

40

41 % SYTO fluo histogram

42

43 bin_edge = exp(linspace(log(min(SYTO61_all(kk,:))),log(max(SYTO61_all(kk,:))),50));

44

45 ydata=histcounts(SYTO61_all(kk,:),bin_edge);

46

47 temp = SYTO61_all(kk,:);

48 temp1 = stage_all(kk,:);

49 troph_fluo = temp(temp1==3);

50 ydata_trogh_fluo=histcounts(troph_fluo,bin_edge);

51

52 % for plot only

53 xnodes = zeros(1,2*length(bin_edge));

54 xnodes(1:2:end) = bin_edge;

55 xnodes(2:2:end) = bin_edge;

56

57 ynodes = zeros(1,2*length(bin_edge));

58 ynodes(2:2:end-1) = ydata;

59 ynodes(3:2:end-1) = ydata;

60

61 ynodes1 = zeros(1,2*length(bin_edge));

62 ynodes1(2:2:end-1) = ydata_trogh_fluo;

63 ynodes1(3:2:end-1) = ydata_trogh_fluo;

64

65 subplot(5,2,2*kk)

66 semilogx(xnodes,ynodes1,’-’,’linewidth’,2,’color’,[0 1 0])

67 hold on

68 semilogx(xnodes,ynodes,’-’,’linewidth’,2,’color’,[0 0 0])

69

70 % xlabel(’SYTO-61 fluorescence’)

71 box off

72

73 set(gca,’xlim’,[1000 100000],’ylim’,[0 400],’fontsize’,14)

74 set(gca,’TickDir’,’out’,’LineWidth’,2)

75 set(gca,’ticklength’,[0.02 0.01])

76end

The model used to simulate in vivo parasite clearance

Here we propose a simple two-compartment model to capture both parasite killing and transition from rings to trophozoites. For a certain age $a$ , we denote the number of live rings by $R$ and the number of live trophozoites by $T$ . Then the equation governing the dynamics of $R$ , $T$ and $a$ are given by

[TABLE]

where the term $\lambda R$ represents the rate of conversion from rings to trophozoites at the population level. Since we assume the transition can only occur after age $A_{r}$ ( $A_{r}=21\ \rm h\ p.i.$ is assumed in the simulation), $\lambda$ takes

[TABLE]

where $C$ is DHA concentration. The hyperbolic relationship is used to model a hypothetical property that a higher drug concentration would lead to a stronger delay in parasite growth. $k_{r}$ and $k_{t}$ are drug-induced parasite killing rates for rings and trophozoites respectively. In order to capture that killing rate for rings is significantly smaller than that for trophozoites [16, 7, 17], we choose

[TABLE]

To model the pharmacokinetic profile where plasma DHA concentration usually follows a biphasic behaviour [16], we use

[TABLE]

where $C_{max}$ is the maximum achievable concentration and $t_{m}$ indicates the time (since drug application) when the maximum concentration is achieved. The in vivo DHA half-life $t_{1/2}=0.9\ \rm h$ [9]. To simulate the case of a single dose of artesunate (2mg/kg), we assume $C_{max}=2820\ \rm nM$ and $t_{m}=1\ \rm h$ [9, 2].

To assess the effect of growth retardation on parasite clearance, we consider two scenarios, drug-dependent $\lambda$ and drug-independent $\lambda$ . The former is modelled by using Eq. S4 while the latter is modelled by fixing $C$ to be zero in Eq. S4 (i.e. assuming constant $\lambda=22/34=0.647\ \rm h^{-1}$ for $a\geq A_{r}$ ).

To simulate the model, we initially distribute $10^{11}$ parasites into 20 age bins from age 1 h p.i. to 20 h p.i. (we only consider integer ages) following a normal distribution with a mean of 10 h p.i. and a standard deviation of 2 h p.i. (see Fig. 12 in the main text). For each age group, we use MATLAB’s bulit-in solver ode15s to solve the model and obtain the time series of $T+R$ . Then the time series of total number of live parasites is given by the sum of all 20 time series (result is shown in Fig. 12 in the main text). Moreover, we can also obtain the age distribution at any time (an example of $t=20\rm\ h$ is given in Fig. 12 in the main text). MATLAB code is provided in the next section.

MATLAB code for solving the model of in vivo parasite clearance

⬇

1clear

2tic

3

4drug_app_time = 10; % time of drug application

5

6th=0.9; % in vivo DHA half-life

7

8init_age=10+2*randn(1e+11,1); % initial ages of the 10^11 parasites

9age_dis = histcounts(init_age,0:1:20); % initial age distribution

10age = 1:20; % age bins

11

12t = 0:0.01:30;

13

14% calculate the time series of DHA concentration

15t1 = t(t<drug_app_time); % time before drug application

16t3 = t(t>=drug_app_time+1); % decreasing phase of PK profile

17t2 = intersect(setdiff(t,t1),setdiff(t,t3)); % increasing phase of PK profile

18C1 = 0*t1;

19C2 = 2820/1*(t2-drug_app_time);

20C3 = C2(end)exp(-log(2)/th(t3-t2(end)));

21C = [C1,C2,C3]; % full drug concentration profile

22

23Numt = ones(length(age),length(t)); % matrix for total number of parasites

24Numt(:,1) = age_dis’; % assign the initial age distribution

25

26wb=waitbar(0,’please wait…’);

27

28for kk=1:length(age)

29 init = [age_dis(kk);0;age(kk);0];

30 for i=2:length(t)

31 [~,Sol] = ode15s(@invivo_sim,[0 0.01],init);

32 init = [Sol(end,1:3)’;C(i-1)];

33 Numt(kk,i) = Sol(end,1)+Sol(end,2);

34 end

35 waitbar(kk/length(age))

36end

37

38close(wb)

39toc

40clear init_age % remove large vector to release memory

where the function $invivo\_sim$ appearing in the command of ode15s is given by

⬇

1function dy = invivo_sim(~,y)

2

3% y = [ring;troph;age;DHA concentration]

4

5kr=0.5*y(4)^2/(y(4)^2+200^2);

6kt=2*y(4)^2/(y(4)^2+200^2);

7

8alt=1; % drug-dependent ring-to-troph transition rate

9% set alt=0 for a drug-independent ring-to-troph transition rate

10

11if y(3)>=21

12 lambda = 22/(34+alt*y(4));

13else

14 lambda = 0;

15end

16

17dy = zeros(4,1);

18

19dy(1) = -kry(1)-lambday(1);

20dy(2) = lambday(1)-kty(2);

21dy(3) = 1;

22dy(4) = 0;

Supplementary figures

Supplementary tables

Bibliography22

The reference list from the paper itself. Each links out to its DOI / PubMed record.

11. World Health Organization. World malaria report 2015. World Health Organization, Geneva, Switzerland. 2015;.
22. Dondorp AM, Nosten F, Yi P, Das D, Phyo AP, Tarning J, et al. Artemisinin resistance in Plasmodium falciparum malaria. N Engl J Med. 2009;361:455–467.
33. Phyo AP, Nkhoma S, Stepniewska K, Ashley EA, Nair S, Mc Gready R, et al. Emergence of artemisinin-resistant malaria on the western border of Thailand: a longitudinal study. Lancet. 2012;379(9830):1960–1966.
44. Ariey F, Witkowski B, Amaratunga C, Beghain J, Langlois AC, Khim N, et al. A molecular marker of artemisinin-resistant Plasmodium falciparum malaria. Nature. 2014;505:50–55.
55. Ashley EA, Dhorda M, Fairhurst RM, Amaratunga C, Lim P, Suon S, et al. Spread of artemisinin resistance in Plasmodium falciparum malaria. N Engl J Med. 2014;371:411–423.
66. Simpson JA, Zaloumis S, De Livera AM, Price RN, Mc Caw JM. Making the most of clinical data: reviewing the role of pharmacokinetic-pharmacodynamic models of anti-malarial drugs. AAPS J. 2014;16(5):962–974.
77. Klonis N, Xie SC, Mc Caw JM, Crespo-Ortiz MP, Zaloumis SG, Simpson JA, et al. Altered temporal response of malaria parasites determines differential sensitivity to artemisinin. Proc Natl Acad Sci USA. 2013;110(3):5157–5162.
88. Dogovski C, Xie SC, Burgio G, Bridgford J, Mok S, Mc Caw JM, et al. Targeting the cell stress response of Plasmodium falciparum to overcome artemisinin resistance. P Lo S Biol. 2015;13(4):e 1002132.