Mathematical models of Plasmodium vivax transmission: A scoping review

Plasmodium vivax is one of the most geographically widespread malaria parasites in the world, primarily found across South-East Asia, Latin America, and parts of Africa. One of the significant characteristics of the P. vivax parasite is its ability to remain dormant in the human liver as hypnozoites and subsequently reactivate after the initial infection (i.e. relapse infections). Mathematical modelling approaches have been widely applied to understand P. vivax dynamics and predict the impact of intervention outcomes. Models that capture P. vivax dynamics differ from those that capture P. falciparum dynamics, as they must account for relapses caused by the activation of hypnozoites. In this article, we provide a scoping review of mathematical models that capture P. vivax transmission dynamics published between January 1988 and May 2023. The primary objective of this work is to provide a comprehensive summary of the mathematical models and techniques used to model P. vivax dynamics. In doing so, we aim to assist researchers working on mathematical epidemiology, disease transmission, and other aspects of P. vivax malaria by highlighting best practices in currently published models and highlighting where further model development is required. We categorise P. vivax models according to whether a deterministic or agent-based approach was used. We provide an overview of the different strategies used to incorporate the parasite’s biology, use of multiple scales (within-host and population-level), superinfection, immunity, and treatment interventions. In most of the published literature, the rationale for different modelling approaches was driven by the research question at hand. Some models focus on the parasites’ complicated biology, while others incorporate simplified assumptions to avoid model complexity. Overall, the existing literature on mathematical models for P. vivax encompasses various aspects of the parasite’s dynamics. We recommend that future research should focus on refining how key aspects of P. vivax dynamics are modelled, including spatial heterogeneity in exposure risk and heterogeneity in susceptibility to infection, the accumulation of hypnozoite variation, the interaction between P. falciparum and P. vivax, acquisition of immunity, and recovery under superinfection.


