J. Math. Biol. DOI 10.1007/s00285-014-0828-1

Mathematical Biology

On the stochastic SIS epidemic model in a periodic environment Nicolas Bacaër

Received: 11 March 2014 / Revised: 11 July 2014 © Springer-Verlag Berlin Heidelberg 2014

Abstract In the stochastic SIS epidemic model with a contact rate a, a recovery rate b < a, and a population size N , the mean extinction time τ is such that (log τ )/N converges to c = b/a − 1 − log(b/a) as N grows to infinity. This article considers the more realistic case where the contact rate a(t) is a periodic function whose average is bigger than b. Then (log τ )/N converges to a new limit C, which is linked to a time-periodic Hamilton–Jacobi equation. When a(t) is a cosine function with small amplitude or high (resp. low) frequency, approximate formulas for C can be obtained analytically following the method used in Assaf et al. (Phys Rev E 78:041123, 2008). These results are illustrated by numerical simulations. Keywords Hamilton–Jacobi equation · Epidemic model · Extinction · Periodic environment Mathematics Subject Classification

35F21 · 60J80 · 92D30

1 Introduction The stochastic SIS epidemic model has been investigated in much detail when the environment is assumed constant, as in the book by Nåsell (2011). With a contact rate a, a recovery rate b < a, and a population size N , the mean extinction time τ (starting for example from one infected person) is such that

N. Bacaër (B) IRD (Institut de Recherche pour le Développement), UMMISCO, Bondy, France e-mail: [email protected] N. Bacaër Université Paris 6, UMMISCO, Campus des Cordeliers, Paris, France

123

N. Bacaër

log τ −→ c = b/a − 1 − log(b/a) > 0 N N →+∞

(1)

(Nåsell 2011, Theorem 12.1). Among other references, see also Andersson and Djehiche (1998), Doering et al. (2005), and Assaf and Meerson (2010). The latter use the method of Wentzel, Kramers and Brillouin [WKB]. The probability Pn (t) of having n ≥ 1 infected people at time t first approaches a quasi-stationary distribution πn . Setting x = n/N so that 0 ≤ x ≤ 1 and letting N go to infinity, −(log πn )/N approaches a continuous function S(x) satisfying the stationary Hamilton–Jacobi  equation H x, ∂∂ xS = 0 with H (x, p) = ax(1 − x)(e p − 1) + bx(e− p − 1) = x(1 − e− p )[a(1 − x)e p − b]

(2)

[Assaf and Meerson 2010, Eq. (12) and §IV.D.3]. More precisely the branch of the level set H = 0 given by a(1 − x)e p − b = 0 leads to S(x) = x log(b/a) + x + (1 − x) log(1 − x) + constant.

(3)

This function has a minimum when x = x ∗ = 1 − b/a, with x ∗ > 0 since b < a. Finally c = S(0) − S(x ∗ ) is the height between the bottom and the border at x = 0 of the potential well S(x). Equivalently the Hamiltonian system dx ∂H = , dt ∂p

dp ∂H =− dt ∂x

(4)

has a heteroclinic orbit connecting (x ∗ , 0) = (1 − ab , 0) to (0, p ∗ ) = (0, log ab ), and 0 c is equal to the action x ∗ p d x along this orbit. This WKB method has been used for a number of other birth-and-death processes from physics or population biology (Ovaskainen and Meerson 2010; Kamenev 2011). For a chemical reaction involving branching and annihilation, Escudero and Rodríguez (2008) studied how a time-periodic environment influenced the heteroclinic orbit that plays a central role in the large N behavior of the mean extinction time. Assaf et al. (2008) investigated the same model more in detail, computing in particular the correction to the mean extinction time due to a periodic perturbation with small amplitude or high (resp. low) frequency. The authors obtained general formulas that can be applied to other birth-and-death processes exhibiting metastability. Some Monte–Carlo simulations of a stochastic SIS model with periodic treatment rate are shown in Billings et al. (2013, Fig. 7). Our goal here is to apply the WKB method used by Assaf et al. (2008) to the SIS epidemic model with a T -periodic contact rate a(t) whose average is bigger than b. Such a model could mimic for example the spread of a bacterial infection that does not confer any immunity in a school with a weekly periodicity due to weekends, or a yearly periodicity due to holidays and seasonality. Of course it is just a first step towards more realistic models.

123

Stochastic SIS epidemic model

In Sect. 2 informal computations suggest that the mean extinction time τ (starting for example from one infected person at time 0) is such that log τ −→ C = min S ∗ (t, 0+ ) − min min S ∗ (t, x). 0≤t≤T 0≤t≤T 0≤x≤1 N N →+∞

(5)

Here S ∗ (t, x) is any T -periodic viscosity solution of the Hamilton–Jacobi equation   ∂S ∂S + H t, x, =0 ∂t ∂x

(6)

for 0 < x < 1 with the mixed Dirichlet-“state constraint” boundary conditions S(t, 0) = 0,

∂S (t, 1) = +∞, ∂x

such that S(t, x) is not identically zero near x = 0. The boundary conditions should be understood in a viscosity sense (Barles 1994) since, e.g., S ∗ (t, 0+ ) may not be equal to 0. The Hamiltonian H (t, x, p) is the same as (2) but with the constant parameter a replaced by a(t). When a(t) = a0 (1 + ε cos(ωt)) with ω = 2π/T, a0 > b, and |ε| ≤ 1, one can set c0 = b/a0 − 1 − log(b/a0 ). Following the methods of Assaf et al. (2008), Sect. 2 shows that C  c0 −

π ω |ε|   ω a0 sinh aπ0 −b

when ε is close to 0, C  c0 − |ε|(1 − b/a0 ) when ω  a0 , and C  c0 −

(a0 − b)2 ε2 (1 + 2b/a0 ) 12 ω2

in the high-frequency limit ω a0 . One can conjecture that C > 0 as long as  1 T T 0 a(t) dt > b. One can also conjecture that C is always less than c0 , as one might expect since seasonal variations are likely to favor disease extinction. This suggests more precisely that a periodic environment leads to an exponential decrease of the mean extinction time for the SIS model similar to that found for the branching and annihilation model (Assaf et al. 2008). Section 3 illustrates these results with numerical simulations. Section 4 adds a few remarks. 2 Analytical computations 2.1 The Hamilton–Jacobi partial differential equation Master equation and Floquet theory. Recall that N is the total population, which is assumed constant. There are S(t) susceptible and I (t) infected people at time t, with

123

N. Bacaër

