Threshold behavior in a stochastic SIR epidemic model with Logistic birth

Physica A xxx (xxxx) xxx

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

Threshold behavior in a stochastic SIR epidemic model with Logistic birth Qun Liu a , Daqing Jiang b,c,d ,

∗

a

School of Mathematics and Statistics, Key Laboratory of Applied Statistics of MOE, Northeast Normal University, Changchun 130024, Jilin Province, PR China b Key Laboratory of Unconventional Oil and Gas Development, China University of Petroleum (East China), Ministry of Education, Qingdao 266580, PR China c College of Science, China University of Petroleum, Qingdao 266580, Shandong Province, PR China d Nonlinear Analysis and Applied Mathematics (NAAM)-Research Group, King Abdulaziz University, Jeddah, Saudi Arabia

article

info

Article history: Received 17 April 2019 Received in revised form 18 October 2019 Available online xxxx Keywords: SIR epidemic model Threshold Stationary distribution Extinction Logistic birth

a b s t r a c t In this paper, we consider a stochastic SIR epidemic model with Logistic birth. By using the stochastic Lyapunov function method, we show that the stochastic basic reproduction number RS0 can be used to determine the threshold dynamics of the stochastic system. If RS0 > 1, we establish sufficient conditions for the existence of a stationary distribution of the positive solutions to the model. While if RS0 < 1, under some extra conditions, we obtain sufficient conditions for extinction of the disease. Finally, some examples and numerical simulations are provided to illustrate the theoretical results. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Since the pioneer work of Kermack and McKendrick [1], mathematical models have been widely used to analyze the spread and control of infectious diseases. During the past few decades, many kinds of epidemic models have been developed [2–13]. For example, Ma et al. [2] considered an SIRS epidemic model with standard incidence. Their work suggests that the basic reproduction number R0 determining whether there exists an endemic equilibrium or not: if R0 < 1 the disease will die out and the entire population will be susceptible; while if R0 > 1, under mild extra conditions, the disease-free equilibrium is unstable and there is a unique positive endemic equilibrium which is globally asymptotically stable. Zhang et al. [3] proposed an SIR epidemic model with Logistic birth which takes on the following form

⎧ ( ) dS N ⎪ ⎪ ⎪ = b − r N − β SI − µS , ⎪ ⎪ K ⎪ ⎨ dt dI

= β SI − δ I , ⎪ ⎪ dt ⎪ ⎪ ⎪ ⎪ dR ⎩ = γ I − µR,

(1.1)

dt where N(t) = S(t) + I(t) + R(t) denotes the total population size at time t. S(t), I(t), R(t) are the numbers of the susceptible individuals, the infectious individuals and the removed individuals at time t, respectively. The parameter β (β > 0) ∗ Corresponding author at: College of Science, China University of Petroleum, Qingdao 266580, Shandong Province, PR China. E-mail address: [email protected] (D. Jiang). https://doi.org/10.1016/j.physa.2019.123488 0378-4371/© 2019 Elsevier B.V. All rights reserved.

Please cite this article as: Q. Liu and D. Jiang, Threshold behavior in a stochastic SIR epidemic model with Logistic birth, Physica A (2019) 123488, https://doi.org/10.1016/j.physa.2019.123488.

2

Q. Liu and D. Jiang / Physica A xxx (xxxx) xxx

denotes the transmission coefficient, b denotes the birth rate, µ is the natural death rate, r = b − µ is the intrinsic rate, δ = α +γ +µ and α is a nonnegative constant representing the disease-caused death rate and γ > 0 denotes the recovery rate of I. They used N as a variable in place of S, then the SIR epidemic model can be expressed as follows

⎧ dI ⎪ ⎪ ⎪ = β I(N − I − R) − δ I , ⎪ ⎪ dt ⎪ ⎨ dR = γ I − µR, ⎪ dt ⎪ ( ) ⎪ ⎪ dN N ⎪ ⎪ ⎩ =r 1− N − αI . dt

(1.2)

K

They proved that the region Γ = {(I , R, N) ∈ R3+ : I + R ≤ N ≤ K } is a positively invariant set of system (1.2). The diseaseβK free equilibrium E0 = (0, 0, K ) of system (1.2) always exists, and if R0 = δ ≤ 1, E0 is globally asymptotically stable in Γ . If R0 > 1, E0 is unstable and there exists a unique endemic equilibrium E ∗ (I ∗ , R∗ , N ∗ ) which is locally asymptotically stable and globally asymptotically stable if α ≤ min{2µ, r /2}. As we know, stochastic epidemic models were advanced in the early 20th century. Early in 1926, McKendrick [14] proposed a stochastic version regarding the general epidemic model. From then on, various forms of stochastic epidemic models have been formulated and studied by many authors [15–20]. For example, May [15] revealed that due to environmental noise, the structures entailed the system of the model to display indiscriminate variations. Lin and Jiang [16] studied a stochastic SIR model which was perturbed by the white noise. They established sufficient conditions for the disease to die out exponentially, existence of a stationary distribution and asymptotic stability. Tornatore et al. [17] suggested a stochastic SIR model with or without allocated delay in time. They gave a satisfactory situation for the asymptotic constancy of the disease-free equilibrium. Gray et al. [18] considered a stochastic SIS epidemic model with constant population size. They obtained conditions for extinction and persistence of the disease. They also proved the existence of a stationary distribution and obtained expressions for its mean and variance. As indicated in [21–25], due to the existence of environmental noise, the parameters involved in system (1.1) will fluctuate around some average values and do not attain fixed values as time evolves [26]. Note that the transmission coefficient is one of the crucial parameters in the epidemic model, and it is more affected by the environmental noise. Then in this paper, we introduce randomness into system (1.1) by perturbing the transmission coefficient β by β + σ ξ (t) to obtain the following stochastic differential equations:

⎧ ( ) N dS ⎪ ⎪ ⎪ = b − r N − (β + σ ξ (t))SI − µS , ⎪ ⎪ K ⎪ ⎨ dt dI

= (β + σ ξ (t))SI − δ I , ⎪ ⎪ dt ⎪ ⎪ ⎪ ⎪ ⎩ dR = γ I − µR,

(1.3)

dt

where ξ (t) denotes a Gaussian white noise and characterized by:

⟨ξ (t)⟩ = 0, ⟨ξ (t)ξ (t ′ )⟩ = δ (t − t ′ ), and ⟨·⟩ denotes ensemble average. Moreover, δ (·) is the Dirac-δ function, σ 2 > 0 denotes the intensity of environmental forcing. Now we can rewrite system (1.3) into the following form of stochastic differential equations:

[( ) ] ⎧ N ⎪ ⎪ N − β SI − µS dt − σ SIdB(t), dS = b−r ⎪ ⎨ K dI = (β SI − δ I)dt + σ SIdB(t), ⎪ ⎪ ⎪ ⎩ dR = (γ I − µR)dt ,

(1.4)

where B(t) is a standard one-dimensional Brownian motion. The relations between the white noise terms and Brownian motion are defined by dB(t) = ξ (t)dt. It is worthy noting that the SDE system (1.4) is the limiting form of the following discrete-time system (1.5). For 0 < h ≪ 1, consider

⎧ [( ) ] √ N(t) ⎪ ⎪ N(t) − β S(t)I(t) − µS(t) dt − σ hS(t)I(t)z(t), ⎨S(t + h) = S(t) + h · b − r K √ ⎪ I(t + h) = I(t) + h · [β S(t)I(t) − δ I(t)]dt + σ hS(t)I(t)z(t), ⎪ ⎩ R(t + h) = R(t) + h · [γ I(t) − µR(t)]dt ,

(1.5)

where for some constant d > 0, {z(t)} is the Bernoulli noise process

P(z(t) = d) = P(z(t) = −d) = 0.5. Please cite this article as: Q. Liu and D. Jiang, Threshold behavior in a stochastic SIR epidemic model with Logistic birth, Physica A (2019) 123488, https://doi.org/10.1016/j.physa.2019.123488.

Q. Liu and D. Jiang / Physica A xxx (xxxx) xxx

3

This system is an appropriate simplification as, in the limit of h → 0, it tends to an Itô equation of the form (1.4). For more details, we refer the readers to [27]. The main aim of this paper is to investigate how environmental fluctuations affect the dynamics of the disease through studying the global dynamics of a stochastic SIR epidemic model with Logistic birth. The organization of this paper is as follows: in Section 2, we present some preliminaries which will be used in our following analysis. In Section 3, we show that there exists a unique global positive solution to system (1.4) with any positive initial value. In Section 4, we establish sufficient conditions for the existence of a stationary distribution of the positive solutions to system (1.4). In Section 5, we obtain sufficient conditions for extinction of the disease. Some examples and numerical simulations are provided to demonstrate our analytical results in Section 6. Finally, some conclusion and future directions are presented to end this paper. 2. Preliminaries In this section, we will present some known results and lemmas which will be used later. Throughout this paper, unless otherwise specified, let (Ω , F , {Ft }t ≥0 , P) be a complete probability space with a filtration {Ft }t ≥0 satisfying the usual conditions (i.e., it is increasing and right continuous while F0 contains all P-null sets) and let B(t) defined on the above complete probability space. Define

Rd+ = {x = (x1 , x2 , . . . , xd ) ∈ Rd : xi > 0, 1 ≤ i ≤ d}. Then we give some basic theories in stochastic differential equations which is introduced in [28]. In general, consider the d-dimensional stochastic differential equation dX (t) = f (X (t), t)dt + g(X (t), t)dB(t) for t ≥ t0 ,

(2.1)

with initial value X (0) = X0 ∈ Rd . B(t) denotes a d-dimensional standard Brownian motion defined on the complete probability space (Ω , F , {Ft }t ≥0 , P). Denote by C 2,1 (Rd × [t0 , ∞); R+ ) the family of all nonnegative functions V (X , t) defined on Rd × [t0 , ∞) such that they are continuously twice differentiable in X and once in t. The differential operator L of Eq. (2.1) is defined by [28] L=

d d ∑ ∂ ∂2 ∂ 1∑ T + fi (X , t) + [g (X , t)g(X , t)]ij . ∂t ∂ Xi 2 ∂ Xi ∂ Xj i,j=1

i=1

If L acts on a function V ∈ C 2,1 (Rd × [t0 , ∞); R+ ), then LV (X , t) = Vt (X , t) + VX (X , t)f (X , t) + where Vt =

∂V ∂t

, VX = ( ∂∂XV , . . . , 1

∂V ), ∂ Xd

1 2

trace[g T (X , t)VXX (X , t)g(X , t)],

