Heterogeneous mosquito exposure increases Plasmodium vivax and Plasmodium falciparum co-infections: a modelling study

(a) Data review

We performed a literature review on PubMed to collect data from field epidemiology studies recording P. vivax, P. falciparum and co-infection prevalence, among the general population. Specifically, we defined those studies as cross-sectional and active infection-detection studies performed in communities, unrestricted by age, sex or symptomatic status. The keywords we used were

Title: (Malaria OR Vivax OR Falciparum) NOT (Health-care facility OR Hospital OR Health centre OR Systematic Review OR Meta*)

Abstract: Prevalence OR Distribution OR Epidemiology OR Cross-section

Text: Vivax AND Falciparum

After title and abstract screening with Rayyan [39], we screened the main text and extracted the P. vivax, P. falciparum and co-infection prevalence from each study, as well as the sample size and diagnostic method used to generate prevalence. We excluded studies that failed to report all those metrics or that did not screen the general population. A spreadsheet of the details of all the studies retained after title and abstract screening is available as a supplementary file. We first evaluated the correlation between P. vivax and P. falciparum prevalence using a Pearson correlation test weighted by study size and evaluated at a significance level $\alpha =0.05$. Then, we defined different metrics related to the coupled dynamics of P. vivax, P. falciparum and co-infections. First, we introduced the notion of expected co-infection prevalence under independence of P. vivax and P. falciparum. That is, we assumed that if P. vivax and P. falciparum transmission are independent, co-infection prevalence should be the product of P. falciparum and P. vivax prevalence. We then defined the proportional excess co-infection prevalence $e$ as the difference between observed and expected co-infection prevalence, ${\displaystyle p_{\mathrm{F}\mathrm{V}}^{\mathrm{O}}}$and ${\displaystyle p_{{\mathrm{F}\mathrm{V}}^{\mathrm{E}}}}$, respectively, normalized by different measures of local transmission intensity: P. vivax prevalence (i), P. falciparum prevalence (ii) or overall malaria prevalence (iii):

where ${\displaystyle p_{\mathrm{V}}^{\mathrm{O}}\ge p_{\mathrm{F}\mathrm{V}}^{\mathrm{O}}}$ and ${\displaystyle p_{\mathrm{F}}^{\mathrm{O}}\ge p_{\mathrm{F}\mathrm{V}}^{\mathrm{O}}}$ are the observed P. vivax and P. falciparum prevalence, respectively.

(b) Transmission model

To capture the dynamics of P. vivax and P. falciparum in humans and in mosquitos, we used two coupled deterministic systems of ordinary differential equations: one for humans and one for mosquitos. For the human system, we considered three binary states of infection: infection with P. falciparum, blood-stage infection with P. vivax and liver-stage infection with P. vivax hypnozoites. All combinations of these states are possible, leading to eight compartments in the model. Infection by a mosquito infected with P. vivax causes humans to move to an infected compartment for P. vivax blood and liver stages at rate ${\displaystyle {\lambda }_{{\mathrm{h}}_{\mathrm{V}}}}$. Infection by a mosquito infected with P. falciparum causes humans to move into an infected compartment for P. falciparum at rate ${\displaystyle {\lambda }_{{\mathrm{h}}_{\mathrm{F}}}}$. Infection by a mosquito infected with both species causes humans to move into the compartment positive for all three infection states at rate ${\displaystyle {\lambda }_{{\mathrm{h}}_{\mathrm{F}\mathrm{V}}}}$. Humans with blood-stage infections recover at rate ${\displaystyle r_{\mathrm{F}}}$ or ${\displaystyle r_{\mathrm{V}}}$, respectively. Hypnozoites activate at rate $f$ and die at rate $\gamma$.

