3.1 Model frameworks

In infectious disease dynamics, modelling frameworks typically involve a combination of mathematical models, statistical analyses, and computer simulations that aim to capture the complex dynamics of disease transmission. The Ross-Macdonald model [57], a compartmental model initially developed to describe malaria transmission dynamics, has been widely used as a modelling framework for P. vivax transmission. This modelling approach has been adapted to investigate a range of vector-borne infectious diseases, and has helped inform public health policies and intervention strategies. The first mathematical model describing P. vivax transmission was introduced by—to the best of our knowledge—Zoysa et al. (1988) [38] in a Ross-Macdonald style modelling approach. Following this, many models have now been developed.

Out of the 48 studies identified that incorporated a P. vivax transmission model, 38 (79%) utilised a deterministic and differential equation (compartmental) framework [4, 20, 24, 28–30, 32–35, 38, 39, 56, 58–82] and nine (19%) used a stochastic and agent-based framework [21, 54, 55, 83–88] (Fig 2). Only one study (2%) used both deterministic and stochastic frameworks [20] to model P. vivax transmission. Robinson et al. (2015) [20] developed the model in a deterministic framework but implemented a stochastic version of the model as a continuous-time Markov chain. For simplicity, we categorise this model as deterministic in Fig 2. Deterministic models are often the first choice amongst modellers due to their simplicity in comparison to stochastic models. Deterministic models are useful for understanding disease dynamics in large populations. Stochastic models provide more realistic and accurate representations of complex systems when dealing with small population sizes or low disease prevalence, as they can account for the randomness and variability observed in real life [89, 90].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. A summary of the 48 P. vivax transmission models currently available in the literature (published as of May 21, 2023) [4, 20, 21, 24, 28–30, 32–35, 38, 39, 54–56, 58–88].

Related models (either modified or motivated by) are connected with a dashed line. Similar/same models are connected with a solid line. The coloured boxes represent key features incorporated in the models (see legend). The hexagonal boxes with the same name represent that the model was also implemented in other frameworks. The timescale (non-linear) is shown on the left.

https://doi.org/10.1371/journal.pcbi.1011931.g002

In contrast to the compartmental differential equation framework, agent-based models represent a system as a collection of individual agents that interact with each other based on a set of rules or behaviours [91, 92]. The main difference between compartmental and agent-based modelling frameworks is that a compartmental model uses aggregate variables or compartments to represent the system, while agent-based models use individuals (agents) [91]. Out of the eight studies that used an agent-based model to capture the dynamics of P. vivax transmission, only two studies modelled both the human and mosquito populations as agents [83, 84]. The other agent-based models modelled the mosquito populations as a deterministic compartmental process, such that they combined ordinary differential equations for mosquitoes with an agent-based model for humans [21, 54, 55, 85–88]. Modelling mosquito dynamics as a deterministic process is an approximate strategy if the size of the mosquito population is very large and P. vivax is not near elimination. In this case, the average behaviour of the stochastic dynamics agrees with those of a deterministic process [92–94]. The actual behaviour of the system depends on the interactions between individuals and mosquitoes, instead of averages. Treating the mosquito compartment as a deterministic process means that modelling elimination is not straightforward, as there will always be some non-zero number of infectious mosquitoes remaining that can trigger an infection in humans again [87].

Environmental features, ecology, and mosquito habitat locations were explicitly included when modelling malaria spread in the agent-based models that modelled both the human and mosquito population as agents [83, 84]. The most recent agent-based models modelling P. vivax dynamics [21, 54, 85, 86, 88] have evolved from a model introduced by White et al. (2018) [55]. The White et al. (2018) [55] model has been adapted to capture disease epidemiology in particular geographical settings [85], and to study the impact of different interventions (drugs or vaccination) [21, 86, 88].

While agent-based models have many advantages, their use poses several challenges. One of the main challenges of agent-based models is the difficulty in parameterising and calibrating the model, given the large number of agents and their interactions [95–97]. For example, despite being an agent-based model, parameterisation is done using an ordinary differential equation system that describes the process in several models [21, 54, 55, 85]. Despite these challenges in parameterisation, agent-based models also often offer a more intuitive representation of epidemiological processes. The computational demands of agent-based models can be challenging [98], although with improving computer technology, this has become less of a concern [99].

