key: cord-0546069-csasum2q
authors: Schreiber, Sebastian J.; Huang, Shuo; Jiang, Jifa; Wang, Hao
title: Extinction and quasi-stationarity for discrete-time, endemic SIS and SIR models
date: 2020-05-17
journal: nan
DOI: nan
sha: 9ca49fc4168fc0b577d5a0c0cf842a82e0da0915
doc_id: 546069
cord_uid: csasum2q

Stochastic discrete-time SIS and SIR models of endemic diseases are introduced and analyzed. For the deterministic, mean-field model, the basic reproductive number $R_0$ determines their global dynamics. If $R_0le 1$, then the frequency of infected individuals asymptotically converges to zero. If $R_0>1$, then the infectious class uniformly persists for all time; conditions for a globally stable, endemic equilibrium are given. In contrast, the infection goes extinct in finite time with probability one in the stochastic models for all $R_0$ values. To understand the length of the transient prior to extinction as well as the behavior of the transients, the quasi-stationary distributions and the associated mean time to extinction are analyzed using large deviation methods. When $R_0>1$, these mean times to extinction are shown to increase exponentially with the population size $N$. Moreover, as $N$ approaches $infty$, the quasi-stationary distributions are supported by a compact set bounded away from extinction; sufficient conditions for convergence to a Dirac measure at the endemic equilibrium of the deterministic model are also given. In contrast, when $R_0<1$, the mean times to extinction are bounded above $1/(1-alpha)$ where $alpha<1$ is the geometric rate of decrease of the infection when rare; as $N$ approaches $infty$, the quasi-stationary distributions converge to a Dirac measure at the disease-free equilibrium for the deterministic model. For several special cases, explicit formulas for approximating the quasi-stationary distribution and the associated mean extinction are given. These formulas illustrate how for arbitrarily small $R_0$ values, the mean time to extinction can be arbitrarily large, and how for arbitrarily large $R_0$ values, the mean time to extinction can be arbitrarily large.

1. Introduction. Infectious disease modeling has been one of the most important topics in mathematical biology (Keeling & Rohani, 2011) . A recent Google scholar search 1 reveals over a million studies referencing SIS (susceptible-infected-susceptible) and SIR (susceptibleinfected-recovered) models. Most of these studies use deterministic, continuous-time equations. However, in seasonal systems or systems where measurements are taking at regular time 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 0 , determines whether a disease can persist or not in a population. If R 0 > 1, persistence often occurs. While if R 0 < 1, the disease-free state (i.e. extinction) often is globally stable and the infection is lost asymptotically as time marches to infinity.

