Modeling the dynamics of Plasmodium falciparum gametocytes in humans during malaria infection

  1. Pengxing Cao  Is a corresponding author
  2. Katharine A Collins
  3. Sophie Zaloumis
  4. Thanaporn Wattanakul
  5. Joel Tarning
  6. Julie A Simpson
  7. James McCarthy
  8. James M McCaw  Is a corresponding author
  1. University of Melbourne, Australia
  2. Radboud University Medical Center, Netherlands
  3. QIMR Berghofer Medical Research Institute, Australia
  4. Mahidol University, Thailand
  5. University of Oxford, United Kingdom
  6. Peter Doherty Institute for Infection and Immunity, Australia

Abstract

Renewed efforts to eliminate malaria have highlighted the potential to interrupt human-to-mosquito transmission — a process mediated by gametocyte kinetics in human hosts. Here we study the in vivo dynamics of Plasmodium falciparum gametocytes by establishing a framework which incorporates improved measurements of parasitemia, a novel gametocyte dynamics model and model fitting using Bayesian hierarchical inference. We found that the model provides an excellent fit to the clinical data from 17 volunteers infected with P. falciparum (3D7 strain) and reliably predicts observed gametocytemia. We estimated the sexual commitment rate and gametocyte sequestration time to be 0.54% (95% credible interval: 0.30–1.00%) per asexual replication cycle and 8.39 (6.54–10.59) days respectively. We used the data-calibrated model to investigate human-to-mosquito transmissibility, providing a method to link within-human host infection kinetics to epidemiological-scale infection and transmission patterns.

https://doi.org/10.7554/eLife.49058.001

Introduction

Malaria is a mosquito-borne parasitic disease caused by protozoan parasites of the Plasmodium genus. It is estimated to have caused approximately 219 million new cases and 435,000 deaths in 2017, primarily due to Plasmodium falciparum (The World Health Organization, 2018). New tools will be required to achieve the ambitious goal of malaria elimination. Among the tools proposed are novel interventions to block transmission from human hosts to vector mosquitoes (The malERA Refresh Consultative Panel on Tools for Malaria Elimination, 2017). P. falciparum malaria is transmitted from humans to the mosquito when terminally differentiated male and female sexual-stages of the parasite, called gametocytes, are taken up by female Anopheles mosquito during a blood meal (Bousema and Drakeley, 2011; Josling and Llinás, 2015). The level of gametocytes in the blood, often referred to as gametocytemia, is highly associated with the probability of human-to-mosquito transmission (Bradley et al., 2018; Churcher et al., 2013). Gametocyte levels below a certain threshold (i.e., <1000 per mL blood) minimize the probability that a mosquito will take up both a male and female gametocyte during a blood-meal, which is necessary to propagate infection (Collins et al., 2018). An accurate understanding of the kinetics of gametocyte development in the human host is essential to predict the probability of transmission. A mathematical model that accurately captures the processes that give rise to observed gametocyte kinetics would be an important predictive tool to facilitate the design and evaluation of effective intervention strategies.

There is significant uncertainty surrounding fundamental aspects of P. falciparum gametocyte dynamics in humans. Parameters such as how many gametocytes are produced during each asexual replication cycle, the period of time in which early gametocytes disappear from the circulation before mature gametocytes appear (referred to as sequestration), and the period in which gametocytes circulate are poorly quantified. These gaps in understanding are due to a range of technical and logistic limitations. The first is the relatively poor sensitivity of the standard diagnostic test, namely microscopic examination of blood-films. Previous in vivo estimates of gametocyte kinetic parameters have been primarily based on historical data from neurosyphilis patients who were treated with so-called malariotherapy (Diebner et al., 2000; Eichner et al., 2001). In these studies, the limit of quantification was approximately 104 parasites/mL blood, at least two orders of magnitude higher than that of current quantitative PCR (qPCR) assays (Rockett et al., 2011). This high limit of quantification prevents an accurate estimation of onset of emergence of both asexual parasites and mature gametocytes in peripheral blood. The second limitation is that the available estimates of gametocyte dynamics parameters based on in vitro cultures (Filarsky et al., 2018; Gebru et al., 2017) may not be applicable to natural infection with P. falciparum gametocytes due to in vitro conditions that may not mimic the human host (Bousema and Drakeley, 2011).

Recent advances in experimental medicine using volunteer infection studies (VIS), otherwise known as controlled human malaria infection (CHMI) studies (Coffeng et al., 2017), allow prospective study design and data collection with the explicit aim of collecting in vivo data (McCarthy et al., 2011), in particular an improved quantification of P. falciparum gametocyte kinetics by qPCR applied in a novel VIS (Collins et al., 2018). Furthermore, the models and fitting methods used in the neurosyphilis patient studies have been superseded for parameter estimation by increasingly sophisticated within-host models (Khoury et al., 2018) and improvements in computational algorithms for Bayesian statistical inference (Piray et al., 2019). Therefore, there is an emerging opportunity to improve our quantitative understanding of the dynamics of P. falciparum gametocytes in human hosts by combining the novel VIS data and advanced modeling approaches.

In this study, we developed a novel mathematical model of gametocyte dynamics, fitted the model to the VIS data and estimated the gametocyte dynamics parameters using a Bayesian hierarchical inference method. We demonstrate that the data-calibrated model can reliably predict the time-course of gametocytemia and thus should form an essential part of modeling studies of malaria transmission.

Results

Model fitting and validation

The outcome variable used in model fitting was the total parasitemia (total circulating asexual parasites and gametocytes per mL blood measured using qPCR) collected from a previously published VIS (Collins et al., 2018), with other measurements from the same study, such as the asexual parasitemia (circulating asexual parasites per mL blood) and gametocytemia (circulating female and male gametocytes per mL blood), used to validate model predictions.

The results of fitting the mathematical model to total parasitemia data for all 17 volunteers are given in Figure 1 where 12 of 17 volunteers experienced recrudescence. The median of posterior predictions and 95% prediction interval (PI) were computed from 5000 model simulations based on 5000 samples from the posterior parameter distribution (see Materials and methods). The results show that the predicted total parasitemia (median and 95% PI) is able to accurately capture the trends of the data through the (visual) posterior predictive check. Furthermore, the narrow 95% PI indicates a low level of uncertainty in predicted total parasitemia.

Figure 1 with 15 supplements see all
Results of model fitting for all 17 volunteers.

Data are presented by circles. The median of posterior predictions (solid line) and 95% prediction interval (PI, shaded area) are generated by 5000 model simulations based on 5000 samples from the posterior parameter distribution as described in the Materials and methods. The histograms showing the posterior distributions of population mean and standard deviation hyperparameters are given in Figure 1—figure supplements 1 and 2. The posterior distribution of each model parameter (see the Materials and methods) for individual volunteers are given in Figure 1—figure supplements 314 Posterior distributions for some biological parameters are given in Figure 1—figure supplement 15, which are generated based on the posterior samples of population mean parameters (see the Materials and methods) and will be used to support the results in Table 1 shown later. The source data and computer code with instructions of implementation to generate Figure 1 and Figure 1—figure supplements 115 are fully publicly available at https://doi.org/10.26188/5cde4c26c8201.

https://doi.org/10.7554/eLife.49058.002
Table 1
Estimates of some key biological parameters and comparison with the literature.

The estimates of the biological parameters (middle column) are shown as the median and 95% credible interval (CI) of the marginal posterior parameter distribution (Figure 1—figure supplement 15). Estimates from the literature (third column) are shown in the format of either ‘mean estimate (95% confidence interval)’ or ‘mean estimate [minimum – maximum estimate]’ or simply ‘a low estimate – a high estimate’. Note some quoted estimates are from either in vivo or in vitro studies of P. falciparum. The source data and computer code with instructions of implementation to generate our model estimates (middle column) in Table 1 are fully publicly available at https://doi.org/10.26188/5cde4c26c8201.

https://doi.org/10.7554/eLife.49058.018
Biological parameters (unit)Median estimate (95% CI)Estimates in the literature
Sexual commitment rate (%/asexual replication cycle)0.54 (0.30–1.00)11 (6.2–15.8)
(Filarsky et al., 2018) (in vitro)
0.64 [0.027–13.5] (Eichner et al., 2001) (in vivo)
Gametocyte sequestration time (days)8.39 (6.54–10.59)7.4 [4 – 12] (Eichner et al., 2001) (in vivo)
Circulating gametocyte lifespan (days)63.5 (12.7–1513.9)16–32 (Gebru et al., 2017) (in vitro)
6.4 [1.3–22.2] (Eichner et al., 2001) (in vivo)
Parasite multiplication factor (per asexual replication cycle)21.8 (17.6–26.9)10–33 (Wockner et al., 2017) (in vivo)
16.4 (15.1–17.8)a (in vivo)
  1. a JS McCarthy, personal communication, May 2019.

Having calibrated the model against total parasitemia, the 5000 posterior parameter sets were used to calculate the median of posterior predictions and 95% PI of the asexual parasitemia and gametocytemia versus time profiles. Model predictions of the asexual parasitemia and gametocytemia for all 17 volunteers are shown in Figure 2 and Figure 3 respectively (curves: median prediction; shaded areas: 95% PI) and are compared to the asexual parasitemia and gametocytemia data (circles). We emphasize that this was not a fitting exercise, rather an independent validation of the calibrated model.

Comparison of model predictions and clinical data for the asexual parasitemia for all 17 volunteers.

Data are presented by circles. The median of posterior predictions (solid curve) and 95% PI (shaded area) are generated by 5000 model simulations based on 5000 samples from the posterior parameter distribution as described in the Materials and methods. The data points with one parasite/mL (i.e., those points which lie on the dotted line) indicate measurements for which no parasites were detected. No data are available for Volunteer 101 and 106 to validate the model predictions. The source data and computer code with instructions of implementation to generate Figure 2 are fully publicly available at https://doi.org/10.26188/5cde4c26c8201.

https://doi.org/10.7554/eLife.49058.019
Comparison of model predictions and clinical data for the gametocytemia for all 17 volunteers.

Data are presented by circles. The median of posterior predictions (solid curve) and 95% PI (shaded area) are generated by 5000 model simulations based on 5000 samples from the posterior parameter distribution as described in the Materials and methods. The data points with one parasite/mL (i.e. those points which lie on the dotted line) indicate measurements for which no parasites were detected. The source data and computer code with instructions of implementation to generate Figure 3 are fully publicly available at https://doi.org/10.26188/5cde4c26c8201.

https://doi.org/10.7554/eLife.49058.020

