A hybrid transmission model for Plasmodium vivax accounting for superinfection, immunity and the hypnozoite reservoir

Malaria is a vector-borne disease that exacts a grave toll in the Global South. The epidemiology of Plasmodium vivax, the most geographically expansive agent of human malaria, is characterised by the accrual of a reservoir of dormant parasites known as hypnozoites. Relapses, arising from hypnozoite activation events, comprise the majority of the blood-stage infection burden, with implications for the acquisition of immunity and the distribution of superinfection. Here, we construct a novel model for the transmission of P. vivax that concurrently accounts for the accrual of the hypnozoite reservoir, (blood-stage) superinfection and the acquisition of immunity. We begin by using an infinite-server queueing network model to characterise the within-host dynamics as a function of mosquito-to-human transmission intensity, extending our previous model to capture a discretised immunity level. To model transmission-blocking and antidisease immunity, we allow for geometric decay in the respective probabilities of successful human-to-mosquito transmission and symptomatic blood-stage infection as a function of this immunity level. Under a hybrid approximation—whereby probabilistic within-host distributions are cast as expected population-level proportions—we couple host and vector dynamics to recover a deterministic compartmental model in line with Ross-Macdonald theory. We then perform a steady-state analysis for this compartmental model, informed by the (analytic) distributions derived at the within-host level. To characterise transient dynamics, we derive a reduced system of integrodifferential equations, likewise informed by our within-host queueing network, allowing us to recover population-level distributions for various quantities of epidemiological interest. In capturing the interplay between hypnozoite accrual, superinfection and acquired immunity—and providing, to the best of our knowledge, the most complete population-level distributions for a range of epidemiological values—our model provides insights into important, but poorly understood, epidemiological features of P. vivax.

• Antiparasite immunity modulates the clearance of blood-stage infection [26], and is typically modelled through accelerated parasite clearance rates and/or reduced parasite densities [13,34].
Transmission-blocking immunity, which modulates the infectivity of sexual blood-stage parasites (gametocytes) to mosquitoes, is also of note, with the mitigation of mosquito-stage development curtailing onward transmission [6,18,25,41]. There is evidence to suggest that the gametocyte circulation is not hampered by clinical and antiparasite immunity [38].
Other forms of immunity are thought to be of comparatively limited consequence in natural transmission settings. Pre-erythrocytic immunity targets sporozoites inoculated through mosquito bites, prior to further liver-stage development. Due to the potential for each sporozoite to develop into a hypnozoite, pre-erythrocytic immune protection has been hypothesised to substantially mitigate the relapse burden [18,29]; exposure to sporozoites in natural transmission settings, however, is believed to be insufficient to induce strong pre-erythrocytic immune protection [18]. Likewise, immune responses targeted towards liver-stage parasites, particularly hypnozoites, are poorly understood [12], but are generally considered to be relatively minor.
The joint dynamics of immunity and the hypnozoite reservoir are of epidemiological interest.
The dichotomisation of both hypnozoite carriage and immune status [7,8,15,20] yields, in some senses, the simplest approach for characterising population-level transmission dynamics. Various models of immunity have been proposed under these dichotomised frameworks, ranging from imperviousness to reinfection (until immunity is lost) [7], to an elevated rate of recovery (antiparasite) [8]; reduced infectiousness to mosquitoes (transmission-blocking) [20]; and necessarily asymptomatic blood-stage infection (clinical) [15]. A slightly extended model of transmissionblocking immunity, superinfection and hypnozoite accrual has been proposed by De Zoysa et al. [5] -with the limitation that each individual can harbour up to two broods of hypnozoites and two overlapping relapses, and a discrete immunity level {0, 1, 2}.
A more comprehensive characterisation of immunity and the hypnozoite reservoir has recently been performed by by White et al. [34]. Under a hypnozoite 'batch' model -whereby hypnozoites are stratified into 'batches', each characterised by a constant rate of relapse over the span of an exponentially-distributed lifetime, with an imposed upper bound K on concurrent batch carriage -White et al. [34] account for the acquisition of both antidisease immunity (which re-duces the probability of symptomatic infection) and antiparasite immunity (which results in an elevated rate of parasite clearance and a reduced probability of detection via light miscroscopy).
In addition to being restricted to short-latency strains, the framework of White et al. [34] ignores size variation in parasite inocula, as we noted in Mehra et al. [49]. Further, White et al. [34] do not explicitly account for superinfection.
Here, we seek to characterise the interplay between hypnozoite accrual, superinfection and acquired immunity, for both short-and long-latency strains. In Mehra [48], we have recently proposed a transmission model for P. vivax that explicitly accounts for superinfection and (short-latency) hypnozoite accrual; our model can be viewed as an extension of pre-existing hypnozoite 'density' models [24,46] to accommodate a rigorous characterisation of superinfection. Under the deterministic model derived in Mehra [48], we can recover population-level distributions for various quantities of epidemiological interest without encountering the computational overheads that have curtailed previous efforts to model explicit hypnozoite densities [34]. The conceptual underpinning of the model detailed in Mehra [48] is the within-host framework we introduced in Mehra et al. [49], which captures hypnozoite and superinfection dynamics as a function of mosquito-to-human transmission intensity, or the force of reinfection (FORI), in a general transmission setting. In the present paper, we adopt an analogous mathematical construction to Mehra [48], extending our previous work to allow for long-latency (temperate) strains and the acquisition of transmission-blocking and antidisease immunity. This paper is structured as follows. Section 2 focuses on the characterisation of within-host dynamics as a function of the FORI. We begin by extending the open network of queues introduced in Mehra et al. [49], which describes the joint dynamics of superinfection and the hypnozoite reservoir, to include a discretised immunity level (Section 2.1); as observed in Mehra et al. [49], this immunity level is governed by a shot noise process, akin to a previous model of antibody dynamics we introduced in Mehra et al. [44]. Rather than solving for the state probabilities based on the Kolmogorov forward differential equations for the queueing network (Section 2.2), we derive a time-dependent PGF as in our previous work [44,49] (Section 2.3). Specific models for antidisease and transmission-blocking immunity are proposed in Sections 2.4 and 2.5 respectively. Section 3 concerns the construction of a hybrid transmission model [19,40], predicated on the coupling of expected host and vector dynamics. To recover the expected dynamics of the vector population, we consider the Kolmogorov forward differential equations for an underlying birth-death process (Section 3.1); while observing that the within-host probability mass function (PMF) can be regarded as the expected population-level frequency distribution of hypnozoite, superinfection and immunity states [40]. We then derive a infinite compartmental model to couple host and vector dynamics (Section 3.2). Steady state analysis -including the identification of a bifurcation parameter governing the existence of endemic equilibria (Section 3.2.1) and a sensitivity analysis of endemic equilibrium solutions (Section 3.2.2) --is performed using the within-host distributions derived in Mehra et al. [49]. To characterise transient population-level dynamics, we adopt the approach detailed in Mehra [48] to derive a reduced system of IDEscomprising an integral equation for the immune-modulated probability of human-to-mosquito transmission (per bloodmeal), and a set of ordinary differential equations (ODEs) governing the number of (un)infected and latent mosquitoes over time (Section 3.3). As a function of the FORI derived under the reduced system of IDEs, we recover population-level distributions for various quantities of epidemiological interest -including the size of the (non)-latent hypnozoite reservoir; superinfection; the prevalence of clinical infection and the relative contribution of relapses to the infection burden -using the distributions derived in Mehra et al. [49]. Concluding remarks are made in Section 4.
2 Within-host human dynamics: hypnozoite accrual, superinfection and immunity We have previously derived the functional dependence between the FORI, and the joint dynamics of blood-stage infection and the hypnozoite reservoir by constructing an open network of infinite server queues [49]. Here, we extend the model detailed in Mehra et al. [49] to allow for the acquisition of immunity. Following the approach detailed in Appendix C.3 of Mehra et al. [49], we assume that the within-host acquisition of immunity is described by a generalised shot noise process such that: • the clearance of each primary infection/relapse elicits a boost of unit magnitude; • the lifetime of each boost is exponentially-distributed with mean 1/w; and • the overall immunity level is given by the cumulative sum of boosts over the entirety of an individual's infection history.
As noted in Mehra et al. [49], this discretised model of immunity can be considered a variation of the antibody model proposed in Mehra et al. [44], in which the clearance of each primary infection/relapse elicits a boost of random magnitude that is then subject to exponential decay at a fixed (deterministic) rate.
In Section 2.1 below, we propose an open network of infinite server queues to capture the within-host dynamics of superinfection, hypnozoite accrual and immune acquisition [49]. The Kolmogorov forward differential equations governing the time evolution of the joint PMF for the network are stated in Section 2.2. Instead of directly solving the the Kolmogorov forward differential equations (which comprise an infinite-dimensional set of ODEs), we derive a joint PGF for the state of the network following a similar approach to Mehra et al. [49] (Section 2.3).
Specific models for transmission-blocking and antidisease immunity are detailed in Sections 2.4 and 2.5 respectively. To elucidate the dynamics captured by our within-host model, we discuss an illustrative sample path in Section 2.6.