Introduction
Malaria remains a significant public health problem, with an estimated 247 million cases and 619,000 deaths reported worldwide in 2021 alone [143].Malaria is most prevalent in the World Health Organisation (WHO) African Region, while the South-East Asia Region has the secondhighest estimated malaria burden globally.Plasmodium vivax is currently the most geographically widespread of the malaria parasites, resulting in significant associated global morbidity and mortality [8,10,19,111].P. vivax has been responsible for approximately 45% of malaria cases in the WHO South-East Asia Region since 2000 and is widely prevalent in countries across Asia, Latin America, and the Pacific Islands [19,111,143].P. vivax has often been overlooked and mistakenly considered as "benign" in the past [111,112].More recent research has produced evidence that, in addition to causing severe illness, P. vivax infection can cause long-term health consequences such as anaemia, impaired cognitive development, and chronic kidney disease [12,23,72,131].The economic impact of P. vivax malaria is also significant, as the disease can lead to decreased productivity, increased healthcare costs, and reduced economic growth in endemic areas [39].
Mathematical modelling is an important tool that allows us to understand dynamic systems in various fields ranging from physics and engineering to social sciences and biology [66].Mathematical modelling can provide valuable insight into infectious disease dynamics and plays an important role in informing public health policy and decision-making [14,57].Infectious disease modelling has been widely used to understand the transmission of malaria, particularly Plasmodium falciparum, and the impact of interventions to control and eliminate malaria [82,124].Modelling of P. vivax transmission differs from P. falciparum modelling, due to the need to account for recurrent infections caused by the activation of hypnozoites, a dormant liver stage of the parasite.P. vivax parasites are introduced into the human body through infectious Anopheles mosquito bites.P. vivax parasites then travel to the liver, where they undergo a series of developmental and replication stages [61,65] before the liver-stage parasites are released into the blood, causing blood-stage infections.Individuals experiencing a blood-stage infection may become symptomatic, with symptoms such as fever and fatigue, or be asymptomatic.One of the significant characteristics of P. vivax infection is that, as part of the parasites' life-cycle, they can remain dormant in the liver for weeks or months [58] as hypnozoites that can cause further blood-stage infections (called relapses) upon reactivation.Importantly, between 79 and 96% of P. vivax cases are due to relapses [2,31,54,116].It can be challenging to distinguish a relapse from other types of recurrent malaria, such as a reinfection (i.e.malaria due to a new infectious bite) or a recrudescence (i.e.recurrence of malaria due to incomplete elimination of blood-stage infections, often associated with treatment failure) [46].Relapse dynamics typically follow temperate or tropical phenotypes, relating to the period between primary infection and hypnozoite activation [78].In tropical regions, relapses occur frequently within a few weeks to a few months, whereas in temperate regions, relapses typically occur between six to 12 months after initial infection.This variation in relapse frequency relates to vector dynamics and the transmission potential of P. vivax.In temperate regions, slower-relapsing hypnozoites may allow the parasites to survive colder months when mosquitoes are less prevalent, whereas, in tropical regions, a faster relapsing frequency may allow the parasite to maximise its transmission potential [18,55].As relapses contribute to the majority of blood-stage infections, it is important to capture these relapse dynamics when modelling P. vivax disease transmission.
The methods of incorporating hypnozoites and their associated relapse dynamics vary across the P. vivax modelling literature.Modellers have often adopted the approach of assuming a binary state (presence or absence) for hypnozoites harboured within an individual [28,42,59,60,116].The P. vivax hypnozoite reservoir (i.e. the number of hypnozoites) is known to be non-binary [138,139].Due to this, more recent P. vivax models have attempted to incorporate the complex hypnozoite dynamics and capture the impact of the hypnozoite reservoir on transmission dynamics [9,10,84,87,138].
The methods used to capture P. vivax immunity also vary across the modelling literature.When individuals are first infected with malaria, they naturally develop some level of immunity.This immunity can be defined as the body's state of resistance to the infection, and, with each subsequent infection, this acquired immunity is enhanced [24].Modellers may consider different types of immunity when modelling P. vivax transmission.This includes immunity against new infections, protection 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.immunity that reduces the probability of parasite transmission to mosquitoes) [41,145,146].
One of the primary reasons for modelling infectious disease transmission is to understand the potential impact of treatment strategies on incidence.In terms of P. vivax, a combination therapy, known as radical cure, is needed to target both the acute infection and the dormant hypnozoite reservoir [103,130,135].The two drugs include: (i) a drug that clears parasites from the blood (such as chloroquine or artemisinin-based combination therapy; and (ii) an 8-aminoquinoline drug that clears hypnozoites from the liver (such as primaquine or tafenoquine).Targeting the hypnozoite reservoir is crucial in controlling or eliminating P. vivax, as transmission can be re-established from the reactivation of hypnozoites [138].Incorporating Glucose-6-phosphate dehydrogenase deficiency (G6PD) testing is recommended before administering primaquine or tafenoquine as these drugs can cause life-threatening haemolysis in individuals with G6PD deficiency, an enzymopathy affecting up to 30% of individuals in malaria-endemic regions [114].
Other interventions that have been modelled include vector control, mass drug administration (MDA), mass screening and testing (MSaT), and P. vivax serological testing and treatment (Pv SeroTAT).Vector control measures are recommended by the WHO in order to achieve elimination [147].MDA is an effective intervention for controlling malaria and was advocated by the WHO in the 1950s to control malaria transmission [53].MDA involves treating the entire population, or a well-defined sub-population, in a geographic location regardless of their infection status [53,96], such that both individuals who are infected and non-infected are treated.In a radical cure MDA intervention, individuals are given artemisinin-based combination therapy to clear blood-stage parasites and primaquine (or tafenoquine) to clear hypnozoites.Due to the risks associated with radical cure treatment in G6PD-deficient individuals, mass administration of radical cure is not recommended by the WHO without first screening for G6PD deficiency [52,134,141].Another strategy for reducing and eliminating malaria is MSaT.This involves identifying and treating infected individuals within a specific geographical location by mass testing of all individuals regardless of their symptom status [121].MSaT is effective in reducing malaria transmission in areas with low to moderate malaria prevalence.However, its success depends on the availability of accurate diagnostic tools, effective antimalarial drugs, and strong community participation [70,142].Pv SeroTAT is a method for identifying individuals with recent blood-stage infections who are potential hypnozoite carriers by measuring antibodies and providing treatment with radical cure [98].This method can identify individuals likely harbouring a hypnozoite reservoir, therefore allowing targeted treatment.Mathematical modelling has been used to understand how these different intervention strategies may impact P. vivax transmission [3,60,137].
In this article, we synthesise the findings of a scoping review of existing mathematical models for population-level P. vivax transmission to provide a comprehensive overview of the modelling frameworks and methods used to characterise P. vivax dynamics.In Section 2, we provide the search and inclusion criteria.We discuss the search results in Section 3 as per the categorical structure in Figure 2 before concluding remarks and open problems are presented in Section 4.

Methods
We conducted a literature search on the 21st of May 2023, using the databases PubMed and Google Scholar to capture all relevant studies using the search terms "hypnozoite", "malaria", "vivax", and "mathematical model" with Boolean operators.We screened the titles, abstracts and full text of articles for the following inclusion criteria: • the paper either applied or described a mathematical model of population-level P. vivax transmission dynamics, and; • the mathematical model of P. vivax incorporated hypnozoite dynamics, as this is a distinguishing feature of P. vivax parasites compared to other Plasmodium spp.
We excluded papers that: • were only concerned with the within-host dynamics of P. vivax.Although within-host models of P. vivax dynamics are important for understanding P. vivax transmission, they were not directly relevant to the aim of our study (i.e. to identify and compare mathematical models of population-level P. vivax transmission).Papers that modelled dynamics at both the withinhost and population level (i.e.multi-scale models) were included.
• only used or described mathematical models of Plasmodium species other than P. vivax (e.g. a mathematical model of P. falciparum infectious disease dynamics).Models that accounted for both P. vivax and another Plasmodium species were included.
• were currently only available as a preprint.
Search terms were conducted in English only, and only literature published in English were considered.No limitations regarding study location, publication status (e.g.accepted, but no preprints), publication type, or publication year were included.To enhance the probability of finding all relevant literature, we screened all references within the articles that met our inclusion and exclusion criteria.Articles were then downloaded to identify key components, which are discussed in Section 3.
We categorised models depending on whether they used a stochastic or deterministic approach, and whether they were compartmental or agent-based.Deterministic models have no random variation and typically utilise a compartmental structure within a population to form differential equations to track the rate of flows between compartments.Stochastic models incorporate random variation and are useful for questions and scenarios where small population numbers or extinction are involved.In terms of P. vivax infectious disease modelling, agent-based models explicitly model P. vivax transmission dynamics at an individual-level, for example, modelling the interaction between humans and vectors and associating respective state variables and parameters to each individual and vector.In our review, we found that almost all stochastic models were also agent-based, so even though these features are not mutually exclusive, we categorised models as (i) deterministic compartmental models or (ii) stochastic agent-based models.

Search results
The initial search yielded 2289 articles, which was reduced to 1005 unique articles after removing duplicates between the two databases.After screening at the title level, a further 901 studies were excluded as they did not fulfil the selection criteria in Section 2. After screening the abstracts, a further 63 studies were excluded due to either (i) no underlying mathematical model being described or (ii) the model was for P. falciparum parasites only.Five additional studies were included from the selected studies' references that were not initially identified.A total of 47 studies were finally selected for review (see Figure 1 for a summary of the selection process).

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 [81], 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 byto the best of our knowledge -Zoysa et al. (1988) [146] in a Ross-Macdonald style modelling approach.Following this, many models have now been developed.
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 [22,132].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) [132].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 [45,102].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 [54,94,95,98,133,136,137]. 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 [22,89,144].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 impossible, as there will always be some non-zero number of infectious mosquitoes remaining that can trigger an infection in humans again [133].
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 [45,102].The most recent agent-based models modelling P. vivax dynamics [54,94,95,98,136] have evolved from a model introduced by White et al. (2018) [137].The White et al. (2018) [137] model has been adapted to capture disease epidemiology in particular geographical settings [95], and to study the impact of different interventions (drugs or vaccination) [54,94,136].
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 [30,37,75].For example, despite being an agent-based model, parameterisation is done using an ordinary differential equation system that describes the process in several models [54,95,98,137].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 [118], although with improving computer technology, this has become less of a concern [38].

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 [43].
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 Figure 2).Few models capture both within-host and population-level impacts (boxes with a strong blue border in Figure 2) [9,54,95,98,[136][137][138].The very first multiscale model for P. vivax transmission was developed by White et al. (2018) [137], and modelled the within-host hypnozoite dynamics using an agentbased model that considered heterogeneity in exposure to mosquito bites.This built on White et al. (2014) [138], 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 [137].The White et al. (2014) [138] 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 [85].The other within-host models that were adapted from White et al. (2018) [137] applied the same assumption regarding batch hypnozoite behaviour [54,95,98,136].Mehra et al. (2020) [85] 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 Figure 2) in White et al. (2014) [138] in analytical form.Later work by Mehra and colleagues embedded the activation-clearance model governing a single hypnozoite in an epidemiological framework [87].This framework accounts for successive mosquito bites, where each bite can simultaneously establish multiple hypnozoites [86,87], and explores the epidemiological consequence of radical cure treatment on a single individual.Anwar et al. (2022) [9] have since developed a multiscale model motivated by White et al. (2014) [138] by embedding the framework of Mehra et al. (2022) [87] 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 withinhost 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) [138] into three compartments and relaxes the artificial truncation needed in White et al. (2014) [138] for numerical simulation.Mehra et al. (2022) [84] 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) [87].In the infinite population size limit, Mehra et al. (2022) [84] 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 [87].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 [84].

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 [145].Since most P. vivax blood-stage infections are due to the reactivation of hypnozoites rather than new primary infections, it is crucial that mathematical models incorporate the size of the hypnozoite reservoir [13,21,32,34,79].Zoysa et al. (1991) [145] 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) [42].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 [139].Despite having the relapse characteristic that makes P. vivax parasites unique, Aldila et al. (2021) [5] did not incorporate relapses in their P. vivax transmission model.In this model, individuals did not harbour hypnozoites when infected with P. vivax and hence did not experience relapse after recovery from blood-stage infection.
Most P. vivax transmission models consider the hypnozoite reservoir as a single compartment, rather than explicitly accounting for a variable number of hypnozoites in the reservoir [1, 3, 4, 11, 27, 35, 44-46, 56, 59, 60, 63, 64, 68, 71, 92, 97, 100, 102, 104-107, 113, 116, 117, 120, 128, 139].Only a handful of models account for the variability in hypnozoite inoculation across mosquito bites (boxes with a bright pink border in Figure 2) [9,10,87,138].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) [138] 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) [137], 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 [54,95,98,136].Furthermore, instead of explicitly modelling hypnozoites independently, they impose the batch hypnozoite model.This assumption means that hypnozoites from a mosquito bite act as a batch, where they all reactivate simultaneously, causing relapse or dying at the same time.This reduces one batch of hypnozoites to a single set of dynamics, which is still truncated at a maximum of k batches.[9] is that the population-level component is considerably simpler than the 2 L max + 1 ordinary differential equations of White et al. (2014) [138].The transmission models proposed by Mehra et al. (2022) [84] likewise account for variation in hypnozoite batch sizes, with Mehra et al. (2022) [84] additionally accommodating long-latency hypnozoite dynamics.The models of Mehra et al. [84] are formulated as systems of integrodifferential equations, informed by the within-host framework of Mehra et al. (2022) [87].The analyses of Anwar and Mehra et al. [9,84] 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.