For the majority of asexual parasitemia data the model predictions (median and 95% PI) can faithfully capture the trends of the data (Figure 2), in particular the accurate predictions for both the recrudescent case where a portion of asexual parasitemia data diverge from the total parasitemia measurement (e.g., Volunteer 103, 104, 105, 201, 203, 304 and 307) and the non-recrudescent case where the posterior-median prediction curve (solid red curve) lies below the limit of detection (one asexual parasite/mL) (e.g., Volunteer 202, 301 and 302). However, there are some discrepant observations. The model under-predicts (Volunteer 204) or over-predicts (Volunteer 303, 305 and 306) a portion of the asexual parasitemia data. Furthermore, for some volunteers such as 202, 301 and 302, the 95% PI (red shaded area) extends into the detectable region again after the asexual parasitemia reaches below the detection limit, indicating that there was a small chance that the patients may have suffered a recrudescence during the observation period (of course, they did not) or after the observation period (although this predication cannot be evaluated because artemisinin combination therapy was given immediately after the period).

Figure 3 shows the data and model predictions for the gametocytemia. Despite some discrepant observations for asexual parasitemia in Figure 2, we found that the model predictions of gametocytemia were able to capture the trends and levels of the gametocytemia data for all 17 volunteers.

Estimation of gametocyte dynamics parameters

The model calibration process provided the joint posterior density for the model parameters, which were used to estimate some key biological parameters governing the dynamics of P. falciparum gametocytes (detailed in the Materials and methods). As shown in Table 1, the sexual commitment rate — the percentage of asexual parasites entering sexual development during each asexual replication cycle — is estimated to be 0.54%/asexual replication cycle (95% credible interval (CI): 0.30–1.00%). This in vivo estimate of 0.54%/asexual replication cycle is much lower than 11%/asexual replication cycle that was estimated from in vitro data (Filarsky et al., 2018). The proportion of committed asexual parasites that survive to become mature gametocytes, calculated by discounting the sexual commitment rate by the probability of survival from the immature (stages I–IV) to mature (stage V) gametocyte life-stage, is 0.52%/asexual replication cycle (95% CI: 0.28–0.95%). The gametocyte sequestration time is the average time that immature gametocytes (stages I–IV) cannot be observed in the peripheral circulation. They re-emerge in the peripheral circulation as mature gametocytes (stage V). It was estimated to be 8.39 days (95% CI: 6.54–10.59 days). The estimate for the circulating gametocyte lifespan is 63.5 days, with a broad 95% CI (12.7–1513.9 days) resulting from the long-tailed posterior distribution (Figure 1—figure supplement 15) and is much longer than the previous in vitro estimate of 16–32 days (Gebru et al., 2017) (note that our lower bound of the 95% CI is lower than the in vitro estimated range). The wide estimate for the circulating gametocyte lifespan, and in particular the high upper bound of the 95% CI, is due to the limited observation time in the VIS which does not enable the lifespan to be accurately determined (explored in more detail in the Discussion).

As shown in Table 1, there are similarities in parameter estimates for P. falciparum between our analysis and the analysis of historical neurosyphilis patient data (Eichner et al., 2001). We found that they exhibited similar in vivo sexual commitment rate (median 0.54%/asexual replication cycle vs. mean 0.64%/asexual replication cycle with overlapping plausible ranges) and gametocyte sequestration time (median 8.39 days vs. mean 7.4 days with overlapping plausible ranges).

Finally, we provided an estimate for the parasite multiplication factor which is the average number of infected red blood cells produced by a single infected red blood cells after one replication cycle. The parasite multiplication factor is an important parameter that quantifies the net growth of asexual parasites and thus influences the rate of gametocyte generation. Our posterior-median estimate is 21.8 parasites per asexual replication cycle (95% CI: 17.6–26.9), consistent with previous estimates which lie in the range 10–33 (Wockner et al., 2017), and a bit larger than an updated estimate calculated from a pooled analysis of parasite counts from 177 volunteers infected with the same P. falciparum strain using a statistical model (16.4 parasites per asexual replication cycle) (JS McCarthy, personal communication, May 2019).

Predicting the impact of gametocyte kinetics on human-to-mosquito transmissibility

Having validated our mathematical model of asexual parasitemia and gametocyte dynamics, we were able to predict how the gametocyte dynamics parameters influence the transmissibility of P. falciparum malaria from humans to mosquitoes in various epidemiological scenarios. In particular, we focused on the early phase of infection where the innate immune response is minimal and treatment has not been administered (in order to avoid complications that our mathematical model was not designed to capture). Two scenarios were considered:

  • Predicting the potential infectiousness of newly hospitalized clinical malaria cases for various values of sexual commitment rate and gametocyte sequestration time. In the model, gametocytemia was assumed to be a surrogate of the potential infectiousness. We further assumed that patients would seek hospital admission when their total parasitemia reached approximately 108 parasites/mL. This choice was based on the microscopic measurements of the total parasitemia (i.e., asexual and sexual parasites) from a study of Cambodia and Thailand hospitalized malaria patients (Saralamba et al., 2011). As illustrated in Figure 4A, we simulated the model (for different sexual commitment rates and gametocyte sequestration times) and looked at the critical gametocytemia level (indicated by Gc) corresponding to the time when the total parasitemia (wave-like black curves) first reached 108 parasites/mL (their associations are indicated by the dotted lines and arrows).

  • Predicting the non-infectious period of malaria patients for various values of sexual commitment rate and gametocyte sequestration time. In the model, the non-infectious period was defined to be time from the inoculation of infected red blood cells to the time when the gametocytemia reached 103 parasites/mL, which is a threshold below which human-to-mosquito transmission was not observed (Collins et al., 2018). Note that this non-infectious period does not include the latent period due to the liver stage which should be considered if the starting time were to be taken from time of mosquito bite. As illustrated in Figure 4B, we simulated the model (for different sexual commitment rates and gametocyte sequestration times) and identified the critical time (indicated by tc) when the gametocytemia (blue curve) first reached 103 parasites/mL (their associations are indicated by the dotted lines and arrows).

Simulation of two scenarios predicting the dependence of human-to-mosquito transmissibility on the sexual commitment rate and gametocyte sequestration time.

(A) illustration of the first scenario: predicting the critical gametocytemia level (indicated by Gc) at the time when the total parasitemia reaches 108 parasites/mL. (B) illustration of the second scenario: predicting the non-infectious period (indicated by tc), which is defined to be time from inoculation of infected red blood cells to the time when the gametocytemia reaches 103 parasites/mL (a threshold below which human-to-mosquito transmission was not observed [Collins et al., 2018]). (C and D) Heatmaps showing the dependence of the critical gametocytemia Gc and the non-infectious period tc on the sexual commitment rate and gametocyte sequestration time. The black dots represent the value obtained by simulating the gametocyte dynamics model using the median estimates of the posterior samples of the population mean parameters as described in the Materials and methods. The red curve in C is the level curve for Gc = 103 parasites/mL. The red curve in D is the level curve for tc = 13.42 days which is the non-infectious period obtained by model simulation using the posterior estimates of the population mean parameters. The source data and computer code with instruction of implementation to generate Figure 4 are fully publicly available at https://doi.org/10.26188/5cde4c26c8201.

https://doi.org/10.7554/eLife.49058.021

A higher sexual commitment rate or a lower gametocyte sequestration time leads to a higher gametocytemia (Gc) at the time of hospitalization (Figure 4C). The red curve in Figure 4C indicates the level curve of 103 gametocytes/mL (i.e., the threshold for infectiousness as mentioned above) dividing the heatmap into two regions. To the left, Gc is below 103 gametocytes/mL, suggesting clinical presentation precedes infectiousness, while to the right Gc is above 103 gametocytes/mL and the converse applies. The Gc value obtained by model simulation using the median estimates of the population mean parameters (indicated by the black dot) is below 103 gametocytes/mL, suggesting that newly hospitalized malaria patients are less likely to be infectious, and thus efforts to identify and treat infections in a timely manner may have a substantial impact in terms of reduced transmission potential. Note that patients from clinical observations of uncomplicated malaria in endemic settings may have higher gametocyte counts at the time of presentation than what our model predicts. For example one study from the TRACII clinical trial reported a range of 16–5120 gametocytes/µL, which is much higher than our prediction of below one gametocytes/µL (or 103 gametocytes/mL) (van der Pluijm et al., 2019). One plausible explanation for the difference is that our model predicted a very fast rise in total parasitemia to 108 parasites/mL while the rise in parasitemia among patients in endemic settings may be slower due to the effect of immunity on the parasite multiplication. Immunity was not considered in our model due to the design of our VIS where only malaria naïve volunteers were recruited.

Figure 4D reinforces the result in Figure 4C using the non-infectious period (tc). As the sexual commitment rate increases or the gametocyte sequestration time decreases, tc decreases. However, for large values of the sexual commitment rate (e.g., >20%), we observed an increase in tc as the sexual commitment rate increases (see the top-right corner of ). This is because an increased sexual commitment rate leads to both a decrease in the rate of asexual parasite growth (due to a direct loss of asexual parasites as they convert to gametocytes) and an increase in the number of sexually committed parasites. For a very high sexual commitment rate, the impact of the former more than counterbalances that of the latter.

Discussion

We have developed a novel mathematical model of gametocyte dynamics that combines an existing multi-state asexual cycle model with a new model for the development of gametocytes. Model parameters were estimated by fitting the model to data from 17 malaria-naïve volunteers inoculated with P. falciparum-infected red blood cells (3D7 strain). Compared to previous studies, our work is distinguished by three novel contributions: (1) the use of a prospectively planned clinical trial to collect more accurate quantitative data of parasite levels measured by qPCR; (2) the development of a novel dynamics mathematical model which allows for robust and biologically-informed extrapolation and hypothesis testing/scenario analysis; and (3) the use of a Bayesian hierarchical inference method for model calibration and parameter estimation.

For gametocyte kinetic parameters, we found that our in vivo estimate of the P. falciparum sexual commitment rate was similar to that found in the neurosyphilis patient data (Eichner et al., 2001) but was much smaller than previous in vitro estimates (Table 1). Importantly, our estimate follows directly from the structure of our mathematical model, and accounts for the fact that some early committed gametocytes may not complete development and thus not emerge in peripheral circulation as mature gametocytes. Novel VIS data using biomarkers specific to early sexual parasites (e.g. AP2-G [Bancells et al., 2019] and PfGEXP5 [Tibúrcio et al., 2015]) would enable a direct (statistical) estimate of the sexual commitment rate, providing an independent validation of our gametocyte dynamics model. Our in vivo estimate for the circulating gametocyte lifespan is imprecise (i.e., has a very wide credible interval) due to the lack of available data for gametocyte clearance (treatment was initiated before gametocyte were naturally cleared in the VIS study). P. falciparum data with gametocytemia measurements over a longer period of time to capture the natural decay of circulating gametocytes, would greatly improve these estimates.