An open network of infinite server queues
To capture the within-host acquisition of immunity, we extend the open network of queues detailed in Mehra et al. [49] to include an additional node I, such that the occupancy of queue • i ∈ {1, . . . , k} to represent hypnozoites that are present in latency compartment i (that is, part of the hidden liver-stage reservoir but unable to activate) • N L to represent non-latent hypnozoites (that is, part of the hidden liver-stage reservoir and able to activate) • A to represent ongoing relapses from activated hypnozoites • D to represent hypnozoites that have died prior to activation • P to represent ongoing primary infections • I to represent cleared blood-stage infections (primary infections or relapses) that have given rise to an immunity increment of unit magnitude • C to represent the loss of immune memory.
As such, the state space for each hypnozoite is S h = {1, . . . , k, N L, A, D, I, C}, while the state space for each primary infection is S p = {P, I, C}.
Arrivals into the network, which represent infective bites, are modelled to follow a non-homogeneous Poisson process with rate λ(t). The consequences of each infective bite are two-fold: • a primary infection is immediately triggered, that is, a single "customer" enters queue P ; • a geometrically-distributed batch of hypnozoites (with mean size ν) is established in the liver, entering latency compartment 1 in the case of long-latency strains (k > 0); and the non-latent compartment N L in the case of short-latency strains (k = 0).
In contrast, non-latent hypnozoites (state N L) undergo death at rate µ, and activation at rate α. As such, service times in queue N L are modelled to be exponentially-distributed with mean 1/(α + µ), with departures routed either into queue A (in which case hypnozoite activation has triggered a relapse) with probability α/(α + µ), or queue D (in which case the hypnozoite has died prior to activation) with probability µ/(α + µ).
The clearance of each blood-stage infection is assumed to be independent, and modelled to occur at some constant rate γ, amounting to exponentially-distributed service times (with mean duration 1/γ) in both queues A and P (representing relapses and primary infections respectively).
To capture the boosting of immunity with exposure, we assume that queue I receives all de-partures from queues A and P (corresponding to cleared blood-stage infections); that is, an immune boost of unit magnitude is acquired upon the clearance of each primary infection or relapse. To capture the waning of immune memory with time, we assume that each immune boost is retained for an exponentially-distributed period of time with mean 1/w -coinciding precisely with the service time in queue I. The number of busy servers in queue I therefore acts as a measure of within-host immunity. All departures from queue I are routed to queue C, where they remain indefinitely.
In a natural generalisation of this queueing network, the the stratification of blood-stage infection and immunity into different compartments could allow us to capture additional stages of the parasite lifecycle and further biological realism.

Modelling correlates of immunity
We can formulate correlates of immune protection as time-dependent functionals of the state of the open network, with the the host immunity level N I (t) mapped to the degree of immune protection at time t. To preserve the independence structure of the queueing network, however, these functionals cannot have any direct feedback into the within-host model. This limits the forms of immunity that are analytically tractable under our model.
A key assumption of the within-host model is that the arrival process is independent of the state of the network. As such, we cannot capture pre-erythrocytic immunity, which modulates the probability of successful mosquito-to-human transmission, whereby the arrival process, comprising mosquito bites, would depend on the number of busy servers N I (t) in node I. The assumption of independence between service rates within each node and the state of queueing network is equally important. We are thus unable to account for the potential modulation of (blood-stage) parasite clearance rates -or equivalently, the service rate for nodes A and Pas a function of the host immunity level N I (t), that is, the state of node I. We could, however, introduce deterministic time variation in the rate of clearance of blood-stage infection γ to model age-related physiological factors.
Immune correlates that are amenable under our within-host framework include: • The probability of exhibiting clinical symptoms or high-density parasitemia, as a manifestation of antidisease immunity (Section 2.4).
• The probability of human-to-mosquito transmission when an uninfected mosquito takes a bloodmeal from a blood-stage infected human host, as a measure of transmission-blocking immunity (Section 2.5).
There is evidence to suggest that these forms of immunity are acquired on different time scales.
Transmission-blocking immune memory, for instance, is believed to be relatively short-lived, with boosting driven largely by successive blood-stage infections in intervals of < 4 months [6]; antidisease immunity, in contrast, is believed to be more robust and longer-lived [18]. By augmenting the rate of decay of the probability of symptomatic blood-stage infection (antidisease) as a function of N I (t), relative to the probability of human-to-mosquito transmission (transmission-blocking), we allow for strong antidisease protection to develop more rapidly than transmission-blocking protection, and be maintained at lower transmission intensities. Before discussing these immune correlates, however, we derive an analytic expression for the distribution of the state of the queueing network at time t [49].

Kolmogorov forward differential equations
Denote by N s (t) the number of hypnozoites/infections in each state s ∈ S := S h ∪ S p at time t and set Then by the Kolmogorov forward differential equations, the time evolution of the state probabilities H i 1 ,...,i k ,i N L ,j,k (t) is governed by the countable system of ODEs reinfection (geometric batch of hypnozoites + primary infection triggered) death of a hypnozoite in the liver (latent or non-latent) prior to activation progression of a latent hypnozoite to the next latency compartment activation of a non-latent hypnozoite, triggering a relapse + γ − jH i 1 ,...,i k ,i N L ,j,k (t) + (j + 1)H i 1 ,...,i k ,i N L ,j+1,k−1 (t) clearance of blood-stage infection + gain of immunity increment + w − kH i 1 ,...,i k ,i N L ,j,k (t) + (k + 1)H i 1 ,...,i k ,i N L ,j,k+1 (t) waning of immune memory . (1) Consider a human population of fixed size P H , with each individual taken to be immune-and infection-naive at time zero. In the absence of demography (that is, birth/death), we can reinterpret the within-host PMF H i 1 ,...,i k ,i N L ,j,k (t) as the expected proportion of humans with i m hypnozoites in state m ∈ {1, . . . , k, N L}; a blood-stage infection comprising j parasite broods and immunity level k. Equation (1) can therefore be thought to govern the expected proportion of humans in each hypnozoite/infection/immunity state [40]. In Section 3.2, we will draw on Equation (1) to construct a hybrid transmission model, comprising a countably infinite system of ODEs. Our aim in the present section, however, is to characterise the within-host PMF.
While the infinite-dimensional system of ODEs given by Equation (1)