Superinfection
Superinfection with malaria is a common phenomenon 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 [122].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 [110,122].
When modelling P. vivax dynamics, it is important to consider the impact of superinfection on recovery and transmission, as the abundance of mosquitoes and the contribution of hypnozoite activation can frequently trigger superinfection.Superinfection can potentially delay recovery from infection [40,123].Most of the literature that incorporates superinfection in P. vivax transmission models (boxes with a brown border in Figure 2) [42,54,59,60,95,100,120,[136][137][138] does so via the recovery rate [42,59,60,138].The superinfection phenomenon was first introduced into malaria models by Macdonald (1950) [80], 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) [40] 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) [42,59,60].Since White et al. (2014) [138] accounts for the variation of hypnozoites, they modified the recovery rate proposed by Dietz et al. (1974) [40] to account for the additional burden of hypnozoites; however, Mehra et al. (2022) [84] 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 [10,84,93] 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 [84] Superinfection was incorporated in later work, where it was assumed that different batches of hypnozoites originated from different mosquito bites [54,98,136,137].Silal et al. (2019) [120] 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 [5], 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 [17].Accordingly, it may not be reasonable to consider this to be a valid model of superinfections.

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 [120,133].For example, in 2019 in Cambodia, co-infection with both P. vivax and P. falciparum accounted for about 17% of malaria cases [29].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 [76,120,125,140].The high risk of P. vivax parasitaemia after P. falciparum infection is possibly related to reactivation of hypnozoites [33,51,129].Hypnozoites may be activated when P. falciparum parasites have been introduced into the body [119] or when the individual is exposed to Anopheles specific proteins [55].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 [7,49,50].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 [36,91].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 [73].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 [88].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 [127].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 [3,5,102,105,106,120].While both P. vivax and P. falciparum species are included in a single model by Aldila et al. (2021) [5], 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 [3,5,102,105,106].
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 [120,125,133].To the best of our knowledge, the first model that accounts for the interaction between both species was developed by Silal et al. (2019) [120].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) [133].This model had reduced complexity compared with Silal et al.'s (2019) co-infection model, but used many of the same parameter values [120] (co-infection models shown with a vivid orange bordered box in Figure 2).

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 [26,83].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 [15,48,77,90].In co-endemic areas, clinical immunity to P. vivax is more rapidly acquired than that due to P. falciparum [90].
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 [145,146].This type of immunity is assumed to be boosted by infection [15].Acquiring some partial immunity (i.e.some degree of protection against malaria) following infection that wanes over time, is most common among published models [3,27,44,45,63,64,[104][105][106][107]128].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 [117,138].The assumption regarding permanent immunity against malaria is not considered valid, as immunity often wanes rapidly when immune adults leave malaria-endemic regions [74].Despite this, one model assumed that recovered individuals become permanently immune to P. vivax [5].Another study assumed that only a fixed proportion of individuals are immune against P. vivax rather than explicitly incorporating immunity into the model [102].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 [74].Repeated exposure to infectious bites will still likely result in infection, though these individuals may be protected against severe malaria or death [74].Silal et al. (2019) [120] applied the opposite assumption and assumed that repeated exposure to infectious bites would likely result in severe infection [120].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 [25,109,137].
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 Figure 2).The only model that explicitly accounts for the acquisition of immunity that increases with new bites was developed by White et al. (2018) [137].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 [137].They also assumed that children acquired immunity through their birth parent's immunity, which then decayed exponentially from birth.Models that were adapted from White et al. (2018) [137] also allow for the acquisition of immunity [54,95,98,136].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 [62,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.

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 [27,42,44,54,56,69,95,100,102,117,120,128,138], but a range of other interventions can contribute to malaria control.For example, ten studies (21%) evaluated the effect of MDA on disease transmission, despite few national control programs considering MDA for P. vivax control [3,10,59,60,87,94,98,116,133,137] (boxes with a dark green border in Figure 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 [47,99,116].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 [116].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% [137].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) [94,98], as well as the effect of chloroquine and primaquine with vector control [95], and the potential effect of three different types of vaccines that target different stages of the P. vivax life cycle [136] in different geographical settings.The mixed-species agent-based model [133] 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) [87].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 [87] and extended the model to study the effect of multiple MDA rounds (up to N rounds) on both within-host and population-level [10].To the best of our knowledge, no other multiscale model has been developed that explicitly accounts for the effect of multiple rounds of MDA.Anwar et al. (2023) [10] also provided optimal intervals if multiple MDA rounds were under consideration.

