Superinfection and the hypnozoite reservoir for Plasmodium vivax: a general framework

A characteristic of malaria in all its forms is the potential for superinfection (that is, multiple concurrent blood-stage infections). An additional characteristic of Plasmodium vivax malaria is a reservoir of latent parasites (hypnozoites) within the host liver, which activate to cause (blood-stage) relapses. Here, we present a model of hypnozoite accrual and superinfection for P. vivax. To couple host and vector dynamics for a homogeneously-mixing population, we construct a density-dependent Markov population process with countably many types, for which disease extinction is shown to occur almost surely. We also establish a functional law of large numbers, taking the form of an infinite-dimensional system of ordinary differential equations that can also be recovered by coupling expected host and vector dynamics (i.e. a hybrid approximation) or through a standard compartment modelling approach. Recognising that the subset of these equations that model the infection status of the human hosts has precisely the same form as the Kolmogorov forward equations for a Markovian network of infinite server queues with an inhomogeneous batch arrival process, we use physical insight into the evolution of the latter process to write down a time-dependent multivariate generating function for the solution. We use this characterisation to collapse the infinite-compartment model into a single integrodifferential equation (IDE) governing the intensity of mosquito-to-human transmission. Through a steady state analysis, we recover a threshold phenomenon for this IDE in terms of a parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0$$\end{document}R0 expressible in terms of the primitives of the model, with the disease-free equilibrium shown to be uniformly asymptotically stable if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0<1$$\end{document}R0<1 and an endemic equilibrium solution emerging if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0>1$$\end{document}R0>1. Our work provides a theoretical basis to explore the epidemiology of P. vivax, and introduces a strategy for constructing tractable population-level models of malarial superinfection that can be generalised to allow for greater biological realism in a number of directions.