Joint PGF for the state of the queue
Rather than solving Equation (1) to yield the probability mass function (PMF) for the state of the queue directly, we derive a joint PGF for the state of the network from first principles, using an argument which is an extension of that in Mehra et al. [49]: The PGF can be viewed as a transformation of the PMF. Since there is a one-to-one correspondence between PGFs and PMFs, the derivation of a joint PGF is sufficient to uniquely characterise the distribution of the queueing network. Under the assumption of geometricallydistributed batch arrivals, we can invert the marginal PGF to recover PMFs for quantities of epidemiological interest, as discussed in Mehra et al. [49].
Treating the dynamics of each hypnozoite/infection to be independent [3], we begin by characterising the probability mass function for a single hypnozoite/primary infection that enters the network at time zero. Here, we extend the activation-clearance model proposed by White et al.
[24] -and discussed in detail in Mehra et al. [42] -to allow for the clearance of blood-stage infection (which was also examined in Mehra et al. [49]) and the gain/loss of immunity (as introduced in the present manuscript). To characterise hypnozoite dynamics, we consider the flow of an arrival into either queue N L (for short-latency strains) or queue 1 (for long-latency strains) through the queueing network shown in Figure 1. Similarly, the dynamics of each pri-mary infection are described by the flow of an arrival into queue P through the network.
Denote by p h,s (t) the probability that a hypnozoite established at time zero is in state s ∈ S h at time t. By the Kolmogorov forward differential equations, it follows that dp h, with the initial condition Following similar reasoning to Mehra et al. [42,49], we can solve the system of ODEs given by Equations (2) to (6) analytically; solutions are given in Equations (32) to (35) in Appendix A.
Since we are not interested in the the distribution of dead hypnozoites or cleared infections over time, we do not provide solutions to the ODEs (7) and (8).
Likewise, we can characterise the probabilistic time course for each primary infection. Denote by p p,s (t) the probability that a primary infection triggered at time zero is in state s ∈ S p at time t. We can solve the Kolmogorov forward equations dp p,P dt = −γp p,P (t) dp p,I dt = −wp p,I (t) + γp p,P (t) dp p,C dt = wp p,I (t) with initial condition p p (0) = (p p,P (0), p p,I (0), p p,C (0)) = (1, 0, 0) to yield p p,P = e −γt p p,I = γ γ − w e −wt − e −γt p p,C (t) = 1 − p p,P (t) − p p,I (t). (10) Embedding these state probabilities in an epidemiological framework, as elucidated in Mehra et al. [49], the joint PGF for is given by In Mehra et al. [49], we recovered analytic expressions for the distributions of several biologicallyrelevant quantities -encompassing the size of the (non)-latent hypnozoite reservoir; the number of parasite broods co-circulating in the bloodstream; the relative contribution of relapses to the infection burden and the cumulative number of recurrences (that is, primary infections and relapses) experienced over time -using the joint PGF given by Equation (11). Formulae relevant to the present manuscript are recapitulated in Appendix B.

