A model for malaria treatment evaluation in the presence of multiple species

To capture the transmission dynamics of both P. falciparum and P. vivax, we use a stochastic agent-based approach for the human population coupled with a deterministic system of ODEs for the mosquito population. Each human agent has their status with respect to both P. falciparum and P. vivax tracked over time, which allows mixed infections to be captured. The agent-based model is implemented by holding rates constant over discrete time-steps for computational efficiency and for ease of coupling to the ODEs that govern the mosquito population.

Each individual’s state is bivariate to specify their state with respect to P. falciparum and to P. vivax. For each Plasmodium species, humans are regarded as being susceptible ($S$), infectious with clinical symptoms ($I$), infectious but asymptomatic ($A$), recovered with no hypnozoites ($R$), recovered with hypnozoites ($L$ for latent: not applicable for P. falciparum), undergoing blood-stage treatment with no radical cure ($T$), or undergoing treatment with radical cure ($G$). Radical cure is defined as low-dose primaquine (3.5mg/kg total) administered over 14 days. A simplified schematic of the human transitions are depicted in Figure 1. Given the large number of connections between states required to describe the transmission and treatment dynamics, the model schematic uses a single line between connectors where multiple transitions apply and does not depict interactions between the species. The full list of possible transition rates and stoichiometries is provided in the Supplementary Table S1.

The dynamics of a human individual infected with type $x$ malaria are briefly described here for $x=f$ (P. falciparum) and $v$ (P. vivax). Individuals susceptible ($S$) to type $x$ malaria are infected at rate $\lambda_{x}$. Upon infection they develop clinical symptoms ($I$) with probability $p_{c,x}$ or are otherwise asymptomatic ($A$). Individuals are symptomatic for a mean duration of $1/\sigma_{x}$, at which point they either become asymptomatic ($A$) or die without treatment with probability $p_{I,x}$. The individuals with clinical symptoms may be treated at rate $c_{x}\tau_{x}$, where $c_{x}$ is the probability that they are able to access healthcare and $\tau_{x}$ is the rate at which medical attention is sought if it is readily available. Individuals with asymptomatic malaria will clear all blood-stage parasites at rate $\alpha_{x}$ and, for P. vivax, will be left with hypnozoites with probability $p_{h,v}$. When an individual with P. vivax seeks treatment they are prescribed radical cure ($G$) with probability $p_{N,x}$, otherwise they receive blood-stage treatment ($T$).

Any infectious individual may additionally be treated at rate $\eta_{x}(t)$ via an intervention program (such as MDA): the form of $\eta_{x}(t)$ will be discussed in Section 2.3. When an individual is treated this way they are prescribed radical cure with probability $p_{M,x}$, otherwise they receive blood-stage treatment. An individual undergoes treatment for an average of $1/\psi$ days (14-day primaquine) at which point they may: die with probability $p_{G,x}$, remain with asymptomatic blood-stage malaria with probability $p_{TfP}$, if $x=v$ they are left with latent hypnozoites ($L$) with probability $p_{P,v}$, otherwise they recover ($R$). Similarly, an individual ends blood-stage treatment after an average of $1/\rho_{x}$ days (3-day ACT) at which point they may: die with probability $p_{T}$, remain with asymptomatic blood-stage malaria with probability $p_{TfA}$, if $x=v$ they are left with latent hypnozoites with probability $p_{A,v}$, otherwise they recover. Latent stage P. vivax infected individuals experience a relapse at rate $\nu_{v}$ or they are reinfected at rate $\lambda_{v}r_{v}$ (where $r_{x}$ represents a possible reduction in susceptibility due to anti-parasite immunity) [24]. Upon relapse or reinfection from $L$ the individual experiences clinical malaria with probability $p_{L,v}$ (where $p_{L,v}<p_{c,v}$). Recovered individuals ($R$) are reinfected with rate $\lambda_{x}r_{x}$ at which point they become a clinical case with probability $p_{R,x}$ (where $p_{R,v}<p_{c,v}$). In addition to the possibility of relapse or reinfection, recovered and latent individuals lose immunity and hypnozoites at rates $\omega_{x}$ and $\kappa_{v}$, respectively, and return to being susceptible.