3.2 Population-level multiscale models

Multiscale disease modelling incorporates at least two interacting scales and provides insights into disease dynamics across these scales that cannot be obtained from a single scale alone [100]. Here we only focus on within-host population models as ‘multiscale models’. For P. vivax, multiscale modelling approaches can incorporate the complex hypnozoite dynamics and their relapse effects on onward disease transmission. Most models in the existing literature only capture the population-level impact of P. vivax (boxes with a light lime green border in Fig 2). Few models capture both within-host and population-level impacts (boxes with a strong blue border in Fig 2) [21, 32, 35, 54, 55, 67, 85, 86]. The very first multiscale model for P. vivax transmission was developed by White et al. (2018) [55], and modelled the within-host hypnozoite dynamics using an agent-based model that considered heterogeneity in exposure to mosquito bites. This built on White et al. (2014) [32], which was the first to develop a within-host model that captured the dynamics of P. vivax hypnozoites. This multiscale model considered the variability in the size of hypnozoite inoculum across each mosquito bite and was subsequently used to parameterise a separate transmission model that captured the entire structure of the hypnozoite reservoir [55]. The White et al. (2014) [32] within-host model for temperate settings assumed collective dormancy. This means that the hypnozoites established by each mosquito bite progress through the dormancy states as a group or batch. This assumption may be biologically unrealistic due to the independence of individual hypnozoite activation and clearance dynamics within liver cells [101]. The other within-host models that were adapted from White et al. (2018) [55] applied the same assumption regarding batch hypnozoite behaviour [21, 54, 85, 86].

Recent work by Mehra et al. (2020) [101] relaxed the collective dormancy assumption. This enabled them to characterise the long-latency period of hypnozoite dynamics (a period of latency prior to hypnozoite activation) modelled (light purple bordered box in Fig 2) in White et al. (2014) [32] in analytical form. Later work by Mehra and colleagues embedded the activation-clearance model governing a single hypnozoite in an epidemiological framework [34]. This framework accounts for successive mosquito bites, where each bite can simultaneously establish multiple hypnozoites [34, 102], and explores the epidemiological consequence of radical cure treatment on a single individual. Anwar et al. (2022) [35] have since developed a multiscale model motivated by White et al. (2014) [32] by embedding the framework of Mehra et al. (2022) [34] for short-latency hypnozoites (deriving the relapse rate by averaging the distribution of hypnozoite burden, which is dependent on the force of reinfection) into a simple population-level model that provides key insights into both within-host level and population level dynamics. The within-host and population models were coupled at each time step (thus producing a multiscale model) to incorporate key parameters that describe the hypnozoite dynamics. This multiscale model can provide the hypnozoite distributions within the population and, more importantly, reduces the infinite compartmental structure of White et al. (2014) [32] into three compartments and relaxes the artificial truncation needed in White et al. (2014) [32] for numerical simulation. Mehra et al. (2022) [36] proposed an alternative approach, constructing a Markov population process to couple host and vector dynamics whilst accounting for (short-latency) hypnozoite accrual and superinfection as per the within-host framework proposed in Mehra et al. (2022) [34]. In the infinite population size limit, Mehra et al. (2022) [36] recovered a functional law of large numbers for this Markov population process, comprising an infinite compartment deterministic model. This infinite compartment model was then reduced into a system of integrodifferential equations based on the expected prevalence of blood-stage infection derived at the within-host scale [34]. This construction yielded population-level distributions of superinfection and hypnozoite burden, and has been generalised to allow for additional complexity, such as long-latency hypnozoites and immunity [36].

3.3 Hypnozoite dynamics and variation

The eradication of P. vivax is challenging due to the presence of the hypnozoite reservoir, which is undetectable and causes new infections long after the initial infection. In developing the first mathematical model for P. vivax, Zoysa et al. (1991) were also the first to model the effect of hypnozoite relapse on P. vivax transmission [39]. Since most P. vivax blood-stage infections are due to the reactivation of hypnozoites rather than new primary infections (in the absence of anti-relapse treatment), it is crucial that mathematical models incorporate the size of the hypnozoite reservoir [103–107]. Zoysa et al. (1991) [39] assumed that the transmission dynamics could be accounted for by modelling a hypnozoite reservoir of size two (to account for up to two relapses). This assumption was later followed by Fujita et al. (2006) [30]. In reality, the average size of the hypnozoite reservoir is likely to be more than two in endemic settings, particularly those with high transmission intensity [33]. Despite having the relapse characteristic that makes P. vivax parasites unique, Nah et al. (2010) [63] and Aldila et al. (2021) [77] did not incorporate relapses in their P. vivax transmission model. In their model, individuals did not harbour hypnozoites when infected with P. vivax and hence did not experience relapse after recovery from blood-stage infection.