N = S(t) + I (t). If I (t) = n, then the probability that one person recovers between t and t + dt is b n dt + o(dt). The probability for one new infection to occur is a(t)n(1 − n/N )dt + o(dt). So a(t) is the contact rate and b the recovery rate. Assume that a(t) is a positive T -periodic continuous function and that 1 r= T



T

a(t) dt − b > 0 or equivalently R0 =

0

1 T

T 0

a(t) dt > 1. b

For a biological interpretation of R0 , see Bacaër and Ait Dads (2012). Hethcote (1973) noticed also that r > 0 (i.e. R0 > 1) is a necessary and sufficient condition for the solution of the mean field equation di/dt = a(t)i(1 −i) − b i to converge to a positive periodic function. Otherwise the solution converges to zero. Let Pn (t) be the probability that I (t) = n. The master equation is d Pn = a(t)(n − 1)(1 − (n − 1)/N )Pn−1 − [a(t)n(1 − n/N ) + b n]Pn dt + b(n + 1)Pn+1

(7)

N for 0 ≤ n ≤ N if one sets P−1 = 0 and PN +1 = 0. Of course n=0 Pn (t) = 1. System (7) may be written as d P/dt = M(t) P, where P(t) is the vector (Pn (t))0≤n≤N and M(t) is the square matrix of size N + 1 ⎛

0 b 0 ⎜ 0 −b − a(t)(1 − 1 ) 2b N ⎜ ⎜ −2b − 2a(t)(1 − a(t)(1 − N1 ) M(t) = ⎜ 0 ⎜. . .. .. ⎝ .. . 0

0

0

0 ··· 0 ··· 2 ) 3b ··· N .. .

0 0 0 .. .

⎞ ⎟ ⎟ ⎟ ⎟. ⎟ ⎠

0 · · · −bN

 0 ∗ , where Q(t) is a square matrix 0 Q(t) of size N . Let X (t) (resp. Y (t)) be the matrix solution of the system d X/dt = M(t) X (resp. dY/dt = Q(t) Y ) with X (0) = I N +1 (resp. Y (0) = I N ), I N being the identity matrix of size N . The set of Floquet multipliers of M(t), i.e., the eigenvalues of the monodromy matrix X (T ), is the union of {μ0 = 1} and of the Floquet multipliers of Q(t). The matrix Q(t) is cooperative: the off-diagonal elements are non-negative. This matrix is also irreducible because the elements just above and below the diagonal are all positive. It follows that all elements of Y (t) are positive for t > 0. The theorem of Perron and Frobenius can be applied. The eigenvalue μ1 of Y (T ) with largest real part is in fact real and positive and the corresponding eigenspace is one-dimensional. Moreover (1 1 . . . 1)Q(t) = (−b 0 0 . . . 0): Q(t) is “leaky” and 0 < μ1 < 1 (Aronsson and Kellogg 1978). Then λ1 = (log μ1 )/T < 0. So the vector (1, 0, 0, . . . , 0) is a steady state, towards which P(t) converges as t → +∞. The goal here is to estimate for N large how close λ1 is to 0. Let v be an eigenvector of X (T ) associated with the eigenvalue μ1 = eλ1 T . One can choose v such that vn > 0 for 1 ≤ n ≤ N . Hence X (T )v = eλ1 T v. As 

This matrix has the block structure M(t) =

123

Stochastic SIS epidemic model

in Floquet theory, set π(t) = e−λ1 t X (t)v. Then dπ dt (t) = −λ1 π(t) + M(t)π(t). Moreover π(T ) = e−λ1 T X (T )v = v = π(0). So π(t) is T -periodic. Write π(t) = (πn (t))0≤n≤N . Then λ1 πn +

dπn = a(t)(n − 1)(1 − (n − 1)/N )πn−1 − [a(t)n(1 − n/N ) + b n]πn dt + b(n + 1)πn+1 . (8)

N d N Summing these equations, one gets λ1 n=0 πn (t) + dt n=0 πn (t) = 0, so

N

N

N

N −λ t 1 =e T -periodic. Therefore n=0 πn n=0 πn (t)

n=0 πn (0). But n=0 πn (t) is

N N (0) = 0 and n=0 πn (t) = 0 for all t. So π0 (t) = − n=1 πn (t). But (8) with n = 0 0 also shows that λ1 π0 (t) + dπ dt = b π1 (t). Integrating over one period and using the periodicity of π0 (t) leads to T λ1 = b

0 T 0

T

π1 (t) dt

π1 (t) dt = −b N 0  T . π0 (t) dt n=1 0 πn (t) dt

(9)

WKB ansatz and the Hamilton–Jacobi equation. When N is large, let us try the WKB ansatz πn (t)  e−N S(t,x) for 1 ≤ n ≤ N , where x = n/N and S(t, x) is a continuous function of t and x for 0 < x < 1, which is T -periodic with respect to t. Then dπn ∂S  −N (t, x) e−N S(t,x) , dt ∂t 1

∂S

∂S

πn+1 (t)  e−N S(t,x+ N )  e−N S(t,x)− ∂ x (t,x) , πn−1 (t)  e−N S(t,x)+ ∂ x (t,x) . Set α(t, x) = a(t)x(1 − x) and β(x) = b x. Then (8) may be written as λ1 πn +

dπn = N α(t, x − 1/N )πn−1 − N [α(t, x)+β(x)]πn + Nβ(x +1/N )πn+1 . dt

Keeping only the highest order terms, one can use α(t, x − 1/N )  α(t, x) and β(x + 1/N )  β(x) to get λ1 πn +

dπn  N α(t, x)[πn−1 − πn ] + Nβ(x)[πn+1 − πn ]. dt

As λ1 is expected to be exponentially small, one can neglect it on the left side. Plugging the WKB ansatz and dividing by N e−N S(t,x) , one gets the Hamilton–Jacobi equation  ∂S   ∂S  ∂S + a(t)x(1 − x) e ∂ x − 1 + b x e− ∂ x − 1 = 0 ∂t

(10)

for 0 < x < 1. This is of the form (6) with the time-periodic Hamiltonian H (t, x, p) given by (2), a(t) replacing a. Equations such as (10) usually have asymptotic solutions of the form S(t, x) = −Et + (t, x) where (t, x) is T -periodic with respect to t

123

N. Bacaër

and E is a constant. Here however only the solutions with E = 0 are of interest: they correspond to the heteroclinic orbits of Sect. 2.2 below. 1 Boundary conditions. As H (t, 0, p) = 0, it follows that ∂∂tS (t, 0) = 0. So S(t, 0) is a constant S0 independent of t. As (10) involves only partial derivatives of S(t, x), its solutions are defined up to an additive constant (recall that the eigenvector v of X (T ) was defined up to a multiplicative constant). So one may choose S0 = 0, hence the Dirichlet condition S(t, 0) = 0 ∀t. (11) Moreover, considering that πn (t) = 0 for n > N and noticing on formula (3) for constant environments that S(1) is finite but dd xS (1) = +∞, one imposes ∂S (t, 1) = +∞ ∀t. ∂x

