Persistence and distribution of a stochastic susceptible–infected–removed epidemic model with varying population size

Accepted Manuscript Persistence and distribution of a stochastic susceptible-infected-removed epidemic model with varying population size Lihong Chen, Fengying Wei PII: DOI: Reference:

S0378-4371(17)30432-6 http://dx.doi.org/10.1016/j.physa.2017.04.114 PHYSA 18216

To appear in:

Physica A

Received date: 9 December 2016 Please cite this article as: L. Chen, F. Wei, Persistence and distribution of a stochastic susceptible-infected-removed epidemic model with varying population size, Physica A (2017), http://dx.doi.org/10.1016/j.physa.2017.04.114 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Persistence and distribution of a stochastic susceptibleinfected-removed epidemic model with varying population size1 Lihong Chen, Fengying Wei

College of Mathematics and Computer Science, Fuzhou University, Fuzhou 350116, PR China

Abstract In this paper, the dynamics of a stochastic susceptible-infected-removed model in a population with varying size is investigated. We firstly show that the stochastic epidemic model has a unique global positive solution with any positive initial value. Then we verify that random perturbations lead to extinction when some conditions are being valid. Moreover, we prove that the solution of the stochastic epidemic model is persistent in the mean by building up a suitable Lyapunov function and using generalized Itˆ o’s formula. Further, the stochastic epidemic model admits a stationary distribution around the endemic equilibrium when parameters satisfy some sufficient conditions. To end this contribution and to check the validity of the main results, numerical simulations are separately carried out to illustrate these results. Keywords: Varying population size; Stochastic SIR model; Extinction; Persistence in the mean; Stationary distribution

1

Introduction

The susceptible-infected-recovered epidemic model is one of the most significant models when epidemiological patterns are taken into account. Kermack and McKendrick (1927) firstly proposed and investigated a classical SIR model. From then on, many papers studied various improved epidemic models and developed some good results and techniques, some of which assume that the populations have constant birth and mortality rates, for example, Zaman et al. (2008). Recently, some contributions about stochastic epidemic models introduced white noises into the corresponding deterministic models for a more real modelling and better understanding. Gray et al. (2011) studied an extension of a classical deterministic SIS epidemic model into a stochastic framework. Zhao et al. (2015) aimed at a stochastic susceptibleinfected-recovered-susceptible model in a population with varying size, in which they considered that the recovered also returned to susceptible individuals and introduced randomness into the contact rate. They proved the threshold of extinction, persistence and discussed the asymptotic behavior around the epidemic proportion equilibrium. Similarly, Zhao (2016) included stochastic perturbations into SIR epidemic model with saturated incidence and investigated the threshold of a stochastic SIR epidemic model and its extensions. Later, Liu and Chen (2016) investigated the dynamics of a stochastic SIR epidemic model, in which they supposed that the removed individuals did not become the susceptible and perturbations of parameters are set in terms of white noises. Zaman et al. (2008) discussed the ordinary differential equations of the deterministic SIR mode in a

1

population with varying size. The system can be expressed as follows: ˙ S(t)

=

˙ I(t)

=

˙ R(t)

=

bN − µS −

βSI , N

βSI − (µ + γ + ε)I, N γI − µR,

(1)

where S(t) is the density of the susceptible individuals, I(t) represents the member who is the infectious individuals and R(t) denotes the number of the individuals who have been removed from the compartment ˙ I˙ and R˙ are the of infection, the total population size is formulated by N (t) = S(t) + I(t) + R(t), and S, corresponding rates of change with respect to time t; b is the birth rate; µ is disease-free death rate; β is effective contact rate of infective individuals and βSI N is a standard incidence rate; γ is remove rate of infected individuals; ε is the rate related to death caused by disease. All parameters in model (1) are assumed to be nonnegative and b > 0, γ > 0. The differential equation of the total population size is governed by N˙ = (b − µ)N − εI. In 1990, Busenberg and Driessche (1990) considered the proportions of individuals in the three epidemiological classes, namely x=

S I R ,y = ,z = . N N N

(2)

That is to say, model (1) becomes the following model when x, y, z are substituted into: x˙ = b − µx − βxy, (3)

y˙ = βxy − (µ + γ + ε)y, z˙ = γy − µz.