Modelling hypnozoite dynamics and relapses can be addressed with varying levels of complexity. The common modelling approaches to embedding hypnozoite dynamics into transmission models are the (i) binary hypnozoite model, (ii) batch hypnozoite model, and (iii) density hypnozoite model. The binary hypnozoite model is the simplest one, assuming that an individual can be either infected with or without hypnozoites [33]. Individuals experience relapse and clear hypnozoites at a constant rate. The batch hypnozoite model assumes that each infectious mosquito bite contributes to a batch of hypnozoites where each batch of hypnozoite is independent of the other. Individuals can experience multiple relapses from each batch, and the overall relapse rate depends on the total batches of hypnozoites present. Hypnozoites are subject to clearance at a constant rate. This reduces one batch of hypnozoites to a single set of dynamics, which is still truncated at a maximum number of batches [55]. The density hypnozoite model is the most complex but most biologically accurate among the approaches. In this approach, hypnozoite numbers within individuals are explicitly modelled, where each infectious mosquito bite can contribute a number of hypnozoites. Hypnozoites can either die or activate at a constant rate, and the dynamics of each hypnozoite are independent and identically distributed [34, 55, 101].

Most P. vivax transmission models consider the hypnozoite reservoir as a single compartment (i.e. binary approach), rather than explicitly accounting for a variable number of hypnozoites in the reservoir [20, 24, 28, 29, 33, 56, 58–66, 68–76, 78–81, 83, 84]. Only a handful of models account for the variability in hypnozoite inoculation across mosquito bites (i.e. hypnozoite density model, boxes with a bright pink border in Fig 2) [4, 32, 34, 35]. If the size of the hypnozoite reservoir is modelled explicitly, the number of compartments in the model increases substantially. The very first model that accounted for the variation in hypnozoites across mosquito bites was introduced by—to the best of our knowledge—White et al. (2014) [32] for a short-latency strain (where hypnozoites can activate immediately after establishment). To account for the variation of hypnozoites across bites, White et al. (2014) modelled a system with an infinite number of compartments to represent individuals with different numbers of hypnozoites. In practice, this is truncated at 2(L_max + 1) ordinary differential equations (for human population only), where L_max is the maximum number of hypnozoites considered. In their model, the hypnozoite reservoir within individuals increases with new infectious bites and decreases with both activation and death of hypnozoites. This infinite compartmental system makes the model very complex, particularly when other important structures must also be incorporated, such as individual heterogeneity in bite exposure. An agent-based model later developed by White et al. (2018) [55], and other models that utilise this agent-based model, consider variation in hypnozoites within individuals, but do not account for the variability in hypnozoites across mosquito bites [21, 54, 85, 86]. Furthermore, instead of explicitly modelling hypnozoites independently, they impose the batch hypnozoite model.

The multiscale model developed by Anwar et al. (2022) [35] accounted for the variation of hypnozoites dynamics across bites. Unlike the White et al. (2014) [32] model, Anwar et al. (2022) only utilised three compartments at the population level by embedding the within-host model (short-latency) developed by Mehra et al. (2022) [34] as a system of integrodifferential equations. This relaxes the artificial truncation for the maximum number of hypnozoites used within the White et al. (2014) [32] model. Under a constant force of reinfection, Anwar et al. (2022) [35] analytically proved that the multiscale model [35] exhibits an identical steady-state hypnozoite distribution as the infinite ordinary differential equation model structure in White et al. (2014) [32]. The advantage of the multiscale model by Anwar et al. (2022) [35] is that the population-level component is considerably simpler than the 2(L_max + 1) ordinary differential equations of White et al. (2014) [32]. The transmission models proposed by Mehra et al. (2022) [36] likewise account for variation in hypnozoite batch sizes, with Mehra et al. (2022) [36] additionally accommodating long-latency hypnozoite dynamics. The models of Mehra et al. [36] are formulated as systems of integrodifferential equations, informed by the within-host framework of Mehra et al. (2022) [34]. The analyses of Anwar and Mehra et al. [35, 36] provided insights into hypnozoite dynamics (e.g. the average size of a hypnozoite reservoir within the population and the average relapse rate), in addition to disease dynamics.