Introduction
Malaria is a parasitic, vector-borne disease with a staggering pulic health burden.The overwhelming majority of malaria cases (98%) are attributed to the parasite Plasmodium falciparum, particularly in Africa [49].Plasmodium vivax, however, exhibits a broader geographical distribution, driving much of the malaria burden in South East Asia, the Americas, the Western Pacific and the Eastern Mediterranean [49].While an estimated 4.5 million malaria cases were attributed to P. vivax in 2020 alone [49], morbidity arising from P. vivax infections remains "obscure and insidious" [46].
The transmission of malaria parasites to humans is mediated by Anopheles mosquito vectors.
During the course of a bloodmeal, an infected mosquito can transmit parasites (sporozoites) to a human host.Following a period of liver-stage development (exoerythrocytic schizogony), parasites are released into the bloodstream, giving rise to a blood-stage infection that is sustained by the replication of parasites in invaded red blood cells [45].A key epidemiological characteristic of malaria is the phenomenon of superinfection.Since the circulation and replication of (pre-existing) parasites in the bloodstream does not preclude further blood-stage infection, an individual can concurrently harbour multiple co-circulating broods of blood-stage parasites.In the context of P. falciparum, we define each blood-stage 'brood' to derive from a single infective bite.The interpretation of a blood-stage 'brood' for P. vivax is more nuanced, in light of its ability to cause relapsing infections following the accrual of a "hypnozoite reservoir" [31,35].
Notably, a P. vivax parasite (sporozoite) injected into a human host has two possible fates: it either gives rise to a primary (blood-stage) infection within approximately 9 days of the bite itself [36], or develops into a hypnozoite [23].Hypnozoites undergo indeterminate latency periods, often lasting weeks or months, during which they are undetectable using standard techniques [48].The activation of each hypnozoite, however, can trigger an additional blood-stage infection, known as a relapse [23].For P. vivax, each primary infection and relapse is defined to comprise a separate blood-stage brood.
We define the multiplicity of broods (MOB) to be the number of co-circulating blood-stage broods in a host at a given point in time, with superinfection taken to be a collective term for blood-stage infections with MOB> 1.As such, superinfection arises from temporally proximate reinfection (that is, infective bites) and, in the case of P. vivax, hypnozoite activation events.
With epidemiological data indicating the preponderance of relapses over primary infections [42] and evidence of P. vivax superinfection even in the absence of reinfection [39], analysis of the statistics of superinfection for P. vivax warrants careful consideration of the hypnozoite reservoir.
The classical framework of malarial superinfection for P. falciparum, proposed initially by Macdonald [1] and formulated mathematically by Bailey [2], assumes independent clearance of each brood, without imposing an upper bound on the MOB.Under this setting, a natural construction to describe the within-host dynamics of superinfection is an infinite-server queue with a time dependent arrival rate given by the intensity of mosquito-to-human transmission [3,8,24,43].Hereafter, we refer to the mosquito-to-human transmission intensity, as quantified by the infective bite rate per human, as the force of reinfection (FORI).For P. vivax, introducing the additional assumption that the dynamics of each hypnozoite are governed by independent stochastic processes [35], we have recently extended this idea to characterise within-host and superinfection dynamics using an open network of infinite-server queues with geometricallydistributed batch arrivals at a time-dependent rate given by the FORI [47,52].In Mehra et al. [52], we derive a time-dependent generating function for the state of the queueing network, which can be inverted analytically to recover marginal distributions for MOB and the hypnozoite burden, amongst other quantities of epidemiological interest, on a within-host scale.
Mathematical modelling of superinfection and hypnozoite dynamics at the population-level can be challenging.As we note in [52], a common approach in the construction of transmission models of P. vivax has been to consolidate hypnozoite carriage into a single state, with an accompanying parametric form for the time to first relapse for hypnozoite-positive individuals [19,28,34,32,37,38].The binarisation of hypnozoite carriage, however, obscures the relationship between transmission intensity and hypnozoite accrual: there is no variation in the risk of relapse based on hypnozoite density, the limitations of which are discussed in [52].Explicit variation in the hypnozoite burden is captured in the deterministic models of White et al. [35,40] and Anwar et al. [50].While the 'batch' model of [40] captures 'broods' of hypnozoites in the liver, it ignores variation in parasite inoculum sizes.On the other hand, both [35,50] take hypnozoite densities into account.To account for superinfection, [35] employs a "pseudoequilibrium approximation" for the (blood-stage) infection recovery rate, adopting the functional form derived by Dietz, Molineaux, and Thomas [3].This functional form, which was derived in the absence of hypnozoite accrual, is embedded in multiple transmission models for P. falciparum as a proxy for superinfection [18,21,22,41].However, we argue that this functional form, as embedded by [35] in a model for P. vivax, is not appropriate in a model that takes hypnozoite accrual into account (see Appendix C for details).The multiscale model of [50] is intended to serve as a re-formulation of the infinite-compartment model proposed in [35].To incorporate hypnozoite accrual in a simple population-level framework, it draws on the relapse rate, conditional on the absence of blood-stage infection, derived in Mehra et al. [52], to derive a numerically tractable system of integrodifferential equations (IDEs).However, while the above-mentioned conditional relapse rate is derived under a framework that explicitly allows for superinfection, the population-level model of [50] does not have separate compartments for different values of the MOB (see Appendix D for details, and a proposed correction to [50] that has been adopted in subsequent work [54]).
To the best of our knowledge, the thesis of Mehra [51] -which forms of the basis of this paper -derives the first model that characterises both superinfection and hypnozoite dynamics for P. vivax.The model assumptions, building on the work of White et al. [35], are detailed in Section 2, where we re-visit the queueing model introduced in Mehra et al. [52] to characterise within-host superinfection and hypnozoite dynamics as a function of the FORI.We then construct a density-dependent Markov population process to couple host and vector dynamics in a homogeneously-mixing population, whilst allowing for superinfection and the accrual of the hypnozoite reservoir in the absence of human demographics in Section 3, proving that disease extinction occurs almost surely (Theorem 3.1).Using the work of Barbour and Luczak [29], we obtain a functional law of large numbers (FLLN), that takes the form of an infinite-dimensional system of ordinary differential equations (ODEs) (Section 3.2) for which the model of Bailey [2] arises as a special case (Appendix E).By drawing on our previous analysis of the within-host model, we show that FLLN can be reduced to a single IDE governing the time evolution of the FORI; the dynamics of superinfection and the hypnozoite reservoir in the human population can be recovered as a function of the FORI solving this IDE (Section 4.2).We then establish a threshold phenomenon for the reduced IDE, with the disease-free equilibrium shown to be uniformly asymptotically stable if R 0 < 1, and an endemic equilibrium solution emerging iff R 0 > 1 (Theorem 4.3, Section 4.3).
A recurrent theme in our analysis of population-level models is the utility of a physical understanding of the within-host model, governing the hypnozoite/MOB burden within a single human as a function of the intensity of mosquito-to-human transmission.In Section 5, we discuss how the ideas presented in this manuscript constitute a general strategy that can be employed to construct tractable population-level models of malarial superinfection, with appropriate constraints on the underlying model structure.
2 Modelling within-host hypnozoite and superinfection dynamics using an open network of infinite server queues with batch arrivals White et al. [35] construct a within-host model for short-latency (tropical) strains of P. vivax, predicated on the following set of assumptions: • Each infective mosquito bite immediately gives rise to a primary (blood-stage) infection and establishes a batch of hypnozoites in the liver.
• Hypnozoite batch sizes N i (for the i th bite) which are independent and identically-distributed (i.i.d.) across bites, are geometrically-distributed with probability mass function for n ∈ Z ≥0 and mean E[N i ] = ν.
• Each hypnozoite in the liver undergoes activation at constant rate α, which gives rise to the host suffering a (blood-stage) relapse; and death at constant rate µ.
• Hypnozoites behave independently; that is, the dynamics of each of the hypnozoites is described by an independent stochastic process.
We extend the framework of White et al. [35] to explicitly allow for superinfection.Specifically, we make the assumption that each blood-stage infection (relapse or primary) is naturally cleared at constant rate γ, with independent dynamics for each hypnozoite and blood-stage infection; as such, the existence of a previous blood-stage infection does not preclude further blood-and liver-stage infections, nor alter the rate of clearance for subsequent blood-stage infections.
We first examine the within-host dynamics of superinfection and the hypnozoite reservoir as a function of the FORI.In Mehra et al. [52], we construct an open network of infinite server queues with batch arrivals to concurrently describe hypnozoite accrual and the burden of blood-stage infection (allowing for superinfection) as a function of mosquito-to-human transmission intensity.The formulation of [52] can be extended to account for long-latency phenotypes, which are predominantly found in temperate regions [35,44], and the administration of drug treatment at a pre-determined sequence of times.Here, we consider the simplest case of the model presented in [52], restricting our attention to short-latency phenotypes (characteristic of tropical transmission settings) [35] in the absence of drug treatment.
We begin by delineating the set of possible states that a single hypnozoite can occupy: • H indicates a hypnozoite that is currently present in the liver; • A indicates a hypnozoite that has activated to give rise to a relapse that is currently in progress; • C indicates a hypnozoite that has previously given rise to a relapse, which has since been cleared; • D indicates a hypnozoite that has died, rather than activating.
Hypnozoites in the liver (state H) undergo activation at rate α and death at rate µ.Relapses (state A) are cleared from the bloodstream at rate γ.This model, introduced in [52], is a simple extension of the short-latency model of [35] that accounts for the clearance of bloodstage infection.A continuous-time Markov chain model that captures the above dynamics has non-zero transition rates q(H, A) = α q(H, D) = µ q(A, C) = γ with absorbing states C and D.
We denote by p s (t) the probability that a hypnozoite is in state s ∈ {H, A, C, D} := S h at time t after inoculation.It is straightforward to establish that when α as in Equations ( 13) to ( 16) of [52].In practice, we would expect the mean duration of hypnozoite carriage 1/(α + µ) to exceed the mean duration of each blood-stage infection 1/γ.
Likewise, we delineate the state space for each primary infection: • P indicates a primary infection that is currently in progress • P C indicates a primary infection that has been cleared from the bloodstream.
We assume that each bite necessarily triggers a primary infection that is cleared naturally from the bloodstream at the constant rate γ, yielding a continuous-time Markov chain model with transition rates q(P, P C) = γ q(P C, P ) = 0.
To embed our model for a single hypnozoite in an epidemiological framework, we construct an open network of infinite server queues, labelled H, A, C, D, P, P C (Figure 1).The arrival process, comprising mosquito bites, is governed by non-homogeneous Poisson process with a time-dependent rate λ(t), such that Each bite leads to the arrival of a single 'individual' in queue P (that is, a primary infection), in addition to a geometrically-distributed batch (with PMF (1)) in queue H (representing the hypnozoite reservoir).Hypnozoites in queue H follow the dynamics described above in moving to queues A, C or D. A primary infection in queue P moves to queue P C at rate γ.
Denote by N s (t) the number of 'individuals' (either hypnozoites or infections) in queue s ∈ {H, A, C, D, P, P C} := S at time t.From Equation (39) of [52], given is guaranteed to converge in the domain z ∈ [0, 1] 6 .
Explicit formulae for quantities of epidemiological interest -including marginal distributions for MOB and the hypnozoite burden; the proportion of recurrences that are expected to be relapses and the cumulative burden of blood-stage infection over time -are provided in [52].
There the marginal distributions for MOB and the hypnozoite burden are expressed in terms of partial exponential Bell polynomials.We can equivalently formulate recurrence relations to compute marginal distributions of interest following the approach of Willmot and Drekic [16], or compute joint distributions using the recurrence relations we derive in Mehra and Taylor [55].