Antidisease immunity
While the number of broods co-circulating in the bloodstream at time t is given by the total number of busy servers in nodes A and P , that is, N A (t) + N P (t), a large proportion of P. vivax infections in endemic settings are asymptomatic, with implications for treatment and elimination strategies [31,33,47]. The relative burden of (a)symptomatic blood-stage infection, which is a function of antidisease immunity, is therefore of epidemiological interest.
Conditional on the presence of blood-stage infection, we assume that the probability of an individual exhibiting clinical symptoms decreases by a factor of p c for each increment of immunity they harbour. As such, the probability of an individual with state N(t) exhibiting clinical symptoms is given by Accounting for stochasticity in within-host hypnozoite and infection dynamics, the probability of an individual exhibiting clinical symptoms at time t can be written using the law of total expectation. Setting z i = 1 for i ∈ S \ I in Equation (11) to recover the marginal PGF for N I (t), we can write By Xekalaki [4], we can recover the unnormalised PGF for N I (t), conditional on the absence of blood-stage infection (that is, N A (t) + N P (t) = 0), by firstly setting z i = 1 for i ∈ S \ {A, P, I} in Equation (11) (to recover the joint PGF for (N A (t), N P (t), N I (t)), and then setting z A = z P = 0 (to exclusively consider the case N A (t) + N P (t) = 0), yielding the expression Substituting Equations (13) and (14) into Equation (12), we recover an expression for p clin (t) as a function of the FORI λ(τ ) in the interval τ ∈ [0, t): where we have used the joint PGF given by Equation (11). In a similar vein, we can introduce analogous models linking the probabilities of (sub)microscopic parasitemia and detectability (through light microscopy versus rapid diagnostic tests versus qPCR assays) to the immunity level N I (t).

Transmission-blocking immunity
Here, we propose a model for transmission-blocking immunity by introducing a functional dependence between the probability of successful human-to-mosquito transmission and the immune status of a blood-stage infected individual.
For an immune-naive, blood-stage infected individual, we set the probability of successful human-to-mosquito transmission to be p 0 . We further assume that the probability of successful humanto-mosquito transmission is reduced by a factor of p tb ∈ [0, 1] for each increment of immunity.
Suppose a mosquito takes a bloodmeal from a human with state N(t) at time t. Conditional on the state of a human host N(t), we thus define the probability of successful human-to-mosquito transmission p h→m (t) to be Under our stochastic epidemiological framework, following similar reasoning to Section 2.4, we can recover the probability of successful human-to-mosquito transmission when a mosquito bites an individual at time t as a function of the FORI λ(τ ) in the interval τ ∈ [0, t) Note that Equation (16) accounts for both the acquisition of immunity and the probability of blood-stage infection.
The quantity p h→m (t) is of particular importance since it underpins the coupling between host and vector dynamics; the time evolution of the expected number of infected mosquitoes (and consequently, the FORI) is dependent only on p h→m (t) and several (known) transmission parameters (see Section 3.1). Equation (16) will be of particular use in Section 3.3, where we construct a reduced hybrid transmission model.

Illustrative sample path
A simulated sample path illustrating temporal variation in the relapse rate, superinfection status (that is, the number of co-circulating parasite broods in the bloodstream) and immunity level, is shown in Figure 2, as an extension of previous simulations presented in Mehra et al. [44,49]. We assume that the individual is both immune-and infection-naive at time zero. Over the course of 10 years, they are subject to 17 infective mosquito bites (shown with vertical dashed lines). Temporal variation in the relapse rate ( Figure 2A), which is proportional to the size of the non-latent hypnozoite reservoir, arises from the interplay between hypnozoite replenishment (through mosquito inoculation) and clearance (through either activation or death). Blood-stage infections include both primary infections (triggered immediately upon mosquito inoculation) and relapses (triggered by hypnozoite activation), with temporally proximate reinfection and/or hypnozoite activation events yielding multiple blood-stage infections ( Figure 2B). The discretised immunity level (shown in Figure 2C) likewise fluctuates, with the clearance of each blood-stage infection eliciting a unit boost that is retained for an exponentially-distributed period of time.
The probability of clinical infection (which serves as a correlate of antidisease immunity) decays geometrically with the discretised immunity level ( Figure 2D). For comparison, in Figure 2E, we illustrate the simplest case of our previous (continuous) model of antibody dynamics [44], whereby the clearance of each blood-stage infection is associated with a unit boost of immunity that then decays exponentially (deterministically). Both the discrete ( Figure 2C) and continuous ( Figure 2E) models of immunity yield qualitatively similar results in this case.

Hybrid transmission models: coupling expected host and vector dynamics
We now construct hybrid transmission models [19,40] to couple host and vector dynamics. We restrict our attention to a homogeneously mixing population of humans and mosquitoes. While the size of the human population P H is taken to be fixed (with no age structure), we allow for time-variation in the size of the mosquito population (e.g. due to climactic variation or the implementation of vector-based control). Vector dynamics are detailed Section 3.1.
In Mehra [48], for a simpler model structure -accounting only for superinfection and shortlatency hypnozoite accrual -we began by constructing a Markov population process (with countably many types) to couple host and vector dynamics. Using the work of Barbour and Luczak [16], we then recovered a deterministic (infinite) compartmental model as the functional law of large numbers; that is, we showed that the sample paths of the Markov process converged to a deterministic sample path in the infinite population size limit, granted the number of mosquitoes per human was held fixed. We noted in Mehra [48], however, that an identical deterministic compartmental model could be recovered under a "hybrid approximation", whereby host and vector dynamics are coupled through expected values, as per the construction of Nåsell [19] and Henry [40]. This hybrid construction is the focus of Section 3.2; by regarding  [27]. We consider long-latency hypnozoites, with k = 2 and δ = 1/100 day −1 . Values for the hypnozoite activation α = 1/334 day −1 and death µ = 1/442 day −1 rates are drawn from White et al. [24]. Blood-stage infections (primary and relapse) are assumed to be cleared at constant rate γ = 1/24 day −1 , as per estimates from White et al. [35]. Under the discretised model of immunity, the lifetime of each immunity boost is assumed to be exponentially-distributed with mean duration 1/w = 250 days, with the probability of clinical symptoms (conditional on the presence of blood-stage infection) assumed to decrease by a factor of p c = 0.65 for each increment of immunity.
(A) The rate of relapse αN N L (t).
(B) The number of parasite broods N A (t) + N P (t) co-circulating in the bloodstream.
(C) The discretised immunity level N I (t).
(D) The probability of clinical infection p The antibody level, as per the model introduced in Mehra et al. [44], whereby the clearance of each blood-stage infection is associated with a unit boost of immunity that then decays exponentially (deterministically) with rate 1/250 day −1 .
the within-host PMF as the expected population-level frequency distribution [40], we recover a compartmental model (comprising an infinite-dimensional system of ODEs) that can be viewed as natural extension of the Ross-Macdonald model to allow for superinfection, hypnozoite accrual and immune acquisition. We characterise endemic equilibria for this compartmental model by drawing on results derived at the within-host level (Section 3.2.1), before performing a sensitivity analysis (Section 3.2.2).
To characterise the transient dynamics of the hybridised system, we adopt the strategy we introduced in Mehra [48]. Specifically, we collapse the infinite-dimensional compartmental model into a reduced system of IDEs -with a set of ODEs governing the time evolution of the number of (un)infected and latent mosquitoes over time; and an integral equation governing the (transmission-blocking) immunity-modulated probability of successful human-to-mosquito tranmsmission (Section 3.3). Based on the time evolution of the number of infected mosquitoes under the reduced system of IDEs, we can recover population-level distributions for several quantities of epidemiological interest using our derived within-host distributions.

Vector dynamics: birth, death and infectivity
Here, we characterise the dynamics of the vector population. We assume mosquito-dynamics are described by a Markovian birth-death process, whereby: • Mosquito births follow a time-dependent rate ω(t) (reflecting, for instance, climactic variation) • Mosquito lifetimes are exponentially-distributed with mean duration 1/g; • Each mosquito bites humans (within a population of fixed size P H ) at a potentially timevarying rate β(t) (reflecting, for instance, the relaxation/intensification of vector-based control measures); • Following successful human-to-mosquito transmission (due to a bloodmeal from a bloodstage infected human), a mosquito undergoes sporogony at rate η; • After sporogony has occured, a mosquito remains infective to humans until death.
A mosquito that is undergoing sporogony following successful human-to-mosquito transmission is hereafter considered to be latent.
Denote by I M (t), L M (t), U M (t) the expected number of infected, latent and uninfected mosquitoes respectively at time t. By the Kolmogorov forward differential equations for the Markovian birth-death process governing the vector population (see Appendix C for details), we obtain the system of coupled ODEs governing the time evolution of I M (t), L M (t), U M (t) as a function of the probability of humanto-mosquito transmission p h→m (t) per bloodmeal.

A countable system of ODEs
Under a hybrid approximation, we seek to couple expected host and vector dynamics [19,40].
Here, we recall two key observations: • As we noted in Section 2.2, for a human population of fixed size P H , the system of ODEs given by Equation (1) governs the expected proportion of humans in each hypnozoite/infection/immunity state as a function of the FORI [40].
• The time evolution of the expected number of (un)infected mosquitoes in the population is governed by Equations (17) and (19) conditional on the probability of successful humanto-mosquito transmission per bloodmeal.
To construct a hybrid transmission model, as per the approach of Nåsell [19], it remains to characterise the FORI, and the probability of human-to-mosquito transmission.
We assume that mosquito-to-human transmission occurs with fixed probability p m→h when a human is bitten by an infected mosquito. We do not allow the parameter p m→h to be modulated by immunity, that is, we do not account for the acquisition of pre-erythrocytic immunity. As a function of the number of infected I M (t) mosquitoes over time, the FORI can therefore be written Likewise, as a function of the proportion of humans H i 1 ,...,i k ,i N L ,j,k (t) with i hypnozoites in state ∈ {1, . . . , k, N L}; j co-circulating broods in the bloodstream; and immunity level k at time t, the probability of successful human-to-mosquito transmission p h→m (t) can be written as per the model of transmission-blocking immunity detailed in Section 2.5.
Then following a similar approach to Nåsell [19], we recover the countable system of ODEs where we have used Equation (1) and Equations (17) to (20). A schematic of this model is provided in Figure 3.
The compartmental model given by Equations (21) and (24) represents a natural extension of the Ross-Macdonald framework to allow for hypnozoite accrual, superinfection and transmissionblocking immunity. We can also view Equations (21) and (24) as an extension of the transmission model proposed by White et al. [24] to allow for long-latency hypnozoites, immunity and explicit superinfection dynamics (as opposed to the approximation of superinfection through an appropriate recovery rate). A summary of model parameters, and their respective interpretations, is provided in Table 1.   [19,40]. Here, the probabilistic distribution of the open network of infinite server queues governing within-host dynamics (Section 2.1) is re-interpreted as the expected proportion of humans in each hypnzoite/superinfection/immunity state. The coupling of host and vector dynamics is predicated on the force of reinfection λ(t) (Equation (20)), which is a function of the number of infected mosquitoes at time t; and the probability of successful human-to-mosquito transmission p h→m (t) per bloodmeal (Equation (16)), which is modulated both by the prevalence of blood-stage infection and the distribution of immunity in the human population at time t.  [27] γ Rate at which each blood-stage infection is cleared White et al. [35] w Rate at which each immune increment/boost is lost Gething et al. [14] ω(t) Mosquito birth rate at time t

The stationary solution
Here, we seek to characterise steady state solutions to the system of ODEs given by Equations (21) to (24). As such, we restrict ourselves to a setting where: • The bite rate per mosquito β(t) = β remains constant over time; and • The mosquito population size I M (t) + L M (t) + U M (t) = P M is fixed, that is, the birth rate ω(t) = g exactly balances the (constant) death rate.
Denote by H * i 1 ,...,i k ,i N L ,j,k , U * M , L * M I * M the stationary solution to the system of IDEs given by Equations (21) to (24), recovered by setting all time derivatives to zero.
We focus on the quantities H * i 1 ,...,i k ,i N L ,j,k and I * M since the overarching effect of sporogony is to introduce a scaling factor g/(g + η) in the fraction of mosquitoes that, in the event of successful mosquito-to-human transmission, survive latency to transition from an uninfected to infected state. Setting the time derivative in Equation (22) to zero and using the assumption of a fixed mosquito population size, we can formulate the number of latent L * M and uninfected U * M at steady state as functions of I * M and P M as follows: We observe that the disease (and immunity) free equilibrium H * 0,...,0,0,0,0 = 1, I * M = 0 always exists, Here, we seek to characterise the existence of endemic equilibrium solutions.
We begin by setting the time derivatives in Equations (22) and (24) to zero to yield We then note that Equation (21) is precisely the set of Kolmogorov forward differential equations for the open network of infinite server queues introduced in Section 2. The PGF for the stationary limiting distribution of this queueing network, given a constant FORI λ(t) = βp m→h I * M /P H , can be recovered by taking the limit t → ∞ in Equation (11). Therefore, using Equation (16) -which we derived from Equation (11) in Section 2.5 -we suggest that Using a simple geometric argument (Appendix D), we can show that Equations (25) and (26) have at most one non-zero intersection (corresponding to an endemic equilibrium solution), and that this intersection exists if and only if Assuming that R 0 > 1 (Equation (27)), an endemic equilibrium solution necessarily exists. As a function of the FORI λ * = βp m→h I * M /P H at the endemic equilibrium, which is a function of the non-trivial solution I * M ∈ (0, P M ] to Equations (25) and (26), we can recover population-level distributions for various quantities of epidemiological interest using the stationary limiting PGF recovered by setting λ(t) = λ * for all t ≥ 0 and taking the limit t → ∞ in Equation (11).
Relevant formulae (based on the derivations presented in Mehra et al. [49])) are provided in Appendix B.