When considering finite populations without external sources of infection, Markov chain models typically predict that the disease goes extinct in finite-time whether R 0 > 1 or < 1. To understand this puzzling difference between the asymptotic behaviors of the deterministic and stochastic models (Bartlett, 1966; Keeling & Rohani, 2011; Diekmann et al., 2012) , one can use the concept of quasi-stationarity that describes the long-term behavior of the stochastic model conditioned on non-extinction (Darroch & Seneta, 1965 , 1967 . For finite state models, the quasi-stationary distribution corresponds to a normalized left eigenvector π of the transition matrix of the Markov chain restricted to the transient states, i.e. the states where the infection persists. In discrete-time, if the disease dynamics follow the quasi-stationary distribution, then the eigenvalue λ associated with this eigenvector corresponds to the probability of 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 extinction 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 between the stochastic and deterministic models, it is natural to ask: How does the intrinsic mean time to extinction increase as the population size gets larger? How is the quasi-stationary distribution related to the asymptotic dynamics of the deterministic model as the population size gets larger? More generally, how do these quantities depend on the parameters such as R 0 ? For continuous-time, stochastic SIS models, there exists a dichotomy in the mean time to extinction when a fixed, positive fraction of the population is infected (Weiss & Dishon, 1971 ; Barbour, 1975; Kryscio & Lefèvre, 1989) . When R 0 > 1, this time increases exponentially with the population size N in the limit of large population sizes. When R 0 < 1, these extinction times remain bounded in the limit of large population size. However, to the best of our knowledge, similar statements for the intrinsic mean extinction times have not been rigorously proven for these continuous-time SIR models. However, in a series of papers (Nåsell, 1996 (Nåsell, , 1999 (Nåsell, , 2001 (Nåsell, , 2002 , Nåsell provided methods to approximate the intrinsic mean extinction times as well as the quasi-stationary distributions. His approximations support the existence of a similar dichotomy for intrinsic mean times to extinction. Moreover, they highlight a remarkable dichotomy about the qualitative behavior of the quasi-stationary distributions. When R 0 > 1, these distributions are well-approximated by a normal distribution centered near the endemic equilibrium. When R 0 < 1, the quasi-stationary distribution is best approximated by a discrete, geometric distribution. Despite these advances, mathematically rigorous results for discrete-time, stochastic SIS and SIR models are lacking.

In this paper, we introduce a new class of discrete-time SIS and SIR deterministic and stochastic models that have several desirable properties including (i) they are derived with individual-based rules and, consequently, preserve non-negativity of all populations, (ii) the deterministic models are the mean field model of the stochastic models, and (iii) the deterministic models converge to the classical continuous-time models in an appropriate limit. For these models, we analyze the global dynamics of deterministic models and then use this analysis in conjunction with large deviation results from Faure & Schreiber (2014) to rigorously characterize the behavior of the intrinsic mean times to extinction and quasi-stationary distributions in the limit of large population sizes. Moreover, for some special cases, we derive explicit approximations for the quasi-stationary distributions and extinction times that apply for all population sizes.

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 quasistationary 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 and intrinsic mean extinction times for the stochastic SIR model. Section 6 discusses our main findings 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 version of the classical susceptible-infected-susceptible (SIS) model where individuals are either susceptible 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 natural mortality with probability e −µ while infected individuals escape recovery with probability e −γ , where µ > 0 and γ > 0. Susceptible individuals from the previous time step who have not died, escape infection with probability e −βIn where β > 0 is the contact and transmission rate. If the population size remains constant, then the disease dynamics are given by

This discrete-time formulation of the SIS model has several advantages. First, it is straightforward to verify that the dynamics of I n remain in the interval [0, 1] provided the initial value I 0 lies in this interval. Second, these dynamics, as described in the next section, correspond to the mean field dynamics of an individual-based model. Finally, if ∆t is the length of a time step, and β =β∆t, γ =γ∆t, µ =μ∆t, then I(t + ∆t) := I n+1 = (1 − I(t))β∆t + I(t) − (μ +γ)∆tI(t) + O(∆t 2 ) where I(t) := I n Hence, in the limit ∆t → 0, we get the classical SIS ordinary differential equation

To understand the dynamics of (2.1), we can linearize at the origin to obtain the per-capita growth rate of the infection at the disease free-equilibrium

The basic reproduction, alternatively, is given by

As α > 1 if and only if R 0 > 1, we can use the basic reproductive number to characterize the global dynamics, as the following theorem shows.

Theorem 2.1. (i)If R 0 ≤ 1, then the origin is globally asymptotically stable.

(ii)If R 0 > 1, then there is a unique positive fixed point in (0, 1] such that it is globally asymptotically stable in (0, 1]. 3. Metastability and extinction in a stochastic SIS model. For the individual-based stochastic model, we require the additional parameter of the total population size N . Given this population size, the state space corresponds to the possible fractions of infected individuals in the population

Let I n ∈ S be the fraction of individuals infected at time step n. To determine the fraction infected in the next time step, we assume that each infected individual remains infected with probability e −µ−γ independent of each other and each susceptible individual lives and becomes infected with probability e −µ (1 − e −βIn ) independent of each other. Hence,

(3.1)

For all 1 ≤ i, j ≤ N , let Q ij = P[X n+1 = j/N |X n = i/N ] be the transition probabilities restricted to the transient states S \ {0} and Q = (Q ij ) be the associated N × N matrix. Unlike the deterministic model, the disease goes extinct in finite time with probability one for the stochastic model. The following proposition follows from standard results in stochastic processes (see e.g. Harier, 2018, Theorem 3.20) .

Proposition 3.1. Assume that µ + γ and β are positive. With probability one, I n = 0 for some n ≥ 1.

Even though extinction is inevitable, it may be preceded by long term transients. To characterize 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 π = (π 1 , . . . , π N ) (depending on N ) such that i π i = 1, π i > 0 for all i, and πQ = λπ. π is the quasi-stationary distribution which satisfies (Darroch & Seneta, 1965) lim

i.e. the probability of having j individuals in the long-term given the disease hasn't gone extinct equals π j . Furthermore,

i.e. λ is the probability the disease persists to the next time step given the process is following the quasi-stationary distribution. Hence, the mean time to extinction, when following the quasi-stationary distribution, is 1 1−λ , what Grimm & Wissel (2004) call the intrinsic mean time to extinction.

Our main result for the stochastic SIS model concerns the behavior of the quasi-stationary distribution and the intrinsic mean time to extinction for large population size N .

Theorem 3.2. Assume µ + γ > 0, β > 0. For each N ≥ 1, let π N be the quasi-stationary distribution and λ N the corresponding eigenvalue. Let α = βe −µ + e −µ−γ . Then (i) If α ≤ 1 (equivalently R 0 ≤ 1), then λ N ≤ α for all N ≥ 1 and lim N →∞ π N = δ 0 where δ 0 is a Dirac measure at 0 and convergence is in the weak* topology on probability measures on [0, 1]. (ii) If α > 1 (equivalently R 0 > 1), then lim N →∞ π N = δ I * where δ I * is the Dirac measure at the unique positive fixed point I * of equation ( 

, the numerically computed intrinsic mean time to extinction (x marks) and the analytical approximation (solid blue curve) 1/(1 − (1 − e −µ ) N ). Parameter values: γ = 0, µ = 1.5, β = 100, and N = 20 for (A).

The first assertion of Theorem 3.2 implies that if α < 1 and the population size is large, then any long-term transients mostly involve low frequencies of infected individuals and the mean time to extinction after these transients is less than 1 1−α . We conjecture that in the limit of large population size, N → ∞, the intrinsic mean time to extinction equals 1 1−α . The second assertion of Theorem 3.2 implies that if α > 1 and the population size is large, then the long-term transients fluctuate around the equilibrium frequency of the deterministic model and the mean time to extinction following these transients increases exponentially with population size, i.e. there exist c 1 , c 2 > 0 such that 1 1−λ N ≥ c 1 e c 2 N for all N ≥ 1. Figure 3 .1 illustrates these conclusions numerically.

Given that Theorem 3.2 describes the effect on increasing population size on the intrinsic mean time to extinction for a fixed value of α, it is natural to ask what effect increasing α has on 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 low recovery rates γ = 0 and very high β 1 contact rates. In the limit of β → ∞, the update rule for I n for I n > 0 is approximately I n+1 ∼ Binom(N, e −µ ). Namely, provided there is at least one infected individual at time step n, all individuals that have not died get infected. In this case, the quasi-stationary distribution is approximately

In particular, the mean intrinsic extinction time is bounded in the limit of α → ∞ due to β → ∞. The accuracy of this approximation is illustrated in Figure 3 .2. Second, in the limit of no recovery and no mortality (i.e. µ = γ = 0), the disease (not surprisingly!) never goes extinct whenever I 0 > 0. Indeed, in this case, I n → 1 as n → ∞ with probability one provided β > 0 and I 0 > 0. These two special cases highlight that the effect of α on the intrinsic mean time to extinction depends in a subtle way as α increases.

4. The dynamics of a deterministic SIR model. As discrete-time counterpart to the classical susceptible-infected-recovered (SIR) model, we assume all individuals escape natural 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 of infected individuals and β > 0 is the contact and transmission rate. If the population size is constant, then the discrete-time dynamics are given by

Like the discrete-time SIS model, this discrete-time formulation of the SIR model has several advantages. First, by adding the two equations of (4.1) together, we obtain that the trajectories of (4.1) remain in the domain

provided the initial value (S 0 , I 0 ) lies in this domain. Furthermore, if we define ∂X 0 := {(S, 0) : 0 ≤ S ≤ 1} and X 0 := X \ ∂X 0 , then X 0 and ∂X 0 are positively invariant. Second, these dynamics, as described in the next section, correspond to the mean field dynamics of an individual-based model. Finally, if ∆t is the length of a time step, and β =β∆t, γ =γ∆t, µ =μ∆t, then

Hence, in the limit ∆t → 0, we get the classical SIR system of ordinary differential equations

For our discrete-time SIR model, the disease-free fixed point is (1, 0). At this fixed point, the per-capita growth rate of the disease still equals α = βe −µ + e −µ−γ and the reproductive number still equals R 0 = βe −µ /(1 − e −µ−γ . We will show that if R 0 > 1, the disease persists and if R 0 ≤ 1, the disease-free equilibrium is globally stable. Furthermore, we will show when the recovery rate γ is sufficiently small, there is a globally stable endemic equilibrium. To state these results precisely, we define the parameter space as P := {λ := (µ, β, γ) : µ > 0, β > 0, γ ≥ 0}. Let C 0 P := {λ = (µ, β, 0) ∈ P : α(λ) > 1} be the parameters corresponding to no recovery (γ = 0) and α > 1. Finally, define C P := {λ ∈ P : α(λ) > 1, (4.1) admits a globally stable hyperbolic fixed point in X 0 }.

Using these definitions, we have the following theorem.

Theorem 4.1. (i) If α ≤ 1, then the disease free fixed point (1, 0) is globally asymptotically stable.

(ii) If α > 1, then F : X 0 → X 0 admits a global and compact attractor contained in the interior of X 0 .

(iii) C P ⊃ C 0 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 > 0 is sufficiently small, and α > 1 (equivalently, R 0 > 1). 5. Metastability and extinction in a stochastic SIR model. As with the stochastic SIS model, the stochastic SIR model requires the additional parameter of the total population size N . For a given population size, the state space S corresponds to the possible fraction of susceptible and infected individuals in the population, i.e.

Let (S n , I n ) ∈ S be the fractions of susceptible and infected individuals at time step n. The fraction of removed individuals at time step n equals R n = 1 − S n − I n . Consistent with the deterministic SIR model, we assume (i) each susceptible individual lives and becomes infected with probability e −µ (1 − e −βIn ) independent of each other, (ii) each infected individual either remains infected, dies and gets replaced with a susceptible individual, or enters the removed class 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 probability of success e −µ (1 − e −βIn ) (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 binomial random variable with N I n −X n+1 trials with probability of success e −γ (i.e. non-dying infected individuals that will not enter the removed class), and Z n+1 be a binomial random variable with N (1 − I n − S n ) trials with probability of success 1 − e −µ (i.e. removed individuals that die and get replaced with a susceptible individual). Under these assumptions, the stochastic SIR model is

As with the stochastic SIS model, the disease goes extinct in finite time with probability one for the stochastic model. The following proposition follows from standard results in stochastic processes (see e.g. Harier, 2018, Theorem 3.20) .

Proposition 5.1. Assume that µ + γ and β are positive. With probability one, I n = 0 for some n ≥ 1.

To characterize metastability and extinction times, define S + = {(x 1 , x 2 ) ∈ S : x 2 > 0} 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 transient states and Q = (Q xy ) x,y∈S + be the associated matrix. Let π = (π x ) x∈S + be the quasistationary distribution with associated persistence probability λ i.e. x∈S + π x = 1, π x > 0 for all x ∈ S and πQ = λπ. Our main result for the stochastic SIR model concerns the behavior of the quasi-stationary distribution and the intrinsic mean time to extinction for large population size N .

Theorem 5.2. Assume µ + γ > 0, β > 0. For each N ≥ 1, let π N be the quasi-stationary distribution and λ N the corresponding eigenvalue for (5.1). Let α = βe −µ + e −µ−γ . Then (i) If α < 1, then λ N ≤ α for all N ≥ 1 and lim N →∞ π N = δ (1,0) where δ (1,0) is a Dirac measure at disease-free equilibrium (1, 0) and convergence is in the weak* topology. (ii) If α > 1, then lim sup N →∞ 1 N log(1−λ N ) < 0 and there exists a compact set K ⊂ (0, 1) 2 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 fixed point (S * , I * ) of equation (4.1).

The first assertion of Theorem 5.2 implies that if α < 1 and the population size is large, then any long-term transient mostly involves low frequencies of infected individuals and the mean time to extinction after these transients is less than 1 1−α . Furthermore, whenever permanence of the deterministic model corresponds to a globally stable equilibrium (see Theorem 4.1), the QSDs concentrate on a Dirac measure at this equilibrium i.e. it supports the only invariant measure in K for the deterministic dynamics. The second assertion of Theorem 5.2 implies that if α > 1 and the population size is large, then the long-term transients fluctuate away from the disease-free equilibrium and the mean time to extinction following these transients increases exponentially with population size i.e. there exist c 1 , c 2 > 0 such that 1 1−λ N ≥ c 1 e c 2 N for all N ≥ 1. Figure 5 .1 illustrates these conclusions.

6. Discussion. This paper formulates and provides a mathematically rigorous analysis of deterministic and stochastic, discrete-time SIS and SIR models. The stochastic models are based on probabilistic, individual-based update rules. The conditional expected change in the fraction of infected and susceptible individuals given the current values of these fractions determines the update rule for the deterministic model and, in this sense, the deterministic model is the mean field model for the stochastic models.

Many earlier discrete-time epidemic models of SIS and SIR dynamics have been derived using numerical approximation schemes for differential equations (Allen, 1994; Ma et al., 2013; Satsuma et al., 2004; Enatsu et al., 2010; Castillo-Chavez & Yakubu, 2001) . These models, 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 escape function, is most similar to higher dimensional, discrete-time epidemiological models that have been used for applications to specific diseases (Emmert & Allen, 2004 , 2006 Allen & van den Driessche, 2008) . Unlike the models based on numerical approximation schemes, our analysis and numerical explorations suggest that our models always approach a global stable equilibrium. When R 0 ≤ 1, we prove that all trajectories asymptotically approach the diseasefree equilibrium for both the SIS and SIR models. When R 0 > 1, we prove the disease persists for both models, approaches a globally stable, endemic equilibrium for the SIS model, and provide 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 is globally stable whenever R 0 > 1. Hence, we conjecture that R 0 > 1 always implies global stability of the endemic equilibrium for the SIR model. Unlike the deterministic model for which the disease persists indefinitely when R 0 > 1 and only goes asymptotically extinct over an infinite time horizon when R 0 ≤ 1, the fraction of infected in the stochastic model becomes zero in finite time for all values of R 0 . To understand this well-know discrepancy (Bartlett, 1966; Keeling & Rohani, 2011; Diekmann et al., 2012) for our model, we characterized the long-term transients using quasi-stationary distributions for finite-state, discrete-time Markov chains (Darroch & Seneta, 1965) . For these characterizations, we used the per-capita growth rate α = βe −µ + e −µ−γ of the infection at the disease free equilibrium. When α < 1 (equivalently, R 0 < 1), the mean time to extinction, when following the quasi-stationary distribution, is ≤ 1/(1 − α) for all population sizes and for both the SIS and SIR models. Indeed, we conjecture that as N → ∞, this mean time to extinction converges to 1/(1 − α). While R 0 < 1 if and only if α < 1, we have α = (1 − e −µ−γ )R 0 + e −µ−γ > R 0 whenever R 0 < 1. Hence, even if R 0 is very small, the mean times to extinction can be arbitrarily long. For example, given any 0 < x < y < 1, we can make R 0 = x and α = y by choosing γ = 0, e −µ = y/2, and β = x(1 − e −µ ).

When R 0 > 1 (equivalently α > 1), we show that the mean extinction times increase exponentially with the population size N and the quasi-stationary distributions concentrate on positive invariant sets for the deterministic model for large N . In particular, coupled with our analysis of the deterministic dynamics, our results imply that the quasi-stationary behavior for large N always concentrates near the globally stable, endemic equilibrium of the 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 consistent with earlier studies of continuous-time Markov models where the analysis was done using diffusion approximations of the individual-based models (Barbour, 1975; Kryscio & Lefèvre, 1989; Nåsell, 1996 Nåsell, , 1999 Andersson & Britton, 2000; Nåsell, 2002; Lindholm & Britton, 2007; Andersson & Lindenstrand, 2011; Clancy & Tjia, 2018) . In contrast, our results apply large deviation methods from (Faure & Schreiber, 2014) directly to the individual-based models. An open question for the stochastic model with R 0 > 1 concerns the asymptotic rate at which the extinction times increase exponentially with N . Specifically, what is the value of α such that the mean time to extinction grows like exp(αN ) for large N ? The diffusion approximations provide one approach to find potential candidates for α.

When R 0 > 1, we found that the mean time to extinction can be arbitrarily large even for a fixed population size. For example, this occurs when recovery and mortality rates approach zero in which case R 0 also increases without bound but α remains bounded above by β + 1. In contrast, increasing contact rates (which increase α and R 0 without bound) leads to extinction times that are constrained by population size, recovery rates and mortality rates.

In addition to the open mathematical questions that we have already raised, future challenges include analyzing extensions of our models. These extensions could include additional compartments such as SEIR models, multi-age group epidemic models, and SIR type models with vaccination (Anderson & May, 1991; Keeling & Rohani, 2011; Kong et al., 2015) . More generally, when the discrete-time system is autonomous, the mathematical approaches used here should be applicable to study quasi-stationary distributions and the intrinsic extinction times. However, when population sizes or transmission rate change stochastically over time (Anderson & May, 1979; Pollicott et al., 2012) , new mathematical approaches are required for studying the impact of these environmentally driven random fluctuations on intrinsic extinction times.

7. Proof of Theorem 2.1. Suppose that α < 1. Then we shall prove that It is easy to see that (7.1) is equivalent to Therefore, g (I) < g (0) = 0, I ∈ (0, 1].

Furthermore, g(I) < g(0) = 0, I ∈ (0, 1].

This shows that (7.2) holds. Fix I 0 ∈ [0, 1] and set I n := F n (I 0 ), J n := L n (I 0 ) = α n I 0 , n = 1, 2, · · · . We claim that (7.3) I n ≤ J n , n = 1, 2, · · · .

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

By mathematical induction, (7.3) is valid for all positive integers. Since α < 1, L n (I 0 ) = α n I 0 → 0 as 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]. Using Feigenbaum's method given in (ii), we can easily prove that 0 is still globally asymptotically stable in this case.

(ii) Suppose that α > 1. Then it is not difficult to prove that F (I) has a unique positive fixed 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 trajectories converge to the fixed point I * .

(iia) Take I 0 ∈ (0, I * ). As depicted in Figure 7 .1(a), the iterating sequence I n increasingly tends to I * . If I 0 ∈ (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 0 > I * c , then the first iteration I 1 ∈ (0, I * ). After the first iteration, 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 (I) ≤ 0, I ∈ (I * c , 1]. By computation,

Therefore, −1 < F (I) ≤ 0, I ∈ (I * c , 1]. First, we prove that I 0 < I 2 < I * if I * c < I 0 < I * . Let K be a line through (I * , F (I * )) whose gradient is −1. Because we have proved F (I) > −1, 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 0 , 0). This line intersects 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 strictly greater than the ordinate of E, the abscissa of C is strictly greater than the abscissa of G.

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 , 1), the ordinate of 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 ordinate of L is strictly greater than the ordinate of A. Then the abscissa of L is strictly greater than the abscissa of A, that is, I 0 ≤ I 2 ≤ I * as shown Figure 7.2(b) . Now, we want to show lim n→∞ F 2n (I 0 ) = I * by contradiction. Suppose it is not true. Then lim n→∞ F 2n (I 0 ) = I * * < I * . Thus we set I * * to be the initial point and repeat the above process, then we have F 2 (I * * ) > I * . This is contradictory to lim n→∞ F 2n (I 0 ) = I * * , therefore, lim n→∞ F 2n (I 0 ) = I * . Similarly, we can prove lim n→∞ F 2n+1 (I 0 ) = I * . This concludes that I n oscillates around I * and converges to I * .

The proof of the case I 0 < I * c is given in the Figure 7 .2(c). 8. Proof of Theorem 3.2. We use results from Faure & Schreiber (2014) to prove Theorem 3.2. We begin by verifying that Standing Hypothesis 1.1 of Faure & Schreiber (2014) holds. In their notation, "ε" corresponds to 1 N in our models i.e. small demographic noise corresponds to large population size N . For all δ > 0 and N ≥ 1, define

where F (x) = e −γ−µ x + e −µ (1 − e −βx )(1 − x) corresponds to the right hand side of the deterministic model in equation (2.1). Standing Hypothesis 1.1 of Faure & Schreiber (2014) requires that lim N →∞ β δ (N ) = 0 for all δ > 0. The following proposition proves something stronger using large deviation estimates.

Proposition 8.1. There exists a function ρ : (0, ∞) → (0, ∞) such that

for all N ≥ 1 and δ > 0.

Proof. While we use standard large deviation estimates, we go through the details to ensure that the estimates can be taken uniformly in x ∈ [0, 1]. Define a = e −µ−γ and b(x) = e −µ (1 − e −βx ). By the exponential Markov inequality, we have for all t

Taking log of equation (8.1) and dividing by N yields

for all t = 0, x ∈ [0, 1], and δ > 0. Similarly, one can show that

for all t = 0, x ∈ [0, 1], and δ > 0. We have

As ψ(0, x) = 0 = ∂ψ ∂t (0, x) = 0, and ψ is strictly convex in t, for all δ > 0, there exists t * (δ) > 0 such that δt * (δ) > ψ(−t * (δ), x) and δt * (δ) > ψ(t * (δ), x) for all x ∈ [0, 1]. Define

for all x ∈ [0, 1] and δ > 0.

To prove the first result of Theorem 3.2, assume that α ≤ 1. Theorem 2.1 implies that 0 is globally stable for the deterministic model I → F (I). Theorem 3.12 of Faure & Schreiber (2014) , which only requires Standing Hypothesis 1.1, implies that lim N →∞ π N = δ 0 . Define R(x) = F (x)/x for x ∈ (0, 1]. Equation (7.1) in the proof of Theorem 2.1 implies that

To prove the second result of Theorem 3.2, assume α > 1 in which case Theorem 2.1 implies that there exists a unique positive globally stable equilibrium I * for the map I → F (I). Assertion (b) of Lemma 3.9 of Faure & Schreiber (2014) , implies that there exists δ > 0 such that 1 − λ N ≤ β δ (N ) for all N ≥ 1. Proposition 8.1 implies that

To complete the proof of the second assertion, we need to verify the assumption in Assertion (b') of Lemma 3.9 of Faure & Schreiber (2014) . Choose η > 0 sufficiently small so that min x∈ [0,η] (1 − e −µ−γ ) x ≥ exp(−ρ(δ)/3) and min x∈[0,η] (1 − e −µ (1 − e −βx )) 1−x ≥ exp(−ρ(δ)/3).

Then min

which verifies the assumption of (b') of Lemma 3.9 of Faure & Schreiber (2014) and implies that lim

i.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 → F (x). As these weak* limit points are supported on [η, 1] and the only invariant measure for x → F (x) on this interval is the Dirac measure δ I * and the unique positive fixed point, it follows that π N converges in the weak* topology to δ I * as claimed.

9. Proof of Theorem 4.1. Let E f = (1, 0). Denote by

Then

DF

It follows that the per-capita growth rate of the disease is still given by (2.2) and the basic reproduction number of (4.1) is still the expression in (2.3). Define the competitive cone

Then it is easy to check that DF (E f ) keeps K invariant (see Wang & Jiang (2001) ). Define

We shall verify that (9.1) F (S, I) ≤ K L(S, I), (S, I) ∈ X.

(9.1) is equivalent to

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, the competitive ordering relation (9.1) holds. Let P 0 := (S 0 , I 0 ) ∈ X, P n := F n (P 0 ) = (S n (P 0 ), I n (P 0 )), Q n := L n (P 0 ). We shall show that

The left inequality is obvious by the definition of competitive order and X. So we will prove the right one. (9.1) deduces that (9.4) holds for n = 1. Suppose that (9.4) holds for n. Then using (9.1) and the order-preserving for DF (E f ) , we obtain

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 with respect to (X 0 , ∂X 0 ), that is, there exists η > 0 such that

It is easy to see that E f is the maximal compact invariant set in ∂X 0 which is positively invariant with respect to F and lies on the stable manifold of E f . Recalling the Hofbauer and So uniform persistence theorem of Hofbauer & So (1989) , the system (4.1) is uniformly persistent if and only if (a) E f is isolated in X, and (b) W s (E f ) ⊂ ∂X 0 . Since the disease free fixed point E f is hyperbolic and a saddle, (a) and (b) follows immediately from the Hartman and Grobman theorem (see Guckenheimer & Holmes (1983) ). This verifies the uniform persistence and hence the system (4.1) admits an attractor in X 0 .

It follows from (9.5) that (4.1) contains a compact attractor A 0 ⊂ {(S, I) ∈ X : I ≥ η}. Besides, by the first equality of (4.1), we have S n (P 0 ) ≥ 1 − e −µ for n ≥ 1 and P 0 ∈ X. This implies that A 0 ⊂ {(S, I) ∈ X : S ≥ 1 − e −µ }. As a result,

This proves (ii).

It remains to prove (iii). First, we consider the system (4.1) with λ 0 = (µ 0 , β 0 , 0) and α(µ 0 , β 0 , 0) > 1. We shall prove that it admits a globally stable fixed point (1 − I * , I * ) in X, where F (I * ) = I * with 0 < I * < 1.

Let ∆ n := S n + I n . Then from (4.1) it follows that (9.6) ∆ n+1 = 1 − e −µ 0 + e −µ 0 ∆ 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 , β = β 0 , γ = 0. Applying Theorem 2.1(ii), we get that the system (2.1) has a globally stable fixed point I * in (0, 1). Thus the system (4.1) admits a globally stable fixed point (1 − 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 Figure  7 .2(a)), we have (9.7) |F (I * )| = e −µ 0 −β 0 I * 1 + β 0 (1 − I * ) < 1.

Next, we will prove the spectral radius of the Jacobian matrix for F (S, I) at the positive fixed point E * (S * , I * ) := (1 − I * , I * ) is strictly less than 1.

An easy calculation yields that

and detDF (E * ) = e −µ 0 F (I * ). This proves that F (I * ) and e −µ 0 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 · the Euclidean norm of R 3 and B C (λ 0 , s) the open ball in P of radius s about the point λ 0 . We fix a λ 0 ∈ C P and an s 0 ∈ (0, µ 0 ) such that α(λ) > 1 for any λ ∈ B C (λ 0 , s 0 ). In order to consider the perturbed systems for parameters, we set F λ 0 (S, I) the given mapping and F λ (S, I) = (S λ (S, I), I λ (S, I)) the mappings for λ ∈ B C (λ 0 , s 0 ). Define

Then R λ (S, I) is continuous on X. By induction, it not difficult to prove that (9.8) F n λ (S, 0) = 1 − e −nµ + e −nµ S, 0 R λ (F n λ (S, 0)) =βe −µ 1 − e −nµ + e −nµ S + e −µ−γ for n = 1, 2, · · · . We claim that there exist an s 1 ∈ (0, s 0 ], an integer N > 0, δ > 0 and ρ > 1, all only depending on λ 0 , such that (9.9) I N λ (S, I) ≥ ρI for all λ ∈ B C (λ 0 , s 1 ) and I ∈ [0, δ],

where F n λ (S, I) = (S n λ (S, I), I n λ (S, I)). From (9.8) it follows that

By the continuity of α(λ), there exists an s 1 ∈ (0, s 0 ] such that α(λ) > α(λ 0 )+1 2 for all λ ∈ B C (λ 0 , s 1 ), and hence

This implies that there is an integer N > 0, only depending on λ 0 , such that (9.10) R λ (F N λ (S, 0)) > α(λ 0 ) + 3 4 for all λ ∈ B C (λ 0 , s 1 ).

Using (9.10) and the uniform continuity of R λ (F N λ (S, I)) on B C (λ 0 , s 1 ) × X, we obtain that there is a δ > 0, only depending on λ 0 , such that R λ (F N λ (S, I)) > α(λ 0 ) + 7 8 := ρ for all λ ∈ B C (λ 0 , s 1 ) and 0 ≤ I ≤ δ, which implies that (9.9) holds, thus the claim is proved. By (9.9), we get that for each λ ∈ B C (λ 0 , s 1 ),

This shows that there exists at least a positive integer m with the property δ, 1] . Then all assumptions of Theorem 2.1 of Smith & Waltman (1999) have been checked. It follows that B C (λ 0 , s 1 ) ⊂ C P . The proof is complete.

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.1 of Faure & Schreiber (2014) holds. For all δ > 0 and N ≥ 1, define β δ (N ) = sup

x,y≥0,x+y≤1 P [ (S n+1 , I n+1 − F (x, y) ∞ ≥ δ|(S n , I n ) = (x, y)]

where F (x, y) = (1 − e −µ + e −µ−βy x, xe −µ (1 − e −βy ) + ye −µ−γ ) corresponds to the right hand side of the deterministic model in equation (4.1) and (x, y) ∞ = max{|x|, |y|} corresponds to the sup norm. Standing Hypothesis 1.1 of (Faure & Schreiber, 2014) requires that lim N →∞ β δ (N ) = 0 for all δ > 0. Like Proposition 8.1 for the stochastic SIS model, the following proposition proves something stronger using large deviation estimates.

Proposition 10.1. There exists a function ρ : (0, ∞) → (0, ∞) such that

for all N ≥ 1 and δ > 0.

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 probability of success a 1 (y) = 1 − e −µ (1 − e −βy ) and X n+1 + Z n+1 has N (1 − x) trials with probability of success b 1 = 1 − e −µ . Using the exponential Markov inequality as in the proof of Proposition 8.1, we get (10.1) 1 N log P[N (S n+1 − F 1 (x, y)) ≥N δ|(S n , I n ) = (x, y)] ≤ −δt + ψ 1 (t, x, y)

1 N log P[N (F 1 (x, y) − S n+1 ) ≥N δ|(S n , I n ) = (x, y)] ≤ −δt + ψ 1 (−t, x, y)

for all t, δ and where ψ 1 (t, x, y) = −tF 1 (x, y) + x log(1 − a 1 (y) + a 1 (y)e t ) + (1 − x) log(1 − b 1 + b 1 e t ).

As ψ 1 (0, x, y) = 0 = ∂ψ ∂t (0, x, y) = 0, and ψ 1 is strictly convex in t, for all δ > 0, there exists t * (δ) > 0 such that δt * (δ) > ψ 1 (−t * (δ), x, y) and δt * (δ) > ψ 1 (t * (δ), x, y) for all x, y ∈ [0, 1] such that x + y ≤ 1. Define ρ 1 (δ) = δt * (δ) − max

x,y∈[0,1],x+y≤1 ψ 1 (t * (δ), x, y) > 0.

Equation (10.1) implies that P[|S n+1 − F 1 (x, y)| ≥ δ|(S n , I n ) = (x, y)] ≤ exp(−N ρ 1 (δ))

for all x, y ∈ [0, 1] with x + y ≤ 1 and δ > 0. W n+1 and Y n+1 conditioned on (S n , I n ) = (x, y) in equation (5.1) are also independent binomial 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 exponential Markov inequality, one can show there exists a function ρ 2 : (0, ∞) → (0, ∞) such that P[|I n+1 − F 2 (x, y)| ≥ δ|(S n , I n ) = (x, y)] ≤ exp(−N ρ 2 (δ)) for all x, y ∈ [0, 1] with x + y ≤ 1 and δ > 0. Setting ρ(δ) = min{ρ 1 (δ), ρ 2 (δ)} completes the proof 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) → F (S, I). Theorem 3.12 of Faure & Schreiber (2014) , which only requires Standing Hypothesis 1.1, implies that lim N →∞ π N = δ (1,0) in the weak* topology. Define R(x, y) = F 2 (x, y)/y for y ∈ (0, 1], x ∈ [0, 1], and x+y ≤ 1. Equation (9.2) in the proof of Theorem 4.1 implies that R(x, y) ≤ α. For N ≥ 1, quasistationarity of π N implies λ N x,y∈S + yπ N x,y =

x,y∈S + E I n+1 (S n , I n ) = (x, y) π N x,y = x,y∈S + F 2 (x, y) π N x,y = x,y∈S + yR (x, y) π N x,y ≤α x,y∈S + yπ N

x,y Since x,y∈S + yπ N x,y > 0, λ N ≤ α for all N ≥ 1 as claimed. To prove the second result of Theorem 5.2, assume α > 1 in which case Theorem 4.1 implies that there exists a global, compact attractor K ⊂ (0, 1) × (0, 1) for (S, I) → F (S, I). Assertion (b) of Lemma 3.9 of Faure & Schreiber (2014) , implies that there exists δ > 0 such that 1 − λ N ≤ β δ (N ) for all N ≥ 1. Proposition 10.1 implies that

Assumption in Assertion (b') of Lemma 3.9 of Faure & Schreiber (2014) holds for an argument similar 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) → F (x, y).

Two applications of urn processes the fringe analysis of search trees and the simulation of quasi-stationary distributions of Markov chains

The basic reproduction number in some discrete-time epidemic models

Some discrete-time SI, SIR, and SIS epidemic models

Comparison of deterministic and stochastic SIS and SIR models in discrete time

Infectious diseases of humans: dynamics and control

Population biology of infectious diseases: Part i

Stochastic epidemics in dynamic populations: quasistationarity and extinction

A stochastic sis epidemic with demography: initial stages and time to extinction

The duration of the closed stochastic epidemic

An Introduction to Stochastic Processes With Special Reference to Methods and Applications

Discrete-time sis models with complex dynamics

Approximating time to extinction for endemic infection models

On quasi-stationary distributions in absorbing discretetime finite markov chains

On quasi-stationary distributions in absorbing continuoustime finite markov chains

Mathematical tools for understanding infectious disease dynamics

Population persistence and extinction in a discrete-time, stage-structured epidemic model

Population extinction in deterministicand stochastic discrete-time epidemic models with periodic coefficients with applications to amphibian populations

Global stability for a class of discrete sir epidemic models

Quasi-stationary distributions for randomly perturbed dynamical systems

The intrinsic mean time to extinction: a unifying approach to analysing persistence and viability of populations

Nonlinear Osillations, Dynamical Systems, and Bifurcations of Vector Fields

Ergodic properties of Markov processes

Uniform persistence and repellors for maps

Modeling infectious diseases in humans and animals

Stage-structured transmission of phocine distemper virus in the dutch 2002 outbreak

The inverse method for a childhood infectious disease model with its application to pre-vaccination and post-vaccination measles data

On the extinction of the SIS stochastic logistic epidemic

Endemic persistence or disease extinction: The effect of separation into sub-communities

Global stability of the endemic equilibrium of a discrete sir epidemic model

The quasi-stationary distribution of the closed endemic sis model

On the quasi-stationary distribution of the stochastic logistic epidemic

Extinction and quasi-stationarity in the Verhulst logistic model

Stochastic models of some endemic infections

Extracting the time-dependent transmission rate from infection data via solution of an inverse ode problem

Extending the sir epidemic model

Perturbation of a globally stable steady state

The general properties of discrete-time competitive dynamical systems

On the asymptotic behavior of the stochastic and deterministic models of an epidemic