Furthermore, hypnozoites have two different relapse patterns: (i) fast relapsing in tropical regions and (ii) slow relapsing in temperate regions [25]. Depending on the setting the model is being developed for, these relapse patterns may also need to be taken into account. In most of the models, the tropical relapse pattern is achieved with a shorter latent period [4, 20, 24, 28–30, 35, 38, 39, 56, 58–62, 64–66, 68–72, 74–76, 78–82]. Only a few models have considered both temperate and tropical relapse patterns [32–34]. White et al. (2014) [32] and Mehra et al. (2022) [34] explicitly accounted for both temperate and tropical relapse patterns by modelling the hypnozoite dynamics. In their temperate model, after hypnozoites were established they undergo a dormancy phase, which is enforced through additional model compartments. In the dormancy phase, the hypnozoites do not activate and can only die. To accommodate the tropical pattern, this dormancy phase is ignored [32, 34]. White et al. (2016) [33] later captured the temperate relapse pattern in a separate transmission model (with a binary hypnozoite model) by considering additional compartments where hypnozoite-positive individuals stay for some duration of time before they experience relapse. Kim et al. (2020) [73] modelled only the temperate relapse pattern in Korea through a survival function to obtain the time to relapse.

3.4 Superinfection

Superinfection with malaria is a common phenomenon, especially in high transmission settings, and can be defined as when an individual has more than one blood-stage infection with the same malaria-causing parasite species at a given time [108]. For P. falciparum malaria, when an infected individual (primary infection) receives a second infectious mosquito bite, they can become infected with two different parasite broods. In reference to P. vivax malaria, individuals can harbour hypnozoites in the liver even after they recover from a primary infection. Therefore, relapsing hypnozoites can provide another pathway to superinfection for individuals infected with P. vivax [108, 109].

When modelling P. vivax dynamics, it is important to consider the impact of superinfection on recovery and transmission, especially in high transmission settings, as the abundance of mosquitoes and the contribution of hypnozoite activation can frequently trigger superinfection. Superinfection can potentially delay recovery from infection [110, 111]. Most of the literature that incorporates superinfection in P. vivax transmission models (boxes with a brown border in Fig 2) [21, 28–30, 32, 55, 72, 80, 85, 86] does so via the recovery rate [28–30, 32]. The superinfection phenomenon was first introduced into malaria models by Macdonald (1950) [112], who assumed “The existence of infection is no barrier to superinfection, so that two or more broods or organisms may flourish side by side”. In the malaria modelling literature, it has been assumed that each brood could be cleared independently at a constant rate. Following this assumption, Dietz et al. (1974) [111] proposed a recovery rate under superinfection for P. falciparum malaria, derived at equilibrium under a constant force of reinfection. This form of the recovery rate was adopted in most studies that included superinfection via the recovery rate. This approach is straightforward when hypnozoites are integrated into the model as a binary state (i.e. an individual either has or does not have hypnozoites) [28–30]. Since White et al. (2014) [32] accounts for the variation of hypnozoites, they modified the recovery rate proposed by Dietz et al. (1974) [111] to account for the additional burden of hypnozoites; however, Mehra et al. (2022) [36] argued that this modified recovery rate does not hold in the presence of hypnozoite accrual.

Generally, there are two approaches when incorporating superinfection: (i) using a corrected recovery rate that explicitly accounts for the history of past infections in the population and hypnozoite accrual dynamics [4, 36, 113] and (ii) coupling the prevalence of blood-stage infection (derived under a within-host model that accounts for superinfection) directly to the proportion of infected mosquitoes [36]. The within-host model of Mehra et al. (2022) [34] included superinfection, with each blood-stage infection (whether primary or relapse) being cleared independently, while the population-level model developed by Anwar et al. (2022) [35], which built on work of Mehra et al. (2022) [34], did not incorporate superinfection. A correction to account for superinfection, based on the recovery rate formulated by Nåsell et al. (2013) [113], was proposed in Mehra et al. (2023) [36] and incorporated in later work by Anwar et al. (2023) [4].