The feasibility region of model (3) is Γ = {(x, y, z) ∈ R3+ |x + y + z = 1}. Obviously, models (1) and (3) are equivalent to each other under the proportions (2). From the recent discussions, in view of results obtained by Zhao et al. (2008), we can easily draw a conclusion that the disease-free steady state and β the threshold of model (3) are respectively P0 (x, y, z) = (1, 0, 0) and R0 = b+γ+ε . We also check that P0 (x, y, z) is globally asymptotically stable in Γ if R0 < 1, and P0 (x, y, z) is unstable and model (3) admits a unique endemic equilibrium P ∗ (x∗ , y ∗ , z ∗ ) if R0 > 1. For a more realistic and reasonable modelling, we introduce random perturbations into SIR epidemic model (1) for further analysis on persistence, extinction and stationary distribution. Now, model (1) becomes the following stochastic framework: βSI
dt + σ1 SdB1 (t), N  βSI  dI(t) = − (µ + γ + ε)I dt + σ2 IdB2 (t), N dR(t) = (γI − µR)dt + σ3 RdB3 (t), dS(t)

=

bN − µS −

(4)

where Bi (t)(i = 1, 2, 3) are independent Brownian motions, σi (i = 1, 2, 3) are the intensities of the white noises. The corresponding equation of the total population size N is governed by:   dN = (b − µ)N − εI dt + σ1 SdB1 (t) + σ2 IdB2 (t) + σ3 RdB3 (t).

2

Let us rewrite model (4) in terms of three proportion variables x, y, z:   dx = b − bx − βxy + εxy − σ12 x2 + x(σ12 x2 + σ22 y 2 + σ32 z 2 ) dt

+σ1 x(1 − x)dB1 (t) − σ2 xydB2 (t) − σ3 xzdB3 (t),   βxy − (b + γ + ε)y + εy 2 − σ22 y 2 + y(σ12 x2 + σ22 y 2 + σ32 z 2 ) dt

dy

=

dz

−σ1 xydB1 (t) + σ2 y(1 − y)dB2 (t) − σ3 yzdB3 (t),  = γy − bz + εyz − σ32 z 2 + z(σ12 x2 + σ22 y 2 + σ32 z 2 )dt

(5)

−σ1 xzdB1 (t) − σ2 yzdB2 (t) + σ3 z(1 − z)dB3 (t).

Due to the relationship x + y + z = 1, we need to discuss the equivalent model (6) as follows:   dy = β(1 − y − z)y − (b + γ + ε)y + εy 2 − σ22 y 2 + y(σ12 (1 − y − z)2 + σ22 y 2 + σ32 z 2 ) dt dz

−σ1 (1 − y − z)ydB1 (t) + σ2 y(1 − y)dB2 (t) − σ3 yzdB3 (t),  = γy − bz + εyz − σ32 z 2 + z(σ12 (1 − y − z)2 + σ22 y 2 + σ32 z 2 )dt

(6)

−σ1 (1 − y − z)zdB1 (t) − σ2 yzdB2 (t) + σ3 z(1 − z)dB3 (t),

with its initial value (y(0), z(0)) ∈ R2+ and y(0) + z(0) < 1. Throughout this paper, let (Ω, {Ft }t≥0 , P) be the complete probability space with a {Ft }t≥0 satisfying the usual conditions, that is, it is a right continuous function and F0 contains all P-null sets. The framework of our paper is organized as follows: we firstly prove the existence and uniqueness of a global positive solution. Then, we verify the extinction and persistence in the mean under some conditions. In the next part, we also prove the existence of the stationary distribution. Consequently, we carry out numerical simulations in support of our analytical results.

2

Existence and uniqueness of positive solution

Theorem 2.1 There is a unique solution (y(t), z(t)) to model (6) on t ≥ 0 for any given initial value (y(0), z(0)) ∈ R2+ and y(0) + z(0) < 1, and the solution will remain in R2+ with probability one, namely (y(t), z(t)) ∈ R2+ for all t ≥ 0 almost surely. Proof Since the coefficients of model (5) are locally Lipschitz continuous, there exists a unique local solution (x(t), y(t), z(t)) on t ∈ [0, τe ), where τe is the explosion time (Mao,2008). To show this assertion, we need to show that τe = ∞ a.s. Let m0 > 0 be sufficiently large such that each component of (x(0), y(0), z(0)) all lies within the interval [ m10 , m0 ]. For each integer m ≥ m0 , we define the stopping time  1 τm = inf t ∈ [0, τe ) : min{x(t), y(t), z(t)} ≤ m or max{x(t), y(t), z(t)} ≥ m . Throughout this paper, we set inf ∅ = ∞. It’s obvious that τm is an increasing function as m → ∞. We also set lim τm = τ∞ , then the fact that τ∞ ≤ τe is valid almost surely. If we can show that τ∞ = ∞ m→∞ a.s., then the desired result τe = ∞ holds for all t ≥ 0. The proof goes by contradiction. There exists a

3

pair of constants T > 0 and ε ∈ (0, 1) such that P{τm ≤ T } ≥ ε for each integer m ≥ m0 . Now we define a C 2 -function from R3+ to R+ : V (x, y, z)