The mosquito system follows a susceptible–exposed–infected dynamic. We considered all combinations of three infection states (S, E and I) for the two species, leading to nine compartments in the model. Infection by feeding on a human infected with blood-stage P. vivax causes mosquitos to move into an exposed compartment for P. vivax at rate ${\displaystyle {\lambda }_{{\mathrm{m}}_{\mathrm{V}}}}$. Infection by feeding on a human infected with P. falciparum causes mosquitos to move into an exposed compartment for P. falciparum at rate ${\displaystyle {\lambda }_{{\mathrm{m}}_{\mathrm{F}}}}$. Infection by feeding on a human infected with blood-stage parasites of both species causes mosquitos to move in the exposed compartment for both at rate ${\displaystyle {\lambda }_{{\mathrm{m}}_{\mathrm{F}\mathrm{V}}}}$. Mosquitos move from exposed to infected at rate ${\displaystyle {\eta }_{\mathrm{V}}}$ or ${\displaystyle {\eta }_{\mathrm{F}}}$. Two compartments account for being simultaneously infected with P. vivax and exposed to P. falciparum or vice-versa. Mosquitos are born and die at rate $\mu$.

Equations for the systems and graphical representation of the compartments and dynamics are shown in the electronic supplementary material, figures S1 and S2.

(c) Heterogeneity

We stratified the human compartments into nine groups accounting for different levels of exposure to infectious mosquito bites. We assumed that the number of bites per day, $X$, are log-normally distributed with mean $1$ and variance ${\sigma }^2$: ${\displaystyle X\sim \mathrm{l}\mathrm{o}\mathrm{g}N\left(1,{\sigma }^2\right)}$. Log-bites, $Z=\log \left(X\right)$, are thus normally distributed, $Z\sim N\left(1,{\sigma }^2\right)$ (electronic supplementary material, figure S3). The larger the variance, the more bites cluster in on smaller parts of the population (electronic supplementary material, figure S4). We discretize the normal distribution using the Gauss–Hermite quadrature approach, described elsewhere [23]. This procedure outputs weights (${\displaystyle w_i}$) used to set the size of the nine strata of the compartmental model, and abscissae ($x_i=\exp \left(z_i\right)$) to scale the force of infection. The subscript $i\in 1,2,\dots 9$ is a population stratum index. The force of infection for P. vivax, P. falciparum and co-inoculation on humans was adjusted with the scaling factors $x_i$ as follows:

The force of infection on mosquitos was adjusted as follows:

Symbol and parameter descriptions for the six equations above can be found in electronic supplementary material, table S1.

(d) Model fitting

We fit the model to the polymerase chain reaction (PCR) data points identified in the literature review to generate parameter estimates for ${\sigma }^2$, $m$ and $\phi$ using ABC. The parameter ${\sigma }^2$ is the variance of the heterogeneity distribution, $m$ is the mosquito density and $\phi$ is the relative vector competence (i.e. the relative ability of a mosquito to transmit one species versus the other). All other parameters had fixed values (electronic supplementary material, table S1), and we tested the dependence of our conclusions on those values with a sensitivity analysis (electronic supplementary material, figure S5). We sampled 100 000 draws from the following prior distributions:

We plugged each of the 100 000 independent draws into the model and simulated P. vivax, P. falciparum and co-infection prevalence by running the model to equilibrium, which is shown in electronic supplementary material, figure S6. For each of the $i$ PCR-data points in parallel, we added measurement error to each of the ${\displaystyle j=1,\dots ,100000}$ prevalence trios by drawing a simulated data point drawn from $Multinom(n_i,p_j)$ where $n_i$ is the sample size used to generate the PCR data point and $p_j$ is the simulated prevalence trio. Then, we computed the Euclidian distance between the PCR data point and each simulated data point on the logit scale and selected its 100 nearest neighbours. We finally assessed how well simulated data points from the fitted heterogeneous model recreate the distribution of the true set of PCR data points, compared with a model without heterogeneity (${\sigma }^2=0$), by plotting an illustrative set of simulated data points from the ensemble of nearest neighbours. The parameter values that generated the $i$ sets of 100 nearest neighbours were combined into a single set of parameter trios, which can be seen as a numerical approximation of the ensemble of posteriors. Parameter estimates used for the following simulations were based on the marginal 25%, 50% and 75% quantiles of the numerical approximation of the posterior ensemble.

(e) Scenario simulations

We performed scenario simulations for three possible epidemiological settings: low transmission with high heterogeneity (25% and 75% quantiles of $m$ and ${\sigma }^2$, respectively), medium transmission with medium heterogeneity (50% quantiles of both ${\displaystyle m\mathrm{a}\mathrm{n}\mathrm{d}{\sigma }^2}$) and high transmission with low heterogeneity (75% and 25% quantiles of $m$ and ${\sigma }^2$, respectively). We chose a single value for $\phi$, the median of the numerical approximation of the posterior ensemble.

