Modelling the population dynamics of Plasmodium falciparum gametocytes in humans during malaria infection

Every year over two hundred million people are infected with the malaria parasite. Renewed efforts to eliminate malaria has highlighted the potential to interrupt transmission from humans to mosquitoes which is mediated through the gametocytes. Reliable prediction of transmission requires an improved understanding of in vivo kinetics of gametocytes. Here we study the population dynamics of Plasmodium falciparum gametocytes in human hosts by establishing a framework which incorporates improved measurements of parasitaemia in humans, a novel mathematical model of gametocyte dynamics, and model validation using a Bayesian hierarchical inference method. We found that the novel mathematical model provides an excellent fit to the available clinical data from 17 volunteers infected with P. falciparum, and reliably predicts observed gametocyte levels. We estimated the P. falciparum’s sexual commitment rate and gametocyte sequestration time in humans to be 0.54% (95% credible interval: 0.30-1.00) per life cycle and 8.39 (6.54-10.59) days respectively. Furthermore, we used the data-calibrated model to predict the effects of those gametocyte dynamics parameters on human-to-mosquito transmissibility, providing a method to link within-human host kinetics of malaria infection to epidemiological-scale infection and transmission patterns.


Abstract 22
Every year over two hundred million people are infected with the malaria parasite. Renewed 23 efforts to eliminate malaria has highlighted the potential to interrupt transmission from humans 24 to mosquitoes which is mediated through the gametocytes. Reliable prediction of transmission 25 requires an improved understanding of in vivo kinetics of gametocytes. Here we study the 26 population dynamics of Plasmodium falciparum gametocytes in human hosts by establishing a 27 framework which incorporates improved measurements of parasitaemia in humans, a novel 28 mathematical model of gametocyte dynamics, and model validation using a Bayesian 29 hierarchical inference method. We found that the novel mathematical model provides an 30 excellent fit to the available clinical data from 17 volunteers infected with P. falciparum, and 31 reliably predicts observed gametocyte levels. We estimated the P. falciparum's sexual 32 commitment rate and gametocyte sequestration time in humans to be 0.54% (95% credible 33 interval: 0.30-1.00) per life cycle and 8.39 (6.54-10.59) days respectively. Furthermore, we used 34 the data-calibrated model to predict the effects of those gametocyte dynamics parameters on 35 human-to-mosquito transmissibility, providing a method to link within-human host kinetics of 36 malaria infection to epidemiological-scale infection and transmission patterns. 37 38 data very well (median of posterior predictions vs. observed data). Furthermore, the narrow 95% 97 PI indicates a low level of uncertainty in predicted total parasitaemia. 98 gametocytaemia data (circles). We emphasise that this was not a fitting exercise, rather an 109 independent validation of the calibrated model. as 202, 301 and 302, the 95% PI (red shaded area) extends into the detectable region again after 133 the asexual parasitaemia reaches below the detection limit, indicating that there was a small 134 chance that the patients may have suffered a recrudescence during the observation period (of 135 course, they did not) or after the observation period (although this predication cannot be 136 evaluated because artemisinin combination therapy was given immediately after the period). 137 Despite the above discrepant observations for asexual parasitaemia, the model predictions of 138 gametocytaemia are very consistent with the data for all 17 volunteers, as shown in Figure 3. 139

Estimation of gametocyte population dynamics parameters 140
The model calibration process provided the joint posterior density for the model parameters, 141 which were used to estimate some key biological parameters governing the population dynamics 142 of P. falciparum gametocytes (see the Material and Methods for details). As shown in Table 1, 143 the sexual commitment rate -the percentage of asexual parasites entering sexual development 144 during each asexual life cycle -is estimated to be 0.54% (95% credible interval (CI): 0.30%, 145 1.00%). This in vivo estimate is much lower than 11% that was estimated from in vitro data 11 . 146 The gametocyte sequestration time is the average time that immature gametocytes (stages I-IV) 147 cannot be observed in the peripheral circulation and was estimated to be 8.39 days (95% CI: 148 6.54, 10.59 days). The estimate for the circulating gametocyte lifespan is 63.5 days, with a broad 149 95% CI (12.7 to 1513.9 days) resulting from the long-tailed posterior distribution (see Figure S1, 150 Supplementary Information) and is much longer than the previous in vitro estimate of 16-32 151 days 12 (note that our lower bound of the 95% CI is lower than the in vitro estimated range). The 152 wide estimate for the circulating gametocyte lifespan, and in particular the high upper bound of 153 the 95% CI, is primarily due to the limited observation time in the VIS as explored in more detail 154 in the Discussion. 155 Table 1: Estimates of some key biological parameters and comparison with the literature. 156 The estimates of the biological parameters (middle column) are shown as the median and 95% 157 credible interval (CI) of the marginal posterior parameter distribution (see Figure S1 in the maximum estimate]" or simply "a low estimate -a high estimate".