= x − 1 − ln x + y − 1 − ln y + z − 1 − ln z = − ln(1 − y − z) − ln y − ln z − 2,

the generalized Itˆo’s formula gives:  (b + γ + ε)(y − x) ε(y − x)y σ22 y(y − x) σ32 (z − x)z dV (x, y, z) = β(y − x) − + − − x x x x (y − x)(σ12 x2 + σ22 y 2 + σ32 z 2 ) r(z − x)y b(z − x) εy(z − x) + + − + x xz x x (z − x)(σ12 x2 + σ22 y 2 + σ32 z 2 ) 1 2 2 2 2 + + σ1 (2x + y + z + 2yz) x 2  2  1 2 y (1 − y)2 − 2(1 − y)zy 2 + y 2 z 2 2 + 1 − 2y + 2y + σ2 2  x2  2 2 2 1 y z − 2yz + 2yz 3 + z 4 − 2z 3 + z 2 2 + 1 − 2z + 2z dt + σ32 2 x2 +(3x − 1)σ1 dB1 (t) + (3y − 1)σ2 dB2 (t) + (3z − 1)σ3 zdB3 (t),

(7)

where the items of having ε in (7)

εy 2 εyz −εy +ε+ − εy + − εy x x x

Therefore, dV (x, y, z) ≤

≤

−εy +ε+ x −εy = +ε+ x = ε − 3εy.

=

εy (y + z) − 2εy x εy (1 − x) − 2εy x



1 1 βy + 2b + γ + ε + σ12 (xy − x2 + xz − x2 + x2 + y 2 + z 2 + yz) 2 2  y3 zy 2 y 2 (1 − y)2 − 2(1 − y)zy 2 + y 2 z 2 y2 2 2 2 +σ2 − +y+ −y + −y + x  x x 2x2 2 2 3 1 yz z z + − y + y 2 + σ32 − z2 − +z+ − z2 2 x x x  y 2 z 2 − 2yz 2 + 2yz 3 + z 4 − 2z 3 + z 2 1 2 + + − z + z dt 2x2 2 +(3x − 1)σ1 dB1 (t) + (3y − 1)σ2 dB2 (t) + (3z − 1)σ3 zdB3 (t) 1 1 1
β + 2b + γ + ε + σ12 + σ22 + σ32 dt + (3x − 1)σ1 dB1 (t) 2 2 2 +(3y − 1)σ2 dB2 (t) + (3z − 1)σ3 zdB3 (t)

:= K0 dt + (3x − 1)σ1 dB1 (t) + (3y − 1)σ2 dB2 (t) + (3z − 1)σ3 zdB3 (t),

For any t ∈ [0, T ] and m ≥ m0 , integrating (7) from 0 to t ∧ τm and taking expectation gives that R τ ∧T EV (x(τm ∧ T ), y(τm ∧ T ), z(τm ∧ T )) ≤ V (x(0), y(0), z(0)) + E 0 m K0 dt ≤ V (x(0), y(0), z(0)) + K0 T < ∞.

We set Ωm = {τm ≤ t} for m ≥ m0 , then P(Ωm ) ≥ ε. Each component of (x(τm ∧T ), y(τm ∧T ), z(τm ∧T )) 1 equals either m or m for all ω ∈ Ωm . Consequently, we have ∞ > V (x(0), y(0), z(0)) + K0 T

1 ≥ P{τm ≤ T } min{m − 1 − ln m, m − 1 + ln m} 1 ≥ ε min{m − 1 − ln m, m − 1 + ln m},

4

Letting m → ∞ leads to the contradiction ∞ > V (x(0), y(0), z(0)) + K0 T = ∞. The proof is complete.

3

Extinction

Two of the most important issues are extinction and persistence when we studied epidemic models. Our aim of this section is to find out the sufficient conditions for the extinction of the disease, and to leave the argument of persistence to the next section. Theorem 3.1 Let (x(t), y(t), z(t)) be a solution of model (5) with initial value (x(0), y(0), z(0)) ∈ Γ. If the coefficients of model (5) satisfy σ12 + σ22 + σ32 ≤ 2(b + γ − β), then the infected individuals will tend to zero exponentially with probability one. Proof Making use of Itˆo’s formula to model (5) and the facts that x < 1, y < 1, z < 1, which yields to d ln y(t)

 1 1 βx − (b + γ + ε) + εy − σ22 y + σ12 x2 + σ22 y 2 + σ32 z 2 − σ12 x2 − σ22 (1 − y)2 2 2  1 − σ32 z 2 dt − σ1 xdB1 (t) + σ2 (1 − y)dB2 (t) − σ3 zdB3 (t) 2 1 1 1 1 ≤ β − (b + γ + ε) + ε + σ12 x2 + σ22 y 2 − σ22 + σ32 z 2 − σ1 xdB1 (t) 2 2 2 2 +σ2 (1 − y)dB2 (t) − σ3 zdB3 (t) 1 ≤ β − (b + γ) + (σ12 + σ22 + σ32 ) − σ1 xdB1 (t) + σ2 (1 − y)dB2 (t) − σ3 zdB3 (t). 2 =