Superinfection was incorporated in later work, where it was assumed that different batches of hypnozoites originated from different mosquito bites [21, 54, 55, 86]. Silal et al. (2019) [72] assumed that superinfection increased the severity of the disease. That is, individuals will transition from lower to higher severity classes with a certain probability due to multiple infections. The only other study incorporating a superinfection-like phenomenon was Aldila et al. (2021), who modelled P. vivax and P. falciparum co-infection and assumed that P. vivax dominates P. falciparum [77], which does not closely resemble the definition of superinfection. This study assumed that if an individual was currently infected with P. falciparum, they would become infected with P. vivax if they received an infectious bite from a mosquito that was infected with P. vivax. The assumption that P. vivax parasites dominate P. falciparum results in the individual being infected with only P. vivax, which is not supported by the empirical biological evidence that shows that the parasitaemic load is much higher for P. falciparum [114]. Accordingly, it may not be reasonable to consider this to be a valid model of superinfections.

3.5 P. vivax and P. falciparum co-infection

Within the Asia-Pacific region, the horn of Africa, and South America, both P. vivax and P. falciparum parasites are common [72, 87]. For example, in 2019 in Cambodia, co-infection with both P. vivax and P. falciparum accounted for about 17% of malaria cases [115]. In co-endemic regions, P. falciparum infections are often followed by P. vivax infection, giving rise to the hypothesis that P. falciparum infections trigger P. vivax hypnozoite activation [72, 116–118]. The high risk of P. vivax parasitaemia after P. falciparum infection is possibly related to reactivation of hypnozoites [119–121]. Hypnozoites may be activated when P. falciparum parasites have been introduced into the body [122] or when the individual is exposed to Anopheles specific proteins [27]. This increased risk of P. vivax blood-stage infection following a P. falciparum infection could alternatively be explained by spatial or demographic heterogeneity in exposure and thus infection risk. Individuals either living in areas where both P. vivax and P. falciparum are highly prevalent or those that engage in an activity bringing them in frequent contact with infected mosquitoes (e.g. forest work) are more likely to be exposed to both parasites than the average person. Having a P. falciparum episode indicates the person has recently been exposed to infectious mosquito bites and is thus likely to have hypnozoites from previous exposure events (that may be triggered or activated spontaneously) or acquire a new primary P. vivax infections following recovery from P. falciparum infection [123–125]. The lack of diagnostics to differentiate primary infections and relapses further complicates determining when an individual is infected with P. vivax hypnozoites. This makes it challenging to disentangle whether P. falciparum infections cause relapses through the reactivation of hypnozoites.

It is also not yet clearly understood how P. vivax and P. falciparum interact, if they compete within the host or if one species causes some, if any, protection against the other [126, 127]. A systematic review and meta-analysis showed that mixed infections (P. falciparum and P. vivax) can often cause a high rate of severe infection regardless of infection order [128]. This evidence was in contrast to a previous study which suggested that severe mixed infections were more likely to happen when P. vivax infection occurred on top of an existing P. falciparum infection (i.e. superinfection), whereas the reverse scenario, P. falciparum infection on top of an existing P. vivax infection, were more likely to result in a lower risk of severe malaria [129]. Furthermore, there is likely ascertainment bias associated with mixed infections in areas with co-circulating parasite strains, as efforts might be biased towards P. falciparum detection [130]. This is likely to be particularly common during episodes of clinical malaria when parasitaemia of one species greatly exceeds the other, and the innate host immune response may suppress both infections. Gaining a better understanding of these cross-species interactions and adjusting accounting for this co-existing phenomenon in the co-endemic region will require multi-species transmission models. Only a handful of mathematical models included both these Plasmodium species [56, 61, 62, 72, 77, 83]. While both P. vivax and P. falciparum species are included in a single model by Aldila et al. (2021) [77], this model did not account for P. vivax relapses. Five studies included both species but used two independent models for each species, which did not allow for interactions between species [56, 61, 62, 77, 83].