We also predicted the effects of altered gametocyte kinetic parameters on the transmissibility from humans to mosquitoes, focusing on two scenarios: the infectiousness of newly hospitalized clinical malaria cases (i.e., the gametocytemia when total parasitemia first reaches a level typically seen upon hospitalization — 108 parasites/mL in the model); and the non-infectious period of malaria patients (i.e., the time from the inoculation of infected red blood cells to the time when the gametocytemia reaches a minimal transmission threshold of 103 parasites/mL in the model). We explored how the sexual commitment rate and gametocyte sequestration time influenced the gametocyte level and the non-infectious period. We would like to emphasize that human-to-mosquito transmissibility is determined by both the level of gametocytemia and the relationship between gametocytemia and the probability of transmission per bite. A reliable prediction of the former is essential but not a sole determinant of transmissibility. Therefore, it is also important to improve our quantitative understanding of the probability of transmission per bite, which may be complicated by and also influenced by the densities and ratios of female and male gametocytes (Bradley et al., 2018; Churcher et al., 2013; Da et al., 2015).

Our study has some limitations. The gametocyte dynamics model, that has been shown to have sufficient complexity to reproduce the clinical observations, is still a rather coarse simplification of the actual biological processes. For example, the model does not assume an adaptive sexual commitment rate (Schneider et al., 2018), nor does it consider the mechanisms of sexual commitment (Bancells et al., 2019). Furthermore, the model assumes a constant gametocyte death rate but does not consider other non-constant formulations as have been previously proposed (Diebner et al., 2000). Another limitation is that we assumed a fixed duration for the asexual replication cycle of 42 hr, while previous work by our group suggests that the replication cycle may be altered by up to a few hours in response to antimalarial drugs (e.g., artemisinin [Cao et al., 2017; Dogovski et al., 2015]), though there is no evidence that piperaquine (which was administered in this VIS) has a similar effect.

In conclusion, we have developed a novel mathematical model of gametocyte dynamics, and demonstrated that it reliably predicts time series data of gametocytemia. The model provides a powerful predictive tool for informing the design of future volunteer infection studies which aim to test transmission-blocking interventions. Furthermore, the within human host transmission model can be incorporated into epidemiological-scale models to refine predictions of the impacts of various antimalarial treatments and transmission interventions.

Materials and methods

Study population and measurements

Request a detailed protocol

The data used in this modeling study are from a previously published VIS (Collins et al., 2018) where 17 malaria-naïve volunteers were inoculated with approximately 2800 P. falciparum-infected red blood cells (3D7 strain). The study was approved by the QIMR Berghofer Human Research Ethics Committee and registered with ClinicalTrials.gov (NCT02431637 and NCT02431650). The volunteers were treated with 480 mg piperaquine phosphate (PQP) on day 7 or 8 post-inoculation to attenuate asexual parasite growth and a second dose of 960 mg PQP was given to any volunteer for treatment of recrudescent asexual parasitemia. All volunteers received a course of artemether/lumefantrine and, if required, a single dose of primaquine (45 mg) to clear all parasites. Parasitemia in the volunteers was monitored approximately daily following inoculation, but with notable variability in the frequency of data collection at later times as described by Collins et al. (2018).

Molecular analysis of parasite levels was carried out throughout the study. The total parasitemia was measured by 18S qPCR (total circulating asexual parasites and gametocytes per mL blood), asexual parasitemia was measured by SBP1 qRT-PCR (circulating asexual parasites per mL blood), and gametocytemia was measured by Pfs25 and PfMGET qRT-PCR (circulating female and male gametocytes per mL blood). Plasma concentrations of PQP were also determined at multiple time points after inoculation. Further details about the VIS are given in Collins et al. (2018). It is important to note that the data used in model fitting is the total parasitemia (from the first measurement to the time before any treatment other than PQP) and the other data, that is asexual parasitemia and gametocytemia (also up to the time of treatment other than PQP) are used to validate the model.

Gametocyte dynamics model

Request a detailed protocol

The mathematical model extends the published models of asexual parasite replication cycle (Saralamba et al., 2011; Zaloumis et al., 2012) by incorporating the development of gametocytes. The model is comprised of three parts describing three populations of parasites: asexual parasites (P), sexually committed parasites (PG) and gametocytes (G). A schematic diagram of the development of those populations based on current knowledge (Bancells et al., 2019; Filarsky et al., 2018) is shown in Figure 5.

Schematic diagram showing the model compartments and transitions.

The model is comprised of three parts describing three populations of parasites: asexual parasites (Pa,t), sexually committed parasites (PGa,t) and gametocytes (Gt). P and PG are functions of asexual parasite age a and time t. Square compartments in the inner loop represent the asexual parasite population which follows a cycle of maturation and replication every aL hours. Sexual commitment occurs from age as and a fraction of asexual parasites become sexually committed (the bigger square compartments in the outer loop) and eventually enter the development of stage I–V gametocytes (G1–G5). The compartments with a dashed boundary are sequestered to tissues and thus not measurable in a blood smear. The notation for each compartment is consistent with those in the model equations and is explained in the main text.

https://doi.org/10.7554/eLife.49058.022

Asexual parasites develop and replicate in the red blood cells (RBCs) until cell rupture at the end of each replication cycle and the released free parasites (merozoites) can initiate new cycles of replication if they successfully invade susceptible RBCs. At the time of inoculation (i.e., t=0 hours in the model), we define the inoculum size to be Pinit and assume the age distribution of inoculated parasites is Gaussian with mean μ and standard deviation σ. As time increments by one hour, the asexual parasites of age a at time t (denoted as P(a,t)) follow the iterative equation:

(1) P(a,t)={P(a1,t1)ekd¯δP,a=2,3,...,aLrPP(aL,t1)ekd¯δP,a=1

where kd¯ represents the average rate of asexual parasite killing by PQP and δP is the rate of asexual parasite death due to processes other than PQP. kd¯ is approximated by the average of kd(t-1) and kd(t) and kdt= kmaxCtγ/(Ctγ+EC50γ) where kmax is the maximum killing rate, EC50 is the PQP concentration at which half maximum killing is achieved, and γ is the Hill coefficient determining the curvature of the dose-response curve. Ct is the PQP concentration which is simulated by a pharmacokinetic model introduced below. aL is the length of each asexual replication cycle and rP is the parasite replication rate indicating the average number of newly infected RBCs attributable to the rupture of a single infected RBC. Note that we distinguish the parasite replication rate rP from the so-called parasite multiplication factor, the latter of which is a 'net replication rate' quantified by the (per cycle) increase in parasite numbers due to replication (rP) and the decrease in parasite numbers due to death or sexual commitment. Sexual commitment is assumed to occur at the first age of the trophozoite stage (denoted to be as) and a fraction (f) of asexual parasites leave the asexual replication cycle and start sexual development in the next hour, which is modeled by

(2) P(as+1,t)=(1f)P(as,t1)ekd¯δP
(3) PG(as+1,t)=fP(as,t1)ekd¯δP.

The first equation describes the proportion of parasites remaining in the asexual replication cycle while the second equation describes the proportion of parasites becoming sexually committed parasites (PG). According to Figure 5, the sexually committed parasites continue the rest of the replication cycle and a part of the next replication cycle (note that they appear indistinguishable from asexual parasites by microscopy) before becoming stage I gametocytes. The process is modeled by

(4) P(a,t)={P(a1,t1)eδP,a=2,3,...,aLexcepta=asanda=as+1rPP(aL,t1)eδP,a=1

Note that we assumed in our model that PQP does not kill gametocytes. Our assumption was based on evidence from both in vitro and in vivo experiments that suggests that PQP has little activity against sexually committed parasites and gametocytes (Collins et al., 2018; Pasay et al., 2016; Bolscher et al., 2015), although we note there is some evidence that PQP might have activity against early-stage I/II gametocytes (Adjalley et al., 2011). The changes of the sequestered stage I–IV gametocytes (G1G4) are governed by difference equations

(5) G1(t)=G1(t1)e(m+δG)+PG(as1,t1)eδP(1e(m+δG))m+δG,
(6) G2(t)=G2(t1)e(m+δG)+mG1(t1)(1e(m+δG))m+δG,
(7) G3(t)=G3(t1)e(m+δG)+mG2(t1)(1e(m+δG))m+δG,
(8) G4(t)=G4(t1)e(m+δG)+mG3(t1)(1e(m+δG))m+δG,

where m is the rate of gametocyte maturation and δG is the death rate of sequestered gametocytes. Stage V gametocytes are circulating in bloodstream (and therefore can be measured from the peripheral blood film) modeled by

(9) G5(t)=G5(t1)eδGm+mG4(t1)(1eδGm)δGm,

where δGm is the death rate of mature circulating gametocytes.

The total parasitemia in the model is given by a=1a=as-1Pa,t+PGa,t+G5(t), which was fitted to the VIS data. After model fitting, we simulated the asexual parasitemia a=1a=as-1P(a,t) and gametocytemia G5(t) and compared them with associated data for model validation. Table 2 presents all the model parameters and their units and descriptions.

Table 2
Details of the gametocyte dynamics model parameters.

The table includes the unit, description and prior distribution for each model parameter. For the uniform prior distributions (U), the lower bounds are non-negative based on the definitions of the model parameters and the upper bounds for the prior distributions were chosen to be sufficiently wide in order to accommodate all biologically plausible values from the literature (Zaloumis et al., 2012). We assumed parasites younger than 25h are circulating and thus fix as to be 25h. For 3D7 strain, the asexual replication cycle is approximately 39–45h (based on in vitro estimates [Duffy and Avery, 2017] and personal communication [JS McCarthy, personal communication, May 2019]) and we fix aL to be 42h.

https://doi.org/10.7554/eLife.49058.023
ParameterUnitDescriptionPrior distribution
Pinitparasites/mLinoculation sizeU(0, 10)
μhmean of the initial parasite age distributionU(0, 35)
σhSD of the initial parasite age distributionU(0, 20)
rP(unitless)parasite replication rateU(0, 100)
kmaxh−1maximum rate of parasite killing by PQPU(0, 1)
EC50ng/mLhalf-maximum effective PQP concentrationU(1, 100)
γ(unitless)Hill coefficient for PQPU(0, 20)
f(unitless)the fraction of parasites entering sexual development per asexual replication cycleU(0, 1)
δPh−1death rate of asexual and sexual parasitesU(0, 0.2)
mh−1maturation rate of gametocytesU(0, 0.1)
δGh−1death rate of sequestered gametocytesU(0, 0.1)
δGmh−1death rate of circulating gametocytesU(0, 0.1)
ashsequestration age of asexual parasitesfixed to be 25
aLhlength of life cycle of asexual parasitesfixed to be 42
  1. SD: standard deviation; h: hour.

Pharmacokinetic model of piperaquine (PQP)

Request a detailed protocol

In the within-host model, the killing rate kd(t) is determined by PQP concentration Ct which was simulated from a pharmacokinetic (PK) model introduced in this section. The PK model, provided by Thanaporn Wattanakul and Joel Tarning (Mahidol-Oxford Tropical Medicine Research Unit, Bangkok), is a three-compartment disposition model with two transit compartments for absorption (see the schematic diagram in Figure 6).

Figure 6 with 1 supplement see all
The pharmacokinetic model of piperaquine (PQP).

The model is a three-compartment disposition model with two transit compartments for absorption. State D represents the dose of PQP. T1 and T2 represent the two transit compartments. C is the central compartment and PQP concentration in this compartment was measured (which are shown in Figure 6—figure supplement 1). P1 and P2 represent two peripheral compartments. kT, q1, q2 and qc are the rates of flow into or out of compartments.

https://doi.org/10.7554/eLife.49058.024

Based on Figure 6, the model is formulated to be a system of ordinary differential equations:

(10) dDdt=kTD,
(11) dT1dt=kTDkTT1,
(12) dT2dt=kTT1kTT2,
(13) dCdt=kTT2+q1P1+q2P2q1Cq2CqcCVc ,
(14) dP1dt=q1Cq1P1V1 ,
(15) dP2dt=q2Cq2P2V2 ,

where kT and q’s are rate constants as shown in Figure 6 and Vc, V1 and V2 are the volume of distribution for the central compartment (in which PQP concentration is C), peripheral compartment 1 (in which PQP concentration is P1) and peripheral compartment 2 (in which PQP concentration is P2) respectively.

Under the sequential pharmacokinetic-pharmacodynamic (PK-PD) approach we have taken, a PQP concentration curve (C(t)) for each volunteer is a required input into the gametocyte dynamics model. The VIS, with its limited sampling of PQP for each volunteer, was not designed to provide this PQP concentration curve directly, so we used a PK model, informed by data from a previous VIS with rich sampling. We drew on an analysis of that previous VIS by Thanaporn Wattanakul and Joel Tarning (unpublished data and estimates). Their analysis provides population-level PQP PK model parameter estimates.

We used MATLAB’s (version 2016b; The MathWorks, Natick, MA) built-in least-squares optimizer lsqcurvefit (with the default setting) to optimize the PK curve for each volunteer in the VIS study. We applied the optimizer to each volunteer’s (limited) PQP data, using the parameter estimates provided by Thanaporn Wattanakul and Joel Tarning as initial values. We applied some further model parameter constraints as specified in Appendix 1. This approach provided us with a data-informed PK curve for each volunteer in the VIS, sufficient for our primary purpose of studying the asexual and sexual parasite dynamics. Of note, Volunteers 202, 301, 302 or 307 had fewer PK data points than PK model parameters, preventing application of this optimization procedure. For these volunteers, their predicted PQP PK curve was derived using the population-level mean PK parameter from Wattanakul and Tarning’s analysis. The MATLAB code (with detailed comments) is publicly available at https://doi.org/10.26188/5cde4c26c8201 The details of the initial conditions, starting point and constraints for the PK curve optimization procedure are provided in Appendix 1. The optimized PK curves and associated parameter values for all volunteers are provided in Figure 6—figure supplement 1 and Appendix 1.

Fitting the model to parasitemia data

Request a detailed protocol

We took a Bayesian hierarchical modeling approach (Gelman et al., 2013) to fit the gametocyte dynamics model to the data from all 17 volunteers. In detail, each volunteer has 12 model parameters (i.e., those in Table 2 except as and aL; also called the individual parameters) and lower and upper bounds of the parameters are given in Table 2. If denoting the individual parameters to be θind and lower and upper bounds to be bL and bU respectively, the following transformations are used to convert the bounded individual parameters to unbounded ones (denoted by φind) in order to in order to improve computational efficiency (Lesaffre et al., 2007; Stan Development Team, 2017):

(16) φind=ln(θindbLbUθind),

oφind beys a multivariate normal distribution 𝒩(φpop, Ωpop) where

(17) φpop=ln(θpopbLbUθpop)

and θpop is a vector containing 12 population mean parameters (hyperparameters) corresponding to the 12 gametocyte dynamics model parameters. Ωpop is the covariance matrix. For more efficient sampling process, φind~𝒩(φpop, Ωpop) was reparameterised to a non-centerd form φind=φpop+ωpopLη, where ωpop is the diagonal standard deviation (SD) matrix whose diagonal elements are the 12 population SD parameters (hyperparameters); L is the lower Cholesky factor of the correlation matrix; η obeys standard multivariate normal distribution. Note that Ωpop= ωpopLLTωpop where LLT is the correlation matrix. The prior distributions for the 12 population mean parameters θpop are given by uniform distributions with bounds given in Table 2. The prior distribution for the 12 population SD parameters is standard half-normal and the prior distribution for the lower Cholesky factor of the correlation matrix L is given by Cholesky LKJ correlation distribution with shape parameter of 2 (Lewandowski et al., 2009; Stan Development Team, 2017). The distribution of the observed parasitemia measurements is assumed to be a log normal distribution with mean given by the model-simulated values and SD parameter with prior distribution of a half-Cauchy distribution with a location parameter of zero and a scale parameter of 5. The distribution for the observed parasitemia measurements was used to calculate the likelihood function and the M3 method (Ahn et al., 2008) was used to penalise the likelihood for data points below the limit of detection for the total parasitemia (10 parasites/mL; Collins et al., 2018).

Model fitting was implemented in R (version 3.2.3) (R Development Core Team, 2017) and Stan (RStan 2.17.3) (Stan Development Team, 2017) using the Hamiltonian Monte Carlo (HMC) optimized by the No-U-Turn Sampler (NUTS) to draw samples from the joint posterior distribution of the parameters including the individual parameters (12 parameters for each volunteer) and population mean parameters (12 hyperparameters). Five chains with different starting points (set by different random seeds) were implemented and 1000 posterior samples retained from each chain after a burn-in of 1000 iterations (in total 5000 samples were drawn from the joint posterior distribution). The marginal posterior and prior distributions of the population mean and SD parameters are shown in Figure 1—figure supplements 1 and 2. The marginal posterior distributions of the individual parameters for all 17 volunteers are shown in Figure 1—figure supplements 314 (using violin plots). For each volunteer, the 5000 sets of individual parameters are used to simulate the gametocyte dynamics model and generate 5000 simulated model outputs (e.g., 5000 time series of total parasitemia, asexual parasitemia or gametocytemia). The posterior prediction and 95% prediction interval (PI) are given by the median and quantiles of 2.5% and 97.5% of the 5000 model outputs at each time respectively (see Figures 13 for example).

The estimates of some key biological parameters (Table 1) were calculated using the 5000 posterior draws of the 12 population mean parameters, that is median and 2.5%- and 97.5%-quantile (95% credible interval). The sexual commitment rate was calculated by fpop×100% (fpop is the population mean parameter for f) and the proportion of committed asexual parasites that survive to become mature gametocytes was calculated by fpop(mpop/(mpop+δGpop))4 where the factor of four arises due to the four sequestered gametocyte stages (I to IV). Circulating gametocyte lifespan was calculated by 1/δGmpop/24  (the factor of 24 converts hours into days). Gametocyte sequestration time was calculated by 4/mpop/24 where 4 indicates four sequestered state (stage I to IV) and 24 converts hours into days. Parasite multiplication factor is calculated by rPpopexp(-δPpopaL)(1-fpop) where the term exp(-δPpopaL)(1-fpop) gives the fraction of surviving asexual parasites after death and sexual conversion per replication cycle.

The gametocyte dynamics model with parameters given by the median estimates of the population mean parameters was used to simulate the two scenarios predicting the dependence of human-to-mosquito transmissibility on the sexual commitment rate and gametocyte sequestration time (Figure 4).

Final analysis and visualization were performed in MATLAB. All computer codes (R, Stan, MATLAB), data and fitting results (CSV and MAT files) and an instruction document (note that reading the document first will make the code much easy to follow) are publicly available at https://doi.org/10.26188/5cde4c26c8201.

Appendix 1

The PK model (see the Materials and methods) describes absorption of administrated drug mass in compartment D and the subsequent kinetics of drug concentration in the central compartment C. The initial conditions for the model simulation were that all the compartments (i.e., D, T1, T2, C, P1 and P2) were zero because PQP was not given until day 7 or 8. When the first dose of 480 mg PQP was given, we set D = 480 and all other compartments to zero. When a second dose of 960 mg was required for some volunteers to treat recrudescence, we again set D to be 960 plus any remaining level of D from the first dose for those volunteers.

The PK data and the PK curves generated using our optimization approach are shown in Figure 6—figure supplement 1. The starting values for optimization are kT=1/2.79, qc=51.4, Vc=804, q1=2020, V1=3010, q2=149 and V2=13300, which are the population-level mean estimates provided by Thanaporn Wattanakul and Joel Tarning (Mahidol-Oxford Tropical Medicine Research Unit, Bangkok). We also assume the parameters are constrained in the optimization procedure by lower bounds of kT=1/3.25, qc=40.6, Vc=469, q1=746, V1=2344, q2=117 and V2=10350 and upper bounds of kT=1/2.32, qc=160, Vc=1342, q1=5034, V1=3986, q2=190 and V2=17546. Note that kT is expressed by the reciprocal of the mean transition time. The lower and upper bounds are the limits of the 95% confidence intervals of the parameter estimate distributions provided by Thanaporn Wattanakul and Joel Tarning except that the upper bound for the clearance rate qc was increased from 62.6 to 160 such that the PK curves could better capture the fast PQP concentration decay observed for some volunteers. The constraints are necessary due to the limited available PK data which prevents identification of the PK parameters. For Volunteer 202, 301, 302, 307 whose numbers of PK data points are less than the number of parameters in the PK model (such that optimization fails), their 'best-fit' parameter values are the starting values given above with some adjustments on qc (e.g., 71.4 for 202, 61.4 for 301 and 81.4 for 302) which are required to allow the simulated curves to visually capture the data (as shown in Figure 6—figure supplement 1 where the predicted PQP concentrations are very close to the observed concentrations).

The PK parameter values obtained from our optimization approach to produce the reasonable PK curves for all 17 volunteers are provided in the Appendix 1—table 1 below.

Appendix 1—table 1
Parameter values used to generate the optimized PK curves for all 17 volunteers.

The units of the model parameters are given in the parentheses. The optimized PK curves are shown in Figure 6—figure supplement 1.

https://doi.org/10.7554/eLife.49058.028
PK parameter (unit)kT(h-1)qc(L/h)Vc(L)q1(L/h)V1(L)q2(L/h)v2(L)
Volunteer 1011/2.4353.78312240395811817177
Volunteer 1021/2.32123.4582849398516310355
Volunteer 1031/3.0761.37562686391012714346
Volunteer 1041/2.32134.7488746398619015007
Volunteer 1051/3.2516013425034398619017546
Volunteer 1061/2.4060.76672992304312316374
Volunteer 2011/2.32160469746398619015557
Volunteer 2021/2.7971.48042020301014913300
Volunteer 2031/2.32141.91243803397319011825
Volunteer 2041/2.3263.38952169332716516205
Volunteer 3011/2.7961.48042020301014913300
Volunteer 3021/2.7981.48042020301014913300
Volunteer 3031/2.45159.41332779397019017355
Volunteer 3041/2.48158.91339762398119017290
Volunteer 3051/2.53120.911521905297912313896
Volunteer 3061/2.32158.8870761398319017488
Volunteer 3071/2.7951.48042020301014913300
https://doi.org/10.7554/eLife.49058.027

Data availability

All data generated or analysed during this study are included in the manuscript, the supporting information and University of Melbourne's repository (publicly available at https://doi.org/10.26188/5cde4c26c8201).

The following data sets were generated
    1. Pengxing Cao
    (2019) Melbourne figshare
    Data and computer code for model fitting and simulation.
    https://doi.org/10.26188/5cde4c26c8201

References

  1. Book
    1. Gelman A
    2. Carlin JBB
    3. Stern HSS
    4. Dunson DB
    5. Vehtari A
    6. Rubin DBB
    (2013)
    Bayesian Data Analysis (Third Edition)
    CRC Press.
    1. Stan Development Team
    (2017)
    Stan Modeling Language Users Guide and Reference Manual
    Stan Modeling Language Users Guide and Reference Manual.

Decision letter

  1. Isabel Rodriguez-Barraquer
    Reviewing Editor; University of California, San Francisco, United States
  2. Neil M Ferguson
    Senior Editor; Imperial College London, United Kingdom

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Modelling the dynamics of Plasmodium falciparum gametocytes in humans during malaria infection" for consideration by eLife. Your article has been reviewed by Neil Ferguson as the Senior Editor, a Reviewing Editor, and two reviewers. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

In this work, the authors have extended a within-host model of a P. falciparum infection to include gametocyte production. The model is calibrated to parasitaemia from 17 volunteers in a controlled human infection study. Although well thought out and well written, a number of concerns remain: (1) several clarifications need to be made in the description of the models; (2) additional justification is needed for the statistical validation of the results; (3) the section on predicting the impact of gametocyte kinetics seems poorly supported.

Essential revisions:

1) Model description

- a and t in the model are never defined (subsection “Study population and measurements”).

- Γ is never defined (subsection “Gametocyte dynamics model”).

- Where do the parameters ranges from Table 2 originate? Why are the death rates so high? For example, the maximum death rate of gametocytes is 0.1/hr, suggesting an average life span of 10 hours.

- Although never stated in this manuscript, it seems that the inoculum from the Collins paper was 2800 parasites. Discuss how the range of 0-10 parasites/mL relates to this.

- Vc, V1 and V2 are never defined (subsection “Pharmacokinetic model of piperaquine (PQP)”).

- How are the initial conditions of the drug model chosen? Why are there two spikes for some patients and only one for others (Figure 6—figure supplement 1)?

- In the fitting of the PK model, was only a single starting set considered? Why is it necessary and valid to increase the upper bound for qc? For different parameters, several patients hit the upper or lower bounds. In particular, Volunteers 105 and 20 hit the bounds on almost all parameters. Could their optimal fit be outside the range?

- Justification for all parameter choices and ranges, specifically when the fitted parameters fall on the bound of these ranges is necessary. If there is a specific biological reason to keep the bounds as they are that should be noted. If there is not a biological reason, the bounds should be widened to show that the parameters fitted, and thus the major conclusions, are not impacted by their choice of bounds.

- Throughout the paper there is discussion of "sexual commitment rate" but reference is always to the percentage commitment. That distinction should be clear and consistent throughout the manuscript. (See Table 1; subsection “Predicting the impact of gametocyte kinetics on human-to-mosquito transmissibility”; Figure 4CD axes labels; Discussion section; Table 2; subsection “Fitting the model to parasitaemia data”; Figure 1—figure supplement 15A axes label; Figure1—figure supplement 1 caption; Figure 1—figure supplement 2 caption; Figure 1—figure supplement 10 caption).

2) Statistical quantification in subsection “Model fitting and validation”, can "very well" be quantified? Can "excellent predictive" be quantified? Also, "very persistent" seems in contradiction of the previous paragraph of discussion of discrepancies. Why is there no discussion in the text of Figure 3? It seems a major point of the paper. How does one validate "visually capture the data" in the caption of Appendix 1?