(12)

Another way of presenting this boundary condition is the “state constraint”   ∂S ∂S + H t, x, ≥0 ∂t ∂x when x = 1. This condition leads for x = 1 to ∂H ∂ p (t, x,

∂H ∂S ∂ p (t, x, ∂ x )

≥ 0 (Soner 1986). But

a(t)x(1 − x)e p − bxe− p , this expression is nonnegative at

p) = since if and only if p = +∞, as in (12).

x =1

Properties of the Hamiltonian. The Hamiltonian H (t, x, p) is convex in p as ∂2 H (t, x, p) = a(t) x(1 − x) e p + b x e− p ≥ 0. Moreover H (t, x, p) → +∞ as ∂ p2 | p| → +∞ provided 0 < x < 1; but this coercivity property is no longer true when x = 0 or x = 1. Notice that H (t, x, 0) = 0. The Lagrangian is L(t, x, v) = max p { pv − H (t, x, p)}. When 0 < x < 1, one has L(t, x, v) = p∗ v − H (t, x, p∗ ) with p∗ the unique solution of v = ∂∂Hp (t, x, p∗ ) = a(t)x(1 − x)e p∗ − bxe− p∗ . This is a second order polynomial equation in e p∗ . It leads to L(t, x, v) = p∗ v − a(t)x(1 − x)(e p∗ − 1) − bx(e− p∗ − 1)    v + v 2 + 4a(t)x(1 − x)bx = v log + a(t)x(1 − x) + bx 2a(t)x(1 − x)  2a(t)x(1 − x)bx v + v 2 + 4a(t)x(1 − x)bx  − . − 2 v + v 2 + 4a(t)x(1 − x)bx For x = 1, one gets L(t, 1, v) = +∞ if v > 0, L(t, 1, 0) = b, and L(t, 1, v) = −v log(−v/b) + v + b if v < 0. For x = 0, one gets L(t, 0, v) = +∞ if v = 0 and L(t, 0, 0) = 0. For x close to 0, notice however that L(t, x, v) ∼ −v log x. So for small η > 0 and any function ξ ∈ C 1 ([θ, t]; [0, 1]) such that ξ(θ ) = 0, one gets

123

Stochastic SIS epidemic model

 θ+η θ

L(s, ξ(s), ξ˙ (s)) ds  −

finite.

 θ+η

dξ ds (s)

θ

log ξ(s) ds = −

 ξ(θ+η) 0

log ξ dξ , which is

Solutions of the Hamilton–Jacobi equation. For a given initial datum S0 (x), the function  S(t, x) = inf

t θ

L(s, ξ(s), ξ˙ (s)) ds + 1θ=0 S0 (ξ(θ )) ; 

0 ≤ θ ≤ t, ξ ∈ C 1 ([θ, t]; [0, 1]), θ = 0 or ξ(θ ) = 0, ξ(t) = x is a viscosity solution of (10) with the mixed Dirichlet-“state constraint” boundary conditions (11)–(12) such that S(0, x) = S0 (x) (Barles 1994). This is the value function of the “exit time” problem at x = 0 with “state constraint” at x = 1. A time-periodic solution S ∗ (t, x) of (10)–(11)–(12) is thus given by a fixed point of the evolution operator above: S ∗ (0, x) = S ∗ (T, x). Roquejoffre (2001) and Mitake (2009) considered similar time-periodic Hamilton-Jacobi equations with Dirichlet boundary conditions. Notice however that there is no uniqueness (for related problems, see Barles and Perthame 1988). Indeed consider the special case where a(t) = a0 is constant. In this case there are two types of stationary viscosity solutions S ∗ (x): on one side the solutions of the form x log(b/a0 ) + x + (1 − x) log(1 − x) + γ with some constant γ ≤ 0, which differ through the constant γ , only that with γ = 0 satisfying the boundary condition at x = 0 in the classical sense; on the other side the solutions of the form min{0, x log(b/a0 ) + x + (1 − x) log(1 − x) + γ } for some constant γ such that 0 < γ ≤ c0 . The latter are identically zero near x = 0 and hence do not yield the right value for C. For the time-periodic equation (10) with the mixed boundary conditions (11)–(12), one can conjecture that it has viscosity solutions S ∗ (t, x) that are T -periodic with respect to t, that are not identically zero near x = 0, and that differ only by a constant (thus yielding the same C). It is such a solution that is chosen for the ansatz. As suggested in Fig. 3 below, the boundary condition at x = 0 has to be understood in a viscosity sense since S ∗ (t, x) may be discontinuous at x = 0. Large N behavior of λ1 . Returning to (9) one gets log b 1 log(−λ1 ) = + log N N N Notice that π1 (t)  e−N S 1 log N



∗ (t,1/N )



T 0

T

0

 N   T 1 π1 (t) dt − log πn (t) dt . N 0 

n=1

 e−N S 

π1 (t) dt

∗ (t,0+ )

for N large. So

−→ − min S ∗ (t, 0+ )

N →+∞

0≤t≤T

because of the method of Laplace for the asymptotic evaluation of integrals. Similarly, ∗ as πn (t)  e−N S (t,n/N ) , one gets

123

N. Bacaër

  N  T 1 πn (t) dt −→ − min min S ∗ (t, x) log N →+∞ 0≤t≤T 0≤x≤1 N 0 n=1

and log(−λ1 ) −→ −C N →+∞ N with C as in (5). One can conjecture that C > 0 if and only if

1 T

T 0

a(t) dt > b.

Mean extinction time. The mean extinction time τn (t) starting from n infected people at time t is a T -periodic solution of the system −1=

dτn + b nτn−1 − (a(t)n(1 − n/N ) + b n)τn + a(t)n(1 − n/N )τn+1 (13) dt

for 1 ≤ n ≤ N , with τ0 (t) = 0. This system involves the adjoint Q ∗ (t) of the matrix Q(t). Set τˆ (t) = (τn (t))1≤n≤N , πˆ (t) = (πn (t))1≤n≤N , and 1 = (1, 1, . . . , 1). Then λ1 πˆ +

d πˆ d τˆ = Q(t)π, ˆ −1 = + Q ∗ (t)τˆ . dt dt

