Extinction and quasi-stationarity for discrete-time, endemic SIS and SIR models
EExtinction and quasi-stationarity fordiscrete-time, endemic SIS and SIR models
Sebastian Schreiber ∗ , Shuo Huang † , Jifa Jiang † , and Hao Wang ‡ Abstract.
Stochastic discrete-time SIS and SIR models of endemic diseases are introduced and analyzed. For thedeterministic, mean-field model, the basic reproductive number R determines their global dynamics.If R ≤
1, then the frequency of infected individuals asymptotically converges to zero. If R > R values. To understand the length of the transient prior to extinctionas well as the behavior of the transients, the quasi-stationary distributions and the associated meantime to extinction are analyzed using large deviation methods. When R >
1, these mean timesto extinction are shown to increase exponentially with the population size N . Moreover, as N approaches ∞ , the quasi-stationary distributions are supported by a compact set bounded awayfrom extinction; sufficient conditions for convergence to a Dirac measure at the endemic equilibriumof the deterministic model are also given. In contrast, when R <
1, the mean times to extinctionare bounded above 1 / (1 − α ) where α < N approaches ∞ , the quasi-stationary distributions converge to a Dirac measure at thedisease-free equilibrium for the deterministic model. For several special cases, explicit formulas forapproximating the quasi-stationary distribution and the associated mean extinction are given. Theseformulas illustrate how for arbitrarily small R values, the mean time to extinction can be arbitrarilylarge, and how for arbitrarily large R values, the mean time to extinction can be arbitrarily large. Key words. infectious diseases, discrete-time SIS model, discrete-time SIR model, times to extinction, quasi-stationary distributions, large deviations
1. Introduction.
Infectious disease modeling has been one of the most important topicsin mathematical biology (Keeling & Rohani, 2011). A recent Google scholar search revealsover a million studies referencing SIS (susceptible-infected-susceptible) and SIR (susceptible-infected-recovered) models. Most of these studies use deterministic, continuous-time equa-tions. However, in seasonal systems or systems where measurements are taking at regulartime intervals (e.g. day or week), discrete-time models play an important role (Anderson& May, 1991; Allen, 1994; Allen & Burgin, 2000; Klepac et al. , 2009; Keeling & Rohani,2011). For many of these models, the basic reproductive number, R , determines whethera disease can persist or not in a population. If R >
1, persistence often occurs. While if R <
1, the disease-free state (i.e. extinction) often is globally stable and the infection is lostasymptotically as time marches to infinity.When considering finite populations without external sources of infection, Markov chain ∗ Department of Evolution and Ecology and Center for Population Biology, University of California, Davis, Cali-fornia 95616, United States ([email protected]) † Mathematics and Science College, Shanghai Normal University, Shanghai 200234, China([email protected], [email protected]) ‡ Department of Mathematical and Statistical Sciences, University of Alberta, Edmonton, Alberta T6G 2G1,Canada ([email protected]) On September 8th, 2019, a search on Google scholar with the search term “disease model SIR OR SIS”returned 1 .
45 million results. a r X i v : . [ q - b i o . P E ] M a y SCHREIBER, HUANG, JIANG, AND HAO models typically predict that the disease goes extinct in finite-time whether R > <
1. Tounderstand this puzzling difference between the asymptotic behaviors of the deterministic andstochastic models (Bartlett, 1966; Keeling & Rohani, 2011; Diekmann et al. , 2012), one canuse the concept of quasi-stationarity that describes the long-term behavior of the stochasticmodel conditioned on non-extinction (Darroch & Seneta, 1965, 1967). For finite state models,the quasi-stationary distribution corresponds to a normalized left eigenvector π of the tran-sition matrix of the Markov chain restricted to the transient states, i.e. the states where theinfection persists. In discrete-time, if the disease dynamics follow the quasi-stationary distri-bution, then the eigenvalue λ associated with this eigenvector corresponds to the probabilityof disease persistence over the next time step (respectively, a time interval of length one).Thus, when the stochastic model follows the quasi-stationary distribution, the time to extinc-tion is exponentially distributed with a mean time of extinction 1 / (1 − λ ). Grimm & Wissel(2004) call 1 / (1 − λ ) the intrinsic mean time to extinction. To understand the link betweenthe stochastic and deterministic models, it is natural to ask: How does the intrinsic meantime to extinction increase as the population size gets larger? How is the quasi-stationarydistribution related to the asymptotic dynamics of the deterministic model as the populationsize gets larger? More generally, how do these quantities depend on the parameters such as R ?For continuous-time, stochastic SIS models, there exists a dichotomy in the mean time toextinction when a fixed, positive fraction of the population is infected (Weiss & Dishon, 1971;Barbour, 1975; Kryscio & Lef`evre, 1989). When R >
1, this time increases exponentiallywith the population size N in the limit of large population sizes. When R <
1, these extinc-tion times remain bounded in the limit of large population size. However, to the best of ourknowledge, similar statements for the intrinsic mean extinction times have not been rigorouslyproven for these continuous-time SIR models. However, in a series of papers (N˚asell, 1996,1999, 2001, 2002), N˚asell provided methods to approximate the intrinsic mean extinction timesas well as the quasi-stationary distributions. His approximations support the existence of asimilar dichotomy for intrinsic mean times to extinction. Moreover, they highlight a remark-able dichotomy about the qualitative behavior of the quasi-stationary distributions. When R >
1, these distributions are well-approximated by a normal distribution centered near theendemic equilibrium. When R <
1, the quasi-stationary distribution is best approximatedby a discrete, geometric distribution. Despite these advances, mathematically rigorous resultsfor discrete-time, stochastic SIS and SIR models are lacking.In this paper, we introduce a new class of discrete-time SIS and SIR deterministic andstochastic models that have several desirable properties including (i) they are derived withindividual-based rules and, consequently, preserve non-negativity of all populations, (ii) thedeterministic models are the mean field model of the stochastic models, and (iii) the deter-ministic models converge to the classical continuous-time models in an appropriate limit. Forthese models, we analyze the global dynamics of deterministic models and then use this analy-sis in conjunction with large deviation results from Faure & Schreiber (2014) to rigorouslycharacterize the behavior of the intrinsic mean times to extinction and quasi-stationary dis-tributions in the limit of large population sizes. Moreover, for some special cases, we deriveexplicit approximations for the quasi-stationary distributions and extinction times that applyfor all population sizes.
N STOCHASTIC DISCRETE-TIME SIS AND SIR MODELS 3
Our paper is structured as follows. Section 2 introduces and analyzes the discrete-time,deterministic SIS model. Section 3 presents mathematical and numerical findings on quasi-stationary distributions and intrinsic mean extinction times for the stochastic SIS model.Section 4 introduces the discrete-time SIR model and proves results about its global attractors.Section 5 presents mathematical and numerical findings on quasi-stationary distributions andintrinsic mean extinction times for the stochastic SIR model. Section 6 discusses our mainfindings and future challenges. The mathematical proofs are given in Sections 7 through 10.
2. The dynamics of a deterministic SIS model.
We begin with a discrete-time versionof the classical susceptible-infected-susceptible (SIS) model where individuals are either sus-ceptible or infected. Let I n denote the fraction of individuals that are infected at time step n ,in which case, the fraction of susceptible individuals equals 1 − I n . Individuals escape natu-ral mortality with probability e − µ while infected individuals escape recovery with probability e − γ , where µ > γ >
0. Susceptible individuals from the previous time step who havenot died, escape infection with probability e − βI n where β > I n +1 = F ( I n ) := e − µ (1 − I n ) (cid:0) − e − βI n (cid:1) + e − µ − γ I n , ≤ I n ≤ . This discrete-time formulation of the SIS model has several advantages. First, it is straight-forward to verify that the dynamics of I n remain in the interval [0 ,
1] provided the initial value I lies in this interval. Second, these dynamics, as described in the next section, correspondto the mean field dynamics of an individual-based model. Finally, if ∆ t is the length of a timestep, and β = ˜ β ∆ t , γ = ˜ γ ∆ t , µ = ˜ µ ∆ t , then I ( t + ∆ t ) := I n +1 = (1 − I ( t )) ˜ β ∆ t + I ( t ) − (˜ µ + ˜ γ )∆ tI ( t ) + O (∆ t ) where I ( t ) := I n Hence, in the limit ∆ t →
0, we get the classical SIS ordinary differential equation dIdt = lim ∆ t → I ( t + ∆ t ) − I ( t )∆ t = (1 − I ) ˜ βI − (˜ µ + ˜ γ ) I. To understand the dynamics of (2.1), we can linearize at the origin to obtain the per-capitagrowth rate of the infection at the disease free-equilibrium(2.2) α = α ( µ, β, γ ) := βe − µ + e − µ − γ . The basic reproduction, alternatively, is given by(2.3) R = βe − µ / (1 − e − µ − γ ) . As α > R >
1, we can use the basic reproductive number to characterize theglobal dynamics, as the following theorem shows.
Theorem 2.1. (i) If R ≤ , then the origin is globally asymptotically stable. (ii) If R > , then there is a unique positive fixed point in (0 , such that it is globallyasymptotically stable in (0 , . SCHREIBER, HUANG, JIANG, AND HAO fraction I infected den s i t y A l l l l l + + Population size N e x t i n c t i on t i m e B0.0 0.2 0.4 0.6 0.8 1.0 fraction I infected den s i t y C l l l l l Population size N e x t i n c t i on t i m e D Figure 3.1.
Quasi-stationary distributions and intrinsic mean extinction times for the stochastic SIS modelfor R > (A,B) and R < (C,D). In (A) and (C), the quasi-stationary distributions approach a Diracdistribution at the equilibrium density (dashed line). In (B), the intrinsic mean time to extinction increasesexponentially with population size for R > . In (D), the intrinsic mean time to extinction saturates at / (1 − α ) (dashed line) as population size increases. Parameter values: γ = 0 . , µ = 0 . , and β = 0 . for(A,B) and β = 0 . for (C,D).
3. Metastability and extinction in a stochastic SIS model.
For the individual-basedstochastic model, we require the additional parameter of the total population size N . Giventhis population size, the state space corresponds to the possible fractions of infected individualsin the population S = (cid:26) , N , N , . . . , N − N , (cid:27) . Let I n ∈ S be the fraction of individuals infected at time step n . To determine the fractioninfected in the next time step, we assume that each infected individual remains infected withprobability e − µ − γ independent of each other and each susceptible individual lives and becomes N STOCHASTIC DISCRETE-TIME SIS AND SIR MODELS 5 infected with probability e − µ (1 − e − βI n ) independent of each other. Hence,(3.1) I n +1 = X n +1 + Y n +1 N where X n +1 ∼ Binom(NI n , e − µ − γ ) and Y n +1 ∼ Binom (cid:16)
N(1 − I n ) , e − µ (cid:0) − e − β I n (cid:1)(cid:17) . For all 1 ≤ i, j ≤ N , let Q ij = P [ X n +1 = j/N | X n = i/N ] be the transition probabilitiesrestricted to the transient states S \ { } and Q = ( Q ij ) be the associated N × N matrix.Unlike the deterministic model, the disease goes extinct in finite time with probability one forthe stochastic model. The following proposition follows from standard results in stochasticprocesses (see e.g. Harier, 2018, Theorem 3.20). Proposition 3.1.
Assume that µ + γ and β are positive. With probability one, I n = 0 forsome n ≥ . Even though extinction is inevitable, it may be preceded by long term transients. Tocharacterize these transients, we use quasi-stationary distributions introduced by Darroch &Seneta (1965). As Q is a sub-stochastic, positive matrix, there exists a dominant eigenvalue λ ∈ (0 ,
1) and associated dominant eigenvector π = ( π , . . . , π N ) (depending on N ) such that (cid:80) i π i = 1, π i > i , and πQ = λπ. π is the quasi-stationary distribution which satisfies(Darroch & Seneta, 1965) lim n →∞ P [ I n = j/N | I n >
0] = π j i.e. the probability of having j individuals in the long-term given the disease hasn’t goneextinct equals π j . Furthermore, N (cid:88) i =1 P [ I n +1 > | I n = i/N ] π i = λ i.e. λ is the probability the disease persists to the next time step given the process is followingthe quasi-stationary distribution. Hence, the mean time to extinction, when following thequasi-stationary distribution, is − λ , what Grimm & Wissel (2004) call the intrinsic meantime to extinction .Our main result for the stochastic SIS model concerns the behavior of the quasi-stationarydistribution and the intrinsic mean time to extinction for large population size N . Theorem 3.2.
Assume µ + γ > , β > . For each N ≥ , let π N be the quasi-stationarydistribution and λ N the corresponding eigenvalue. Let α = βe − µ + e − µ − γ . Then(i) If α ≤ (equivalently R ≤ ), then λ N ≤ α for all N ≥ and lim N →∞ π N = δ where δ is a Dirac measure at and convergence is in the weak* topology on probabilitymeasures on [0 , .(ii) If α > (equivalently R > ), then lim N →∞ π N = δ I ∗ where δ I ∗ is the Dirac measureat the unique positive fixed point I ∗ of equation (2.1) and lim sup N →∞ N log(1 − λ N ) < . SCHREIBER, HUANG, JIANG, AND HAO . . . fraction infected p r obab ili t y numericalanalytic A 0 20 40 60 80 + + Population size N e x t i n c t i on t i m e B Figure 3.2.
Quasi-stationary distributions and intrinsic mean extinction times for the stochastic SIS modelfor high contact rates and low recovery rates. In (A), the numerically computed quasi-stationary distribution(x marks) and the analytical approximation (solid blue curve) π i = (cid:0) Ni (cid:1) e − µi (1 − e − µ ) N − i / (1 − e − µ ) N . In (B),the numerically computed intrinsic mean time to extinction (x marks) and the analytical approximation (solidblue curve) / (1 − (1 − e − µ ) N ) . Parameter values: γ = 0 , µ = 1 . , β = 100 , and N = 20 for (A). The first assertion of Theorem 3.2 implies that if α < − α . We conjecture that in thelimit of large population size, N → ∞ , the intrinsic mean time to extinction equals − α . The second assertion of Theorem 3.2 implies that if α > c , c > − λ N ≥ c e c N for all N ≥ . Figure 3.1illustrates these conclusions numerically.Given that Theorem 3.2 describes the effect on increasing population size on the intrinsicmean time to extinction for a fixed value of α , it is natural to ask what effect increasing α hason these extinction times for a fixed population size. In general, this is a challenging question.However, we can answer this question for two special cases. First, we consider the case of lowrecovery rates γ = 0 and very high β (cid:29) β → ∞ , the updaterule for I n for I n > I n +1 ∼ Binom(N , e − µ ). Namely, provided there is atleast one infected individual at time step n , all individuals that have not died get infected.In this case, the quasi-stationary distribution is approximately π i = (cid:0) Ni (cid:1) e − µi (1 − e − µ ) N − i /λ N for i = 1 , , . . . , N with the persistence eigenvalue λ N = (1 − e − µ ) N . In particular, the meanintrinsic extinction time is bounded in the limit of α → ∞ due to β → ∞ . The accuracyof this approximation is illustrated in Figure 3.2. Second, in the limit of no recovery and nomortality (i.e. µ = γ = 0), the disease (not surprisingly!) never goes extinct whenever I > . Indeed, in this case, I n → n → ∞ with probability one provided β > I > . These two special cases highlight that the effect of α on the intrinsic mean time to extinctiondepends in a subtle way as α increases.
4. The dynamics of a deterministic SIR model.
As discrete-time counterpart to theclassical susceptible-infected-recovered (SIR) model, we assume all individuals escape natural
N STOCHASTIC DISCRETE-TIME SIS AND SIR MODELS 7 mortality with probability e − µ , infected individuals escape recovery with probability e − γ ,and susceptible individuals escape infection with probability e − βI where I is the frequency ofinfected individuals and β > (cid:40) S n +1 = 1 − e − µ + S n e − µ − βI n ,I n +1 = e − µ S n (cid:0) − e − βI n (cid:1) + e − µ − γ I n . Like the discrete-time SIS model, this discrete-time formulation of the SIR model has sev-eral advantages. First, by adding the two equations of (4.1) together, we obtain that thetrajectories of (4.1) remain in the domain X := { ( S, I ) : S ≥ , I ≥ , S + I ≤ } provided the initial value ( S , I ) lies in this domain. Furthermore, if we define ∂X := { ( S,
0) : 0 ≤ S ≤ } and X := X \ ∂X , then X and ∂X are positively invariant. Second,these dynamics, as described in the next section, correspond to the mean field dynamics ofan individual-based model. Finally, if ∆ t is the length of a time step, and β = ˜ β ∆ t , γ = ˜ γ ∆ t , µ = ˜ µ ∆ t , then S ( t + ∆ t ) := S n +1 = ˜ µ ∆ t + S ( t ) − (˜ µ + ˜ βI ( t ))∆ tS ( t ) + O (∆ t ) where S ( t ) := S n ,I ( t + ∆ t ) := I n +1 = S ( t ) I ( t ) ˜ β ∆ t + I ( t ) − (˜ µ + ˜ γ ) I ( t )∆ t + O (∆ t ) where I ( t ) := I n . Hence, in the limit ∆ t →
0, we get the classical SIR system of ordinary differential equations dSdt = lim ∆ t → S ( t + ∆ t ) − S ( t )∆ t = ˜ µ − (˜ µ + ˜ βI ) SdIdt = lim ∆ t → I ( t + ∆ t ) − I ( t )∆ t = ˜ βIS − (˜ µ + ˜ γ ) I. For our discrete-time SIR model, the disease-free fixed point is (1 , α = βe − µ + e − µ − γ and the reproductivenumber still equals R = βe − µ / (1 − e − µ − γ . We will show that if R >
1, the disease persistsand if R ≤
1, the disease-free equilibrium is globally stable. Furthermore, we will show whenthe recovery rate γ is sufficiently small, there is a globally stable endemic equilibrium. Tostate these results precisely, we define the parameter space as P := { λ := ( µ, β, γ ) : µ > , β > , γ ≥ } . Let C P := { λ = ( µ, β, ∈ P : α ( λ ) > } be the parameters corresponding to norecovery ( γ = 0) and α >
1. Finally, define C P := { λ ∈ P : α ( λ ) > , (4.1) admits a globally stable hyperbolic fixed point in X } . Using these definitions, we have the following theorem.
Theorem 4.1. (i) If α ≤ , then the disease free fixed point (1 , is globally asymptoticallystable. (ii) If α > , then F : X → X admits a global and compact attractor contained in theinterior of X . (iii) C P ⊃ C P is a non-empty open subset in P . In addition we conjecture that there is a globally stable endemic equilibrium when µ =˜ µ ∆ t , β = ˜ β ∆ t , γ = ˜ γ ∆, ∆ t > α > R > SCHREIBER, HUANG, JIANG, AND HAO susceptible frequency i n f e c t ed f r equen cy . . . A l l l l l l population size N m ean e x t i n c t i on t i m e Bsusceptible frequency i n f e c t ed f r equen cy . . . C l l l l l l population size N m ean e x t i n c t i on t i m e D Figure 5.1.
Quasi-stationary distributions and mean intrinsic extinction times for the stochastic SIRmodel. For parameters with α > , the quasi-stationary distributions, estimated numerically using the methodof Aldous et al. (1988), concentrate on the stable equilibrium as the population size goes from N = 100 (A)to N = 10 , (C). For this α > , the associated intrinsic mean extinction times increase exponentially withpopulation size in (B). For parameters with α < , the mean extinction times are bounded by / (1 − α ) in (D).Parameter values: µ = 0 . , γ = 0 . and β = 0 . for (A)-(C) and β = 0 . for (D).
5. Metastability and extinction in a stochastic SIR model.
As with the stochastic SISmodel, the stochastic SIR model requires the additional parameter of the total populationsize N . For a given population size, the state space S corresponds to the possible fraction ofsusceptible and infected individuals in the population, i.e. S = { ( i/N, j/N ) : i, j ∈ { , , . . . , N } , i + j ≤ N } ⊂ X. Let ( S n , I n ) ∈ S be the fractions of susceptible and infected individuals at time step n . Thefraction of removed individuals at time step n equals R n = 1 − S n − I n . Consistent with thedeterministic SIR model, we assume (i) each susceptible individual lives and becomes infectedwith probability e − µ (1 − e − βI n ) independent of each other, (ii) each infected individual eitherremains infected, dies and gets replaced with a susceptible individual, or enters the removedclass with probabilities e − µ − γ , 1 − e − µ , or e − µ (1 − e − γ ) independent of each other, and (iii)each removed individual dies and creates a new susceptible individual with probability 1 − e − µ .To account for these transitions, let W n +1 be a binomial random variable with N S n trials and N STOCHASTIC DISCRETE-TIME SIS AND SIR MODELS 9 probability of success e − µ (1 − e − βI n ) (i.e. susceptible individuals that will become infected), X n +1 be a binomial random variable with N I n trails and probability of success 1 − e − µ (i.e.infected individuals that die and get replaced by a susceptible individual), Y n +1 be a binomialrandom variable with N I n − X n +1 trials with probability of success e − γ (i.e. non-dying infectedindividuals that will not enter the removed class), and Z n +1 be a binomial random variablewith N (1 − I n − S n ) trials with probability of success 1 − e − µ (i.e. removed individuals thatdie and get replaced with a susceptible individual). Under these assumptions, the stochasticSIR model is(5.1) S n +1 = 1 N ( N S n − W n +1 + X n +1 + Z n +1 ) ,I n +1 = 1 N ( W n +1 + Y n +1 ) , where W n +1 ∼ Binom(NS n , e − µ (cid:0) − e − β I n (cid:1) ) ,X n +1 ∼ Binom(NI n , − e − µ ) ,Y n +1 ∼ Binom(NI n − X n+1 , e − γ ) , and Z n +1 ∼ Binom(N(1 − I n − S n ) , − e − µ ) . As with the stochastic SIS model, the disease goes extinct in finite time with probabilityone for the stochastic model. The following proposition follows from standard results instochastic processes (see e.g. Harier, 2018, Theorem 3.20).
Proposition 5.1.
Assume that µ + γ and β are positive. With probability one, I n = 0 forsome n ≥ . To characterize metastability and extinction times, define S + = { ( x , x ) ∈ S : x > } to be all the states where the disease persists. For all pairs of states x, y ∈ S + , let Q xy = P [( S n +1 , I n +1 ) = y | ( S n , I n ) = x ] be the transition probabilities restricted to the transientstates and Q = ( Q xy ) x,y ∈S + be the associated matrix. Let π = ( π x ) x ∈S + be the quasi-stationary distribution with associated persistence probability λ i.e. (cid:80) x ∈S + π x = 1, π x > x ∈ S and πQ = λπ. Our main result for the stochastic SIR model concerns thebehavior of the quasi-stationary distribution and the intrinsic mean time to extinction forlarge population size N . Theorem 5.2.
Assume µ + γ > , β > . For each N ≥ , let π N be the quasi-stationarydistribution and λ N the corresponding eigenvalue for (5.1) . Let α = βe − µ + e − µ − γ . Then(i) If α < , then λ N ≤ α for all N ≥ and lim N →∞ π N = δ (1 , where δ (1 , is a Diracmeasure at disease-free equilibrium (1 , and convergence is in the weak* topology.(ii) If α > , then lim sup N →∞ N log(1 − λ N ) < and there exists a compact set K ⊂ (0 , such that π ∗ ( K ) = 1 for every weak* limit point π ∗ of π N as N → ∞ , and where π ∗ is invariant for the deterministic model (4.1) . Moreover, if ( µ, β, γ ) ∈ C P , then lim N →∞ µ N = δ ( S ∗ ,I ∗ ) where δ ( S ∗ ,I ∗ ) is the Dirac measure at the unique positive fixedpoint ( S ∗ , I ∗ ) of equation (4.1) . The first assertion of Theorem 5.2 implies that if α < mean time to extinction after these transients is less than − α . Furthermore, whenever per-manence of the deterministic model corresponds to a globally stable equilibrium (see Theo-rem 4.1), the QSDs concentrate on a Dirac measure at this equilibrium i.e. it supports theonly invariant measure in K for the deterministic dynamics. The second assertion of Theo-rem 5.2 implies that if α > c , c > − λ N ≥ c e c N for all N ≥ . Figure 5.1 illustrates these conclusions.
6. Discussion.
This paper formulates and provides a mathematically rigorous analysis ofdeterministic and stochastic, discrete-time SIS and SIR models. The stochastic models arebased on probabilistic, individual-based update rules. The conditional expected change inthe fraction of infected and susceptible individuals given the current values of these fractionsdetermines the update rule for the deterministic model and, in this sense, the deterministicmodel is the mean field model for the stochastic models.Many earlier discrete-time epidemic models of SIS and SIR dynamics have been derivedusing numerical approximation schemes for differential equations (Allen, 1994; Ma et al. ,2013; Satsuma et al. , 2004; Enatsu et al. , 2010; Castillo-Chavez & Yakubu, 2001). Thesemodels, including the one-dimensional ones, can exhibit oscillatory dynamics. In contrast,our model, which is based on individual-based update rules and uses an exponential escapefunction, is most similar to higher dimensional, discrete-time epidemiological models thathave been used for applications to specific diseases (Emmert & Allen, 2004, 2006; Allen &van den Driessche, 2008). Unlike the models based on numerical approximation schemes, ouranalysis and numerical explorations suggest that our models always approach a global stableequilibrium. When R ≤
1, we prove that all trajectories asymptotically approach the disease-free equilibrium for both the SIS and SIR models. When R >
1, we prove the disease persistsfor both models, approaches a globally stable, endemic equilibrium for the SIS model, andprovide sufficient conditions for global stability of the endemic equilibrium for the SIR model.Extensive numerical simulations suggest that this endemic equilibrium of the SIR model isglobally stable whenever R >
1. Hence, we conjecture that R > R > R ≤
1, the frac-tion of infected in the stochastic model becomes zero in finite time for all values of R . Tounderstand this well-know discrepancy (Bartlett, 1966; Keeling & Rohani, 2011; Diekmann et al. , 2012) for our model, we characterized the long-term transients using quasi-stationarydistributions for finite-state, discrete-time Markov chains (Darroch & Seneta, 1965). For thesecharacterizations, we used the per-capita growth rate α = βe − µ + e − µ − γ of the infection atthe disease free equilibrium. When α < R < ≤ / (1 − α ) for all population sizesand for both the SIS and SIR models. Indeed, we conjecture that as N → ∞ , this meantime to extinction converges to 1 / (1 − α ). While R < α <
1, we have α = (1 − e − µ − γ ) R + e − µ − γ > R whenever R <
1. Hence, even if R is very small, the meantimes to extinction can be arbitrarily long. For example, given any 0 < x < y <
1, we can
N STOCHASTIC DISCRETE-TIME SIS AND SIR MODELS 11 make R = x and α = y by choosing γ = 0, e − µ = y/
2, and β = x (1 − e − µ ).When R > α > N and the quasi-stationary distributions concentrateon positive invariant sets for the deterministic model for large N . In particular, coupledwith our analysis of the deterministic dynamics, our results imply that the quasi-stationarybehavior for large N always concentrates near the globally stable, endemic equilibrium ofthe SIS model. We provide sufficient conditions for the same conclusion for the SIR model,and conjecture that this always occurs for the SIR model. These conclusions are consistentwith earlier studies of continuous-time Markov models where the analysis was done usingdiffusion approximations of the individual-based models (Barbour, 1975; Kryscio & Lef`evre,1989; N˚asell, 1996, 1999; Andersson & Britton, 2000; N˚asell, 2002; Lindholm & Britton, 2007;Andersson & Lindenstrand, 2011; Clancy & Tjia, 2018). In contrast, our results apply largedeviation methods from (Faure & Schreiber, 2014) directly to the individual-based models. Anopen question for the stochastic model with R > N . Specifically, what is the value of α such thatthe mean time to extinction grows like exp( αN ) for large N ? The diffusion approximationsprovide one approach to find potential candidates for α. When R >
1, we found that the mean time to extinction can be arbitrarily large even fora fixed population size. For example, this occurs when recovery and mortality rates approachzero in which case R also increases without bound but α remains bounded above by β + 1. Incontrast, increasing contact rates (which increase α and R without bound) leads to extinctiontimes that are constrained by population size, recovery rates and mortality rates.In addition to the open mathematical questions that we have already raised, future chal-lenges include analyzing extensions of our models. These extensions could include additionalcompartments such as SEIR models, multi-age group epidemic models, and SIR type modelswith vaccination (Anderson & May, 1991; Keeling & Rohani, 2011; Kong et al. , 2015). Moregenerally, when the discrete-time system is autonomous, the mathematical approaches usedhere should be applicable to study quasi-stationary distributions and the intrinsic extinctiontimes. However, when population sizes or transmission rate change stochastically over time(Anderson & May, 1979; Pollicott et al. , 2012), new mathematical approaches are requiredfor studying the impact of these environmentally driven random fluctuations on intrinsic ex-tinction times.
7. Proof of Theorem 2.1.
Suppose that α <
1. Then we shall prove that(7.1) F ( I ) := e − µ (1 − I ) (cid:0) − e − βI (cid:1) + e − µ − γ I < (cid:0) βe − µ + e − µ − γ (cid:1) I =: L ( I ) , I ∈ (0 , . It is easy to see that (7.1) is equivalent to(7.2) (1 − I ) (cid:0) − e − βI (cid:1) < βI, I ∈ (0 , . Let g ( I ) := (1 − I ) (cid:0) − e − βI (cid:1) − βI . Then g (0) = 0, g (cid:48) ( I ) = − (1 + β ) + e − βI + β (1 − I ) e − βI , g (cid:48) (0) = 0 , and g (cid:48)(cid:48) ( I ) = − βe − βI (cid:0) β (1 − I ) (cid:1) < , I ∈ [0 , . Therefore, g (cid:48) ( I ) < g (cid:48) (0) = 0 , I ∈ (0 , . Furthermore, g ( I ) < g (0) = 0 , I ∈ (0 , . This shows that (7.2) holds.Fix I ∈ [0 ,
1] and set I n := F n ( I ) , J n := L n ( I ) = α n I , n = 1 , , · · · . We claim that(7.3) I n ≤ J n , n = 1 , , · · · . For n = 1, (7.3) follows immediately from (7.1). Assume that (7.3) holds for n . Then by (7.1)and the increasing of L , we have that I n +1 = F ( I n ) ≤ L ( I n ) ≤ L ( J n ) = J n +1 . By mathematical induction, (7.3) is valid for all positive integers. Since α < L n ( I ) = α n I → n →
0, i.e. 0 is globally asymptotically stable.Suppose that α = 1. Then it follows from (7.1) that F ( I ) < I, I ∈ (0 , I ]. UsingFeigenbaum’s method given in (ii), we can easily prove that 0 is still globally asymptoticallystable in this case.(ii) Suppose that α >
1. Then it is not difficult to prove that F ( I ) has a unique positivefixed point, denoted by I ∗ , and F ( I ) has a unique positive critical point I ∗ c . We will divide(ii) into two cases: (iia) I ∗ ≤ I ∗ c , (iib) I ∗ > I ∗ c . We use Feigenbaum’s method by depicting graphs to show that all nontrivial trajectoriesconverge to the fixed point I ∗ .(iia) Take I ∈ (0 , I ∗ ). As depicted in Figure 7.1(a), the iterating sequence I n increasinglytends to I ∗ . If I ∈ ( I ∗ , I ∗ c ] with I ∗ < I ∗ c , then Figure 7.1(b) shows that the iterating sequence I n decreasingly tends to I ∗ . If I > I ∗ c , then the first iteration I ∈ (0 , I ∗ ). After the firstiteration, I n increasingly tends to I ∗ , see Figure 7.1(c).(iib) It is easy to see that F ( I ) is decreasing when I > I ∗ c , that is, F (cid:48) ( I ) ≤ , I ∈ ( I ∗ c , F (cid:48) ( I ) ≥ e − µ (cid:0) − e − βI + e − γ (cid:1) > − . Therefore, − < F (cid:48) ( I ) ≤ , I ∈ ( I ∗ c , I < I < I ∗ if I ∗ c < I < I ∗ . Let K be a line through ( I ∗ , F ( I ∗ ))whose gradient is −
1. Because we have proved F (cid:48) ( I ) > − S = F ( I ) is between K and S = I as shown in Figure 7.2(a).Plotting a line which is perpendicular to I -axis through ( I , S = I, S = F ( I ) and K at A, E and B , respectively. Through B and E , we add lines parallel to I -axis, which intersect S = I at C and G , respectively. Because the ordinate of B is strictlygreater than the ordinate of E , the abscissa of C is strictly greater than the abscissa of G . N STOCHASTIC DISCRETE-TIME SIS AND SIR MODELS 13
Sketch the lines perpendicular to I -axis through C and G , which intersect S = F ( I ) at J and H , respectively. Because S = F ( I ) is monotonically decreasing on ( I ∗ c , J are strictly less than the ordinate of H . Now we extend the segment CJ such that it intersects K at D . The ordinate of D is strictly less than the ordinate of H . Sketch the lines parallel to I -axis through D and H . Then they intersect S = I at A and L , respectively. The ordinateof L is strictly greater than the ordinate of A . Then the abscissa of L is strictly greater thanthe abscissa of A , that is, I ≤ I ≤ I ∗ as shown Figure 7.2(b).Now, we want to show lim n →∞ F n ( I ) = I ∗ by contradiction. Suppose it is not true. Thenlim n →∞ F n ( I ) = I ∗∗ < I ∗ . Thus we set I ∗∗ to be the initial point and repeat the aboveprocess, then we have F ( I ∗∗ ) > I ∗ . This is contradictory to lim n →∞ F n ( I ) = I ∗∗ , therefore,lim n →∞ F n ( I ) = I ∗ . Similarly, we can prove lim n →∞ F n +1 ( I ) = I ∗ . This concludes that I n oscillates around I ∗ and converges to I ∗ .The proof of the case I < I ∗ c is given in the Figure 7.2(c). Figure 7.1. I ∗ ≤ I ∗ c . Figure 7.2. I ∗ > I ∗ c .
8. Proof of Theorem 3.2.
We use results from Faure & Schreiber (2014) to prove The-orem 3.2. We begin by verifying that Standing Hypothesis 1.1 of Faure & Schreiber (2014)holds. In their notation, “ ε ” corresponds to N in our models i.e. small demographic noisecorresponds to large population size N . For all δ > N ≥
1, define β δ ( N ) = sup x ∈ [0 , P [ | I n +1 − F ( x ) | ≥ δ | I n = x ] where F ( x ) = e − γ − µ x + e − µ (1 − e − βx )(1 − x ) corresponds to the right hand side of thedeterministic model in equation (2.1). Standing Hypothesis 1.1 of Faure & Schreiber (2014)requires that lim N →∞ β δ ( N ) = 0 for all δ > . The following proposition proves somethingstronger using large deviation estimates.
Proposition 8.1.
There exists a function ρ : (0 , ∞ ) → (0 , ∞ ) such that β δ ( N ) ≤ exp( − N ρ ( δ )) for all N ≥ and δ > . Proof.
While we use standard large deviation estimates, we go through the details toensure that the estimates can be taken uniformly in x ∈ [0 , . Define a = e − µ − γ and b ( x ) = e − µ (1 − e − βx ). By the exponential Markov inequality, we have for all t (8.1) P [ N ( I n +1 − F ( x )) ≤ N δ | I n = x ] ≤ e − tNδ E [ e t ( N ( I n +1 − F ( x ))) | I n = x ]= e − tF ( x ) N − tδN (cid:0) − a + ae t (cid:1) Nx (cid:0) − b ( x ) + b ( x ) e t (cid:1) N (1 − x ) Define ψ ( x, t ) = − tF ( x ) + x log(1 − a + ae t ) + (1 − x ) log(1 − b ( x ) + b ( x ) e t ) . Taking log of equation (8.1) and dividing by N yields(8.2) 1 N log P [ N ( I n +1 − F ( x )) ≥ N δ | I n = x ] ≤ − δt + ψ ( t, x )for all t (cid:54) = 0, x ∈ [0 , δ >
0. Similarly, one can show that(8.3) 1 N log P [ N ( F ( x ) − I n +1 ) ≥ N δ | I n = x ] ≤ − δt + ψ ( − t, x )for all t (cid:54) = 0, x ∈ [0 , δ >
0. We have ∂ψ∂t = − F ( x ) + xae t − a + ae t + (1 − x ) b ( x ) e t − b ( x ) + b ( x ) e t and ∂ ψ∂t = ax e t (1 − a )(1 − a + ae t ) + (1 − x ) b ( x ) e t (1 − b ( x ))1 − b ( x ) + b ( x ) e t > . As ψ (0 , x ) = 0 = ∂ψ∂t (0 , x ) = 0, and ψ is strictly convex in t , for all δ >
0, there exists t ∗ ( δ ) > δt ∗ ( δ ) > ψ ( − t ∗ ( δ ) , x ) and δt ∗ ( δ ) > ψ ( t ∗ ( δ ) , x ) for all x ∈ [0 , . Define ρ ( δ ) = δt ∗ ( δ ) − max x ∈ [0 , ψ ( t ∗ ( δ ) , x ) . Equations (8.2)–(8.3) imply that P [ | I n +1 − F ( x ) | ≥ δ | I n = x ] ≤ exp( − N ρ ( δ ))for all x ∈ [0 ,
1] and δ > . N STOCHASTIC DISCRETE-TIME SIS AND SIR MODELS 15
To prove the first result of Theorem 3.2, assume that α ≤
1. Theorem 2.1 implies that 0is globally stable for the deterministic model I (cid:55)→ F ( I ). Theorem 3.12 of Faure & Schreiber(2014), which only requires Standing Hypothesis 1.1, implies that lim N →∞ π N = δ . Define R ( x ) = F ( x ) /x for x ∈ (0 , R ( x ) ≤ α for x ∈ (0 , . For N ≥
1, quasi-stationarity of π N implies λ N N (cid:88) i =1 iN π Ni = N (cid:88) i =1 E (cid:20) I n +1 (cid:12)(cid:12)(cid:12) I n = iN (cid:21) π Ni = N (cid:88) i =1 F (cid:18) iN (cid:19) π Ni = N (cid:88) i =1 iN R (cid:18) iN (cid:19) π Ni ≤ α N (cid:88) i =1 iN π Ni . Since (cid:80) Ni =1 iN π Ni > λ N ≤ α for all N ≥ α > I ∗ for the map I (cid:55)→ F ( I ).Assertion (b) of Lemma 3.9 of Faure & Schreiber (2014), implies that there exists δ > − λ N ≤ β δ ( N ) for all N ≥
1. Proposition 8.1 implies thatlim sup N →∞ N log(1 − λ N ) ≤ − ρ ( δ ) . To complete the proof of the second assertion, we need to verify the assumption in As-sertion (b’) of Lemma 3.9 of Faure & Schreiber (2014). Choose η > x ∈ [0 ,η ] (1 − e − µ − γ ) x ≥ exp( − ρ ( δ ) /
3) and min x ∈ [0 ,η ] (1 − e − µ (1 − e − βx )) − x ≥ exp( − ρ ( δ ) / . Then min x ∈ [0 ,η ] P [ I n +1 = 0 | I n = x ] ≥ exp( − N ρ ( δ ) / N →∞ β δ ( N )min x ∈ [0 ,η ] P [ I n +1 = 0 | I n = x ] ≤ lim N →∞ exp( − N ρ ( δ ))exp( − N ρ ( δ ) /
3) = 0which verifies the assumption of (b’) of Lemma 3.9 of Faure & Schreiber (2014) and impliesthat lim N →∞ (cid:88) i/N ≤ η π Ni = 0i.e. for any weak* limit point π ∗ of π N , π ∗ ([0 , η ]) = 0. As λ N →
1, Proposition 3.11 of Faure& Schreiber (2014) implies that any weak* limit point of π N is invariant for the dynamics x (cid:55)→ F ( x ). As these weak* limit points are supported on [ η,
1] and the only invariant measurefor x (cid:55)→ F ( x ) on this interval is the Dirac measure δ I ∗ and the unique positive fixed point, itfollows that π N converges in the weak* topology to δ I ∗ as claimed.
9. Proof of Theorem 4.1.
Let E f = (1 , F ( S, I ) := (cid:18) − e − µ + Se − µ − βI e − µ S (cid:0) − e − βI (cid:1) + e − µ − γ I (cid:19) . Then DF ( E f ) := (cid:18) e − µ − βe − µ βe − µ + e − µ − γ (cid:19) . It follows that the per-capita growth rate of the disease is still given by (2.2) and the basicreproduction number of (4.1) is still the expression in (2.3). Define the competitive cone K := { ( u, v ) : u ≤ , v ≥ } . Then it is easy to check that DF ( E f ) keeps K invariant (see Wang & Jiang (2001)). Define L ( S, I ) := (cid:18) (cid:19) + (cid:18) e − µ − βe − µ βe − µ + e − µ − γ (cid:19) (cid:18) S − I (cid:19) . We shall verify that(9.1) F ( S, I ) ≤ K L ( S, I ) , ( S, I ) ∈ X. (9.1) is equivalent to(9.2) (cid:40) − e − µ + Se − µ − βI ≥ e − µ ( S − − βIe − µ ,e − µ S (cid:0) − e − βI (cid:1) + e − µ − γ I ≤ (cid:0) βe − µ + e − µ − γ (cid:1) I on X . By simple computation, (9.2) is equivalent to(9.3) S (1 − e − βI ) ≤ βI, ( S, I ) ∈ X. Since (
S, I ) ∈ X , S ≤ (1 − I ). Thus (9.3) follows immediately from (7.2), that is, thecompetitive ordering relation (9.1) holds.Let P := ( S , I ) ∈ X, P n := F n ( P ) = ( S n ( P ) , I n ( P )) , Q n := L n ( P ). We shall showthat(9.4) (1 , τ ≤ K P n ≤ K Q n , n = 1 , , · · · . The left inequality is obvious by the definition of competitive order and X . So we will provethe right one. (9.1) deduces that (9.4) holds for n = 1. Suppose that (9.4) holds for n . Thenusing (9.1) and the order-preserving for DF ( E f ) , we obtain P n +1 = F ( P n ) ≤ K L ( P n ) ≤ K L ( Q n ) = Q n +1 . N STOCHASTIC DISCRETE-TIME SIS AND SIR MODELS 17
By mathematical induction, (9.4) holds.Let α <
1. Then Q n → (1 ,
0) as n → ∞ . Therefore, (i) is proved by (9.4).Suppose α >
1. Then we shall prove that the system (4.1) is uniformly persistent withrespect to ( X , ∂X ), that is, there exists η > n →∞ I n ( P ) ≥ η, for all P = ( S , I ) ∈ X . It is easy to see that E f is the maximal compact invariant set in ∂X which is positivelyinvariant with respect to F and lies on the stable manifold of E f . Recalling the Hofbauerand So uniform persistence theorem of Hofbauer & So (1989), the system (4.1) is uniformlypersistent if and only if(a) E f is isolated in X , and(b) W s ( E f ) ⊂ ∂X .Since the disease free fixed point E f is hyperbolic and a saddle, (a) and (b) follows immediatelyfrom the Hartman and Grobman theorem (see Guckenheimer & Holmes (1983)). This verifiesthe uniform persistence and hence the system (4.1) admits an attractor in X .It follows from (9.5) that (4.1) contains a compact attractor A ⊂ { ( S, I ) ∈ X : I ≥ η } .Besides, by the first equality of (4.1), we have S n ( P ) ≥ − e − µ for n ≥ P ∈ X . Thisimplies that A ⊂ { ( S, I ) ∈ X : S ≥ − e − µ } . As a result, A ⊂ { ( S, I ) ∈ X : S ≥ − e − µ , I ≥ η, S + I ≤ } . This proves (ii).It remains to prove (iii). First, we consider the system (4.1) with λ = ( µ , β ,
0) and α ( µ , β , >
1. We shall prove that it admits a globally stable fixed point (1 − I ∗ , I ∗ ) in X ,where F ( I ∗ ) = I ∗ with 0 < I ∗ < n := S n + I n . Then from (4.1) it follows that(9.6) ∆ n +1 = 1 − e − µ + e − µ ∆ n . It is easy to see that (9.6) has the positive fixed point 1 and all positive orbits tend to 1 as n → + ∞ . Therefore, the system (4.1) is reduced to the system (2.1) with µ = µ , β = β , γ =0. Applying Theorem 2.1(ii), we get that the system (2.1) has a globally stable fixed point I ∗ in (0 , − I ∗ , I ∗ ) in X , where F ( I ∗ ) = I ∗ with 0 < I ∗ <
1. Recalling the proof of Theorem 2.1(ii)(see Figure 7.1 and Figure7.2(a)), we have(9.7) | F (cid:48) ( I ∗ ) | = (cid:12)(cid:12) e − µ − β I ∗ (cid:0) β (1 − I ∗ ) (cid:1)(cid:12)(cid:12) < . Next, we will prove the spectral radius of the Jacobian matrix for F ( S, I ) at the positivefixed point E ∗ ( S ∗ , I ∗ ) := (1 − I ∗ , I ∗ ) is strictly less than 1.An easy calculation yields that DF ( E ∗ ) := (cid:18) e − µ − β I ∗ − β S ∗ e − µ − β I ∗ e − µ (1 − e − β I ∗ ) β S ∗ e − µ − β I ∗ + e − µ (cid:19) , (cid:18) e − µ − β I ∗ − β S ∗ e − µ − β I ∗ e − µ (1 − e − β I ∗ ) β S ∗ e − µ − β I ∗ + e − µ (cid:19) (cid:18) − (cid:19) = F (cid:48) ( I ∗ ) (cid:18) − (cid:19) and det DF ( E ∗ ) = e − µ F (cid:48) ( I ∗ ). This proves that F (cid:48) ( I ∗ ) and e − µ are two eigenvalues of DF ( E ∗ ), that is, the spectral radius of DF ( E ∗ ) is strictly less than 1.In what follows, we shall use Theorem 2.1 of Smith & Waltman (1999) to prove that C P is open in the parameter space P .For this purpose, denote by (cid:107)·(cid:107) the Euclidean norm of R and B C ( λ , s ) the open ball in P of radius s about the point λ . We fix a λ ∈ C P and an s ∈ (0 , µ ) such that α ( λ ) > λ ∈ B C ( λ , s ). In order to consider the perturbed systems for parameters, we set F λ ( S, I )the given mapping and F λ ( S, I ) = ( S λ ( S, I ) , I λ ( S, I )) the mappings for λ ∈ B C ( λ , s ). Define R λ ( S, I ) = (cid:40) I λ ( S,I ) I , if I > ,βe − µ S + e − µ − γ , if I = 0 . Then R λ ( S, I ) is continuous on X . By induction, it not difficult to prove that(9.8) F nλ ( S,
0) = (cid:0) − e − nµ + e − nµ S, (cid:1) R λ ( F nλ ( S, βe − µ (cid:0) − e − nµ + e − nµ S (cid:1) + e − µ − γ for n = 1 , , · · · . We claim that there exist an s ∈ (0 , s ], an integer N > δ > ρ > λ , such that(9.9) I Nλ ( S, I ) ≥ ρI for all λ ∈ B C ( λ , s ) and I ∈ [0 , δ ] , where F nλ ( S, I ) = ( S nλ ( S, I ) , I nλ ( S, I )).From (9.8) it follows that R λ ( F nλ ( S, ≥ α ( λ ) − βe − ( n +1) µ . By the continuity of α ( λ ), there exists an s ∈ (0 , s ] such that α ( λ ) > α ( λ )+12 for all λ ∈ B C ( λ , s ), and hence R λ ( F nλ ( S, ≥ α ( λ ) + 12 − ( β + s ) e − ( n +1)( µ − s ) for all λ ∈ B C ( λ , s ) . This implies that there is an integer
N >
0, only depending on λ , such that(9.10) R λ ( F Nλ ( S, > α ( λ ) + 34 for all λ ∈ B C ( λ , s ) . Using (9.10) and the uniform continuity of R λ ( F Nλ ( S, I )) on B C ( λ , s ) × X , we obtain thatthere is a δ >
0, only depending on λ , such that R λ ( F Nλ ( S, I )) > α ( λ ) + 78 := ρ for all λ ∈ B C ( λ , s ) and 0 ≤ I ≤ δ, N STOCHASTIC DISCRETE-TIME SIS AND SIR MODELS 19 which implies that (9.9) holds, thus the claim is proved.By (9.9), we get that for each λ ∈ B C ( λ , s ), I mNλ ( S, I ) ≥ ρ m I if F kNλ ( S, I ) ∈ [0 , × (0 , δ ] for k = 0 , , · · · , m − . This shows that there exists at least a positive integer m with the property F mNλ ( S, I ) ∈ [0 , × ( δ,
1] if (
S, I ) ∈ [0 , × (0 , δ ] . Let U = X , Λ = B C ( λ , s ) , B λ = [0 , × [ δ, B C ( λ , s ) ⊂ C P . The proof iscomplete.
10. Proof of Theorem 5.2.
As in the proof of Theorem 3.2, we use results from Faure &Schreiber (2014) to prove Theorem 5.2. We begin by verifying that Standing Hypothesis 1.1of Faure & Schreiber (2014) holds. For all δ > N ≥
1, define β δ ( N ) = sup x,y ≥ ,x + y ≤ P [ (cid:107) ( S n +1 , I n +1 − F ( x, y ) (cid:107) ∞ ≥ δ | ( S n , I n ) = ( x, y )]where F ( x, y ) = (1 − e − µ + e − µ − βy x, xe − µ (1 − e − βy ) + ye − µ − γ ) corresponds to the righthand side of the deterministic model in equation (4.1) and (cid:107) ( x, y ) (cid:107) ∞ = max {| x | , | y |} cor-responds to the sup norm. Standing Hypothesis 1.1 of (Faure & Schreiber, 2014) requiresthat lim N →∞ β δ ( N ) = 0 for all δ > . Like Proposition 8.1 for the stochastic SIS model, thefollowing proposition proves something stronger using large deviation estimates.
Proposition 10.1.
There exists a function ρ : (0 , ∞ ) → (0 , ∞ ) such that β δ ( N ) ≤ exp( − N ρ ( δ )) for all N ≥ and δ > . Proof.
We begin by observing that
N S n − W n +1 and X n +1 + Z n +1 in (5.1) conditioned on( S n , I n ) = ( x, y ) are independent binomials where N S n − W n +1 has N x trials with probabilityof success a ( y ) = 1 − e − µ (1 − e − βy ) and X n +1 + Z n +1 has N (1 − x ) trials with probabilityof success b = 1 − e − µ . Using the exponential Markov inequality as in the proof of Proposi-tion 8.1, we get(10.1) 1 N log P [ N ( S n +1 − F ( x, y )) ≥ N δ | ( S n , I n ) = ( x, y )] ≤ − δt + ψ ( t, x, y )1 N log P [ N ( F ( x, y ) − S n +1 ) ≥ N δ | ( S n , I n ) = ( x, y )] ≤ − δt + ψ ( − t, x, y )for all t, δ and where ψ ( t, x, y ) = − tF ( x, y ) + x log(1 − a ( y ) + a ( y ) e t ) + (1 − x ) log(1 − b + b e t ) . As ψ (0 , x, y ) = 0 = ∂ψ∂t (0 , x, y ) = 0, and ψ is strictly convex in t , for all δ >
0, there exists t ∗ ( δ ) > δt ∗ ( δ ) > ψ ( − t ∗ ( δ ) , x, y ) and δt ∗ ( δ ) > ψ ( t ∗ ( δ ) , x, y ) for all x, y ∈ [0 , x + y ≤
1. Define ρ ( δ ) = δt ∗ ( δ ) − max x,y ∈ [0 , ,x + y ≤ ψ ( t ∗ ( δ ) , x, y ) > . Equation (10.1) implies that P [ | S n +1 − F ( x, y ) | ≥ δ | ( S n , I n ) = ( x, y )] ≤ exp( − N ρ ( δ ))for all x, y ∈ [0 ,
1] with x + y ≤ δ > .W n +1 and Y n +1 conditioned on ( S n , I n ) = ( x, y ) in equation (5.1) are also independentbinomial random variables where W n +1 has N x trials with probability of success e − µ (1 − e − βy )and Y n +1 has N y trials with probability of success e − γ . Therefore using the exponentialMarkov inequality, one can show there exists a function ρ : (0 , ∞ ) → (0 , ∞ ) such that P [ | I n +1 − F ( x, y ) | ≥ δ | ( S n , I n ) = ( x, y )] ≤ exp( − N ρ ( δ ))for all x, y ∈ [0 ,
1] with x + y ≤ δ > . Setting ρ ( δ ) = min { ρ ( δ ) , ρ ( δ ) } completes theproof of the proposition.To prove the first result of Theorem 5.2, assume that α ≤
1. Theorem 4.1 implies that(1 ,
0) is globally stable for the deterministic model (
S, I ) (cid:55)→ F ( S, I ). Theorem 3.12 of Faure& Schreiber (2014), which only requires Standing Hypothesis 1.1, implies that lim N →∞ π N = δ (1 , in the weak* topology. Define R ( x, y ) = F ( x, y ) /y for y ∈ (0 , , x ∈ [0 , x + y ≤ R ( x, y ) ≤ α . For N ≥
1, quasi-stationarity of π N implies λ N (cid:88) x,y ∈S + yπ Nx,y = (cid:88) x,y ∈S + E (cid:104) I n +1 (cid:12)(cid:12)(cid:12) ( S n , I n ) = ( x, y ) (cid:105) π Nx,y = (cid:88) x,y ∈S + F ( x, y ) π Nx,y = (cid:88) x,y ∈S + yR ( x, y ) π Nx,y ≤ α (cid:88) x,y ∈S + yπ Nx,y
Since (cid:80) x,y ∈S + yπ Nx,y > λ N ≤ α for all N ≥ α > K ⊂ (0 , × (0 ,
1) for (
S, I ) (cid:55)→ F ( S, I ).Assertion (b) of Lemma 3.9 of Faure & Schreiber (2014), implies that there exists δ > − λ N ≤ β δ ( N ) for all N ≥
1. Proposition 10.1 implies thatlim sup N →∞ N log(1 − λ N ) ≤ − ρ ( δ ) . Assumption in Assertion (b’) of Lemma 3.9 of Faure & Schreiber (2014) holds for an argumentsimilar to the proof of Theorem 3.2 and consequently any weak* limit point π ∗ of π N satisfies π ∗ ( K ) = 1. As λ N →
1, Proposition 3.11 of Faure & Schreiber (2014) implies that any weak*limit point of π N is invariant for the dynamics ( x, y ) (cid:55)→ F ( x, y ). N STOCHASTIC DISCRETE-TIME SIS AND SIR MODELS 21
References.
Aldous, D., Flannery, B. & Palacios, J.L. (1988). Two applications of urn processes thefringe analysis of search trees and the simulation of quasi-stationary distributions of Markovchains.
Probability in the Engineering and Informational Sciences , 2, 293–307.Allen, L. & van den Driessche, P. (2008). The basic reproduction number in some discrete-timeepidemic models.
Journal of Difference Equations and Applications , 14, 1127–1147.Allen, L.J.S. (1994). Some discrete-time SI, SIR, and SIS epidemic models.
MathematicalBiosciences , 124, 83–105.Allen, L.J.S. & Burgin, A.M. (2000). Comparison of deterministic and stochastic SIS and SIRmodels in discrete time.
Mathematical Biosciences , 163, 1–33.Anderson, R.M. & May, R. (1991).
Infectious diseases of humans: dynamics and control .Oxford University Press.Anderson, R.M. & May, R.M. (1979). Population biology of infectious diseases: Part i.
Nature ,280, 361–367.Andersson, H. & Britton, T. (2000). Stochastic epidemics in dynamic populations: quasi-stationarity and extinction.
Journal of Mathematical Biology , 41, 559–580.Andersson, P. & Lindenstrand, D. (2011). A stochastic sis epidemic with demography: initialstages and time to extinction.
Journal of Mathematical Biology , 62, 333–348.Barbour, A.D. (1975). The duration of the closed stochastic epidemic.
Biometrika , 62, 477–482.Bartlett, M. (1966).
An Introduction to Stochastic Processes With Special Reference to Meth-ods and Applications . Cambridge University Press, London-New York.Castillo-Chavez, C. & Yakubu, A.A. (2001). Discrete-time sis models with complex dynamics.
Nonlinear Analysis, Theory, Methods and Applications , 47, 4753–4762.Clancy, D. & Tjia, E. (2018). Approximating time to extinction for endemic infection models.
Methodology and Computing in Applied Probability , 20, 1043–1067.Darroch, J.N. & Seneta, E. (1965). On quasi-stationary distributions in absorbing discrete-time finite markov chains.
Journal of Applied Probability , 2, 88–100.Darroch, J.N. & Seneta, E. (1967). On quasi-stationary distributions in absorbing continuous-time finite markov chains.
Journal of Applied Probability , 4, 192–196.Diekmann, O., Heesterbeek, H. & Britton, T. (2012).
Mathematical tools for understandinginfectious disease dynamics . vol. 7. Princeton University Press.Emmert, K.E. & Allen, L.J. (2004). Population persistence and extinction in a discrete-time,stage-structured epidemic model.
Journal of Difference Equations and Applications , 10,1177–1199.Emmert, K.E. & Allen, L.J. (2006). Population extinction in deterministicand stochasticdiscrete-time epidemic models with periodic coefficients with applications to amphibianpopulations.
Natural Resource Modeling , 19, 117–164.Enatsu, Y., Nakata, Y. & Muroya, Y. (2010). Global stability for a class of discrete sirepidemic models.
Math. Biosci. Eng , 7, 347–361.Faure, M. & Schreiber, S.J. (2014). Quasi-stationary distributions for randomly perturbeddynamical systems.
The Annals of Applied Probability , 24, 553–598.Grimm, V. & Wissel, C. (2004). The intrinsic mean time to extinction: a unifying approachto analysing persistence and viability of populations.
Oikos , 105, 501–511.
Guckenheimer, J. & Holmes, P. (1983).
Nonlinear Osillations, Dynamical Systems, and Bi-furcations of Vector Fields . Springer-Verlag.Harier, M. (2018). Ergodic properties of Markov processes.Hofbauer, J. & So, J.W.H. (1989). Uniform persistence and repellors for maps.
Proc. Amer.Math. Soc. , 107, 1137–1142.Keeling, M.J. & Rohani, P. (2011).
Modeling infectious diseases in humans and animals .Princeton University Press.Klepac, P., Pomeroy, L.W., Bjørnstad, O.N., Kuiken, T., Osterhaus, A.D. & Rijks, J.M.(2009). Stage-structured transmission of phocine distemper virus in the dutch 2002 out-break.
Proceedings of the Royal Society B: Biological Sciences , 276, 2469–2476.Kong, J.D., Jin, C. & Wang, H. (2015). The inverse method for a childhood infectious diseasemodel with its application to pre-vaccination and post-vaccination measles data.
Bulletinof Mathematical Biology , 77, 2231–2263.Kryscio, R.J. & Lef`evre, C. (1989). On the extinction of the SIS stochastic logistic epidemic.
Journal of Applied Probability , 26, 685–694.Lindholm, M. & Britton, T. (2007). Endemic persistence or disease extinction: The effect ofseparation into sub-communities.
Theoretical Population Biology , 72, 253–263.Ma, X., Zhou, Y. & Cao, H. (2013). Global stability of the endemic equilibrium of a discretesir epidemic model.
Advances in Difference Equations , 2013, 42.N˚asell, I. (1996). The quasi-stationary distribution of the closed endemic sis model.
Advancesin Applied Probability , 28, 895–932.N˚asell, I. (1999). On the quasi-stationary distribution of the stochastic logistic epidemic.
Mathematical Biosciences , 156, 21–40.N˚asell, I. (2001). Extinction and quasi-stationarity in the Verhulst logistic model.
Journal ofTheoretical Biology , 211, 11–27.N˚asell, I. (2002). Stochastic models of some endemic infections.
Mathematical Biosciences ,179, 1–19.Pollicott, M., Wang, H. & Weiss, H. (2012). Extracting the time-dependent transmission ratefrom infection data via solution of an inverse ode problem.
Journal of Biological Dynamics ,6, 509–523.Satsuma, J., Willox, R., Ramani, A., Grammaticos, B. & Carstea, A. (2004). Extending thesir epidemic model.
Physica A: Statistical Mechanics and its Applications , 336, 369–375.Smith, H.L. & Waltman, P. (1999). Perturbation of a globally stable steady state.
Proc.Amer. Math. Soc. , 127, 447–453.Wang, Y. & Jiang, J. (2001). The general properties of discrete-time competitive dynamicalsystems.
J. Diff. Eqns. , 24, 470–493.Weiss, G.H. & Dishon, M. (1971). On the asymptotic behavior of the stochastic and deter-ministic models of an epidemic.