3) Predicting likelihood of transmission. Although a reference is given for the choice of 108 parasite/mL for newly hospitalized cases, this number differs from the values seen in other important references (Eichner et al.,). For the non-infectious period, the 103 parasites/mL is listed as a value below which there is no transmission. As written in the manuscript it sounds like transmission is likely above this value. Furthermore, Figure 4 is confusing. What is the line in D? The only mention refers to when Gc=103, which is all of the values in D. The scaling (log vs linear) in these figures is confusing. Why use log on the fraction of sexual commitment? Why use linear on the Gc value?

4) Subsection “Predicting the impact of gametocyte kinetics on human-to-mosquito transmissibility” seem to suggest that, when determining the point during the infection at which the patient is hospitalised, total parasite density is a better indicator that asexual parasite density. This statement surprised me: a reference to Saralambda et al., 2010 is provided, but I can't find any mention of this. Have I understood the statement correctly, and is there any evidence of this? This is an important point, as the results presented in Figure 4 depend on the determination of the time of hospitalisation (which presumably is a proxy for the patient becoming febrile). The asexual parasite population is responsible for the rupturing of red blood cells every 48 hours, which is often linked to symptomatic malaria. The statement in subsection “Predicting the impact of gametocyte kinetics on human-to-mosquito transmissibility” seems to be contradicted in the Discussion section, which is a bit confusing.

5) In subsection “Gametocyte dynamics model” the authors state that piperaquine does not kill immature or mature gametocytes, although there is no reference for this. in vitro evidence (S.H. Adjalley et al., 2011) suggests that the drug does have some effect on very young gametocytes (stages I and II). I think it is reasonable to neglect these effects in the model, but less reasonable to state that there aren't any (I appreciate that I'm being a bit fussy here). Neglecting this drug effect may lead to an increased estimate for the death rate of immature gametocytes, but this is just speculation on my part and I'm not suggesting that the model should be re-fitted at this stage.

6) In Figure 2, several of the panels contain data points that look to be below the stated limit of detected (e.g. all data points after day 14 for Volunteer 202). I imagine that these represent zeros (data points for which the parasitaemia was below the limit of detection), this should be stated somewhere if so.

7) I'm a bit confused by the red line in Figure 4D. The caption of the figure states that, "The red curves indicate the cases corresponding to gametocytaemia of 103 parasites/mL". But my understanding of Figure 4D is that at every point (of the 2D surface) gametocytaemia has reached 103 gametocytes / mL. Have I understood this correctly? It could be that the red line on this panel is a contour of constant tc, which should be clearly stated if this is the case. *The results presented in Figure 4C suggest that newly hospitalised malaria patients are unlikely to be infectious. How do these results compare with clinical trial data? Baseline gametocytaemia is routinely recorded in clinical trials of uncomplicated malaria. It would be interesting to see how the results compare (and why they might be different).

8) One limitation of the study is that the circulation time of mature gametocytes cannot be estimated with any accuracy, due to the lack of patient follow-up. The authors do acknowledge this, but I think the statement (subsection “Estimation of gametocyte dynamics parameters”), "… we found that the circulating gametocyte lifespan… was much longer than that estimated from the neurosyphilis patient data…" is too strong and should be adjusted. In particular, using a constant hazard for gametocyte death (which was not found to provide the best model fit to the neurosyphilis patient data) in the case where an adequate follow up period was not available will overestimate the circulation time of the gametocytes.

9) In Table 2, I wonder if the prior on parameter f should be 0-100%, not 0-1%. Figure 1—figure supplement 1 and Figure 1—figure supplement 10 suggests that the relevant marginal posteriors extend beyond 0.01 (here f not converted to a percentage). Furthermore, in Table 1 the authors compare the commitment rate to the much higher values obtained elsewhere, particularly in vitro. It would be a curious comparison, if higher values were excluded from the model by using a uniform prior between 0-1%.

10) In Figure 1—figure supplement 1, I was initially confused by the high values for rp. Toward the end of the Methods section, the authors do explain how to convert this value to a net multiplication rate, which does clarify the matter. This parameter could be called (e.g.) the raw multiplication rate, but I leave this to the authors' discretion.

https://doi.org/10.7554/eLife.49058.032

Author response

Summary:

In this work, the authors have extended a within-host model of a P. falciparum infection to include gametocyte production. The model is calibrated to parasitaemia from 17 volunteers in a controlled human infection study. Although well thought out and well written, a number of concerns remain: (1) several clarifications need to be made in the description of the models; (2) additional justification is needed for the statistical validation of the results; (3) the section on predicting the impact of gametocyte kinetics seems poorly supported.

We have made a number of changes to improve the manuscript based on the above comments. We address each of them in detail below and very much appreciate the editors and reviewers’ time for providing the helpful suggestions.

Essential revisions:

1) Model description

- a and t in the model are never defined (subsection “Study population and measurements”).

The variables a and t were defined in the originally submitted manuscript (Line 350). However, we note that Figure 5, which introduced the notation, was referred to before these definitions. We have revised the caption of Figure 5 in the revised manuscript to ensure that notation is defined on first use:

“The model is comprised of three parts describing three populations of parasites: asexual parasites (Pa,t), sexually committed parasites (PGa,t) and gametocytes (Gt). P and PG are functions of asexual parasite age a and time t. Square compartments in the inner loop represent the asexual parasite population which follows a cycle of maturation and replication every aL hours. Sexual commitment occurs from age as and a fraction of asexual parasites become sexually committed (the bigger square compartments in the outer loop) and eventually enter the development of stage I–V gametocytes (G1–G5). The compartments with a dashed boundary are sequestered to tissues and thus not measurable in a blood smear. The notation for each compartment is consistent with those inthe model equations and is explained in the main text.”

Note that the text in subsection “Study population and measurements” of the originally submitted manuscript has not been changed.

- γ is never defined (subsection “Gametocyte dynamics model”).

γ is now defined in the revised manuscript (subsection “Gametocyte dynamics model”):

“…, and γ is the Hill coefficient determining the curvature of the dose-response curve”