Sensitivity analysis for endemic equilibrium solutions
We now perform a sensitivity analysis for the endemic equilibrium solutions. In Section 3.2.2.1, we examine endemic equilibrium solutions in the absence of immunity. Endemic equilibria, allowing for transmission-blocking and antidisease immunity, are detailed in Section 3.2.2.2.

Short-latency vs long-latency strains in the absence of transmission-blocking immunity
We begin by examining steady state solutions for both short-latency (k = 0) and long-latency (k > 0) strains in the absence of transmission-blocking immunity (p tb = 1). Figure 4A depicts a sensitivity analysis for the bifurcation parameter R 0 (Equation (27)) with respect to the hypnozoite activation α, death µ and latency k parameters. Recall that an endemic equilibrium exists, and is unique, if R 0 > 1; if R 0 < 1, only the disease-free equilibrium exists. The bifurcation boundary R 0 = 1 for parameters (α, µ) is shown in white in Figure 4A.
In the absence of hypnozoite accrual (ν = 0), no endemic equilibria exist for the set of transmission parameters considered here; the existence of endemic equilibrium solutions is therefore contingent on the relapse burden. The interplay between hypnozoite activation α and death We set P M /P H = 1.2, p 0 = 0.25 and parameters γ, ν, g = ω(t), β, p m→h as per Table 1.
µ rates governs the expected number of relapses per bite α/(α + µ) [24]. As such, when the hypnozoite activation rate α is low relative to the hypnozoite death rate µ, there are insufficient relapses to sustain transmission and no endemic equilibrium solution exists, that is, R 0 < 1 ( Figure 4A).
In the case of short-latency strains (k = 0), we further observe that excessively high activation rates α preclude the existence of endemic equilibrium solutions (that is, yield R 0 < 1) ( Figure   4A); similar observations have been posited by White et al. [27] and Anwar et al. [46], albeit in the absence of superinfection. Without an enforced dormancy period, an elevated activation rate α gives rise to a high risk of relapse immediately after each infective bite. The rapid depletion of the hypnozoite reservoir following each bite -driven by temporally proximate hypnozoite activation events -leads to a divergence in relapse risk conditional on status of blood-stage (super)infection ( Figure 4B2). To justify why a high relapse rate for (blood-stage) superinfected individuals is a weaker driver of onward transmission than a high relapse rate for blood-stage uninfected individuals, we observe that the expected time to clearance for j parasite broods is γ(1 + 1/2 + · · · + 1/j); as such, an additional relapse for an individual with m pre-existing broods in their bloodstream increases the expected time to (blood-stage) recovery by an increment of γ/(m + 1). We deduce that the stratification of relapse risk, conditional on the status of blood-stage infection, is driven by the time to the most recent infective bite in the case of fast-activating short-latency hypnozoites: while recently-inoculated individuals experience a high burden of both liver-and blood-stage infection, there is a limited burden of liver-and blood-stage infection between successive mosquito bites, yielding a population-level reduction in the overall burden of blood-stage infection ( Figure 4B1). Hence, for elevated activation rates α, there is a limited window of time following each infective bite for which an individual remains blood-stage infected, and therefore infective to mosquitoes; if the activation rate α is sufficiently high, then these windows of human-to-mosquito infectivity may be insufficient to sustain transmission in the steady state, in which case R 0 < 1 and no endemic equilibrium solution exists ( Figure 4A).
For long-latency strains (k > 0), stochasticity in the enforced dormancy period prevents excessive overlap between hypnozoite activation events arising from the same bite, thereby reducing the sensitivity of the endemic equilibrium burden of blood-stage infection to elevated hypnozoite activation rates α ( Figure 4B1). Decreasing the variance of the dormancy period k/δ 2 , whilst fixing the expected duration k/δ, would presumably increase the sensitivity of endemic equilibria to the hypnozoite activation rate α, since hypnozoites would be more likely to emerge from dormancy at similar times. We observe that the assumption of independent dormancy, introduced in Mehra et al. [42], underpins this observation for long-latency strains; the collective dormancy assumption of White et al. [24] -whereby synchronicity in the latency phase means that hypnozoites established through the same infective bite emerge collectively from dormancy -leads to greater sensitivity of endemic equilibrium solutions to elevated hypnozoite activation rates α. Indeed, under a 'binary' hypnozoite model predicated implicitly on the assumption of collective dormancy, White et al. [27] predict stronger constraints on the hypnozoite activation rate α than we predict in Figure 4 under the assumption of independent dormancy. Elevated hypnozoite activation rates α, however, give rise to a stratification of relapse risk by superinfection status, even in the case of long-latency strains (k > 0) ( Figure 4B2). In the absence of hypnozoite death (that is, µ = 0), the burden of blood-stage infection at the endemic equilibrium is maximised for low hypnozoite activation rates α ( Figure 4B1), which yield broad temporal relapse distributions for each infective bite, and a population-level relapse risk that does not vary strongly by superinfection status ( Figure 4B2). For non-zero death rates µ > 0, however, hypnozoite death during the enforced dormancy period -during which activation is prohibited -serves as a key constraint. As such, elevated activation rates α (up to a point) yield an increasing burden of blood-stage infection for long-latency strains (k > 0, Figure 4B1), even as the risk of relapse stratified by superinfection status becomes more unbalanced and proportionately higher for individuals with pre-existing blood-stage infections ( Figure 4B2).

