Modeling the dynamics of Plasmodium falciparum gametocytes in humans during malaria infection
Abstract
Renewed efforts to eliminate malaria have highlighted the potential to interrupt humantomosquito 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 datacalibrated model to investigate humantomosquito transmissibility, providing a method to link withinhuman host infection kinetics to epidemiologicalscale infection and transmission patterns.
https://doi.org/10.7554/eLife.49058.001Introduction
Malaria is a mosquitoborne 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 sexualstages 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 humantomosquito 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 bloodmeal, 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 bloodfilms. Previous in vivo estimates of gametocyte kinetic parameters have been primarily based on historical data from neurosyphilis patients who were treated with socalled malariotherapy (Diebner et al., 2000; Eichner et al., 2001). In these studies, the limit of quantification was approximately 10^{4} 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 withinhost 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 datacalibrated model can reliably predict the timecourse 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.
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.
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 nonrecrudescent case where the posteriormedian 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 underpredicts (Volunteer 204) or overpredicts (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 lifestage, 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 reemerge 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 longtailed 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 posteriormedian 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 humantomosquito 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 10^{8} 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 G_{c}) corresponding to the time when the total parasitemia (wavelike black curves) first reached 10^{8} parasites/mL (their associations are indicated by the dotted lines and arrows).
Predicting the noninfectious period of malaria patients for various values of sexual commitment rate and gametocyte sequestration time. In the model, the noninfectious period was defined to be time from the inoculation of infected red blood cells to the time when the gametocytemia reached 10^{3} parasites/mL, which is a threshold below which humantomosquito transmission was not observed (Collins et al., 2018). Note that this noninfectious 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 t_{c}) when the gametocytemia (blue curve) first reached 10^{3} parasites/mL (their associations are indicated by the dotted lines and arrows).
A higher sexual commitment rate or a lower gametocyte sequestration time leads to a higher gametocytemia (G_{c}) at the time of hospitalization (Figure 4C). The red curve in Figure 4C indicates the level curve of 10^{3} gametocytes/mL (i.e., the threshold for infectiousness as mentioned above) dividing the heatmap into two regions. To the left, G_{c} is below 10^{3} gametocytes/mL, suggesting clinical presentation precedes infectiousness, while to the right G_{c} is above 10^{3} gametocytes/mL and the converse applies. The G_{c} value obtained by model simulation using the median estimates of the population mean parameters (indicated by the black dot) is below 10^{3} 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 10^{3} 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 10^{8} 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 noninfectious period (t_{c}). As the sexual commitment rate increases or the gametocyte sequestration time decreases, t_{c} decreases. However, for large values of the sexual commitment rate (e.g., >20%), we observed an increase in t_{c} as the sexual commitment rate increases (see the topright 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 multistate asexual cycle model with a new model for the development of gametocytes. Model parameters were estimated by fitting the model to data from 17 malarianaïve volunteers inoculated with P. falciparuminfected 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 biologicallyinformed 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. AP2G [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 — 10^{8} parasites/mL in the model); and the noninfectious 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 10^{3} parasites/mL in the model). We explored how the sexual commitment rate and gametocyte sequestration time influenced the gametocyte level and the noninfectious period. We would like to emphasize that humantomosquito 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 nonconstant 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 transmissionblocking interventions. Furthermore, the within human host transmission model can be incorporated into epidemiologicalscale models to refine predictions of the impacts of various antimalarial treatments and transmission interventions.
Materials and methods
Study population and measurements
Request a detailed protocolThe data used in this modeling study are from a previously published VIS (Collins et al., 2018) where 17 malarianaïve volunteers were inoculated with approximately 2800 P. falciparuminfected 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 postinoculation 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 qRTPCR (circulating asexual parasites per mL blood), and gametocytemia was measured by Pfs25 and PfMGET qRTPCR (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 protocolThe 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 (P_{G}) 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.
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 ${P}_{init}$ and assume the age distribution of inoculated parasites is Gaussian with mean μ and standard deviation $\sigma $. As time increments by one hour, the asexual parasites of age $a$ at time $t$ (denoted as $P(a,t)$) follow the iterative equation:
where $\overline{{k}_{d}}$ represents the average rate of asexual parasite killing by PQP and ${\delta}_{P}$ is the rate of asexual parasite death due to processes other than PQP. $\overline{{k}_{d}}$ is approximated by the average of ${k}_{d}(t1)$ and ${k}_{d}\left(t\right)$ and ${k}_{d}\left(t\right)={k}_{max}{C\left(t\right)}^{\gamma}/({C\left(t\right)}^{\gamma}+{{EC}_{50}}^{\gamma})$ where ${k}_{max}$ is the maximum killing rate, ${EC}_{50}$ is the PQP concentration at which half maximum killing is achieved, and $\gamma $ is the Hill coefficient determining the curvature of the doseresponse curve. $C\left(t\right)$ is the PQP concentration which is simulated by a pharmacokinetic model introduced below. ${a}_{L}$ is the length of each asexual replication cycle and ${r}_{P}$ 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 ${r}_{P}$ from the socalled parasite multiplication factor, the latter of which is a 'net replication rate' quantified by the (per cycle) increase in parasite numbers due to replication (${r}_{P}$) 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 ${a}_{s}$) and a fraction ($f$) of asexual parasites leave the asexual replication cycle and start sexual development in the next hour, which is modeled by
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 (${P}_{G}$). 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
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 earlystage I/II gametocytes (Adjalley et al., 2011). The changes of the sequestered stage I–IV gametocytes (${G}_{1}$–${G}_{4}$) are governed by difference equations
where m is the rate of gametocyte maturation and ${\delta}_{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
where ${\delta}_{Gm}$ is the death rate of mature circulating gametocytes.
The total parasitemia in the model is given by ${\sum}_{a=1}^{a={a}_{s}1}\left[P\left(a,t\right)+{P}_{G}\left(a,t\right)\right]+{G}_{5}\left(t\right)$, which was fitted to the VIS data. After model fitting, we simulated the asexual parasitemia ${\sum}_{a=1}^{a={a}_{s}1}P(a,t)$ and gametocytemia ${G}_{5}\left(t\right)$ and compared them with associated data for model validation. Table 2 presents all the model parameters and their units and descriptions.
Pharmacokinetic model of piperaquine (PQP)
Request a detailed protocolIn the withinhost model, the killing rate ${k}_{d}\left(t\right)$ is determined by PQP concentration $C\left(t\right)$ which was simulated from a pharmacokinetic (PK) model introduced in this section. The PK model, provided by Thanaporn Wattanakul and Joel Tarning (MahidolOxford Tropical Medicine Research Unit, Bangkok), is a threecompartment disposition model with two transit compartments for absorption (see the schematic diagram in Figure 6).
Based on Figure 6, the model is formulated to be a system of ordinary differential equations:
where ${k}_{T}$ and $q$’s are rate constants as shown in Figure 6 and ${V}_{c}$, ${V}_{1}$ and ${V}_{2}$ are the volume of distribution for the central compartment (in which PQP concentration is C), peripheral compartment 1 (in which PQP concentration is P_{1}) and peripheral compartment 2 (in which PQP concentration is P_{2}) respectively.
Under the sequential pharmacokineticpharmacodynamic (PKPD) 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 populationlevel PQP PK model parameter estimates.
We used MATLAB’s (version 2016b; The MathWorks, Natick, MA) builtin leastsquares 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 datainformed 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 populationlevel 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 protocolWe 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 ${a}_{s}$ and ${a}_{L}$; 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 ${\theta}_{ind}$and lower and upper bounds to be ${b}_{L}$ and ${b}_{U}$ respectively, the following transformations are used to convert the bounded individual parameters to unbounded ones (denoted by ${\phi}_{ind}$) in order to in order to improve computational efficiency (Lesaffre et al., 2007; Stan Development Team, 2017):
o${\phi}_{ind}$ beys a multivariate normal distribution $\mathcal{\mathcal{N}}$(${\phi}_{pop}$, ${\mathrm{\Omega}}_{pop}$) where
and ${\theta}_{pop}$ is a vector containing 12 population mean parameters (hyperparameters) corresponding to the 12 gametocyte dynamics model parameters. ${\mathrm{\Omega}}_{pop}$ is the covariance matrix. For more efficient sampling process, ${\phi}_{ind}~\mathcal{\mathcal{N}}$(${\phi}_{pop}$, ${\mathrm{\Omega}}_{pop}$) was reparameterised to a noncenterd form ${\phi}_{ind}={\phi}_{pop}+{\omega}_{pop}L\eta $, where ${\omega}_{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; $\eta $ obeys standard multivariate normal distribution. Note that ${\mathrm{\Omega}}_{pop}={\omega}_{pop}L{L}^{T}{\omega}_{pop}$ where $L{L}^{T}$ is the correlation matrix. The prior distributions for the 12 population mean parameters ${\theta}_{pop}$ are given by uniform distributions with bounds given in Table 2. The prior distribution for the 12 population SD parameters is standard halfnormal 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 modelsimulated values and SD parameter with prior distribution of a halfCauchy 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 NoUTurn 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 burnin 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 3–14 (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 1–3 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 ${f}_{pop}\times 100\%$ (${f}_{pop}$ is the population mean parameter for $f$) and the proportion of committed asexual parasites that survive to become mature gametocytes was calculated by ${f}_{pop}{({m}_{pop}/({m}_{pop}+{{\delta}_{G}}_{pop}\left)\right)}^{4}$ where the factor of four arises due to the four sequestered gametocyte stages (I to IV). Circulating gametocyte lifespan was calculated by $1/{{\delta}_{Gm}}_{pop}/24$ (the factor of 24 converts hours into days). Gametocyte sequestration time was calculated by $4/{m}_{pop}/24$ where 4 indicates four sequestered state (stage I to IV) and 24 converts hours into days. Parasite multiplication factor is calculated by ${{r}_{P}}_{pop}\mathrm{e}\mathrm{x}\mathrm{p}\left({{\delta}_{P}}_{pop}{a}_{L}\right)(1{f}_{pop})$ where the term $\mathrm{e}\mathrm{x}\mathrm{p}\left({{\delta}_{P}}_{pop}{a}_{L}\right)(1{f}_{pop})$ 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 humantomosquito 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, T_{1}, T_{2}, C, P_{1} and P_{2}) 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 ${k}_{T}=1/2.79$, ${q}_{c}=51.4$, ${V}_{c}=804$, ${q}_{1}=2020$, ${V}_{1}=3010$, ${q}_{2}=149$ and ${V}_{2}=13300$, which are the populationlevel mean estimates provided by Thanaporn Wattanakul and Joel Tarning (MahidolOxford Tropical Medicine Research Unit, Bangkok). We also assume the parameters are constrained in the optimization procedure by lower bounds of ${k}_{T}=1/3.25$, ${q}_{c}=40.6$, ${V}_{c}=469$, ${q}_{1}=746$, ${V}_{1}=2344$, ${q}_{2}=117$ and ${V}_{2}=10350$ and upper bounds of ${k}_{T}=1/2.32$, ${q}_{c}=160$, ${V}_{c}=1342$, ${q}_{1}=5034$, ${V}_{1}=3986$, ${q}_{2}=190$ and ${V}_{2}=17546$. Note that ${k}_{T}$ 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 ${q}_{c}$ 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 'bestfit' parameter values are the starting values given above with some adjustments on ${q}_{c}$ (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.
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).

Melbourne figshareData and computer code for model fitting and simulation.https://doi.org/10.26188/5cde4c26c8201
References

Likelihood based approaches to handling data below the quantification limit using NONMEM VIJournal of Pharmacokinetics and Pharmacodynamics 35:401–421.https://doi.org/10.1007/s1092800890944

A combination of new screening assays for prioritization of transmissionblocking antimalarials reveals distinct dynamics of marketed and experimental drugsJournal of Antimicrobial Chemotherapy 70:1357–1366.https://doi.org/10.1093/jac/dkv003

A mechanistic model quantifies artemisinininduced parasite growth retardation in bloodstage Plasmodium falciparum infectionJournal of Theoretical Biology 430:117–127.https://doi.org/10.1016/j.jtbi.2017.07.017

The power of malaria vaccine trials using controlled human malaria infectionPLOS Computational Biology 13:e1005255.https://doi.org/10.1371/journal.pcbi.1005255

A controlled human malaria infection model enabling evaluation of transmissionblocking interventionsJournal of Clinical Investigation 128:1551–1562.https://doi.org/10.1172/JCI98012

Modelling the transition of asexual blood stages of Plasmodium falciparum to gametocytesJournal of Theoretical Biology 202:113–127.https://doi.org/10.1006/jtbi.1999.1041

Plasmodium falciparum in vitro continuous culture conditions: A comparison of parasite susceptibility and tolerance to antimalarial drugs throughout the asexual intraerythrocytic life cycleInternational Journal for Parasitology: Drugs and Drug Resistance 7:295–302.https://doi.org/10.1016/j.ijpddr.2017.07.001

Genesis, sequestration and survival of Plasmodium falciparum gametocytes: parameter estimates from fitting a model to malariatherapy dataTransactions of the Royal Society of Tropical Medicine and Hygiene 95:497–501.https://doi.org/10.1016/S00359203(01)900161

Sexual development in plasmodium parasites: knowing when it's time to commitNature Reviews Microbiology 13:573–587.https://doi.org/10.1038/nrmicro3519

Withinhost modeling of bloodstage malariaImmunological Reviews 285:168–193.https://doi.org/10.1111/imr.12697

Generating random correlation matrices based on vines and extended onion methodJournal of Multivariate Analysis 100:1989–2001.https://doi.org/10.1016/j.jmva.2009.04.008

Hierarchical bayesian inference for concurrent model fitting and comparison for group studiesPLOS Computational Biology 15:e1007043.https://doi.org/10.1371/journal.pcbi.1007043

SoftwareR: A language and environment for statistical computingR Found Stat Comput, Vienna, Austria.

Stan Modeling Language Users Guide and Reference ManualStan Modeling Language Users Guide and Reference Manual.
Article and author information
Author details
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).
Version history
 Received: June 5, 2019
 Accepted: October 15, 2019
 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,324
 views

 320
 downloads

 34
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
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)
Further reading

 Computational and Systems Biology
 Medicine
