Optimal interruption of P. vivax malaria transmission using mass drug administration

TL;DR

This study develops a multiscale mathematical model to evaluate the effectiveness of mass drug administration targeting hypnozoites in controlling Plasmodium vivax malaria, highlighting the temporary impact and challenges in achieving elimination.

Contribution

The paper introduces a novel multiscale integro-differential equation model for P. vivax hypnozoite dynamics and optimizes MDA timing considering seasonality and initial prevalence.

Findings

MDA effects are temporary and depend on initial prevalence.

Radical cure alone may not achieve elimination.

Optimal MDA intervals vary with seasonality and prevalence.

Abstract

\textit{Plasmodium vivax} is the most geographically widespread malaria-causing parasite resulting in significant associated global morbidity and mortality. One of the factors driving this widespread phenomenon is the ability of the parasites to remain dormant in the liver. Known as hypnozoites, they reside in the liver following an initial exposure, before activating later to cause further infections, referred to as relapses. As around 79-96$\%$ of infections are attributed to relapses, we expect it will be highly impactful to apply treatment to target the hypnozoite reservoir to eliminate \textit{P. vivax}. Treatment with a radical cure to target the hypnozoite reservoir is a potential tool to control or eliminate \textit{P. vivax}. We have developed a multiscale mathematical model as a system of integro-differential equations that captures the complex dynamics of \textit{P. vivax} hypnozoites and the effect of hypnozoite relapse on disease transmission. Here, we use our model to study the anticipated effect of radical cure treatment administered via a mass drug administration (MDA) program. We implement multiple rounds of MDA with a fixed interval between rounds, starting from different steady-state disease prevalences. We then construct an optimisation model to obtain the optimal MDA interval. We also incorporate mosquito seasonality in our model to study its effect on the optimal treatment regime. We find that the effect of MDA interventions is temporary and depends on the pre-intervention disease prevalence (and choice of model parameters) as well as the number of MDA rounds under consideration. We find radical cure alone may not be enough to lead to \textit{P. vivax} elimination under our mathematical model (and choice of model parameters) since the prevalence of infection eventually returns to pre-MDA levels.

Tables

Table 1. Table 1 : Definitions, values and sources for model parameters.

Symbol	Definition	Value/s	Source
$a$	Biting rate of mosquitoes	80 year-1	[18]
$b$	Transmission probability: mosquito to human	0.5	[51]
$c$	Transmission probability: human to mosquito	0.23	[8]
$b_m\left(t\right)$	Mosquito birth rate (seasonal)	Time-varying
$g$	Baseline mosquito birth rate	0.1 day-1	[19]
$\eta$	Seasonal amplitude	0.1	Assumed
$\varphi$	Seasonal phase	0	Assumed
$m$	Number of mosquitoes per human	Varried
$n$	Rate of mosquito sporogony	1/12 days-1	[19]
$\gamma$	Blood-stage infection clearance rate	1/60 day-1 [1/156 1/46]	[13]
$\alpha$	Hypnozoite activation rate	1/332 day-1 [1/700 1/110]	[57]
$\mu$	Hypnozoite death rate	1/425 day-1 [1/750 1/155]	[57]
$\nu$	Average number of hypnozoites per mosquito bite	8.5 [1.5 20.5]	[57]
$\lambda \left(t\right)$	Force of reinfection	Calculated	$\lambda \left(t\right)=m_{0}ab{I}_m$
$p_1\left(t\right)$	Probability that a blood-stage infected individual has no hypnozoites and $MOI=1$	Time-varying	Calculated
$p_2\left(t\right)$	Probability that a blood-stage infected individual has hypnozoites and $MOI=1$	Time-varying	Calculated
$k_i\left(t\right)$	Probability that liver-stage infected individual has $i$ hypnozoites within liver	Time varying	Calculated, [37]
$k_T\left(t\right)$	Average number of hypnozoites within liver for liver-stage infected individuals	Time varying	Calculated, [37]
$p_{blood}$	Probability that ongoing blood-stage infections are cleared instantaneously	0.9	Assumed
$p_{rad}$	Probability that hypnozoites dies instantaneously	0.9	Assumed
$D_b\left(t\right)$	Clearance rate of blood-stage parasite	Calculated
$D_l\left(t\right)$	Clearance rate of liver-stage parasite (hypnozoite)	Calculated
$N$	Total number of MDA rounds	varied
$s_j\left(j=1,2,\dots ,N\right)$	MDA administration time	Estimated with optimisation
$x_j\left(j=1,2,\dots ,N-1\right)$	MDA intervals	Estimated with optimisation

Equations

$$b_m(t)=b_m(0)\left(1+\eta \cos \left(\frac{2\pi t}{365}+\varphi \right)\right),$$

$$\frac{\mathrm{d}S}{\mathrm{d}t}=-\lambda S+\mu k_1(t)L+p_1(t)\gamma I+D_l(l)L+D_b(t)p(t)I,$$

$$\displaystyle\frac{\mathrm{d}I}{\mathrm{d}t}=\lambda(S+I)+\alpha k_{T}(t)L-\gamma\big{(}p_{1}(t)+p_{2}(t)\big{)}I-D_{b}(t)I,$$

$$\displaystyle\frac{\mathrm{d}L}{\mathrm{d}t}=-\lambda L-\mu k_{1}(t)L-\alpha k_{T}(t)L+\gamma p_{2}(t)I-D_{l}(t)L+D_{b}(t)\big{(}1-p(t)\big{)}I,$$

$$\frac{\mathrm{d}{S}_m}{\mathrm{d}t}=b_m(t)-acI{S}_m-b_m(t)S_m,$$

$$\frac{\mathrm{d}{E}_m}{\mathrm{d}t}=acI{S}_m-\left(b_m(t)+n\right)E_m,$$

$$\frac{\mathrm{d}{I}_m}{\mathrm{d}t}=n{E}_m-b_m(t)I_m,$$

$$\lambda =$$

$$D_b(t)=$$

$$D_l(t)=$$

$$=$$

$$p_H^r(t,s_1,\dots ,s_N),p_A^r(t,s_1,\dots ,s_N),p_C^r(t,s_1,\dots ,s_N),p_D^r(t,s_1,\dots ,s_N)$$

$$\frac{d{p}_H^r}{dt}=$$

$$\frac{d{p}_A^r}{dt}=$$

$$\frac{d{p}_C^r}{dt}=$$

$$\frac{d{p}_D^r}{dt}=$$

$$p_H^r(t,s_1,s_2,s_3,\dots ,s_N)=$$

$$p_A^r(t,s_1,s_2,s_3,\dots ,s_N)=$$

$$\displaystyle+(1-p_{\text{rad}})^{N}\big{(}p_{A}(t)-e^{-\gamma(t-s_{N})}p_{A}(s_{N})\big{)},$$

$$p_H(t)=$$

$$p_A(t)=$$

$$\mathbf{N}(t)=(N_H(t),N_A(t),N_C(t),N_D(t),N_P(t),N_{PC})$$

$$G(t,z_H,z_A,z_C,z_D,z_P,z_{PC}):$$

$$=\text{exp}\left\{-q(t)+{\int }_0^t\frac{\lambda (\tau )\left(z_P{e}^{-\gamma (t-\tau )}+(1-e^{-\gamma (t-\tau )})z_{PC}\right)}{1+\nu \left(1-{\sum }_{f\in F}{z}_f{p}_f(t-\tau )\right)}d\tau \right\},$$

$$q(t)={\int }_0^t\lambda (\tau )d\tau .$$

$$G^{s_1,s_2,\dots s_N}(t,z_H,z_A,z_C,z_D,z_P,z_{PC}):=\mathbb{E}\left[\prod _{f\in F^{\prime }}{z}_f^{N_s^{s_1,s_2,\dots s_N}(t)}\right]$$