Integrating (8) from 0 to t and dividing by t on both sides, we have ln y(t) 1 y(0) M1 (t) ≤ β − (b + γ) + (σ12 + σ22 + σ32 ) + + , t 2 t t where M1 (t) = −σ1

Rt 0

x(s)dB1 (s) + σ2

Rt 0

(1 − y(s))dB2 (s) − σ3

By strong law of large numbers for martingales lim

t→∞

M1 (t) = 0 a.s. t

Thus, lim sup t→∞

ln y(t) 1 ≤ β − (b + γ) + (σ12 + σ22 + σ32 ) < 0 a.s. t 2

That is to say, lim y(t) = 0 a.s. The proof is complete. t→∞

5

Rt 0

z(s)dB3 (s).

(8)

4

Persistence in the mean

Definition 4.1 (see [8]) In this part, we establish a condition of model (6) for persistence in the mean if Z Z 1 t 1 t lim inf y(s)ds > 0, lim inf z(s)ds > 0 a.s t→∞ t 0 t→∞ t 0 For simplicity and convenience, we introduce the notation: Z 1 t hx(t)i = x(s)ds. t 0

Theorem 4.1 Let (x(t), y(t), z(t)) be a solution of model (5) with initial value (x(0), y(0), z(0)) ∈ Γ. If βb 1 1 + ε − σ12 − σ22 , 2ε 2 2 then the following properties are valid 3b + γ <

lim inf hy(t)i ≥ y˜∗ , lim inf hz(t)i ≥ t→∞

t→∞

γ − k(b + γ + ε + σ22 ) ∗ y˜ > 0. a.s. b + σ32

where 0
γ H βb 1 1 > 0, H = − (3b + γ + σ12 + σ22 − ε). , y˜∗ = b + γ + ε + σ22 β 2ε 2 2

Proof By the similar discussion, applying the Itˆ o’s formula to model (5) leads to d(ln x + ln y)

=

1 b 1 − (2b + γ + ε + σ12 + σ22 ) + + βx + (2ε − β)y + σ12 x2 + σ22 y 2 2 2 x  +σ32 z 2 dt + (σ1 − 2σ1 x)dB1 (t) + (σ2 − 2σ2 y)dB2 (t) − 2σ3 zdB3 (t). 

Integrating (9) from 0 to t and dividing by t on both sides, which implies that ln x(t) ln y(t) + t t

1 1 1 = −(2b + γ + ε + σ12 + σ22 ) + bh i + βhxi + (2ε − β)hyi 2 2 x M2 (t) ln x(0) ln y(0) +hσ12 x2 + σ22 y 2 + σ32 z 2 i + + + t t t 1 2 1 2 z ≥ −(2b + γ + ε + σ1 + σ2 ) + 2εhyi − βhyi + bh i + βhxi 2 2 x M (t) ln x(0) ln y(0) 2 +hσ12 x2 + σ22 y 2 + σ32 z 2 i + + + . t t t Z t Z t Z t where M2 (t) = (σ1 − 2σ1 x(s))dB1 (s) + (σ2 − 2σ2 y(s))dB2 (s) − σ3 z(s)dB3 (s). We denote

0

0

0

z z F (x, y, z) = 2εhyi + bh i + βhxi = h2εy + b + βxi =: hLa (x, y, z)i. x x Using Lagrange’s Method of Multipliers, then z La (x, y, z) = 2εy + b + βx + λ(x + y + z − 1) x  βb  b(β − 2ε) βb ≥ min L(x, y, z) = 2ε 1 − 2 + + x+y+z=1 4ε 2ε 2ε βb = 2ε + − b. 2ε 6

(9)

Hence F (x, y, z) ≥

min hLa (x, y, z)i = 2ε +

x+y+z=1

βb − b. 2ε

Therefore, ln x(t) ln y(t) + t t

≥

1 1 βb −(2b + γ + ε + σ12 + σ22 ) + 2ε + −b 2 2 2ε M2 (t) ln x(0) ln y(0) −βhyi + + + . t t t