- Where do the parameters ranges from Table 2 originate? Why are the death rates so high? For example, the maximum death rate of gametocytes is 0.1/hr, suggesting an average life span of 10 hours.

The parameter ranges shown in Table 2 are for the (Bayesian) prior distributions on the model’s parameters. They should not be interpreted as the biological estimates for the parameters. The biological estimates for the parameters, informed by the VIS data, are given by the Bayesian posterior distributions. The ranges for the prior distributions are conventionally chosen such that they are sufficiently wide to accommodate the estimates in the literature (e.g. for our case some parameter ranges were chosen to be wider than the in vivo values given in Zaloumis et al., 2012) but also confined by biological plausibility (e.g. for our model the lower bounds were all non-negative because of the definitions of the model parameters). Since there is no gold standard for determining how wide the ranges of the (uniform) priors should be, we chose relatively large values for the upper bounds. For example, as the reviewers pointed out, the upper bound of gametocyte death rate was chosen to be 0.1/h (i.e. a mean life span of 10 hours), which is much higher than the maximum of our posterior estimates ~0.009/h (see Figure 1—figure supplement 1). This is actually a clear demonstration of how the different ranges between the prior distribution and the posterior distribution indicate that the (VIS) data holds considerable information on the value for this parameter. We have adjusted the caption of Table 2 in the revised manuscript:

“The table includes the unit, description and prior distribution for each model parameter. For the uniform prior distributions (U), the lower bounds are non-negative based on the definitions of the model parameters and the upper bounds for the prior distributions were chosen to be sufficiently wide in order to accommodate all biologically plausible values from the literature (Zaloumis et al., 2012). We assumed parasites younger than 25h are circulating and thus fix as to be 25h. For 3D7 strain, the asexual replication cycle is approximately 39–45h (based on in vitro estimates (Duffy and Avery, 2017) and personal communication (JS McCarthy, personal communication, May 2019)) and we fix aL to be 42h.”

- Although never stated in this manuscript, it seems that the inoculum from the Collins paper was 2800 parasites. Discuss how the range of 0-10 parasites/mL relates to this.

Once again, we note that our range of 0–10 is for the prior distribution, not an estimate for the range of the actual value. Furthermore, there is no simple conversion from the inoculum of 2800 infected red blood cells to the initial parasite density due to a number of factors, such as blood volume, stochastic variation in the proportion of viable parasites, heterogeneity of the parasites in the blood, inter-individual variation in parasite growth rates etc. We chose an upper bound on the prior of 10 because the initial total parasite load was below the detection limit of 10 parasites/mL. Note that even on day 4 where the total parasitaemia should have been at least two orders of magnitude more than the initial parasitaemia, none of the volunteers had a parasitaemia level > 100 parasites/mL, suggesting that the initial parasite load must be far less than 10. As above, the credible range of the posterior distribution is what constitutes the biological contribution of our analysis. The posterior for the initial parasite load Pinit was found to be 0.0232–0.2948, as seen in Figure 1—figure supplement 1, indicating that the chosen upper bound for the prior was both appropriate and had an absolutely negligible effect on the posterior distribution.

- Vc, V1 and V2 are never defined (subsection “Pharmacokinetic model of piperaquine (PQP)”).

Thank you for identifying this oversight. We have explained them in subsection “Pharmacokinetic model of piperaquine (PQP)” in the revised manuscript:

“…, where kT and q’s are rate constants as shown in Figure 6 and Vc, V1 and V2 are the volume of distribution for the central compartment (in which PQP concentration is C), peripheral compartment 1 (in which PQP concentration is P1) and peripheral compartment 2 (in which PQP concentration is P2) respectively.”

- How are the initial conditions of the drug model chosen?

Thank you for identifying this oversight. The initial conditions for PK model simulation were that all the compartments (i.e. D, T1, T2, C, P1and P2) were zeroes because PQP was not given until day 7 or 8. When the first dose of 480mg PQP was given, we set D = 480 and zero for other compartments. When a second dose of 960mg was required for some volunteers to treat recrudescence, we set D to be 960 plus any remaining level of D from the first dose for those volunteers. Note our PK model (as described in the Section “Pharmacokinetic model of piperaquine (PQP)”) explicitly models the absorption of administrated drug mass in compartment D, and how that drives drug concentration in the central compartment C, and so the second dose is modelled as a step change in D, not a step change in C. We now describe how the initial conditions are determined in the Appendix 1:

“The PK model (see the Materials and methods section) describes absorption of administrated drug mass in compartment D and the subsequent kinetics of drug concentration in the central compartment C. The initial conditions for the model simulation were that all the compartments (i.e. D, T1, T2, C, P1 and P2) were zero because PQP was not given until day 7 or 8. When the first dose of 480mg PQP was given, we set D = 480 and all other compartments to zero. When a second dose of 960mg was required for some volunteers to treat recrudescence, we again set D to be 960 plus any remaining level of D from the first dose for those volunteers.”

and in the main text in the revised manuscript (subsection “Fitting the model to parasitaemia data”):

“The details of the initial conditions, starting point and constraints for the PK curve optimisation procedure are provided in Appendix 1.”

Why are there two spikes for some patients and only one for others (Figure 1—figure supplement 1)?

A second dose of PQP was given to the volunteers showing the recrudescence. While this was mentioned in the first section of the Materials and methods section in the originally submitted manuscript, we now also mention it in the caption of Figure 6—figure supplement 1:

“Some volunteers have two peaks of PQP concentrations because they had recrudescent asexual parasitaemia (see Figure 1 in the main text) and were treated with a second dose of 960mg PQP.”

- In the fitting of the PK model, was only a single starting set considered? Why is it necessary and valid to increase the upper bound for qc? For different parameters, several patients hit the upper or lower bounds. In particular, Volunteers 105 and 20 hit the bounds on almost all parameters. Could their optimal fit be outside the range?

These questions all concern the PK parameters used for each volunteer to generate PQP concentration data for our gametocyte dynamics model fitting. The Volunteer Infection Study (VIS) on which our paper is based was not designed to accurately determine the PK parameters for each trial participant. Nonetheless, we did have some PK data for each volunteer. For our gametocyte dynamics model fitting, our primary aim was to determine an appropriate PK curve for each volunteer (rather than a unique estimate of their PK parameters) based on the limited data. That is, we have taken a sequential PK-PD approach, where predicted PQP concentration (i.e. C(t)) for each individual is included as a variable in the PD model. To achieve this requirement, we drew on an analysis of rich PK data for PQP from a different VIS, where the PQP’s PK model parameter estimates were determined by our collaborators at Mahidol-Oxford Tropical Medicine Research Unit (Bangkok), Thanaporn Wattanakul and Joel Tarning.

The approach we took, as described (and now further clarified) in the Section “Pharmacokinetic model of piperaquine (PQP)” and the Appendix 1, was to use the MATLAB built-in least-squares optimiser lsqcurvefit (with the default setting) to optimise the PK curve for each volunteer in the VIS study. We applied the optimiser to each volunteer’s (limited) PQP data, using the parameter estimates provided by Thanaporn Wattanakul and Joel Tarning as initial values. We applied some further model parameter constraints as specified in the Appendix 1. This approach provided us with the required data-informed PK curve for each volunteer in the VIS, sufficient for our sequential PK-PD study focussing on the asexual and sexual parasite dynamics. Of note, Volunteers 202, 301, 302 or 307 had fewer PK data points than PK model parameters, preventing application of this optimisation procedure. For these volunteers, their predicted PQP PK curve was derived using the population-level mean parameter estimates from Wattanakul and Tarning’s analysis. We have modified the manuscript to clarify the method for producing appropriate PK curves (subsection “Pharmacokinetic model of piperaquine (PQP)”):

“Under the sequential pharmacokinetic-pharmacodynamic (PK-PD) approach we have taken, a PQP concentration curve (C(t)) for each volunteer is a required input into the gametocyte dynamics model. The VIS, with its limited sampling of PQP for each volunteer, was not designed to provide this PQP concentration curve directly, so we used a PK model, informed by data from a previous VIS with rich sampling. We drew on an analysis of that previous VIS by Thanaporn Wattanakul and Joel Tarning (unpublished data and estimates). Their analysis provides population-level PQP PK model parameter estimates.

We used MATLAB’s (version 2016b; The MathWorks, Natick, MA) built-in least-squares optimiser lsqcurvefit (with the default setting) to optimise the PK curve for each volunteer in the VIS study. We applied the optimiser to each volunteer’s (limited) PQP data, using the parameter estimates provided by Thanaporn Wattanakul and Joel Tarning as initial values. We applied some further model parameter constraints as specified in Appendix 1. This approach provided us with a data-informed PK curve for each volunteer in the VIS, sufficient for our primary purpose of studying the asexual and sexual parasite dynamics. Of note, Volunteers 202, 301, 302 or 307 had fewer PK data points than PK model parameters, preventing application of this optimisation procedure. For these volunteers, their predicted PQP PK curve was derived using the population-level mean PK parameter from Wattanakul and Tarning’s analysis. The MATLAB code (with detailed comments) is publicly available at https://doi.org/10.26188/5cde4c26c8201. The details of the initial conditions, starting point and constraints for the PK curve optimisation procedure are provided in Appendix 1. The optimised PK curves and associated parameter values for all volunteers are provided in Figure 6-figure-supplement 1 and Appendix 1.”

We increased qc in order to generate reasonable PK curves. Because some volunteers (e.g. 102, 104, 105, 201, 306, etc.) exhibited a fast drug decay, it was necessary to relax the upper bound of qc to improve the fit to the (limited) available data. We have explained this in Appendix 1:

“… except that the upper bound for the clearance rate qc was increased from 62.6 to 160 such that the PK curves could better capture the fast PQP concentration decay observed for some volunteers.”

- Justification for all parameter choices and ranges, specifically when the fitted parameters fall on the bound of these ranges is necessary. If there is a specific biological reason to keep the bounds as they are that should be noted. If there is not a biological reason, the bounds should be widened to show that the parameters fitted, and thus the major conclusions, are not impacted by their choice of bounds.

Our previous response to the questions of parameter estimation – and that what we require for our analysis of the VIS data is reasonable PK curves, not PK parameter estimates – has partly addressed this concern. For the PK model, the bounds were chosen to be the 95% confidence interval of the PK parameter estimates provided by Thanaporn Wattanakul and Joel Tarning, with some relaxation on the upper bound of the clearance rate 𝑞𝑐 as mentioned above.

We have adjusted the text in Appendix 1:

“The lower and upper bounds are the limits of the 95% confidence intervals of the parameter estimate distributions provided by Thanaporn Wattanakul and Joel Tarning except that the upper bound for the clearance rate 𝑞𝑐 was increased from 62.6 to 160 such that the PK curves could better capture the fast PQP concentration decay observed for some volunteers. The constraints are necessary due to the limited available PK data which prevents identification of the PK parameters.”