For individuals with mixed infections, there are several transitions in the model which depend on the individual’s state with respect to both species: we refer to these dependencies as species “interactions”. When an individual with a mixed infection is treated, they move from states in $\{I,A,L\}$ to a state in $\{T,G\}$ for both species (depending on the treatment); this is referred to as “treatment entanglement”. Similarly, when a patient stops treatment with respect to one malaria species, they are moved to one of the post-treatment states with respect to the other species. Death with respect to one species will cause a transition to death with respect to the other. The model allows treatment efficacies to vary for mixed infections; however, in this work we have assumed antimalarial efficacy against each species to be equivalent to the efficacy against mono-infections. We model P. vivax relapses triggered by the recovery of P. falciparum (“triggering”) by setting the relapse rate for a person that is recovered ($R$) from P. falciparum  but has latent stage ($L$) P. vivax to be $\hat{\nu}_{fv}=z_{f}\nu_{v}$, where $z_{f}>1$. The model also allows blood-stage P. vivax to be masked by blood-stage P. falciparum (“masking”) by treating a mixed infection as though it were a P. falciparum infection only with probability $h_{v}$. Explicitly, for an individual with P. falciparum and P. vivax both in states $I$ or $A$, the probability of receiving radical cure is $p_{N,fv}=h_{v}p_{N,f}+(1-h_{v})p_{N,v}$.

The dynamics of the mosquito population are governed by a system of ODEs (presented in the Supplementary §2). The mosquitoes follow standard Susceptible-Exposed-Infectious (SEI) dynamics with the addition of a seasonally varying death rate and the ability for mosquitoes to carry and spread mixed infection in a single bite (known as simultaneous inoculation). Asymptomatic individuals tend to have a lower peripheral parasitaemia and therefore were assumed to be less infectious to mosquitoes than individuals with clinical malaria with a relative infectiousness of $\zeta_{A,x}=0.1$.

We consider model parameters that would be representative of a Cambodia-like context, where both P. falciparum and P. vivax circulate, and the mosquito species present (An. dirus, An. minimus, An. maculatus, and An. barbirostris) are able to transmit both parasite species [25] allowing mosquitoes to be modelled as a single population that spreads both P. falciparum and P. vivax.

We simulate three treatment scenarios: current practice, accelerated radical cure, and unified radical cure. In each scenario we assume that when a person tests positive for malaria the species is always identified correctly for mono-infections, since specialised RDTs have been shown to have high sensitivity and specificity, particularly for P. falciparum (see, for example, [12]). The three treatment practices considered are:

1. 1.

   Current practice: Under this scenario, P. falciparum and most P. vivax cases are prescribed a blood-stage treatment but only 16% of P. vivax cases are prescribed radical cure. The low coverage of radical cure was chosen to match the rates reported from a study in Cambodia where radical cure was prescribed conservatively to 16% of detected P. vivax cases [26]. That is, the probability that an individual receives radical cure when being treated for malaria $x$, is

   $$p_{N,x}=\begin{cases}0,&\text{for }x=f\\ 0.16,&\text{for }x=v\\ (1-h_{v})0.16,&\text{for }x=fv,\end{cases}$$

   where $h_{v}$ is the probability that P. vivax is masked by P. falciparum when the individual is tested.

2. 2.

   Accelerated radical cure: Under this scenario, any eligible person who is diagnosed with P. vivax and returns a G6PD normal RDT is prescribed radical cure with a low-dose 14-day course of primaquine (total dose 3.5 mg/kg) alongside a 3-day course of blood-stage treatment. Anyone who is $>$6 months old and is not pregnant or lactating is considered eligible for radical cure. We assume that the G6PD RDTs have a sensitivity of 94% and a specificity of 91% [27], 6% of the population have G6PD enzyme activity $<$30%, 2% of the population are pregnant and/or lactating and 2% of the population are $<$6 months old. When combined, this means that 18% of people diagnosed with P. vivax will be ineligible to receive radical cure (further details are given Supplementary §3). The probability of receiving radical cure under this scenario is

   $$p_{N,x}=\begin{cases}0,&\text{for }x=f\\ 0.82,&\text{for }x=v\\ (1-h_{v})0.82,&\text{for }x=fv.\end{cases}$$