ln x(t) ln y(t) + ≥ H − βhyi + ϕ(t), where t t βb 1 1 1 βb H = 2ε + − b − (2b + γ + ε + σ12 + σ22 ) = − (3b + γ + σ12 + 2ε 2 2 2ε 2 M2 (t) x(0) y(0) + + . ϕ(t) = t t t It follows from the large number theory for martingales that lim ϕ(t) = 0. Due So we have

t→∞

−∞ < ln y(t) < 0, it implies that lim inf hy(t)i ≥ t→∞

1 2 σ − ε) > 0, 2 2

to −∞ < ln x(t) < 0,

H := y˜∗ > 0, a.s.. β

From model (6) and generalized Itˆo’s formula, we have  kdy + dz = kβxy − k(b + γ + ε)y + kεy 2 − kσ22 y 2 + ky(σ12 x2 + σ22 y 2 + σ32 z 2 ) + γy  −bz + εyz − σ32 z 2 + z(σ12 x2 + σ22 y 2 + σ32 z 2 ) dt − (kσ1 xy + σ1 xz)dB1 (t) +[kσ2 y(1 − y) − σ2 yz]dB2 (t) + [σ3 z(1 − z) − kσ3 yz]dB3 (t)   γ − k(b + γ + ε + σ22 ) y − (b + σ32 )z dt − (kσ1 xy + σ1 xz)dB1 (t)

≥

(10)

+[kσ2 y(1 − y) − σ2 yz]dB2 (t) + [σ3 z(1 − z) − kσ3 yz]dB3 (t).

Let 0 < k < that k where

γ , integrating (10) from 0 to t and dividing by t on both sides, which derives b + γ + ε + σ22

 M3 (t) y(t) − y(0) z(t) − z(0)  + ≥ γ − k(b + γ + ε + σ22 ) hyi − (b + σ32 )hzi + , t t t

M3 (t)

=

− +

Z

t

kσ1 x(s)y(s) + σ1 x(s)z(s) dB1 (s) +

0

Z

0

t




Z

t

0

  σ3 z(s)(1 − z(s)) − kσ3 y(s)z(s) dB3 (s).

  kσ2 y(s)(1 − y(s)) − σ2 y(s)z(s) dB2 (s)

According to strong law of large numbers for martingale, then lim t→∞ on both sides yields that lim inf hz(t)i ≥ t→∞

M3 (t) = 0. Taking the inferior limit t

γ − k(b + γ + ε + σ22 ) ∗ y˜ > 0. b + σ32

Due to the relationship x + y + z = 1 and the above discussion, we have shown the persistence in the mean of y(t) and z(t). The proof is complete.

7

5

Stationary distribution and ergodicity

Assumption A (see [7]) There exists a bounded domain U ⊂ El with regular boundary, such that: (A1) In the open domain U and some neighborhood thereof, the smallest eigenvalue of the diffusion matrix A(x) is bounded away from zero; (A2) If x ∈ El \U , the mean time τ at which a path issuing from x reaches the set U is finite, and sup Ex τ < ∞

x∈K

for every compact subset K ⊂ El . Lemma 5.1 (see [12]) If Assumption (A) holds, then the Markov process x(t) has a stationary distribution µ(·), and Z Z n o 1 T Px lim f (x(t))dt = f (x)µ(dx) = 1, T →∞ T 0 El

where f is a function integrable with respect to the measure µ. Remark 5.1 (i) To verify (A1), it is sufficient to prove that there exists a positive constant M > 0 such that (see Strang, 2006) l X

i,j=1

aij (x)ξi ξj ≥ M |ξ|2 , x ∈ U, ξ ∈ El .

(ii) To proof (A2), it suffices to show that there exists a non-negative C 2 -function V and a neighborhood U such that, LV is negative for any El \U (see Zhu and Yin, 2007). Theorem 5.1 Let (x(t), y(t), z(t)) be a solution of model (5) with initial value (x(0), y(0), z(0)) ∈ Γ. If p 1 3b + γ < 2 βb − ε − (σ12 + σ22 + σ32 ), 2 then there is a stationary distribution and the solution is ergodic. Proof We define V (x(t), y(t), z(t)) = − ln(x + y) − ln x − ln y − ln z := V1 + V2 + V3 + V4 . Generalized Itˆ o’s formula, we have !   1 1 b − bx + (ε − β)xy − σ12 x2 + x(σ12 x2 + σ22 y 2 + σ32 z 2 ) − − LV1 (x, y, z) = x+y x+y βxy − (b + γ + ε)y + (ε − σ22 )y 2 + y(σ12 x2 + σ22 y 2 + σ32 z 2 )    1 1 σ1 x(1 − x) −σ1 xy 1    (x + y)2 (x + y)2  + trace  −σ2 xy σ2 y(1 − y)    1 1 2 2 2 −σ3 xz −σ3 yz (x + y) (x + y) ! σ1 x(1 − x) −σ2 xy −σ3 xz −σ1 xy σ2 y(1 − y) −σ3 yz =