Open questions and conclusion
Mathematical modelling is a powerful tool for understanding, analyzing, and predicting complex real-world phenomena, as well as simulating different scenarios, testing hypotheses, and making informed decisions based on the results.Mathematical models have proven useful to characterise P. vivax transmission in different parts of the world and provide insights into the effect of different strategies to achieve elimination, including treatment, vaccination, and vector control.In this work, we provided a review of the existing mathematical models that capture P. vivax disease progression and transmission.P. vivax transmission dynamics are particularly challenging to model given the difficulties discerning relapses from reinfections and recrudescences.The choice of transmission model framework comes down to the research question at hand.While mathematical models can provide key insights without the expense of large trials or epidemiological studies, it is important to recognize that mathematical models are not perfect representations of reality and are always subject to limitations, uncertainties, and assumptions.Therefore, using mathematical models in conjunction with empirical data, expert knowledge, and critical thinking is essential to obtain meaningful and reliable results.
Across the different approaches of mathematical modelling of P. vivax, there were varying assumptions regarding parasite dynamics and acquisition of immunity.Some models were motivated to capture realistic biological aspects of the parasite [9,138,146], or epidemiological and public health aspects [3,10,27,35,42,45,54,59,60,87,95,98,102,116,117,120,128,133,136,137,139,145], whereas some models were motivated to construct a novel or extended mathematical model of P. vivax dynamics, i.e., focusing on the mathematical aspects of P. vivax dynamics [1, 5,44,46,56,63,64,68,71,92,97,100,[104][105][106][107]113].As the dynamics of these type of models are well established, we argue that more importance should be placed on using these models to address the current hurdles and setbacks in achieving P. vivax elimination.For example, the effect of new drugs, emerging drug resistance, and the potential effect of vaccination (when it becomes available).Modelling different scenarios with the available tools under the current recommendations is crucial to inform decision-making regarding malaria elimination.Furthermore, given that some of the biological aspects of P. vivax are well understood, we argue that researchers should shift their focus to modelling these important aspects.
The spatial distribution of P. vivax transmission is heterogeneous, and the number of hypnozoites that an individual harbours might vary significantly; this contributes directly to the risk of hypnozoite reactivation and P. vivax relapse [108,126].This heterogeneity can be partially captured by modelling individuals' movement using metapopulations and including parasite movement between different sub-populations [95,137].However, none of the current models explicitly consider this spatial heterogeneity.Given the high degree of heterogeneity of P. vivax risk in almost all populations, future model development should address this.
As more than 80% of P. vivax infections may be due to relapse, and multiple hypnozoites can be established from each infectious bite, modelling the dynamics of hypnozoite variation and activation is crucial [9,87,116].Another important aspect that requires more detailed attention is the interaction between multiple species of Plasmodium, particularly in areas where P. falciparum and P. vivax are co-endemic (Asia, the Horn of Africa, and the Americas).Studies show that there is a high risk of P. vivax parasitaemia after P. falciparum infection that is possibly related to reactivation of hypnozoites [33,51,129].This is in line with the hypothesis that P. falciparum infection might trigger underlying P. vivax infection [76,120,125,140].Hence, we argue that this hypothetical triggering phenomenon should be investigated when modelling P. falciparum and P. vivax interactions.Future P. vivax modelling efforts should also account for superinfections.Where mosquito abundance is high, transmission intensity is also likely to be high if malaria parasites are present [16,67,110,115].In these scenarios, infected individuals are likely to experience multiple episodes of infection at once (i.e.superinfection).Superinfection can significantly delay recovery time, leaving ample opportunity for onward transmission from the infected individual to susceptible mosquitoes.P. vivax models should hence account for the transmission dynamics associated with superinfection.Immunity against P. vivax strongly correlates to past exposure; therefore, focus should also be placed on modelling the acquisition (and waning) of immunity related to superinfection, as multiple concurrent exposures may boost immunity more than singular exposures [84,137].Furthermore, as parasites from relapse are either genetically identical or related to a previous primary infection, they are more efficiently targeted by naturally acquired immune responses previously developed from the primary infection than further, genetically unrelated primary infections.As a consequence, relapses are less likely to be associated with (severe) clinical symptoms.[62,140].This interplay between immunity and relapse has not been fully addressed in any models developed to date.Given these important biological aspects, we suggest that future modelling should focus on developing the above-mentioned key areas: (i) spatial heterogeneity in exposure risk, (ii) accumulation of hypnozoites variation, (iii) P. falciparum and P. vivax interactions, (iv) acquisition of immunity, and (v) recovery under superinfection.Different modelling communities have recently started focusing on these areas recently, for example, modelling hypnozoite dynamics [85,87], multispecies interactions (P.falciparum and P. vivax ) [120,133], bite exposure immunity [137] and superinfection [10,84,87].No model currently includes all of the above factors that play a role in P. vivax transmission due to the complexity the resulting model would have, and not all of the factors may need to be modelled to answer the research questions at hand.Therefore, when developing models to explore P. vivax disease progression with a focus on answering specific research questions, mathematical epidemiologists and modellers should consider relevant aspects within the context of existing recommendations.
To address the outstanding research questions identified here, a suitably skilled interdisciplinary team is required.We hope that this review can contribute to developing the common language needed for communication between different scientists by highlighting the progress of P. vivax transmission models to date.