Short-latency strains with transmission-blocking immunity
We now perform a sensitivity analysis for endemic equilibrium solutions allowing for the acquisition of transmission-blocking and antidisease immunity ( Figure 5). Here, we restrict ourselves to short-latency strains (k = 0). For a fixed set of hypnozoite activation and death rates (α, µ), long-latency strains (k > 0) yield similar qualitative patterns as a function of the immunity parameters w, p tb and p c . In the absence of transmission-blocking and clinical immunity (that is, p tb = p c = 1), we revert to the setting examined in Section 3.2.2.1; we highlight this case with closed circles in Figure 5.
A two-way sensitivity analysis, with respect to parameters w and (1−p tb ), is shown in Figure 5A.
Recall that the probability of human-to-mosquito transmission (when an uninfected mosquito takes a bloodmeal from a blood-stage infected human) decays geometrically, with factor p tb , as a function of an individual's immunity level (Section 2.5); since protection rises as p tb → 0, we Here, we set p 0 = 0.65 and P M /P H = 1.2, and parameters α, µ, γ, ν, g = ω(t), β, p m→h as per Table 1. think of (1−p tb ) as a transmission-blocking protection parameter. In contrast, the parameter 1/w governs the time scale for which immunity is retained, with the limiting case w = 0 corresponding to a scenario where immunity is never lost. The endemic equilibrium prevalence of blood-stage infection ( Figure 5A1) decreases both as: • immunity becomes longer-lived (that is, w → 0), whereby a larger subset of an individual's infection history is expected to contribute to their current immunity level; and • the protective effect associated with each cleared blood-stage infection is augmented (that is, (1 − p tb ) → 1).
Mitigation of the blood-stage infection burden in light of transmission-blocking immunity, however, necessarily limits exposure; reductions in the the prevalence of blood-stage infection are therefore accompanied by reductions in the population-level distribution of immunity. The mean immunity level at the endemic equilibrium therefore decreases, even as the rate of immune decay w decreases, and immunity accrues over a larger time scale ( Figure 5A2). Likewise, augmenting the transmission-blocking effect of each immunity increment (1 − p tb ) -whereby the probability of human-to-mosquito is suppressed strongly, even at low immunity levels -leads to a reduction in the mean immunity level at the endemic equilibrium ( Figure 5A2).
In particular, we see a substantially reduced burden of (blood-stage) superinfection at the endemic equilibrium as the transmission-blocking effect of each immune increment (1 − p tb ) increases ( Figure 5B). The suppression of superinfection explains the stronger decay in the mean immunity level ( Figure 5A2), relative to the prevalence of blood-stage infection ( Figure 5A1), as the transmission-blocking protection parameter (1 − p tb ) is strengthened: since the clearance of each primary infection and relapse yields an immunity boost, irrespective of temporal overlap with other blood-stage infections, superinfection is an important driver of acquired immunity.
We also observe a trade-off between transmission-blocking and antidisease immunity. The mitigation of transmission as the transmission-blocking protection parameter (1 − p tb ) is augmented leads to a lower population-level distribution of immunity at the endemic equilibrium ( Figure   5A2). If the distribution of immunity at the endemic equilibrium is sufficiently reduced, then an increasingly strong transmission-blocking effect per immune increment (1 − p tb ) can give rise to an increasing prevalence of clinical infection at the endemic equilibrium ( Figure 5C), even as the burden of blood-stage infection continues to decline ( Figure 5B).

A reduced system of IDEs to study transient behaviour
The hybrid transmission model given by Equations (21) to (24) yields population-level dynamics of superinfection, the hypnozoite reservoir and acquired immunity. However, the countable system of ODEs given by Equations (21) to (24) is not necessarily readily amenable to numerical solution; truncating the system at reasonable endpoints could yield thousands of coupled ODEs, since the size of the hypnozoite reservoir (and, by extension, the immunity level) could reasonably be expected to range up to 30 in moderate to high transmission settings (see Figure   4 of White et al. [24] and Figure 6 of Anwar et al. [46]).
Here, we propose a reduced system of integrodifferential equations (IDEs) to couple host and vector dynamics, following the approach detailed in Mehra [48]. In particular, we observe that: • The dependence of the human population on vector dynamics can be distilled into the FORI, which is proportional to the number of infected mosquitoes in the population I M (t) (see Equation (20)).
• The dependence of the vector population on the state of the human population can be distilled into the probability of successful human-to-mosquito transmission p h→m (t) when an uninfected mosquito bites any human in the population; note that the quantity p h→m (t) accounts for both the prevalence of blood-stage infection and the distribution of (transmission-blocking) immunity within the human population.
At time t = 0, we make the assumption that each individual in the human population (of fixed size P H ) has immunity level zero; carries no hypnozoites; and harbours no ongoing blood-stage infections, whereby p h→m (0) = 0. As such, we consider the introduction of a number of infected mosquitoes into an otherwise infection-and immune-naive human population. As a function of the FORI, the probability of successful human-to-mosquito transmission p h→m (t) is then governed by the integral given in Equation (16). Likewise, as a function of the probability of successful human-to-mosquito transmission p h→m (t) per bloodmeal, the time evolution of the expected number of infected mosquitoes over time I M (t) is governed by the coupled system of ODEs given by Equations (22) to (24), which also captures the time evolution of the expected number of latent L M (t) and uninfected U M (t) mosquitoes over time.
Coupling expected host and vector dynamics under a hybrid approximation thus yields the system of IDEs with initial condition I M (0), L M (0), U M (0) ≥ 0, where we have used Equations (16) and (22) to (24). Recall that • p h,A (x) denotes the probability that a hypnozoite has activated to give rise to a relapse that is ongoing time x after inoculation (Equation (34)).
• p h,I (x) denotes the probability that immune memory has been gained (following the clearance of a relapse) time x after a hypnozoite is established in the liver (Equation (35)).
• p p,A (x) denotes the probability that a primary infection is ongoing time x after onset (Equation (10)).
• p p,I (x) denotes the probability that immune memory has been gained time x after the onset of a primary infection (Equation (10)).
Interpretations for each transmission/within-host parameter are detailed in Table 1.
The system of IDEs given by Equations (28) to (31)  Observe that the integral equation governing p h→m (t) (Equation (31)) -which we derived using the within-host PGF given by Equation (11) -satisfies the Kolmogorov forward differential equations for the within-host queueing structure, which in turn constitutes the expected frequency distribution of the human population (Equation (21)), granted an individual is both immune-and infection-naive at time zero and the FORI at time t is β(t)p m→h I M (t)/P H . As such, under the assumption that the human population is both immune-and infection-naive at time zero, the time evolution of the FORI and the probability of human-to-mosquito transmission per bloodmeal are equivalent under the system of IDEs given by Equations (28) to (31), and the countable system of ODEs given by Equations (21) to (24).
While the quantities I M (t), L M (t), U M (t) and p h→m (t) are sufficient to couple host and vector dynamics, we ultimately seek to characterise population-level distributions for quantities of epidemiological interest. We note, however, that the complete population-level distribution of superinfection, immunity and hypnozoite states can be recovered conditional on the FORI using the results derived in Mehra et al. [49]. To reiterate the premise of the hybrid approximation, the population-level transmission models discussed here have been constructed by casting the within-host probabilistic distribution as the population-level frequency distribution [40]. Given the time evolution of the FORI λ(t) = β(t)p m→h I M (t)/P H derived from the system of IDEs given by Equations (28) to (31), we can recover population-level distributions for quantities of epidemiological interest using the formulae derived in Mehra et al. [49].
The methodology adopted here uses an integral system, within which we can enforce timedependence in the bite rate per mosquito β(t) and the mosquito birth rate ω(t). In Section 3.3.1 below, we present illustrative results for two scenarios: a constant transmission setting, where all transmission parameters are fixed (Section 3.3.1.1); and a seasonal transmission setting, with a sinusoidal mosquito birth rate ω(t) (Section 3.3.1.2). Vector-based control interventions represent a natural extension [11,13,34], but are not presented here.