- Throughout the paper there is discussion of "sexual commitment rate" but reference is always to the percentage commitment. That distinction should be clear and consistent throughout the manuscript. (See Table 1; subsection “Predicting the impact of gametocyte kinetics on human-to-mosquito transmissibility”; Figure 4CD axes labels; Discussion section; Table 2; subsection “Fitting the model to parasitaemia data”; Figure 1—figure supplement 15A axes label; Figure1—figure supplement 1 caption; Figure 1—figure supplement 2 caption; Figure 1—figure supplement 10 caption).

Thank you for identifying this ambiguity in nomenclature. By rate, we refer to a “per replication cycle rate”. The revised manuscript uses a new terminology: the sexual commitment rate is now consistently expressed as “percentage/asexual replication cycle” (e.g. 0.54%/asexual replication cycle). The parameter f indicates the fraction and thus lies between 0 and 1. The conversion between the sexual commitment rate and parameter f is “sexual commitment rate = f×100%”. We have corrected any ambiguous expressions and provided the conversion method (subsection “Fitting the model to parasitaemia data”) in the revised manuscript.

2) Statistical quantification in subsection “Model fitting and validation”, can "very well" be quantified? Can "excellent predictive" be quantified? How does one validate "visually capture the data" in the caption of Supplementary file 1?

In the Bayesian analysis literature, much of model validation is qualitative, based on a combination of measures and analysis including the posterior predictive check. In our context, “very well” or "excellent predictive capability" means the predicted parasitaemia (i.e. the median and 95% prediction interval) is able to accurately capture the detailed dynamics/trajectory of the data. This was primarily judged by visual check rather than based on quantities measuring the goodness of fit, such as information criteria, which are often less meaningful (i.e. no better than visual check) unless used for model comparison purposes. We have revised the manuscript (Introduction):

“The results show that the predicted total parasitaemia (median and 95% PI) is able to accurately capture the trends of the data through the (visual) posterior predictive check.”

and (subsection “Model fitting and validation”):

“For the majority of asexual parasitaemia data the model predictions (median and 95% PI) can faithfully capture the trends of the data (Figure 2)”.

For determining patient-specific PQP concentrations at the times when parasitaemia was measured (which is an essential step of our sequential PK-PD approach), we aimed to generate PK curves which are visually close to the limited PK data (note that the difficulty of fitting to the limited PK data was explained in detail in our response to comment (1)), and one can validate the quality of the predicted PQP concentrations by looking at Figure 6—figure supplement 1 where both the predicted PQP profiles and available data are provided for all the volunteers. We have replaced “visually capture the data” by “visually capture the data (as shown in Figure 6—figure supplement 1 where the predicted PQP concentrations are very close to the observed concentrations)” in the Appendix 1 to minimise the ambiguity.

In the subsection “Model fitting and validation”, "very persistent" seems in contradiction of the previous paragraph of discussion of discrepancies. Why is there no discussion in the text of Figure 3? It seems a major point of the paper.

We begin by noting that presumably the reviewers are asking about the use of the phrase “very consistent” rather than (as typed) “very persistent”. Figure 3 was discussed, albeit briefly, in subsection “Model fitting and validation” of the originally submitted manuscript (which is the text being queried by the reviewers).

In the originally submitted manuscript, Figure 3 shows that the predicted gametocytaemia for all volunteers is very close to the measured data for those volunteers, hence our description as “very consistent”. Note that we acknowledged that the model’s predictions of the asexual parasitaemia (not the gametocytaemia) for some volunteers were less consistent (subsection “Model fitting and validation” which refers to Figure 2).

We have modified the manuscript to ensure that these very important results in Figure 3 are clearly emphasised for the reader. We have modified the discussion of Figure 3 and highlighted it as a new paragraph in the revised manuscript (subsection “Model fitting and validation”):

“Figure 3 shows the data and model predictions for the gametocytaemia. Despite some discrepant observations for asexual parasitaemia in Figure 2, we found that the model predictions of gametocytaemia were able to capture the trends and levels of the gametocytaemia data for all 17 volunteers.”

3) Predicting likelihood of transmission. Although a reference is given for the choice of 108 parasite/mL for newly hospitalized cases, this number differs from the values seen in other important references (Eichner et al.,). For the non-infectious period, the 103 parasites/mL is listed as a value below which there is no transmission. As written in the manuscript it sounds like transmission is likely above this value. Furthermore, Figure 4 is confusing. What is the line in D? The only mention refers to when Gc=103, which is all of the values in D. The scaling (log vs linear) in these figures is confusing. Why use log on the fraction of sexual commitment? Why use linear on the Gc value?

We used a unit of “parasites per milli litre” while Eichner et al. used “parasites per micro litre”. A unit conversion of 108 parasites/mL gives 105 parasites/uL which matches the level shown in Eichner et al.

For the threshold of 103, we have adjusted subsection “Predicting the impact of gametocyte kinetics on human-to-mosquito transmissibility” in the revised manuscript:

“…, which is a threshold below which human-to-mosquito transmission was not observed (Collins et al., 2018).”

We have provided the definition of the red curve in Figure 4D in the revised manuscript (subsection “Predicting the impact of gametocyte kinetics on human-to-mosquito transmissibility”):

“The red curve in C is the level curve for Gc = 103 parasites/mL. The red curve in D is the level curve for tc= 13.42 days which is the non-infectious period obtained by model simulation using the posterior estimates of the population mean parameters.”

We used a log scale for the sexual commitment rate because the plausible range of the rate spanned a few orders of magnitude (i.e. from order 0.1% to 10%; our in vivo estimate was about 0.5% while the in vitro estimate is about 10%). For Gc, although parasitaemia is typically shown on a log scale, the variation in Gc was roughly linear (see the colour bar in Figure 4C) and therefore we used a linear scale for Gc.. Based on these reviewers’ queries, we have considered alternative plots but have decided that our current visualisations are the most appropriate.

4) Subsection “Predicting the impact of gametocyte kinetics on human-to-mosquito transmissibility” seem to suggest that, when determining the point during the infection at which the patient is hospitalised, total parasite density is a better indicator that asexual parasite density. This statement surprised me: a reference to Saralambda et al., [2010] is provided, but I can't find any mention of this. Have I understood the statement correctly, and is there any evidence of this? This is an important point, as the results presented in Figure 4 depend on the determination of the time of hospitalisation (which presumably is a proxy for the patient becoming febrile). The asexual parasite population is responsible for the rupturing of red blood cells every 48 hours, which is often linked to symptomatic malaria. The statement in subsection “Predicting the impact of gametocyte kinetics on human-to-mosquito transmissibility” seems to be contradicted in the Discussion section, which is a bit confusing.

Thank you for pointing out the confusing statement (Subsection “Predicting the impact of gametocyte kinetics on human-to-mosquito transmissibility”) and a typographical error in the Discussion section. First of all, “asexual parasitaemia” should be “total parasitaemia”, as the total parasitaemia was consistently used in the simulation of the clinical scenario and in the Results section. We have changed “asexual” to “total” in the revised manuscript (Discussion section).

The reference to Saralamba et al., 2010 was to support the choice of 108 parasites/mL as a realistic value for the total parasitaemia of newly hospitalised patients. Since the parasitaemia data in Saralamba et al., 2010 was measured under microscope where both asexual and sexual rings were counted, we chose total parasitaemia to simulate the clinical scenario. We have modified the statement (Lines 199-200 in the originally submitted manuscript) in the revised manuscript to be (subsection “Predicting the impact of gametocyte kinetics on human-to-mosquito transmissibility”):

“We further assumed that patients would seek hospital admission when their total parasitaemia reached approximately 108 parasites/mL. This choice was based on the microscopic measurements of the total parasitaemia (i.e. asexual and sexual parasites) from a study of Cambodia and Thailand hospitalised malaria patients (Saralamba et al., 2010).”

5) In subsection “Gametocyte dynamics model” the authors state that piperaquine does not kill immature or mature gametocytes, although there is no reference for this. in vitro evidence (S.H. Adjalley et al., 2011) suggests that the drug does have some effect on very young gametocytes (stages I and II). I think it is reasonable to neglect these effects in the model, but less reasonable to state that there aren't any (I appreciate that I'm being a bit fussy here). Neglecting this drug effect may lead to an increased estimate for the death rate of immature gametocytes, but this is just speculation on my part and I'm not suggesting that the model should be re-fitted at this stage.

We agree that our statement was inaccurate and have modified the manuscript (subsection “Gametocyte dynamics model”):

“Note that we assumed in our model that PQP does not kill gametocytes. Our assumption was based on evidence from both in vitro and in vivo experiments that suggests that PQP has little activity against sexually committed parasites and gametocytes (Bolscher et al., 2014; Collins et al., 2018; Pasay et al., 2016), although we note there is some evidence that PQP may have activity against early-stage I/II gametocytes (Adjalley et al., 2011).”

6) In Figure 2, several of the panels contain data points that look to be below the stated limit of detected (e.g. all data points after day 14 for Volunteer 202). I imagine that these represent zeros (data points for which the parasitaemia was below the limit of detection), this should be stated somewhere if so.

This interpretation is correct in that the points are those with a read-out of below the limit of detection, but we cannot ever say that the values are actually zero (as in no parasites at all). The data points placed at one parasite/mL (or zero on a log scale) represent the measurements for which no parasites were detected. We now clearly state this in the caption of Figure 2 in the revised manuscript:

“The data points with one parasite/mL (i.e. those points which lie on the dotted line) indicate measurements for which no parasites were detected.”

7) I'm a bit confused by the red line in Figure 4D. The caption of the figure states that, "The red curves indicate the cases corresponding to gametocytaemia of 103 parasites/mL". But my understanding of Figure 4D is that at every point (of the 2D surface) gametocytaemia has reached 103 gametocytes / mL. Have I understood this correctly? It could be that the red line on this panel is a contour of constant tc, which should be clearly stated if this is the case.

Thank you for raising this issue of interpretation, which was also raised in comment (3) above. We have addressed it in the revised manuscript (Figure 4 caption):

“The red curve in C is the level curve for Gc = 103 parasites/mL. The red curve in D is the level curve for tc= 13.42 days which is the non-infectious period obtained by model simulation using the posterior estimates of the population mean parameters.”

The results presented in Figure 4C suggest that newly hospitalised malaria patients are unlikely to be infectious. How do these results compare with clinical trial data? Baseline gametocytaemia is routinely recorded in clinical trials of uncomplicated malaria. It would be interesting to see how the results compare (and why they might be different).