Let ·, · be the usual scalar product of real vectors. Then d πˆ d τˆ d πˆ , τˆ  =  , τˆ  + πˆ ,  = Q(t)π, ˆ τˆ  − λ1 πˆ , τˆ  − πˆ , 1 − πˆ , Q ∗ (t)τˆ . dt dt dt The terms involving Q(t) and Q ∗ (t) cancel. Integrating over one period and using the fact that both πˆ (t) and τˆ (t) are T -periodic, one gets T

πˆ , 1 . −λ1 =  T0 ˆ , τˆ  dt 0 π This suggests that the mean time to extinction τ , starting for example from one infected person at time 0, is of the same order of magnitude as −1/λ1 : log(τ ) −→ C. N N →+∞ One can conjecture that this mainly informal derivation can be proved rigorously, as for the SIS model in a constant environment (Nåsell 2011). 2.2 The heteroclinic orbit General case. Recall that the Hamilton–Jacobi equation (6) may be solved at least locally by ray tracing, i.e., by solving the Hamiltonian system (4) together with the equation

123

Stochastic SIS epidemic model

∂H dz = p(t) (t, x(t), p(t)) − H (t, x(t), p(t)) dt ∂p and the initial conditions x(0) = x0 , p(0) = ∂∂ xS (0, x0 ), z(0) = S(0, x0 ), so that z(t) = S(t, x(t)). Following Assaf et al. (2008), consider more closely the Hamiltonian system (4). In the present case, ∂H (t, x, p) = a(t)x(1 − x)e p − bxe− p , ∂p ∂H (t, x, p) = a(t)(1 − 2x)(e p − 1) + b(e− p − 1). ∂x

(14)

Let us first look for a nontrivial T -periodic solution such that x ≡ 0 and dp ∂H =− (t, 0, p) = −(a(t) − b e− p )(e p − 1). dt ∂x Setting p = log(1 + q), one gets a Bernoulli differential equation that can easily be solved. This leads to the T -periodic solution ⎛



p ∗ (t) = log ⎝1 +

t e−bt+ 0 a(s) ds

ep

∗ (0)

−1

+

t

a(s) e

t −b(t−s)+ s a(u) du

−1 ⎞ ds

⎠,

0

where ⎛



T

p ∗ (0) = log ⎝1 +  T 0

1 − e−bT +

a(s) ds

0

T

a(s) e−b(T −s)+

s

a(u) du

⎠. ds

˜ and p(t) = The periodic solution (0, p ∗ (t)) is unstable. Indeed, setting x(t) = x(t) ˜ and linearizing the equations, one gets p ∗ (t) + p(t) 

d x/dt ˜ d p/dt ˜



 =





0 a(t)e p (t) − be− p (t) ∗ ∗ ∗ 2a(t)(e p (t) − 1) −a(t)e p (t) + be− p (t)

  x˜ . p˜

T ∗ ∗ The Floquet multipliers are f = exp 0 [a(t)e p (t) − be− p (t) ]dt and 1/ f , hence the instability. The instability may also be seen as a consequence of Liouville’s theorem concerning the invariance of the “volume” in the phase space (x, p) under a Hamiltonian flow. Secondly let us look for a nontrivial T -periodic solution such that p ≡ 0 and ∂H dx = (t, x, 0) = a(t)x(1 − x) − bx. dt ∂p

123

N. Bacaër

This is the mean field equation for the SIS model. The unique nonzero T -periodic solution is given by x ∗ (t) =



t 1 bt− 0 a(s) ds e + x ∗ (0)

with



t

t

a(u) eb(t−u)−

x (0) =  T 0

1 − ebT − a(u) e

a(s) ds

−1 du

.

0 T



u

0

a(s) ds

T b(T −u)− u a(s) ds

.

(15)

du

(x ∗ (t), 0)

The periodic solution is also unstable. Indeed, setting x(t) = x ∗ (t) + x(t) ˜ and p(t) = p(t) ˜ and linearizing the equations, one gets 

d x/dt ˜ d p/dt ˜



 =

a(t)(1 − 2x ∗ (t)) − b a(t)x ∗ (t)(1 − x ∗ (t)) + bx ∗ (t) 0 −a(t)(1 − 2x ∗ (t)) + b

  x˜ . p˜

The Floquet multipliers are once again inverse of each other, hence the instability. Recall from Sect. 1 that in a constant environment there is a heteroclinic orbit in the plane (x, p) from the stationary point (x ∗ , 0) = (1 − b/a, 0) to the stationary point (0, p ∗ ) = (0, log(b/a)) when a > b. Recall also that Escudero and Rodríguez (2008) found a similar heteroclinic orbit for the time-periodic model involving branching and annihilation of identical particles, at least for a small amplitude periodic perturbation. So one can also expect the existence of a heteroclinic orbit (x(t), ˆ p(t)) ˆ from the periodic solution (x ∗ (t), 0) to the periodic solution (0, p ∗ (t)); the existence can probably be proved using a variational approach, as in (Rabinowitz 1994). This special orbit may be found numerically by a shooting method. According to Assaf et al. (2008, Eq. (20)),  +∞  ∂H p(t) ˆ (t, x(t), ˆ p(t)) ˆ − H (t, x(t), ˆ p(t)) ˆ dt. (16) C= ∂p −∞ This integral can be estimated numerically. Perturbation method. When a(t) is a constant a0 , call (xˆ0 (t), pˆ 0 (t)) the heteroclinic orbit connecting the stationary points (x ∗ , 0) = (1 − b/a0 , 0) and (0, p ∗ ) = (0, log(b/a0 )). This orbit is such that a0 (1 − x)e p − b = 0, as can be seen from (2). Using this equation to express p as a function of x and inserting the result in the first equation in (4), one gets ddtx = b x − a0 x(1 − x). The solution is  x(t) =

1 (a0 −b)(t−t0 ) a0 e (1 − e(a0 −b)(t−t0 ) ) + x(t0 ) a0 − b

−1

.

Choosing for example x(t0 ) = (1 − b/a0 )/2, one gets xˆ0 (t) =

123

1 − b/a0 1 + e(a0 −b)(t−t0 ) . and p ˆ (t) = log 0 1 + e(a0 −b)(t−t0 ) 1 + e(a0 −b)(t−t0 ) a0 /b

Stochastic SIS epidemic model

Assume now that a(t) = a0 (1+ε φ(t)) with a0 > b, ε small, and φ(t) a T -periodic T function such that 0 φ(t) dt = 0. The Hamiltonian can be written as H (t, x, p) = H0 (x, p) + ε H1 (t, x, p), where H0 (x, p) is the same as (2) with a replaced by a0 , and where H1 (t, x, p) = a0 φ(t) x(1 − x)(e p − 1). Set c0 = b/a0 − 1 − log(b/a0 ). Assaf et al. (2008, Eq. (24)) showed that C  min Γ (t0 ), Γ (t0 ) = c0 − ε t0

