An epidemic model for an evolving pathogen with strain-dependent immunity
AAn epidemic model for an evolving pathogen withstrain-dependant immunity
Adam Griffin a,b, ∗ , Gareth O. Roberts a , Simon E.F. Spencer a a Department of Statistics, University of Warwick, Coventry, UK b Centre for Ecology & Hydrology, Benson Lane, Wallingford, Oxfordshire, OX10 8BB, UK
Abstract
Between pandemics, the influenza virus exhibits periods of incremental evolutionvia a process known as antigenic drift. This process gives rise to a sequence ofstrains of the pathogen that are continuously replaced by newer strains, prevent-ing a build up of immunity in the host population. In this paper, a parsimoniousepidemic model is defined that attempts to capture the dynamics of evolvingstrains within a host population. The ‘evolving strains’ epidemic model hasmany properties that lie in-between the Susceptible-Infected-Susceptible andthe Susceptible-Infected-Removed epidemic models, due to the fact that indi-viduals can only be infected by each strain once, but remain susceptible toreinfection by newly emerged strains. Coupling results are used to identify keyproperties, such as the time to extinction. A range of reproduction numbers areexplored to characterise the model, including a novel quasi-stationary reproduc-tion number that can be used to describe the re-emergence of the pathogen intoa population with ‘average’ levels of strain immunity, analogous to the begin-ning of the winter peak in influenza. Finally the quasi-stationary distributionof the evolving strains model is explored via simulation.
Keywords:
Epidemiology, Probabilistic models, Quasistationary distributions ∗ Corresponding author
Email addresses: [email protected] (Adam Griffin), [email protected] (Gareth O. Roberts), [email protected] (Simon E.F. Spencer)
Preprint submitted to Mathematical Biosciences August 18, 2020 a r X i v : . [ q - b i o . P E ] A ug . Introduction Epidemic models have become important tools for understanding, predict-ing and developing mitigation strategies for public health planners dealing withinfectious diseases. Recent advances in genetic epidemiology have greatly accel-erated our understanding of the complex interactions between host immunityand pathogen evolution, and emphasised the important role that pathogen evo-lution can have on the dynamics of infection. However, it remains extremelychallenging to combine together these two interacting processes within the samemathematical framework [1]. In this paper we develop a parsimonious epidemicmodel that describes the transmission dynamics of a multi-strain pathogen withevolutionary dynamics similar to the influenza A virus evolving via antigenicdrift.Multi-strain models have become increasingly popular due to the rise inavailability of pathogen genetic analyses. Many models have been based onordinary differential equations (ODE), despite the fact that stochastic effectsplay an important role in mutation [see 1, for a review]. Bichara et al. [2] developan epidemic model with competition between finitely many pathogen strains,and include vertical transmission and immunity from maternal antibodies inthe infection dynamics. Meehan et al. [3] analyse multi-strain epidemic modelswith mutation between strains within an ODE framework. However since theirfocus is on drug-resistance, they do not consider the effect of immunity. In themulti-strain models discussed in Gog et al. [4], there is assumed to be a finitenumber of possible strains, and each individual may be infected with one or moreof such strains. Evolution was been modelled by a random jump process on afinite strain space using a nearest neighbour jump process. Models involving acountable number of infectious statuses have been discussed in the past [5], butthese typically only use the previously mentioned nearest-neighbour evolution.In [5] this is expressed as a model for parasitic infections where the “type” ofan individual is defined by the quantity of parasites in a host. Despite themany modelling papers on multi-strain epidemics, the methodology required to2t these models to data is only just emerging [6].Between pandemics, the 4 main sub-types of the influenza virus evolve ac-cording to a process called antigenic drift [7]. Antigenic drift arises due to thefact that infection with a particular strain of influenza provides the host with along-lasting immunity to future infection by the same strain. Once immunity toa particular strain has built up in the population, there is a selection advantageto strains that do not elicit the same immune response. To capture within amathematical model the complex processes driving the evolution of the influenzavirus is extremely challenging due to the interactions between host immunityand viral evolution [8]. Nonetheless, simple models can give rise to surprisinglycomplex dynamics [9, 10]. The H3N2 subtype of influenza A, in particular, ex-hibits a narrow spread in its evolutionary tree, with all strains a short geneticdistance from a single branch [11, 12]. Each strain persists for a relatively shortamount of time before being replaced.In this paper we define a novel epidemic model with countably infinite,evolving strains that sits between the traditional susceptible-infected-susceptible(SIS) and susceptible-infected-removed (SIR) epidemic models, in that each in-dividual may be infected many times with the pathogen, but only once by astrain. The model is designed to reflect the linear pattern of evolution observedin pathogens undergoing antigenic drift, such as seasonal influenza. By intro-ducing an equivalence relation on the state space, we are able to describe theequilibrium behaviour of the model prior to elimination of the pathogen. InSection 5, coupling arguments are used to make precise the relationship be-tween our new model and the traditional SIS and SIR models and to explorethe large population limit. In Section 6 we discuss three reproduction num-bers for the novel model. Finally in Section 7 we explore simulations from thequasi-equilibrium distributions. 3 . Definition of the model
Consider a closed population of N individuals which are classified as suscep-tible , infective or removed . For each time t , we denote the number of susceptiblesby S ( t ), the number of infectives by I ( t ), and the number of removed individ-uals by R ( t ). An infective remains in this class for a random period of timeknown as their infectious period, after which they become removed. Similarly,removed individuals become susceptible again after their immune period, duringwhich they cannot be infected by any strain (even a new one). We assume thedurations of infectious periods are i.i.d. draws from L I ∼ Exp( γ ) and immuneperiods are i.i.d draws from L R ∼ Exp( δ ).To capture dynamics of competing and evolving strains, every individualhas a strain index k ∈ Z which denotes the most recent strain with which anindividual was (or is currently) infected. We denote the number of susceptibles,infectives and removed individuals respectively with strain index k by S k ( t ), I k ( t ) and R k ( t ). Finally we denote by K ∗ ( t ) the largest strain index observedup to time t and use K ∗ := K ∗ ( t − ) where the time is clear from context.As in the standard SIRS model [13], we assume homogeneous mixing ofindividuals, and so each pair of individuals makes contact at the points of aPoisson process with rate βN >
0. New strains are introduced into the populationin the following way. Each time an infective makes contact with a susceptibleindividual, we assume that with some probability θ ∈ [0 , K ∗ ( t − ) + 1. With probability 1 − θ , the original strain in the infectiveattempts to infect the susceptible; the success of this infection depends on thestrain index of the susceptible. For simplicity, we assume that immunity iscumulative: a susceptible with strain index k is immune to all strains withindex j ≤ k . Removed and susceptible individuals retain the strain index of thestrain they have most recently recovered from. All contact processes, mutationevents, infectious periods and immune periods are assumed to be independent4rom each other.To summarise, the epidemic proceeds according to the following events. • Infection without mutation: ( S j ( t ) , I k ( t )) (cid:55)→ ( S j ( t ) − , I k ( t ) + 1)for all j < k , with rate β (1 − θ ) N − S j ( t ) I k ( t ). • Infection with mutation: ( S j ( t ) , I K ∗ +1 ( t ) = 0) (cid:55)→ ( S j ( t ) − , I K ∗ +1 ( t ) = 1) for j ∈ Z with rate βθN − S j ( t ) I ( t ). • Recovery: ( I k ( t ) , R k ( t )) (cid:55)→ ( I k ( t ) − , R k ( t ) + 1) for k ∈ Z with rate γ I k ( t ). • Loss of global immunity: ( R k ( t ) , S k ( t )) (cid:55)→ ( R k ( t ) − , S k ( t ) + 1) for k ∈ Z with rate δ R k ( t ).The state space of the E-SIRS model is given by Ω (cid:48) = { ( s , i , r ) : (cid:80) k ∈ Z (s k +i k + r k ) = N } , with s , i , r being infinite sequences taking values in { , , . . . , N } .A natural initial condition might be I (0) = 1 , I k (0) = 0 for k (cid:54) = 1, S (0) = N − S k (0) = 0 for k (cid:54) = 0, and R k (0) = 0 for all k . This equates to a singlecurrently infective individual infected with a strain to which all other individualsare susceptible, and to which no-one is currently recovering. Consider a second model where following an infectious period, an individualbecomes immediately susceptible, corresponding to the E-SIRS model where δ = ∞ . In this model there are no periods of immunity and so recovery eventsgenerate susceptibles: • Recovery: ( I k ( t ) , S k ( t )) (cid:55)→ ( I k ( t ) − , S k ( t ) + 1) for k ∈ Z with rate γ I k ( t ).5lus infection transitions as above. The E-SIS model evolves over the sub-space { ( s , i ) : (cid:80) k ∈ Z (s k + i k ) = N } . We will refer to both spaces by Ω (cid:48) , themeaning will always be clear from context. Consider the E-SIS model with θ = 1. All contacts are mutation contactsand hence successful, and so ( S ( t ) , I ( t )), the total numbers of susceptibles andinfectives, follow a traditional single-strain SIS model as defined in [14]. Wecan also perform a similar identification between ( S ( t ) , I ( t ) , R ( t )), the numberof susceptibles, infectives and immune individuals in the E-SIRS model and thesingle-strain SIRS model as defined in [13].On the other hand, consider the E-SIS model with θ = 0. Since no contactsare mutations, no individual can be infected more than once. If the popula-tion starts with strain index 0 except for the initial infectives with strain 1,( S ( t ) , I ( t ) , S ( t )) behaves as a traditional single-strain SIR model as definedin [14].
3. Equivalence relation
We wish to study the long-term average behaviour of characteristics suchas the levels of immunity and pathogen diversity, however the constant emer-gence and extinction of strains means that the evolving epidemic process hasno steady-state. To counter this we introduce an equivalence relation to fix theprocess against the most recently emerged strain.
Definition 1.
The active strain set of a state ( s , i , r ) ∈ Ω (cid:48) is given by K = { k ∈ Z : i k > } . Let elements of this set be indexed from 1 to | K | in ascendingorder, so for k a , k b ∈ K , we have k a < k b whenever a < b . Definition 2.
Two states ( s , i , r ) and ( s (cid:48) , i (cid:48) , r (cid:48) ) ∈ Ω (cid:48) are equivalent if and onlyif the following conditions hold. The total numbers of susceptibles, infectives and removed individuals areequal: | s | = (cid:80) k ∈ Z s k = (cid:80) k ∈ Z s (cid:48) k = | s (cid:48) | , and similarly | i | = | i (cid:48) | and | r | = | r (cid:48) | . The numbers of active strains are equal: | K | = | K (cid:48) | . Each active strain has the same number of infectives: i k a = i (cid:48) k (cid:48) a , for a =1 , . . . , | K | . The numbers of individuals that are susceptible to the a th active strain areequal: (cid:80) k Definition 3. The representative of the equivalence class containing ( s , i , r ) ,denoted ( s ∗ , i ∗ , r ∗ ) is defined as follows. If K (cid:54) = ∅ , denote the active strains forthe representative by K ∗ = { − | K | , . . . , } . Let φ : K → K ∗ be a bijectiondefined by φ ( k a ) = a − | K | for a = 1 . . . , | K | . Then the representative ( s ∗ , i ∗ , r ∗ ) is given by: ∗ k = i φ − ( k ) for k ∈ K ∗ , otherwise. s ∗ k = (cid:80) φ − ( k +1) − j = φ − ( k ) s j for k ∈ { − | K | , . . . , − } , (cid:80) ∞ j = φ − (0) s j for k = 0 , (cid:80) φ − (1 −| K | ) − j = −∞ s j for k = −| K | , otherwise. r ∗ k = (cid:80) φ − ( k +1) − j = φ − ( k ) r j for k ∈ { − | K | , . . . , − } , (cid:80) ∞ j = φ − (0) r j for k = 0 , (cid:80) φ − (1 −| K | ) − j = −∞ r j for k = −| K | , otherwise. If K = ∅ then i ∗ k = 0 for all k ∈ Z and s ∗ = (cid:80) j ∈ Z s j and s ∗ k = 0 for k (cid:54) = 0 , andsimilarly for r . Let { } denote the set of all these absorbing states. In the rest of this paper, the process of representatives on the space ofequivalence classes (states described with starred states as in Definition 3) willbe referred to as the normalised process, and will be denoted by ( S ∗ , I ∗ , R ∗ ) or( S ∗ , I ∗ ) as appropriate. Definitions 2 and 3 remove all strains with no infectiveindividuals, and give index 0 to the most recent strain to have emerged and haveinfectives. All susceptibles and removed individuals are given the strain indexone less than the nearest infective above them in strain order. Any individualsimmune to all existing strains are given strain 0, as though they just recoveredfrom the most recently emerged strain. Example 4. Consider the state ( s , i , r ) ∈ Ω (cid:48) given by (s , . . . , s ) = (0 , , , , , , , . . . , i ) = (0 , , , , , , , . . . , r ) = (1 , , , , , , here all remaining terms of s , i , r are zero. Strains 2 and 6 are active so K = { , } ⇒ K ∗ = {− , } and the representative under the equivalence relation isgiven by (s ∗− , . . . , s ∗ ) = (0 , , , (i ∗− , . . . , i ∗ ) = (0 , , , (r ∗− , . . . , r ∗ ) = (1 , , . Notation. Recall that without the equivalence relation, the state space of theepidemic process was Ω (cid:48) . We denote the state space of the normalised processover the set of equivalence class representatives by Ω.We will use x = ( s , i , r ) ∈ Ω with, for example, s = (s k ) k ∈ { , . . . , N } Z todenote a typical element of the state space. We will also use, for example, | s | = (cid:80) k ∈ Z s k to denote the total number of susceptibles. A random variable writtenin calligraphic type, e.g. R k ( t ), refers to a process without the equivalencerelation. The corresponding variable written in roman type, e.g. R k ( t ), refersto the normalised process evolving over representatives from the equivalenceclasses. 4. Quasi-stationarity and absorbing states Like many infectious disease models, the E-SIRS and E-SIS models defined inSection 2 have an absorbing class of states that corresponds to the populationcontaining no infected individuals, i = . For finite population models, theabsorbing state is reached with certainty in finite time, and so the limitingdistribution is degenerate with no mass in non-absorbing states. However, likethe single strain SIS model, these processes may not go extinct for a long time(individuals can be reinfected) and the transient quasi-stable behaviour is ofinterest. The quasi-stationary distribution and limiting conditional distributionconditioned on the epidemic not going extinct, represent the long-term averagebehaviour for an endemic disease. 9 .1. Properties of quasi-stationary distributions for epidemics In this section and the rest of this paper, P u [ A ] = P [ A | X (0) ∼ u ]. Definition 5. Let X = ( X ( t )) t ≥ be a Markov process on a countable statespace Ω with absorbing state from which it cannot escape. Then a distribution u on Ω \ { } is a quasi-stationary distribution (QSD) if P u [ X ( t ) ∈ A | X ( t ) (cid:54) =0] = u ( A ) for all t ≥ .Given initial condition v on S = Ω \ { } , u is a v - limiting conditionaldistribution (LCD) if lim t →∞ P v [ X ( t ) ∈ A | X ( t ) (cid:54) = 0] = u ( A ) . Note that, forprocesses where S is a single communicating class, every QSD u is a u -LCDand every LCD is a QSD. Related to the QSD on irreducible state spaces is the notion of the decayparameter which describes the rate of decay of the transition probabilities. Definition 6. Let X = ( X ( t )) t ≥ be an irreducible Markov process on a count-able state space Ω with absorbing state . Let u be a QSD associated to X .Then the decay parameter α is given by α = inf { a ≥ P ij ( t ) = o ( e − at ) } . for i, j ∈ Ω \ { } The absorption parameter α is given by α = inf { a ≥ (cid:90) ∞ P i [ T > t ] e at dt = ∞} , for i ∈ Ω \ { } where T is the extinction time of X starting from state i . Notethat for irreducible processes, α is independent of i, j and α is independent of i . According to Theorem 6 of [15], a necessary condition for the existence of aQSD is that α > Theorem 7. Conditional on non-absorption, the following hold. The QSD for the number of infectives in the single-strain SIS model existsuniquely and gives weight to all states { , . . . , N } . The QSD for the number of susceptibles and infectives in the single-strainSIR model exists uniquely and gives full weight to the state { ( S, I ) =(0 , } . The QSD for the number of susceptibles and infectives in the single-strainSIRS model exists uniquely, and gives weight to all non-absorbing states.Proof. Theorem 1 of [16] states that QSDs exist and are unique on finite irre-ducible state spaces, and so there is a unique QSD for the SIS model and forthe SIRS model conditional on { I > } , and non-zero weight is given to all non-absorbing states. For reducible processes, Theorem 8 of [16] states that QSDswill give full weight to the communicating class with the longest expect time toleave and any states accessible from this “slowest” communicating class. Thischaracterises the QSD for the SIR model.Further work on characterising the QSD for the standard SIS model can beseen in [17, 18, 19] making use of recurrent processes and normal approximations. Here we will summarise the existence and uniqueness results for the E-SISand E-SIRS processes. Theorem 8. Let the E-SIS model be defined as in Subsection 2.2 with param-eters β, γ > and θ ∈ (0 , . Then for the normalised process, conditionalon the events { I ( t ) > } , there exists a unique QSD which equals the uniqueLCD of the process and gives weight to all non-absorbing (i.e. transient) states, { ( s , i ) ∈ Ω : | i | > } . If θ = 0 and the process begins with a single infective,then there exists a unique LCD which gives full weight to the state with a singleinfective with strain index , and N − susceptibles with strain index 0.Proof. For θ ∈ (0 , S =Ω \{ } is a single finite communicating class, which immediately gives existence,uniqueness and equality of the QSD and LCD by Theorem 3 from [15].Under the equivalence relation, there can be at most N different strain in-dices present in the population. This implies that every individual must appear11n one of the states s − N , . . . s or i − N , . . . , i . Therefore we can bound abovethe size of S , the set of transient (i.e. non-absorbing) states, by (2 N ) N .One can see that the transient states form a single communicating class bynoting that one can get from a single infective of strain index 0 with N − θ = 0 we consider the E-SIS model starting with a single infectiveof strain 0 and susceptibles of strain index − 1. If θ = 0, then mutation isimpossible. As a result, once an individual has become infected and recovered,they join the s class and cannot be reinfected. In this way { S ( t ) } behavesidentically to the { R ( t ) } class in the SIR model, and we identify the two modelsin this way. Point 2 in Theorem 7 then gives the required LCD. Theorem 9. Let the E-SIRS model be defined as in Subsection 2.1 with pa-rameters β, γ, δ > and θ ∈ (0 , . Then, conditional on having at least oneinfective, there exists a unique QSD. If θ = 0 and there is initially one infective,a QSD still exists and gives full weight to the state with one infective with strainindex and N − susceptibles with strain index .Proof. For θ ∈ (0 , N ) N , since individuals may alsoreside in classes r − N , . . . , r . The fact that the transient states form a singlecommunicating class also follows as in Theorem 8. For θ = 0, we see that eachstate that can be reached is a transient communicating class; there is no way toreturn to a state once left. As such, we need to consider the decay parameteron leaving each state, which equals the exponential rate of leaving such a state: β s − i /N + γ i + δ r . The decay parameter for the process is therefore theminimal such value across all non-absorbing states. Therefore s − = 0, i = 1and r = 0 minimise the decay parameter. According to Theorem 1 of [16] thisforces the QSD to have full mass on this state where i = 1 , s = N − 1, sincethe only state accessible from this state is an absorbing one.12 .3. Sampling the quasi-stationary distribution Samples from quasi-stationary distributions can be produced using the se-quential Monte Carlo (SMC) sampler with stopping time resampling methodsdeveloped in [20]. In brief, M realisations of the model (referred to as particles)are simulated forward in time. Absorbed particles (with no infected individu-als) are given weight zero and non-absorbed particles are given weight 1 initially.The distribution of weights converges to the limiting conditional distribution.Once the total weight drops below a proportion λ of the initial weight, theparticles are replenished via a resampling step. Combine-split resampling [20]was used, which prevents any occupied states from being lost and has the ad-vantage that the distribution of weights remains unchanged after resampling.This resampling method combines particles in the same location together, drawsnew particle locations uniformly from the existing locations and equalizes theweight on particles in the same location. In our implementation, after a burn-inof T b = 1, weighted samples were drawn every T d = 1 time units until time T max = 100.Figure 1 shows the expected number of individuals in each class under theQSD. In this example we used M = 1000 particles and a resampling thresholdof λ = 0 . 6. The Figure shows that when δ is smaller there is a larger proportionof globally immune individuals in the removed classes, and so the populationcan support fewer strains. 5. Limiting behaviour One aspect of interest is how the E-SIS and E-SIRS processes relate to thosewithout mutation. To this end, we consider the limits of the times to extinctionof the processes as θ tends to 0 or 1, and the limit, for fixed θ , of the time toextinction as the population size N tends to infinity. Large population limits canbe used to justify the use of infinite population models as approximations for,for example, the decay parameters for the relevant processes which we cannotobtain analytically. 13 24 −20 −16 −12 −8 −4 0RemovedInfectivesSusceptiblesStrain Index E x pe c t ed N u m be r o f I nd i v i dua l s −24 −20 −16 −12 −8 −4 0RemovedInfectivesSusceptiblesStrain Index E x pe c t ed N u m be r o f I nd i v i dua l s . . . . . . . Figure 1: Comparisons of expected population composition under E-SIRS QSDs with (a) β = 2, θ = 0 . δ = 0 . N = 25; (b) β = 2, θ = 0 . δ = 2, N = 25. Theorem 10. Let T θ be the time to extinction of the E-SIS model, and T the time to extinction of the standard SIS model, both starting from a singleinfective (nominally of strain index 1) in a population of N individuals. Then T θ → T in distribution as θ → .Proof. We make use of a coupling of ( T θ : 0 < θ < 1) and T . Firstly, let X θ ( t ) = ( S θ ( t ) , I θ ( t )) be the E-SIS model. We assume the process to be definedover a population indexed by n = 1 , . . . , N . • For each individual n , define a sequence of i.i.d. infectious periods { L ( n ) m ∼ Exp( γ ) : m ∈ N } . • For each ordered pair of individuals ( n, n (cid:48) ), define a homogeneous Poissonprocess { A ( n,n (cid:48) ) ( t ) } on [0 , ∞ ) with rate β/N . • For each ordered pair ( n, n (cid:48) ) define a sequence of uniform random variables U ( n,n (cid:48) ) l ∼ Unif[0 , 1] for l ∈ N . • Let all L ( n ) m , A ( n,n (cid:48) ) and U ( n,n (cid:48) ) l be independent.14ow let ( E, F , P ) be the product space of these random processes and variables.Using these building blocks, we construct an E-SIS model { X θ ( t ) } and SISmodel { Y ( t ) = ( S ( t ) , I ( t )) } as follows. Set S θ − (0) = N − , I θ (0) = 1 for theE-SIS model and set S (0) = N − 1, and I (0) = 1 for Y (0). Assume the initialinfective individual has index n = 1 without loss of generality. In both processesinfectious individual n generates contacts with each susceptible individual n (cid:48) atpoints of the Poisson process { A ( n,n (cid:48) ) ( u ) } , where u denotes the cumulative timethat individual n has been infectious and n (cid:48) has been susceptible. In other wordsthe Poisson processes are stopped whenever it is not possible for individual n toinfect individual n (cid:48) . At each contact event in the SIS model an infection occurs.The newly infected individual n (cid:48) stays infected for a period of length L ( n (cid:48) ) m ( n (cid:48) ,t )+1 ,where m ( n (cid:48) , t ) is the number of infections individual n (cid:48) has recovered from upto the contact time t . In the E-SIS model the i th contact event in { A ( n,n (cid:48) ) ( u ) } results in a mutation if and only if U ( n,n (cid:48) ) i ≤ θ , in which case individual n (cid:48) isinfected with a new strain and given the lowest unused strain index. Howeverif U ( n,n (cid:48) ) i > θ then individual n attempts to infect n (cid:48) with their current strainand the infection is successful if the strain index of individual n (cid:48) is strictly lessthan the stain index of n . As in the SIS model, successful infections in the E-SIS model are given infectious period length L ( n (cid:48) ) m ( n (cid:48) ,t )+1 . Notice that under thiscoupling non-mutation contacts of n (cid:48) by n are only successful if the strain indexof n is strictly greater than that of n (cid:48) in the E-SIS model. However mutationcontacts and all contacts in the SIS model are always successful.Fix ω ∈ E , our probability space. For the SIS model, we have T < ∞ almost surely. On the interval [0 , T ( ω )), there are two possibilities for the E-SISmodel. At each infective-susceptible contact we compare the strain indices. Thefirst possibility is that every contact will always lead to a successful infection,arising from a sequence of infection events which always contact a susceptibleof a lower index. In this case, we have T θ ( ω ) = T ( ω ) for all θ ∈ [0 , U ( n,n (cid:48) ) i > θ . Since we must have a finite number of such events15ccurring in [0 , T ), we can find θ ∗ such that U ( n,n (cid:48) ) i ≤ θ ∗ for all such U ( n,n (cid:48) ) i corresponding to potentially unsuccessful events. This means that for θ ≥ θ ∗ we must have T θ ∗ ( ω ) = T ( ω ). So for every ω ∈ E , there exists θ ∗ ∈ (0 , T θ ( ω ) = T ( ω ) for all θ ≥ θ ∗ . Hence T θ ( ω ) → T ( ω ) as θ → ω ∈ E , and hence T θ → T in distribution by the SkorohodDudley theorem from, for example, Theorem 3 of [21]. Theorem 11. Let T be the time to extinction of the standard SIR model. Then T θ → T in distribution as θ → . Intuitively, one can think of identifying the S -class for the E-SIS modeland the R -class of the SIR model. As mutation events get rarer, the chanceof mutation happening before extinction becomes smaller and smaller, and sothe two processes are more likely to coincide under a suitable coupling untilextinction. Proof. This proof follows in a similar fashion to Theorem 10. In this version,the coupling is constructed between the E-SIS and the SIR model. The onlydifferences are that in the SIR model individuals enter the removed categoryafter their infectious period and the Poisson process { A ( n,n (cid:48) ) ( u ) } progressesduring any time for which n is infective and n (cid:48) is susceptible in the E-SIS model(as before); but when n (cid:48) is susceptible or removed in the SIR model.In the E-SIS model infectious contacts between n and n (cid:48) are only successfulif the event is a mutation or n (cid:48) is of a strictly lower strain index than n . In theSIR model, only the first infectious contact is successful. This means that thetwo epidemics must be identical up to the time of the first repeat contact, whenone identifies the {S , S , S , . . . } classes in the E-SIS model with the R classof the SIR model.Similar to the proof of Theorem 10, for each ω ∈ E one can find a value of θ ∗ so that T θ ∗ ( ω ) = T ( ω ) for all θ ≤ θ ∗ , and so T θ → T almost surely as θ → T θ → T in distribution by the Skorohod Dudley theorem of[21]. 16 .2. Large population limits In order to obtain some large population limit results, we will consider an“infinite” population model. We will refer to this as an Evolving Birth-DeathProcess (E-BDP). More precisely we assume that infected individuals are anegligible part of an infinite population of individuals that are not immune toany strains at the start of the epidemic, and so all infections will be successfulalmost surely. This implies infections from a given strain k and recoveries fromthat strain behave as a linear birth-death process with birth rate β and deathrate γ . Additionally, at the point of each infection, with probability θ ∈ [0 , K ∗ + 1.The possible events comprise: • Infection with mutation: I K ∗ +1 ( t ) = 0 (cid:55)→ I K ∗ +1 ( t ) = 1 with rate βθ I ( t ). • Infection without mutation: I k ( t ) (cid:55)→ I k ( t ) + 1 for k ∈ Z with rate β (1 − θ ) I k ( t ). • Recovery: I k ( t ) (cid:55)→ I k ( t ) − k ∈ Z with rate γ I k ( t ).After it emerges, each strain behaves according to a linear birth-death pro-cess with birth rate β (1 − θ ) and death rate γ . The total number of infectives I ( t ) also behaves according a birth-death process with birth rate β and deathrate γ .The time to extinction of the E-SIS model converges to that of the E-BDPmodel, noting that under a suitable coupling, the time to extinction of the E-BDP equals the Linear BDP without mutation. This leads us to the followingresult. Theorem 12. Let T θ,N be the time to extinction of the E-SIS model, and T thetime to extinction of the E-BDP model. Then we have T θ,N → T in distributionas N → ∞ when β < γ . If β ≥ γ , then on the event { T < ∞} , a region ofprobability − γ/β , we also have T θ,N → T in distribution as N → ∞ roof. Using Theorems 10 and 11 we can conclude that for any fixed N that T ,N is the time to extinction for the standard SIR model, and T ,N is equal tothe time to extinction for the standard SIS epidemic model. Furthermore, fromthese theorems we can construct a coupling of the SIS, E-SIS and SIR modelsusing two sets of Poisson processes and mutation indicator variables such that,for any θ ∈ [0 , T ,N ( ω ) ≤ T θ,N ( ω ) ≤ T ,N ( ω ) (1)for almost every ω ∈ E . From [22], we know that if β < γ then T ,N convergesin distribution to T , the time to extinction of a Linear BDP with the sameparameters β and γ . From [23] we obtain that the same thing happens for SISmodels, i.e. T ,N → T in distribution as N → ∞ . Using the bounds in Equation(1), we obtain that T θ,N → T as N → ∞ for all θ ∈ [0 , β ≥ γ we note that on a set of probability 1 − γ/β , the timeto extinction of the linear BDP is infinite, as discussed in Chapter 3.2 of [24].From [23], we know that T ,N → T almost surely (and hence in distribution)on the event { T < ∞} . From [22] we know that on this event, T ,N → T indistribution. Therefore we must have that T θ,N → T as N → ∞ for all θ ∈ [0 , { T = ∞} we don’t have T ,N →∞ . Instead T ,N converges to an extreme-value distribution as mentioned inTheorem 8.1 [14].Next we show existence of a QSD for the E-BDP model. Theorem 13. Let X ( t ) be the E-BDP with parameters γ > β > and θ ∈ [0 , .Then, under the equivalence relation described in Section 3 and conditional onthe event { I ( t ) > } , there exists a unique QSD.Proof. To prove existence, we first show that the state space we are interestedin is countable. To do this we use the following construction. Starting witha single infective of strain 0, we can define a method of constructing the statespace. By having a birth in strain 0, or a mutation event, one can systematically18rrive at any state in the state space. Given these two possible events, one canencode each state according to a finite binary sequence, which corresponds toa unique integer which we can use to enumerate the space. Given that thereexists a lower bound l ∈ Z such that I k = 0 for all k ≤ l , we construct the stateas follows. Starting with the lowest non-zero strain index l + 1 consider I l +1 strain 0 within-strain-infection events. Then for each higher strain k , we choosea mutation event followed by I k − X ( t ) =( X j ( t )) j ∈ Z be the E-BDP. Let α X be the decay parameter for X ( t ). Let( Y ( t )) t ≥ be the process defined on the same probability space, given by Y ( t ) = (cid:80) j ∈ Z X j ( t ). Since the mutations do not affect whether infections are successfulor not, Y ( t ) can be seen to be a single-strain linear BDP with birth rate β , anddeath rate γ . As discussed in Example 1 of [25], Y ( t ) has the decay parameter α Y = γ − β . Let T X be the extinction time of X ( t ) and T Y for Y ( t ).Letting α Y be the absorption parameter for Y ( t ), and α X for X ( t ), we alsoknow that α Y = γ − β . Since T X = T Y under the coupling, we use the definitionof the decay parameter to deduce that α X = γ − β , and hence α X ≥ α X > α X = γ − β since there is only one state from whichextinction can occur: one must have 1 infective before extinction, which mustbe of strain 0 under the equivalence relation. This leads to the uniqueness ofthe QSD. 6. Reproduction numbers To characterise the dynamics of the models, we look to a number of keystatistics which are related to the commonly used basic reproduction number , R ,that illustrates whether or not an epidemic is likely to infect a large proportionof the population. The basic reproduction number is defined as the number of19ndividuals infected by a single typical infective in a large, otherwise susceptiblepopulation [27]. In the E-SIRS model, we still have R = β/γ . One issuewith R is that it fails to take into account the likely immunities present in thepopulation, or how much the pathogen evolves during the opening phase of theepidemic. R ∗ In [28], an epidemic is considered which spreads through a population groupedinto households, such that individuals in the same household make contact at adifferent rate to individuals in different households. The households reproduc-tion number R ∗ is shown in [28, Section 2.3] to be equal to R ∗ = µR H where µ is the expected number of individuals infected in a single household epidemic(including the initial infective), and R H is the the mean number of contacts aninfective individual makes with individuals in other households during a singleinfectious period.For the E-SIRS model, we consider each strain as a “household” which hascountably many individuals, and mutations are considered contacts betweenhouseholds. In this case µ is the expected total population of a birth deathprocess with birth rate β (1 − θ ) and death rate γ , including the initial infective.One can use the branching property to compute E [ Z ] = γ/ ( γ − β (1 − θ )) andnote that the between household reproduction rate R H = βθ/γ and so, R ∗ = βθγ − β (1 − θ ) β (1 − θ ) < γ ∞ β (1 − θ ) ≥ γ To recontextualise this in terms of strains and mutations, one can think of R H as the expected number of new strains originating from a single individualduring one infectious period, and µ as the expected number of individuals thatever get infected by a specific strain.One could consider R to be the “intra-strain” reproduction number, and µ to be the “inter-strain” reproduction number. With these we obtain one ofthree regimes: 20 If R = β/γ < 1, then the whole population would die out with certainty,and no large epidemic would occur. • If R ≥ µ < ∞ then a large epidemic occurs with positive probabil-ity, but each individual strain dies out quickly. • If R ≥ µ = ∞ , then each strain has a positive probability ofproducing a large outbreak.Figure 2 shows realisations of the genetic trees under the E-SIS model under thetwo supercritical regimes. For small θ , we obtain only a small number of strains,and the epidemic is more likely to die out. Moreover, in a finite population, thislow θ leads to high immunity in the population and hence shorter epidemics.The trees highlight how only a small number of the strains survive for a longtime, particularly in Fig. 2(a). Figure 2 also show similarities to the tree forH3N2 in Extended Data Figure 9(c) of [12], a paper which specifically looks tomodel influenza. Time E m e r g i ng S t r a i n s R � H � ∗ � ∞ (a) Time E m e r g i ng S t r a i n s R � H � ∗ � (b) Figure 2: Comparisons of emergence of strains under different R , R H under the E-SIS modelwith N = 25, γ = 1, β = 1 . θ = 0 . 1; (b) θ = 0 . .2. Quasi-stationary reproduction number R Q One drawback to the R is that it only usefully describes the initial behaviourof an epidemic in a naive population and doesn’t take into account the build upof immunity in the E-SIRS model. One alternative is to consider the effectivereproduction number , denoted R t , defined as R t = R S ( t ) N in a population ofsize N . Much work has been done in trying to evaluate R t for specific infectionssuch as influenza by [29] and Ebola by [30]. However, R t is time-dependent andcan therefore be difficult to compute and interpret. At a quasi-stable equilib-rium the number of new infections balances the recoveries and so R t ≈ 1, andhence R t is not informative about the disease characteristics. Ideally, we wouldlike a reproduction number that adjusts for the build-up of immunity in thepopulation, but remains informative about the infectivity of a disease.We offer an alternative reproduction number, based on the QSD, which aimsto describe the infectiousness of strains of an endemic disease in a populationwith ‘average’ levels of historical immunity. The quasi-stationary reproductionnumber ( R Q ) is the average number of secondary infections caused by a singletypical infective introduced into an otherwise uninfected (S status) populationwith levels of immunity (strain indexes) drawn from the quasi-stationary dis-tribution, so each other individual may or may not be immune to the currentstrain of the infective. By typical infective, we mean an individual with strainindex sampled from the distribution of strain indexes of infectives in the QSD.Under the E-SIRS model, the total number of infectives is always less than theSIS model without evolving strains, and so R Q ≤ R .The quasi-stationary reproduction number provides a measure of the abilityof a pathogen to re-invade a population from which it has been eradicated. Fordiseases like seasonal influenza which have greatly reduced incidence during thesummer months, R Q measures the reproduction number at the beginning of thenext influenza season after accounting for the residual immunity left over fromlast year.More precisely, we draw the single infective from the marginal number ofinfectives in the QSD u I ( k ): the probability that given an individual is infective,22t is of strain index k . For QSD u this is given by u I ( k ) = (cid:88) ( s , i , r ) ∈ Ω u ( s , i , r ) i k | i | . Under the equivalence relation described in Section 3, we can have a maximumof N strains in a population of size N , and so the strain index ranges over k ∈ K ∗ = { , − N } . The susceptible population is drawn from the total strainmarginals of the QSD u K ( k ): the probability that under the QSD that a givenindividual is of strain index k . u K ( k ) = (cid:88) ( s , i , r ) ∈ Ω u ( s , i , r ) i k + r k + s k N Finally, we require the probability that a randomly chosen individual drawnfrom the strain marginal will be susceptible to strain k (i.e. will have a strainindex lower than k ): u L ( k ) = k − (cid:88) j =1 − N u K ( j ) . During their infectious period the infective makes infectious contact with eachindividual at the points of a Poisson process with rate β/N . For large popu-lations the infective is unlikely to contact the same individual twice (or them-selves), and so the expected number of contacts is β/γ . With probability θ thecontacts are mutations and are successful infections. With probability (1 − θ )the contacts are non-mutations and are only successful if the individual con-tacted has a lower strain index, which occurs with probability u L ( k ) when theinfective has strain index k . To calculate R Q we condition on the strain indexof the initial infective, hence R Q = βγ (cid:88) k =1 − N ( θ + (1 − θ ) u L ( k )) u I ( k )= βγ ( θ + (1 − θ ) u TL u I ) . (2)Since u I and u L are both probability mass functions, 0 ≤ u TL u I ≤ R H ≤ R Q ≤ R . As θ → R Q → β/γ = R , as does R H .23 .0 0.2 0.4 0.6 0.8 1.0 . . . . . . q R ep r odu c t i on N u m be r R R * R Q q R ep r odu c t i on N u m be r R R * R Q Figure 3: Comparisons of R , R ∗ and R Q under varying θ . with (a) β = 0 . γ = 1, N = 100;(b) β = 2, γ = 1, N = 100. The three notions of a reproduction number in this section describe threedifferent facets of the epidemic model, and can be compared in Figure 3. Itshows that R Q is always less than R due to the effects of immunity, and R ∗ depends greatly on θ ; the unplotted points for R ∗ are in the regions where itis infinite, namely where β (1 − θ ) ≥ γ . Values of R Q were calculated using theSMC sampler described in Section 4.3 with M = 100 particles and a resamplingthreshold of λ = 0 . 7. Simulation Study To further explore the E-SIRS model we use the SMC sampler describedin Section 4.3 to investigate numerically features of the QSD which we cannotobtain analytically. We wish to observe how various key properties behave as wevary parameters of the model. To this end we look at the following expectationsover the QSD. For brevity we omit time indices and conditioning, and denoteexpectations under the QSD by E Q . • The expected total number of infectives E Q [ I ] and immune individuals E Q [ R ] in the QSD, where I = (cid:80) k = −∞ I k , R = (cid:80) k = −∞ R k .24 The expected total number of active strains E Q [ K ] in the QSD where K = |{ k : I k > }| . • We also look at how varying the model parameters affects strain diversityin infectives and the whole population.We will focus on the E-SIS model, but also discuss for each statistic how theaddition of an immune period, as in the E-SIRS models, changes the numberof infectives and strain diversity. Unless otherwise stated, all expectations overthe QSD were produced with the SMC sampler described in 4.3 with M = 100particles and resampling threshold λ = 0 . Figure 4 shows a heatmap of the expected number of infectives in the pop-ulation under quasi-stationarity, E Q [ I ], and how this depends on the contactrate and the probability of mutation. Increasing the contact rate β or muta-tion probability θ increases the expected number of infectives. However, for afixed population size (in this case N = 100), the number of infectives increaseslinearly in β when E Q [ I ] is much smaller than N . This can be observed inFigure 5(a), which shows that the number of infectives grows more slowly as β increases, especially when θ is small and so the probability of failed infections ishigh. As was also noted in Section 4.3, increased levels of global immunity resultin fewer infectives under quasi-stationarity, due to the increased possibility offailed infections.Figure 5(b) shows that as N increases the expected proportion of infectives( E Q [ I ] /N ) decreases in the case where β < γ , whereas in the supercritical casewe see that E Q [ I ] /N remains fairly constant. In the E-SIRS model, we seethat E Q [ I ] /N is decreased by the introduction of transient global immunity.Furthermore, as δ gets smaller the transient immunity lasts longer and E Q [ I ] /N further decreases. 25 q b Figure 4: Expected number of infectives as β and θ change in E-SIS with γ = 1, N = 100. . . . . b E x pe c t ed P r opo r t i on o f I n f e c t i v e s SIS , q = , q = d = , q = d = , q = d = , q = d = , q = . . . . . . N E x pe c t ed P r opo r t i on o f I n f e c t i v e s SIS , b = , b = d = , b = d = , b = d = , b = d = , b = Figure 5: Expected proportion of infectives in the E-SIS and E-SIRS models: (a) as β varieswith γ = 1 , N = 100; (b) as N varies with γ = 1, θ = 0 . .2. Expected number of strains We investigated what happens to the expected number of active strains, E Q [ K ] (strains held by infectives) as the parameters change. Under our models,the number of strains is always less than the number of infectives due to theabsence of super-infectivity (infection of an individual by multiple strains duringa single infectious period). As such, much of the behaviour is similar to thatof the expected number of infectives in the previous subsection. For example,the expected number of strains increases linearly with β when E Q [ K ] is muchless than N . This follows since we already know that for θ = 1 every infectivebegins a new strain and so E Q [ I ] = E Q [ K ]. At the other end of the scale, weautomatically have that E Q [ K ] = 1 if θ = 0.Figure 6 shows the expected number of strains for fixed βθ (mutation con-tact rate) and β (1 − θ ) (non-mutation contact rate), as β and θ vary. Note thatfor Figure 6(b), both θ and β increase from left to right, whereas, to maintainfixed βθ , β decreases as θ increases. In the case when βθ is high, one mightexpect E Q [ I ] and E Q [ K ] to be closer in value since there is a high probabilityof mutation leading to a high number of co-circulating strains. This is demon-strated in Figure 6(a), where we see that for fixed βθ , the number of strainsis larger (and therefore closer to the number of infectives) for the βθ = 2 linethan for the βθ = 0 . 05 line. In Figure 6(a) there is a maximum point for thenumber of strains as β increases, after which the number of strains decreases.As β (1 − θ ) increases in Figure 6(b), the number of strains becomes more linearin θ . In Figure 7 we investigate the distribution of immunity across the activestrains. The figure shows the expected proportion of infectives E Q [ I/N ] andtotal individuals for each strain index E Q [ I k + S k + R k ], relative to the mostrecently emerged strain. The expectations taken over the quasi-stationary dis-tribution were calculated using the SMC sampler described in Section 4.3, with M = 100 particles, resampling threshold λ = 0 . T b = 11.27 .0 0.2 0.4 0.6 0.8 1.0 q E x pe c t ed N u m be r o f S t r a i n s l l l b q = b q = b q = b q = l b = b = b = b = q E x pe c t ed N u m be r o f S t r a i n s ll b ( - q ) = b ( - q ) = b ( - q ) = b ( - q ) = Figure 6: Expected number of strains under quasi-stationarity for (a) fixed βθ ; (b) fixed β (1 − θ ); with γ = 1 and N = 100. Figure 7(a) illustrates that larger values of β greatly increase the diversity ofstrains in the infectives and the variation in immunity in the population, sincethere are more infectives and so more chances for mutation contacts. Anotherpoint of interest is the lag of the strain diversity: the number of strains betweenthe mode of the infective strains and the mode of the total population. Thelag is fairly consistent for the different values of β , but does increase slowlyin β . Figure 7(b) shows the change in strain diversity in θ . As θ increases,the number of strains present increases, so the strain diversity curve flattensout. For high values of θ , a larger lag is observed between the infectives andthe whole population, due to the higher diversity. In Figure 7(c), the effectof population size is explored. As N increases, we observe a wider numberof strains, as one would expect given E Q [ K ]’s behaviour. However, unlike thebehaviour as β changes, the peak moves away from 0 but the lag between theinfectives and the rest of the population appears more consistent. For the E-SIRS model explored in Figure 7(d), the immune period reduces strain diversityby reducing the expected number of infectives.In applications, one might wish to look further into the lag between thestrain distribution of the infectives and the immunity in the population. If28 40 −30 −20 −10 0 . . . . . . . . Strain Index E x pe c t ed P r opo r t i on b = b = b = b = b = b = . . . . . Strain Index E x pe c t ed P r opo r t i on q = q = q = q = q = q = . . . . . Strain Index E x pe c t ed P r opo r t i on N = 50 − InfectivesN = 50 − Total IndivsN = 100 − InfectivesN = 100 − Total IndivsN = 200 − InfectivesN = 200 − Total Indivs −50 −40 −30 −20 −10 0 . . . . . Number of Infectives E x pe c t ed P r opo r t i on E−SIS − InfectivesE−SIS − Total IndivsE−SIRS, d = d = d = Figure 7: Strain diversity as a function of the parameters in the E-SIS and E-SIRS models:(a) β varies; (b) θ varies; (c) N varies; (d) δ varies. Unless otherwise stated β = 2, γ = 1, N = 100, θ = 0 . 29 pathogen has a long lag, then vaccination can be effective in updating theimmunity present in the population. However, if the lag is short, then a vaccinebased on a recent strain will have little effect in increasing the levels of immunityin the population, as the most immunity profiles in the population will alreadyrepresent the currently circulating pathogen. 8. Conclusions In this paper we defined an epidemic model for a pathogen undergoing ge-netic drift, that lies between the well-studied SIR and SIS epidemic models. Themodel appears to capture some qualitative aspects of the evolution of strainsof influenza A, despite depending on just 4 parameters. Compared to modelsused by [12] and [31], which require the storage of a whole antigenic history,our model is much simpler, which makes simulation, computation and inferencemuch easier. Despite these simplifications, the simulated genetic trees in Fig-ure 2 show similarities to the tree for H3N2 in Extended Data Figure 9(c) of[12]. The relative simplicity of our model enables analytical insights into modelbehaviour, such as the relationship between our models and the SIS and SIRmodels discussed in Theorems 11 and 10. A simulation study showed that thereis a nonlinear tradeoff between mutation and infectivity when trying to estimatethe number of co-circulating strains under quasistationarity. The developmentof a quasistationary reproduction number R Q also allows summary the expectedbehaviour of an epidemic under quasistationarity, by comparing it to a house-hold epidemic model. Clearly this work could be reduced to a finite state spaceof strains, but also include more complex strain evolution models, accountingfor similarity of strains conferring some amount of partial immunity. Acknowledgements Funding: AG was supported by EPSRC Grant Number EP/HO23364/1;GOR and SEFS were supported by EPSRC grant EP/R018561/1; SEFS wassupported by MRC grant MR/P026400/1.30 eferences [1] A. J. Kucharski, V. Andreasen, J. R. Gog, Capturing the dynamics ofpathogens with many strains, Journal of Mathematical Biology 72 (1-2)(2015) 1–24.[2] D. Bichara, A. Iggidr, G. Sallet, Global analysis of multi-strains SIS, SIRand MSIR epidemic models, Journal of Applied Mathematics and Comput-ing 44 (1-2) (2013) 273–292.[3] M. T. Meehan, D. G. Cocks, J. M. Trauer, E. S. McBryde, Coupled, multi-strain epidemic models of mutating pathogens, Mathematical Biosciences296 (2018) 82–92.[4] J. Gog, B. Grenfell, Dynamics and selection of many-strain pathogens,Proceedings of the National Academy of Sciences 99 (26) (2002) 17209–17214.[5] S. Moy, Extensions of a limit theorem of Everett, Ulam and Harris onmultitype branching processes to a branching process with countably manytypes, The Annals of Mathematical Statistics 38 (4) (1967) 992–999.[6] P. Touloupou, B. Finkenst¨adt, N. French, S. Spencer, Bayesian inferencefor multi-strain epidemics with application to