=

b σ 2 x2 (γ + ε)y σ2 y2 + b − εy + 1 + + 2 − (σ12 x2 + σ22 y 2 + σ32 z 2 ) x+y x+y x+y x+y σ 2 x2 z 2 σ2 y2 z2 σ 2 z 2 (x + y)2 + 2 + 3 + 1 2 2 2(x + y) 2(x + y) 2(x + y)2 2 2 −b + (γ + ε)y σ1 x + σ22 y 2 (σ 2 x2 + σ22 y 2 )z 2 b − εy + + + 1 x+y x+y 2(x + y)2 1 −(σ12 x2 + σ22 y 2 + σ32 z 2 ), 2

−

8

LV2 (x, y, z)

=

= LV3 (x, y, z)

=

LV4 (x, y, z) = −

 1 b − bx − βxy + εxy − σ12 x2 + x(σ12 x2 + σ22 y 2 + σ32 z 2 ) x  1  + 2 σ12 x2 (1 − x)2 + σ22 x2 y 2 + σ32 x2 z 2 2x 1 1 1 1 b − + b + βy − εy − ( σ12 x2 + σ22 y 2 + σ32 z 2 ) + σ12 , x 2 2 2 2 −

1 1 1 1 −βx + (b + γ + ε) − εy − ( σ12 x2 + σ22 y 2 + σ32 z 2 ) + σ22 , 2 2 2 2 γy 1 1 1 1 + b − εy − ( σ12 x2 + σ22 y 2 + σ32 z 2 ) + σ32 . z 2 2 2 2

Hence, LV (x, y, z)

= LV1 (x, y, z) + LV2 (x, y, z) + LV3 (x, y, z) + LV4 (x, y, z)

b γy 1 − βx + βy − 4εy 4b + γ + ε + (σ12 + σ22 + σ32 ) − − 2 x z (σ 2 x2 + σ22 y 2 )z 2 −b + (γ + ε)y σ12 x2 + σ22 y 2 + + 1 + x+y x+y 2(x + y)2 5 2 25 2 2 + − ( σ1 x σ2 y + 2σ32 z 2 ) 2 2 b γy 1 b − − − βx − 4εy + (4b + 2γ + 2ε + β + 2σ12 + 2σ22 + σ32 ) ≤ − x+y x z 2 b bz γy ≤ − − − − 4εy + K, x+y x z =

1 where K = 4b + 2γ + 2ε + β + 2σ12 + 2σ22 + σ32 . Let the region 2 D = {(x, y) : 0 < x < 1, 0 < y < 1, 0 < x + y < 1},

Uδ0 ,δ1 ,δ2 = {(x, y) ∈ D : δ0 ≤ x ≤ 1 − δ0 − δ2 , δ0 ≤ y ≤ 1 − δ0 − δ2 , δ1 ≤ x + y ≤ 1 − δ2 }, where δ0 , δ1 , δ2 always lie within the interval (0, 1) of having small values such that 2δ0 < δ1 , δ0 + δ2 <

1 1 δ22 5 , + − < −1, 2 1 − δ2 2(1 − δ2 )2 2

−b −bδ2 −γδ0 −b + K < −1, + K < −1, + K < −1, + K < −1. δ1 δ0 δ2 1 − δ2

Case 1 For all (x, y) ∈ Dδ10 ,δ1 ,δ2 = {(x, y) ∈ D : 0 < x + y < δ1 }, (11) gives that LV ≤

−b −b +K ≤ + K < −1. x+y δ1

Case 2 For all (x, y) ∈ Dδ20 ,δ1 ,δ2 = {(x, y) ∈ D : 0 < x < δ0 , δ1 ≤ x + y ≤ 1 − δ2 }, (11) derives that LV ≤

−bz −bδ2 +K ≤ + K < −1. x δ0

Case 3 For all (x, y) ∈ Dδ30 ,δ1 ,δ2 = {(x, y) ∈ D : δ0 ≤ y < 1, 1 − δ2 < x + y < 1}, (11) implies that LV ≤

−γy −γδ0 +K ≤ + K < −1. z δ2 9

(11)

Case 4 For all (x, y) ∈ Dδ40 ,δ1 ,δ2 = {(x, y) ∈ D : 0 < y < δ0 , δ1 < x + y < 1 − δ2 }, (11) derives that LV ≤

−b −b +K ≤ + K < −1. x+y 1 − δ2

