1 Introduction

Plasmodium vivax is a parasite that causes malaria, responsible for 4.5 million cases in 2020 (World Health Organization 2021). 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’ (White 2014; Ricardo Águaset al. 2012). Because of the relapse characteristics, P. vivax has become the most globally widespread parasite and is responsible for significant morbidity and mortality (Antinori 2012; Battle 2019). 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 (Price 2020).

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 (Mehra et al. 2022b; Price 2020). An estimated 79–96% of the total vivax cases are due to hypnozoite activation (Robinson 2015; Huber 2021; Adeshina and Adekunle 2015; Commons 2020). Thus, targeting the hypnozoite reservoir with treatment is an important element of any program for P. vivax elimination (Campo 2015). Mass drug administration (MDA) is an effective intervention for controlling many diseases and was advocated by the World Health Organization (WHO) in the 1950 s to control malaria transmission (Hsiang 2013). MDA involves treating the entire population, or a well-defined sub-population, in a geographic location regardless of their infection status (Newby 2015; Hsiang 2013). 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 (Timothy et al. 2010; Taylor 2019; Poespoprodjo 2022). 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 (World Health Organization 2021; Howes Rosalind 2012; Watson 2018).

The effect of radical cure treatment on P. vivax transmission has been explored in a number of mathematical models (Ishikawa 2003; Ricardo Águaset al. 2012; Chamchod and Beier 2013; Roy Manojit 2013; White 2014, 2016, 2018). However, most of these models do not consider the variation in hypnozoite number per mosquito bite (Ishikawa 2003; Ricardo Águaset al. 2012; Chamchod and Beier 2013; Roy Manojit 2013; White 2016). Mehra et al. (2022b) 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 (Anwar 2022) 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 (Phommasone 2020). 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 (Kaehler 2019; Jambulingam 2016). 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 (Herdicho et al. 2021; Buonomo and Marca 2018; Kabirul and Nobuko 2014; Galardo 2009). 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 (Ishikawa 2003; Chamchod and Beier 2013; White 2014; Silal 2019; Mehra et al. 2022a), few have also captured the rich dynamics of hypnozoites (Mehra et al. 2022a) and none have considered the impact of multiple drug administration explicitly on hypnozoite dynamics. 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 (Smith 2012; Dietz et al. 1974). However, only a few P. vivax mathematical transmission models account for superinfection (White 2014, 2018; Silal 2019; Mehra 2022; Mehra et al. 2022a, b).

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 Sect. 2, we extend our existing multiscale model (Anwar 2022) 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 (Mehra et al. 2022b) under multiple rounds of MDA and obtain the recovery rate under superinfection. In Sect. 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 Sect. 4.

2 Methods

A multiscale mathematical model that accounts for hypnozoite variation within individuals without treatment has already been developed (Anwar 2022). 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 (Anwar 2022), 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 (2022). 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 (blood-stage infected individuals may or may not have hypnozoites), either become susceptible (S) or liver-stage infected (L). Following the work of Mehra (2022), 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 )\).

Fig. 1
figure 1

Schematic illustration of the multiscale model with treatment. S, I and L represent the fractions of the human population that are susceptible with no hypnozoites, blood-stage infected, and liver-stage infected, respectively. Note that, blood-stage infected individuals may or may not have hypnozoites, but liver-stage infected individuals have at least 1 hypnozoite. The left (top and bottom) part of the schematic demonstrates the transmission dynamics between the human and mosquito populations while the right part of the schematic demonstrates how the within-host model has been embedded within the population scale model. The within-host model takes into account the history of infective bites and calculates the probability of blood-stage infected individuals having 0 hypnozoites and one blood-stage infection (\(p_1(t)\)), blood-stage infected individuals having more than 0 hypnozoites and one blood-stage infection (\(p_2(t)\)), liver-stage infected individuals having 1 hypnozoite (\(k_1(t)\)), the expected size of the hypnozoite reservoir (\(k_T(t)\)), and the probability of blood-stage infected individuals having 0 hypnozoites (p(t)) at any given time t as a function of the force of infection, \(\lambda (t)\). The red area on the right part of the schematic indicates the force of infection from time \(t=0\) to t. The functions \(D_b(t)\) and \(D_l(t)\) capture the effect of treatment when implemented. Other parameters are defined in Table 1 in “Appendix”