Erectile dysfunction (ED) affects a significant proportion of men aged 40–70 and is caused by cavernous tissue dysfunction. Presently, the most common treatment for ED is phosphodiesterase 5 inhibitors; however, this is less effective in patients with severe vascular disease such as diabetic ED. Therefore, there is a need for development of new treatment, which requires a better understanding of the cavernous microenvironment and cellcell communications under diabetic condition. Pericytes are vital in penile erection; however, their dysfunction due to diabetes remains unclear. In this study, we performed singlecell RNA sequencing to understand the cellular landscape of cavernous tissues and cell typespecific transcriptional changes in diabetic ED. We found a decreased expression of genes associated with collagen or extracellular matrix organization and angiogenesis in diabetic fibroblasts, chondrocytes, myofibroblasts, valverelated lymphatic endothelial cells, and pericytes. Moreover, the newly identified pericytespecific marker, Limb BudHeart (Lbh), in mouse and human cavernous tissues, clearly distinguishing pericytes from smooth muscle cells. Cellcell interaction analysis revealed that pericytes are involved in angiogenesis, adhesion, and migration by communicating with other cell types in the corpus cavernosum; however, these interactions were highly reduced under diabetic conditions. Lbh expression is low in diabetic pericytes, and overexpression of LBH prevents erectile function by regulating neurovascular regeneration. Furthermore, the LBHinteracting proteins (Crystallin Alpha B and Vimentin) were identified in mouse cavernous pericytes through LCMS/MS analysis, indicating that their interactions were critical for maintaining pericyte function. Thus, our study reveals novel targets and insights into the pathogenesis of ED in patients with diabetes.

 Computational and Systems Biology
Largescale microbiome studies are progressively utilizing multiomics designs, which include the collection of microbiome samples together with host genomics and metabolomics data. Despite the increasing number of data sources, there remains a bottleneck in understanding the relationships between different data modalities due to the limited number of statistical and computational methods for analyzing such data. Furthermore, little is known about the portability of general methods to the metagenomic setting and few specialized techniques have been developed. In this review, we summarize and implement some of the commonly used methods. We apply these methods to real data sets where shotgun metagenomic sequencing and metabolomics data are available for microbiome multiomics data integration analysis. We compare results across methods, highlight strengths and limitations of each, and discuss areas where statistical and computational innovation is needed.