Case 5 For all (x, y) ∈ Dδ50 ,δ1 ,δ2 = {(x, y) ∈ D : 0 < y < δ0 , 1 − δ2 ≤ x + y < 1}, (11) yields that 1 b γy b (γ + ε)y 4b + γ + ε + (σ12 + σ22 + σ32 ) − − − βx + βy − 4εy − + 2 x z x+y x+y 5 2 2 5 2 2 (σ12 x2 + σ22 y 2 )z 2 σ12 x2 + σ22 y 2 − ( σ1 x + σ2 y + 2σ32 z 2 ) + + x+y 2(x + y)2 2 2 2 2 2 2 b γy (γ + ε)δ0 σ x + σ2 y (σ 2 x2 + σ22 y 2 )δ22 ≤ − − βx − − 4εy + + 1 + 1 x z 1 − δ2 1 − δ2 2(1 − δ22 )2 5 5 1 +βδ0 − ( σ12 x2 + σ22 y 2 + 2σ32 z 2 ) + 3b + γ + ε + (σ12 + σ22 + σ32 ) 2 2 2 p  (γ + ε)δ0  1 5 δ22 2 2 2 2 − + βδ0 ≤ −2 βb + + (σ1 x + σ2 y ) + 2 1 − δ2 1 − δ2 2(1 − δ2 ) 2 1 2 + (σ1 + σ22 + σ32 ) + 3b + γ + ε 2 p 1 ≤ −2 βb + (σ12 + σ22 + σ32 ) + 3b + γ + ε < 0. 2 √ We denote −ψ = max{−1, −2 βb + 12 (σ12 + σ22 + σ32 ) + 3b + γ + ε} < 0, and take U = Uδ0 ,δ1 ,δ2 , for ¯ , we can get LV ≤ −ψ < 0, which proves the condition (A2) of Assumption A. arbitrary x ∈ D\U For any (x, y, z) ∈ D, there is a positive number  M = min (σ12 (1 − x)2 + σ22 y 2 + σ32 z 2 )x2 , (σ12 x2 + σ22 (1 − y)2 + σ32 z 2 )y 2 , (σ12 x2 + σ22 y 2 + σ32 (1 − z)2 )z 2 > 0 LV (x + y + z)

=

such that for any ξ ∈ R3 , 3 X

i,j=1

aij ξi ξj

= (σ12 (1 − x)2 + σ22 y 2 + σ32 z 2 )x2 ξ12 + (σ12 x2 + σ22 (1 − y)2 + σ32 z 2 )y 2 ξ22 +(σ12 x2 + σ22 y 2 + σ32 (1 − z)2 )z 2 ξ32 ≥ M kξk2 .

We also verify the condition (A1) of Assumption A. Hence, the stochastic model (5) has a stable stationary distribution and the solution is ergodic. The proof is complete. Remark 5.2 We can draw an interesting conclusion after comparison of the conditions of Theorem 4.1 and Theorem 5.1. In other words, if model (5) has a stationary distribution, then the solution of model (5) is also persistent in the mean. According to the following fact p βb 1 p 1 2 βb − = − ( βb − 2ε)2 + 2ε < 2ε + σ32 2ε 2ε 2

is true without any further proof, we therefore have

p 1 βb 1 1 + ε − σ12 − σ22 . 2 βb − ε − (σ12 + σ22 + σ32 ) < 2 2ε 2 2

On the contrary, if model (5) is just persistent, we can not draw such a conclusion in which model (5) admits a stationary distribution.

10

6

Examples and simulations

Example 6.1 We assume that the parameters of model (5) are listed in the Table 1: Parameters Values

σ1 0.1

σ2 0.15

σ3 0.2

b 0.1

γ 0.1

β 0.15

ε 0.15

and the initial value is (x(0), y(0), z(0)) = (0.2, 0.7, 0.1). Here 0.0725 = σ12 + σ22 + σ32 < 2(b + γ − ε) = 0.1. According to Theorem 3.1, the condition is satisfied, thus the infected individuals tend to extinction (see Figure 1).

1 0.9 0.8 0.7

y(t)

0.6 0.5 0.4 0.3 0.2 0.1 0

0

0.5

1 time t

11

1.5

2 4

x 10

0.8 0.7 0.6

y(t)

0.5 0.4 0.3 0.2 0.1 0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x(t)

Figure 1: The extinction of the infected individuals of model (5). Example 6.2 Let the parameters of model (5) be listed in the Table 2: Parameters Values

σ1 0.01

σ2 0.008

σ3 0.02

b 0.04

γ 0.01

β 0.5

ε 0.15

and the initial value is (x(0), y(0), z(0)) = (0.5, 0.2, 0.3). Here 0.13 = 3b + γ <

βb 1 1 + ε − σ12 − σ22 = 0.21658. 2ε 2 2