3. 3.

   Unified radical cure: Under this scenario, radical cure is prescribed to any eligible person in whom peripheral parasitaemia is detected, with eligibility as defined and calculated in the accelerated radical cure scenario. The probability of receiving radical cure in this scenario is

   $$p_{N,x}=\begin{cases}0.82,&\text{for }x=f\\ 0.82,&\text{for }x=v\\ 0.82,&\text{for }x=fv.\end{cases}$$

The three treatment scenarios and the probability of receiving radical cure is summarised in Table 1 with an assumed probability of masking of $h_{v}=0.5$.

For each of the three treatment scenarios we also consider the impact of an MDA intervention where a proportion of the population are given a blood-stage treatment irrespective of infective status, thus allowing a proportion of clinical and asymptomatic blood-stage infections to be treated. We assume that a proportion, $p$, of the population are prescribed a blood-stage treatment over a fixed period of time, $\Delta t=t_{2}-t_{1}$, so that the treatment rate of an individual with species $x$ due to MDA is:

We assume that people are not screened for G6PD status nor prescribed radical cure during MDA, based on concerns about haemolytic risks outweighing the benefits in patients who do not have malaria [28]. That is, the probability that an individual receives radical cure under MDA, given that they are treated for malaria type $x$ is $p_{M,x}$= 0 for all $x$.

We present the impact of different treatment and intervention strategies on the number of malaria cases and deaths from 2021 to 2030, where 2030 is the regional target for malaria elimination.

For each treatment and intervention scenario we run 50 model simulations and record the model compartments over time. Given the relatively short time frame, we ignore background human demographic dynamics in our model to reduce computational complexity.

For the current practice scenarios we assume that radical cure is prescribed conservatively and does not increase the risk of haemolysis (as a best case scenario). To account for heamolytic risk in the accelerated radical cure and unified radical cure scenarios we increase the probability of death by $5\times 10^{-6}$ for patients administered primaquine and the probability of radical cure failure by $5\times 10^{-4}$ (details on these values are given in Supplementary §3).

For the MDA intervention, we assume that half of the population receive blood-stage treatment over a 30 day period. The MDA is rolled out twice yearly, before and after the annual peak in transmission.

Initial conditions are set to resemble an endemic setting in eastern Cambodia, because the prevalence of clinical and asymptomatic infections in the region is well documented [29] and levels of immunity in eastern Cambodia were studied in a 2005 sero-survey [30]. Some parameter values were based on expert elicitation or a limited evidence base, and some parameter estimates vary greatly between studies. As such, the scenarios presented here are indicative of population dynamics of multi-species infections over time that accommodates interaction between P. falciparum and P. vivax infections and expected trends in the impacts of different interventions. Accordingly, the scenarios here should not be interpreted as forecasts of malaria cases and deaths in Cambodia and other similar malaria-endemic regions over the next 10 years. All parameters and initial conditions are given in Supplementary Tables S2 and S3.

We performed a sensitivity analysis on the model, in which each model parameter was modified separately and the relative change in model outputs recorded. To implement the sensitivity analysis, we considered the baseline value of each parameter (given in Tables 2 and 3) and ran simulations with the parameter scaled down to 80% and up to 120% while all other parameters remained fixed. If scaling a probability parameter up to 120% compared to baseline led to a probability being greater than 1 the value was instead held at 1. Similarly, the relative susceptibility of partially-immune individuals compared to susceptible individuals was not scaled up, so as not to exceed 1. For each set of parameter values, 50 repeats of the simulation were conducted and various outputs were recorded, including: the total number of P. falciparum infections, P. vivax infections, clinical P. falciparum infections, clinical P. vivax infections, P. vivax relapses, deaths, blood-stage treatments administered and radical cure treatments administered.