(f) Hypnozoite activation-by-fever

We explored the hypothesis that P. falciparum-induced fevers could activate P. vivax hypnozoites. In addition to the spontaneous P. vivax relapse rate $f$ originating from hypnozoite-infected compartments, we introduced an additional relapse rate triggered by P. falciparum fevers. A proportion ${\displaystyle {\pi }_{\mathrm{F}}{\pi }_{\mathrm{T}}}$ of the force of infection ${\displaystyle {\lambda }_{\mathrm{F}}}$ exiting from hypnozoite-infected compartments is deviated towards the compartment infected with hypnozoites, and blood-stage P. falciparum and blood-stage P. vivax, while the remainder ${\displaystyle 1-\left({\pi }_{\mathrm{F}}{\pi }_{\mathrm{T}}\right)}$ continues as usual towards the compartment infected with hypnozoites and P. falciparum only. These dynamics are illustrated in the electronic supplementary material, figure S7. The parameter ${\displaystyle {\pi }_{\mathrm{F}}}$ represents the proportion of P. falciparum infections that are febrile, and ${\displaystyle {\pi }_{\mathrm{T}}}$ represents the hypothetical proportion of P. falciparum infections that would trigger a P. vivax relapse. Given the unknown nature of the parameter ${\displaystyle {\pi }_{\mathrm{T}}}$, we performed a sensitivity analysis by setting the parameter to 0.1, 0.5 or 1. As seen in electronic supplementary material, table S1, we set ${\displaystyle {\pi }_{\mathrm{F}}=0.5}$.

(g) Treatment

We implemented treatment through a partial instantaneous blocking of a given force of infection. A proportion $\pi \theta$ of the force of infection is deviated towards a susceptible compartment, and the residual proportion $1-\pi \theta$ continues towards an infected compartment as usual. The parameter $\pi$ is the proportion of infections that are febrile, and the parameter $\theta$ is a composite parameter, encompassing the various factors that govern successful treatment (health care access, adherence, treatment effectiveness, etc.). We differentiated between ${\displaystyle {\pi }_{\mathrm{F}}}$ and ${\displaystyle {\pi }_{\mathrm{V}}}$ as the febrile proportions for P. falciparum and P. vivax, respectively. To simulate the standard treatment guidelines, we considered that a new P. vivax infection would be treated with radical cure, clearing hypnozoites and blood-stage parasites, and a new P. falciparum infection would be treated with blood-stage-only treatment. When a new infection is treated, any concurrent infection would also be cured due to the treatment. For a new co-inoculation, a proportion ${\displaystyle {\pi }_{\mathrm{V}}\theta }$ is blocked with radical cure, and ${\displaystyle {\pi }_{\mathrm{F}}\left(1-{\pi }_{\mathrm{V}}\right)\theta }$ with only blood-stage-only treatment. The treatment dynamics are illustrated in the electronic supplementary material, figure S8. As seen in electronic supplementary material, table S1, we set ${\displaystyle {\pi }_{\mathrm{F}}=0.5}$, ${\displaystyle {\pi }_{\mathrm{V}}=0.1}$ and $\theta =0.25$.

(a) Data review

We collected data from field epidemiology studies recording population prevalence of P. vivax, P. falciparum and co-infection. Our search terms identified 759 articles on PubMed. After abstract screening, we obtained 125 articles, among which 37 reported three-way prevalence. From the 37 articles, we extracted 104 data points (each a prevalence trio). The 104 data points included 58 data points that were generated using PCR, which is more precise than other methods such as microscopy or rapid diagnostic tests. When establishing the presence of a species by microscopy, one is likely to stop the test once a first species is observed thereby discarding the potential presence of a second or third species [6]. This bias would likely lead to a strong underestimation of co-infection prevalence. We therefore present these data in the electronic supplementary material, figure S9 but exclude them from the main analyses.