Full size image

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 (Sect. 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 Sect. 2.2. The model schematic is depicted in Fig. 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 (Mehra et al. 2022b). That is, on administration, all ongoing blood-stage infections are instantaneously cleared with probability \(p_{\textrm{blood}}\), and each hypnozoite in the liver dies instantaneously with probability \(p_{\textrm{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_{\textrm{blood}}\). We define p(t) as the probability of blood-stage infected individuals having no hypnozoites in their liver (Anwar 2022). 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 h (Birgersson 2016) and primaquine, which kills hypnozoites have a relatively short half-life of approximately 3.7\(-\)9.6 h (Jittamala 2015). 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 (Schlagenhauf 2019). Since the number of mosquitoes in the environment influences P. vivax dynamics dramatically (Herdicho et al. 2021; Buonomo and Marca 2018), it is important to account for seasonal environmental effects on the mosquito population (see, for example, Kabirul and Nobuko 2014; Galardo 2009). 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:

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

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):

$$\begin{aligned}&\frac{\textrm{d}S}{\textrm{d}t}=-\lambda S+\mu k_1(t)L+p_1(t)\gamma I+D_l(l)L+D_b(t)p(t)I, \end{aligned}$$
(1)
$$\begin{aligned}&\frac{\textrm{d}I}{\textrm{d}t}=\lambda (S+I)+\alpha k_T(t)L-\gamma \big (p_1(t)+p_2(t)\big ) I-D_b(t)I, \end{aligned}$$
(2)
$$\begin{aligned}&\frac{\textrm{d}L}{\textrm{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, \end{aligned}$$
(3)
$$\begin{aligned}&\frac{\textrm{d}S_m}{\textrm{d}t}=b_m(t)-acIS_m-b_m(t)S_m, \end{aligned}$$
(4)
$$\begin{aligned}&\frac{\textrm{d}E_m}{\textrm{d}t}=acIS_m-\left( b_m(t)+n\right) E_m, \end{aligned}$$
(5)
$$\begin{aligned}&\frac{\textrm{d}I_m}{\textrm{d}t}=nE_m-b_m(t)I_m, \end{aligned}$$
(6)

where,

$$\begin{aligned} \lambda =m_0ab I_m \text {exp}\left\{ \frac{365g\eta }{2\pi } \sin \left( \frac{2\pi t}{365}+\phi \right) \right\} , \end{aligned}$$

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:

$$\begin{aligned} D_b(t)&=\ln {\big ((1-p_{\textrm{blood}})^{-1}\big )}\sum _{j=1}^N \delta _D(t-s_j),\\ D_l(t)&=\big \{k_1(t)\ln {\big ((1-p_{\textrm{rad}})^{-1}}\big )+k_2(t)\ln {\big ((1-p_{\textrm{rad}}^2)^{-1}}\big )+\ldots \big \}\sum _{j=1}^N \delta _D(t-s_j)\\&=\sum _{i=1}^\infty \ln {\Big ({1-p_{\textrm{rad}}^i}}\Big )^{-k_i(t)} \sum _{j=1}^N \delta _D(t-s_j), \end{aligned}$$

where \(\delta _D(\cdot )\) is the Dirac delta function. That is, any blood-stage parasite will be instantaneously cleared with probability \(p_{\textrm{blood}}\) every time the radical cure is administered (Mehra et al. 2022b) 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 (Eq. 20) and \(p_2(t)\) which is the probability that an individual experiencing only one infection has hypnozoites in the liver given blood-stage infection (Eq. 21), 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_{\textrm{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 in “Appendix”.

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. (2022b). 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)

$$p^r_H(t,s_1,\ldots ,s_N),\ p^r_A(t,s_1,\ldots ,s_N),\ p^r_C(t,s_1,\ldots ,s_N),\ p^r_D(t,s_1,\ldots ,s_N)$$

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

$$\begin{aligned} \frac{\textrm{d}p^r_H}{\textrm{d}t}=&-(\alpha +\mu )p^r_H-\ln {\big ((1-p_{\text {rad}})\big )^{-1}}\sum _{j=1}^N \delta _D(t-s_j)p^r_H, \end{aligned}$$
(7)
$$\begin{aligned} \frac{\textrm{d}p^r_A}{\textrm{d}t}=&-\gamma p^r_A+\alpha p^r_H-\ln {\big ((1-p_{\text {blood}})\big )^{-1}}\sum _{j=1}^N \delta _D(t-s_j)p^r_A,\end{aligned}$$
(8)
$$\begin{aligned} \frac{\textrm{d}p^r_C}{\textrm{d}t}=&\gamma p^r_A+\ln {\big ((1-p_{\text {blood}})\big )^{-1}}\sum _{j=1}^N \delta _D(t-s_j)p^r_A,\end{aligned}$$
(9)
$$\begin{aligned} \frac{\textrm{d}p^r_D}{\textrm{d}t}=&-\mu p^r_H+\ln {\big ((1-p_{\text {rad}})\big )^{-1}}\sum _{j=1}^N \delta _D(t-s_j)p^r_H, \end{aligned}$$
(10)

where the parameters \(\alpha ,\ \gamma ,\) and \(\mu \) are as per Table 1. Since our population model in Eqs. (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)\)) (Anwar 2022); hence we solve Eqs. (7)–(8) for N MDA rounds to give:

$$\begin{aligned} p^r_H(t,s_1,s_2,s_3,\ldots , s_N)&=(1-p_{\text {rad}})^N p_H(t), \end{aligned}$$
(11)
$$\begin{aligned} p^r_A(t,s_1,s_2,s_3,\ldots , s_N)&=(1-p_{\text {blood}})e^{-\gamma (t-s_N)}p_A^r(s_N,s_1,s_2,\ldots ,s_{N-1})\nonumber \\&+(1-p_{\text {rad}})^N\big (p_A(t)-e^{-\gamma (t-s_N)}p_A(s_N)\big ), \end{aligned}$$
(12)

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:

$$\begin{aligned} p_H(t)&=e^{-(\alpha +\mu )t},\\ p_A(t)&=\frac{\alpha }{(\alpha +\mu )-\gamma }\left( e^{-\gamma t}-e^{-(\alpha +\mu )t}\right) . \end{aligned}$$

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_{\textrm{blood}}=p_{\textrm{rad}}=0.5\)) is illustrated in Fig. 2A, 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.

Fig. 2
figure 2

Effect of radical cure (three rounds of MDA) on a single hypnozoite. A Probability of hypnozoite establishment as per Eq. (11) and B probability of hypnozoite activation as per Eq. (12). For each subplot, blue represents the probability without considering treatment and orange represents treatment assuming 50% efficacy \((p_{\textrm{blood}}=p_{\textrm{rad}}=0.5)\) of the drugs. The vertical lines indicate the times of drug administration. Other parameters are as in Table 1

Full size image

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':=\left\{ H,A,C,D,P,PC\right\} \), the probability generating function (PGF) for

$$\begin{aligned} \textbf{N} (t)=(N_H(t),N_A(t),N_C(t),N_D(t),N_P(t),N_{PC}) \end{aligned}$$

with \(\textbf{N} (0)=\textbf{0}\) can be written following from Equation (30) in Mehra et al. (2022b) (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):

$$\begin{aligned}&G(t,z_H,z_A,z_C,z_D,z_P,z_{PC}):=\mathbb {E}\left[ \displaystyle \prod _{f\in F'} z_f^{N_f(t)}\right] \nonumber \\&\quad =\text {exp}\left\{ -q(t)+\int _0^t \frac{\lambda (\tau )\left( z_Pe^{-\gamma (t-\tau )}+(1-e^{-\gamma (t-\tau )})z_{PC}\right) }{1+\nu \left( 1-\sum _{f\in F} z_fp_f(t-\tau )\right) }\textrm{d}\tau \right\} , \end{aligned}$$
(13)

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

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

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. (2022b). 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:

$$\begin{aligned}&G^{s_1,s_2,\ldots s_N}(t,z_H,z_A,z_C,z_D,z_P,z_{PC}):=\mathbb {E}\left[ \displaystyle \prod _{f\in F'} z_f^{N^{s_1,s_2,\ldots s_N}_s(t)}\right] \nonumber \\&\quad =\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) }\textrm{d}\tau \nonumber \\&\qquad +\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^r_s(t-\tau ,s_1-\tau )\right) }\textrm{d}\tau \nonumber \\&\qquad +\int _{s_1}^{s_2} \lambda (\tau )\frac{(1-p_{\text {blood}})^2e^{-\gamma (t-\tau )}z_P+(1-(1-p_{\text {blood}})^2e^{-\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 )\right) }\textrm{d}\tau \nonumber \\&\qquad +\int _{s_2}^{s_3} \lambda (\tau )\frac{(1-p_{\text {blood}})^3e^{-\gamma (t-\tau )}z_P+(1-(1-p_{\text {blood}})^3e^{-\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 )\right) }\textrm{d}\tau \nonumber \\&\qquad +\ldots +\int _{s_{N-1}}^{s_N} \lambda (\tau )\frac{(1-p_{\text {blood}})^Ne^{-\gamma (t-\tau )}z_P+(1-(1-p_{\text {blood}})^Ne^{-\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) }\textrm{d}\tau \Bigg \}. \end{aligned}$$
(14)

We now use the PGF in Eq. (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 (2022) we have:

$$\begin{aligned} p(t)&=P\big (N_H(t)=0|N_A(t)>0 \cup N_P(t)>0\big )\nonumber \\&=\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 )}. \end{aligned}$$
(15)

where the probability that an individual has an empty hypnozoite reservoir at time t, \(P(N_H(t)=0)\), is given by:

$$\begin{aligned}&P(N_H(t)=0)=G^{t,s_1,s_2,\ldots ,s_N}(t,z_H=0, z_A=1, z_C=1, z_D=1, z_P=1, z_{PC}=1)\nonumber \\&= {\left\{ \begin{array}{ll} \text {exp}\Bigl \{-q(t)+\int _0^t \frac{\lambda (\tau )}{1+\nu p_H(t-\tau )}\textrm{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 )}\textrm{d}\tau +\int _0^{s_1} \frac{\lambda (\tau )}{1+\nu p_H^r(t-\tau ,s_1-\tau )}\textrm{d}\tau \\ \qquad +\int _{s_1}^{s_2} \frac{\lambda (\tau )}{1+\nu p_H^r(t-\tau ,s_1-\tau ,s_2-\tau )}\textrm{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 )}\textrm{d}\tau \Bigr \}&{} \text {if}\ t\ge s_N, \end{array}\right. } \end{aligned}$$
(16)

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:

$$\begin{aligned}&P\big (N_A(t)+N_P(t)=0\big )=G^{t,s_1,s_2,\ldots ,s_N}(t,z_H=1, z_A=0, z_C=1, z_D=1, z_P=0, z_{PC}=1) \nonumber \\&= {\left\{ \begin{array}{ll} \text {exp}\Bigl \{-q(t)+\int _0^t \frac{\lambda (\tau )(1-e^{-\gamma (t-\tau )})}{1+\nu p_A(t-\tau )}\textrm{d}\tau \Bigr \} &{} \text {if}\ t< s_1\\ \text {exp}\Bigl \{-q(t)+\int _{s_N}^t \frac{\lambda (\tau )(1-e^{-\gamma (t-\tau )})}{1+\nu p_A(t-\tau )}\textrm{d}\tau \\ \qquad +\int _0^{s_1} \frac{\lambda (\tau )(1-(1-p_{\textrm{blood}})e^{-\gamma (t-\tau )})}{1+\nu p_A^r(t-\tau ,s_1-\tau )}\textrm{d}\tau \\ \qquad +\int _{s_1}^{s_2} \frac{\lambda (\tau )(1-(1-p_{\textrm{blood}})^2e^{-\gamma (t-\tau )})}{1+\nu p_A^r(t-\tau ,s_1-\tau ,s_2-\tau )}\textrm{d}\tau \\ \qquad +\ldots +\int _{s_{N-1}}^{s_N} \frac{\lambda (\tau )(1-(1-p_{\textrm{blood}})^Ne^{-\gamma (t-\tau )})}{1+\nu p_A^r(t-\tau ,s_1-\tau ,\ldots ,s_N-\tau )}\textrm{d}\tau \Bigr \} &{} \text {if}\ t\ge s_N. \end{array}\right. } \end{aligned}$$
(17)

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:

$$\begin{aligned} P\big (N_H(t)&=N_A(t)=N_P(t)=0\big )\nonumber \\&=G^{t,s_1,s_2,\ldots ,s_N}(t,z_H=0, z_A=0, z_C=1, z_D=1, z_P=0, z_{PC}=1)\nonumber \\&={\left\{ \begin{array}{ll} \text {exp}\Bigl \{-q(t)+\int _0^t \frac{\lambda (\tau )(1-e^{-\gamma (t-\tau )})}{1+\nu (p_H(t-\tau )+p_A(t-\tau )}\textrm{d}\tau \Bigr \} &{} \text {if}\ t< s_1\\ \text {exp}\Bigl \{-q(t)+\int _{s_N}^t \frac{\lambda (\tau )(1-e^{-\gamma (t-\tau )})}{1+\nu (p_H(t-\tau )+p_A(t-\tau )}\textrm{d}\tau \\ \qquad +\int _0^{s_1} \frac{\lambda (\tau )(1-(1-p_{\textrm{blood}})e^{-\gamma (t-\tau )})}{1+\nu (p_H^r(t-\tau ,s_1-\tau )+p_A^r(t-\tau ,s_1-\tau )}\textrm{d}\tau \\ \qquad +\int _{s_1}^{s_2} \frac{\lambda (\tau )(1-(1-p_{\textrm{blood}})^2e^{-\gamma (t-\tau )})}{1+\nu (p_H^r(t-\tau ,s_1-\tau ,s_2-\tau )+p_A^r(t-\tau ,s_1-\tau ,s_2-\tau )}\textrm{d}\tau \\ \qquad +\ldots +\int _{s_{N-1}}^{s_N} \frac{\lambda (\tau )(1-(1-p_{\textrm{blood}})^Ne^{-\gamma (t-\tau )})}{1+\nu (p_H^r(t-\tau ,s_1-\tau ,\ldots ,s_N-\tau )+p_A^r(t-\tau ,s_1-\tau ,\ldots ,s_N-\tau )}\textrm{d}\tau \Bigr \} &{} \text {if}\ t\ge s_N. \end{array}\right. } \end{aligned}$$
(18)

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:

$$\begin{aligned} p_1(t)=\frac{P\big (N_A(t)+N_P(t)=1|N_H(t)=0\big )P(N_H(t)=0)}{1-P(N_A(t)+N_P(t)=0)}. \end{aligned}$$
(19)

The expression for \(P(N_H(t)=0)\) and \(P(N_H(t)+N_P(t)=0)\) follows from Eqs. (16) and (17). The expression for \(P\big (N_A(t)+N_P(t)=1|N_H(t)=0\big )\) can be obtained from Eq. (B-27) (see “Appendix B” for details) which is

$$\begin{aligned} P(N_A(t)+N_P(t)&=1|N_H(t)=0)=\text {exp}\left\{ h(0,t)-h(1,t)\right\} \frac{\partial h(0,t)}{\partial z}\\&=\frac{G(t,z_H=0,\ z_A=0,\ z_C=1,\ z_D=1,\ z_P=0\, z_{PC}=1)}{G(t,z_H=0,\ z_A=1,\ z_C=1,\ z_D=1,\ z_P=1\, z_{PC}=1)}\frac{\partial h(0,t)}{\partial z},\\&=\frac{P(N_H(t)=N_A(t)=N_P(t)=0)}{P(N_H(t)=0)}\frac{\partial h(0,t)}{\partial z}. \end{aligned}$$

Finally, from Eq. (19),

$$\begin{aligned} p_1(t)&=\frac{P\big (N_A(t)+N_P(t)=1|N_H(t)=0\big )P(N_H(t)=0)}{1-P(N_A(t)+N_P(t)=0)},\nonumber \\&=\frac{P(N_H(t)=N_A(t)=N_P(t)=0)}{1-P(N_A(t)=N_P(t)=0)}\frac{\partial h(0,t)}{\partial z}, \end{aligned}$$
(20)

where

$$\begin{aligned} \frac{\partial h(0,t)}{\partial z}&={\left\{ \begin{array}{ll} \int _0^t \lambda (\tau )\frac{e^{-\gamma (t-\tau )}\big (1+\nu p_H(t-\tau )\big )+\nu p_A(t-\tau )}{\big [1+\nu \big (p_A(t-\tau )+p_H(t-\tau )\big )\big ]^2}\textrm{d}\tau &{} \text {if}\ t< s_1\\ \int _{s_N}^t \lambda (\tau )\frac{e^{-\gamma (t-\tau )}\big (1+\nu p_H(t-\tau )\big )+\nu p_A(t-\tau )}{\big [1+\nu \big (p_A(t-\tau )+p_H(t-\tau )\big )\big ]^2}\textrm{d}\tau \\ +\int _0^{s_1} \lambda (\tau )\frac{(1-p_{\textrm{blood}})e^{-\gamma (t-\tau )}\big (1+\nu p_H^r(t-\tau ,s_1-\tau )\big )+\nu p_A^r(t-\tau ,s_1-\tau )}{\big [1+\nu \big (p_A^r(t-\tau ,s_1-\tau )+p_H^r(t-\tau ,s_1-\tau )\big )\big ]^2}\textrm{d}\tau \\ \int _{s_{1}}^{s_2} \lambda (\tau )\frac{(1-p_{\textrm{blood}})^2e^{-\gamma (t-\tau )}\big (1+\nu p_H^r(t-\tau ,s_1-\tau ,s_2-\tau )\big )+\nu p_A^r(t-\tau ,s_1-\tau ,s_2-\tau )}{\big [1+\nu \big (p_A^r(t-\tau ,s_1-\tau ,s_N-\tau )+p_H^r(t-\tau ,s_1-\tau ,s_N-\tau )\big )\big ]^2}\textrm{d}\tau \\ +\ldots +\int _{s_{N-1}}^{s_N} \lambda (\tau )\frac{(1-p_{\textrm{blood}})^Ne^{-\gamma (t-\tau )}\big (1+\nu p_H^r(t-\tau ,s_1-\tau ,\ldots ,s_n-\tau )\big )+\nu p_A^r(t-\tau ,s_1-\tau ,\ldots ,s_N-\tau )}{\big [1+\nu \big (p_A^r(t-\tau ,s_1-\tau ,\ldots ,s_N-\tau )+p_H^r(t-\tau ,s_1-\tau ,\ldots ,s_N-\tau )\big )\big ]^2}\textrm{d}\tau &{} \text {if}\ t\ge s_N. \end{array}\right. } \end{aligned}$$

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

$$\begin{aligned} p_2(t)&=P(N_A(t)+N_P(t)=1,N_H(t)>0|P(N_A(t)+N_P(t)>0),\nonumber \\&=\frac{P\big (N_A(t)+N_P(t)=1\big )}{1-P(N_A(t)+N_P(t)=0)}-\frac{P\big (N_A(t)+N_P(t)=1|N_H(t)=0\big )P(N_H(t)=0)}{1-P(N_A(t)+N_P(t)=0)},\nonumber \\&=\frac{P\big (N_A(t)+N_P(t)=1\big )}{1-P(N_A(t)=N_P(t)=0)}-p_1(t). \end{aligned}$$
(21)

The expression \(P(N_A(t)+N_P(t)=0)=P(N_A(t)=N_P(t)=0)\) is given by Eq. (17). The expression for \(P(N_A(t)+N_P(t)=1)\) follows from Equation (81) in Mehra et al. (2022b) and is given by

$$\begin{aligned} P(N_A(t)+N_P(t)=1)=P(N_A(t)=N_P(t)=0)\frac{\partial f(0,t)}{\partial z}, \end{aligned}$$

where,

$$\begin{aligned} \frac{\partial f(0,t)}{\partial z}=&{\left\{ \begin{array}{ll} \int _0^t \lambda (\tau )\frac{e^{-\gamma (t-\tau )}+\nu p_A(t-\tau )}{[1+\nu p_A(t-\tau )]^2}\textrm{d}\tau &{} \text {if}\ t< s_1\\ \int _{s_n}^t \lambda (\tau )\frac{e^{-\gamma (t-\tau )}+\nu p_A(t-\tau )}{[1+\nu p_A(t-\tau )]^2}\textrm{d}\tau \\ +\int _0^{s_1} \lambda (\tau )\frac{(1-p_{\textrm{blood}})e^{-\gamma (t-\tau )}+\nu p_A^r(t-\tau ,s_1-\tau )}{[1+\nu p_A^r(t-\tau ,s_1-\tau )]^2}\textrm{d}\tau \\ +\int _{s_{1}}^{s_2} \lambda (\tau )\frac{(1-p_{\textrm{blood}})^2e^{-\gamma (t-\tau )}+\nu p_A^r(t-\tau ,s_1-\tau ,s_2-\tau )}{[1+\nu p_A^r(t-\tau ,s_1-\tau ,s_2-\tau )]^2}\textrm{d}\tau \\ +\ldots +\int _{s_{N-1}}^{s_N} \lambda (\tau )\frac{(1-p_{\textrm{blood}})^Ne^{-\gamma (t-\tau )}+\nu p_A^r(t-\tau ,s_1-\tau ,\ldots ,s_n-\tau )}{[1+\nu p_A^r(t-\tau ,s_1-\tau ,\ldots ,s_n-\tau )]^2}\textrm{d}\tau &{} \text {if}\ t\ge s_N. \end{array}\right. } \end{aligned}$$

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:

$$\begin{aligned} k_1(t)&=P(N_H(t)=1|N_A(t)=N_P(t)=0,N_H(t)>0)\nonumber \\&=\frac{P(N_H(t)=1|N_A(t)=N_p(t)=0)}{1-P(N_H(t)=0|N_A(t)=N_P(t)=0)}.\nonumber \\&=\frac{\text {exp}\left\{ g(0,t)-g(1,t)\right\} }{1-P(N_H(t)=0|N_A(t)=N_P(t)=0)} \frac{\partial g(0,t)}{\partial z}\nonumber \\&=\frac{P(N_H(t)=N_A(t)=N_P(t)=0)}{\big (1-P(N_H(t)=0|N_A(t)=N_P(t)=0)\big )P(N_A(t)=N_P(t)=0)} \frac{\partial g(0,t)}{\partial z}, \end{aligned}$$
(22)

where

$$\begin{aligned} \frac{\partial h(0,t)}{\partial z}=&{\left\{ \begin{array}{ll} \int _0^t \frac{\lambda (\tau )\nu p_H(t-\tau )(1-e^{-\gamma (t-\tau )})}{[1+\nu (p_H(t-\tau )+p_A(t-\tau ))]^2}\textrm{d}\tau &{} \text {if}\ t< s_1\\ \int _{s_N}^t \frac{\lambda (\tau )\nu p_H(t-\tau )(1-e^{-\gamma (t-\tau )})}{[1+\nu (p_H(t-\tau )+p_A(t-\tau ))]^2}\textrm{d}\tau \\ +\int _0^{s_1} \frac{\lambda (\tau )\nu p_H^r(t-\tau ,s_1-\tau )(1-(1-p_{\textrm{blood}})e^{-\gamma (t-\tau )})}{[1+\nu (p_H^r(t-\tau ,s_1-\tau )+p_A^r(t-\tau ,s_1-\tau ))]^2}\textrm{d}\tau \\ +\int _{s_{1}}^{s_2} \frac{\lambda (\tau )\nu p_H^r(t-\tau ,s_1-\tau ,s_2-\tau )(1-(1-p_{\textrm{blood}})^2e^{-\gamma (t-\tau )})}{[1+\nu (p_H^r(t-\tau ,s_1-\tau ,s_2-\tau )+p_A^r(t-\tau ,s_1-\tau ,s_2-\tau ))]^2}\textrm{d}\tau \\ +\ldots +\int _{s_{N-1}}^{s_N} \frac{\lambda (\tau )\nu p_H^r(t-\tau ,s_1-\tau ,\ldots ,s_N-\tau )(1-(1-p_{\textrm{blood}})^Ne^{-\gamma (t-\tau )})}{[1+\nu (p_H^r(t-\tau ,s_1-\tau ,\ldots ,s_N-\tau )+p_A^r(t-\tau ,s_1-\tau ,\ldots ,s_N-\tau ))]^2}\textrm{d}\tau &{} \text {if}\ t\ge s_N. \end{array}\right. } \end{aligned}$$
(23)

The expression for \(P(N_H(t)=1|N_A(t)=N_p(t)=0)\) follows from Equation (78) in Mehra et al. (2022b) and \(P(N_H(t)=0|N_A(t)=N_P(t)=0)\) is obtained by dividing Eq. (18) by Eq. (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:

$$\begin{aligned} k_T=\sum _{i=1}^\infty ik_i&= \Big (\frac{\mathbb {E}\left[ N_H(t)|N_A(t)=N_P(t)=0\right] }{1-P(N_H(t)=0|N_A(t)=N_P(t)=0)}\Big ) \end{aligned}$$

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:

$$\begin{aligned}&\mathbb {E}\left[ N_H(t)|N_A(t)=N_P(t)=0\right] \nonumber \\&\quad ={\left\{ \begin{array}{ll} \int _0^t \frac{\nu p_H(t-\tau )\lambda (\tau )\big (1-e^{-\gamma (t-\tau )}\big ) }{[1+\nu p_A(t-\tau )]^2 }\textrm{d}\tau &{} \text {if}\ t< s_1\\ \int _{s_{n}}^t \frac{\nu p_H(t-\tau )\lambda (\tau )\big (1-e^{-\gamma (t-\tau )}\big ) }{[1+\nu p_A(t-\tau )]^2 }\textrm{d}\tau \\ + \int _{0}^{s_1} \frac{\nu p_H^r(t-\tau ,s_1-\tau )\lambda (\tau )\big (1-(1-p_{\textrm{blood}})e^{-\gamma (t-\tau )}\big ) }{[1+\nu p_A^r(t-\tau ,s_1-\tau )]^2 }\textrm{d}\tau \\ + \int _{s_{1}}^{s_2} \frac{\nu p_H^r(t-\tau ,s_1-\tau ,s_2-\tau )\lambda (\tau )\big (1-(1-p_{\textrm{blood}})^2e^{-\gamma (t-\tau )}\big ) }{[1+\nu p_A^r(t-\tau ,s_1-\tau ,s_2-\tau )]^2 }\textrm{d}\tau \\ +\ldots + \int _{s_{N-1}}^{s_N} \frac{\nu p_H^r(t-\tau ,s_1-\tau ,\ldots ,s_n-\tau )\lambda (\tau )\big (1-(1-p_{\textrm{blood}})^Ne^{-\gamma (t-\tau )}\big ) }{[1+\nu p_A^r(t-\tau ,s_1-\tau ,\ldots ,s_N-\tau )]^2 }\textrm{d}\tau &{} \text {if}\ t\ge s_N \end{array}\right. } \end{aligned}$$
(24)

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 4th-order Runge–Kutta method, as described by Algorithm 1 in Anwar (2022). 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:

$$\begin{aligned} \begin{aligned}&{\mathop {{{\,\textrm{minimise}\,}}}\limits _{s_1, s_2,\ldots , s_N}} \quad \ Z\\&\quad \quad \quad \text {s.t.}\quad 0<s_1<s_2<\ldots <s_N,\\ \end{aligned} \end{aligned}$$

where Z is the objective function to be minimised. Based on public health relevance, 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 \Big )\),

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 and an estimated 79–96% of the total vivax infections are due to relapse, it is important to target the hypnozoite reservoir to disrupt the effect of relapse. Therefore, 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 Fig. 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:

$$\begin{aligned} \begin{aligned}&{\mathop {{{\,\textrm{minimise}\,}}}\limits _{x_1, x_2,\ldots , x_{N-1}}} \quad Z \\&\quad \quad \quad \quad \text {s.t.} \quad x_1,\ x_2,\ \ldots ,\ x_{N-1}>0\ \text {and}\ \sum x_i \le t_{\max }\\ \end{aligned} \end{aligned}$$
(25)

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 (Eq. 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:

$$\begin{aligned} \begin{aligned} {\mathop {{{\,\textrm{minimise}\,}}}\limits _{x_0, x_1,\ldots , x_{N-1}}} \quad&Z\\ \text {s.t.} \quad&\ x_1,\ \ldots ,\ x_{N-1}>0,\ x_0\ge 0 \ \text {and}\ \sum x_i \le t_{\max }\\ \end{aligned} \end{aligned}$$
(26)

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 Sect. 3.1. The effect of drug efficacy (varying \(p_{\textrm{rad}}\)) with one MDA round is presented in Sect. 3.2. We then present numerical results on the effect of multiple MDA rounds on disease prevalence (Sect. 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 (Sect. 3.4).

3.1 The Effect of a Single Round of MDA (with \(p_{\textrm{blood}}=0.9,\ p_{\textrm{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 (Lydeamore 2019). 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_{\textrm{blood}}=0.9\), \(p_{\textrm{rad}}=0.9\)); the hypnozoite reservoir size just before the MDA (Fig. 3B) becomes smaller in size (Fig. 3C; mode is 0) 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 (Fig. 3A), as we assume that the drug clears any ongoing blood-stage infections with \(90\%\) efficacy (\(p_{\textrm{blood}}=0.9\)). For liver-stage infected individuals, as the drugs are assumed to kill each hypnozoite with probability \(p_{\textrm{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 (Fig. 3D) as before MDA (Fig. 3B) eventually.

Fig. 3
figure 3

Results from multiscale model under radical cure treatment (\(p_{\textrm{blood}}=0.9\), \(p_{\textrm{rad}}=0.9\)) with a single round of MDA without seasonality. Parameters are as per Table 1. Subplot A depicts the proportion of blood-stage and liver-stage infected humans over time under treatment. The colored dashed lines indicate the times at which the hypnozoite distribution is quantified in Subplots (B–D). Hypnozoite distribution in population (obtained as per Equations (74)–(75) in Mehra et al. (2022b)) before MDA (time is indicated by the blue dashed line in Subplot A is depicted in Subplot (B). Subplot (C) and Subplot D depict the hypnozoite distribution in the population just after and 10 years after the MDA, respectively (times are indicated by green and purple dashed lines in Subplot A)

Full size image
Fig. 4
figure 4

Effect of radical cure on liver-stage infected individuals without seasonality. Subplot A depicts the proportion of liver-stage infected for different hypnozoitocidal efficacy levels (\(p_{\textrm{rad}}\)). Yellow, blue, and green lines corresponds to \(p_{\textrm{rad}}=0.9\), \(p_{\textrm{rad}}=0.95\), and \(p_{\textrm{rad}}=1\), respectively. Here \(p_{\textrm{blood}}=0.9\) for all scenarios. Subplots B, C, and D show the hypnozoite distribution within the population just after the MDA program when \(p_{\textrm{rad}}=0.95\), \(p_{\textrm{rad}}=0.99\), and \(p_{\textrm{rad}}=1\), respectively (obtained as per Equations (74)–(75), in Mehra et al. (2022b)). Other parameters are as in Table 1

Full size image
Fig. 5
figure 5

Effect of multiple rounds of MDA without seasonality. Subplot A shows the effect of four MDA rounds on prevalence over time whereas Subplot B is a snapshot from Subplot A that shows the transient dynamics during the MDA rounds. A sensitivity analysis of up to \(N=6\) MDA rounds over different steady state prevalences is illustrated in Subplots (C–D). Subplot C shows the disease prevalence 5 years after the first MDA round for up to six MDA rounds and Subplot D shows the disease prevalence 15 years after the first MDA round. The green and red asterisks in Subplots C and D are the prevalences corresponding to the green and red asterisks in Subplot (A), respectively. The intervals between MDA rounds are fixed at 30 days. Other parameters are as in Table 1

Full size image

3.2 The Effect of a Single Round of MDA, Varying \({p_{\textrm{rad}}}\)

The effect of the drug on disease transmission and the hypnozoite reservoir also changes with the efficacy of the drug (Fig. 4). Figure 4A illustrates the effect of varying efficacies of the hypnozoicidal drug (i.e. \(p_{\textrm{rad}}\)) on liver-stage infected proportions. Figure 4B illustrates the hypnozoite distribution just after the application of MDA when \(p_{\textrm{blood}}=p_{\textrm{rad}}=0.9\) and Fig. 4C illustrates the hypnozoite distribution when \(p_{\textrm{blood}}=0.9,\ p_{\textrm{rad}}=0.95\). If the hypnozoicidal drug were 100% effective (that is, \(p_{\textrm{rad}}=1\)) then all liver-stage infected individuals would recover (green line in Fig. 4A). In the case of \(p_{\textrm{rad}}=1\), the hypnozoite reservoir within the human population would be completely cleared (Fig. 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. (Fig. 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 (Fig. 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 Fig. 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 (Newby 2015). The effect of four successive MDA rounds is clearly visible in Fig. 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, D illustrate the sensitivity analysis of up to six MDA rounds at different assumed prevalences (20–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 Eq. (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 Fig. 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 Fig. 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 Fig. 6D. The dashed vertical lines in Fig. 6C, 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 (Fig. 6B). The left vertical axis in Fig. 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–60%), the optimisation problem with the objective function \(Z_1\) suggests an interval of around 480 days between the two MDA rounds.

Fig. 6
figure 6

Effect of two rounds of optimally timed MDA. Subplot A shows the impact of varying intervals between two MDAs on the objective function \(Z_1\) (using a starting steady-state disease prevalence of 20%) where the minimum objective value obtained from the optimisation problem (Eq. 25) is when the interval is 34.7 days. Subplot C depicts the objective function, \(Z_1\), over time, before and after the MDAs when using the optimal interval of \(\approx 35\) days. Subplot D illustrates the transient disease dynamics corresponding to two optimally timed MDAs (time for the first MDA is arbitrary), separated in time by \(\approx 35\) days. Finally, subplot B illustrates the optimal interval for different disease prevalences (right vertical axis) by varying the mosquito ratio (left vertical axis), m corresponding to the objective function \(Z_1\). These optimal intervals are for two MDA rounds found by solving Eq. (25) where the red rectangle shows the \(\approx 35\) day optimal interval for the 20% steady state disease prevalence used in Subplot (A), (C), and (D). All parameters are as in Table 1

Full size image
Fig. 7
figure 7

Sensitivity analysis of two and three rounds of optimal MDA intervals over different disease prevalence (right vertical axis) without seasonality. Prevalence is varied by varying mosquito ratio, m (left vertical axis). The optimal interval between two and three rounds of MDA for the objective function \(Z_1\) (Subplot A, D, respectively), \(Z_2\) with \(W_h=1,\ W_m=0\) (Subplot B, E, respectively) and \(Z_2\) with \(W_h=1,\ W_m=1\) (Subplot C, F, respectively) are shown. All parameters are as in Table 1

Full size image

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 (Fig. 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\) (Fig. 7B, C).

The optimal intervals for three MDA rounds depend on both m, hence prevalence, and the choice of the objective function (Fig. 7D–F). If three MDA rounds are considered, the optimisation problem (Eq. 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 Fig. 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 (Fig. 8D–F). The seasonal amplitude, \(\eta \), is thought to play an important role in intervention strategies (Selvaraj et al. 2018); 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 Fig. 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 Fig. 8D–F indicate the optimal intervals corresponding to Fig. 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)\) (Fig. 8E), while with \(Z_2\, (W_h=W_m=1)\), the algorithm suggests a similar interval as for low prevalences (Fig. 8F). With objective function \(Z_1\), the interval between the two MDA rounds is also very similar for different annual mean prevalences (Fig. 8D).

Fig. 8
figure 8

Effect of two rounds of optimally timed MDA with mosquito seasonality. Subplots A–C depict the impact of optimal MDA on disease prevalence (annual mean disease prevalence before MDA of 54.9%) with objective function \(Z_1\), \(Z_2\, (W_h=1,\ W_m=0)\), and \(Z_2\, (W_h=W_m=1)\), respectively. The solid vertical line indicates the time when the prevalence reaches a peak before the initial MDA. The dashed vertical lines indicate the optimal times for the MDA rounds. Subplots D–F depict the sensitivity analysis over different annual mean disease prevalences with the objective function \(Z_1\), the objective function \(Z_2\) with \(W_h=1, W_m=0\), and objective function \(Z_2\) with \(W_h=W_m=1\), respectively. All parameters are as in Table 1

Full size image

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 Fig. 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 Fig. 9D–F indicate the optimal intervals corresponding to Fig. 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.

Fig. 9
figure 9

Effect of three rounds of optimally timed MDA with mosquito seasonality. Subplots A–C depict the impact of three optimally timed MDAs on objective functions \(Z_1\), \(Z_2\,(W_h=1,\ W_m=0)\), and \(Z_2\, (W_h=W_m=1)\), respectively (note the change in X-axis between Subplots A, B, C) (annual mean disease prevalence before MDA of 49.9%). Subplots D-F depict a sensitivity analysis over different annual mean disease prevalences (right vertical axis) obtained by varying initial mosquito ratio (left vertical axis) with objective function \(Z_1\), the objective function \(Z_2\) with \(W_h=1, W_m=0\), and the objective function \(Z_2\) with \(W_h=W_m=1\), respectively. All parameters are as in Table 1

Full size image
Fig. 10
figure 10

Effect of change in model parameters on the optimal interval (without seasonality) between 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. Subplots A–E illustrate the impact of varying m (number of mosquitoes per human), \(\alpha \) (hypnozoite activation rate), \(\mu \) (hypnozoite death rate), \(\gamma \) (natural recovery rate), and \(\nu \) (average number of hypnozoites per bite) on the optimal interval, respectively. The colorbars in each subplot illustrate the equilibrium prevalence corresponding to the parameters before the first MDA was implemented. The red arrows in Subplots A–E indicate the baseline parameters in Table 1 and optimal interval when prevalence is around \(40\%\) (see Fig. 6B as a reference). Parameter ranges for Subplot A–E are as in Table 1

Full size image

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 Figs. 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 Fig. 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 depicts the effect of varying \(m,\ \alpha ,\ \mu ,\ \gamma \), and \(\nu \) on the optimal intervals, respectively. The color of the bars in Fig. 10A–E depicts the equilibrium prevalence corresponding to the parameter value. The red arrows in Fig. 10A–E depict the baseline parameters in Table 1 that generate a prevalence of \(40\%\) as shown in Fig. 6B. As illustrated in Fig. 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_0abI_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\%\) (Fig. 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 Figs. 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_{\textrm{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 Eq. C-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_TL\)). 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 (Figs. 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\) (Fig. 12A) and for objective function \(Z_2\) with \(W_h=1, W_m=0\) (Fig. 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.

Fig. 11
figure 11

Effect of change in two model parameters (the hypnozoite activation rate, \(\alpha \), and the number of mosquitoes per human, m) on the optimal interval (without seasonality) between two MDA rounds with the objective functions \(Z_1\) and \(Z_2\) with \(W_h=1, W_m=0\). Subplot A is a violin plot that illustrates the optimal interval corresponding to the objective function \(Z_1\), whereas Subplot B illustrates the optimal interval corresponding to the objective function \(Z_2\) with \(W_h=1, W_m=0\). In both cases, for a given value of \(\alpha \) (outer color of the scatter points), we choose the parameter value m (inner color of the scatter points) so that the steady-state prevalence is \(30\%\). The colorbars on the bottom illustrate the value of m and \(\alpha \), respectively. All other parameters are as in Table 1

Full size image
Fig. 12
figure 12

Effect of change in two model parameters (the hypnozoite activation rate, \(\alpha \), and the initial number of mosquitoes per human, \(m_0\)) on the optimal interval (with seasonality) between two MDA rounds with the objective functions \(Z_1\) and \(Z_2\) with \(W_h=1, W_m=0\). Subplot A illustrates the optimal intervals \(x_0\) and \(x_1\) corresponding to the objective function \(Z_1\), whereas Subplot B illustrates the optimal interval corresponding to the objective function \(Z_2\) with \(W_h=1, W_m=0\). In both cases, for a given value of \(\alpha \) (outer color of the scatter points), we choose the parameter value \(m_0\) (inner color of the scatter points) so that the annual mean prevalence is \(\approx 30\%\). The colorbars on the bottom illustrate the value of \(m_0\) and \(\alpha \), respectively. All other parameters are as in Table 1

Full size image

4 Conclusions and 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 (Hsiang 2013). 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 (Anwar 2022) by integrating treatment into the model with multiple MDA rounds accounting for superinfection. We have extended Mehra et al. (2022b) 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 Mehra (2022). 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 (Fig. 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, Figs. 10 and 11). However, when there is seasonal variation in the mosquito population in the environment, the optimal intervals are very similar (Fig. 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 (Maude 2012; Greenwood 2008; Hsiang 2013), treatment coverage (the proportion of the population who are treated) drives the overall effectiveness of the MDA (Finn Timothy 2020; Dyson 2017). The higher the coverage of MDA, the more transmission will be reduced and the closer we will reach towards elimination (Slater 2014). For simplicity, here we assume 100% treatment coverage in our model which is difficult to achieve in reality due to various factors (Agboraw 2021; Finn Timothy 2020). However, if and when this assumption of 100% coverage is relaxed and each MDA round has a certain coverage, the correlation of coverage will be important as systematic non-adherence can greatly undermine the success of the MDA program and be ineffective in disrupting onward transmission (Dyson 2017; Rock 2015; Plaisier 2000). We have also assumed that all drugs (both blood-stage and liver-stage) are 90% effective, unless specified otherwise. 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 (Huber 2021; Nelwan 2015; Llanos-Cuentas 2014). 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 (Figs. 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 h to 28 days (Jittamala 2015; Schlagenhauf 2019) 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 (Allen and Burgin 2000). 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 (Huber 2021). 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 (Mehra et al. 2022b). 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_{\textrm{blood}}\).

Although being an effective intervention strategy, MDA has some disadvantages, especially in terms of drug resistance (Zuber and Takala-Harrison 2018; Commons 2018). 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 (World Health Organization 2020). 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 (Howes Rosalind 2012; Watson 2018). 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 (Keeling and Rohani 2011; Greenhalgh et al. 2015). 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 (Le Menach 2007; Price 2020; Newby 2015). 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 (Zuber and Takala-Harrison 2018). 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.