Whether it is important to model species interaction depends on the particular geographical setting. If both parasites are co-endemic in a setting, and the research question being considered relates to both species, then it may be important to use a model that can capture the interactions between the parasite species [72, 87, 116]. To the best of our knowledge, the first model that accounts for the interaction between both species was developed by Silal et al. (2019) [72]. In this study, a separate model (deterministic, meta-population) for both species was proposed, and these two models were entangled at each time step to incorporate interactions between the species, including treatment, triggering, and masking (non-P. falciparum infections are misdiagnosed as P. falciparum). Following this work, the first agent-based model transmission model accounting for both P. vivax and P. falciparum infections and treatment was developed by Walker et al. (2023) [87]. This model had reduced complexity compared with Silal et al.’s (2019) co-infection model, but used many of the same parameter values [72] (co-infection models shown with a vivid orange bordered box in Fig 2).

3.6 Immunity

Immunity against disease acquired through infection is usually referred to as adaptive immunity, and the primary function of adaptive immunity is to destroy foreign pathogens [131, 132]. Naturally acquired immunity to malaria is characterised by relatively rapid acquisition of immunity against severe disease and a more gradual establishment of immunity against uncomplicated malaria, while sterile immunity against infections is never achieved [133–136]. In co-endemic areas, clinical immunity to P. vivax is more rapidly acquired than that due to P. falciparum [135].

How immunity is accounted for in mathematical models of malaria varies since different models consider different types of immunity. For example, immunity against new infections, immunity against severe malaria, anti-parasite immunity (i.e. the ability to control parasite density upon infection), clinical immunity (i.e. protection against clinical disease), and transmission-blocking immunity (i.e. reducing the probability of parasite transmission to mosquitoes). Immunity against new infections and severe malaria is assumed to be acquired through infection. This reduces the probability of reinfection from an infectious mosquito bite and has been modelled using up to two immunity levels [38, 39]. This type of immunity is assumed to be boosted by infection [133]. Acquiring some partial immunity (i.e. some degree of protection against malaria) following infection that wanes over time, is most common among published models [56, 58–62, 64, 69, 70, 79, 84]. Some assumptions regarding immunity include that, if treated, individuals acquire some level of immunity that reduces the probability of reinfection (i.e. gain immunity against new infection) and that this wanes over time [32, 65]. The assumption regarding permanent immunity against malaria is not considered valid, as immunity often wanes rapidly when immune adults leave malaria-endemic regions [137]. Despite this, two models assumed that recovered individuals become permanently immune to P. vivax [68, 77]. Another study assumed that only a fixed proportion of individuals are immune against P. vivax rather than explicitly incorporating immunity into the model [83]. Strategies for incorporating immunity into P. vivax transmission models thus widely vary, where some assumptions are more realistic and appropriate than others.

Individuals who have not previously experienced malaria infection almost invariably become infected when first exposed to infectious mosquito bites, as immunity against malaria has not yet developed [137]. Repeated exposure to infectious bites will still likely result in infection, though these individuals may be protected against severe malaria or death [137]. Silal et al. (2019) [72] applied the opposite assumption and assumed that repeated exposure to infectious bites would likely result in severe infection however, with less probability per bite. With increasing exposure, naturally acquired immunity will also give some level of protection against symptomatic malaria. Adults living in endemic areas are more likely to have developed protective immunity compared to children due to repeated exposure over their lifetime. Adults living in endemic areas are likely to have experienced substantially more infectious mosquito bites compared to children due to age (and therefore lengthened opportunity to acquire infectious mosquito bites), greater skin surface area, and more time spent outside in environments with a higher prevalence of mosquitoes [55, 138, 139].

Immunity should be considered when using mathematical models to capture underlying disease dynamics. The assumption regarding immunity varies among models (boxes with a purple border in Fig 2). The only model that explicitly accounts for the acquisition of immunity that increases with new bites was developed by White et al. (2018) [55]. The assumption in regards to both anti-parasite immunity (ability to reduce parasite density upon infection) and clinical immunity (protection against clinical disease) depends on age and exposure to mosquito bites which is modelled using partial differential equations [55]. They also assumed that children acquired immunity through their birth parent’s immunity, which then decayed exponentially from birth. Models that were adapted from Whiteet al. (2018) [55] also allow for the acquisition of immunity [21, 54, 85, 86]. However, the immunity acquired from a primary infection may protect more strongly against relapses (which are genetically related to the primary infection) than against a new, genetically distinct primary infection. That is, hypnozoites established from an infectious bite, when reactivated, may be less likely to cause clinical infection. This is because the parasites could be genetically identical or related, which could elicit a more protective immune response due to familiarity with the primary infection [117, 140]. Thus, relapses from the same batch of hypnozoites may only cause asymptomatic infections. Despite this, no models to date have fully accounted for the relationship between relapse and immunity. Model assumptions regarding the acquisition of immunity may be too simple to capture the true underlying biology and dynamics.