There is a positive linear relationship between P. vivax and P. falciparum PCR-prevalence (figure 1a). Figure 1b illustrates the difference between co-infection prevalence expected under independence, defined as the product of P. vivax and P. falciparum prevalence, and observed co-infection prevalence. Most of the data points exist above or at the line of equality, signifying more co-infections than expected under independence. To explore whether that excess is related to transmission intensity, we computed the proportional excess in co-infection prevalence and plotted it against species-specific prevalence (figure 1c,1d). Departure from independence varies more when prevalence is lower.

• Download figure
• Open in new tab
• Download PowerPoint

(b) Model fitting and validation

We estimated the values of three model parameters ${\sigma }^2$, $m$ and $\phi$, by fitting our model of co-infection with heterogeneity to the 58 PCR data points. The marginal 25%, 50% and 75% quantiles, respectively, of the parameters used to generate the combined 5 800 nearest-neighbour simulated points are 0.67, 1.53 and 2.4 (${\sigma }^2)$; 0.27, 0.34 and 0.63 ($m$); and 0.85, 1.01 and 1.15 ($\phi$). While for most data points, the pairwise marginal posteriors are well behaved (electronic supplementary material, figure S11a,b), for some points where any of the three prevalence values approaches 0, a range of parameter pairs generate the same prevalence triplet, indicating issues of identifiability (electronic supplementary material, figure S11c).

We evaluated our model’s capacity to generate simulated data that resemble real data by plotting 58 simulated points drawn uniformly at random from the 5 800 nearest neighbours (figure 2). This illustrative set of 58 simulated data points resembles the distribution of the true set of 58 PCR data points, indicating that our simple model of co-infection with heterogeneity captures departure from independence in observed co-infection prevalence. Draws from the model without heterogeneity (${\sigma }^2 = 0$) show that some excess co-infections are still expected but not to the extent observed in the data (electronic supplementary material, figure S12).

• Download figure
• Open in new tab
• Download PowerPoint

(c) Scenario simulations

We used the aforementioned parameter estimates to explore the activation-by-fever hypothesis and the radical cure proposition under three epidemiologically plausible conditions. The simulations for all other combinations of the quartiles of $m$ and ${\sigma }^2$ are shown in the electronic supplementary material, figures S13, S14.

In figure 3, we explore the epidemiological consequences of triggering of P. vivax relapses by febrile P. falciparum infections in the absence of treatment. Allowing P. falciparum-induced fever to trigger relapses has the greatest effect in the setting with high transmission and low heterogeneity (figure 3c). In this setting, allowing all P. falciparum-induced fevers to activate hypnozoites (${\displaystyle {\pi }_{\mathrm{T}}=1}$), the scenario most favourable to the hypothesis, would increase P. vivax prevalence from 39.9% to 42.0%, that is by at most 5.3% and co-infection prevalence from 20.7% to 22.3%, that is by 7.5%. While this model does not provide evidence on the likelihood of febrile P. falciparum triggering P. vivax relapses, we can conclude that if it does occur, then it causes very limited increases in P. vivax transmission at the population level.

• Download figure
• Open in new tab
• Download PowerPoint

In all settings, P. falciparum treatment reduced P. vivax prevalence, and P. vivax treatment reduced P. falciparum prevalence (figure 4a–f,m–o). Higher levels of heterogeneity were associated with more frequent P. vivax and P. falciparum co-infection, and a greater effect of treating one parasite on the other. For instance, P. falciparum treatment reduced P. vivax prevalence by 49.6% versus 2.9% in the high versus low heterogeneity setting, respectively (figure 4m versus 4o); meanwhile, P. vivax treatment reduced P. falciparum prevalence by 53.9% versus 6.2% in the high versus low heterogeneity setting (figure 4m versus 4o).

• Download figure
• Open in new tab
• Download PowerPoint

Treating P. vivax with radical cure and P. falciparum with blood-stage only treatment leads to a 77.0%, 26.1% and 17.8% relative reduction in P. vivax prevalence in the low, medium and high transmission settings, respectively (figure 4g–i,m–o). Compared with that, treating both P. vivax and P. falciparum with radical cure leads to a 77.8%, 48.8% and 26.1% relative reduction in P. vivax prevalence in the low, medium and high transmission settings, respectively (figure 4j–l,m–o).