Figure 2 depicts the total number of infectious individuals in the population over time (the median, minimum and maximum of the 50 simulations) and figures showing all model compartments through time are presented in Supplementary Figures S2 and S3. The different treatment strategies had little effect on the prevalence of P. falciparum, with the unified treatment scenario resulting in a marginally lower prevalence of P. falciparum. For P. vivax  increasing the coverage of radical cure has a large impact on the prevalence of P. vivax, which appears to be approaching elimination towards the end of the time period. Meanwhile, the MDA intervention greatly reduced the prevalence of P. falciparum but had less of an effect on P. vivax. No scenarios led to elimination of malaria over the ten year period by 2030, although elimination may have been achieved over a longer time-frame.

Figure 3 shows boxplots of the cumulative number of infections and deaths by species over the 10 year period. The unified treatment strategy with G6PD testing of all individuals resulted in fewer infections and fewer deaths due to P. falciparum and P. vivax, despite the increased risk of haemolysis from radical cure. The accelerated radical cure approach approximately halved P. vivax infections but resulted in a small increase in P. falciparum infections, this increase is likely due to a reduction in the prevalence of mixed infections, which are more likely to be detected than mono-infections, inadvertently causing a reduction in the detection rate of P. falciparum.

In Figure 4, the results of the sensitivity analysis are presented in terms of the ten most influential parameters on cumulative infections (including clinical and asymptomatic cases). Sensitivity analyses with respect to all parameters and other outcomes are given in Supplementary Figure S4. These figures present the mean, minimum and maximum relative outcome compared to the baseline over 50 simulations, for each parameter set and orders them based on their relative sensitivity (in terms of absolute difference between the 80% and 120% scenarios).

In Figure 4 we see that P. falciparum and P. vivax infections were most sensitive to many of the vector-related parameters, including: the bite rate ($b$), the death rate of mosquitoes ($\delta_{0}$), the probability of transmission given an infectious bite (from humans to mosquitoes and vice versa, $\epsilon_{M,x}$ and $\epsilon_{H,x}$) and the rate at which exposed mosquitoes become infectious ($\gamma_{x}$). The spread of mosquito-borne infectious diseases are well known to be sensitive to these parameters [31].

Aside from the vector-related parameters, we identified several important human-related parameters, including: the relative infectiousness of asymptomatic carriers ($\zeta_{A,x}$), the relative susceptibility of partially-immune individuals ($r$), the rate at which asymptomatic infections are cleared ($\alpha_{x}$), the rate of treatment seeking ($\tau_{x}$) and accessibility of treatment ($c$). The parameters $\zeta_{A,x}$ and $\alpha_{x}$ determine the expected number of secondary infections generated by asymptomatic individuals. The parameter $r$ is related to anti-parasite immunity, it represents a possible lower rate of infection in recovered individuals. The parameters $\tau_{x}$ and $c$ determine the rate at which clinical cases seek treatment and the probability that they receive treatment for malaria.

In addition to the parameters that were influential on cumulative infections of both species, P. vivax infections were sensitive to the probability that hypnozoites are cleared when blood-stage parasites are cleared without treatment ($p_{h}$). This emphasises the contribution of relapses in the overall P. vivax malaria burden.

Other model outcomes (clinical infections, total deaths, total relapses, total treatments received) were broadly sensitive to the same model parameters without any additions or omissions (see Supplement Figure 4 for details). Small differences in relative sensitivity between measuring cumulative infections versus clinical infections were largely centred on the parameters having to do with the proportions immune expected to develop clinical malaria upon reinfection (that is, $p_{R}$ for P. falciparum and $p_{c}$ for P. vivax).

ACREME funded the salary of RIH, and contributed to the costs of data cleaning and organisation by PN and the CNM.