Illustrative results for the reduced system of IDEs
To recover transient host and vector dynamics, we solve the system of IDEs given by Equations (28) and (31) numerically, using Euler's method (for the ODEs given by Equations (28) to (30)) and the trapezoidal rule (for the integral given by Equation (31)) with a fixed time step; this procedure is a variation of the algorithm proposed by Anwar et al. [46]. As a function of the FORI derived from Equations (28) and (31), we recover the time evolution of several quantities of epidemiological interest. Relevant formulae (as per Mehra et al. [49]) are provided in Appendix B, including: • the mean and variance for the size of the (non)-latent hypnozoite reservoir (Equations (36) and (37)); • the PMF for the number of co-circulating blood-stage broods (Equations (38) to (40)); • the relapse rate conditional on the blood-stage infection status (Equations (45) to (49)); • the distribution of immunity, as quantified by the mean and variance of the discrete immunity levels (Equations (50) and (51)); • the prevalence of clinical infection (Equation (15)).

Non-seasonal transmission
We begin by considering host and vector dynamics in the absence of seasonality (Figure 6).
At time zero, we consider the introduction of several infected mosquitoes into an (blood-and liver-stage) infection and immune naive human population. Predicted endemic equilibrium solutions, given by the unique non-trivial solution to Equations (26) and (25), are shown with dashed blue lines for the FORI (Figures 6A1/B1) and the immunity-modulated probability of human-to-mosquito transmission ( Figures 6A5/B5).
Illustrative dynamics for short latency strains (k = 0) are shown in Figure 6A. The low initial level of infection in the mosquito population constrains the transmission intensity at early time points. Prior to the acquisition of extensive transmission-blocking immunity -with relatively low immunity levels harboured for a year following the introduction of infected mosquitoes into an immune-naive human population ( Figure 6A6) -human-to-mosquito transmission remains comparatively unconstrained, leading to a sustained increase in the FORI ( Figure 6A1), and consequently, the blood-stage infection burden ( Figure 6A2) and the size of the hypnozoite reservoir ( Figure 6A7). A pronounced rise in the prevalence of blood-stage infection during this early period leads to an increase in both the probability of human-to-mosquito transmission ( Figure   6A5) and the prevalence of clinical infection ( Figure 6A3). Intensified transmission, however, is accompanied by the sustained acquisition of immunity (Figures 6A6), which eventually mitigates the probability of human-to-mosquito transmission ( Figure 6A5), leading to a reduction in the FORI ( Figure 6A5), as well as a slight reduction in the burden of (clinical) blood-stage infection ( Figures 6A2, A3). These transient effects eventually subside, and for the set of parameters considered here, the predicted endemic equilibrium (obtained by numerically solving Equations (25) and (26), and indicated with dotted lines blue lines) is reached within approximately 4 years.
Analogous results for long-latency strains (k > 0) are shown in Figure 6B. P H and (A5/B5) probability of human-to-mosquito tranmission p h→m (t) are governed by the system of IDEs given by Equations (28) and (31). Endemic equilibrium solutions for the FORI βp m→h I M * P H and human-to-mosquito transmission probability p * h→m , given by the non-trivial solution to Equations (26) and (25), are indicated with dashed blue lines. All other quantities are calculated as a function of the numerical solution for the FORI Model parameters α, µ, γ, ν, w, g, η, β, p m→h , p tb as per Table 1, with p 0 = 0.65. size for approximately one year ( Figure 6B7); as such, single-brood primary infections dominate the infection burden for an extended period of time ( Figures 6B2 and 6B4) relative to short-latency strains. In tandem with the emergence of hypnozoites from dormancy, relapses eventually drive up the burden of (clinical) blood-stage infection (Figures 6B2 and 6B4). As for short-latency strains (k = 0), the acquisition of transmission-blocking immunity eventually mitigates onward human-to-mosquito transmission ( Figure 6B5), leading to a slight reduction in transmission intensity before the predicted endemic equilibrium (obtained by numerically solving Equations (25) and (26), and indicated with dotted lines blue lines) is reached six years after the introduction of infected mosquitoes into an infection-and immune-naive human population.