According to Theorem 4.1, the condition is valid, thus densities of the susceptible and infected individuals approach persistence in the mean (see Figure 2). Example 6.3 We suppose that the parameters of model (5) are listed in the Table 3: Parameters Values

σ1 0.01

σ2 0.008

σ3 0.02

b 0.04

γ 0.01

β 0.5

ε 0.15

and the initial value is (x(0), y(0), z(0)) = (0.5, 0.2, 0.3). Here p 1 βb 1 1 + ε − σ12 − σ22 . 0.13 = 3b + γ < 2 βb − ε − (σ12 + σ22 + σ32 ) = 0.13256 < 0.21658 = 2 2ε 2 2

From Theorem 5.1 and Assumption A, the condition is satisfied, thus model (5) admits a stationary distribution and it is ergodic. Further, the solution of model (5) is persistent in the mean (see Figure 3).

12

0.5 x y

0.45

population densities of x and y

0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0

0

2

4

6 time t

8

10 4

x 10

Figure 2: The persistence in the mean of the susceptible individuals and the infected individuals of model (5).

Acknowledgments The authors would like to thank the reviewers for their good comments. This work is supported by the National Natural Science Foundation of China (Grant No. 11201075), Natural Science Foundation of Fujian Province of China (Grant No. 2016J01015) and Scholarship under the Education Department of Fujian Province.

References [1] Kermack W.O., McKendrick A.G., Contributions to the mathematical theory of epidemics (Part I) Proceedings of the Royal Society A, 1927, 115: 700-721. [2] Zaman G., Kang Y.H., Jung I.H., Stability analysis and optimal vaccination of an SIR epidemic model[J]. BioSystems, 2008, 93(3): 240-249. [3] Gray A., Greenhalgh D., Hu L. et al., A stochastic differential equation SIS epidemic model[J]. SIAM Journal on Applied Mathematics, 2011, 71(3): 876-902. [4] Zhao Y., Jiang D., Mao X.R., et al., The threshold of a stochastic SIRS epidemic model in a population with varying size[J]. Discrete Continuous Dynamical Systems-Series B, 2015, 20: 12771295. [5] Zhao D., Study on the threshold of a stochastic SIR epidemic model and its extensions[J]. Communications in Nonlinear Science and Numerical Simulation, 2016, 38: 172-177.

13

The histogram of x(t) 600

500

Frequency

400

300

200

100

0

0.28

0.3

0.32 0.34 x at time 4000

0.36

0.38

The histogram of y(t)

250

Frequency

200

150

100

50

0

0.18

0.2

0.22 0.24 y at time 4000

14

0.26

0.28

The histogram of z(t)

250

Frequency

200

150

100

50

0 0.25

0.3

0.35

0.4 0.45 z at time 4000

0.5

0.55

0.6

Figure 3: The histogram of the susceptible individuals, the infected individuals and the removed individuals of model (5). [6] Liu Q., Chen Q., Dynamics of a stochastic SIR epidemic model with saturated incidence[J]. Applied Mathematics and Computation, 2016, 282: 155-166. [7] Liu J.M., Wei F.Y., Dynamics of stochastic SEIS epidemic model with varying population size[J]. Physica A: Statistical Mechanics and its Applications, 2016, 464: 241-250. [8] Zhao Y.N., Jiang D.Q., The threshold of a stochastic SIRS epidemic model with saturated incidence[J]. Applied Mathematics Letters, 2014, 34: 90-93. [9] Mao X.R., Stochastic Differential Equations and Applications (2nd Edition)[M]. Horwood, Chichester, 2008. [10] Zhou J., Hethcote W.H., Population size dependent incidence in models for diseases without immunity[J]. Journal of Mathematical Biology, 199, 32(8): 809-834. [11] Busenberg S., Driessche P. van den, Analysis of a disease transmission model in a population with varying size[J]. Journal of mathematical biology, 1990, 28(3): 257-270. [12] Hasminskii R.Z., Stochastic Stability of Differential Equations[M]. Netherlands, Sijithoff and Noordhoff, Alphen aan den Rijin, 1980. [13] Strang G., Linear Algebra and its Applications[M]. America, Thomson Learning Inc, 2006. [14] Zhu C., Yin G., Asymptotic properties of hybrid diffusion systems[J]. SIAM Journal on Control and Optimization, 2007, 46(4): 1155-1179.

15

The highlights of this paper are presented below: • A stochastic susceptible-infected-removed model is formulated when population size varies with time. • Extinction of solution is derived when intensities of white noises are controlled by parameters of model. • Persistence, stationary distribution and ergodicity are followed by careful computation via constructing suitable functions. • Three examples and numerical simulations are separately carried out to check the validity of the main results.