This is a good suggestion. Patients from field clinical observations of uncomplicated malaria may carry more gametocytes than what our model predicts at the time of presentation, for example one study from the TRACII clinical trial reported a range of 16–5120 gametocytes/µL, much higher than our prediction of below 1 gametocytes/µL (or 103 gametocytes/mL) (Pluijm et al., 2019). One explanation for the difference would be that our model predicted a very fast rise in parasitaemia to 108 parasites/mL while patients in clinical trials with naturally occurring malaria may take a longer time to reach such a high parasitaemia level due to immunity. Immunity was, of course, not considered in our model due to the design of our Volunteer Infection Studies where only malaria naïve volunteers are recruited. Thus, we have added the following in the revised manuscript (subsection “Predicting the impact of gametocyte kinetics on human-to-mosquito transmissibility”) to address this comment:

“Note that patients from clinical observations of uncomplicated malaria in endemic settings may have higher gametocyte counts at the time of presentation than what our model predicts. For example one study from the TRACII clinical trial reported a range of 16–5120 gametocytes/µL, which is much higher than our prediction of below 1 gametocytes/µL (or 103 gametocytes/mL) (van der Pluijm et al., 2019). One plausible explanation for the difference is that our model predicted a very fast rise in total parasitaemia to 108 parasites/mL while the rise in parasitaemia among patients in endemic settings may be slower due to the effect of immunity on the parasite multiplication. Immunity was not considered in our model due to the design of our VIS where only malaria naïve volunteers were recruited.”

8) One limitation of the study is that the circulation time of mature gametocytes cannot be estimated with any accuracy, due to the lack of patient follow-up. The authors do acknowledge this, but I think the statement (subsection “Estimation of gametocyte dynamics parameters”), "… we found that the circulating gametocyte lifespan… was much longer than that estimated from the neurosyphilis patient data…" is too strong and should be adjusted. In particular, using a constant hazard for gametocyte death (which was not found to provide the best model fit to the neurosyphilis patient data) in the case where an adequate follow up period was not available will overestimate the circulation time of the gametocytes.

We fully accept this critique of our original comment. Noting that our comments in the Discussion section emphasise that we cannot make a meaningful estimate of this parameter based on the VIS data, we have decided to simply delete this statement from the revised manuscript.

9) In Table 2, I wonder if the prior on parameter f should be 0-100%, not 0-1%. Figure 1—figure supplement 1 and Figure 1—figure 10 suggests that the relevant marginal posteriors extend beyond 0.01 (here f not converted to a percentage). Furthermore, in Table 1 the authors compare the commitment rate to the much higher values obtained elsewhere, particularly in vitro. It would be a curious comparison, if higher values were excluded from the model by using a uniform prior between 0-1%.

Yes, the reviewers understood correctly that the prior on parameter f was defined to be 0-1 (i.e. 0-100%), which allowed the range of posterior distribution to extend beyond 0.01 (i.e. 1%) shown in Figure 1—figure supplement 1 and Figure 1—figure supplement 10. In order to avoid any confusion between the parameter f (which indicates the fraction) and the sexual commitment rate (which was expressed as percentage), as for our response to an earlier comment, we have made adjustments throughout the manuscript.

10) In Figure 1—figure supplement 1, I was initially confused by the high values for rp. Toward the end of the Methods section, the authors do explain how to convert this value to a net multiplication rate, which does clarify the matter. This parameter could be called (e.g.) the raw multiplication rate, but I leave this to the authors' discretion.

We appreciate that the reviewers have raised this matter which we believe many readers familiar with the details of malaria kinetics may also note. We have added the following text immediately following the introduction of the parasite replication rate rP in the revised manuscript (subsection “Gametocyte dynamics model”):

“Note that we distinguish the parasite replication rate rP from the so-called parasite multiplication factor, the latter of which is a “net multiplication rate” quantified by the (per cycle) increase in parasite numbers due to replication (rP) and the decrease in parasite numbers due to death or sexual commitment.”

https://doi.org/10.7554/eLife.49058.033

Article and author information

Author details

  1. Pengxing Cao

    School of Mathematics and Statistics, University of Melbourne, Melbourne, Australia
    Contribution
    Conceptualization, Formal analysis, Visualization, Methodology, Writing—original draft
    For correspondence
    pengxing.cao@unimelb.edu.au
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3721-9850
  2. Katharine A Collins

    1. Department of Medical Microbiology, Radboud University Medical Center, Nijmegen, Netherlands
    2. QIMR Berghofer Medical Research Institute, Brisbane, Australia
    Contribution
    Data curation, Formal analysis, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7080-2215
  3. Sophie Zaloumis

    Centre for Epidemiology and Biostatistics, Melbourne School of Population and Global Health, University of Melbourne, Melbourne, Australia
    Contribution
    Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8253-8896
  4. Thanaporn Wattanakul

    1. Mahidol-Oxford Tropical Medicine Research Unit, Faculty of Tropical Medicine, Mahidol University, Bangkok, Thailand
    2. Centre for Tropical Medicine and Global Health, Nuffield Department of Medicine, University of Oxford, Oxford, United Kingdom
    Contribution
    Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7570-4665
  5. Joel Tarning

    1. Mahidol-Oxford Tropical Medicine Research Unit, Faculty of Tropical Medicine, Mahidol University, Bangkok, Thailand
    2. Centre for Tropical Medicine and Global Health, Nuffield Department of Medicine, University of Oxford, Oxford, United Kingdom
    Contribution
    Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4566-4030
  6. Julie A Simpson

    Centre for Epidemiology and Biostatistics, Melbourne School of Population and Global Health, University of Melbourne, Melbourne, Australia
    Contribution
    Conceptualization, Formal analysis, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2660-2013
  7. James McCarthy

    QIMR Berghofer Medical Research Institute, Brisbane, Australia
    Contribution
    Conceptualization, Data curation, Formal analysis, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6596-9718
  8. James M McCaw

    1. School of Mathematics and Statistics, University of Melbourne, Melbourne, Australia
    2. Centre for Epidemiology and Biostatistics, Melbourne School of Population and Global Health, University of Melbourne, Melbourne, Australia
    3. Epidemiology, Peter Doherty Institute for Infection and Immunity, Parkville, Australia
    Contribution
    Conceptualization, Formal analysis, Funding acquisition, Writing—review and editing
    For correspondence
    jamesm@unimelb.edu.au
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2452-3098

Funding

Australian Research Council (Discovery project DP170103076)

  • James M McCaw

National Health and Medical Research Council (1025319)

  • Julie A Simpson

National Health and Medical Research Council (Centre for Research Excellence 1134989)

  • Julie A Simpson
  • James McCarthy
  • James M McCaw

Australian Research Council (Early Career Researcher Award DE170100785)

  • Sophie Zaloumis

National Health and Medical Research Council (Senior Research Fellowship 1104975)

  • Julie A Simpson

Bill and Melinda Gates Foundation

  • Thanaporn Wattanakul
  • Joel Tarning

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We acknowledge useful conversations with David S Khoury, Deborah Cromer and Miles P Davenport (Kirby Institute, UNSW Australia, Sydney, Australia) and assistance in drafting the manuscript from Laura Cascales (QIMR Berghofer, Brisbane, Australia). We thank Jörg J Möhrle, Head of Translation from Medicines for Malaria Venture for his support and permission to use the trial data. This research was supported by use of the Nectar Research Cloud, a collaborative Australian research platform supported by the National Collaborative Research Infrastructure Strategy (NCRIS).

Senior Editor

  1. Neil M Ferguson, Imperial College London, United Kingdom

Reviewing Editor

  1. Isabel Rodriguez-Barraquer, University of California, San Francisco, United States

Publication history

  1. Received: June 5, 2019
  2. Accepted: October 15, 2019
  3. Version of Record published: October 29, 2019 (version 1)

Copyright

© 2019, Cao et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,069
    Page views
  • 289
    Downloads
  • 22
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Pengxing Cao
  2. Katharine A Collins
  3. Sophie Zaloumis
  4. Thanaporn Wattanakul
  5. Joel Tarning
  6. Julie A Simpson
  7. James McCarthy
  8. James M McCaw
(2019)
Modeling the dynamics of Plasmodium falciparum gametocytes in humans during malaria infection
eLife 8:e49058.
https://doi.org/10.7554/eLife.49058

Further reading

    1. Computational and Systems Biology
    Valentina Baldazzi, Delphine Ropers ... Hidde de Jong
    Research Article

    Different strains of a microorganism growing in the same environment display a wide variety of growth rates and growth yields. We developed a coarse-grained model to test the hypothesis that different resource allocation strategies, corresponding to different compositions of the proteome, can account for the observed rate-yield variability. The model predictions were verified by means of a database of hundreds of published rate-yield and uptake-secretion phenotypes of Escherichia coli strains grown in standard laboratory conditions. We found a very good quantitative agreement between the range of predicted and observed growth rates, growth yields, and glucose uptake and acetate secretion rates. These results support the hypothesis that resource allocation is a major explanatory factor of the observed variability of growth rates and growth yields across different bacterial strains. An interesting prediction of our model, supported by the experimental data, is that high growth rates are not necessarily accompanied by low growth yields. The resource allocation strategies enabling high-rate, high-yield growth of E. coli lead to a higher saturation of enzymes and ribosomes, and thus to a more efficient utilization of proteomic resources. Our model thus contributes to a fundamental understanding of the quantitative relationship between rate and yield in E. coli and other microorganisms. It may also be useful for the rapid screening of strains in metabolic engineering and synthetic biology.

    1. Computational and Systems Biology
    2. Neuroscience
    Kai J Sandbrink, Pranav Mamidanna ... Alexander Mathis
    Research Article

    Biological motor control is versatile, efficient, and depends on proprioceptive feedback. Muscles are flexible and undergo continuous changes, requiring distributed adaptive control mechanisms that continuously account for the body's state. The canonical role of proprioception is representing the body state. We hypothesize that the proprioceptive system could also be critical for high-level tasks such as action recognition. To test this theory, we pursued a task-driven modeling approach, which allowed us to isolate the study of proprioception. We generated a large synthetic dataset of human arm trajectories tracing characters of the Latin alphabet in 3D space, together with muscle activities obtained from a musculoskeletal model and model-based muscle spindle activity. Next, we compared two classes of tasks: trajectory decoding and action recognition, which allowed us to train hierarchical models to decode either the position and velocity of the end-effector of one's posture or the character (action) identity from the spindle firing patterns. We found that artificial neural networks could robustly solve both tasks, and the networks'units show tuning properties similar to neurons in the primate somatosensory cortex and the brainstem. Remarkably, we found uniformly distributed directional selective units only with the action-recognition-trained models and not the trajectory-decoding-trained models. This suggests that proprioceptive encoding is additionally associated with higher-level functions such as action recognition and therefore provides new, experimentally testable hypotheses of how proprioception aids in adaptive motor control.