$$\displaystyle=\text{exp}\Bigg{\{}-q(t)+\int_{s_{N}}^{t}\lambda(\tau)\frac{e^{-\gamma(t-\tau)}z_{P}+(1-e^{-\gamma(t-\tau)})z_{PC}}{1+\nu\left(1-\sum_{f\in F}z_{f}.p_{f}(t-\tau)\right)}d\tau$$

$$+{\int }_0^{s_1}\lambda (\tau )\frac{(1-p_{\text{blood}})e^{-\gamma (t-\tau )}{z}_P+(1-(1-p_{\text{blood}})e^{-\gamma (t-\tau )})z_{PC}}{1+\nu \left(1-{\sum }_{f\in F}{z}_f.p_s^r(t-\tau ,s_1-\tau )\right)}d\tau$$

$$+{\int }_{s_1}^{s_2}\lambda (\tau )\frac{(1-p_{\text{blood}}{)}^2{e}^{-\gamma (t-\tau )}{z}_P+(1-(1-p_{\text{blood}}{)}^2{e}^{-\gamma (t-\tau )})z_{PC}}{1+\nu \left(1-{\sum }_{f\in F}{z}_f.p_s^r(t-\tau ,s_1-\tau ,s_2-\tau )\right)}d\tau$$

$$+{\int }_{s_2}^{s_3}\lambda (\tau )\frac{(1-p_{\text{blood}}{)}^3{e}^{-\gamma (t-\tau )}{z}_P+(1-(1-p_{\text{blood}}{)}^3{e}^{-\gamma (t-\tau )})z_{PC}}{1+\nu \left(1-{\sum }_{f\in F}{z}_f.p_s^r(t-\tau ,s_1-\tau ,s_2-\tau ,s_3-\tau )\right)}d\tau$$

$$\displaystyle+\ldots+\int_{s_{N-1}}^{s_{N}}\lambda(\tau)\frac{(1-p_{\text{blood}})^{N}e^{-\gamma(t-\tau)}z_{P}+(1-(1-p_{\text{blood}})^{N}e^{-\gamma(t-\tau)})z_{PC}}{1+\nu\left(1-\sum_{f\in F}z_{f}.p^{r}_{s}(t-\tau,s_{1}-\tau,s_{2}-\tau,s_{3}-\tau,\ldots,s_{N}-\tau)\right)}d\tau\Bigg{\}}.$$

$$p(t)$$

$$\displaystyle=\frac{P\big{(}N_{H}(t)=0)-P(N_{H}(t)=N_{A}(t)=N_{P}(t)=0\big{)}}{1-P\big{(}N_{A}(t)=N_{P}(t)=0\big{)}}.$$