A Markov population process with countably many types
To couple vector and human dynamics, whilst accounting for superinfection and the hypnozoite reservoir in the absence of human demographics (that is, without allowing for births/deaths in the human population), we construct a density-dependent Markov population process with countably many types.We consider a closed, homogeneously mixing population, with fixed numbers of humans P H and mosquitoes P M , where: • Each mosquito bites humans at constant rate β.
• When an uninfected mosquito bites a blood-stage infected human (that is, a human with at least one ongoing relapse or primary infection), human-to-mosquito transmission occurs with probability q.
• When an infected mosquito bites any human in the population, mosquito-to-human transmission occurs with probability p.This involves the establishment of: a primary (blood-stage) infection, which increases the MOB of the human host by one; and a batch of hypnozoites in the liver, where the batch sizes are geometrically-distributed with PMF (1) and are i.i.d.across bites.
• Within a human, each hypnozoite and primary infection is governed by an independent stochastic processes.
• Within a human, each hypnozoite undergoes activation at rate α and death at rate µ.
• Within a human, each primary infection or relapse (triggered by a hypnozoite activation event) is cleared independently at constant rate γ.
• Each infected mosquito is replaced by an uninfected mosquito at rate g.
We can formulate the state space of the Markov chain in several ways.One approach is to label each human with a hypnozoite/MOB state, and each mosquito with a binary infection state.
Specifically, we can denote the state at time t as (u(t), v(t)) where • u(t) is a P H -dimensional vector whose r th component itself is a vector giving the number of 'individuals' in compartment H, and the sum of individuals in compartments A and P in Figure 1 for the r th human; and • v(t) is a is a P M -dimensional vector whose s th component is 1 if the s th mosquito is infected and zero otherwise, yielding the state space Alternatively, we can formulate the state space to count the number of humans and mosquitoes in each respective infection state.At time t, define • h i,j (t), i, j ∈ Z ≥0 to be the number of humans with a hypnozoite reservoir of size i and MOB precisely j, and to be the number of uninfected and infected mosquitoes respectively.
We can thus denote the state of the Markov chain at time t to be (h(t), m(t)) where The state space is then where the 1 2 (i + j + 1)(i + j) + i + 1 th term of h(t) corresponds to h i,j (t).For notational convenience, we denote by e i,j the 1 2 (i + j + 1)(i + j) To obtain a density-dependent Markov population process, we proceed with the state description (h(t), m(t)).The transition rates are which can be understood as follows: • q (h,m),(h−e i,j +e i,j−1 ,m) : a human with hypnozoite reservoir size i and MOB j clears a single brood.
• q (h,m),(h−e i,j +e i−1,j ,m) : a single hypnozoite dies in a human with hypnozoite reservoir size i and MOB j.
• q (h,m),(h−e i,j +e i−1,j+1 ,m) : a single hypnozoite activates, thereby giving rise to a relapse (that is, a blood-stage infection with one additional brood), in a human with hypnozoite reservoir size i and MOB j.
• q (h,m),(h−e i,j +e i+k,j+1 ,m) : a human with hypnozoite reservoir size i and MOB j gains an additional batch of k hypnozoites, in addition to a primary infection (that is, a blood-stage infection with one additional brood) through the bite of an infected mosquito.
• q (h,m),(h,m−e 0 +e 1 ) : an uninfected mosquito is infected by taking a bloodmeal from an infected human.
We set the ratio of the human and mosquito population sizes to be v := P M P H , allowing us to write the transition rates given by Equations ( 7) to (12) in the form where ω J h ,Jm : R → R are given by ω (e i−1,j −e i,j ,0) (H, M) = µiH i,j , i ≥ 1, j ≥ 0 ( 14) 16) As such, we have defined a density-dependent Markov process with size parameter P M .