Furthermore, as repeated exposure to infectious bites or relapse from the same batch of hypnozoites might cause only asymptomatic infections, this reservoir of asymptomatic infections in addition to the hypnozoite reservoir might greatly affect the overall transmission dynamics. This is recognised as an obstacle in malaria elimination [141, 142]. Asymptomatic malaria infections have received less research attention than symptomatic infections in major studies [143, 144] and even when developing mathematical models. Only a few models have explicitly considered asymptomatic infections [55, 72, 74, 87].

3.7 Effect of interventions for malaria control

In most of the models included in this review, it was assumed that treatment would be targeted towards infected individuals [21, 30, 32, 64–66, 70, 79, 80, 83, 85, 145], but a range of other interventions can contribute to malaria control. The impact of a pre-erythrocytic vaccine on reducing P. vivax transmission has been studied and hypothesised that even with low efficacy, it can substantially reduce transmission [67]. The impact of radical cure treatment on the infected population has also been studied in models [21, 30, 32, 64, 85] where only a handful of models incorporate G6PD testing [21, 64, 85]. Furthermore, 11 studies (23%) evaluated the effect of MDA on disease transmission, despite few national control programs considering MDA for P. vivax control [4, 20, 28, 29, 34, 54–56, 72, 87, 88] (boxes with a dark green border in Fig 2). Since MDA is recommended as an important tool to reduce asymptomatic P. falciparum infection, it is also likely to be of great importance for P. vivax elimination [20, 146, 147]. One study examined the effect of multiple MDAs and MSaTs (up to two rounds) with different drug combinations (blood-stage drug only, blood-stage drug and primaquine, or blood-stage drug and tafenoquine), finding that MDA with tafenoquine following G6PD screening could significantly reduce transmission compared to MSaT, given that no tools were available at the time to identify individuals with hypnozoites [20]. The effect of long-lasting insecticide nets along with MDA was studied using an agent-based model in Papua New Guinea, where the model predicted that MDA could reduce P. vivax transmission by between 58% and 86% [55]. The same agent-based model was later used to investigate the effect of multiple treatment strategies, including MDA, MSaT with light microscopy detection of blood-stage parasitemia, and P. vivax serological test and treatment (PvSeroTAT) [54, 88], as well as the effect of chloroquine and primaquine with vector control [85], and the potential effect of three different types of vaccines that target different stages of the P. vivax life cycle [86] in different geographical settings. The impact of different intervention scenarios, including five annual rounds of MDA on both P. falciparum and P. vivax across the Asia-Pacific, has been studied using a metapopulation model [72]. The only mixed-species agent-based model [87] was used to investigate different treatment scenarios, including current practice, accelerated radical cure, and unified radical cure provided with and without MDA (radical cure was with 14 days of primaquine and a G6PD test while the MDA was with blood-stage treatments only).

The only within-host model that accounted for the effect of MDA on each of the hypnozoites and infections was proposed by Mehra et al. (2022) [34]. This model provided base analytical expressions for the effect of multiple rounds of MDA on hypnozoite dynamics and provided the epidemiological impact of one round of MDA on a single individual. Anwar et al. (2023) recently embedded Mehra et al.’s work [34] and extended the model to study the effect of multiple MDA rounds (up to N rounds) on both within-host and population-level [4]. To the best of our knowledge, no other multiscale model has been developed that explicitly accounts for the effect of multiple rounds of MDA. The interval between MDA rounds when multiple rounds are considered, varies across different models, mostly with either one week [28, 29] or six months [20, 54] with different level of coverage. The only model that provides optimal intervals if multiple MDA rounds were under consideration is by Anwar et al. (2023) [4].