+∞

−∞

H1 (t, xˆ0 (t), pˆ 0 (t)) dt

for ε close to 0. In the present case, (1 − xˆ0 )e pˆ0 = b/a0 . So Γ (t0 ) = c0 − ε a0

+∞ −∞

  φ(t) xˆ0 (t) b/a0 − 1 + xˆ0 (t) dt

= c0 + ε(1 − b/a0 )



+∞ −∞

φ(t0 + u/(a0 − b))

eu du. (1 + eu )2

T Thus Γ (t0 ) is a T -periodic function of t0 such that 0 Γ (t0 ) dt0 = 0. Consider the

k iωt , where ω = 2π/T , φ = 0 Fourier decomposition of φ(t), φ(t) = +∞ 0 k=−∞ φk e ¯ since the average of φ(t) is zero, and φ−k = φk (the complex conjugate). Then

Γ (t0 ) = c0 + ε(1 − b/a0 ) = c0 + ε(1 − b/a0 )

+∞  k=−∞ +∞  k=−∞

φk ek iωt0 φk ek iωt0

+∞ −∞

k iωu

e a0 −b

kπ ω a0 −b



sinh

kπ ω a0 −b

eu du (1 + eu )2 

(see Appendix). In particular if φ(t) = cos(ωt), then φ±1 = 1/2 and φk = 0 otherwise. So π ω cos(ωt0 ) .  Γ (t0 ) = c0 + ε (17) ω a0 sinh aπ0 −b Recall from Escudero and Rodríguez (2008) and Assaf et al. (2008) that the perturbed system is of the form ∂ H0 ∂ H1 dx = +ε , dt ∂p ∂p

dp ∂ H0 ∂ H1 =− −ε , dt ∂x ∂x

(18)

123

N. Bacaër

and that xˆ0 (t) and pˆ 0 (t) depend only on t − t0 ; so the Melnikov function M(t0 ) is +∞ 

 ∂ H1 ∂ H0 ∂ H1 ∂ H0 M(t0 ) = − + (t, xˆ0 (t), pˆ 0 (t)) dt ∂x ∂p ∂p ∂x −∞  +∞  ∂ H1 d xˆ0 ∂ H1 d pˆ 0 − − (t, xˆ0 (t), pˆ 0 (t)) dt = ∂ x dt ∂ p dt −∞  +∞  ∂ H1 d xˆ0 ∂ H1 d pˆ 0 1 dΓ = + . (t, xˆ0 (t), pˆ 0 (t)) dt = − ∂ x dt0 ∂ p dt0 ε dt0 −∞

Using (17), one gets

M(t0 ) =

π ω sin(ωt0 ) .  ω a0 sinh aπ0 −b

So M(t0 ) crosses 0 for t0 = kπ/ω (k integer). It follows that the heteroclinic orbit exists at least for ε small. The minimum of Γ (t0 ) in (17) is obtained for t0 = T /2 if ε > 0 and for t0 = 0 if ε < 0: in both cases, one gets C  c0 −

π ω|ε|   ω a0 sinh aπ0 −b

(19)

for ε close to 0, as announced in the introduction. Notice that (19) resembles equation (4.76) of Kamenev (2011) derived starting from a slightly different Hamiltonian. When ω is small (or T is large) so that ω  a0 , then (19) shows that C  c0 − |ε|(1 − b/a0 ) ,

(20)

which is independent of ω. This formula is the same as that obtained from (1) with a = a0 (1 − |ε|): b b b b − 1 − log = − 1 − log − |ε|(1 − b/a0 ) + o(ε) a0 (1 − |ε|) a0 (1 − |ε|) a0 a0 when ε is close to 0. As in the “adiabatic approximation” of Assaf et al. (2008, Sect. 4), formula (20) is expected to hold not just for small ε but as long as ω  a0 and the right side of (20) is positive. As sinh(x) ≥ x for all x ≥ 0, one may notice that the approximate value of C given by (20) is always less than that given by (19).

123

Stochastic SIS epidemic model

High-frequency limit. Let us assume now that ω a0 , still with φ(t) = cos(ωt). System (18) may be written as ∂ H0 dx = (x, p) + a0 ε cos(ωt)x(1 − x)e p dt ∂p ∂ H0 dp =− (x, p) − a0 ε cos(ωt)(1 − 2x)(e p − 1). dt ∂x Following the Kapitsa method (Assaf et al. 2008, §III.B), write x(t) = X (t) + ξ(t) and p(t) = Y (t) + η(t), where X and Y are slow variables, while ξ and η are small but fast oscillations. The fast oscillating terms balance each other: dξ  a0 ε cos(ωt)X (1 − X )eY , dt

dη  −a0 ε cos(ωt)(1 − 2X )(eY − 1). dt

Considering X and Y as constant during the short period T = 2π/ω, one gets ξ(t) 

a0 ε a0 ε sin(ωt)X (1 − X )eY , η(t)  − sin(ωt)(1 − 2X )(eY − 1). ω ω

This suggests the transformation a0 ε sin(ωt)X (1 − X )eY ω a 2 ε2 a0 ε sin(ωt)(1 − 2X )(eY − 1) + 0 2 Φ(t, X, Y ) , p=Y− ω ω

x=X+

where Φ(t, X, Y ) is chosen for the transformation to be almost canonical, i.e., the Poisson brackets must satisfy the condition {x, p} =

∂x ∂p ∂x ∂p − = 1 + o(a02 /ω2 ) . ∂ X ∂Y ∂Y ∂ X

(21)

Since     2 ε2 a a0 ε ε ∂Φ a 0 {x, p} = 1 + sin(ωt)(1 − 2X )eY 1 − sin(ωt)(1 − 2X )eY + 0 2 ω ω ω ∂Y    a ε a ε a 2 ε2 ∂Φ 0 0 sin(ωt)X (1 − X )eY 2 sin(ωt)(eY − 1) + 0 2 − , ω ω ω ∂X condition (21) may be written as a02 ε2 ∂Φ a02 ε2 2 2 2Y sin (ωt)(1 − 2X ) e + ω2 ω2 ∂Y 2 2 a ε − 2 0 2 sin2 (ωt)X (1 − X )eY (eY − 1) + o(a02 /ω2 ) = 1 + o(a02 /ω2 ). ω

{x, p} = 1 −

123

N. Bacaër