163
As shown in Table 1, there are both similarities and differences in parameter estimates for P. 164 falciparum based on our analysis and the analysis of historical neurosyphilis patient data 9 . We 165 found that they exhibited very similar in vivo sexual commitment rate (median 0.54% vs. mean 166 0.64% with overlapping plausible ranges) and gametocyte sequestration time (median 8.39 days 167 vs. mean 7.4 days with overlapping plausible ranges). On the other hand, we found that the 168 circulating gametocyte lifespan for P. falciparum was much longer than that estimated from the 169 neurosyphilis patient data 9 (median 63.5 days vs. 6.4 days with very different ranges). 170 Finally, we provided an estimate for the parasite multiplication factor, an important parameter 171 that quantifies the net replication rate of asexual parasites and also influences the rate of Predicting the impact of gametocyte kinetics on human-to-mosquito 177 transmissibility 178 Having validated our mathematical model of asexual parasitaemia and gametocyte dynamics, we 179 were able to predict how the gametocyte dynamics parameters influence the transmissibility of P. 180 falciparum malaria from humans to mosquitoes in various epidemiological scenarios. In 181 particular, we focused on the early phase of infection where the innate immune response is 182 minimal and treatment has not been administered (in order to avoid complications that our 183 mathematical model was not designed to capture). Two scenarios were considered: 184 • Predicting the potential infectiousness of newly hospitalised malaria patients for various 185 values of sexual commitment rate and gametocyte sequestration time. In the model, 186 gametocytaemia was assumed to be a surrogate of the potential infectiousness and newly 187 hospitalised patients were assumed to exhibit a total parasitaemia of approximately 10 8 188 parasites/mL based on published clinical data from Cambodia and Thailand 18 . As illustrated 189 in Figure 4A, we simulated the model (for different sexual commitment rates and 190 gametocyte sequestration times) and looked at the critical gametocytaemia level (indicated 191 by Gc) corresponding to the time when the total parasitaemia (wave-like black curves) first 192 reached 10 8 parasites/mL (their associations are indicated by the dotted lines and arrows). 193 • Predicting the non-infectious period of malaria patients for various values of sexual 194 commitment rate and gametocyte sequestration time. In the model, the non-infectious period 195 was defined to be time from the inoculation of infected red blood cells to the time when the 196 gametocytaemia reach 10 3 parasites/mL, which is a threshold allowing for human-to-197 mosquito transmission 7 . Note that this non-infectious period does not include the latent 198 period due to the liver stage which should be considered if the starting time were taken as 199 the inoculation by mosquito. As illustrated in Figure 4B, we simulated the model (for Materials and Methods for details). The red curves indicate the cases corresponding to 217 gametocytaemia of 10 3 parasites/mL. 218 Figure 4C shows that a higher sexual commitment rate or a lower gametocyte sequestration time 219 leads to a higher gametocytaemia (Gc) at the time of hospitalisation. The red curve in Figure 3C  220 indicates the level curve of 10 3 gametocytes/mL (i.e. the threshold for infectiousness as 221 mentioned above) dividing the heatmap into two regions. To the left, Gc is below 10 3 222 gametocytes/mL, suggesting clinical presentation precedes infectiousness, while to the right Gc is 223 above 10 3 gametocytes/mL and the converse applies. The Gc value obtained by model simulation 224 using the median estimates of the population mean parameters (indicated by the black dot) is 225 below 10 3 gametocytes/mL, suggesting that newly hospitalized malaria patients are less likely to 226 be infectious, and thus efforts to identify and treat infections in a timely manner may have a 227 substantial impact in terms of reduced transmission potential. 228 Figure 4D reinforces the result in Figure 4C using the non-infectious period (tc). As the sexual 229 commitment rate increases or the gametocyte sequestration time decreases, tc decreases. 230 However, for relatively large values of the sexual commitment rate (e.g. > 20%), we observed an 231 increase in tc as the sexual commitment rate increases (see the top-right corner of Figure 4D). 232 This is because an increased fraction of sexual conversion can lead to both a decrease in the rate 233 of asexual parasite growth and an increase in the number of sexual parasites, and the impact of 234 the former more than counterbalances that of the latter. For gametocyte kinetic parameters, we found that our in vivo estimate of the P. falciparum 247 sexual commitment rate was similar to that found in the neurosyphilis patient data 9 but was much 248 smaller than previous in vitro estimates (see Table 1). Novel VIS data using biomarkers specific 249 to early sexual parasites (e.g. PfGEXP5 19 ) may help to validate and improve our in vivo estimate. 250 Furthermore, our in vivo estimate of the circulating gametocyte lifespan is substantially larger 251 than previous in vitro estimates and is also larger than that found in the neurosyphilis patient 252 data 9 (see Table 1). One way to improve our in vivo estimate would be to use new P. falciparum 253 data with gametocytaemia measurements over a longer period of time to capture the natural 254 decay of circulating gametocytes, which was not observed in the current VIS data (see Figure 3). 255 We also predicted the effects of altered gametocyte kinetic parameters on the transmissibility 256 from human to mosquito, focusing on two scenarios: the infectiousness of newly hospitalised 257