Seasonal transmission
To allow for seasonality, arising, for instance, from external climactic variation, we impose sinusoidal forcing (with period one year) on the mosquito birth rate. Illustrative dynamics for short-latency (k = 0) and long-latency (k > 0) strains are shown in Figures 7A and 7B respectively. With the imposition of seasonal forcing, infection levels within both humans and mosquitoes exhibit oscillations (with period one year) that eventually stabilise around a steady mean. The nature of these oscillations within a season, however, varies between short-and long-latency strains. Oscillations in the FORI ( Figure  Model parameters α, µ, γ, ν, w, g, η, β, p m→h , p tb as per Table 1, with p 0 = 0.65.

Discussion
The interplay between the hypnozoite reservoir, superinfection and acquired immunity is a key aspect of the epidemiology of P. vivax. Here, we propose a hybrid transmission model for P. vivax, accounting for hypnozoite accrual, (blood-stage) superinfection and the acquisition of transmission-blocking and antidisease immunity. To capture within-host dynamics as a function of the FORI, we extend the open network of infinite server queues constructed in Mehra et al. [49] to embed a discretised version of the antibody model we introduced in Mehra et al. [44].
By deriving the joint PGF for the state of the queueing network (Equation (11)), we obtain an analytic description of within-host dynamics in a general transmission setting. To couple host and vector dynamics, we adopt the hybrid approximation of Nåsell [19] under which probabilistic within-host distributions are cast as expected population-level proportions [40]. We thus recover a deterministic compartmental model, comprising a countably infinite system of ODEs (Equations (21) and (24)), which can be viewed as a natural extension of the Ross-Macdonald framework. For a simpler system with countably many states, we demonstrated the equivalence of the hybrid approximation to the functional law of large numbers [16] for an appropriate Markov population process in Mehra [48].
We draw on distributions derived at the within-host level [49] to characterise both the transient and steady state behaviour of this compartmental model. In particular, following the approach we developed in Mehra [48], we derive a reduced system of IDEs governing the time evolution of the number of (un)infected and latent mosquitoes; and the immunity-modulated probability of human-to-mosquito transmission (Equations (28) to (31)). As a function of the FORI predicted under this reduced system of IDEs -which is equivalent to the complete compartmental model granted the human population is initially immune-and infection-naive -we recover complete population-level distributions for various quantities of epidemiological interest, using the formulae derived in Mehra et al. [49] (see Appendix B). By drawing on the within-host queueing models we introduced in Mehra et al. [44,49], and the construction developed in Mehra [48], we circumvent the practical constraints that have previously limited the tractability of hypnozoite density models [34].
Our model, to the best of our knowledge, provides the most complete description of superinfection, immunity and hypnozoite dynamics for P. vivax thus far, while remaining readily amenable to numerical solution and analysis. In Mehra [48], we developed a framework to capture the dynamics of (short-latency) hypnozoite accrual and superinfection, addressing a gap in the lit-erature with respect to the rigorous analysis of superinfection; we have extended the framework of Mehra [48] in the present manuscript to allow for greater biological realism, namely, acquired immunity and long-latency phenotypes. The joint population-level dynamics of the hypnozoite reservoir and acquired immunity have been previously examined by White et al. [34]. The construction of White et al. [34] is predicated on a continuous age-and exposure-dependent immunity level, which is subsequently mapped (using Hill functions) to correlates of antidisease immunity (that is, a reduced probability of clinical infection) and antiparasite immunity (including accelerated parasite clearance and mitigated parasite densities, manifesting in a reduced probability of detection via light microscopy). Here, we instead consider a discretised exposure-dependent immunity level, which we map to correlates of antidisease immunity and transmission-blocking immunity (that is, a reduced probability of human-to-mosquito transmission) under the assumption of geometric decay. While White et al. [34] explicitly account for treatment, and age structure and heterogeneity in the human population, we restrict our attention to a homogeneous human population in the absence of treatment and age structure.
Unlike White et al. [34], however, we monitor hypnozoite densities rather than broods (thereby capturing variation in parasite inocula across bites), in addition to the population-level distribution of superinfection; our framework also holds for long-latency hypnozoite strains, unlike the framework of White et al. [34] which is restricted to short-latency strains.
While most previous hypnozoite 'batch' and 'density' models [24,34,46] have relied on numerical solution to characterise steady state solutions, our sensitivity analyeses are informed by the within-host distributions derived in Mehra et al. [49]. We recover a threshold phenomenon for the hybrid model, deriving a bifurcation parameter (Equation (27)) governing the existence of endemic equilibria. In the absence of transmission-blocking immunity (p tb = 1) and mosquito latency (1/η = 0), we were able to perform an asymptotic stability analysis in Mehra [48] for the first-order IDE governing the time-evolution of the FORI using the stability criterion of Brauer [2]; the imposition of transmission blocking immunity (p tb < 1) or mosquito latency (1/η > 0) yields a higher-order system of IDEs governing the FORI, for which we are unaware of asymptotic stability criteria.
The transient and stationary results presented in this manuscript are underpinned by analyticity at the within-host level, which, in turn, is predicated on the assumption that each hypnozoite/infection is governed by an independent stochastic process. The assumption of independent, spontaneous hypnozoite activation in line with the 'genetic clock' hypothesis, as implemented by White et al. [24], is critical to our construction: external triggers of hypnozoite activation (e.g. febrile illness, arising from parasitic or bacterial infections [21], particularly P. falciparum [37]) necessarily introduce synchronicity between activating hypnozoites, thereby violating the assumption of independent hypnozoite behaviour. Our model does not readily accommodate interactions between concurrent hypnozoites/infections, for instance, competition between co-circulating parasite broods [9]. Antiparasite immunity (manifest in the modulation parasite clearance rates [30]) and pre-erythrocytic immunity [18], which render hypnozoite/infection dynamics dependent on the infection history, are likewise intractable. On a population-level, our model is constrained by the assumption of homogeneity for the human population. We do not account for age structure or demography within the human population; heterogeneity in the risk of relapse and immunity levels across individuals is driven purely by stochastic fluctuations, rather than differences in the time over which the hypnozoite reservoir has been accrued and immunity has been acquired.
Our formulation of immunity, moreover, is non-mechanistic and subject to a number of simplifying assumptions. Adopting a discretised version of the model presented in Mehra et al. [44], we assume that the clearance of each blood-stage infection is accompanied by an immunity boost with unit magnitude and an exponentially-distributed lifetime. Empirical characterisation of antibody titres, however, has revealed substantial heterogeneity in the magnitude of antibody boosts across successive infections, and the temporal distribution of antibody boosts associated with different antigens [23]. A key omission in our model is strain specificity; while homologous challenge yields a strong immune response, immune protection following heterologous challenge is contingent on the degree of cross-reactivity between strains [18]. As such, the discretised immunity level considered here largely serves as a proxy for 'recent' exposure to blood-stage infection, with the immune decay parameter w governing the time scale on which immunity is retained.
Nonetheless, in capturing the interplay between hypnozoite accrual, superinfection and acquired immunity -and providing, to the best of our knowledge, the most complete population-level distributions for a range of epidemiological values -our model provides insights into important, but poorly understood, epidemiological features of P. vivax, with natural extensions to explore the consequences of control and elimination strategies.

B Population-level distributions for quantities of epidemiological interest
As a function of the FORI λ(t) in the interval τ ∈ [0, t), we can recover distributions for various quantities of epidemiological interest at time t using the results derived in Mehra et al. [49].
In the case of the open network of infinite server queues described in Section 2, we assume a functional form for λ(t), wherein we interpret these quantities as probabilistic distributions at the within-host level. In the context of the hybrid transmission models constructed in Section 3, where we capture the time evolution of λ(t) by coupling host and vector dynamics, these quantities are instead interpreted as population-level proportions [40]. We adopt the same notation for both situations, with N s (t) denoting the marginal distribution for the number of hypnozoites/infections in state s ∈ S at time t on either the within-host or population-level.

B.3 Relapse rate conditional on blood-stage infection status
The joint PGF for the size of the non-latent hypnozoite reservoir N N L (t) and the number of co-circulating parasite broods M I (t) := N A (t) + N P (t) can be written where we have used Equation (11).
The (unconditional) relapse rate, which is proportional to the expected size of the non-latent hypnozoite reservoir, can be written as in Equation (39) of Mehra et al. [49].
Conditional on the absence of blood-stage infection (M I (t) = 0), the relapse rate is given by as in Equation (78) of Mehra et al. [49].
Likewise, by Xekalaki [4], we recover the relapse rate conditional on a single-brood blood-stage infection (M I (t) = 1), as well as a double-brood blood-stage infection (M I (t) = 2) (1 + νp h,A (t − τ )) 4 dτ By the law of total expectation, we can recover the relapse rate conditional on a blood-stage using Equations (38) to (40) and (43) to (47).  (40) and (41) of Mehra et al. [49]. The relative contribution of relapses to the infection burden (that is, the proportion of blood-stage infections that encompass at least one relapse) follows readily from the quantities P (N P (t) = 0), P (N A (t) = 0) and P (N A (t) = N P (t) = 0) (that is, the probability that each individual has neither a primary infection, nor a relapse, at time t).

B.5 Distribution of immunity
From the joint PGF given by Equation (11), we recover the PGF governing the population-level distribution of immunity I(t) at time t E z N I (t) ] = G(t, z 1 = 1, . . . , z k = 1, z N L = 1, z A = 1, z D = 1, z C = 1, z I = z, z P = 1) We can thus compute the expected immunity level The complete population-level distribution of immunity can be recovered using Faa di Bruno's formula and Leibiniz's integral rule, using a similar approach to Mehra et al. [49] (see Equations (82) and (83) of Mehra et al. [49], which have an analogous functional form to that consided here).

C Expected vector dynamics
Here, we characterise the dynamics of the vector population, as a function of the probability of successful human-to-mosquito transmission p h→m (t). The structure of the Markovian birthdeath process governing the state of the vector population is described in Section 3.1.
Define the generating function Using the ODEs given by Equation (52), we recover a PDE for the generating function F ∂F ∂t Then for s ∈ {i, l, u}, by differentiating the PDE given by Equation (53) with respect to z s and evaluating the resultant expression at z i = z l = z u = 1, we recover precisely the system of ODEs given by Equations (17) to (19)  From Equations (13) and (14), we observe that On the other hand, observing that F 2 , F 2 are continuous; F 2 (0) > 0, F 2 (0) < 0 and As such, R 0 > 1 is a sufficient and necessary condition for the existence of an endemic equilibrium; when R 0 > 1, the endemic equilibrium solution is unique.