Modeling the impact of Plasmodium falciparum sexual stage immunity on the composition and dynamics of the human infectious reservoir for malaria in natural settings

Malaria transmission remains high in Sub-Saharan Africa despite large-scale implementation of malaria control interventions. A comprehensive understanding of the transmissibility of infections to mosquitoes may guide the design of more effective transmission reducing strategies. The impact of P. falciparum sexual stage immunity on the infectious reservoir for malaria has never been studied in natural settings. Repeated measurements were carried out at start-wet, peak-wet and dry season, and provided data on antibody responses against gametocyte/gamete antigens Pfs48/45 and Pfs230 as anti-gametocyte immunity. Data on high and low-density infections and their infectiousness to anopheline mosquitoes were obtained using quantitative molecular methods and mosquito feeding assays, respectively. An event-driven model for P. falciparum sexual stage immunity was developed and fit to data using an agent based malaria model infrastructure. We found that Pfs48/45 and Pfs230 antibody densities increased with increasing concurrent gametocyte densities; associated with 55–70% reduction in oocyst intensity and achieved up to 44% reduction in proportions of infected mosquitoes. We showed that P. falciparum sexual stage immunity significantly reduces transmission of microscopic (p < 0.001) but not submicroscopic (p = 0.937) gametocyte infections to mosquitoes and that incorporating sexual stage immunity into mathematical models had a considerable impact on the contribution of different age groups to the infectious reservoir of malaria. Human antibody responses to gametocyte antigens are likely to be dependent on recent and concurrent high-density gametocyte exposure and have a pronounced impact on the likelihood of onward transmission of microscopic gametocyte densities compared to low density infections. Our mathematical simulations indicate that anti-gametocyte immunity is an important factor for predicting and understanding the composition and dynamics of the human infectious reservoir for malaria.

The method analyzes the change in antibody titers at the population level between the peak-wet (time point i) and end-wet (time point j) transmission season assuming alterations of antibody kinetics by frequent boosting (e.g. reinfection, persistent gametocyte carriage) unlikely. The time-lag between sampling time points was estimated to be 111 days as from the last sampling day at the peak-wet season (August 23 th ) to the last sampling day at the end-wet season (December 12 th ). For Pfs230, the antibody half-life was estimated 46.4 days (Median OD = 0.897 for ni= 104 samples, Median OD = 0.258 for nj = 24 samples). For Pfs48/45, the antibody half-life was estimated 142.0 days (Median OD = 0.8315 for ni= 41, Median OD = 0.4837 for nj = 10 samples).

Section 2.1 Model initialization
P. falciparum sexual stage immune model was initialized using simulated gametocyte data ( Figure 1) as previously configured [2] and antibody half-lives as determined by a statistical model (section 1 -equation 1). The dynamical sexual stage immune model can be described as follow: (2) where: Equation (2) represents the antibody titer at any time (day) of life (ti). A(ti) = Antibody level at day ti. r = log (2)/antibody half-life = exponential antibody decay rate. β(γ(tᵢ), d(tᵢ)) = Boost size at day ti. γ(ti) = gametocyte density at day ti. d(ti)= Duration of exposure in days at day ti. A0 = Initial antibody level from maternal immunity at day 0.

Section 2.2 Model fitting to estimate best anti-Pfs230 and anti-Pfs48/45 antibody half-lives.
P. falciparum sexual stage immune model was then fitted to Pfs48/45 and Pfs230 antibody data under the following structural assumptions: higher gametocyte density will boost more strongly than lower density -two categories were considered here (submicroscopic and microscopic gametocyte density) -and age-dependent exposure should impact the magnitude of the boost and the inter-individual variation within a smoothness constraint. Parameters were estimated by using observed features of antibody titers distributions: median titers and inter-quartile ranges by age, season, and current gametocyte density.
Multiple parameters sets were run to estimate the best fit allowing some stochasticity within a 95% confidence interval such that: αi and A(ti) = Antibody titers of individuals at day ti from the field data and simulated data respectively, whose durations of exposure (d(ti)) fall into the same bin. ɛ (set to 0.1) = absolute difference (margin of error) between field data and simulated data when comparing means of antibody titers of individuals whose durations of exposure (d(ti)) fall into the same bin.

Section 2.2.1. Fitting outputs and interpretation
In simulation b (Table 1, Fig. 2), the model is initialized and run with Pfs230 antibody half-life equals to 46.4 days as determined in section 1 -equation (1). In simulation b, the model significantly diverges from the data likely due to slow antibody decay and / or overestimation of boost sizes. In simulation c, the model is run as in b, except that boost size is reduced (half-lowered), which seems not to better fit the model. In simulation d, the model is run as in c, except that the half-life of the antibody response is extended to 56 days. Note that antibody half-life was extended to 56 days by simply doubling our best fit antibody half-life of 28 days (simulation e). This part of the process aimed to test the impact of extended antibody half-life against lowered boost sizes. At first glance, an antibody half-life of 56 days well fits the model at the population level ( Fig. 2 -2d), except that the resolution in antibody response (median and inter-quartile ranges) between microscopic and submicroscopic gametocyte densities as observed in the field data is no picked up by the model (Fig. 2 -1a versus Fig. 2 -3d). In simulation e, the model seems to best fit the data balancing boost size and antibody half-life to account for fine resolution in immune boost and decay.
Two advantages were taken from the fitting procedure: i) Best estimates in boost sizes and antibody half-lives (28 days for both Pfs230 and Pfs48/45) fit the present data (S1 Fig) and data from and independent cohort study (26 and 32 days for Pfs230 and Pfs48/45 respectively, S1 Fig), ii) The present dynamical model estimates demonstrates that half-life of sexual stage antibodies can be as short as estimated at the population level (Equation (1), at least for anti-Pfs230 antibody half-life) and as reported in previous studies [3][4][5][6].  Table 2 and 3) multiplied by a factor of 0.5.
| where δ is set to be ≤ ε.