So   ∂Φ = sin2 (ωt) (1 − 2X )2 e2Y + 2X (1 − X )eY (eY − 1) . ∂Y For Φ(t, X, 0) = 0 to hold, one must chose   Φ(t, X, Y ) = sin2 (ωt) (1 − 2X )2 (e2Y − 1)/2 + X (1 − X )(eY − 1)2 . The generating function F2 (t, x, Y ) of this transformation, such that o(a02 /ω2 ) and ∂∂Fx2 = p + o(a02 /ω2 ), is given by

∂ F2 ∂Y

= X +

a0 ε sin(ωt)x(1 − x)(eY − 1) ω a 2 ε2 + 0 2 sin2 (ωt)x(1 − x)(1 − 2x)(e2Y − 1). 2ω

F2 (t, x, Y ) = xY −

Set H (t, x, y) = h(t, X, Y ). The new Hamiltonian is h(t, X, Y ) + ∂∂tF2 . Taking the average of this Hamiltonian over one period T = 2π/ω, the second term cancels out T T as 0 ∂∂tF2 dt = 0 and only the effective Hamiltonian H¯ (X, Y ) = T1 0 h(t, X, Y ) dt T remains. A tedious computation, using that T1 0 sin2 (ωt) dt = 1/2, leads to  H¯ (X, Y ) X (1 − e−Y ) a0 (1 − X )eY − b +

a02 ε2  −a0 X (1 − X )2 e2Y 2ω2



+ b(1 − X )(1 − 2X )e − bX (1 − X )(e − 1) − b(1 − 2X ) Y

Y

2

.

One gets the perturbed heteroclinic orbit by setting the term in brackets to zero. It links (X ε∗ , 0) to (0, Yε∗ ), where X ε∗

 b(a0 − b) ε2 a0 (a0 − b)ε2 ∗ , Y  (1 − b/a0 ) 1 −  log(b/a ) + . 0 ε 2ω2 2ω2 

0 The action along this heteroclinic orbit is C = X ∗ Y d X . Another tedious computation ε finally leads to (a0 − b)2 ε2 C  c0 − (1 + 2b/a0 ) , (22) 12 ω2 as announced in the introduction. Since the function z → (1 − z)2 (1 + 2z) is less than a 2 ε2

1 on the interval 0 < z < 1, the correction term for C is always less than 120 ω2 , which is small since ω a0 by assumption. Not unexpectedly, a population subject to a high frequency perturbation does not depend much on the amplitude ε of this perturbation.

123

Stochastic SIS epidemic model

3 Numerical computations Floquet multipliers. λ1 may be estimated directly by computing the Floquet multipliers of the master equation using a software (such as Scilab) that can solve ordinary differential equations and compute matrix eigenvalues numerically. eλ1 T is the eigenvalue with a real part that is second largest, the largest being 1. One can then plot − log(−λ1 ) as a function of N . The slope of this curve gives an estimate of C. Heteroclinic orbit. As in Assaf et al. (2008), the orbit connecting (x ∗ (t), 0) to (0, p ∗ (t)) can be obtained by a shooting method, taking as initial condition x ∗ (0) given by (15) and a very small negative value for p(0). This value is varied until getting a solution (x(t), p(t)) that tends to become periodic, i.e., with x(t) approaching 0 and p(kT ) approaching p ∗ (0) for k large (but not too large to avoid numerical instability). One can then use the integral (16) to compute C numerically. PDE method. It is also possible to compute a periodic solution S ∗ (t, x) of the Hamilton–Jacobi partial differential equation (PDE) (6) using the numerical methods of the theory of viscosity solutions. For example, let Δt be the time step and Δx be the spatial step. Let S mj be an approximation of S(mΔt, jΔx), where j and m are integers such that m ≥ 0 and 0 ≤ j ≤ J with J = 1/Δx. One can use the following Godunov type scheme S m+1 − S mj j Δt



S mj − S mj−1 S mj+1 − S mj + H mΔt, jΔx, , Δx Δx

 = 0,

where the numerical Hamiltonian H(t, x, p − , p + ) is given by H(t, x, p − , p + ) =



min{H (t, x, p); p − ≤ p ≤ p + } if p − < p + , max{H (t, x, p); p + ≤ p ≤ p − } if p + ≤ p −

(Osher and Shu 1991). As H (t, x, p) is convex in p, the second expression involving a maximum is equal to max{H (t, x, p + ), H (t, x, p − )}. As for the first expression involving a minimum, notice from (14) that H (t, x, p) has a minimum with respect b . So to p when ∂∂Hp = 0, i.e., when p = p  = 21 log a(t)(1−x) ⎧ ⎨ H (t, x, p + ) if p − < p + ≤ p  , − + min{H (t, x, p); p ≤ p ≤ p } = H (t, x, p − ) if p  ≤ p − < p + , ⎩ H (t, x, p  ) if p − ≤ p  ≤ p + . m For the boundary conditions, one sets S0m = 0 and (S m J − S J −1 )/Δx = K with some large value for K . The time step Δt must be small enough compared to Δx for the Courant–Friedrichs–Lewy (CFL) condition to hold. As initial condition we took S(0, x) = x log(b/a0 ) + x + (1 − x) log(1 − x), the smooth stationary solution when a(t) is replaced by its time average. One could also take a constant function S(0, x) = σ with σ sufficiently negative, but the convergence to the periodic regime is then slower.

123

N. Bacaër Fig. 1 Computing the Floquet multipliers of the master equation: − log(−λ1 ) as a function of N for ε = 0.2, 0.5 or 0.8 and N = 10, 20, …, 60. C is the slope of these lines. Parameter values: T = 1, a0 = 20, b = 5

60 50 40 30 20 10 0 0

10 20 30 40 50 60 70 80 90 100