Study population and measurements 286
The data used in this modelling study are from a previously published VIS 7 where 17 malaria-287 naïve volunteers were inoculated with P. falciparum-infected red blood cells (3D7 strain). The  Parasitaemia in the volunteers were monitored at least daily since the inoculation. 295 The data contain the measurements of total parasitaemia (total circulating asexual parasites and 296 gametocytes per mL blood), asexual parasitaemia (circulating asexual parasites per mL blood), 297 gametocytaemia (circulating female and male gametocytes per mL blood) and plasma 298 concentration of PQP, at multiple time points after inoculation (see Figure 1-3 and Figure S2 in 299 the Supplementary Information). Further details about the VIS are given in 7 . We emphasize that 300 the data used in model fitting is the total parasitaemia (from the first measurement to the time before any treatment other than PQP) and the other data (i.e. asexual parasitaemia and 302 gametocytaemia) are used to validate the model. 303   Research Unit, Bangkok), is a three-compartment disposition model with two transit 373 compartments for absorption (see schematic diagram in Figure S3, Supplementary Information). 374

Gametocyte dynamics model
Based on Figure S3, the model is formulated to be a system of ordinary differential equations: 375 The model parameters were obtained by fitting the model output ( ) to the PQP concentration 382 data from the VIS in MATLAB (version 2016b; The MathWorks, Natick, MA) using lsqcurvefit. 383 The MATLAB code (with detailed comments) are publicly available at 384 https://doi.org/10.26188/5cde4c26c8201. Note that for volunteers whose number of data points is 385 less than the number of parameters in the PK model (such that optimization fails), the parameter 386 values are given by the estimates from another VIS (provided by Thanaporn Wattanakul and Joel 387 Tarning). The best-fit curve and associated parameter values for all volunteers are provided in 388 Figure S2 and Table S1 in the Supplementary Information. 389 Fitting the model to parasitaemia data 390 We took a Bayesian hierarchical modelling approach 27 to fit the gametocyte dynamics model to 391 the data from all 17 volunteers. In detail, each volunteer has 12 model parameters (i.e. those in 392 Table 2 except % and $ ; also called the individual parameters) and lower and upper bounds of 393 the parameters are given in Table 2. If denoting the individual parameters to be )*B and lower 394 and upper bounds to be $ and x respectively, the following transformations are used to convert 395 the bounded individual parameters to unbounded ones (denoted by )*B ) in order to in order to 396 improve computational efficiency 28,29 : 397 Note that Ω~•~= ~•~q~•~ where q is the correlation matrix. The prior distributions for 407 the 12 population mean parameters ~•~ are given by uniform distributions with bounds given in 408 Table 2. The prior distribution for the 12 population SD parameters is standard half-normal and 409 the prior distribution for the lower Cholesky factor of the correlation matrix is given by 410 Cholesky LKJ correlation distribution with shape parameter of 2 29,30 . The distribution of the 411 observed parasitaemia measurements is assumed to be a lognormal distribution with mean given 412 by the model-simulated values and SD parameter with prior distribution of a half-Cauchy 413 distribution with a location parameter of zero and a scale parameter of 5. The distribution for the 414 observed parasitaemia measurements was used to calculate the likelihood function and the M3 415 method 31 was used to penalise the likelihood for data points below the limit of detection for the 416 total parasitaemia (10 parasites/mL 7 ).   The posterior prediction and 95% prediction interval (PI) are given by the median and quantiles 431 of 2.5% and 97.5% of the 5000 model outputs at each time respectively (see Figure 1-3 for 432 example). 433 The estimates of some key biological parameters (Table 1)  The gametocyte dynamics model with parameters given by the median estimates of the 443 population mean parameters was used to simulate the two scenarios predicting the dependence of 444 human-to-mosquito transmissibility on the sexual commitment rate and gametocyte sequestration 445 time ( Figure 4). 446 Final analysis and visualization were performed in MATLAB. All computer codes (R, Stan, 447 MATLAB), data and fitting results (CSV and MAT files) and an instruction document (note that 448 reading the document first will make the code much easy to follow) are publicly available at 449 https://doi.org/10.26188/5cde4c26c8201. 450