$$P($$

$$\displaystyle=\begin{cases}\text{exp}\Bigl{\{}-q(t)+\int_{0}^{t}\frac{\lambda(\tau)}{1+\nu p_{H}(t-\tau)}d\tau\Bigr{\}}&\text{if}\ t<s_{1}\\ \text{exp}\Bigl{\{}-q(t)+\int_{s_{N}}^{t}\frac{\lambda(\tau)}{1+\nu p_{H}(t-\tau)}d\tau+\int_{0}^{s_{1}}\frac{\lambda(\tau)}{1+\nu p_{H}^{r}(t-\tau,s_{1}-\tau)}d\tau\\ \qquad+\int_{s_{1}}^{s_{2}}\frac{\lambda(\tau)}{1+\nu p_{H}^{r}(t-\tau,s_{1}-\tau,s_{2}-\tau)}d\tau+\ldots+\int_{s_{N-1}}^{s_{N}}\frac{\lambda(\tau)}{1+\nu p_{H}^{r}(t-\tau,s_{1}-\tau,\ldots,s_{N}-\tau)}d\tau\Bigr{\}}&\text{if}\ t\geq s_{N},\end{cases}$$

$$\displaystyle P\big{(}N_{A}(t)+N_{P}(t)=0\big{)}=$$

$$=$$

$$\displaystyle P\big{(}N_{H}(t)=$$

$$=$$

$$p_1(t)=$$

$$P(N_A(t)+N_P(t)=$$

Full text

Keywords: P. vivax dynamics, multi-scale model, superinfection, mass drug administration, radical cure

1 Introduction

Plasmodium vivax is a parasite that causes malaria, responsible for 4.5 million cases in 2020 [42]. After an infective mosquito bite, the P. vivax parasite triggers a primary infection and can remain dormant (known as a ‘hypnozoite’) within the human liver for a prolonged period before causing a secondary infection known as a ‘relapse’ [57, 2]. Because of the relapse characteristics, P. vivax has become the most globally widespread parasite and is responsible for significant morbidity and mortality [4, 7]. The reason for hypnozoite activation is still not clear, and the number of hypnozoites established per infective mosquito bite and the recurrence time vary geographically [45].

When it comes to P. vivax for control or elimination, the biological characteristics of P. vivax make it more challenging than other malaria parasites because P. vivax transmission can be re-established from hypnozoite activation [37, 45]. An estimated 79–96% of the total vivax cases are due to hypnozoite activation [46, 25, 1, 14]. Thus, targeting the hypnozoite reservoir with treatment is an important element of any program for P. vivax elimination [11]. Mass drug administration (MDA) is an effective intervention for controlling many diseases and was advocated by the World Health Organization (WHO) in the 1950s to control malaria transmission [24]. MDA involves treating the entire population, or a well-defined sub-population, in a geographic location regardless of their infection status [39, 24]. Most of the antimalarial drugs currently used to treat malaria only clear blood-stage parasites. Drugs that clear hypnozoites from the liver are referred to as ‘radical cure’, examples of which are primaquine and tafenoquine [55, 53, 44]. In a radical cure MDA intervention, individuals are given a combination of two such drugs: artemisinin combination therapy (ACT) for clearing blood-stage parasites and primaquine to clear hypnozoites. However, because of the risk of haemolysis in glucose 6 phosphate dehydrogenase (G6PD)–deficient individuals, radical cure is not recommended by the WHO without screening for G6PD deficiency [40, 23, 54].

The effect of radical cure treatment on P. vivax transmission has been explored in a number of mathematical models [26, 2, 12, 47, 57, 58, 56]. However, most of these model do not consider the variation in hypnozoite number per mosquito bite [26, 2, 12, 47, 58]. Mehra et al. [37] have developed a within-host model capturing hypnozoite dynamics and variation across infected mosquito bites that explicitly models the effect of radical cure treatment on the hypnozoite dynamics and reservoir. We have previously developed a multiscale model [5] by embedding Mehra et al.’s within-host model without treatment, which only uses three compartments at the population level while considering hypnozoite dynamics and the effect of the hypnozoite reservoir on disease transmission. The effect of three rounds of MDA with radical cure on P. vivax prevalence has been studied in a randomised controlled trial [43]. However, as the MDA implementations are expensive, and the empirical evidence remains unclear as to the overall impact they are expected to have, mathematical modelling is well suited to explore the overall impact and establish efficient designs before the actual implementation of the MDAs [29, 27]. The impact of multiple rounds ($>3$) of MDA on P. vivax transmission has not been explored with a mathematical model as far as we are aware. As P. vivax is transmitted by mosquitoes, overall disease transmission is greatly affected by the mosquito population distribution in a region, which in turn is influenced by climate factors [22, 10, 6, 17]. Thus, the effect of treatment can be influenced by an abundance of mosquitoes and, hence, by seasonality. However, while a few P. vivax transmission models have considered the role of seasonality in mosquito population dynamics [26, 12, 50, 36], few have also captured the rich dynamics of hypnozoites [36] and none have considered the impact of drug administration. Also, the abundance of mosquitoes and the contribution of hypnozoite activation can frequently trigger superinfection (reinfection of individuals that are already infected) which can potentially delay recovery from infection ([52, 16]). However, only a few P. vivax mathematical transmission models account for superinfection [57, 56, 50, 35, 36].

In this article, we study the impact of multiple rounds of radical cure treatment within an MDA program on disease transmission by incorporating hypnozoite dynamics into an epidemic transmission framework. We account for superinfection and consider the impact of seasonal mosquito population changes (which we refer to as “seasonality” throughout). In Section 2, we extend our existing multiscale model [5] to incorporate the effect of radical cure treatment and seasonality. We then obtain some key parameters for the population model from the within-host model [37] under multiple rounds of MDA and obtain the recovery rate under superinfection. In Section 3, we provide illustrative results for both the within-host scale and transmission setting. We construct an optimisation problem to determine the optimal interval between MDA rounds with and without accounting for seasonality before concluding remarks are presented in Section 4.

2 Methods

A multiscale mathematical model that accounts for hypnozoite variation within individuals without treatment has already been developed [5]. In our previous work, we did not account for superinfection in the population level model. Here, we extend the population level model to account for superinfection and allow treatment (via MDA) with a radical cure. The inclusion of superinfection in the model is important since for high transmission settings, overlapping blood-stage infections are frequent due to exposure to multiple infectious bites as well as the activation of hypnozoites.

2.1 Population transmission model with treatment

Let $S$, $I$ and $L$ represent the fraction of the human population who are susceptible with no hypnozoites, blood-stage infected and liver-stage infected, respectively. Individuals in both $S$ and $L$ compartments are susceptible to infective mosquito bites and become blood-stage infected ($I$) at the rate $\lambda(t)=mabI_{m}$, where $m$ is the human-to-mosquito ratio, $a$ is mosquito biting rate, and $b$ is the transmission probability from mosquito to human. Recovery without accounting for superinfection is straightforward. In our previous model [5], we did not account for superinfection in the population-level model while the within-host framework permits superinfection. Here, we consider superinfection at the population level by following the work of Mehra [35]. To do that we need to consider the multiplicity of infection (MOI), defined as the number of distinct parasites co-circulating within a blood-stage infected individual (for P. vivax, either from a new infectious bite or hypnozoite activation). When considering superinfection, an individual might experience multiple blood-stage infections, and recovery from the blood-stage infection is conditioned upon how many infections (MOI) they are currently experiencing. Those blood-stage infected individuals who are experiencing only one infection will recover and move out of $I$ and, depending on the hypnozoite reservoir size, either become susceptible ($S$) or liver-stage infected ($I$). Following the work of Mehra [35], we define two parameters $p_{1}(t)$ and $p_{2}(t)$ where $p_{1}(t)$ is the probability that a blood-stage infected individual only experiencing one infection (MOI=1) has no hypnozoites in their liver at time $t$ and $p_{2}(t)$ is the probability that a blood-stage infected individual only experiencing one infection (MOI=1) has hypnozoites in their liver at time $t$. Hence, after recovery from blood-stage infection, individuals become susceptible ($S$) at rate $p_{1}(t)\gamma$ and become liver-stage infected ($L$) at rate $p_{2}(t)\gamma$, where $\gamma$ is the natural recovery rate. Thus, the probability of staying blood-stage infected at time $t$ is \big{(}1-(p_{1}(t)+p_{2}(t))\big{)}.

Individuals suffer relapses from hypnozoite activation, the rate at which depends on the hypnozoite reservoir size and the baseline activation rate for each hypnozoite, $\alpha$. We define $k_{T}(t)$ to be the average hypnozoite reservoir size given liver-stage infected. That is $\alpha k_{T}(t)$ is the relapse rate. Individuals from the $L$ compartment become susceptible without experiencing a relapse if they have only one hypnozoite (with probability $k_{1}(t)$) and the hypnozoite dies naturally before activation, at rate $\mu$. For the mosquito population, we define $S_{m},\ E_{m}$, and $I_{m}$ to be the fraction of susceptible, exposed, and infectious mosquitoes, respectively. Susceptible mosquitoes become exposed when they take a blood meal from an infected individual at the rate $acI$, where $c$ is the transmission probability from human to mosquito. After the incubation period (mean $1/n$ days), they become infectious and can transmit parasites to humans. The time-dependent parameters $p_{1}(t),\ p_{2}(t),\ k_{1}(t)$, and $k_{T}(t)$ capture the dynamics of hypnozoites and are obtained from the within-host model (Section 2.2) as an integral function of the force of infection, $\lambda(t)$, which makes the model a system of integro-differential equations (IDEs). The definition and derivation of these time-dependent parameters are discussed in Section 2.2. The model schematic is depicted in Figure 1.

Suppose that drug treatment is administered successively at times $s_{1},s_{2},\ldots,s_{N}$, where $N$ is the total number of MDA rounds. The effect of blood-stage and liver-stage radical cure treatments are captured by the time-dependent functions $D_{b}(t)$ and $D_{l}(t)$, respectively. To implement the effect of radical cure in the population level model, we assume that radical cure has an instantaneous effect [37]. That is, on administration, all ongoing blood-stage infections are instantaneously cleared with probability $p_{blood}$, and each hypnozoite in the liver dies instantaneously with probability $p_{rad}$. Without treatment, blood-stage infections are cleared one at a time, but with treatment, blood-stage infections will all be cleared with probability $p_{blood}$. We define $p(t)$ as the probability of blood-stage infected individuals having no hypnozoites in their liver [5]. Therefore, at the time when radical cure is administered, individuals that were blood-stage infected either become susceptible with probability $p(t)$ or become liver-stage infected with probability $(1-p(t))$. Blood-stage infected individuals who are not cured following treatment undergo the same dynamics as those who receive no treatment. Liver-stage infected individuals whose hypnozoites have not been fully cleared following treatment will undergo the same dynamics as if without treatment but starting with the reduced hypnozoite reservoir. Hence, if radical cure is administered at time $t=s_{j}$, where $j$ is the number of MDA rounds, the drug has an effect only on the ongoing infections and hypnozoites established from time $t=s_{j-1}$ until $t=s_{j}$. Hypnozoites that are established after time $t=s_{j}$ or any blood-stage infections caused by either hypnozoite activation or infectious mosquito bites after $t=s_{1}$ will undergo dynamics as if without treatment (until the next time of MDA application). Since we are concerned with disease dynamics over a time scale of years, the assumption of an instantaneous effect of the radical cure is appropriate, as drugs such as artemisinin, which clears blood-stage parasites, have a half-life of 1.93 hours [9] and primaquine, which kills hypnozoites have a relatively short half-life of approximately 3.7–9.6 hours [28]. Another drug, tafenoquine, that also kills hypnozoites has a half-life of approximately 14–28 days which is short compared to a time scale of years [48]. Since the number of mosquitoes in the environment influences P. vivax dynamics dramatically [22, 10], it is important to account for seasonal environmental effects on the mosquito population (see, for example, [6, 17]). To incorporate mosquito seasonality, we consider that the mosquito birth rate at time t, $b_{m}(t)$, is regulated by a cosine function with a period of 1 year as follows:

[TABLE]

where $b_{m}(0)=g$ is the baseline mosquito birth rate, $\eta\in[0\ 1)$ is the seasonal amplitude and $\phi$ is the seasonal phase (taken to be 0). Note that if $b_{m}(t)=b_{m}(0)=g$, then the mosquito population is constant, that is, without seasonality. With all the assumptions outlined above, the system of IDEs that describe the dynamics is (see Appendix A for a detailed derivation of the model):

[TABLE]

where,

[TABLE]

is the force of reinfection, and $m_{0}=\frac{N_{m}(0)}{N_{h}}$ is the initial mosquito ratio. Here $D_{b}(t)$ and $D_{l}(t)$ are blood-stage parasite and liver-stage parasite (hypnozoite) clearance rates respectively and are given by:

[TABLE]

where $\delta_{D}(\cdot)$ is the Dirac delta function. That is, any blood-stage parasite will be instantaneously cleared with probability $p_{blood}$ every time the radical cure is administered [37] and depending on the parameter $p_{1}(t)$ which is the probability that an individual experiencing only one infection has no hypnozoites in the liver given blood-stage infection (Equation 2.2.2) and $p_{2}(t)$ which is the probability that an individual experiencing only one infection has hypnozoites in the liver given blood-stage infection (Equation 2.2.3), move to the susceptible compartment ($S$) at rate $p_{1}(t)D_{b}(t)$ and to the liver-stage infected compartment ($L$) at rate $p_{2}(t)D_{b}(t)$, respectively. As each hypnozoite is cleared with probability $p_{rad}$, the liver-stage clearance rate $D_{l}(t)$ depends on how many hypnozoites are present in the liver. That is, $D_{l}(t)$ depends on $k_{1}(t),\ k_{2}(t),\ \ldots,k_{i}(t)$, where $k_{i}(t)$ is the probability that a liver-stage infected individual has $i$ hypnozoites. All model parameters are defined in Table 1.

2.2 Within-host model with treatment

A within-host model considering the effect of radical cure on hypnozoite dynamics was introduced by Mehra et al. [37]. They developed the framework considering $N$ MDA rounds but explored analytically and numerically considering one MDA round. Here, we solve the necessary equations for $N$ MDA rounds. First, the dynamics of a single hypnozoite under treatment were modelled, then a fixed number of hypnozoites introduced by a single mosquito bite before accounting for continuous mosquito inoculation where each mosquito bite contributes an average of $\nu$ hypnozoites to the reservoir. The within-host model also assumes that radical cure has an instantaneous effect.

For the short-latency case (in which hypnozoites can immediately activate after establishment without going through a latency phase), a hypnozoite can be in one of four different states. Let $H$, $A$, $C$, and $D$ represent the state of establishment, activation, clearance and death for a single hypnozoite, respectively. Suppose that drug treatment is administered successively at times $s_{1},s_{2},\ldots,s_{N}$. We denote the state of the hypnozoite at time $t$ with $X_{r}(t,\ s_{1},\ s_{2},\ \ldots,s_{N})\in(H,A,C,D)$ with corresponding probability mass function (PMF)

[TABLE]

respectively. The governing equations for the state probabilities under treatment are given by Equations (17)–(22) from Mehra et al. [37]:

[TABLE]

where the parameters $\alpha,\ \gamma,$ and $\mu$ are as per Table 1. Since our population model in Equations (2)–(6) uses the parameters $p_{1}(t),\ p_{2}(t),\ k_{1}(t)$ and $k_{T}(t)$, we seek to obtain expressions for these parameters from the within-host model under multiple rounds of MDA. Evaluating the parameters $p_{1}(t),\ p_{2}(t),\ k_{1}(t)$ and $k_{T}(t)$ in the population model requires the probability of hypnozoite establishment ($p^{r}_{H}(t)$) and the probability of hypnozoite activation ($p^{r}_{A}(t)$) [5]; hence we solve Equations (7)–(8) for $N$ MDA rounds to give:

[TABLE]

where $p_{H}(t)$ and $p_{A}(t)$ are the probability of establishment and activation of a hypnozoite without treatment, respectively, and are given by:

[TABLE]

Figure 2 shows the effect of three rounds of MDA on the dynamics of a single hypnozoite. The probability of hypnozoite establishment ($p^{r}_{H}(t)$) and hypnozoite activation ($p^{r}_{A}(t)$) under 3 rounds of MDA (with $p_{blood}=p_{rad}=0.5$) are illustrated in Figure 2(A) and Figure 2(B), respectively. The drug is administered for the first time 200 days after the hypnozoite is established and the interval between each MDA round is fixed at 30 days.

We now define two additional states, $P$ and $PC$, to denote an ongoing primary infection from infective mosquito bites and a cleared primary infection, respectively. Let $N_{f}(t)$ denote the number of hypnozoites in states $f\in\{H,A,C,D\}:=F$ at time $t$ and $N_{P}(t),\ N_{PC}(t)$ denote the number of ongoing and cleared primary infections, respectively, at time $t$. Defining the state space $F^{\prime}:=\left\{H,A,C,D,P,PC\right\}$, the probability generating function (PGF) for

[TABLE]

with $\mathbf{N}(0)=\mathbf{0}$ can be written following from Equation (30) in Mehra et al. [37] (for short-latency case ($k=0$) with probability of a blood-stage infection after an infectious bite, $p_{prim}=1$) (by the law of total expectation):

[TABLE]

where $q(t)$ is the mean number of infective bites in the interval $(0,t]$ and is given by:

[TABLE]

All parameters are as per Table 1. The expression for the joint PGF with drug administration at time $t=s_{1}$ is given by Equation (31) in Mehra et al. [37]. Following a similar analysis, if the drug is administered at $N$ successive times ($s_{1},\ s_{2},\ldots,s_{N}$) then the joint PGF for the number of hypnozoites/infections in each state is:

[TABLE]

We now use the PGF in Equation (14) to derive expressions for the population-level parameters $p(t)$, $p_{1}(t),\ p_{2}(t),\ k_{1}(t)$, and $k_{T}(t)$ under multiple MDA rounds.

2.2.1 Probability blood-stage infected individual has no hypnozoites (under $N$ rounds of MDA)

With $p(t)$ defined as the probability that an individual has an empty hypnozoite reservoir conditional on an ongoing blood-stage infection (i.e. primary infection or relapse) from Equation (13) of Anwar et al. [5] we have:

[TABLE]

where the probability that an individual has an empty hypnozoite reservoir at time $t$, $P(N_{H}(t)=0)$, is given by:

[TABLE]

the probability that an individual is neither experiencing a relapse nor a primary infection at time $t$, P\big{(}N_{A}(t)+N_{P}(t)=0\big{)} (i.e., no blood-stage infection), is given by:

[TABLE]

and the probability that an individual is neither experiencing an infection nor has any hypnozoites in their liver at time $t$, P\big{(}N_{H}(t)=N_{A}(t)=N_{P}(t)=0\big{)}, is given by:

[TABLE]

2.2.2 Probability of blood-stage infected individual having one infection and no hypnozoites (under $N$ rounds of MDA)

With $p_{1}(t)$ defined as the probability that an individual has one infection ($N_{A}(t)+N_{P}(t)=1$) and an empty hypnozoite reservoir ($N_{H}(t)=0$) conditional on an ongoing blood-stage infection (i.e. primary infections or relapse, $N_{A}(t)+N_{P}(t)>1$), we have:

[TABLE]

The expression for $P(N_{H}(t)=0)$ and $P(N_{H}(t)+N_{P}(t)=0)$ follows from Equations (16) and (17). The expression for P\big{(}N_{A}(t)+N_{P}(t)=1|N_{H}(t)=0\big{)} can be obtained from Equation (27) (see appendix B for details) which is

[TABLE]

Finally, from Equation (19),

[TABLE]

where

[TABLE]

2.2.3 Probability of blood-stage infected individual having one infection and non-zero hypnozoites (under $N$ rounds of MDA)

The probability that a blood-stage infected individual experiencing only one infection ($N_{A}(t)+N_{P}(t)=1$) and has hypnozoites ($N_{H}(t)>0$) at time $t$, $p_{2}(t)$, is

[TABLE]

The expression $P(N_{A}(t)+N_{P}(t)=0)=P(N_{A}(t)=N_{P}(t)=0)$ is given by Equation (17). The expression for $P(N_{A}(t)+N_{P}(t)=1)$ follows from Equation (81) in Mehra et al. [37] and is given by

[TABLE]

where,

[TABLE]

2.2.4 Probability liver-stage infected individual has 1 hypnozoite in liver (under $N$ rounds of MDA)

The probability that a liver-stage infected individual has 1 hypnozoite in the liver at time $t$ (that is, the conditional probability for $N_{H}(t)$ given an individual does not have an ongoing blood-stage infection at time $t$) under $N$ MDA rounds is:

[TABLE]

where

[TABLE]

The expression for $P(N_{H}(t)=1|N_{A}(t)=N_{p}(t)=0)$ follows from Equation (78) in [37] and $P(N_{H}(t)=0|N_{A}(t)=N_{P}(t)=0)$ is obtained by dividing Equation (18) by Equation (17).

2.2.5 Average number hypnozoites within liver-stage infected individuals (under $N$ rounds of MDA)

The average number of hypnozoites within liver-stage infected individuals, $k_{T}(t)$, is defined by:

[TABLE]

where $\mathbb{E}\left[N_{H}(t)|N_{A}(t)=N_{P}(t)=0\right]$ is the expected size of the hypnozoite reservoir in an uninfected (no blood-stage infection) individual under $N$ rounds of MDA and is given by:

[TABLE]

The time-dependent parameters $p_{1}(t),\ p_{2}(t),\ k_{1}(t)$, and $k_{T}(t)$ that characterise the hypnozoite dynamics at the population level, account for all the infective bites received throughout time and change instantaneously with MDA because of the assumption of the instantaneous effect of the drug.

As these parameters involve numerical integration, we implement our own integro differential equation (IDE) solver using a $4^{th}$-order Runge-Kutta method, as described by Algorithm 1 in Anwar et al. [5]. Considering treatment at times $s_{1},\ s_{2},\ \ldots,\ s_{N}$, the parameters $p_{1}(t),\ p_{2}(t),\ k_{1}(t)$, and $k_{T}(t)$ are first obtained from the within-host model at each time $t$ to then obtain the solution of the population-level model at time $t$.

2.3 Optimisation model for the MDA intervals

In this section, we construct a mathematical optimisation model to obtain the optimal timing for each MDA round. Suppose $s_{1},s_{2},\ldots,s_{N}$ are the $N$ MDA administration times. We want to optimise the MDA intervention times so that the outcome of the MDA implementation is optimised. We construct the optimisation problem as:

[TABLE]

where $Z$ is the objective function to be minimised. Based on good public health, we investigate two objective functions:

• •

Z_{1}=\min_{t}\Big{(}I(t)+k_{T}(t)L(t)\Big{)},

• •

Z_{2}=\min_{t}\Big{(}\big{(}I(t)+L(t))W_{h}+(E_{m}(t)+I_{m}(t)\big{)}W_{m} ),

where $W_{h}$, $W_{m}$ are weighting factors for the human and mosquito population, respectively and $t\in[s_{1}\ t_{max}]$. That is, $Z_{1}$ is the minimum of the sum of the blood-stage infected proportion and the average hypnozoite burden in liver-stage infected individuals for $t\in[s_{1}\ t_{max}]$ and $Z_{2}$ is the minimum of the weighted sum of the proportion of infected humans (both blood-stage and liver-stage) and infected mosquitoes (exposed and infectious) for $t\in[s_{1}\ t_{max}]$. Since the P. vivax transmission is mainly dominated by hypnozoite dynamics, it is worth exploring the optimum effect of the drugs on disease prevalence and hypnozoite burden with the objective function, $Z_{1}$. As mosquito populations are an integral part in P. vivax transmission, we explore the potential effect of infected (exposed and infectious) mosquitoes along with infected humans with the objective function $Z_{2}$. By setting $W_{m}=0$, we can also investigate the optimal effect on only the human infected proportions (see Figure 7, for example).

2.3.1 Without seasonality

When seasonality is not considered, the time of the first MDA, $s_{1}$, can be considered arbitrary (as long as the dynamics have reached an equilibrium). In this case, we can fix $s_{1}=0$ (without loss of generality) and then the remaining MDA implementation times are optimised. Here, we minimise the objective function $Z$ over the time period of $[s_{1}\ t_{max}]$. We reconstruct the optimisation problem in terms of time intervals between MDA rounds. Let $x_{1},\ x_{2},\ \ldots,\ x_{N-1}$ be the time intervals between the first and second rounds of MDA, second and third rounds of MDA, and so on, respectively. That is, $x_{1}=s_{2}-s_{1},\ x_{2}=s_{3}-s_{2},\ \ldots,x_{N-1}=s_{N}-s_{N-1}$. Then the optimisation problem becomes:

[TABLE]

2.3.2 With seasonality

When considering seasonality in the mosquito population, the time of the first MDA round is no longer arbitrary as the dynamics are periodic oscillations around the mean annual prevalence. As the periodic function that governs the mosquito birth rate, $b_{m}(t)$, has an oscillation period of one year (assumed), the dynamics of human and mosquito populations have a peak within each year. Our optimisation problem without seasonality (Equation (25)) is constructed in terms of MDA intervals $x_{1},\ x_{2},\ \ldots,x_{N-1}$; for the case when we consider seasonality, we set a range of two years (starting from the time when prevalence is at a peak) for the optimisation algorithm to find the first MDA time, $s_{1}$. Here we define $x_{0}=s_{1}-\theta$ where $\theta$ is the peak prevalence time. That is, $x_{0}$ represents the interval between the prevalence peak time and the initial MDA time. The remaining times are obtained similarly without seasonality. Hence, the optimisation problem with seasonality in the mosquito population is:

[TABLE]

3 Results

In this section, we present some numerical results. First, we consider the effect of MDA rounds if there were no seasonality. We explore the effect of one MDA round on disease prevalence (as a function of human to mosquito ratio, $m$), liver-stage infected proportions, and the hypnozoite reservoir in Section 3.1. The effect of drug efficacy (varying $p_{rad}$) with one MDA round is presented in Section 3.2. We then present numerical results on the effect of multiple MDA rounds on disease prevalence (Section 3.3) by varying mosquito ratio where we present the rebound (e.g., minimum) disease prevalence obtained after 5 and 15 years for varying MDA rounds (up to $N=6$ rounds) with varying pre-MDA prevalence ($20\%-60\%$). Finally, we explore the effect of optimal MDA intervals on different disease prevalence by varying mosquito ratios for the two different objective functions constructed in the previous section, both with and without seasonality (Section 3.4).

3.1 The effect of a single round of MDA (with $p_{blood}=0.9,\ p_{rad}=0.9$)

To quantify the effect of radical cure MDA, we first assume that one round of MDA is applied when the system is at a steady state (see Appendix C for detail on the steady-state derivation). Treatment coverage plays a significant role in the effect of an MDA program [33]. To study the model behaviour, we assume that $100\%$ of the population is covered by the MDA scheme and that there is $90\%$ drug efficacy. Figure 3 shows the results from our multiscale model under one round of MDA. The drugs were assumed to have an instantaneous effect (with $p_{blood}=0.9$, $p_{rad}=0.9$); the hypnozoite reservoir size just before the MDA (Figure 3B) becomes smaller in size (Figure 3C; mode is [math]) as a result of the radical cure. That is, just immediately following MDA, most individuals will have no hypnozoites within their liver (with probability $\approx 0.7$). Disease prevalence drops significantly at the time of radical cure (Figure 3A), as we assume that the drug clears any ongoing blood-stage infections with $90\%$ efficacy ($p_{blood}=0.9$). For liver-stage infected individuals, as the drugs are assumed to kill each hypnozoite with probability $p_{rad}=0.9$, the overall effect of the drug depends on the size of the hypnozoite reservoir. If the size of the hypnozoite reservoir is substantial before the treatment, the overall effect would be insignificant, and vice versa. As individuals are still exposed to infectious mosquito bites, and each infective bite contributes to an average of $\nu$ number of hypnozoites that activate at a constant rate $\alpha$, both blood-stage and liver-stage proportions reach the same equilibrium state (Figure 3D) as before MDA (Figure 3B) eventually.

3.2 The effect of a single round of MDA, varying $p_{blood}$

The effect of the drug on disease transmission and the hypnozoite reservoir also changes with the efficacy of the drug (Figure 4). Figure 4A illustrates the effect of varying efficacies of the hypnozoicidal drug (i.e., $p_{rad}$) on liver-stage infected proportions. Figure 4B illustrates the hypnozoite distribution just after the application of MDA when $p_{blood}=p_{rad}=0.9$ and Figure 4C illustrates the hypnozoite distribution when $p_{blood}=0.9,\ p_{rad}=0.95$. If the hypnozoicidal drug were 100% effective (that is, $p_{rad}=1$) then all liver-stage infected individuals would recover (green line in Figure 4A). In the case of $p_{rad}=1$, the hypnozoite reservoir within the human population would be completely cleared (Figure 4D). In other words, immediately following drug administration, no individuals would be liver-stage infected. However, the disease will eventually reach the same equilibrium state as if no treatment were administered. (Figure 4A).

3.3 The effect of multiple MDA rounds

We also examined the impact of multiple MDA rounds on transmission and hypnozoite dynamics in the absence of seasonality (Figure 5). Figure 5A depicts the long-term behaviour of the transmission dynamics under four MDA rounds where the transient behaviour over the time of the MDA rounds (150 days) is depicted in Figure 5B. Here, we assumed a fixed interval (30 days) between MDA rounds, although intervals between MDA rounds among studies vary widely, from weeks to several months [39]. The effect of four successive MDA rounds is clearly visible in Figure 5A. Disease prevalence was driven down to approximately zero after the fourth round. However, as we model the system as a deterministic process and the effect of the drug is temporary, over time the disease reaches the same endemic steady state as before treatment. The overall effect of radical cure MDA treatment also depends on the disease prevalence; the lower the prevalence, the more effective the MDA in reducing the disease prevalence and hypnozoite-positive proportions. Figure 5C and D illustrate the sensitivity analysis of up to six MDA rounds at different assumed prevalences ($20\%\mbox{--}60\%$) obtained by varying mosquito ratio, $m$, showing the rebound prevalence 5 years and 15 years after the first MDA round was applied. The interval between each MDA round was again fixed at 30 days. If the prevalence before MDA is high, the dynamics reach the equilibrium state faster than when the prevalence is low before MDA.

3.4 Optimal MDA programs

3.4.1 Without seasonality

To obtain the optimal interval between MDA rounds, we use the optimisation problem defined in Equation (25). We used the MATLAB optimisation tool ‘Multistart’ (with 80 different initial starting points) with fmincon (SQP algorithm) to generate global optimal solutions. The results of two optimally timed MDA rounds are illustrated in Figure 6 for a steady state disease prevalence (see Appendix C for the derivation of the steady state disease prevalence) of 20% with the objective function $Z_{1}$. With our choice of parameter values (see Table 1), the optimisation problem gives an optimal interval of $x_{1}=s_{2}-s_{1}=34.7$ days, as illustrated in Figure 6A. Figure 6C depicts the sum of blood-stage infected population proportion and the hypnozoite burden on liver-stage infected population over time, $I(t)+k_{T}(t)L(t)$, before and after the MDA rounds using the optimal interval of $x_{1}=34.7$ days. The effect of the optimally timed MDA rounds on disease prevalence (20%) is depicted in Figure 6D. The dashed vertical lines in Figure 6C and D indicate the optimal time for the MDA rounds ($x_{1}=34.7$). When no seasonality is considered, the time of the first MDA can be at any arbitrary time (after an equilibrium has been reached). The equilibrium disease prevalence (obtained by varying mosquito ratio, $m$) greatly affects the optimum intervals (Figure 6B). The left vertical axis in Figure 6B illustrates the mosquito ratio, $m$, and the values on the right vertical axis depict the prevalence corresponding to each green bar. For higher prevalence ($25\%\mbox{--}60\%$), the optimisation problem with the objective function $Z_{1}$ suggests an interval of around 480 days between the two MDA rounds.

Figure 7 shows the optimum interval for two (first row) and three (second row) MDA rounds for different equilibrium disease prevalences (right vertical axis) obtained through changing the mosquito ratio, $m$, with three different choices of the objective function. The first, second, and third columns represent the objective function $Z_{1}$, $Z_{2}$ with $W_{h}=1,\ W_{m}=0$, and $Z_{2}$ with $W_{h}=1,\ W_{m}=1$, respectively. In contrast with the objective function $Z_{2}$ with $W_{h}=1,\ W_{m}=0$ and $Z_{2}$ with $W_{h}=W_{m}=1$, the optimisation problem suggests a longer interval for the second MDA round for higher prevalences ($>35\%$) with $Z_{1}$ when only two rounds of MDA are used (Figure 7B). The interval between the two MDA rounds is very similar for different prevalences for the objective function $Z_{2}$ with $W_{h}=1\ W_{m}=0$ and $Z_{2}$ with $W_{h}=W_{m}=1$ (Figure 7B, C).

The optimal intervals for three MDA rounds depend on both $m$, hence prevalence, and the choice of the objective function (Figure 7D-F). If three MDA rounds are considered, the optimisation problem (Equation (25)) suggests a similar interval for all of the MDA rounds with all three choices of the objective function for low prevalence ($<50\%$). But for higher prevalences ($>55\%$), for $Z_{1}$ and $Z_{2}$ with $W_{h}=1,\ W_{m}=0$, the optimisation routine suggests an immediate implementation of the third round of MDA after a long delay in between. However, for $Z_{2}$ with $W_{h}=1,\ W_{m}=0$, the interval $x_{2}$ becomes shorter as $m$, hence prevalence gets higher (but remains the same).

3.4.2 With seasonality

The effect of two optimally timed MDA rounds (including the first round, which was not required to be considered when seasonality was not considered) is illustrated in Figure 8. The optimal time for the first MDA round is approximately the same for different annual mean disease prevalences (right vertical axis, obtained by varying initial mosquito ratio, $m_{0}$) and the objective functions (Figure 8D–F). The seasonal amplitude, $\eta$, is thought to play an important role in intervention strategies [49]; here we have assumed $\eta=0.1$. Figure 8A–C shows the impact of two MDA rounds on disease prevalence for all the objective functions when there is a 54.9% annual mean disease prevalence before MDA for demonstrative purposes. The vertical solid line indicates the time when the pre-MDA prevalence reaches a peak and the vertical dashed lines indicate the time of the MDA implementations. When the annual mean disease prevalence is 54.9%, the optimisation problem with our choice of parameters as per Table 1 along with the objective function $Z_{1}$, provides the interval $x_{0}=103.4$ days and $x_{1}=26.7$ days for two MDA rounds. The intervals with $Z_{2}\,(W_{h}=1,W_{m}=0)$ are $x_{0}=132.1$ days, $x_{1}=34.3$ days and with $Z_{2}\,(W_{h}=W_{m}=1)$ are $x_{0}=133.7$ days and $x_{1}=33.7$ days. The sensitivity analysis for optimal interval time with different annual mean prevalences (right vertical axis, corresponding to each bar) is illustrated in Figure 8D–F with $Z_{1}$, $Z_{2}\,(W_{h}=1,\ W_{m}=0)$ and $Z_{2}\,(W_{h}=W_{m}=1)$, respectively. The red rectangles in Figure 8D–F indicate the optimal intervals corresponding to Figure 8A–C. With respect to all objective functions, $Z_{1}$, $Z_{2}\,(W_{h}=1,\ W_{m}=0)$, and $Z_{2}\,(W_{h}=W_{m}=1)$, the optimal intervals are very similar when the mean annual prevalence is low ($<50\%$). However, the optimisation algorithm suggests an immediate implementation for the two MDA rounds for higher annual mean prevalence with the objective function $Z_{2}\,(W_{h}=1,\ W_{m}=0)$ (Figure 8E), while with $Z_{2}\,(W_{h}=W_{m}=1)$, the algorithm suggests a similar interval as for low prevalences (Figure 8F). With objective function $Z_{1}$, the interval between the two MDA rounds is also very similar for different annual mean prevalences (Figure 8D).

Figure 9 shows the optimal intervals when three MDA rounds are considered for each objective function. Figure 9A–C demonstrates the effect of three optimally timed MDA rounds on the objective function $Z_{1}$, $Z_{2}\,(W_{h}=1,\ W_{m}=0)$ and $Z_{2}\,(W_{h}=W_{m}=1)$, respectively over time where the vertical solid line indicates the time when the prevalence reaches a peak and the three subsequent vertical dashed lines indicate the optimal time for the three MDA rounds. The sensitivity analysis for optimal interval time with different annual mean prevalences (right vertical axis) is illustrated in Figure 9D–F with $Z_{1}$, $Z_{2}\,(W_{h}=1,\ W_{m}=0)$ and $Z_{2}\,(W_{h}=W_{m}=1)$, respectively where the violet rectangles in Figure 9D–F indicate the optimal intervals corresponding to Figure 9A–C. With respect to all objective functions, $Z_{1}$, $Z_{2}\,(W_{h}=1,\ W_{m}=0)$, and $Z_{2}\,(W_{h}=W_{m}=2)$, the optimal timing for the second and third MDA round, that is the interval between the last two rounds is almost identical throughout all different prevalences.

The choice of model parameters can significantly influence the optimal MDA intervals. In order to obtain an equilibrium disease prevalence (without seasonality, see Appendix C) for Figures 5 and 6, we only varied the human-to-mosquito ratio ($m$) and kept all other parameter values as per Table 1. We note that there are (possibly) many other combinations of model parameters that could generate the same equilibrium prevalence (see Figure 10). Hence, we performed a sensitivity analysis for the parameters $m,\ \alpha,\ \mu,\ \gamma$, and $\nu$ on the optimal interval (without seasonality) for two MDA rounds with the objective function $Z_{1}$ which is the minimum of the sum of the blood-stage infected proportion and the average hypnozoite burden in liver-stage infected individuals at time $t$. Figure 10A–E depict the effect of varying $m,\ \alpha,\ \mu,\ \gamma$, and $\nu$ on the optimal intervals, respectively. The color of the bars in Figure 10A–E depicts the equilibrium prevalence corresponding to the parameter value. The red arrows in Figure 10A–E depict the baseline parameters in Table 1 that generate a prevalence of $40\%$ as shown in Figure 6B. As illustrated in Figure 10A, the abundance of mosquitoes can drastically influence the disease equilibrium as the force of reinfection (that is, the probability of reinfection per unit time) increases with $m$ $(\lambda=m_{0}abI_{m})$. We see that the optimal intervals can be different for the same equilibrium prevalence generated with a different combination of the parameters $m,\ \alpha,\ \mu,\ \gamma$, and $\nu$. That is, the optimal interval without seasonality depends on the input model parameters. To investigate this further, we vary the values of $\alpha$ and set a value of $m$ such that the steady-state disease prevalence is $30\%$ (Figure 11). Figure 11A illustrates the distribution of the optimal interval for two MDA rounds for different values of $\alpha$ and $m$. The optimal interval varies from around 47 days to 446 days for the different combinations of $\alpha$ and $m$ (objective function $z_{1}$). Figure 11B depicts the distribution of optimal intervals for the same set of parameters but with the objective function $Z_{2}$ with $W_{h}=1,W_{m}=0$. In this case, the optimal interval varies from around 37 days to 172 days. The results illustrate that prevalence alone is not sufficient to determine an optimal MDA interval when there is no seasonality in the mosquito population.

The jump in optimal interval seen in Figure 8A, 10, and 11 as we vary model parameters is related to the choice of the objective function. Regardless of the choice of model parameters (and hence disease prevalence), the effect of the drug on the blood-stage infected population ($I$) is to cause an instantaneous reduction at the time of the MDA corresponding to the effectiveness of the drug ($p_{blood}$). Hence, to minimise the blood-stage infected population alone, the optimisation will always suggest the immediate implementation of the second MDA round. Similarly, the effect of the drug on the hypnozoite reservoir (which has an average size, $k_{T}$) is always to reduce its size regardless of model parameters (and disease prevalence). However, since the effect of the drug on $k_{T}$ will be more for larger hypnozoite reservoir sizes, to minimize the hypnozoite reservoir alone, the optimisation would suggest a longer interval (regardless of disease prevalence) so that the reservoir has time to build up before the next MDA round. In contrast, the effect of the drug on the liver-stage infected population ($L$) does vary with model parameters (and disease prevalence). For low disease prevalence, the average hypnozoite reservoir size will be smaller (see Equation (34)) in which case $L$ will decrease at the time of the first MDA application. For higher disease prevalence, the average hypnozoite reservoir size will be larger and it is possible that $L$ will increase at the time of the first MDA application since those in $I$ have their blood-stage infection cleared but not all of their hypnozoites due to the large average hypnozoite reservoir size. The objective function, $Z_{1}$, considers minimising both blood-stage infections ($I$) and hypnozoite burden within liver-stage infected fractions ($k_{T}L$). This increase in the liver-stage infected fractions when the first MDA is applied under higher disease prevalence means that a longer interval for the second MDA will be optimal to reduce the burden (Figure 8A, 10) while when prevalence is low the liver-stage infected fraction will decrease with the first MDA and a second MDA round within a short interval will be optimal. Furthermore, if we consider seasonality in the mosquito populations, the results are quite different. Figure 12 illustrates the distribution of the optimal intervals ($x_{0}$ and $x_{1}$) when the annual mean prevalence is approximately $30\%$ for objective function $z_{1}$ (Figure 12A) and for objective function $Z_{2}$ with $W_{h}=1,W_{m}=0$ (Figure 12B). The distribution of the optimal interval is consistent for both objective functions. The results indicate that when there are fluctuations in the abundance of mosquitoes in the environment, the optimal interval can be identified by measuring the prevalence.

4 Discussion

Targeting the hypnozoite reservoir is the most crucial action in any P. vivax elimination strategy, as hypnozoites dominate P. vivax transmission dynamics. In malaria elimination efforts around the world, interest in MDA using primaquine or tafenoquine has grown, as these are the only available drugs to treat liver-stage P. vivax infections [24]. In this paper, we have developed a multiscale model that captures hypnozoite dynamics and the effect of the hypnozoite reservoir on disease transmission under radical cure treatment as a method of MDA. This model extends our previous work [5] by integrating treatment into the model with multiple MDA rounds accounting for superinfection. We have extended Mehra et al.’s [37] within-host model by obtaining key parameters regarding hypnozoite dynamics under multiple MDA rounds and embedding these into a population-level transmission model that considers superinfection based on [35]. We have also included mosquito seasonality in our model to study the impact of MDA treatment when there is a seasonal effect on mosquitoes in the environment. According to our model and choice of parameters, MDA with radical cure can significantly reduce disease burden at the time the program is administered and maintain it at low levels when prevalence before the MDA intervention is low and if multiple MDA rounds are implemented (Figure 5). Our model results are sensitive to some parameters, especially for parameter regimes where superinfection is likely. However, we found that the optimal MDA intervals for a specific objective depend on the parameter values (without seasonality), especially the ones that have more influence on the transmission dynamics (mosquitoes per human, hypnozoite activation rate, hypnozoite death rate, natural recovery rate, and average hypnozoite per mosquito bite). That is, even where different combinations of the model parameters correspond to the same equilibrium prevalence, the optimal intervals are not necessarily the same (without seasonality, Figure 11) and 10). However, when there is seasonal variation in the mosquito population in the environment, the optimal intervals are very similar (Figure 12) for different combinations of the model parameters that correspond to the same annual mean prevalence. Hence, prevalence alone should not be considered a reliable measure when determining optimal intervals between rounds of MDA, especially in regions where seasonal variation in the mosquito population is negligible.

Although the optimal interval, frequency, and population coverage with MDA are not clear in practice [34, 21, 24], here we assume 100% treatment coverage and that all drugs (both blood-stage and liver-stage) are 90% effective. This assumption about the effectiveness of the radical cure drug is realistic, as studies show radical cure efficacy varies between 57·7% and 95% depending on the combination of drugs [25, 38, 32]. According to our model, the optimal intervals between MDA rounds vary with the prevalence before MDA, the number of MDA rounds under consideration, and the choice of the objective function (Figures 7, 8, 9). However, regardless of the objective and number of MDA rounds, the overall effect of the drug is only temporary under our model assumptions. This temporary effect is due to the assumption of the instantaneous effect of the drugs. This assumption is appropriate given that available drugs have half-lives varying from 3.7 hours to 28 days [28, 48] which is short compared to the time frame of interest (years). Hence, in the long term, the dynamical system does not observe any drug effect and the system returns to its pre-MDA state, which is the expected outcome from a deterministic framework such as ours. A deterministic framework is useful to understand the disease dynamics for a large population size however for a small population size, it will be important to use a stochastic model to study disease-extinction scenarios [3]. Currently, prophylaxis is not taken into account in our model. Accounting for prophylaxis might vary model outcomes, as a longer duration of prophylaxis leads to greater measured efficacy, especially in higher transmission settings [25]. Furthermore, given the mosquito population has a shorter lifespan, for a longer duration of prophylaxis a reasonable proportion of infectious mosquitoes may die out and disrupt the chains of transmission. The assumption of blood-stage infection clearance in the presence of superinfection is slightly different in the population model in comparison to the within host model. The within host model assumes that each blood-stage infection is cleared independently for analytical tractability [37]. However, since we are not aware of any study that suggests that the blood-stage drugs act differently on each blood-stage infection, we assumed that the clearance of all blood-stage infections (regardless of how many there are) depends only on the efficacy of the drug, $p_{blood}$.

Although being an effective intervention strategy, MDA has some disadvantages, especially in terms of drug resistance [59, 15]. Because of the extensive use of antimalarial drugs, the parasite has developed resistance to some drugs, particularly chloroquine. However, chloroquine is still effective in most parts of the world for P. vivax [41]. Another challenge with MDA is the use of the anti-hypnozoicidal drugs primaquine and tafenoquine, as these can cause blood hemolysis in individuals with G6PD deficiency and problems in pregnant women [23, 54]. We do not consider G6PD deficiency in our model, but it could easily be extended to do so. We also do not consider drug resistance, immunity, or heterogeneity in bite exposure.

Since our model is deterministic, disease fade-out is not possible, but a disease in a real-life setting may undergo stochastically driven fade-out when the disease prevalence is sufficiently low [30, 20]. The primary purpose of this work is to optimise the implementation of the timing of the rounds of MDAs. However, disease elimination could be investigated with our multiscale model by approximating the elimination probability as a Binomial random variable. As P. vivax parasites are transmitted through infectious mosquito bites, contributing to hypnozoites in the liver, it is as important to reduce mosquito-bite exposure or the abundance of mosquitoes as it is to clear hypnozoites from the liver [31, 45, 39]. Insecticide-treated nets, indoor residual spraying, and long-lasting insecticide–treated nets are some of the standard vector-control interventions for controlling malaria transmission and are necessary additional interventions alongside MDA as per the WHO guidelines [59]. Including vector-control interventions with MDA and stochasticity in the model to obtain the probability of disease eradication is an avenue for potential future work.

To our knowledge, ours is the first multiscale model to provide a framework for studying the effect of multiple MDA rounds in both the within-host and population scale for P. vivax transmission. The results from the model demonstrate the effect of several MDA rounds delivered at optimal intervals on both the transmission setting and hypnozoite dynamics. According to our model, P. vivax transmission can only be interrupted for a certain period (the duration of which depends on the prevalence before MDA) when using MDA. That is, MDA alone is not sufficient to progress us towards sustained P. vivax elimination under our model. While our model has not been parameterised for any particular geographical setting, it has the potential to aid policymakers in MDA control strategy decision-making.

5 Funding

M.N. Anwar is supported by a Melbourne Research Scholarship. J.M. McCaw’s research is supported by the Australian Research Council (DP170103076, DP210101920) and the NHMRC Australian Centre of Research Excellence in Malaria elimination (ACREME). J.A. Flegg’s research is supported by the Australian Research Council (DP200100747, FT210100034).

6 Data availability

Data sharing is not applicable to this article as no datasets were generated or analysed during the current study.

Appendix A Model derivation with mosquito seasonality

Let $X$, $Y$, and $Z$ represent the number of susceptible, blood-stage and liver-stage infected individuals and $U$, $V$, and $W$ represent the number of susceptible, exposed and infectious mosquitoes. Let $N_{h}=X+Y+Z$ be the total human population and $N_{m}(t)=U+V+W$ be the total mosquito population at time $t$, respectively. With mosquito seasonality, the model equations for the number of individuals in each compartment are:

[TABLE]

where

[TABLE]

All model parameters are defined in Table 1. Here, $\mathrm{d}(N_{h})/\mathrm{d}t=\mathrm{d}(X+Y+Z)/\mathrm{d}t=0$ so that the human population is constant in size over time. But for the mosquito population:

[TABLE]

where $N_{m}(0)$ is the initial mosquito population size.

We now convert the above transmission model into a model on the proportion scale for consistency with the model given in Equations (2)-(6). Let $S=X/N_{h},\ I=Y/N_{h},\ L=Z/N_{h},\ S_{m}=U/N_{m}(t),\ E_{m}=V/N_{m}(t),\ I_{m}=W/N_{m}(t)$. Therefore

[TABLE]

The equations for the human population on the proportion scale become:

[TABLE]

And the equations for the mosquitoes on the proportion scale become:

[TABLE]

The model dynamics with mosquito seasonality in comparison with no seasonality are depicted in Figure 13.

Appendix B Multiple infections given blood-stage infected

An individual might experience multiple blood-stage infections at the same time either due to bites from infectious mosquitoes or relapses from hypnozoite activation. We define the multiplicity of infection (MOI) as the number of distinct parasites co-circulating within a blood-stage infected individual. Thus, the multiplicity of infection (MOI) is given by the total number of bloods-stage infections (infections from mosquito bites and relapses) at time $t$: $M_{I}(t)=N_{A}(t)+N_{P}(t)$. Now, multiplicity of infection given empty hypnozoite reservoir: $M_{I}(t)|N_{H}(t)=0$ can be obtained from the PGF given by Equation (13) that holds for before treatment and by Equation (14) which holds following treatment at times $s_{1},\ s_{2},\ \ldots,s_{N}$ as

[TABLE]

where

[TABLE]

Now, the probability mass function for $M_{I}(t)|N_{H}(t)=0$ is

[TABLE]

where

[TABLE]

Appendix C Steady state analysis (without seasonality)

To obtain the steady state of the system without treatment and seasonality, first assume that the time-dependent parameters $p_{1}(t),\ p_{2}(t),\ k_{1}(t)$, and $k_{T}(t)$ are in steady state. We define

[TABLE]

Therefore at steady state, we have:

[TABLE]

From Equations (30), (31) and (33), at steady state we have:

[TABLE]

Substituting the value of $S_{m}$ and $E_{m}$ in Equation (32), we get

[TABLE]

We have a constant human population, therefore

[TABLE]

Now, substituting the value of $S$ and $L$ in Equation (28) gives:

[TABLE]

To obtain a steady state prevalence, $I^{*}$, we need $\lambda=\lambda_{SS}$ to be at a steady state. Each of the parameters $\bar{p_{1}},\ \bar{p_{2}},\ \bar{k_{1}}$, and $\bar{k_{T}}$ can be expressed as a function of $\lambda_{SS}$ from the Equations (2.2.2), (2.2.3), (2.2.4), and (2.2.5) respectively. That is,

[TABLE]

which makes Equation (34) a nonlinear function of $\lambda_{SS}$. That is

[TABLE]

Therefore, given a fixed $I^{*}$, we can solve Equation (35) to obtain $\lambda_{SS}$ and hence can find the other human and mosquito proportions at steady state. Note that the steady state disease prevalence, $I^{*}$, is obtained as a function of the human to mosquito ratio, $m$. That is, by varying $m$, we can vary $\lambda_{SS}$ and hence $I^{*}$. The steady state solution following this analysis is illustrated in Figure 14 and captures the long-term behaviour of the model based on numerical simulation.