The constant σ has to be chosen sufficiently negative to avoid the non-uniqueness problem already mentioned in Sect. 2. Once the solution of the time-dependent problem has reached a periodic regime, one can estimate C = mint S ∗ (t, 0+ )−mint,x S ∗ (t, x). Monte Carlo method. The mean extinction time can also be estimated by a Monte Carlo method, averaging stochastic simulations. Notice however that the Gillespie algorithm taking advantage of exponential waiting times cannot be used because a(t) is periodic in time. As N increases, time to extinction can become very large. We do not show any results using this method. Examples. Assume that a(t) = a0 (1+ε cos(2π t/T )) with T = 1 week. First consider the case where a0 = 20 per week and b = 5 per week. The mean infectious period is 1/b = 1.4 days. Thus R0 = a0 /b = 4 > 1 and c0 = b/a0 − 1 − log(b/a0 )  0.636. Figure 1 shows − log(−λ1 ) as a function of N for ε = 0.2, 0.5, or 0.8 and N = 10, 20, …, 60, computed using the Floquet multipliers. The lines correspond to a linear regression of the last 3 points N = 40, 50, 60. The slopes of these lines—estimates of C—were 0.524, 0.364, and 0.225 for ε = 0.2, 0.5, and 0.8 respectively. Notice in this example that a0 and ω = 2π/T are of the same order, a case one can call intermediate frequency. So one expects (19) to give a good approximation for C when ε is small. Figure 2 shows the following curves as a function of ε for 0 ≤ ε ≤ 1: a plain line for C using the heteroclinic orbit and a dashed line with long dashes for C using the PDE method with Δx = 0.002 and Δt = 0.0002 (these first two lines are almost indistinguishable); three dots representing the values of C obtained in Figure 1 (notice how they fall on the two previous lines); a dashed line with short dashes for the approximate formula (19); a dash-dotted line for the low-frequency approximation (20). One can see that the approximate formula (19) is close to C even when ε is only moderately small. Figure 3 shows a time-periodic solution S ∗ (t, x) of the Hamilton–Jacobi PDE, plotted as a function of x for different values of t, when ε = 0.5. Notice the discontinuity of the solution at x = 0. A zoom near x = 0 would show that S ∗ (t, 0+ ) is indeed time-periodic, so that the boundary condition S ∗ (t, 0) = 0 can be satisfied only in a generalized sense.

123

Stochastic SIS epidemic model Fig. 2 Intermediate frequency: C computed using the heteroclinic orbit (plain line) or the PDE method (dashed line with long dashes) (the two lines are almost indistinguishable), the Floquet multipliers as in Fig. 1 (dots), the approximate formula (19) (dashed line with short dashes), and the low-frequency formula (20) (dash-dotted line), as a function of ε. Same parameter values as in Fig. 1

0.7

C

0.6 0.5 0.4 0.3 0.2 0.1 0 0

Fig. 3 A time-periodic solution S ∗ (t, x) of the Hamilton–Jacobi PDE, plotted as a function of x for t = 0 (plain line), t = T /4 (long dashes), t = T /2 (short dashes), and t = 3T /4 (dash-dotted line). Same parameter values as in Fig. 1 and ε = 0.5

0.2

0.4

0.6

0.8

1

0.6

0.8

1

0 −0.1 −0.2 −0.3 −0.4 −0.5 −0.6 −0.7 0

0.1935

0.2

0.4

C

0.193 0.1925 0.192 0.1915 0.191 0.1905 0.19 0.1895 0.189 0.1885 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Fig. 4 High-frequency regime: C computed using the heteroclinic orbit (plain line) and the high-frequency formula (22) (dashed line) as a function of ε. Parameter values: T = 1, a0 = 2, b = 1

Figure 4 considers a high-frequency example: a0 = 2 per week and b = 1 per week. Thus R0 = 2 and c0  0.1931. In this case ω  6.28 per week is somewhat bigger than a0 . The number C is computed using the heteroclinic orbit and the high-frequency

123

N. Bacaër Fig. 5 The components t  → x(t) ˆ and t  → p(t) ˆ of the heteroclinic orbit (x(t), ˆ p(t)) ˆ connecting the two periodic solutions (0, p ∗ (t)) and (x ∗ (t), 0). Same parameter values as in Fig. 4 and ε = 0.1

0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 0

2

4

6

8

10 12 14 16 18 20

formula (22) as a function of ε for 0 ≤ ε ≤ 1. There is a good agreement over the entire ε range. Finally Fig. 5 shows the orbit connecting (x ∗ (t), 0) to (0, p ∗ (t)) for the same parameter values when ε = 0.1. 4 Remarks – It may be possible to get more precise estimates with the refined WKB ansatz πn (t)  e−N S0 (t,n/N )−S1 (t,n/N ) . Inserting πn+1 (t)  e

−N S0 (t,n/N )−

2 ∂ S0 1 ∂ S0 1 ∂ S1 ∂ x (t,n/N )− 2N ∂ x 2 (t,n/N )−S1 (t,n/N )− N ∂ x (t,n/N )

and a similar expression for πn−1 (t) in (8), and separating the highest-order terms, one gets the Hamilton–Jacobi equation (10) for S0 (t, x) and the transport equation ∂ S0 ∂ S0  ∂ S1 ∂ S1  + a(t)x(1 − x)e ∂ x − bxe− ∂ x ∂t ∂x     ∂ S0 ∂S x(1 − x) ∂ 2 S0 x ∂ 2 S0 − ∂ x0 ∂ x 1 − 2x + +be −1 + = a(t) e 2 ∂x2 2 ∂x2

for S1 (t, x). If a(t) is periodic then both S0 (t, x) and S1 (t, x) have to be computed numerically. Recall from Assaf and Meerson (2010) that if a(t) = a is constant then the pre-exponential factor for λ1 may be obtained by the method of matched asymptotic expansions as follows. First S0 (x) is given by (3) and the equation for √ −N S0 (n/N ) S1 (x) leads to S1 (x) = log(x 1 − x) + const. So πn  κ en √ n for some constant κ. For small n, πn 

N

κ N −N S0 (0)−nS  (0) κ N −N S0 (0)  a n 0 e e = . n n b

1− N

(23)

But for small n, system (8) may also be approximated by the recurrence a(n − n 1)πn−1 − n(a + b)πn + b(n + 1)πn+1  0 for n ≥ 1, giving πn  πn1 1−(a/b) 1−a/b .

123

Stochastic SIS epidemic model (a/b) As n increases, this expression is equivalent to πn1 a/b−1 . This coincides with (23) a −N S (0) 0 ( b − 1). Finally (9) gives only if π1  κ N e n

$ √ % ∗ 1 − x ∗ N S  (x ∗ ) (b − a)x 0 −b π1 /(κ N ) (a − b)2 N λ1   1 −N S0 (x)  =− √ ∗ e√ a ecN 2π e N [S0 (0)−S0 (x )] 2π 0 x 1−x d x (Assaf and Meerson 2010, Eq. (71)). If for example a = 20, b = 5, and N = 50, then this estimate is just 2 % above the value of λ1 obtained with some software computing eigenvalues of large matrices as in Fig. 1. The method of matched asymptotics can probably be extended to the periodic case but will unlikely lead to such an explicit asymptotic formula.