VXX = ( ∂ ∂X ∂VX )d×d . In view of Itô’s formula [28], if X (t) ∈ Rd , then 2 i

j

dV (X (t), t) = LV (X (t), t)dt + VX (X (t), t)g(X (t), t)dB(t). Consider the integral equation

∫

t

b(s, X (s))ds +

X (t) = X (t0 ) + t0

k ∫ ∑ r =1

t

σr (s, X (s))dBr (s).

(2.2)

t0

Lemma 2.1 ([29]). Suppose that the coefficients of (2.2) are independent of t and satisfy the following conditions for some constant B:

|b(s, x) − b(s, y)| +

k ∑

|σr (s, x) − σr (s, y)| ≤ B|x − y|, |b(s, x)| +

r =1

k ∑

|σr (s, x)| ≤ B(1 + |x|),

(2.3)

r =1

in DR for every R > 0, and that there exists a nonnegative C 2 -function V (x) in Rd such that LV (x) ≤ −1 outside some compact set. Then there exists a solution of system (2.2) which is a stationary Markov process. Remark 2.1. The condition (2.3) in Lemma 2.1 can be replaced by the global existence of the solutions to (2.2) according to Remark 5 of Xu et al. [30]. Please cite this article as: Q. Liu and D. Jiang, Threshold behavior in a stochastic SIR epidemic model with Logistic birth, Physica A (2019) 123488, https://doi.org/10.1016/j.physa.2019.123488.

4

Q. Liu and D. Jiang / Physica A xxx (xxxx) xxx

3. Existence and uniqueness of the global positive solution In this section, we show that there exists a unique global positive solution to system (1.4) which is a prerequisite for analyzing the long-time behavior of system (1.4). Theorem 3.1. For any initial value (S(0), I(0), R(0)) ∈ R3+ , there is a unique solution (S(t), I(t), R(t)) to system (1.4) on t ≥ 0 and the solution will remain in R3+ with probability one, namely, (S(t), I(t), R(t)) ∈ R3+ for all t ≥ 0 almost surely (a.s.) Proof. Since the coefficients of system (1.4) satisfy the local Lipschitz condition, then for any initial value (S(0), I(0), R(0)) ∈ R3+ there exists a unique local solution (S(t), I(t), R(t)) ∈ R3+ on t ∈ [0, τe ), where τe denotes the explosion time [28]. To show this solution is global, we only need to prove τe = ∞ a.s. To this end, let n0 > 0 be sufficiently large such that S(0), I(0) and R(0) all lie within the interval [ n1 , n0 ]. For each integer n ≥ n0 , define the stopping time as in [28] 0

{

τn = inf t ∈ [0, τe ) : min{S(t), I(t), R(t)} ≤

}

1

or max{S(t), I(t), R(t)} ≥ n ,

n

where throughout this paper we set inf ∅ = ∞ (as usual ∅ denotes the empty set). Obviously, τn is increasing as n → ∞. Set τ∞ = limn→∞ τn , whence τ∞ ≤ τe a.s. If τ∞ = ∞ a.s. is true, then τe = ∞ a.s. and (S(t), I(t), R(t)) ∈ R3+ a.s. for all t ≥ 0. That is to say, to finish the proof we only need to show τ∞ = ∞ a.s. If this assertion is not true, then there is a pair of constants T > 0 and ϵ ∈ (0, 1) such that

P{τ∞ ≤ T } > ϵ. Therefore, there exists an integer n1 ≥ n0 such that

P{τn ≤ T } ≥ ϵ, ∀n ≥ n1 . Moreover, for t ≤ τn , we can obtain that d(S + I + R) =dN

[( =

b−r

)

N K

]