Steady state behaviour
The disease-free state is necessarily absorbing.However, since we have a density-dependent Markov population process with countably many types, absorption in the disease-free state is not guaranteed a priori since the mean MOB and/or hypnozoite reservoir size in the population can drift to infinity [15].
In Theorem 3.1 below, we show that the process does, in fact, reach the disease free state with probability one.To do this, we adopt the state description (u(t), v(t)) whereby each human is labelled with a hypnzoite/MOB state, and each mosquito is assigned a binary infection state.
If the arrivals of blood stage infections and batches of hypnozoites to the humans and infections to the mosquitoes were independent Poisson processes, then the the process (u(t), v(t)) could be regarded as a modelling P H independent batch arrival infinite server queueing processes with a structure as in Figure 1 and P M independent on-off processes whose s th component tells us whether the sth mosquito is infected or not.
However this is not quite the case.An 'arrival' of an infection to a human depends on the number of infected mosquitoes and conversely, an 'arrival' of an infection to a mosquito depends on the number of infected humans.
To overcome this, we couple the actual process (u(t), v(t)) with a second process (u in which humans (mosquitoes) become infected at constant rates independent of the number of mosquitoes (humans) that are infected.Specifically, in the process (u ′ (t), v ′ (t)), in modelling the rate at which humans become infected, we assume that all the mosquitoes are infected all of the time and, in modelling the rate at which mosquitoes become infected, we assume that all the humans are infected all of the time.
This amount of infection in the process (u ′ (t), v ′ (t)) can be shown to dominate the amount of infection in (u(t), v(t)).Furthermore, we can recognise (u ′ (t), v ′ (t)) as an independent network of infinite server queues for which a stability criterion is known, and hence establish that this criterion must apply to (u(t), v(t)) as well.
The details are given below.
Theorem 3.1.With probability one, disease is eventually eliminated, that is, and the time to disease elimination has finite expectation.
Proof.To show that absorption in the disease-free state (u(t), v(t)) = (0, 0) occurs almost surely, with a finite expected time to absorption, we adopt a coupling argument.
Consider an ensemble of P H independent networks of infinite server queues, as defined in Section Under this setting, we necessarily have In particular, and as such, the hitting time yields an upper bound for the time to absorption in the disease-free state under the process Using Corollary 4.1.1 of Mehra and Taylor [55], since the expected network occupation time for each hypnozoite/primary infection is finite and the hypnozoite batch size has finite mean, each component {u ′ r (t) : t ≥ 0}, r = 1, . . ., P H is ergodic.Further, since each Markov chain {v ′ s (t) : t ≥ 0}, s = 1, . . ., P M is irreducible and possesses a finite state space, it is also ergodic.
Each component {u ′ r (t) : t ≥ 0} and {v ′ s (t) : t ≥ 0} is positive recurrent and thus possesses a stationary distribution.The components evolve independently, therefore the stationary distribution for the multidimensional product {(u ′ (t), v ′ (t)) : t ≥ 0} should be given by the product of sta- it follows that it must be positive recurrent.This establishes that each state (u 0 , v 0 ) ∈ χ ′ is positive recurrent and the return time to the disease-free state has finite expectation E[T (0,0) ] < ∞.
Any state (u 0 , v 0 ) ∈ χ ′ can be reached from the disease-free state (0, 0) with positive probability p(u 0 , v 0 ) > 0 prior to return to the disease-free state through a concerted series of mosquito-inoculation events, with appropriate constraints on the time course of each infection and mosquito lifetimes.Consequently, that is, the expected hitting time for the disease-free state under the process {(u ′ (t), v ′ (t)) : t ≥ 0} has finite mean, irrespective of the initial condition (u 0 , v 0 ) ∈ χ ′ .
Since the time to absorption in the disease free state under the process {(u(t), v(t)) : t ≥ 0} with intial condition (u 0 , v 0 ) ∈ χ ′ is bounded above by T (u 0 ,v 0 ) , it immediately follows that and that the time to disease elimination, moreover, has finite expectation.■

A functional law of large numbers
We now use the work of Barbour and Luczak [29] to obtain a FLLN for the density-dependent Markov population process.In Theorem 3.2 below, we show that the sample paths of the continuous time Markov chain in converge to the solution of in the limit P M → ∞.Here, we define the η-norm and the set To show that the the semilinear problem (19) with initial condition (H(0), M(0)) ∈ R η has a unique mild solution in the interval [0, ∞), we use the results of Barbour and Luczak [29] (which draw on Theorem 1.4, Chapter 6 of Pazy [30]), bounding the η-norm of (H, M) using the expected occupancy of nodes H, A and P of the network introduced in Section 2 (Appendix A.2).
To show the convergence of sample paths to the semilinear problem (19), we use Theorem 4.7 of Barbour and Luczak [29] after verifying a series of semigroup and transition rate assumptions (Appendices A.1 and A.3).
Theorem 3.2.The semilinear problem (19) with initial condition (H(0), M(0)) ∈ R η has a unique mild solution in the interval [0, ∞).For each T > 0, ∃K T , K T , K T such that if T log P M P M for P M large enough, then T log P M P M .

A deterministic population-level model
We devote the present section to the analysis of the semilinear system given by Equation (19), which arises as the FLLN limit [29] for the density-dependent Markov population process detailed in Section 3.

A countable system of ODEs
It is instructive to write out the components of the vector equation (19), rescaled to consider the proportion of humans in each hypnozoite/MOB state.At time t, define • H i,j (t) to be the proportion of humans with hypnozoite reservoir size i and MOB j • I M (t) to be the proportion of infected mosquitoes.
Then we obtain the the infinite-dimensional system of ODEs where, for notational convenience, we have set The case ν = 0 yields the classical model of superinfection for P. falciparum, as formulated by Bailey [2] (Appendix E).
Equation (20)  Rather than solving the system of ODEs (20) directly for some time dependent function I M (t), we can draw on a physical understanding of the queueing network to derive a generating function for the time dependent state probabilites H i,j (t) as a function of I M (τ ), τ ∈ [0, t).In particular, from Equation ( 6), given the absence of liver-and blood-stage infection in the human population at time zero, that is, where Equation ( 22) is guaranteed to converge in the domain (x, y) ∈ [0, 1] 2 for any fixed t.
We can exploit the functional dependence (22) between the generating function for the state probabilities H i,j (t) of the human population at time t, and the proportion of infected mosquitoes I M (τ ) for τ ∈ [0, t) to aid analysis of the infinite-compartment model given by Equations ( 20) and ( 21).

A reduced integrodifferential equation
Upon examination of the deterministic model structure, we note that the coupling between human and mosquito dynamics is completely determined by the time evolution of • the prevalence of blood-stage infection in the human population (that is, humans with MOB at least one), which can be written The utility of our observation in Section 4.1 emerges in the characterisation of the functional dependence between I H (t) and the FORI aI M (τ ) for τ ∈ [0, t).Specifically, given H 0,0 (0) = 1, H i,j (0) = 0 for all i + j > 0, (23) we can readily write using the generating function (22).
Substituting the integral equation ( 24) into the ODE (21) yields the IDE When the human population is initially blood-and liver-stage infection-naive (that is, initial condition (23) holds), the time evolution of the FORI aI M (t) is equivalent under the IDE (25), and the infinite-dimensional system of ODEs given by Equations ( 20) and (21).
As a function of the FORI aI M (τ ), τ ∈ [0, t) which solves the IDE (25), we can use the generating function (22) to recover the population-level distribution of hypnozoite and superinfection states in the human population at time t.The recurrence relations provided in Mehra and Taylor [55] can be used to obtain the proportion of humans H i,j (t) in each hypnozoite/MOB state.The time-evolution of various quantities of biological interest can be recovered using the formulae provided in Mehra et al. [52].

The steady state distribution
Steady state analysis for the countably infinite system of ODEs given by Equations ( 20) and ( 21) is not straightforward.Instead, we characterise the existence and asymptomatic stability of equilibria for the IDE (25).
Upon inspection, we find that an equilibrium solution I * M of the IDE (25) must satisfy the integral equation The disease-free equilibrium I * M = 0 clearly satisfies Equation (26).On the domain I * M ∈ [0, 1], we note that the LHS of ( 26) is a monotonically increasing, convex function of I * M that approaches infinity in the limit I * M → 1 while the RHS is monotonically increasing, concave.The existence of a second solution to Equation ( 26) therefore depends on the limiting slope as I M → 0. It exists if and only if the derivative of the LHS is less than the derivative of the RHS at I * M = 0 (see Figure 2).This allows us to identify a bifurcation parameter with R 0 > 1 a necessary and sufficient condition for the existence of an endemic equilibrium We can readily interpret the integral expression in Equation (27).The probability of a current blood-stage infection time τ after an infective bite is given by the integrand no ongoing primary infection therefore yields the expected cumulative duration of blood-stage infection attributable to a single infective bite.The functional form of the bifurcation parameter R 0 is analogous to that of a reproduction number [10].
Using the stability criterion of Brauer [5], whereby the IDE ( 25) is linearised about each equilibrium solution, we can also establish sufficient conditions for equilibrium solutions I * M to be uniformly asymptotically stable.The threshold behaviour of the IDE (25) is summarised in Theorem 4.3 below, with a proof provided in Appendix B. 1.If R 0 < 1, then the disease-free equilibrium I * M = 0 is uniformly asymptotically stable and no endemic equilibrium solution exists.
2. If R 0 > 1, there exists precisely one endemic equilibrium I * M > 0, given by the non-trivial solution to the equation This endemic equilibrium I * M > 0 is uniformly asymptotically stable if

A general strategy for constructing tractable models of malarial superinfection
Our conceptual approach for modelling superinfection and hypnozoite dynamics for vivax malaria has been three-fold.
To establish a functional relationship between the intensity of mosquito-to-human transmission, and the dynamics of superinfection and hypnozoite accrual on the within-host scale (Section 2), we construct an open network of infinite server queues that accounts for stochasticity in: • the temporal sequence of mosquito-to-human transmission events, modelled with a nonhomogeneous Poisson process; • parasite inocula and consequently hypnozoite batch sizes, assumed to be geometricallydistributed for each bite [35]; and • the time course of each hypnozoite and primary infection, governed by independent stochastic processes.
The independence structure of the queueing network enables us to characterise the occupancy distribution through relatively straightforward physical arguments [52].
We model population level dynamics by constructing a Markov population process to addressing the dependence between the burden of blood-stage infection in the human population, and the intensity of mosquito-to-human transmission (Section 3).To characterise the steady state behaviour, we couple the Markov population process to an independent ensemble of queueing networks, with the same structure as the within host model and a known stability criterion (Theorem 3.1).Re-formulation of the state space to count the number of humans and mosquitoes in each permissible infection state, followed by a rescaling argument, yields a density-dependent Markov population process for which a FLLN limit, taking the form of a countably-infinite system of ODEs, is recovered using the work of Barbour and Luczak [29] (Theorem 3.2).
We pay particular attention to the FLLN limit, recognising that it takes a form identical to what could be derived using a deterministic compartmental model of the type that are ubiquitous in mathematical epidemiology (Section 4).By recognising that those parts of the compartmental model that relate to the infection status of the human population are identical to the Kolmogorov forward differential equations for a network of infinite-server queues, we propose, to the best of our knowledge, a novel reduction to collapse the countable system of ODEs into a reduced IDE that is more amenable to analysis, and for which a threshold phenomenon can be characterised (Theorem 4.3).
A common vein in the analysis of population-level models through the course of this manuscript is the utility of a physical understanding of the underlying within-host model.We now seek to illustrate how the ideas presented in this manuscript can be adopted more generally to construct tractable models of malarial superinfection, with appropriate constraints on the model structure.

From a countable system of ODEs to a reduced integrodifferential equation
We arrived at the countable system of ODEs given by Equations ( 20) and ( 21) by recovering a FLLN for an appropriately-scaled density dependent Markov population process.Following a standard compartment modelling approach, however, we can view this system as a natural extension of the Ross-Macdonald framework to allow for superinfection [2] and hypnozoite accrual [35].
We can also derive Equations ( 20) and ( 21) under the "hybrid approximation", as delineated by Nåsell [33] and Henry [43].Equation (20), which captures the time evolution of the within-host PMF, can also be interpreted as the governing equation for the expected frequency distribution in the human population, as a function of the FORI.Similarly, Equation ( 21) is precisely the Kolmogorov forward differential equation equation governing the time evolution of the probability that each mosquito in the population is infected -which can likewise be interpreted as the expected proportion of infected mosquitoes -as a function of the prevalence of blood-stage infection in the human population.The premise of the "hybrid approximation" is to allow for stochasticity within the human and mosquito populations respectively, but to couple host and vector dynamics through expected values [33,43]; the hybrid approximation to the present model of hypnozoite accrual and superinfection yields precisely Equations ( 20) and (21).Parallels between the FLLN and hybrid approximation have been noted previously in the literature: according to Hadeler and Dietz [7], Rost [6] proves that hybrid and mean-field approximations are equivalent1 , while Lewis [4] establishes similar results through simulation.
Irrespective of the derivation of the infinite-compartment model it is the observation that Equa-tion (20) constitutes the Kolmogorov forward differential equations for an open network of infinite server queues (Section 4.1) that constitutes the crux of our analysis.Rather than solving the system of ODEs (20) directly, we use a physical understanding of the queueing structure to derive the PGF (22).From the PGF, we extract the prevalence of blood-stage infection as a function of the FORI -which underpins the coupling between host and vector dynamics, and enables the derivation of a single IDE governing the FORI; in addition to quantities of biological interest, the time evolution of which can be recovered as a function of the FORI solving the aforementioned IDE.
Challenges posed by infinite-dimensional systems of ODEs for simulation and analysis have prompted the development of methods like the pseudoequilibrium approximation to yield finitedimensional systems of ODEs in the presence of superinfection [3,43].Our IDE reduction offers an alternative construction that yields an identical FORI to the original system of ODEs, under the assumption that the human population is initially uninfected.
Relative to the countable system of ODEs which served as our starting point, the reduced IDE is more amenable to steady state analysis.For compartment models with finitely many types (comprising finite-dimensional systems of ODEs), the next generation matrix method of Diekmann, Heesterbeek, and Roberts [26] is routinely used to characterise threshold phenomena: the basic reproduction number R 0 -which quantifies the number of additional cases stemming from an index case -is computed as the spectral radius of a "next generation operator", with the key implication that R 0 < 1 is a necessary and sufficient condition for the stability of the trivial (disease-free) equilibrium [20].While the notion of a reproduction number and spectral bound generalises to epidemic models with countably many types [25], computation is not straightforward.In the present manuscript, we apply the stability criterion of Brauer [5] to perform an asymptotic stability analysis for the reduced IDE, but are unaware of readilyverifiable stability criteria applicable to the countable system of ODEs.

Analyticity at the within-host scale
The strategy delineated above is underpinned by the fact that the within-host model is analyticallytractable.The critical assumption is that consequences of each mosquito-to-human transmission event are governed by independent stochastic processes, for which a time-dependent PGF can be derived.We can accommodate multiple 'types' of incoming parasites/infections (here, we distinguish hypnozoites from immediately developing parasite forms, represented by a primary infection), and occupancy times for each parasite/lifecycle stage do not necessarily have to be Markovian (that is, we can allow for general service time distributions).Tractable extensions to the within-host model that permit the IDE reduction proposed here encompass the class of models analysed in Mehra and Taylor [55], comprising open networks of infinite server queues with nonhomogeneous multivariate batch Poisson arrivals.They are applicable when the evolution of 'individuals' (hypnozoites or blood-stage infections in this case) occurs independently.
Importantly, this is not the case when there is competition or another form of interaction between parasite broods.Further, while service time distributions and arrival rates can vary deterministically with time, they should not be dependent on the within-host state, barring the consideration of certain forms of immunity.
Nonetheless, our work provides a theoretical basis to explore the epidemiology of malarial superinfection, with natural extensions to allow for greater biological realism.A generalisation of the present framework, allowing for long-latency hypnozoites [35,44]; time-varying transmission parameters, and the acquisition of transmission-blocking and antidisease immunity for P. vivax, is presented in Mehra et al. [53].There, we derive an infinite compartment model under a "hybrid approximation" [33,43], comprising an infinite-dimensional system of ODEs, from which we derive a reduced system of IDEs using the approach detailed in Section 4. Analyses of both steady state solutions and transient dynamics under this generalised deterministic model are tractable using the within-host distributions derived in Mehra et al. [52].

A A functional law of large numbers
We now verify the conditions for the FLLN Barbour and Luczak [29] of the Markov population process governing hypnozoite and superinfection dynamics, as described in Section 3. The scaled transition rates ω J are given by Equations ( 13) to (18).We define in addition to the subset = 0 for all but finitely many i, j}.
Let J be the set of jumps with non-zero transition rates.Jumps are of "bounded influence" in the sense of [29], in that at most two compartments are affected.Further, Therefore, Assumptions 1.2 and 1.3 of Barbour and Luczak [29] are satisfied, with the implication that Markov process with transition rates given by Equations ( 13) to ( 18) is a "pure jump process, at least for some non-zero length of time" [29].
In the notation of Barbour and Luczak [29], we write the proposed limit as a semilinear equation where A is a constant (that is, not dependent on time) matrix, with non-zero terms while F has components A.1 Assumption 2.1 of Barbour and Luczak [29] We recapitulate and verify Assumption 2.1 of Barbour and Luczak [29] below, with changes in notation/exposition in line with the current manuscript.
We begin by introducing the set By assigning a 'size' to each human and mosquito infection state through a fixed vector η = we can construct an empirical moment analogue To bound these empirical moments, Barbour and Luczak [29] introduce the quantites .
In a similar vein, we set Assumption 2.1 of Barbour and Luczak [29] then reads as follows.
Assumption 1.There exists η ∈ R satisfying condition (29) and r (1) max ≥ 1 such that the We begin by noting that In the case r = 1, we can simplify these expressions substantially to yield the bounds For r ≥ 1 and (H, M) ∈ R with |H| 1 = 1/v, |M| 1 = 1, we obtain the bounds where H b is a geometrically-distributed random variable with mean ν and state space Z ≥0 , noting that E[H p b ] < ∞ for all p ≥ 0.
Likewise, for a given value (H, M) ∈ X + , for all r > 0 since the RHS constitutes a finite sum.We thus see that Assumption 1 is satisfied for arbitrarily large r max ≥ 1. ■ A.2 Existence of a unique weak solution to the semilinear equation (19) We now establish the existence and uniqueness of weak solutions to the semilinear equation (19) (which can be written explicitly as the system of ODEs given by Equations ( 20) and ( 21) after re-scaling H i,j by a factor of v).
We are then required to fix µ ∈ R with µ ≥ 1, with the constraint that: 1. there exists r ≤ r max such that sup i,j≥0 {µ 2. A T µ ≤ wµ for some w ≥ 0 (Assumption 3.2 of [29]) Since Assumption 1 is satisfied for arbitrarily large r max , it suffices to set µ = η, in which case we can choose w = max{α + µ, γ, g}.We thus define the η-norm The conditions of Theorem 3.1 of Barbour and Luczak [29] are thus satisfied.The implications of this theorem are as follows.As in Barbour and Luczak [29], denote by P (•) "the semigroup of Markov transition matrices corresponding to the minimal process associated with Q", where Then R is a strongly continuous semigroup on R η , with A = R ′ (0) and R ′ (t) = R(t)A for all t ≥ 0.
(22), we see that for all t ≥ 0 Hence, A.
≥ 1 such that max and some b = b(ζ) ≥ 1, and that Proof.For our choice of ζ ∈ R, we have where H b is a geometrically-distributed random variable with mean ν and state space Z ≥0 . Further, ) is monotonically decreasing on the interval [0, M 0 ) and monotonically increasing on the interval (M 0 , 1).
Therefore, an endemic equilibrium solution exists iff and given R 0 > 1, this endemic equilibrium is unique.
We now characterise the asymptotic stability of the disease-free and endemic equilibria.By Theorem 2 of Brauer [5], an equilibrium solution I M (t) = I * M is uniformly asymptotically stable if the trivial solution u(t) = 0 to the first-order linear IDE is uniformly asymptotically stable.By Theorem 1 of Brauer [5], since P (τ ) is continuous and non-negative for all τ ∈ [0, ∞), it follows that u(t) = 0 is uniformly asymptotically stable if and Equation ( 33) serves as a sufficient condition for a steady state solution I * M to be uniformly asymptotically-stable.Noting that any steady state solution I * M must solve Equation ( 32), we derive the equivalent condition Using Equation (34), it immediately follows that the disease-free equilibrium I * M = 0 is uniformly asymptotically stable in the case where R 0 < 1.
When R 0 > 1, we have h(0) > 0 and h(1) < 0, so there exists precisely one root h(x) = 0 in the domain x ∈ (0, 1) and serves as a sufficient condition for the endemic equilibrium to be uniformly asymptotically stable.