N Pn (t) z n with 0 ≤ z ≤ 1. Then – Consider the generating function φ(t, z) = n=0 φ(t, 1) = 1 for all t. A simple computation starting from (7) shows that   a(t)z ∂φ a(t) 2 ∂φ ∂ 2φ = (1 − z) b + − a(t)z + z (1 − z) 2 ∂t N ∂z N ∂z for 0 < z < 1. In the quasi-stationary regime, one expects that φ(t, z)  1 + eλ1 t ψ(t, z) with ψ(t, z) periodic with respect to t, ψ(t, 1) = 0, and λ1 ψ +

  ∂ψ ∂ 2ψ a(t)z a(t) 2 ∂ψ = (1 − z) b + − a(t)z + z (1 − z) 2 . ∂t N ∂z N ∂z

Thus λ1 is also the largest nonzero eigenvalue of this parabolic problem. This might be a way to prove more rigorously the asymptotic results concerning λ1 for N large. – Setting Pn (t) = P(t, x) where x = n/N and performing a second-order Taylor expansion of the master equation (7), one gets the Fokker-Planck or diffusion equation     ∂  1 ∂ 2  ∂P =− a(t)x(1 − x) − bx P + a(t)x(1 − x) + bx P . ∂t ∂x 2N ∂ x 2 Similarly, setting for the mean time to extinction τn (t) = τ (t, x) where x = n/N , (13) leads to the adjoint problem −1 =

 ∂τ  ∂ 2τ 1  ∂τ  + a(t)x(1 − x) − bx + a(t)x(1 − x) + bx . ∂t ∂x 2N ∂x2

However it is already known from the case of time-independent coefficients that these equations usually do not give the correct value for C (Doering et al. 2005); the value tends to be correct only when R0 is close to 1. The same problem occurs for other similar epidemic models. For example Diekmann et al. (2012, pp. 112113) used a diffusion approximation to estimate C for an SIR model despite having

123

N. Bacaër

R0  10 in their numerical application. For a WKB analysis of the stochastic SIR model, see Kamenev and Meerson (2008). – The number 1/C may be called the critical community size (Diekmann et al. 2012). Acknowledgments This article was stimulated by a discussion with Hans Metz concerning effective population sizes in periodic environments at the Conference on Stochastic Models in Ecology, Evolution and Genetics (Angers, December 2013).

Appendix Let us prove that

+∞

−∞

eiλu

eu πλ . du = u 2 (1 + e ) sinh(π λ)

First of all, eu /(1 + eu )2

= 1/(4 cosh2 (u/2)) is an even function. This fact combined with an integration by parts shows that

+∞

−∞

e

iλu

eu du = 2 (1 + eu )2



+∞

eu cos(λu) du (1 + eu )2 0 ∞  ∞ − cos(λu) λ sin(λu) =2 −2 du 1 + eu 0 1 + eu 0 ∞ −u e sin(λu) = 1 − 2λ du. 1 + e−u 0

Developing 1/(1 + e−u ) in a power series, one gets

+∞

−∞

eiλu



 eu du = 1 − 2λ (−1)n u 2 (1 + e ) = 1 + 2λ

n=0 ∞  2 n=0





e−(n+1)u sin(λu) du

0

(−1)n+1 . λ2 + (n + 1)2

The sum of this series may be computed by taking z = iπ λ in Euler’s formula ∞

1  1 2z = + (−1)n 2 , sin z z z − n2π 2 n=1

which holds for all complex number z such that z = nπ (n integer). As sin(iπ λ) = i sinh(π λ), one gets ∞

 (−1)n πλ = 1 + 2λ2 sinh(π λ) λ2 + n 2 n=1

and the result follows.

123

Stochastic SIS epidemic model

References Andersson H, Djehiche B (1998) A threshold limit theorem for the stochastic logistic epidemic. J Appl Probab 35:662–670 Aronsson G, Kellogg RB (1978) On a differential equation arising from compartmental analysis. Math Biosci 38:113–122 Assaf M, Kamenev A, Meerson B (2008) Population extinction in a time-modulated environment. Phys Rev E 78:041123 Assaf M, Meerson B (2010) Extinction of metastable stochastic populations. Phys Rev E 81:021116 Bacaër N, Ait Dads E (2012) On the biological interpretation of a definition for the parameter R0 in periodic population models. J Math Biol 65:601–621 Barles G (1994) Solutions de viscosité des équations de Hamilton–Jacobi. Springer, Berlin Barles G, Perthame B (1988) Exit time problems in optimal control and vanishing viscosity method. SIAM J Control Optim 26:1133–1148 Billings L, Mier-y-Teran-Romero L, Lindley B, Schwartz IB (2013) Intervention-based stochastic disease eradication. PLoS ONE 8(8):e70211 Diekmann O, Heesterbeek H, Britton T (2012) Mathematical tools for understanding infectious disease dynamics. Princeton University Press, Princeton Doering CR, Sargsyan KV, Sander LM (2005) Extinction times for birth-death processes: exact results, continuum asymptotics, and the failure of the Fokker-Planck approximation. Multiscale Model Simul 3:283–299 Escudero C, Rodríguez JA (2008) Persistence of instanton connections in chemical reactions with timedependent rates. Phys Rev E 77:011130 Hethcote H (1973) Asymptotic behavior in a deterministic epidemic model. Bull Math Biol 35:607–614 Kamenev A (2011) Field theory of non-equilibrium systems. Cambridge University Press, Cambridge Kamenev A, Meerson B (2008) Extinction of an infectious disease: a large fluctuation in a nonequilibrium system. Phys Rev E 77:061107 Mitake H (2009) Large time behavior of solutions of Hamilton–Jacobi equations with periodic boundary data. Nonlinear Anal 71:5392–5405 Nåsell I (2011) Extinction and quasi-stationarity in the stochastic logistic SIS model. Springer, Berlin Osher S, Shu CW (1991) High-order essentially nonoscillatory schemes for HamiltonJacobi equations. SIAM J Numer Anal 28:907–922 Ovaskainen O, Meerson B (2010) Stochastic models of population extinction. Trends Ecol Evol 25:643–652 Rabinowitz PH (1994) Heteroclinics for a reversible Hamiltonian system. Ergod Theor Dyn Syst 14:817– 829 Roquejoffre JM (2001) Convergence to steady states or periodic solutions in a class of Hamilton–Jacobi equations. J Math Pures Appl 80:85–104 Soner HM (1986) Optimal control with state-space constraint I. SIAM J Control Optim 24:552–561

123

On the stochastic SIS epidemic model in a periodic environment.

In the stochastic SIS epidemic model with a contact rate a, a recovery rate b < a, and a population size N, the mean extinction time τ is such that (l...
446KB Sizes 0 Downloads 6 Views