Figure 2: Fitting outputs (Mean OD of antibody) under different parameters sets in comparison to field data (a).
Median antibody titers from the study field data are shown in 1a and 2a. Simulations (b, c, d and e) are run with parameters sets (Table 1) to fit the data. Outputs 1, 2, 3 and 4 from each of simulations b, c, d and e represent: Antibody response dynamics at the individual level (1), mean antibody response at the population level (2), median antibody response per age at the population level under two gametocyte density categories (3) and median antibody response per season at the population level (4).

Section 3. Oocyst simulation and effect of P. falciparum sexual stage immunity on transmission.
Regression coefficients (β values) that associate with anti-Pfs48/45 and anti-Pfs230 antibodies were transformed into odd ratios ( ) and the effect of the corresponding antibody on oocyst intensity defined as: where TRAoocyst is defined as transmission reducing activity on oocyst intensity. The following TRA function: = (3) where Y = residual mean oocyst intensity X = observed mean oocyst intensity, = Adjusted effect of immunity on oocyst intensity, was incorporated into the EMOD malaria model to assess the effect of sexual stage immunity on simulated oocyst intensity. Simulated oocyst intensity (X) per individual was sampled based on gametocyte density carriage using a kernel density estimation. Oocyst data from which are drawn the model's oocyst count (Fig 3) were obtained from individuals seronegative for both anti-Pfs48/45 and anti-Pfs230 antibodies in the field data and stratified by gametocyte density (submicroscopic, <16, 17-100, ≥100 gametocytes/µL). Oocyst sampling was simulated 40 times in a closed loop mathematical draw to mimic a sample size of at least 20 fed mosquitoes required for the membrane feeding assay in the field. The arithmetic mean of oocyst intensity was then estimated per simulated individual. Residual oocyst intensity of a seropositive individual was obtained using equation (3) and the corresponding mosquito's infection rate was inferred using a polynomial interpolation over a reference curve [7] of the oocyst intensity versus oocyst prevalence as illustrated in Fig 4. The reduction in mosquito's infection rate (TRAmosquito) achieved in seropositive individuals was estimated using the following conventional TRA equation [7]: Distribution of simulated mosquito's oocysts against measured field data. N = 15000 observations from a birth cohort of 100 simulated individuals (over 50 years' simulation was performed to store data across all ages as in the field data) with oocyst intensity collected at different time points (start-wet, peak-wet and dry season). Field data consisted of oocyst intensity collected at different time points (start-wet, peak-wet and dry season) from a sample size of N = 291 blood samples fed to mosquitoes. Figure 4: Illustration of the relationship between oocyst prevalence and intensity in standard membrane feeds data. Data extracted from [7].

Section 4. Accounting for gametocyte density measurement uncertainty by QT-NASBA
Because we quantified gametocyte density by QT-NASBA, it appeared relevant that our simulations account for gametocyte measurement uncertainty by the QT-NASBA method. QT-NASBA is a RNA to DNA based pathogen detection method. The specificity of the assay [8] depends on primers that are specific for the RNA target. The reaction relies on enzymes AMV-RT, Rnase H and T7 RNA polymerase. An isothermal condition combined to T7 ensures that primers only anneal to single stranded target RNA and product cDNA preventing amplification of genomic DNA. Moreover, QT-NASBA uses a real-time molecular beacon technology combined to time to positivity (TTP) measurement as a proxy and a calibration curve of 10 0 -10 6 gametocytes/µL dilution series to determine the parasite density in the test ample. Despite its high sensitivity and specificity, the method is subject to uncertainties because of random errors of analyzers (e.g. precision of TTP) or users (e.g. parallax type error). By averaging random error data extracted from a series of calibration curves [9], we established an uncertainty range within which our model's true gametocyte densities were randomly sampled, thus incorporating measurement uncertainty. Fig. 5 shows the 95% Confidence Interval of measurement uncertainty which increases with decreasing true gametocyte density. When comparing to survey data quantified by QT-NASBA, the simulated true gametocyte densities were transformed by a log-normal smearing and fixed detection threshold (Fig. 5) to account for the method uncertainty described in [9]. . Relationship between gametocyte density measurement uncertainty and simulated true gametocyte density. N = 15000 observations from a birth cohort of 100 simulated individuals with data on gametocyte density collected at different time points (start-wet, peak-wet and dry season) over 50 years' simulation.