■
C The pseudeoquilibrium approximation to the superinfection recovery rate proposed by White et al. [35] White et al. [35] construct a deterministic compartmental model, consisting of an infinitedimensional system of ODEs monitoring the hypnozoite burden and the absence or presence of blood-stage infection.Superinfection is captured through a "pseudoequilibrium approximation" for the rate of recovery from blood-stage infection, conditional on the size of the hypnozoite reservoir.Below, we characterise this construction in detail.
The pseudoequilibrium approximation to the rate of recovery from (super)infection was initially proposed by Dietz, Molineaux, and Thomas [3].For an M/M/∞ queue with homogeneous arrival rate λ and service rate γ -which describes within-host superinfection dynamics for P.
falciparum in a constant transmission setting -Dietz, Molineaux, and Thomas [3] note that the expected duration of each busy period, or 'superinfection' at steady state is Therefore, as a function of the FORI λ(t) = λ at time t, Dietz, Molineaux, and Thomas [3] argue that the rate of recovery from superinfection at time t can be approximated to be This argument essentially replaces a rate of recovery, which necessarily depends on the MOB, with its expectation taken under steady-state conditions with the arrival rate taken to be constant.While this might be a reasonable thing to do if there is no better option, as we show in this paper, the model is still tractable when the MOB is included in the state space.It might be that the pseudoequilibrium approximation yields reasonable results when changes in the FORI occur on a slower time scale than the clearance of blood-stage infection; it has, however, been critiqued on the basis of its inability to capture MOB dynamics following abrupt changes in transmission intensity [24].
To account for superinfection, White et al. [35] suggest that the recovery rate from blood-stage infection, conditional on the hypnozoite burden i should be + 1 (36) where λ represents the FORI.For an individual with a hypnozoite reservoir of size i, additional blood-stage infections are acquired at rate λ + iα -with primary infections occurring at rate λ due to mosquito inoculation, and each hypnozoite in the liver activating independently at rate α to give rise to a relapse.Equation ( 36) is a modification of Equation ( 35) accounting for the additional contribution of hypnozoite activation to the blood-stage infection burden.However, like the approximation (35) it still arises from replacing a recovery rate that depends on the current value of the MOB with its expectation under steady-state conditions, in which case the MOB would have a Poisson distribution with parameter (λ + iα)/γ.
If a single hypnozoite were established through each infective bite, then the distribution of MOB conditional on the size of the hypnozoite reservoir would be Poisson-distributed; for a network of infinite server queues with single arrivals, the numbers of busy servers in each queue follow independent Poisson distributions at each time t [14].The batch structure in hypnozoites established through each infective bite, however, breaks this construction; indeed, Boucherie and Taylor [12] show that the only open networks of queues that have transient product form distributions are networks of M t /M/∞ queues.
The distribution of MOB, conditional on the size of the hypnozoite reservoir, can be recovered using the within-host queueing model introduced in Mehra et al. [52]  where H denotes the generating function given by Equation (22).Following similar reasoning to Nåsell [33], given a hypnozoite reservoir of size i, the rate of recovery from (super)infection can be written Explicit formulae for the rate of recovery from superinfection, in terms of partial Bell polynomials, can be recovered using Leibniz's integral rule and Faa di Bruno's formula [27]; the resultant expressions, however, are unwieldy and are therefore omitted here.A simpler frameworkformulated in terms of a recovery rate corrected for superinfection -is detailed in Appendix D.  [50].
However, rather than formulating a model in terms of the relapse rate and the rate of recovery from (super)infection, a simpler but equivalent approach is to directly consider the (integral) expression for the expected prevalence of blood-stage infection in the human population over time under the within-host framework of [52], as we have done in this paper.
While the proposed correction to the framework of [50] is somewhat subtle, it has important implications.As discussed in [50], the omission of superinfection at the population-level leads to discrepancies relative to the framework of [35], the re-formulation of which was the primary aim; these discrepancies are evident at high transmission intensities, where superinfections would be expected.The deterministic model introduced in Section 4, in contrast, serves as an extension of the framework of [35] to allow for superinfection.The omission of superinfection at the population-level in [50] means that a number of key epidemiological parameters that are analytically tractable at the within-host level do not make sense at the population-level.A key example is the expected contribution of relapse vs primary infection to the overall burden of recurrence (which includes the case of overlapping relapses and primary infections, as in [52]),

Figure 1 :
Figure 1: Schematic of the open network of infinite server queues governing the within-host hypnozoite and MOB burden, as a function of the intensity of mosquito-to-human transmission.Adapted from Figure 3 of Mehra et al. [52] as a special case (short-latency hypnozoites, no drug treatment).

Case 1 :Figure 2 :
Figure 2: Schematic of the geometric argument to establish the existence of a non-zero solution to Equation (26).
H ) − γI H , which appeared in