N − µN − α I dt

) ] [( N N − µN dt ≤ b−r K ( ) N =r 1 − Ndt K

and so N(t) ≤ K , if N(0) ≤ K . Define a C 2 -function V : R3+ → R+ ∪ {0} by V (S , I , R) = (S − 1 − ln S) + (I − 1 − ln I) + (R − 1 − ln R). The nonnegativity of this function can be seen from u − 1 − ln u ≥ 0, ∀u > 0. Let n ≥ n0 and T > 0 be arbitrary. Applying Itô’s formula to V , we obtain dV (S , I , R) = LV (S , I , R)dt + σ (I − S)dB(t), where LV : R3+ → R is defined by

(

LV (S , I , R) = 1 −

1

)[(

S

b−r

N K

)

]

N − β SI − µS +

σ2 2

2

I +

( 1−

1 I

)

(β SI − δ I) +

σ2 2

)

K

(

2

1−

S +

N σ2 2 σ2 2 N N =N b − r − µN − α I + 2µ + δ + β I + I + S − βS − b−r K 2 2 S K ( ) N σ2 2 σ2 2 ≤r 1 − N + 2µ + δ + β I + I + S K 2 2 ( ) N σ2 2 σ2 2 ≤r 1 − N + 2µ + δ + β K + K + K

(

(

2

) −

1 R

)

(γ I − µR)

γI R

2

:=K1 . Here K1 is a positive constant. The rest of the proof follows that of Dalal et al. [31] and hence is omitted. This completes the proof. Please cite this article as: Q. Liu and D. Jiang, Threshold behavior in a stochastic SIR epidemic model with Logistic birth, Physica A (2019) 123488, https://doi.org/10.1016/j.physa.2019.123488.

Q. Liu and D. Jiang / Physica A xxx (xxxx) xxx

5

Remark 3.1. Theorem 3.1 shows that for any initial value (S(0), I(0), R(0)) ∈ R3+ , there exists a unique global solution (S(t), I(t), R(t)) ∈ R3+ a.s. to system (1.4). So d(S + I + R) =dN

) ] [( N N − µN − α I dt = b−r K [( ) ] N ≤ b−r N − µN dt K ( ) N =r 1 − Ndt . K

If N(0) < K , then N(t) < K a.s. Therefore, the region

Γ ∗ = {(S , I , R) ∈ R3+ : S + I + R < K } is a positively invariant set of system (1.4) on Γ ∗ . From now on, we always assume that the initial value (S(0), I(0), R(0)) ∈ Γ ∗ . 4. Existence of stationary distribution In this section, we will establish sufficient conditions for the existence of a stationary distribution of the positive solutions to system (1.4). The existence of a stationary distribution can be regarded as the stochastic weak stability of the system and means that the disease will prevail in the long term. βK δ+ 12 σ 2 K 2

Theorem 4.1. If RS0 := Markov process.

> 1, then there exists a positive solution (S(t), I(t), R(t)) to system (1.4) which is a stationary

Proof. By Lemma 2.1, we only need to construct a nonnegative C 2 -function V (S , I , R) and a closed set D ⊂ Γ ∗ such that LV (S , I , R) ≤ −1 for any (S , I , R) ∈ Γ ∗ \ D. In view of system (1.4), we have

( L −

I N

) =

β SI N

−

≤β I − ≤β I −

δI

δI N δI N

N

−

rI(1 −

N ) K

N

+

αI

+

αI

+

αI 2 N2

2

N2

(4.1)

N I

=β I − (µ + γ ) , N

( L(− ln N) = − r 1 −

N K

) +

αI N

.

(4.2)

From (4.1) and (4.2) it follows that

( L − ln N −

α I µ+γ N

)

( ) N ≤−r 1− + K

αβ I. µ+γ

(4.3)

Moreover, we have

( L − ln I +

) β 1 βγ R = − β S + δ + σ 2S2 + I − βR µ 2 µ βγ 1 ≤ − β S + δ + σ 2K 2 + I − βR 2 µ 1

= − β (N − I − R − K ) − β K + δ + σ 2 K 2 + 2

βγ I − βR µ

(4.4)

β (µ + γ ) 1 = − β (N − K ) + I − β K + δ + σ 2K 2. µ 2 Define V1 (S , I , R) = − ln I +

( ) β βK α I β r(µ + γ )2 + αβ 2 µK R− ln N + + I, µ r µ+γ N µr δ (µ + γ )

Please cite this article as: Q. Liu and D. Jiang, Threshold behavior in a stochastic SIR epidemic model with Logistic birth, Physica A (2019) 123488, https://doi.org/10.1016/j.physa.2019.123488.

6

Q. Liu and D. Jiang / Physica A xxx (xxxx) xxx

then in view of (4.3) and (4.4), we get

β (µ + γ ) αβ 2 K β r(µ + γ )2 + αβ 2 µK I+ I+ (β SI − δ I) 2 µ r(µ + γ ) µ r δ (µ + γ ) 1 β 2 r(µ + γ )2 + αβ 3 µK = − β K + δ + σ 2K 2 + SI 2 µr δ (µ + γ ) ) ( β 2 r(µ + γ )2 + αβ 3 µK 1 SI , = − δ + σ 2 K 2 (RS0 − 1) + 2 µr δ (µ + γ )

LV1 = − β K + δ +

1

σ 2K 2 +

(4.5)

where RS0 :=

βK . δ + 21 σ 2 K 2

Next, define V2 (S) = − ln S , V3 (S , I , R) = − ln(K − N), V4 (S , I , R) = − ln R, then making use of Itô’s formula, we have LV2 = −

(

N

b−r

S

2

rN(1 − r

2

N ) K

+ βI + µ +

σ2 2

I2 (4.6)

K2

σ2

≤β K + µ +

LV3 =

K

σ2

≤β I + µ +

)

N

K 2,

− αI

K −N αI

(4.7)

= N− K

K −N αI ≤− + r, K −N

LV4 = −

γI R

+ µ.

(4.8)

Define a C 2 -function Q : Γ ∗ → R as follows Q (S , I , R) = MV1 (S , I , R) + V2 (S) + V3 (S , I , R) + V4 (S , I , R) + V5 (R), where M > 0 is a sufficiently large number satisfying the following condition

(

1

−M δ + σ K 2

2

)

2

(RS0 − 1) + β K + 2µ + r +

σ2 2

K 2 ≤ −2.

(4.9)

In addition, Q (S , I , R) is not only continuous, but also tends to ∞ as (S , I , R) approaches the boundary of Γ ∗ . Therefore, it must be lower bounded and achieve this lower bound at a point (S0 , I0 , R0 ) in the interior of Γ ∗ . Then we define a C 2 -function V : Γ ∗ → R+ ∪ {0} as follows V (S , I , R) =Q (S , I , R) − Q (S0 , I0 , R0 )

=MV1 (S , I , R) + V2 (S) + V3 (S , I , R) + V4 (S , I , R) + V5 (R) − Q (S0 , I0 , R0 ). Then in view of (4.5)–(4.8), we have

(

LV ≤ − M δ +

1 2

σ K 2

2

)

(RS0 − 1) +

M(β 2 r(µ + γ )2 + αβ 3 µK )

µ r δ (µ + γ )

SI −

αI γI σ2 2 − + β K + 2µ + r + K . K −N R 2 (4.10)

Now we are in the position to construct a compact subset Dϵ such that the condition A2 in Lemma 2.1 holds. Define the following bounded closed set Dϵ = {(S , I , R) ∈ Γ ∗ : ϵ ≤ S , ϵ ≤ I , ϵ 2 ≤ R, S + I + R ≤ K − ϵ 2 }, Please cite this article as: Q. Liu and D. Jiang, Threshold behavior in a stochastic SIR epidemic model with Logistic birth, Physica A (2019) 123488, https://doi.org/10.1016/j.physa.2019.123488.

Q. Liu and D. Jiang / Physica A xxx (xxxx) xxx

7

where 0 < ϵ < 1 is a sufficiently small constant. In the set Γ ∗ \ Dϵ , we can choose ϵ sufficiently small such that the following conditions hold

ϵ≤

µr δ (µ + γ ) MK (β 2 r(µ + γ )2 + αβ 3 µK )

,

(4.11)

−

MK 2 (β 2 r(µ + γ )2 + αβ 3 µK ) σ2 2 γ + + β K + 2µ + r + K ≤ −1, ϵ µ r δ (µ + γ ) 2

(4.12)

−

α MK 2 (β 2 r(µ + γ )2 + αβ 3 µK ) σ2 2 + + β K + 2µ + r + K ≤ −1. ϵ µr δ (µ + γ ) 2

(4.13)

For convenience, we can divide Γ ∗ \ Dϵ into four domains, D1 = {(S , I , R) ∈ Γ ∗ : S < ϵ}, D2 = {(S , I , R) ∈ Γ ∗ : I < ϵ}, D3 = {(S , I , R) ∈ Γ ∗ : I ≥ ϵ, R < ϵ 2 }, D4 = {(S , I , R) ∈ Γ ∗ : K − ϵ 2 < S + I + R, I ≥ ϵ}. Clearly, Dcϵ = D1 ∪ D2 ∪ D3 ∪ D4 . Next, we will prove that LV (S , I , R) ≤ −1 for any (S , I , R) ∈ Dcϵ , which is equivalent to proving it on the above four domains, respectively. Case 1. For any (S , I , R) ∈ D1 , by (4.10), we have

) M(β 2 r(µ + γ )2 + αβ 3 µK ) σ2 2 σ 2 K 2 (RS0 − 1) + SI + β K + 2µ + r + K 2 µr δ (µ + γ ) 2 ) ( MK (β 2 r(µ + γ )2 + αβ 3 µK ) σ2 2 1 S + β K + 2µ + r + K ≤ − M δ + σ 2 K 2 (RS0 − 1) + 2 µr δ (µ + γ ) 2 ( ) 1 MK (β 2 r(µ + γ )2 + αβ 3 µK ) σ2 2 ≤ − M δ + σ 2 K 2 (RS0 − 1) + ϵ + β K + 2µ + r + K 2 µr δ (µ + γ ) 2 ≤−2+1 = − 1, (

LV ≤ − M δ +

1

(4.14)

which follows from (4.9) and (4.11). Therefore LV ≤ −1 for any (S , I , R) ∈ D1 . Case 2. For any (S , I , R) ∈ D2 , according to (4.10), we obtain

(

LV ≤ − M δ +

1

σ K 2

2

)

(RS0 − 1) +

M(β 2 r(µ + γ )2 + αβ 3 µK )

SI + β K + 2µ + r +

σ2

K2

µr δ (µ + γ ) 2 ) ( 2 2 3 MK (β r(µ + γ ) + αβ µK ) 1 σ2 2 ≤ − M δ + σ 2 K 2 (RS0 − 1) + I + β K + 2µ + r + K 2 µr δ (µ + γ ) 2 ( ) 1 MK (β 2 r(µ + γ )2 + αβ 3 µK ) σ2 2 ≤ − M δ + σ 2 K 2 (RS0 − 1) + ϵ + β K + 2µ + r + K 2 µr δ (µ + γ ) 2 ≤−2+1 = − 1, 2

(4.15)

which follows from (4.9) and (4.11). So LV ≤ −1 for any (S , I , R) ∈ D2 . Case 3. For any (S , I , R) ∈ D3 , it follows from (4.10) that LV ≤ −

γI

+

M(β 2 r(µ + γ )2 + αβ 3 µK )

SI + β K + 2µ + r +

σ2

K2

µr δ (µ + γ ) 2 2 2 2 3 γI MK (β r(µ + γ ) + αβ µK ) σ2 2 ≤− + + β K + 2µ + r + K R µr δ (µ + γ ) 2 γϵ MK 2 (β 2 r(µ + γ )2 + αβ 3 µK ) σ2 2 ≤− 2 + + β K + 2µ + r + K ϵ µr δ (µ + γ ) 2 γ MK 2 (β 2 r(µ + γ )2 + αβ 3 µK ) σ2 2 =− + + β K + 2µ + r + K ϵ µ r δ (µ + γ ) 2 ≤ − 1, R

(4.16)

Please cite this article as: Q. Liu and D. Jiang, Threshold behavior in a stochastic SIR epidemic model with Logistic birth, Physica A (2019) 123488, https://doi.org/10.1016/j.physa.2019.123488.

8

Q. Liu and D. Jiang / Physica A xxx (xxxx) xxx

which follows from (4.12). Thus LV ≤ −1 for any (S , I , R) ∈ D3 . Case 4. For any (S , I , R) ∈ D4 , in view of (4.10), we have M(β 2 r(µ + γ )2 + αβ 3 µK ) σ2 2 αI + SI + β K + 2µ + r + K K −N µr δ (µ + γ ) 2 αI MK 2 (β 2 r(µ + γ )2 + αβ 3 µK ) σ2 2 ≤− + + β K + 2µ + r + K K −N µ r δ (µ + γ ) 2 αϵ MK 2 (β 2 r(µ + γ )2 + αβ 3 µK ) σ2 2 ≤− 2 + + β K + 2µ + r + K ϵ µr δ (µ + γ ) 2 MK 2 (β 2 r(µ + γ )2 + αβ 3 µK ) σ2 2 α + β K + 2µ + r + K =− + ϵ µr δ (µ + γ ) 2 ≤ − 1,

LV ≤ −

(4.17)

which follows from (4.13). Consequently, we have LV ≤ −1 for any (S , I , R) ∈ D4 . Therefore, from (4.14)–(4.17), we can obtain that for a sufficiently small ϵ , LV ≤ −1 for all (S , I , R) ∈ Γ ∗ \ Dϵ . By Lemma 2.1, we obtain that there is a positive solution (S(t), I(t), R(t)) to system (1.4) which is a stationary Markov process. This completes the proof. 5. Extinction of the disease In this section, we will establish sufficient conditions for extinction of the disease. Theorem 5.1. Let (S(t), I(t), R(t)) be a solution to system (1.4) with any initial value (S(0), I(0), R(0)) ∈ Γ ∗ . If

σ 2 > max

{

β2 β , 2δ K

} (5.1)

or RS0 =

β βK < 1 and σ 2 < , 1 2 2 K δ + 2σ K

(5.2)

then the disease will go to extinction exponentially a.s., i.e., lim I(t) = 0 a.s.

t →∞

Proof. Applying Itô’s formula to ln I, we have

( d ln I =

βS − δ −

σ2 2

)

S 2 dt + σ SdB(t).

(5.3)

Integrating (5.3) from 0 to t and then dividing by t on both sides, we obtain ln I(t) − ln I(0) t

= :=

1 t 1 t

∫ t( 0 ∫ t( 0

β S(τ ) − δ − β S(τ ) − δ −

σ2 2

σ2 2

∫t

)

S (τ ) d τ + 2

)

S 2 (τ ) d τ +

1

t

∫

t

σ S(τ )dB(τ ) 0

G(t) t

(5.4)

,

∫t

where G(t) := 0 σ S(τ )dB(τ ) which is a continuous local martingale whose quadratic variation is ⟨G, G⟩t = 0 σ 2 S 2 (τ )dτ ≤ σ 2 K 2 t. Applying the strong law of large numbers for local martingales [28], we get limt →∞ G(t) = 0 a.s. Moreover, under t the condition (5.1), we obtain

σ2 2 φ (τ ) :=β S(τ ) − δ − S (τ ) 2 ( )2 σ2 β β2 =− S(τ ) − 2 + −δ 2 σ 2σ 2 β2 ≤ 2 − δ. 2σ Please cite this article as: Q. Liu and D. Jiang, Threshold behavior in a stochastic SIR epidemic model with Logistic birth, Physica A (2019) 123488, https://doi.org/10.1016/j.physa.2019.123488.

Q. Liu and D. Jiang / Physica A xxx (xxxx) xxx

9

Substituting the inequality above into (5.4) and taking the limit on both sides, we have lim

ln I(t) t

t →∞

≤

β2 − δ < 0 a.s. 2σ 2

Now we consider the case σ 2 <

φ (τ ) :=β S(τ ) − δ −

σ

2

2

β K

(5.5)

. In this case, we get

S 2 (τ )

σ2 2 ≤β K − K −δ 2 ( ) 2 σ 2 = δ+ K (RS0 − 1). 2

In view of (5.2), we have lim

ln I(t) t

t →∞

(

≤ δ+

σ2 2

K

2

)

(RS0 − 1) < 0 a.s.

(5.6)

From (5.5) and (5.6) it follows that lim I(t) = 0 a.s.

t →∞

This completes the proof. Remark 5.1. From Theorem 5.1, we can find that if RS0 < 1 and the white noise is not large, then the disease will go to extinction exponentially a.s. In addition, note that the expressions of RS0 and R0 , we can find that RS0 < R0 and if σ = 0, then RS0 = R0 . In other words, the conditions for the disease to die out in system (1.4) are much weaker than in the corresponding deterministic system (1.1). 6. Examples and numerical simulations In this section, we will introduce some numerical simulations to illustrate our theoretical results. We numerically simulate the solution to system (1.4) with initial value (S(0), I(0), R(0)) = (0.8, 0.6, 0.4). For the numerical simulations, we use Milstein’s Higher Order Method mentioned in [32] to obtain the discretization equations of system (1.4)

] ) [( ⎧ √ σ2 2 2 Sk + Ik + Rk ⎪ ⎪ (Sk + Ik + Rk ) − β Sk Ik − µSk ∆t − σ Sk Ik ∆t νk + Sk Ik (νk − 1)∆t , Sk+1 = Sk + b−r ⎪ ⎨ K 2 2 √ σ 2 ⎪ S Ik (νk2 − 1)∆t , Ik+1 = Ik + [β Sk Ik − (α + γ + µ)Ik ]∆t + σ Sk Ik ∆t νk + ⎪ ⎪ 2 k ⎩ Rk+1 = Rk + [γ Ik − µRk ]∆t , where the time increment ∆t > 0, σ 2 > 0 denotes the intensity of the white noise, νk are the independent Gaussian random variables which follow the distribution N(0, 1) for k = 1, 2, . . . , n. Example 6.1. In order to obtain the existence of a stationary distribution numerically, in Fig. 1, we choose parameters b = 0.6, r = 0.4, K = 0.8, β = 0.8, µ = 0.2, γ = 0.1, α = 0.2 and σ 2 = 0.4. Direct calculation leads to . βK RS0 = 1 2 2 = 1.0191 > 1. In other words, the condition in Theorem 4.1 is satisfied. From Theorem 4.1 it follows that δ+ 2 σ K

there is a stationary distribution of system (1.4) which implies that all the individuals are persistent a.s. Fig. 1 confirms this. Example 6.2. In order to examine the extinction of the disease, in Fig. 2, we choose parameters b = 0.6, r = 0.4, K = 0.8, β = 0.8, µ = 0.2, γ = 0.1 and α = 0.2. Case i. In this case, we choose the intensity of environmental white noise σ 2 = 1.2. Direct calculation leads to

} { 2 β β σ 2 = 1.2 > 1 = max , . 2δ K In view of the condition (i) in Theorem 5.1, we can obtain that the disease goes to extinction exponentially with probability one. Fig. 2(a) illustrates this. Case ii. In this case, we choose the intensity of environmental white noise σ 2 = 0.5. By a simple calculation, we have RS0 =

βK β . = 0.9697 < 1 and σ 2 = 0.5 < 1 = . 1 2 2 K δ + 2σ K

By the condition (ii) in Theorem 5.1, we can derive that the disease dies out with probability one. Fig. 2(b) shows this. Please cite this article as: Q. Liu and D. Jiang, Threshold behavior in a stochastic SIR epidemic model with Logistic birth, Physica A (2019) 123488, https://doi.org/10.1016/j.physa.2019.123488.

10

Q. Liu and D. Jiang / Physica A xxx (xxxx) xxx

Fig. 1. The left column shows the paths of S(t), I(t) and R(t) of system (1.4) with initial value (S(0), I(0), R(0)) = (0.8, 0.6, 0.4) under the noise intensity σ 2 = 0.4. The blue lines represent the solution to system (1.4) and the red lines represent the solution to the corresponding undisturbed system (1.1). The right column displays the histogram of the probability density functions of S, I, R populations with values of b = 0.6, r = 0.4, K = 0.8, β = 0.8, µ = 0.2, γ = 0.1, α = 0.2 and σ 2 = 0.4.

Comparing Figs. 1 with 2(a), we can find that the smaller white noise can make the disease persist while the larger white noise can make the disease extinct. This assertion shows that environmental white noise has important effects on controlling the spread of the disease. 7. Conclusion and future directions In this paper, we analyzed a stochastic SIR epidemic model with Logistic birth. By using the stochastic Lyapunov function method, we show that the stochastic basic reproduction number RS0 can be used to determine the threshold dynamics of system (1.4). If RS0 > 1, we establish sufficient conditions for the existence of a stationary distribution of the positive solutions to model (1.4) (i.e., Theorem 4.1). Whereas if RS0 < 1, under extra conditions, we obtain sufficient conditions for extinction of the disease (i.e., Theorem 5.1). This means that environmental fluctuations play positive effects in the control of infectious diseases. Some interesting topics deserve further investigation. On the one hand, one may propose some more realistic but complex models, such as considering the effects of impulsive perturbations on system (1.4). On the other hand, our model Please cite this article as: Q. Liu and D. Jiang, Threshold behavior in a stochastic SIR epidemic model with Logistic birth, Physica A (2019) 123488, https://doi.org/10.1016/j.physa.2019.123488.

Q. Liu and D. Jiang / Physica A xxx (xxxx) xxx

11

Fig. 2(a). The paths of S(t), I(t) and R(t) of system (1.4) with initial value (S(0), I(0), R(0)) = (0.8, 0.6, 0.4) under the noise intensity σ 2 = 1.2.

Fig. 2(b). The paths of S(t), I(t) and R(t) of system (1.4) with initial value (S(0), I(0), R(0)) = (0.8, 0.6, 0.4) under the noise intensity σ 2 = 0.5.

Please cite this article as: Q. Liu and D. Jiang, Threshold behavior in a stochastic SIR epidemic model with Logistic birth, Physica A (2019) 123488, https://doi.org/10.1016/j.physa.2019.123488.

12

Q. Liu and D. Jiang / Physica A xxx (xxxx) xxx

is autonomous, it is interesting to introduce the colored noise, such as continuous-time Markov chain, into model (1.4). These problems will be the subject of our future work. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments This work is supported by the National Natural Science Foundation of China (No. 11871473) and Shandong Provincial Natural Science Foundation, China (No. ZR2019MA010). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32]

W.O. Kermack, A.G. McKendrick, Contributions to the mathematical theory of epidemics-I, Proc. R. Soc. A 115 (1927) 701–721. Z. Ma, Y. Zhou, J. Wu, Modeling and Dynamics of Infectious Diseases, Higher Education Press, 2009. J. Zhang, J.Q. Li, Z.E. Ma, Global analysis of SIR epidemic models with population size dependent contact rate, J. Eng. Math. 21 (2004) 259–267. Y. Cai, Y. Kang, M. Banerjee, W. Wang, A stochastic SIRS epidemic model with infectious force under intervention strategies, J. Differential Equations 259 (2015) 7463–7502. Y. Cai, Y. Kang, M. Banerjee, W. Wang, A stochastic epidemic model incorporating media coverage, Commun. Math. Sci. 14 (2016) 1–18. Y. Cai, Y. Kang, W. Wang, A stochastic SIRS epidemic model with nonlinear incidence rate, Appl. Math. Comput. 305 (2017) 221–240. Y. Cai, K. Wang, W. Wang, Global transmission dynamics of a Zika virus model, Appl. Math. Lett. 92 (2019) 190–195. Y. Cai, J. Jiao, Z. Gui, Y. Liu, W. Wang, Environmental variability in a stochastic epidemic model, Appl. Math. Comput. 329 (2018) 210–226. W. Guo, Y. Cai, Q. Zhang, W. Wang, Stochastic persistence and stationary distribution in an SIS epidemic model with media coverage, Physica A 492 (2018) 2220–2236. B. Yang, Y. Cai, K. Wang, W. Wang, Global threshold dynamics of a stochastic epidemic model incorporating media coverage, Adv. Difference Equ. 2018 (2018) 462. L.J.S. Allen, A.M. Burgin, Comparison of deterministic and stochastic SIS and SIR models in discrete time, Math. Biosci. 163 (2000) 1–33. I. Nåsell, On the time to extinction in recurrent epidemics, J. R. Stat. Soc. Ser. B Stat. Methodol. 61 (1999) 309–330. I. Nåsell, On the quasi-stationary distribution of the stochastic logistic epidemic, Math. Biosci. 156 (1999) 21–40. A.G. McKendrick, Applications of mathematics to medical problems, Proc. Edinburgh Math. Soc. 44 (1926) 98–130. R.M. May, Stability and Complexity in Model Ecosystems, Princeton University, 1973. Y. Lin, D. Jiang, Long-time behaviour of a perturbed SIR model by white noise, Discrete Contin. Dyn. Syst. Ser. B 18 (2013) 1873–1887. E. Tornatore, S.M. Buccellato, P. Vetro, Stability of a stochastic SIR system, Physica A 354 (2005) 111–126. A. Gray, D. Greenhalgh, L. Hu, X. Mao, J. Pan, A stochastic differential equation SIS epidemic model, SIAM J. Appl. Math. 71 (2011) 876–902. Q. Liu, D. Jiang, Stationary distribution and extinction of a stochastic SIR model with nonlinear perturbation, Appl. Math. Lett. 73 (2017) 8–15. Y. Zhang, K. Fan, S. Gao, S. Chen, A remark on stationary distribution of a stochastic SIR epidemic model with double saturated rates, Appl. Math. Lett. 76 (2018) 46–52. C. Xu, S. Yuan, An analogue of break-even concentration in a simple stochastic chemostat model, Appl. Math. Lett. 48 (2015) 62–68. Q. Liu, D. Jiang, Stationary distribution of a stochastic SIS epidemic model with double diseases and the Beddington–DeAngelis incidence, Chaos 27 (2017) 083126. T. Ma, X. Meng, Z. Chang, Dynamics and optimal harvesting control for a stochastic one-predator-two-prey time delay system with jumps, Complexity 2019 (2019) 5342031, 19. F. Zhu, X. Meng, T. Zhang, Optimal harvesting of a competitive n-species stochastic model with delayed diffusions, Math. Biosci. Eng. 16 (2019) 1554–1574. G. Liu, Z. Chang, X. Meng, Asymptotic analysis of impulsive dispersal predator–prey systems with Markov switching on finite-state space, J. Funct. Space 2019 (2019) 8057153, 18. C. Xu, S. Yuan, T. Zhang, Sensitivity analysis and feedback control of noise-induced extinction for competition chemostat model with mutualism, Physica A 505 (2018) 891–902. X. Mao, G. Marion, E. Renshaw, Environmental Brownian noise suppresses explosions in population dynamics, Stochastic Process. Appl. 97 (2002) 95–110. X. Mao, Stochastic Differential Equations and Applications, Horwood Publishing, Chichester, 1997. R. Khasminskii, Stochastic Stability of Differential Equations, Sijthoff and Noordhoff, Alphen aan den Rijn, The Netherlands, 1980. D. Xu, Y. Huang, Z. Yang, Existence theorems for periodic Markov process and stochastic functional differential equations, Discrete Contin. Dyn. Syst. 24 (2009) 1005–1023. N. Dalal, D. Greenhalgh, X. Mao, A stochastic model for internal HIV dynamics, J. Math. Anal. Appl. 341 (2008) 1084–1101. D.J. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Rev. 43 (2001) 525–546.

Please cite this article as: Q. Liu and D. Jiang, Threshold behavior in a stochastic SIR epidemic model with Logistic birth, Physica A (2019) 123488, https://doi.org/10.1016/j.physa.2019.123488.