Figure 1 :
Figure 1: Summary of the article selection process, illustrating papers included and excluded at each stage of the review process.

Figure 2 :
Figure 2: A summary of the 47 P. vivax transmission models currently available in the literature (published as of May 21, 2023).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.
The multiscale model developed by Anwar et al. (2022) [9] accounted for the variation of hypnozoites dynamics across bites.Unlike the White et al. (2014) [138] 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) [87] as a system of integrodifferential equations.This relaxes the artificial truncation for the maximum number of hypnozoites used within the White et al. (2014) [138] model.Under a constant force of reinfection, Anwar et al. (2022) [9] analytically proved that the multiscale model [9] exhibits an identical steady-state hypnozoite distribution as the infinite ordinary differential equation model structure in White et al. (2014) [138].The advantage of the multiscale model by Anwar et al. (2022) . The within-host model of Mehra et al. (2022) [87] 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) [9], which built on work of Mehra et al. (2022)[87], did not incorporate superinfection.A correction to account for superinfection, based on the recovery rate formulated by Nåsell et al. (2013) [93], was proposed in Mehra et al. (2023) [84] and incorporated in later work by Anwar et al. (2023) [10].

FundingL.
Smith is supported by the National Health and Medical Research Council (NHMRC) (GNT2016726) and the Department of Foreign Affairs and Trade Australia through the project Strengthening Preparedness in the Asia-Pacific Region through Knowledge (SPARK). A. Devine's research is supported through the NHMRC (2019152).E. Conway, and I. Mueller's research are supported by the NHMRC (GNT2016726) and the Department of Foreign Affairs and Trade Australia through the project Strengthening Preparedness in the Asia-Pacific Region through Knowledge (SPARK).J.M. McCaw's research is supported by the Australian Research Council (DP210101920) and the NHMRC Australian Centre of Research Excellence in Malaria Elimination (ACREME) (APP1134989).J.A. Flegg's research is supported by the Australian Research Council (DP200100747, FT210100034) and the NHMRC (APP2019093).The contents of the published material are solely the responsibility of the individual authors and do not reflect the views of NHMRC.