- Email: [email protected]
Survival and stationary distribution of a SIR epidemic model with stochastic perturbations
Applied Mathematics and Computation 244 (2014) 118–131
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Survival and stationary distribution of a SIR epidemic model with stochastic perturbations Yanli Zhou a,b,⇑, Weiguo Zhang a, Sanling Yuan a a b
University of Shanghai for Science and Technology, Shanghai 200093, China Shanghai Medical Instrumentation College, Shanghai 200093, China
a r t i c l e
i n f o
Keywords: Stochastic SIR epidemic model Stationary distribution Itô formula Extinction
a b s t r a c t In this paper, the dynamics of a SIR epidemic model is investigated. First, we show that the system admits a unique positive global solution starting from the positive initial value. Then, when R0 > 1, we show that the stochastic model has a stationary distribution under certain parametric restrictions. In particular, we show that random effects may lead the disease to extinction in scenarios where the deterministic model predicts persistence. When R0 6 1, a result on fluctuation of the solution around the disease-free equilibrium of deterministic system is established under suitable conditions. Finally, numerical simulations are carried out to illustrate the theoretical results. Ó 2014 Elsevier Inc. All rights reserved.
1. Introduction In 1927, Kermack and McKendrick establish the classical deterministic SIR (susceptible-infected-removed) model [1]. Since then, many people have studied the SIR disease model. Regarding research on the deterministic SIR model and its generalizations, the reader can refer to [2–6]. If SðtÞ; IðtÞ; RðtÞ denote the number of the individuals susceptible to the disease, of infected members and members who have been removed from the possibility of infection through full immunity, respectively, then the differential equations which describe the spread of the disease are:
8 _ > < SðtÞ ¼ K bSðtÞIðtÞ lSðtÞ; _ ¼ bSðtÞIðtÞ ðc þ l þ eÞIðtÞ; IðtÞ > :_ RðtÞ ¼ cIðtÞ lRðtÞ: The parameters in the model are summarized in the following list:
K: the influx of individuals into the susceptible; b: transmission coefficient between compartments S and I; l: natural death rate of S; I; R compartments; e: the disease induced mortality rate; c: the rate of recovery from infection; K; b; l; e; c are all positive. ⇑ Corresponding author at: University of Shanghai for Science and Technology, Shanghai 200093, China. E-mail address: [email protected] (Y. Zhou). http://dx.doi.org/10.1016/j.amc.2014.06.100 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.
ð1:1Þ
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
119
In system (1.1), R0 ¼ lðlbþKeþcÞ is the threshold of the system for an epidemic to persist or disappear. If R0 1, system (1.1) always has a disease-free equilibrium E0 ¼ ðKl ; 0; 0Þ and moreover, it is globally stable in int C, where
C ¼ fðS; I; RÞ : S > 0; I P 0; R P 0; S þ I þ R 6 Klg, this means that the disease dies out. If R0 > 1, then E0 is unstable and there exists a globally asymptotically stable endemic equilibrium E ¼ lþbeþc ; lþKeþc lb ; lðlKþceþcÞ bc , this means that the disease
remains. These above conclusions of system (1.1) can be obtained from [6]. Those important and useful works on deterministic models provide a great insight into the effect of the epidemic model, but in the real world, the models of population dynamics of diseases are inevitably affected by random fluctuations. Introduction of stochasticity into population dynamics models of diseases can bring to light new insights. May [7] pointed out that due to continuous fluctuation in the environment, the birth rates, death rates, transmission coefficient and all other parameters involved with the model exhibit random fluctuations to a greater or lesser extent. As a result the equilibrium distribution never attains a steady value, but fluctuates randomly around some average value [7]. Therefore stochastic differential equation models play a significant role in various branches of applied sciences including the population dynamics of diseases. As related work, the reader can refer to [8–12]. In [10], Jiang et al. assumed that stochastic perturbations are of a white noise type which are directly proportional _ _ _ to SðtÞ; IðtÞ; RðtÞ, influenced on the SðtÞ; IðtÞ; RðtÞ in system (1.1). Accordingly, system (1.1) was modified to the following form:
8 dSðtÞ ¼ ðK bSðtÞIðtÞ lSðtÞÞdt þ r1 SðtÞdB1 ðtÞ; > > < dIðtÞ ¼ ðbSðtÞIðtÞ ðc þ l þ eÞIðtÞÞdt þ r2 IðtÞdB2 ðtÞ; > > : dRðtÞ ¼ ðcIðtÞ lRðtÞÞdt þ r3 RðtÞdB3 ðtÞ;
ð1:2Þ
where Bi ðtÞ ði ¼ 1; 2; 3Þ is a standard Brownian motion, ri ði ¼ 1; 2; 3Þ is the intensity of environmental white noise. In [10], the long time behavior of the stochastic system was studied. In [13], Lin and Jiang assumed that the disease transmission coefficient b is subject to the environmental white noise, that _ is b ! b þ rBðtÞ. Then system (1.1) becomes Itô SDE
8 dSðtÞ ¼ ðK bSðtÞIðtÞ lSðtÞÞdt rSðtÞIðtÞdBðtÞ; > > < dIðtÞ ¼ ðbSðtÞIðtÞ ðc þ l þ eÞIðtÞÞdt þ rSðtÞIðtÞdBðtÞ; > > : dRðtÞ ¼ ðcIðtÞ lRðtÞÞdt;
ð1:3Þ
where BðtÞ is a standard Brownian motion, r is the intensity of environmental white noise. For system (1.3), Lin and Jiang presented sufficient conditions for the disease to go extinct exponentially and found the support of the invariant density. In particular, they obtained the existence of a stationary distribution of system (1.3) and its asymptotic stability by use of the Markov semigroup theory. Motivated by the works of Jiang et al. [10] and Lin and Jiang [13], in this paper, we consider the effects of random fluctuations from environmental white noise. The stochastic version based on the deterministic system (1.1) takes the following form:
8 dSðtÞ ¼ ðK bSðtÞIðtÞ lSðtÞÞdt þ r1 SðtÞdB1 ðtÞ r4 SðtÞIðtÞdB4 ðtÞ; > > < dIðtÞ ¼ ðbSðtÞIðtÞ ðc þ l þ eÞIðtÞÞdt þ r2 IðtÞdB2 ðtÞ þ r4 SðtÞIðtÞdB4 ðtÞ; > > : dRðtÞ ¼ ðcIðtÞ lRðtÞÞdt þ r3 RðtÞdB3 ðtÞ;
ð1:4Þ
where Bi ðtÞ ði ¼ 1; 2; 3; 4Þ is a standard Brownian motion with Bi ð0Þ ¼ 0 ði ¼ 1; 2; 3; 4Þ and ri ði ¼ 1; 2; 3; 4Þ is the intensity of environmental white noise. We are not aware of any literatures analytically concerning the existence of a stationary distribution of system (1.4). The main aim of this paper is to study the existence of a stationary distribution of system (1.4). The paper is organized as follows. In Section 2, the global existence and positivity of the solution of system (1.4) are studied. In Section 3, by choosing appropriate Lyapunov functions, we show that there is a stationary distribution for system (1.4) under certain parametric restrictions. In Section 4, sufficient conditions are derived for the extinction of disease. In Section 5, how the solution spirals around the disease-free equilibrium of deterministic system (1.1) is established. Finally, numerical simulations are carried out to illustrate the theoretical results. Throughout of this paper, unless otherwise specified, let ðX; F; fFt gtP0 ; PÞ be a complete probability space with a filtration satisfying the usual conditions (i.e., it is right continuous and F0 contains all P-null sets). Let Bi ðtÞ ði ¼ 1; 2; 3; 4Þ denote the independent standard Brownian motions defined on this probability space. We also denote
R3þ ¼ fðx; y; zÞ 2 R3 : x > 0; y > 0; z > 0g:
120
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
2. Existence and uniqueness of the positive solution Since the SðtÞ; IðtÞ; RðtÞ in Eq. (1.4) are the sizes of the susceptible individuals, the infected individuals and the removed individuals at time t, respectively, they should be non-negative. It is obvious that we are only interested in positive solutions when we study the epidemic systems. In order to get a unique global solution (i.e., no explosion in a finite time) for any given initial value, the coefficients of the stochastic differential equation are generally required to satisfy the linear growth condition and local Lipschitz condition [14]. However, the coefficients of Eq. (1.4) do not satisfy the linear growth condition, though they are locally Lipschitz continuous, which may include the solution of Eq. (1.4) to explode at a finite time. In what follows, we shall prove that the solution of system (1.4) is positive and global under given conditions. Theorem 2.1. For any given initial condition ðSð0Þ; Ið0Þ; Rð0ÞÞ 2 R3þ , there is a unique positive solution ðSðtÞ; IðtÞ; RðtÞÞ to system (1.4) on t P 0, and the solution remains in R3þ with probability one, namely, ðSðtÞ; IðtÞ; RðtÞÞ 2 R3þ for all t P 0 almost surely. Proof. Since the coefficients of the equation are locally Lipschitz continuous, for any initial value ðSð0Þ; Ið0Þ; Rð0ÞÞ 2 R3þ , there is a unique local solution ðSðtÞ; IðtÞ; RðtÞÞ on t 2 ½0; se Þ, where se is the explosion time (i.e., the moment at which the solution tends to infinity). Next, we show that se ¼ þ1 almost surely. Let m0 > 0 be sufficiently large so that Sð0Þ; Ið0Þ and Rð0Þ lie h i within the interval m10 ; m0 . For each integer m P m0 , define the stopping time
sm ¼ inf t 2 ½0; se Þ : minfSðtÞ; IðtÞ; RðtÞg 6
1 or maxfSðtÞ; IðtÞ; RðtÞg P m ; m
where throughout this paper we set inf / ¼ 1 (as usual / denotes the empty set). According to the definition, sm is increasing as m ! 1. Set s1 ¼ limm!1 sm , whence s1 6 sm a.s. If we can show that s1 ¼ 1 a.s., then se ¼ 1 and ðSðtÞ; IðtÞ; RðtÞÞ 2 R3þ a.s. for all t P 0. If this statement is false, then there exist a pair of constants T > 0 and 2 ð0; 1Þ such that
Pðs1 6 TÞ P : Hence there is an integer m1 P m0 such that
Pfsm 6 Tg P for all m P m1 :
ð2:1Þ
Consider the function VðtÞ by
SðtÞ þ ðIðtÞ 1 log IðtÞÞ þ ðRðtÞ 1 log RðtÞÞ; VðSðtÞ; IðtÞ; RðtÞÞ ¼ SðtÞ c c log c where c is a positive constant to be defined later. The nonnegativity of this function can be seen from the fact that for any fixed constant v > 0; f ðuÞ ¼ u v logðu=v Þ P 0; for all u > 0. With Itô formula, we can get
c c 1 1 1 2 2 dSðtÞ þ dIðtÞ þ dRðtÞ dVðSðtÞ; IðtÞ; RðtÞÞ ¼ 1 ðdSðtÞÞ þ 1 ðdIðtÞÞ þ 1 SðtÞ IðtÞ RðtÞ 2SðtÞ2 2IðtÞ2 þ
1 2RðtÞ2
2
ðdRðtÞÞ ;
then
dVðSðtÞ; IðtÞ; RðtÞÞ ¼ LVðSðtÞ; IðtÞ; RðtÞÞdt þ ðSðtÞ cÞðr1 dB1 ðtÞ r4 IðtÞdB4 ðtÞÞ þ ðIðtÞ 1Þðr2 dB2 ðtÞ þ r4 SðtÞdB4 ðtÞÞ þ ðRðtÞ 1Þr3 dB3 ðtÞ; where LVðSðtÞ; IðtÞ; RðtÞÞ :
R3þ
LVðSðtÞ; IðtÞ; RðtÞÞ ¼
ð2:2Þ
! Rþ by
c cr2 þ r24 I2 1 ðK bSðtÞIðtÞ lSðtÞÞ þ 1 ðbSðtÞIðtÞ ðl þ c þ eÞIðtÞÞ þ 1 SðtÞ IðtÞ 2 r2 þ r24 S2 1 r2 þ 2 ðcI lRÞ þ 3 þ 1 RðtÞ 2 2 1
¼ K þ cl þ l þ l þ e þ c þ c
cr21 þ r24 I2 ðtÞ cK þ ½cb ðl þ eÞIðtÞ ðl þ bÞSðtÞ lRðtÞ 2 SðtÞ
IðtÞ r22 þ r24 S2 ðtÞ r23 þ : þ RðtÞ 2 2
e Choose c ¼ lþ such that ½cb ðl þ eÞIðtÞ ¼ 0, then b
ð2:3Þ
121
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
LVðtÞ 6 K þ cl þ 2l þ c þ e þ
cr21 þ r22 þ r24 I2 ðtÞ þ r24 S2 ðtÞ þ r23 ¼: K: 2
ð2:4Þ
Therefore,
dVðSðtÞ;IðtÞ;RðtÞÞ 6 KdtðSðtÞcÞðr1 dB1 ðtÞ r4 IðtÞdB4 ðtÞÞþðIðtÞ1Þðr2 dB2 ðtÞþ r4 SðtÞdB4 ðtÞÞþðRðtÞ1Þr3 dB3 ðtÞ: ð2:5Þ Integrating both sides of (2.5) from 0 to
Z sm ^T
dVðSðtÞ; IðtÞ; RðtÞÞ 6
Z sm ^T
0
sm ^ T
Kdt þ
0
Z sm ^T 0
þ r4 SðtÞdB4 ðtÞÞ þ
ðSðtÞ cÞðr1 dB1 ðtÞ r4 IðtÞdB4 ðtÞÞ þ
Z sm ^T
Z sm ^T
ðIðtÞ 1Þðr2 dB2 ðtÞ
0
ðRðtÞ 1Þr3 dB3 ðtÞ:
ð2:6Þ
0
Taking the expectations on both sides of (2.6), we obtain
EVðSðsm ^ TÞ; Iðsm ^ TÞ; Rðsm ^ TÞÞ 6 VðSð0Þ; Ið0Þ; Rð0ÞÞ þ KEðsm ^ TÞ; so
EVðSðsm ^ TÞ; Iðsm ^ TÞ; Rðsm ^ TÞÞ 6 VðSð0Þ; Ið0Þ; Rð0ÞÞ þ KT:
ð2:7Þ
Let Xm ¼ sm ^ T for m P m1 . Then, by (2.1), we known that PðXm Þ P . Note that, for every x 2 Xm , there is at least one of Sðsm ^ TÞ; Iðsm ^ TÞ and Rðsm ^ TÞ that equals either m or m1 . Then
VðSð0Þ; Ið0Þ; Rð0ÞÞ þ KT P E½IXm ðxÞVðSðsm ^ TÞ; Iðsm ^ TÞ; Rðsm ^ TÞÞ ¼ E½IXm ðxÞVðSðsm ; xÞ; Iðsm ; xÞ; Rðsm ; xÞÞ 1 1 þ log m ; P ðm 1 log mÞ ^ m where IXm ðxÞ is the indicator function of Xm ðxÞ. Letting m ! 1, it leads to the contradiction
1 > VðSð0Þ; Ið0Þ; Rð0ÞÞ þ KT ¼ 1; which is a contradiction, therefore we have s1 ¼ 1 a.s. Therefore, it implies SðtÞ; IðtÞ and RðtÞ will not explode in a finite time with probability one. The proof is therefore completed. h
3. Existence of the stationary distribution An important aspect on studying the disease is the long term prevalence of the disease or whether it will ultimately vanish from the population. In this subsection, we show that there is a stationary distribution when the white noise is small. The stationary distribution proofs in this paper are similar to [15]. Before proving the main theorem of this section we first cite a known result from Has’miskii [16] which will be useful to prove the theorem. Let XðtÞ be a regular time-homogeneous Markov process in El (the l dimensional Euclidean space) described by stochastic equation
dXðtÞ ¼ bðXÞdt þ
k X
f r dBr ðtÞ:
ð3:1Þ
r¼1
The diffusion matrix is defined as follows:
AðxÞ ¼ ðaij ðxÞÞ;
aij ðxÞ ¼
k1 X
fri ðxÞfrj ðxÞ:
r¼1
There exists a bounded domain U El with regular boundary C, having the following properties: (B.1) In the domain U and some neighborhood thereof, the smallest eigenvalue of the diffusion matrix AðxÞ is bounded away from zero. (B.2) If x 2 El n U, the mean time s at which a path emerging from x reaches the set U is finite, and supx2K Ex s < 1 for every compact subset K El .
Lemma 3.1. If above assumptions (B.1) and (B.2) hold, then the Markov process XðtÞ has a stationary distribution lðÞ. Let f ðÞ be a integrable function with respect to the measure lðÞ. Then
122
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
( Px lim
T*1
Z
T
f ðXðtÞÞdt ¼
Z
0
) f ðxÞlðdxÞ
¼1
El
for all x 2 El . Remark 3.1. Lemma 3.1 asserts the existence of a stationary distribution with suitable density function. Its proof can be found in [16]. To validate (B.1), it suffices to prove V is uniformly elliptical in U, where Vu ¼ bðxÞ:ux þ ½trðAðxÞuxx =2, i.e. there is a positive number M such that k X ai;j ni nj P Mjnj2 ;
x 2 U; n 2 Rk
i;j¼1
(see [17] and Rayleigh’s principle in [18]). To verify (B.2), it is sufficient to show that there exists some neighborhood U and a non-negative C 2 -function such that LV is negative at every point x 2 El n U like in [19]. Theorem 3.1. Consider the stochastic system (1.4) with initial condition in R3þ . Assume that R0 > 1 and 0 < d < minðm1 S2 ; m2 I2 ; m3 R2 Þ, then there exists a stationary distribution lðÞ for system (1.4). Here ðS ; I ; R Þ is the unique endemic equilibrium of system (1.1),
m1 ¼ l r21
r24 ð2l þ c þ eÞI b
m2 ¼ l þ c r22
;
2
c
2l
;
l
r23 ; 2l þ c þ e 2 r24 I ð2l þ e þ cÞ 2 d ¼ r21 S2 þ I2 þ I r2 þ S þ r23 R2 : b 2b
m3 ¼
ð3:2Þ
2
Especially, we have
1 lim sup E t t!1
Z
t
2
2
2
½m1 ðSðuÞ S Þ þ m2 ðIðuÞ I Þ þ m3 ðRðuÞ R Þ du 6 d:
ð3:3Þ
0
Proof. Since R0 > 1, then there exists a positive equilibrium of system (1.1), and
K ¼ bS I þ lS ; bS I ¼ ðc þ l þ eÞI ; cI ¼ lR :
ð3:4Þ
Let us consider the function
VðSðtÞ; IðtÞ; RðtÞÞ ¼ V 1 ðtÞ þ
2l þ c þ e V 2 ðtÞ þ V 3 ðtÞ; b
ð3:5Þ
where
V 1 ðtÞ ¼
1 2 ðSðtÞ þ IðtÞ S I Þ ; 2
V 2 ðtÞ ¼ IðtÞ I I ln
IðtÞ ; I
2
V 3 ðtÞ ¼
ðRðtÞ R Þ : 2
ð3:6Þ
Using Itô formula, we have
dVðtÞ ¼ dV 1 ðtÞ þ
2l þ c þ e dV 2 ðtÞ þ dV 3 ðtÞ; b
in detail,
1 2 dV 1 ðtÞ ¼ ðSðtÞ þ IðtÞ S I ÞðdSðtÞ þ dIðtÞÞ þ ðdSðtÞ þ dIðtÞÞ 2 ¼ LV 1 ðtÞdt þ ðSðtÞ þ IðtÞ S I Þðr1 SðtÞdB1 ðtÞ þ r2 IðtÞdB2 ðtÞÞ; I I 2 dV 2 ðtÞ ¼ 1 dIðtÞ þ 2 ðdIðtÞÞ ¼ LV 2 dt þ ðIðtÞ I Þr2 dB2 ðtÞ þ ðIðtÞ I Þr4 SðtÞdB4 ðtÞ; IðtÞ 2I ðtÞ 2
dV 3 ðtÞ ¼ ðRðtÞ R ÞdRðtÞ þ ðdRðtÞÞ ¼ LV 3 ðtÞdt þ ðRðtÞ R ÞRðtÞr3 dB3 ðtÞ;
ð3:7Þ
123
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
where
LV 1 ðtÞ ¼ ðSðtÞ þ IðtÞ S I Þ½K lSðtÞ ðl þ c þ eÞIðtÞ þ
r21 S2 ðtÞ þ r22 I2 ðtÞ 2
LV 2 ðtÞ ¼ ðIðtÞ I Þ½bSðtÞ ðl þ c þ eÞ þ
;
I ðr2 þ r24 S2 ðtÞÞ; 2 2
ð3:9Þ
1 LV 3 ðtÞ ¼ ðRðtÞ R ÞðcIðtÞ lRðtÞÞ þ r23 R2 ðtÞ: 2 2
ð3:8Þ
ð3:10Þ
2
Using (3.4) and the inequality ða þ bÞ 2a2 þ 2b , we can have from (3.8) that
1 2 LV 1 ðtÞ ¼ ðSðtÞ S þ IðtÞ I Þ½lS þ ðl þ c þ eÞI lSðtÞ ðl þ c þ eÞIðtÞ þ ½r21 ðSðtÞ S þ S Þ 2 2 2 þ r22 ðIðtÞ I þ I Þ 6 ðSðtÞ S þ IðtÞ I Þ½lðSðtÞ S Þ ðl þ c þ eÞðIðtÞ I Þ þ r21 ðS S Þ 2
2
þ r21 S2 þ r22 ðI I Þ þ r22 I2 ¼ ðl r21 ÞðSðtÞ S Þ ð2l þ c þ eÞðSðtÞ S ÞðIðtÞ I Þ 2
ðl þ c þ e r22 ÞðIðtÞ I Þ þ r21 S2 þ r22 I2 :
ð3:11Þ
Similarly, we can have from (3.9) and (3.10) that
LV 2 ðtÞ ¼ ðIðtÞ I ÞðbSðtÞ bS Þ þ 6 bðSðtÞ S ÞðIðtÞ I Þ þ
I 2 r24 I 2 r þ ðSðtÞ S þ S Þ 2 2 2
I 2 r þ r24 I ðSðtÞ S Þ2 þ r24 I S2 2 2
ð3:12Þ
and
LV 3 ðtÞ ¼ ðRðtÞ R ÞðcIðtÞ lRðtÞ cI þ lR Þ þ
1 2 r ðRðtÞ R þ R Þ2 2 3
2
2
6 cðRðtÞ R ÞðIðtÞ I Þ lðRðtÞ R Þ þ r23 ðRðtÞ R Þ þ r23 R pffiffiffiffi c 2 2 ¼ pffiffiffiffi ðIðtÞ I Þ lðRðtÞ R Þ lðRðtÞ R Þ þ r23 ðRðtÞ R Þ þ r23 R2
l
2
c
6
2l
2
ðIðtÞ I Þ
l 2
2
2
ðRðtÞ R Þ þ r23 ðRðtÞ R Þ þ r23 R2 :
ð3:13Þ
It follows from (3.7) that
LVðtÞ ¼ LV 1 ðtÞ þ
2l þ c þ e LV 2 ðtÞ þ LV 3 ðtÞ: b
ð3:14Þ
Substituting (3.11)–(3.13) into (3.14), we obtain that
l r2 ð2l þ c þ eÞI c2 2 2 2 ðSðtÞ S Þ l þ c r22 ðIðtÞ I Þ LVðtÞ 6 l r21 4 r23 ðRðtÞ R Þ b 2l 2 2l þ c þ e 2 r24 I ð2l þ e þ cÞ 2 2 2 þ r21 S2 þ I2 þ I r2 þ S þ r23 R2 , m1 ðSðtÞ S Þ m2 ðIðtÞ I Þ b 2b 2
m3 ðRðtÞ R Þ þ d;
ð3:15Þ
where m1 ; m2 ; m3 and d are defined respectively in (3.2). Therefore, 2
2
2
dVðtÞ 6 m1 ðSðtÞ S Þ m2 ðIðtÞ I Þ m3 ðRðtÞ R Þ þ d þ ðSðtÞ þ IðtÞ S I Þðr1 SðtÞdB1 ðtÞ
þ r2 IðtÞdB2 ðtÞÞ þ ðIðtÞ I Þr2 dB2 ðtÞ þ ðIðtÞ I Þr4 SðtÞdB4 ðtÞ þ ðRðtÞ R ÞRðtÞr3 dB3 ðtÞ:
ð3:16Þ
Integrating both sides of (3.16) from 0 to t, yields
VðtÞ Vð0Þ 6
Z
t
2
2
2
½m1 ðSðuÞ S Þ m2 ðIðuÞ I Þ m3 ðRðuÞ R Þ du þ dt þ MðtÞ;
ð3:17Þ
0
where MðtÞ is a martingale defined by
MðtÞ ¼
Z 0
þ
t
ðSðuÞ þ IðuÞ S I Þðr1 SðuÞdB1 ðuÞ þ r2 IðuÞdB2 ðuÞÞ þ
Z 0
Z 0
t
ðRðuÞ R ÞRðuÞr3 dB3 ðuÞ:
t
ððIðuÞ I Þr2 dB2 ðuÞ þ ðIðuÞ I Þr4 SðuÞdB4 ðuÞÞ
124
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
Taking expectation on both sides of (3.17), yields
EVðtÞ EVð0Þ 6 E
Z
t
2
2
2
½m1 ðSðuÞ S Þ m2 ðIðuÞ I Þ m3 ðRðuÞ R Þ du þ dt:
0
Hence,
1 lim sup E t t!1
Z
t
2
2
2
½m1 ðSðuÞ S Þ þ m2 ðIðuÞ I Þ þ m3 ðRðuÞ R Þ du 6 d:
0
That is, we obtain (3.3). Noticing that if 0 < d < minðm1 S2 ; m2 I2 ; m3 R2 Þ, then the ellipsoid 2
2
2
m1 ðSðtÞ S Þ m2 ðIðtÞ I Þ m3 ðRðtÞ R Þ þ d ¼ 0 lies entirely in R3þ . We can take U to be any neighborhood of the ellipsoid such that U # El ¼ R3þ , so for x 2 U n El ; LVðtÞ 6 C (C is a positive constant), which implies condition (B.2) is satisfied. On the other hand, we can write system (1.4) as the form of system (3.1):
0
1 0 1 0 1 0 1 K bSðtÞIðtÞ lSðtÞ þ cIðtÞ r1 SðtÞ 0 0 B C B C B C B C B C d@ IðtÞ A ¼ @ bSðtÞIðtÞ ðl þ e þ cÞIðtÞ Adt þ @ 0 AdB1 ðtÞ þ @ r2 IðtÞ AdB2 ðtÞ þ @ 0 AdB3 ðtÞ RðtÞ 0 0 cIðtÞ lRðtÞ r3 RðtÞ 0 1 r4 SðtÞIðtÞ B C þ @ r4 SðtÞIðtÞ AdB4 ðtÞ: 0 SðtÞ
1
0
Here the diffusion matrix is
0
r2 S2 þ r24 S2 I2 r24 S2 I2 B 1 A¼B r22 I2 þ r24 S2 I2 r24 S2 I2 @ 0
0 0
r23 R2
0
1 C C: A
There is a M ¼ minfr21 S2 ; r22 I2 ; r23 R2 g > 0, such that for all ðx1 ; x2 ; x3 Þ 2 U and n 2 R3 , 3 X aij ni nj ¼ r21 S2 þ r24 S2 I2 n21 þ r22 I2 þ r24 S2 I2 n22 2r24 S2 I2 n1 n2 þ r23 R2 n23 i;j¼1
¼ r21 S2 n21 þ r22 I2 n22 þ r23 R2 n23 þ r24 S2 I2 ðn1 n2 Þ2 P r21 S2 n21 þ r22 I2 n22 þ r23 R2 n23 P min
n
o
r21 S2 ; r22 I2 ; r23 R2 jnj2
¼ Mjnj2 ; which shows that condition (B.1) is also satisfied. Therefore, we can conclude that stochastic system (1.4) has a stationary distribution lðÞ. h Next, we can further have the following almost sure persistence of the total population result of system (1.4). Set
XðtÞ ¼ SðtÞ þ IðtÞ þ RðtÞ
ð3:18Þ
Theorem 3.2. If K > 0, then for any given initial ðSð0Þ; Ið0Þ; Rð0ÞÞ 2 R3þ , the solution of model (1.4) obeys
0 < lim inf XðtÞ 6 lim sup XðtÞ < þ1: t!1
ð3:19Þ
t!1
Proof. It follows from (3.18) that
dXðtÞ ¼ ðK lXðtÞ eIðtÞÞdt þ r1 SðtÞdB1 ðtÞ þ r2 IðtÞdB2 ðtÞ þ r3 RðtÞdB3 ðtÞ: 1 Define a function VðXðtÞÞ ¼ XðtÞ and by the Itô formula, we have
d
1 1 1 2 ¼ 2 dXðtÞ þ 3 ðdXðtÞÞ XðtÞ X ðtÞ X ðtÞ ¼
1 2
X ðtÞ
ðK lXðtÞ eIðtÞÞdt þ
þ r3 RðtÞdB3 ðtÞÞ:
1 3
X ðtÞ
ðr21 S2 þ r22 I2 þ r23 R2 Þ
1 2
X ðtÞ
ðr1 SðtÞdB1 ðtÞ þ r2 IðtÞdB2 ðtÞ
ð3:20Þ
125
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
Since
SðtÞ XðtÞ
d
IðtÞ 1; XðtÞ 1 and
"
RðtÞ XðtÞ
1, then
# 1 t K l þ e r21 þ r22 þ r23 et e ¼ 2 þ þ þ 1 et dt 2 ðr1 SðtÞdB1 ðtÞ þ r2 IðtÞdB2 ðtÞ þ r3 RðtÞdB3 ðtÞÞ: XðtÞ XðtÞ X ðtÞ XðtÞ X ðtÞ
It is easy to obtain that
1 1 t ¼ e þ et XðtÞ Xð0Þ
Z
t
eu HðXðuÞÞdu et
"Z
0
0
t
eu X 2 ðuÞ
r1 SðuÞdB1 ðuÞ þ
Z 0
t
eu X 2 ðuÞ
r2 IðuÞdB2 ðuÞ þ
Z
t
0
eu X 2 ðuÞ
#
r3 RðuÞdB3 ðuÞ : ð3:21Þ
Obviously, if K > 0, then we have HðXðtÞÞ 6 M1 for all XðtÞ > 0. In fact, we only need to prove the next two equations
Pðlim inf Xðx; tÞ ¼ 0Þ ¼ 0
ð3:22Þ
Pðlim sup Xðx; tÞ ¼ þ1Þ ¼ 0:
ð3:23Þ
t!1
and t!1
We now begin to prove assertion (3.22). If it is not true, then there is a sufficiently small d 2 ð0; 1Þ such that
PðX1 Þ > d;
ð3:24Þ
where X1 ¼ fxj lim inf Xðx; tÞ ¼ 0g. Define the stopping time
t!1
1 k
sk ¼ inf t P 0j; Xðx; tÞ 6 ; x 2 X1 ; k 2 Z þ : According to the definition,
sk is increasing and sk ! 1 as k ! 1. It then follows from (3.21) that
Z sk 1 1 sk 1 þ M1 : E I X1 ¼ E e þ esk eu HðXÞdu 6 Xðsk Þ Xð0Þ Xð0Þ 0
ð3:25Þ
However,
1 E IX1 P kEðIX1 Þ P kd: Xðsk Þ Letting k ! 1, leads to yield
1 E IX1 P kEðIX1 Þ P kd ! 1; Xðsk Þ
as k ! 1:
But this contradicts (3.25). We therefore must have the desired assertion (3.22), namely, lim inf t!1 XðtÞ > 0 a.s. On the other hand, it follows from (3.20) that Z t Z t Z t Z t XðtÞ ¼ Xð0Þelt þ elt elu ðK eIðuÞÞdu þ elt elu r1 SðuÞdB1 ðuÞ þ elt elu r2 IðuÞdB2 ðuÞ þ elt elu r3 RðuÞdB3 ðuÞ: 0
0
0
0
In the same way as in the proof of (3.22), we find that
PðX2 Þ ¼ 0; where X2 ¼ fxj lim supt!1 Xðx; tÞ ¼ þ1g a.s.. That is, lim supt!1 XðtÞ < þ1. The proof is hence complete.
h
4. Extinction of the disease In the deterministic models, the value of the basic reproductive number R0 determines persistence or extinction of the disease. If R0 1, the disease is eliminated, whereas if R0 > 1, the disease persists in the population. However, we will show in this section that if the noise is sufficiently large, the disease will become extinct for stochastic system (1.4), although it may be persistent for its deterministic version (1.1). The following theorem gives a condition for the extinction of a disease expressed in terms of intensities of noise and system parameters. Theorem 4.1. Consider stochastic system (1.4) with initial condition in R3þ . We have
lim sup t!1
If
b2 2r24
ln IðtÞ b2 r2 < ðc þ l þ eÞ 2 a:s: t 2 2r24
< ð c þ l þ eÞ þ
r22 2
holds, then I (t) will go to zero exponentially with probability one.
126
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
Before proving the main theorem of this section, we first cite a Lemma from Mao [20] which will be useful to prove the theorem. Lemma 4.1. Let MðtÞ : t P 0 be a continuous local martingale with value Mð0Þ ¼ 0. Let hMðtÞi be its quadratic variation. Let d > 1 be a number, and let mn and sn be two sequences of positive numbers. Then, for almost all x 2 X, there exists a random integer n0 ¼ n0 ðxÞ such that, for all n P n0 ,
MðtÞ 6
1 d ln n mn hMðtÞ; MðtÞi þ ; 2 mn
0 6 t 6 sn :
ð4:1Þ
Proof. Define VðIðtÞÞ ¼ ln IðtÞ and by Itô formula, we have
# Z t" Z t r22 r24 S2 ðuÞ du þ r2 B2 ðtÞ þ ln IðtÞ ¼ ln Ið0Þ þ bSðuÞ ðc þ l þ eÞ r4 SðuÞdB4 ðtÞ 2 2 0 0
Z t r2 r2 ¼ ln Ið0Þ þ 4 S2 ðuÞ þ bSðuÞ ðc þ l þ eÞ 2 du þ r2 B2 ðtÞ þ MðtÞ; 2 2 0 where MðtÞ ¼
Rt 0
r4 SðuÞdB4 ðuÞ is a continuous local martingale whose quadratic variation is
hMðtÞ; MðtÞi ¼ r24
1 2 mr 2 4
Z
t
S2 ðuÞdu:
0
mn ¼ m > 0 and sn ¼ n in (4.1). For almost all x 2 X, there is a random integer n0 ðxÞ such that, for n > n0 ,
Choosing d ¼ 2;
MðtÞ 6
ð4:2Þ
Z
t
S2 ðuÞdu þ
0
2
m
ln n for all 0 6 t 6 n:
ð4:3Þ
Substituting (4.3) into (4.2), we obtain that
ln IðtÞ < ln Ið0Þ þ
Z t Z t r2 r2 1 2 4 S2 ðuÞ þ bSðuÞ ðc þ l þ eÞ 2 du þ r2 B2 ðtÞ þ mr24 S2 ðuÞdu þ ln n: 2 2 2 m 0 0
ð4:4Þ
Clearly, we have that
2 1 2b 1 b b2 b2 ð1 mÞr24 S2 ðtÞ 2 SðtÞ ¼ ð1 mÞr24 SðtÞ þ 6 : 2 2 ð1 mÞr24 r4 2ð1 mÞr24 2ð1 mÞr24
ð4:5Þ
It follows from (4.4) and (4.5) that
" ln IðtÞ < ln Ið0Þ þ
# b2 r22 2 ð c þ l þ e Þ t þ r2 B2 ðtÞ þ ln n 2 m 2ð1 mÞr24
ð4:6Þ
Thus, for n 1 6 t 6 n, we obtain
ln IðtÞ ln Ið0Þ b2 r2 r2 B2 ðtÞ 2 ln n < þ ðc þ l þ eÞ 2 þ þ : t t t 2 m n1 2ð1 mÞr24 By the large number theorem for martingales [14], we have
lim
t!1
B2 ðtÞ ¼ 0: t
Hence, let n ! 1, so t ! 1, we obtain
lim sup t!1
Let
ln IðtÞ b2 r2 < ðc þ l þ eÞ 2 : 2 t 2 2ð1 mÞr4
m ! 0, we obtain lim sup t!1
ln IðtÞ b2 r2 < ðc þ l þ eÞ 2 : 2 t 2 2r4
The proof is therefore completed. h Remark 4.1. When R0 > 1, the positive solution converges to the endemic equilibrium of deterministic system (1.1). However, the disease will die2 out exponentially regardless of the magnitude of R0 provided that r2 and r4 are big enough such 2 r that 2br2 < ðc þ l þ eÞ þ 22 . That is to say, large noises can lead to the extinction of disease. 4
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
127
5. Asymptotic behavior around the disease-free equilibrium E0 As mentioned in system (1.1), if R0 1, then system (1.1) always has a globally asymptotically stable disease-free equilibrium E0 , which means the disease will die out with the advancement of time. Noting that E0 is not an equilibrium of stochastic system (1.4), it is natural to ask whether the disease will go to extinction in the population. In this section we mainly use the way of estimating the oscillation around E0 to reflect how the solution of system (1.4) spirals closely around E0 . We have the following theorem. Theorem 5.1. Assume R0 1. Then for any solution ðSðtÞ; IðtÞÞ of system (1.4) with initial value ðSð0Þ; Ið0Þ; Rð0ÞÞ 2 R3þ , if r21 < l; r22 < 2ðl þ eÞ þ c and r23 < 2l c, we have
lim sup t!1
1 E t
# Z t " K 2 r 2 K2 SðuÞ þ I2 ðuÞ þ R2 ðuÞ du 6 1 2 ; l ml 0
ð5:1Þ
n o r2 r2 where m ¼ min l r21 ; l þ e þ 2c 22 ; l 2c 23 . Proof. Define a function V : R3þ ! Rþ by
VðtÞ ¼ V 1 ðtÞ þ
2l þ e þ c V 2 ðtÞ þ V 3 ðtÞ ¼ b
2 SðtÞ Kl þ IðtÞ 2
þ
2l þ e þ c 1 IðtÞ þ R2 ðtÞ: 2 b
ð5:2Þ
Along the trajectories of system (1.4), we have
dVðtÞ ¼ dV 1 ðtÞ þ
2l þ e þ c dV 2 ðtÞ þ dV 3 ðtÞ: b
ð5:3Þ
In detail,
K 1 K 2 dV 1 ðtÞ ¼ SðtÞ þ IðtÞ ðdSðtÞ þ dIðtÞÞ þ ðdSðtÞ þ dIðtÞÞ ¼ LV 1 ðtÞdt þ SðtÞ þ IðtÞ ðr1 SðtÞdB1 ðtÞ þ r2 IðtÞdB2 ðtÞÞ; 2 l l dV 2 ðtÞ ¼ LV 2 ðtÞdt þ r2 IðtÞdB2 ðtÞ þ r4 SðtÞIðtÞdB4 ðtÞ; and
dV 3 ðtÞ ¼ LV 3 ðtÞdt þ r3 R2 ðtÞdB3 ðtÞ; where
K r2 S2 ðtÞ þ r22 I2 ðtÞ LV 1 ðtÞ ¼ SðtÞ þ IðtÞ ðK lSðtÞ ðl þ e þ cÞIðtÞÞ þ 1 ; 2 l
ð5:4Þ
LV 2 ðtÞ ¼ bSðtÞIðtÞ ðl þ c þ eÞIðtÞ;
ð5:5Þ
LV 3 ðtÞ ¼ cRðtÞIðtÞ lR2 ðtÞ þ
r
2 2 3 R ðtÞ
2
ð5:6Þ
:
From (5.3), we have that
2
r2 SðtÞ K þ K þ r2 I2 ðtÞ 1 2 l l K K ðl þ e þ cÞIðtÞ þ LV 1 ðtÞ ¼ SðtÞ þ IðtÞ l SðtÞ l l 2 2 K 2 K r2 K IðtÞ l þ e þ c 2 I2 ðtÞ þ r21 ð2l þ e þ cÞ SðtÞ ; 6 ðl r21 Þ SðtÞ l l 2 l 2
ð5:7Þ
2
where in the last step we have used the inequality ða þ bÞ 2a2 þ 2b . Similarly, we can have from (5.4) and (5.5) that
K K K K IðtÞ ðl þ c þ eÞIðtÞ ¼ b SðtÞ IðtÞ þ b ðl þ c þ eÞ IðtÞ LV 2 ðtÞ ¼ b SðtÞ þ
l
l
l
l
and
LV 3 ðtÞ 6
c 2
Noting that R0 ¼
½R2 ðtÞ þ I2 ðtÞ lR2 ðtÞ þ Kb lðlþeþcÞ
r23 R2 ðtÞ 2
c r2 c ¼ l 3 R2 ðtÞ þ I2 ðtÞ: 2 2 2
1, it follows from above that
K IðtÞ: LV 2 ðtÞ 6 b SðtÞ
l
ð5:8Þ
128
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
It follows from (5.6), (5.7) and (5.2) that
2 K 2 c r2 c r2 K LVðtÞ 6 ðl r21 Þ SðtÞ l þ e þ 2 I2 ðtÞ l 3 R2 ðtÞ þ r21 : l 2 2 2 2 l
ð5:9Þ
Integrating from 0 to t on both sides of (5.3), then taking expectation, and considering inequality (5.9), yields
EVðtÞ Vð0Þ 6 ðl r21 ÞE þ r21
2 K
l
Z t Z t Z t K 2 c r2 c r2 SðuÞ du l þ e þ 2 E I2 ðuÞdu l 3 E R2 ðuÞdu l 2 2 2 2 0 0 0
t:
Hence, we have that
1 lim sup E t t!1 Let m ¼ min
# 2 Z t" K 2 c r2 c r2 K ðl r21 Þ SðuÞ þ l þ e þ 2 I2 ðuÞ þ l 3 R2 ðuÞ du 6 r21 : l 2 2 2 2 l 0
n
2
2
ð5:10Þ
o
l r21 ; l þ e þ 2c r22 ; l 2c r23 , it follows from (5.9) that
1 lim sup E t t!1
# Z t " K 2 r 2 K2 2 2 SðuÞ þ I ðuÞ þ R ðuÞ du 6 1 2 : l ml 0
The proof of Theorem 5.1 is thus completed. h Notice that when r1 ¼ 0; E0 is also the disease-free equilibrium of system (1.4) and (5.9) is reduced to the following form
K 2 c r2 c r2 LVðtÞ 6 l SðtÞ l þ e þ 2 I2 ðtÞ l 3 R2 ðtÞ; l 2 2 2 2 which means that LVðtÞ is negative definite provided that r22 < 2ðl þ eÞ þ c and r23 < 2l c. Therefore, we can have the following result [14]. Corollary 5.1. Assume that R0 1 and r1 ¼ 0. Then the disease-free equilibrium E0 of system (1.4) is stochastically asymptotically stable in the large provided r22 < 2ðl þ eÞ þ c and r23 < 2l c. 6. Discussions and numerical simulations In order to make the readers understand our results better, we will perform some numerical simulations to illustrate our theoretical results. The numerical simulations are given by the following Milstein scheme [21]. Consider the discretization of model (1.4) for t ¼ 0; Mt; 2Mt; ; nMt: Group(a)
S(t)
Relative frequency density
2
100
1
50
0
0
100
200
300
400
0 −1
t I(t)
Group(b)
100
0.5
50 0
100
200
300
400
0 −1
t R(t)
Group(c)
100
1
50 0
100
200
t
0 1 I(t) at time 150
2
Relative frequency density
2
0
2
Relative frequency density
1
0
0 1 S(t) at time 150
300
400
0 −1
0 1 R(t) at time 150
Fig. 1. The solution of stochastic system (1.4) and its histogram with Sð0Þ ¼ 0:62; Ið0Þ ¼ 0:54; Rð0Þ ¼ 0:54; K ¼ 0:4; b ¼ 0:8; r1 ¼ 0:1; r2 ¼ 0:1; r3 ¼ 0:1; r4 ¼ 0:1.
2
l ¼ 0:2; e ¼ 0:1; c ¼ 0:2;
129
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
Group(a)
S(t)
Relative frequency density
2
4
1
2
0
0
100
200
300
0 −1
400
4
1
2 0
100
200
300
0 −1
400
0 1 I(t) at time 150
t R(t)
Group(c)
2
Relative frequency density
2
0
0 1 S(t) at time 150
t I(t)
Group(b)
2
Relative frequency density
4
5
2 0
0
100
200
300
0 −1
400
0 1 R(t) at time 150
t Fig. 2. The solution of stochastic system (1.4) and e ¼ 0:1; c ¼ 0:2; r1 ¼ 0:2; r2 ¼ 0:2; r3 ¼ 0:2; r4 ¼ 0:2.
histogram
with
Sð0Þ ¼ 0:62; Ið0Þ ¼ 0:54; Rð0Þ ¼ 0:54; K ¼ 0:4; b ¼ 0:8;
pffiffiffiffiffiffiffi 1 1 2 2
pffiffiffiffiffiffiffi 1 2 2 pffiffiffiffiffiffi 1 2 2 Ikþ1 ¼ Ik þ ½bSk Ik ðc þ l þ ÞIk Mt þ Ik r2 n2;k Mt þ r2 ðn21;k 1ÞMt þ Sk Ik r4 n4;k 4t þ r4 ðn4;k 1ÞMt ; 2 2
pffiffiffiffiffiffi 1 2 2 Rkþ1 ¼ Rk þ ðcIk lRk ÞMt þ Rk r3 n3;k Mt þ r3 ðn3;k 1ÞMt ; 2
Skþ1 ¼ Sk þ ðK bSk Ik lSk ÞMt þ Sk
its
2
l ¼ 0:2;
pffiffiffiffiffiffi
r1 n1;k Mt þ r21 ðn21;k 1ÞMt þ Sk Ik r4 n4;k 4t þ r24 ðn24;k 1ÞMt ;
ð6:1Þ
where Mt is time increment n1k ; n2k ; n3k and n4k ðk ¼ 1; 2; . . . ; nÞ are Nð0; 1Þ-distributed independent random variables. In the figures, the blue lines and the red lines represent solutions of deterministic system (1.1) and stochastic system (1.4) respectively. In Figs. 1–3, we choose the parameters K ¼ 0:4; b ¼ 0:8; l ¼ 0:2; e ¼ 0:1; c ¼ 0:2 and the initial point Sð0Þ ¼ 0:62; Ið0Þ ¼ 0:54; Rð0Þ ¼ 0:54, for numerical simulations of stochastic system (1.4). We use different values of r1 ; r2 ; r3 and r4 in order to understand their role on the dynamics.
3
1.4
1.2
2.5
1 2
I(t)
S(t)
0.8 1.5
0.6 1 0.4 0.5
0
0.2
0
50
100
t
150
0
0
50
100
150
t
Fig. 3. Deterministic and stochastic trajectories of system for initial condition and parameter values: Sð0Þ ¼ 0:62; Ið0Þ ¼ 0:54; Rð0Þ ¼ 0:54; K ¼ 0:4; b ¼ 0:8; l ¼ 0:2; e ¼ 0:1; c ¼ 0:2; r1 ¼ 0:1; r2 ¼ 0:7; r3 ¼ 0:1; r4 ¼ 0:6.
130
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
Group(a) 1.5 0.2
I(t)
R(t)
0.2
S(t)
1 0.5
0
50
100
0
150
0.1
0
50
t
Group(b) 1.5
100
0.3
0.3
0.2
0.2
0.1
50
50
100
0
150
100
150
100
150
t
I(t)
S(t) 0.5
0
0
t
1
0
0
150
R(t)
0
0.1
0.1
0
50
t
100
0
150
0
50
t
t
Fig. 4. Deterministic and stochastic trajectories of system for initial condition and parameter values: Sð0Þ ¼ 0:7; Ið0Þ ¼ 0:2; Rð0Þ ¼ 0:1; K ¼ 0:2; b ¼ 0:4; l ¼ 0:2; e ¼ 0:1; c ¼ 0:2; Group (a): r1 ¼ 0:1; r2 ¼ 0:1; r3 ¼ 0:1, r4 ¼ 0:1; Group (b): r1 ¼ 0:05; r2 ¼ 0:05; r3 ¼ 0:05; r4 ¼ 0:05.
For a clear understanding of this idea we can look at the stationary distribution for parameter values as mentioned in Section 3 and for two different sets of forcing intensities r1 ¼ 0:2; r2 ¼ 0:2; r3 ¼ 0:2; r4 ¼ 0:2 and r1 ¼ 0:1; r2 ¼ 0:1; r3 ¼ 0:1; r4 ¼ 0:1. Parameter values chosen above are consistent with the conditions required for the existence of stationary distribution (see Theorem 3.1). The distributions of the population obtained through numerical simulation for two different sets of forcing intensities are reported in Figs. 1 and 2. These figures report the results of one simulation run and stationary distribution about the populations of susceptible individuals, infected individuals and removed individuals. After some initial transients the population densities fluctuate around the deterministic steady-state values S ¼ 0:625; I ¼ 0:55 and R ¼ 0:55. From the Figs. 1 and 2, with intensities increasing, the amplitude of fluctuation is larger. In Fig. 3, we choose r2 ¼ 0:7; r4 ¼ 0:6, other parameters take the same values as Fig. 1. As a result, the population of infected individuals goes to extinction. And, we can verify that the conditions in Theorem 4 are satisfied. That is to say, large noises can lead the disease to extinction, which is a phenomenon different from its corresponding deterministic system (1.1). Group(a) 1.5 0.2
I(t)
R(t)
0.2
S(t)
1 0.5
0
50
100
0
150
0.1
0
50
t
Group(b) 1.5
100
0.3
0.3
0.2
0.2
0.1
50
100
t
50
150
100
150
100
150
t
I(t)
S(t) 0.5
0
0
t
1
0
0
150
R(t)
0
0.1
0
0.1
0
50
100
t
150
0
0
50
t
Fig. 5. Deterministic and stochastic trajectories of system for initial condition and parameter values: Sð0Þ ¼ 0:7; Ið0Þ ¼ 0:2; Rð0Þ ¼ 0:1; K ¼ 0:2; b ¼ 0:4; l ¼ 0:2; e ¼ 0:1; c ¼ 0:2; Group (a): r1 ¼ 0:1; r2 ¼ 0:2; r3 ¼ 0:3; r4 ¼ 0:2; Group (b): r1 ¼ 0; r2 ¼ 0:1; r3 ¼ 0:1; r4 ¼ 0:1.
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
131
In the following, we consider the long time behavior of system (1.4) in the case of R0 < 1. In Figs. 4 and 5, we always choose initial value Sð0Þ ¼ 0:7; Ið0Þ ¼ 0:2, Rð0Þ ¼ 0:1 and parameters K ¼ 0:2; b ¼ 0:4; l ¼ 0:2; e ¼ 0:1; c ¼ 0:2 with different intensities of white noise which satisfy the conditions in Theorem 5.1. In Fig. 4(a), we choose r1 ¼ 0:1; r2 ¼ 0:1; r3 ¼ 0:1; r4 ¼ 0:1 and In Fig. 4(b), we choose r1 ¼ 0:05; r2 ¼ 0:05; r3 ¼ 0:05; r4 ¼ 0:05. As we can see, the curves of system (1.4) always oscillate around the curves of system (1.1) in time, and moreover, the larger the intensities of the white noises are, the larger the fluctuations of the solutions will be. Furthermore, if we take r1 the same value as in Fig. 5(a), but r2 ; r3 and r4 take larger values r2 ¼ 0:2; r3 ¼ 0:3; r4 ¼ 0:2. We can see from Fig. 5(a) that the fluctuation of the solution will be very smaller. That is to say, the intensities of r2 ; r3 and r4 have little effects on the solution. When r1 ¼ 0, then E0 becomes the disease-free equilibrium of system (1.4). We can see from Fig. 5(b) that solution of system (1.4) is stochastically asymptotically stable in the large, which supports the conclusion of Corollary 5.1. Acknowledgments The authors thank the editors and reviewers for their valuable comments. This research is supported by National Natural Science Foundation of China (11271260, 11071164) and the Innovation Program of Shanghai Municipal Education Commission (13ZZ116). Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/ j.amc.2014.06.100. References [1] W.O. Kermack, A.G. McKendrick, Contributions to the mathematical theory of epidemics (Part I), Proc. Soc. London Ser. A 115 (1927) 700–721. [2] J.M. Tchuenche, A. Nwagwo, R. Levins, Global behaviour of an SIR epidemic model with time delay, Math. Methods Appl. Sci. 30 (2007) 733–749. [3] T.L. Zhang, Z.D. Teng, Permanence and extinction for a nonautonomous SIRS epidemic model with time delay, Appl. Math. Model. 33 (2009) 1058– 1071. [4] F.P. Zhang, Z.Z. Li, F. Zhang, Global stability of an SIR epidemic model with constant infectious period, Appl. Math. Comput. 199 (2008) 285–291. [5] A. Korobeinikov, Lyapunov functions and global stability for SIR and SIRS epidemiological models with non-linear transmission, Bull. Math. Biol. 68 (2006) 615–626. [6] R.M. Anderson, R.M. May, Population biology of infectious diseases. Part I, Nature 280 (1979) 361–367. [7] R.M. May, Stability and Complexity in Model Ecosystems, Princeton University Press, New Jersey, 1973. [8] C.Y. Ji, D.Q. Jiang, N.Z. Shi, The behavior of an SIR epidemic model with stochastic perturbation, Stoch. Anal. Appl. 30 (2012) 755–773. [9] A. Gray, D. Greenhalgh, L. Hu, X.R. Mao, J.F. Pan, A stochastic differential equation SIS epidemic model, SIAM J. Appl. Math. 71 (2011) 876–902. [10] D.Q. Jiang, J.J. Yu, C.Y. Ji, N.Z. Shi, Asymptotic behavior of global positive solution to a stochastic SIR model, Math. Comput. Model. 54 (2011) 221–232. [11] J.R. Artalejo, A. Economou, M.J. Lopez-Herrero, On the number of recovered individuals in the SIS and SIR stochastic epidemic models, Math. Biosci. 228 (2010) 45–55. [12] Q. Lu, Stability of SIRS system with random perturbations, Physica A 388 (2009) 3677–3686. [13] Y.G. Lin, D.Q. Jiang, Long-time behavior of perturbed SIR model by white noise, Discrete Contin. Dyn. Syst. Ser. B 18 (2013) 1873–1887. [14] X.R. Mao, Stochastic Differential Equations and Applications, 2nd Edition., Horwood, Chichester, UK, 2007. [15] C.Y. Ji, D.Q. Jiang, Dynamics of a stochastic density dependent predator-prey system with Beddington–DeAngelis functional response, J. Math. Anal. Appl 381 (2011) 441–453. [16] R.Z. Hasminskii, Stochastic Stability of Differential Equations, Sijthoff and Noordhoff, 1980. [17] T.C. Gard, Introduction to Stochastic Differential Equations, Marcel Dekker, New York, 1987. [18] G. Strang, Linear Algebra and Its Applications, Thomson Learning Inc, 1988. [19] C. Zhu, G. Yin, Asymptotic properties of hybrid diffusion systems, SIAM J. Control Optim. 46 (2007) 1155–1179. [20] X.R. Mao, Almost sure asymptotic bounds for a class of stochastic differential equations, Stoch. Stoch. Rep. 41 (1992) 57–69. [21] D.J. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Rev. 43 (2001) 525–546.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Survival and stationary distribution of a SIR epidemic model with stochastic perturbations Yanli Zhou a,b,⇑, Weiguo Zhang a, Sanling Yuan a a b
University of Shanghai for Science and Technology, Shanghai 200093, China Shanghai Medical Instrumentation College, Shanghai 200093, China
a r t i c l e
i n f o
Keywords: Stochastic SIR epidemic model Stationary distribution Itô formula Extinction
a b s t r a c t In this paper, the dynamics of a SIR epidemic model is investigated. First, we show that the system admits a unique positive global solution starting from the positive initial value. Then, when R0 > 1, we show that the stochastic model has a stationary distribution under certain parametric restrictions. In particular, we show that random effects may lead the disease to extinction in scenarios where the deterministic model predicts persistence. When R0 6 1, a result on fluctuation of the solution around the disease-free equilibrium of deterministic system is established under suitable conditions. Finally, numerical simulations are carried out to illustrate the theoretical results. Ó 2014 Elsevier Inc. All rights reserved.
1. Introduction In 1927, Kermack and McKendrick establish the classical deterministic SIR (susceptible-infected-removed) model [1]. Since then, many people have studied the SIR disease model. Regarding research on the deterministic SIR model and its generalizations, the reader can refer to [2–6]. If SðtÞ; IðtÞ; RðtÞ denote the number of the individuals susceptible to the disease, of infected members and members who have been removed from the possibility of infection through full immunity, respectively, then the differential equations which describe the spread of the disease are:
8 _ > < SðtÞ ¼ K bSðtÞIðtÞ lSðtÞ; _ ¼ bSðtÞIðtÞ ðc þ l þ eÞIðtÞ; IðtÞ > :_ RðtÞ ¼ cIðtÞ lRðtÞ: The parameters in the model are summarized in the following list:
K: the influx of individuals into the susceptible; b: transmission coefficient between compartments S and I; l: natural death rate of S; I; R compartments; e: the disease induced mortality rate; c: the rate of recovery from infection; K; b; l; e; c are all positive. ⇑ Corresponding author at: University of Shanghai for Science and Technology, Shanghai 200093, China. E-mail address: [email protected] (Y. Zhou). http://dx.doi.org/10.1016/j.amc.2014.06.100 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.
ð1:1Þ
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
119
In system (1.1), R0 ¼ lðlbþKeþcÞ is the threshold of the system for an epidemic to persist or disappear. If R0 1, system (1.1) always has a disease-free equilibrium E0 ¼ ðKl ; 0; 0Þ and moreover, it is globally stable in int C, where
C ¼ fðS; I; RÞ : S > 0; I P 0; R P 0; S þ I þ R 6 Klg, this means that the disease dies out. If R0 > 1, then E0 is unstable and there exists a globally asymptotically stable endemic equilibrium E ¼ lþbeþc ; lþKeþc lb ; lðlKþceþcÞ bc , this means that the disease
remains. These above conclusions of system (1.1) can be obtained from [6]. Those important and useful works on deterministic models provide a great insight into the effect of the epidemic model, but in the real world, the models of population dynamics of diseases are inevitably affected by random fluctuations. Introduction of stochasticity into population dynamics models of diseases can bring to light new insights. May [7] pointed out that due to continuous fluctuation in the environment, the birth rates, death rates, transmission coefficient and all other parameters involved with the model exhibit random fluctuations to a greater or lesser extent. As a result the equilibrium distribution never attains a steady value, but fluctuates randomly around some average value [7]. Therefore stochastic differential equation models play a significant role in various branches of applied sciences including the population dynamics of diseases. As related work, the reader can refer to [8–12]. In [10], Jiang et al. assumed that stochastic perturbations are of a white noise type which are directly proportional _ _ _ to SðtÞ; IðtÞ; RðtÞ, influenced on the SðtÞ; IðtÞ; RðtÞ in system (1.1). Accordingly, system (1.1) was modified to the following form:
8 dSðtÞ ¼ ðK bSðtÞIðtÞ lSðtÞÞdt þ r1 SðtÞdB1 ðtÞ; > > < dIðtÞ ¼ ðbSðtÞIðtÞ ðc þ l þ eÞIðtÞÞdt þ r2 IðtÞdB2 ðtÞ; > > : dRðtÞ ¼ ðcIðtÞ lRðtÞÞdt þ r3 RðtÞdB3 ðtÞ;
ð1:2Þ
where Bi ðtÞ ði ¼ 1; 2; 3Þ is a standard Brownian motion, ri ði ¼ 1; 2; 3Þ is the intensity of environmental white noise. In [10], the long time behavior of the stochastic system was studied. In [13], Lin and Jiang assumed that the disease transmission coefficient b is subject to the environmental white noise, that _ is b ! b þ rBðtÞ. Then system (1.1) becomes Itô SDE
8 dSðtÞ ¼ ðK bSðtÞIðtÞ lSðtÞÞdt rSðtÞIðtÞdBðtÞ; > > < dIðtÞ ¼ ðbSðtÞIðtÞ ðc þ l þ eÞIðtÞÞdt þ rSðtÞIðtÞdBðtÞ; > > : dRðtÞ ¼ ðcIðtÞ lRðtÞÞdt;
ð1:3Þ
where BðtÞ is a standard Brownian motion, r is the intensity of environmental white noise. For system (1.3), Lin and Jiang presented sufficient conditions for the disease to go extinct exponentially and found the support of the invariant density. In particular, they obtained the existence of a stationary distribution of system (1.3) and its asymptotic stability by use of the Markov semigroup theory. Motivated by the works of Jiang et al. [10] and Lin and Jiang [13], in this paper, we consider the effects of random fluctuations from environmental white noise. The stochastic version based on the deterministic system (1.1) takes the following form:
8 dSðtÞ ¼ ðK bSðtÞIðtÞ lSðtÞÞdt þ r1 SðtÞdB1 ðtÞ r4 SðtÞIðtÞdB4 ðtÞ; > > < dIðtÞ ¼ ðbSðtÞIðtÞ ðc þ l þ eÞIðtÞÞdt þ r2 IðtÞdB2 ðtÞ þ r4 SðtÞIðtÞdB4 ðtÞ; > > : dRðtÞ ¼ ðcIðtÞ lRðtÞÞdt þ r3 RðtÞdB3 ðtÞ;
ð1:4Þ
where Bi ðtÞ ði ¼ 1; 2; 3; 4Þ is a standard Brownian motion with Bi ð0Þ ¼ 0 ði ¼ 1; 2; 3; 4Þ and ri ði ¼ 1; 2; 3; 4Þ is the intensity of environmental white noise. We are not aware of any literatures analytically concerning the existence of a stationary distribution of system (1.4). The main aim of this paper is to study the existence of a stationary distribution of system (1.4). The paper is organized as follows. In Section 2, the global existence and positivity of the solution of system (1.4) are studied. In Section 3, by choosing appropriate Lyapunov functions, we show that there is a stationary distribution for system (1.4) under certain parametric restrictions. In Section 4, sufficient conditions are derived for the extinction of disease. In Section 5, how the solution spirals around the disease-free equilibrium of deterministic system (1.1) is established. Finally, numerical simulations are carried out to illustrate the theoretical results. Throughout of this paper, unless otherwise specified, let ðX; F; fFt gtP0 ; PÞ be a complete probability space with a filtration satisfying the usual conditions (i.e., it is right continuous and F0 contains all P-null sets). Let Bi ðtÞ ði ¼ 1; 2; 3; 4Þ denote the independent standard Brownian motions defined on this probability space. We also denote
R3þ ¼ fðx; y; zÞ 2 R3 : x > 0; y > 0; z > 0g:
120
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
2. Existence and uniqueness of the positive solution Since the SðtÞ; IðtÞ; RðtÞ in Eq. (1.4) are the sizes of the susceptible individuals, the infected individuals and the removed individuals at time t, respectively, they should be non-negative. It is obvious that we are only interested in positive solutions when we study the epidemic systems. In order to get a unique global solution (i.e., no explosion in a finite time) for any given initial value, the coefficients of the stochastic differential equation are generally required to satisfy the linear growth condition and local Lipschitz condition [14]. However, the coefficients of Eq. (1.4) do not satisfy the linear growth condition, though they are locally Lipschitz continuous, which may include the solution of Eq. (1.4) to explode at a finite time. In what follows, we shall prove that the solution of system (1.4) is positive and global under given conditions. Theorem 2.1. For any given initial condition ðSð0Þ; Ið0Þ; Rð0ÞÞ 2 R3þ , there is a unique positive solution ðSðtÞ; IðtÞ; RðtÞÞ to system (1.4) on t P 0, and the solution remains in R3þ with probability one, namely, ðSðtÞ; IðtÞ; RðtÞÞ 2 R3þ for all t P 0 almost surely. Proof. Since the coefficients of the equation are locally Lipschitz continuous, for any initial value ðSð0Þ; Ið0Þ; Rð0ÞÞ 2 R3þ , there is a unique local solution ðSðtÞ; IðtÞ; RðtÞÞ on t 2 ½0; se Þ, where se is the explosion time (i.e., the moment at which the solution tends to infinity). Next, we show that se ¼ þ1 almost surely. Let m0 > 0 be sufficiently large so that Sð0Þ; Ið0Þ and Rð0Þ lie h i within the interval m10 ; m0 . For each integer m P m0 , define the stopping time
sm ¼ inf t 2 ½0; se Þ : minfSðtÞ; IðtÞ; RðtÞg 6
1 or maxfSðtÞ; IðtÞ; RðtÞg P m ; m
where throughout this paper we set inf / ¼ 1 (as usual / denotes the empty set). According to the definition, sm is increasing as m ! 1. Set s1 ¼ limm!1 sm , whence s1 6 sm a.s. If we can show that s1 ¼ 1 a.s., then se ¼ 1 and ðSðtÞ; IðtÞ; RðtÞÞ 2 R3þ a.s. for all t P 0. If this statement is false, then there exist a pair of constants T > 0 and 2 ð0; 1Þ such that
Pðs1 6 TÞ P : Hence there is an integer m1 P m0 such that
Pfsm 6 Tg P for all m P m1 :
ð2:1Þ
Consider the function VðtÞ by
SðtÞ þ ðIðtÞ 1 log IðtÞÞ þ ðRðtÞ 1 log RðtÞÞ; VðSðtÞ; IðtÞ; RðtÞÞ ¼ SðtÞ c c log c where c is a positive constant to be defined later. The nonnegativity of this function can be seen from the fact that for any fixed constant v > 0; f ðuÞ ¼ u v logðu=v Þ P 0; for all u > 0. With Itô formula, we can get
c c 1 1 1 2 2 dSðtÞ þ dIðtÞ þ dRðtÞ dVðSðtÞ; IðtÞ; RðtÞÞ ¼ 1 ðdSðtÞÞ þ 1 ðdIðtÞÞ þ 1 SðtÞ IðtÞ RðtÞ 2SðtÞ2 2IðtÞ2 þ
1 2RðtÞ2
2
ðdRðtÞÞ ;
then
dVðSðtÞ; IðtÞ; RðtÞÞ ¼ LVðSðtÞ; IðtÞ; RðtÞÞdt þ ðSðtÞ cÞðr1 dB1 ðtÞ r4 IðtÞdB4 ðtÞÞ þ ðIðtÞ 1Þðr2 dB2 ðtÞ þ r4 SðtÞdB4 ðtÞÞ þ ðRðtÞ 1Þr3 dB3 ðtÞ; where LVðSðtÞ; IðtÞ; RðtÞÞ :
R3þ
LVðSðtÞ; IðtÞ; RðtÞÞ ¼
ð2:2Þ
! Rþ by
c cr2 þ r24 I2 1 ðK bSðtÞIðtÞ lSðtÞÞ þ 1 ðbSðtÞIðtÞ ðl þ c þ eÞIðtÞÞ þ 1 SðtÞ IðtÞ 2 r2 þ r24 S2 1 r2 þ 2 ðcI lRÞ þ 3 þ 1 RðtÞ 2 2 1
¼ K þ cl þ l þ l þ e þ c þ c
cr21 þ r24 I2 ðtÞ cK þ ½cb ðl þ eÞIðtÞ ðl þ bÞSðtÞ lRðtÞ 2 SðtÞ
IðtÞ r22 þ r24 S2 ðtÞ r23 þ : þ RðtÞ 2 2
e Choose c ¼ lþ such that ½cb ðl þ eÞIðtÞ ¼ 0, then b
ð2:3Þ
121
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
LVðtÞ 6 K þ cl þ 2l þ c þ e þ
cr21 þ r22 þ r24 I2 ðtÞ þ r24 S2 ðtÞ þ r23 ¼: K: 2
ð2:4Þ
Therefore,
dVðSðtÞ;IðtÞ;RðtÞÞ 6 KdtðSðtÞcÞðr1 dB1 ðtÞ r4 IðtÞdB4 ðtÞÞþðIðtÞ1Þðr2 dB2 ðtÞþ r4 SðtÞdB4 ðtÞÞþðRðtÞ1Þr3 dB3 ðtÞ: ð2:5Þ Integrating both sides of (2.5) from 0 to
Z sm ^T
dVðSðtÞ; IðtÞ; RðtÞÞ 6
Z sm ^T
0
sm ^ T
Kdt þ
0
Z sm ^T 0
þ r4 SðtÞdB4 ðtÞÞ þ
ðSðtÞ cÞðr1 dB1 ðtÞ r4 IðtÞdB4 ðtÞÞ þ
Z sm ^T
Z sm ^T
ðIðtÞ 1Þðr2 dB2 ðtÞ
0
ðRðtÞ 1Þr3 dB3 ðtÞ:
ð2:6Þ
0
Taking the expectations on both sides of (2.6), we obtain
EVðSðsm ^ TÞ; Iðsm ^ TÞ; Rðsm ^ TÞÞ 6 VðSð0Þ; Ið0Þ; Rð0ÞÞ þ KEðsm ^ TÞ; so
EVðSðsm ^ TÞ; Iðsm ^ TÞ; Rðsm ^ TÞÞ 6 VðSð0Þ; Ið0Þ; Rð0ÞÞ þ KT:
ð2:7Þ
Let Xm ¼ sm ^ T for m P m1 . Then, by (2.1), we known that PðXm Þ P . Note that, for every x 2 Xm , there is at least one of Sðsm ^ TÞ; Iðsm ^ TÞ and Rðsm ^ TÞ that equals either m or m1 . Then
VðSð0Þ; Ið0Þ; Rð0ÞÞ þ KT P E½IXm ðxÞVðSðsm ^ TÞ; Iðsm ^ TÞ; Rðsm ^ TÞÞ ¼ E½IXm ðxÞVðSðsm ; xÞ; Iðsm ; xÞ; Rðsm ; xÞÞ 1 1 þ log m ; P ðm 1 log mÞ ^ m where IXm ðxÞ is the indicator function of Xm ðxÞ. Letting m ! 1, it leads to the contradiction
1 > VðSð0Þ; Ið0Þ; Rð0ÞÞ þ KT ¼ 1; which is a contradiction, therefore we have s1 ¼ 1 a.s. Therefore, it implies SðtÞ; IðtÞ and RðtÞ will not explode in a finite time with probability one. The proof is therefore completed. h
3. Existence of the stationary distribution An important aspect on studying the disease is the long term prevalence of the disease or whether it will ultimately vanish from the population. In this subsection, we show that there is a stationary distribution when the white noise is small. The stationary distribution proofs in this paper are similar to [15]. Before proving the main theorem of this section we first cite a known result from Has’miskii [16] which will be useful to prove the theorem. Let XðtÞ be a regular time-homogeneous Markov process in El (the l dimensional Euclidean space) described by stochastic equation
dXðtÞ ¼ bðXÞdt þ
k X
f r dBr ðtÞ:
ð3:1Þ
r¼1
The diffusion matrix is defined as follows:
AðxÞ ¼ ðaij ðxÞÞ;
aij ðxÞ ¼
k1 X
fri ðxÞfrj ðxÞ:
r¼1
There exists a bounded domain U El with regular boundary C, having the following properties: (B.1) In the domain U and some neighborhood thereof, the smallest eigenvalue of the diffusion matrix AðxÞ is bounded away from zero. (B.2) If x 2 El n U, the mean time s at which a path emerging from x reaches the set U is finite, and supx2K Ex s < 1 for every compact subset K El .
Lemma 3.1. If above assumptions (B.1) and (B.2) hold, then the Markov process XðtÞ has a stationary distribution lðÞ. Let f ðÞ be a integrable function with respect to the measure lðÞ. Then
122
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
( Px lim
T*1
Z
T
f ðXðtÞÞdt ¼
Z
0
) f ðxÞlðdxÞ
¼1
El
for all x 2 El . Remark 3.1. Lemma 3.1 asserts the existence of a stationary distribution with suitable density function. Its proof can be found in [16]. To validate (B.1), it suffices to prove V is uniformly elliptical in U, where Vu ¼ bðxÞ:ux þ ½trðAðxÞuxx =2, i.e. there is a positive number M such that k X ai;j ni nj P Mjnj2 ;
x 2 U; n 2 Rk
i;j¼1
(see [17] and Rayleigh’s principle in [18]). To verify (B.2), it is sufficient to show that there exists some neighborhood U and a non-negative C 2 -function such that LV is negative at every point x 2 El n U like in [19]. Theorem 3.1. Consider the stochastic system (1.4) with initial condition in R3þ . Assume that R0 > 1 and 0 < d < minðm1 S2 ; m2 I2 ; m3 R2 Þ, then there exists a stationary distribution lðÞ for system (1.4). Here ðS ; I ; R Þ is the unique endemic equilibrium of system (1.1),
m1 ¼ l r21
r24 ð2l þ c þ eÞI b
m2 ¼ l þ c r22
;
2
c
2l
;
l
r23 ; 2l þ c þ e 2 r24 I ð2l þ e þ cÞ 2 d ¼ r21 S2 þ I2 þ I r2 þ S þ r23 R2 : b 2b
m3 ¼
ð3:2Þ
2
Especially, we have
1 lim sup E t t!1
Z
t
2
2
2
½m1 ðSðuÞ S Þ þ m2 ðIðuÞ I Þ þ m3 ðRðuÞ R Þ du 6 d:
ð3:3Þ
0
Proof. Since R0 > 1, then there exists a positive equilibrium of system (1.1), and
K ¼ bS I þ lS ; bS I ¼ ðc þ l þ eÞI ; cI ¼ lR :
ð3:4Þ
Let us consider the function
VðSðtÞ; IðtÞ; RðtÞÞ ¼ V 1 ðtÞ þ
2l þ c þ e V 2 ðtÞ þ V 3 ðtÞ; b
ð3:5Þ
where
V 1 ðtÞ ¼
1 2 ðSðtÞ þ IðtÞ S I Þ ; 2
V 2 ðtÞ ¼ IðtÞ I I ln
IðtÞ ; I
2
V 3 ðtÞ ¼
ðRðtÞ R Þ : 2
ð3:6Þ
Using Itô formula, we have
dVðtÞ ¼ dV 1 ðtÞ þ
2l þ c þ e dV 2 ðtÞ þ dV 3 ðtÞ; b
in detail,
1 2 dV 1 ðtÞ ¼ ðSðtÞ þ IðtÞ S I ÞðdSðtÞ þ dIðtÞÞ þ ðdSðtÞ þ dIðtÞÞ 2 ¼ LV 1 ðtÞdt þ ðSðtÞ þ IðtÞ S I Þðr1 SðtÞdB1 ðtÞ þ r2 IðtÞdB2 ðtÞÞ; I I 2 dV 2 ðtÞ ¼ 1 dIðtÞ þ 2 ðdIðtÞÞ ¼ LV 2 dt þ ðIðtÞ I Þr2 dB2 ðtÞ þ ðIðtÞ I Þr4 SðtÞdB4 ðtÞ; IðtÞ 2I ðtÞ 2
dV 3 ðtÞ ¼ ðRðtÞ R ÞdRðtÞ þ ðdRðtÞÞ ¼ LV 3 ðtÞdt þ ðRðtÞ R ÞRðtÞr3 dB3 ðtÞ;
ð3:7Þ
123
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
where
LV 1 ðtÞ ¼ ðSðtÞ þ IðtÞ S I Þ½K lSðtÞ ðl þ c þ eÞIðtÞ þ
r21 S2 ðtÞ þ r22 I2 ðtÞ 2
LV 2 ðtÞ ¼ ðIðtÞ I Þ½bSðtÞ ðl þ c þ eÞ þ
;
I ðr2 þ r24 S2 ðtÞÞ; 2 2
ð3:9Þ
1 LV 3 ðtÞ ¼ ðRðtÞ R ÞðcIðtÞ lRðtÞÞ þ r23 R2 ðtÞ: 2 2
ð3:8Þ
ð3:10Þ
2
Using (3.4) and the inequality ða þ bÞ 2a2 þ 2b , we can have from (3.8) that
1 2 LV 1 ðtÞ ¼ ðSðtÞ S þ IðtÞ I Þ½lS þ ðl þ c þ eÞI lSðtÞ ðl þ c þ eÞIðtÞ þ ½r21 ðSðtÞ S þ S Þ 2 2 2 þ r22 ðIðtÞ I þ I Þ 6 ðSðtÞ S þ IðtÞ I Þ½lðSðtÞ S Þ ðl þ c þ eÞðIðtÞ I Þ þ r21 ðS S Þ 2
2
þ r21 S2 þ r22 ðI I Þ þ r22 I2 ¼ ðl r21 ÞðSðtÞ S Þ ð2l þ c þ eÞðSðtÞ S ÞðIðtÞ I Þ 2
ðl þ c þ e r22 ÞðIðtÞ I Þ þ r21 S2 þ r22 I2 :
ð3:11Þ
Similarly, we can have from (3.9) and (3.10) that
LV 2 ðtÞ ¼ ðIðtÞ I ÞðbSðtÞ bS Þ þ 6 bðSðtÞ S ÞðIðtÞ I Þ þ
I 2 r24 I 2 r þ ðSðtÞ S þ S Þ 2 2 2
I 2 r þ r24 I ðSðtÞ S Þ2 þ r24 I S2 2 2
ð3:12Þ
and
LV 3 ðtÞ ¼ ðRðtÞ R ÞðcIðtÞ lRðtÞ cI þ lR Þ þ
1 2 r ðRðtÞ R þ R Þ2 2 3
2
2
6 cðRðtÞ R ÞðIðtÞ I Þ lðRðtÞ R Þ þ r23 ðRðtÞ R Þ þ r23 R pffiffiffiffi c 2 2 ¼ pffiffiffiffi ðIðtÞ I Þ lðRðtÞ R Þ lðRðtÞ R Þ þ r23 ðRðtÞ R Þ þ r23 R2
l
2
c
6
2l
2
ðIðtÞ I Þ
l 2
2
2
ðRðtÞ R Þ þ r23 ðRðtÞ R Þ þ r23 R2 :
ð3:13Þ
It follows from (3.7) that
LVðtÞ ¼ LV 1 ðtÞ þ
2l þ c þ e LV 2 ðtÞ þ LV 3 ðtÞ: b
ð3:14Þ
Substituting (3.11)–(3.13) into (3.14), we obtain that
l r2 ð2l þ c þ eÞI c2 2 2 2 ðSðtÞ S Þ l þ c r22 ðIðtÞ I Þ LVðtÞ 6 l r21 4 r23 ðRðtÞ R Þ b 2l 2 2l þ c þ e 2 r24 I ð2l þ e þ cÞ 2 2 2 þ r21 S2 þ I2 þ I r2 þ S þ r23 R2 , m1 ðSðtÞ S Þ m2 ðIðtÞ I Þ b 2b 2
m3 ðRðtÞ R Þ þ d;
ð3:15Þ
where m1 ; m2 ; m3 and d are defined respectively in (3.2). Therefore, 2
2
2
dVðtÞ 6 m1 ðSðtÞ S Þ m2 ðIðtÞ I Þ m3 ðRðtÞ R Þ þ d þ ðSðtÞ þ IðtÞ S I Þðr1 SðtÞdB1 ðtÞ
þ r2 IðtÞdB2 ðtÞÞ þ ðIðtÞ I Þr2 dB2 ðtÞ þ ðIðtÞ I Þr4 SðtÞdB4 ðtÞ þ ðRðtÞ R ÞRðtÞr3 dB3 ðtÞ:
ð3:16Þ
Integrating both sides of (3.16) from 0 to t, yields
VðtÞ Vð0Þ 6
Z
t
2
2
2
½m1 ðSðuÞ S Þ m2 ðIðuÞ I Þ m3 ðRðuÞ R Þ du þ dt þ MðtÞ;
ð3:17Þ
0
where MðtÞ is a martingale defined by
MðtÞ ¼
Z 0
þ
t
ðSðuÞ þ IðuÞ S I Þðr1 SðuÞdB1 ðuÞ þ r2 IðuÞdB2 ðuÞÞ þ
Z 0
Z 0
t
ðRðuÞ R ÞRðuÞr3 dB3 ðuÞ:
t
ððIðuÞ I Þr2 dB2 ðuÞ þ ðIðuÞ I Þr4 SðuÞdB4 ðuÞÞ
124
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
Taking expectation on both sides of (3.17), yields
EVðtÞ EVð0Þ 6 E
Z
t
2
2
2
½m1 ðSðuÞ S Þ m2 ðIðuÞ I Þ m3 ðRðuÞ R Þ du þ dt:
0
Hence,
1 lim sup E t t!1
Z
t
2
2
2
½m1 ðSðuÞ S Þ þ m2 ðIðuÞ I Þ þ m3 ðRðuÞ R Þ du 6 d:
0
That is, we obtain (3.3). Noticing that if 0 < d < minðm1 S2 ; m2 I2 ; m3 R2 Þ, then the ellipsoid 2
2
2
m1 ðSðtÞ S Þ m2 ðIðtÞ I Þ m3 ðRðtÞ R Þ þ d ¼ 0 lies entirely in R3þ . We can take U to be any neighborhood of the ellipsoid such that U # El ¼ R3þ , so for x 2 U n El ; LVðtÞ 6 C (C is a positive constant), which implies condition (B.2) is satisfied. On the other hand, we can write system (1.4) as the form of system (3.1):
0
1 0 1 0 1 0 1 K bSðtÞIðtÞ lSðtÞ þ cIðtÞ r1 SðtÞ 0 0 B C B C B C B C B C d@ IðtÞ A ¼ @ bSðtÞIðtÞ ðl þ e þ cÞIðtÞ Adt þ @ 0 AdB1 ðtÞ þ @ r2 IðtÞ AdB2 ðtÞ þ @ 0 AdB3 ðtÞ RðtÞ 0 0 cIðtÞ lRðtÞ r3 RðtÞ 0 1 r4 SðtÞIðtÞ B C þ @ r4 SðtÞIðtÞ AdB4 ðtÞ: 0 SðtÞ
1
0
Here the diffusion matrix is
0
r2 S2 þ r24 S2 I2 r24 S2 I2 B 1 A¼B r22 I2 þ r24 S2 I2 r24 S2 I2 @ 0
0 0
r23 R2
0
1 C C: A
There is a M ¼ minfr21 S2 ; r22 I2 ; r23 R2 g > 0, such that for all ðx1 ; x2 ; x3 Þ 2 U and n 2 R3 , 3 X aij ni nj ¼ r21 S2 þ r24 S2 I2 n21 þ r22 I2 þ r24 S2 I2 n22 2r24 S2 I2 n1 n2 þ r23 R2 n23 i;j¼1
¼ r21 S2 n21 þ r22 I2 n22 þ r23 R2 n23 þ r24 S2 I2 ðn1 n2 Þ2 P r21 S2 n21 þ r22 I2 n22 þ r23 R2 n23 P min
n
o
r21 S2 ; r22 I2 ; r23 R2 jnj2
¼ Mjnj2 ; which shows that condition (B.1) is also satisfied. Therefore, we can conclude that stochastic system (1.4) has a stationary distribution lðÞ. h Next, we can further have the following almost sure persistence of the total population result of system (1.4). Set
XðtÞ ¼ SðtÞ þ IðtÞ þ RðtÞ
ð3:18Þ
Theorem 3.2. If K > 0, then for any given initial ðSð0Þ; Ið0Þ; Rð0ÞÞ 2 R3þ , the solution of model (1.4) obeys
0 < lim inf XðtÞ 6 lim sup XðtÞ < þ1: t!1
ð3:19Þ
t!1
Proof. It follows from (3.18) that
dXðtÞ ¼ ðK lXðtÞ eIðtÞÞdt þ r1 SðtÞdB1 ðtÞ þ r2 IðtÞdB2 ðtÞ þ r3 RðtÞdB3 ðtÞ: 1 Define a function VðXðtÞÞ ¼ XðtÞ and by the Itô formula, we have
d
1 1 1 2 ¼ 2 dXðtÞ þ 3 ðdXðtÞÞ XðtÞ X ðtÞ X ðtÞ ¼
1 2
X ðtÞ
ðK lXðtÞ eIðtÞÞdt þ
þ r3 RðtÞdB3 ðtÞÞ:
1 3
X ðtÞ
ðr21 S2 þ r22 I2 þ r23 R2 Þ
1 2
X ðtÞ
ðr1 SðtÞdB1 ðtÞ þ r2 IðtÞdB2 ðtÞ
ð3:20Þ
125
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
Since
SðtÞ XðtÞ
d
IðtÞ 1; XðtÞ 1 and
"
RðtÞ XðtÞ
1, then
# 1 t K l þ e r21 þ r22 þ r23 et e ¼ 2 þ þ þ 1 et dt 2 ðr1 SðtÞdB1 ðtÞ þ r2 IðtÞdB2 ðtÞ þ r3 RðtÞdB3 ðtÞÞ: XðtÞ XðtÞ X ðtÞ XðtÞ X ðtÞ
It is easy to obtain that
1 1 t ¼ e þ et XðtÞ Xð0Þ
Z
t
eu HðXðuÞÞdu et
"Z
0
0
t
eu X 2 ðuÞ
r1 SðuÞdB1 ðuÞ þ
Z 0
t
eu X 2 ðuÞ
r2 IðuÞdB2 ðuÞ þ
Z
t
0
eu X 2 ðuÞ
#
r3 RðuÞdB3 ðuÞ : ð3:21Þ
Obviously, if K > 0, then we have HðXðtÞÞ 6 M1 for all XðtÞ > 0. In fact, we only need to prove the next two equations
Pðlim inf Xðx; tÞ ¼ 0Þ ¼ 0
ð3:22Þ
Pðlim sup Xðx; tÞ ¼ þ1Þ ¼ 0:
ð3:23Þ
t!1
and t!1
We now begin to prove assertion (3.22). If it is not true, then there is a sufficiently small d 2 ð0; 1Þ such that
PðX1 Þ > d;
ð3:24Þ
where X1 ¼ fxj lim inf Xðx; tÞ ¼ 0g. Define the stopping time
t!1
1 k
sk ¼ inf t P 0j; Xðx; tÞ 6 ; x 2 X1 ; k 2 Z þ : According to the definition,
sk is increasing and sk ! 1 as k ! 1. It then follows from (3.21) that
Z sk 1 1 sk 1 þ M1 : E I X1 ¼ E e þ esk eu HðXÞdu 6 Xðsk Þ Xð0Þ Xð0Þ 0
ð3:25Þ
However,
1 E IX1 P kEðIX1 Þ P kd: Xðsk Þ Letting k ! 1, leads to yield
1 E IX1 P kEðIX1 Þ P kd ! 1; Xðsk Þ
as k ! 1:
But this contradicts (3.25). We therefore must have the desired assertion (3.22), namely, lim inf t!1 XðtÞ > 0 a.s. On the other hand, it follows from (3.20) that Z t Z t Z t Z t XðtÞ ¼ Xð0Þelt þ elt elu ðK eIðuÞÞdu þ elt elu r1 SðuÞdB1 ðuÞ þ elt elu r2 IðuÞdB2 ðuÞ þ elt elu r3 RðuÞdB3 ðuÞ: 0
0
0
0
In the same way as in the proof of (3.22), we find that
PðX2 Þ ¼ 0; where X2 ¼ fxj lim supt!1 Xðx; tÞ ¼ þ1g a.s.. That is, lim supt!1 XðtÞ < þ1. The proof is hence complete.
h
4. Extinction of the disease In the deterministic models, the value of the basic reproductive number R0 determines persistence or extinction of the disease. If R0 1, the disease is eliminated, whereas if R0 > 1, the disease persists in the population. However, we will show in this section that if the noise is sufficiently large, the disease will become extinct for stochastic system (1.4), although it may be persistent for its deterministic version (1.1). The following theorem gives a condition for the extinction of a disease expressed in terms of intensities of noise and system parameters. Theorem 4.1. Consider stochastic system (1.4) with initial condition in R3þ . We have
lim sup t!1
If
b2 2r24
ln IðtÞ b2 r2 < ðc þ l þ eÞ 2 a:s: t 2 2r24
< ð c þ l þ eÞ þ
r22 2
holds, then I (t) will go to zero exponentially with probability one.
126
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
Before proving the main theorem of this section, we first cite a Lemma from Mao [20] which will be useful to prove the theorem. Lemma 4.1. Let MðtÞ : t P 0 be a continuous local martingale with value Mð0Þ ¼ 0. Let hMðtÞi be its quadratic variation. Let d > 1 be a number, and let mn and sn be two sequences of positive numbers. Then, for almost all x 2 X, there exists a random integer n0 ¼ n0 ðxÞ such that, for all n P n0 ,
MðtÞ 6
1 d ln n mn hMðtÞ; MðtÞi þ ; 2 mn
0 6 t 6 sn :
ð4:1Þ
Proof. Define VðIðtÞÞ ¼ ln IðtÞ and by Itô formula, we have
# Z t" Z t r22 r24 S2 ðuÞ du þ r2 B2 ðtÞ þ ln IðtÞ ¼ ln Ið0Þ þ bSðuÞ ðc þ l þ eÞ r4 SðuÞdB4 ðtÞ 2 2 0 0
Z t r2 r2 ¼ ln Ið0Þ þ 4 S2 ðuÞ þ bSðuÞ ðc þ l þ eÞ 2 du þ r2 B2 ðtÞ þ MðtÞ; 2 2 0 where MðtÞ ¼
Rt 0
r4 SðuÞdB4 ðuÞ is a continuous local martingale whose quadratic variation is
hMðtÞ; MðtÞi ¼ r24
1 2 mr 2 4
Z
t
S2 ðuÞdu:
0
mn ¼ m > 0 and sn ¼ n in (4.1). For almost all x 2 X, there is a random integer n0 ðxÞ such that, for n > n0 ,
Choosing d ¼ 2;
MðtÞ 6
ð4:2Þ
Z
t
S2 ðuÞdu þ
0
2
m
ln n for all 0 6 t 6 n:
ð4:3Þ
Substituting (4.3) into (4.2), we obtain that
ln IðtÞ < ln Ið0Þ þ
Z t Z t r2 r2 1 2 4 S2 ðuÞ þ bSðuÞ ðc þ l þ eÞ 2 du þ r2 B2 ðtÞ þ mr24 S2 ðuÞdu þ ln n: 2 2 2 m 0 0
ð4:4Þ
Clearly, we have that
2 1 2b 1 b b2 b2 ð1 mÞr24 S2 ðtÞ 2 SðtÞ ¼ ð1 mÞr24 SðtÞ þ 6 : 2 2 ð1 mÞr24 r4 2ð1 mÞr24 2ð1 mÞr24
ð4:5Þ
It follows from (4.4) and (4.5) that
" ln IðtÞ < ln Ið0Þ þ
# b2 r22 2 ð c þ l þ e Þ t þ r2 B2 ðtÞ þ ln n 2 m 2ð1 mÞr24
ð4:6Þ
Thus, for n 1 6 t 6 n, we obtain
ln IðtÞ ln Ið0Þ b2 r2 r2 B2 ðtÞ 2 ln n < þ ðc þ l þ eÞ 2 þ þ : t t t 2 m n1 2ð1 mÞr24 By the large number theorem for martingales [14], we have
lim
t!1
B2 ðtÞ ¼ 0: t
Hence, let n ! 1, so t ! 1, we obtain
lim sup t!1
Let
ln IðtÞ b2 r2 < ðc þ l þ eÞ 2 : 2 t 2 2ð1 mÞr4
m ! 0, we obtain lim sup t!1
ln IðtÞ b2 r2 < ðc þ l þ eÞ 2 : 2 t 2 2r4
The proof is therefore completed. h Remark 4.1. When R0 > 1, the positive solution converges to the endemic equilibrium of deterministic system (1.1). However, the disease will die2 out exponentially regardless of the magnitude of R0 provided that r2 and r4 are big enough such 2 r that 2br2 < ðc þ l þ eÞ þ 22 . That is to say, large noises can lead to the extinction of disease. 4
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
127
5. Asymptotic behavior around the disease-free equilibrium E0 As mentioned in system (1.1), if R0 1, then system (1.1) always has a globally asymptotically stable disease-free equilibrium E0 , which means the disease will die out with the advancement of time. Noting that E0 is not an equilibrium of stochastic system (1.4), it is natural to ask whether the disease will go to extinction in the population. In this section we mainly use the way of estimating the oscillation around E0 to reflect how the solution of system (1.4) spirals closely around E0 . We have the following theorem. Theorem 5.1. Assume R0 1. Then for any solution ðSðtÞ; IðtÞÞ of system (1.4) with initial value ðSð0Þ; Ið0Þ; Rð0ÞÞ 2 R3þ , if r21 < l; r22 < 2ðl þ eÞ þ c and r23 < 2l c, we have
lim sup t!1
1 E t
# Z t " K 2 r 2 K2 SðuÞ þ I2 ðuÞ þ R2 ðuÞ du 6 1 2 ; l ml 0
ð5:1Þ
n o r2 r2 where m ¼ min l r21 ; l þ e þ 2c 22 ; l 2c 23 . Proof. Define a function V : R3þ ! Rþ by
VðtÞ ¼ V 1 ðtÞ þ
2l þ e þ c V 2 ðtÞ þ V 3 ðtÞ ¼ b
2 SðtÞ Kl þ IðtÞ 2
þ
2l þ e þ c 1 IðtÞ þ R2 ðtÞ: 2 b
ð5:2Þ
Along the trajectories of system (1.4), we have
dVðtÞ ¼ dV 1 ðtÞ þ
2l þ e þ c dV 2 ðtÞ þ dV 3 ðtÞ: b
ð5:3Þ
In detail,
K 1 K 2 dV 1 ðtÞ ¼ SðtÞ þ IðtÞ ðdSðtÞ þ dIðtÞÞ þ ðdSðtÞ þ dIðtÞÞ ¼ LV 1 ðtÞdt þ SðtÞ þ IðtÞ ðr1 SðtÞdB1 ðtÞ þ r2 IðtÞdB2 ðtÞÞ; 2 l l dV 2 ðtÞ ¼ LV 2 ðtÞdt þ r2 IðtÞdB2 ðtÞ þ r4 SðtÞIðtÞdB4 ðtÞ; and
dV 3 ðtÞ ¼ LV 3 ðtÞdt þ r3 R2 ðtÞdB3 ðtÞ; where
K r2 S2 ðtÞ þ r22 I2 ðtÞ LV 1 ðtÞ ¼ SðtÞ þ IðtÞ ðK lSðtÞ ðl þ e þ cÞIðtÞÞ þ 1 ; 2 l
ð5:4Þ
LV 2 ðtÞ ¼ bSðtÞIðtÞ ðl þ c þ eÞIðtÞ;
ð5:5Þ
LV 3 ðtÞ ¼ cRðtÞIðtÞ lR2 ðtÞ þ
r
2 2 3 R ðtÞ
2
ð5:6Þ
:
From (5.3), we have that
2
r2 SðtÞ K þ K þ r2 I2 ðtÞ 1 2 l l K K ðl þ e þ cÞIðtÞ þ LV 1 ðtÞ ¼ SðtÞ þ IðtÞ l SðtÞ l l 2 2 K 2 K r2 K IðtÞ l þ e þ c 2 I2 ðtÞ þ r21 ð2l þ e þ cÞ SðtÞ ; 6 ðl r21 Þ SðtÞ l l 2 l 2
ð5:7Þ
2
where in the last step we have used the inequality ða þ bÞ 2a2 þ 2b . Similarly, we can have from (5.4) and (5.5) that
K K K K IðtÞ ðl þ c þ eÞIðtÞ ¼ b SðtÞ IðtÞ þ b ðl þ c þ eÞ IðtÞ LV 2 ðtÞ ¼ b SðtÞ þ
l
l
l
l
and
LV 3 ðtÞ 6
c 2
Noting that R0 ¼
½R2 ðtÞ þ I2 ðtÞ lR2 ðtÞ þ Kb lðlþeþcÞ
r23 R2 ðtÞ 2
c r2 c ¼ l 3 R2 ðtÞ þ I2 ðtÞ: 2 2 2
1, it follows from above that
K IðtÞ: LV 2 ðtÞ 6 b SðtÞ
l
ð5:8Þ
128
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
It follows from (5.6), (5.7) and (5.2) that
2 K 2 c r2 c r2 K LVðtÞ 6 ðl r21 Þ SðtÞ l þ e þ 2 I2 ðtÞ l 3 R2 ðtÞ þ r21 : l 2 2 2 2 l
ð5:9Þ
Integrating from 0 to t on both sides of (5.3), then taking expectation, and considering inequality (5.9), yields
EVðtÞ Vð0Þ 6 ðl r21 ÞE þ r21
2 K
l
Z t Z t Z t K 2 c r2 c r2 SðuÞ du l þ e þ 2 E I2 ðuÞdu l 3 E R2 ðuÞdu l 2 2 2 2 0 0 0
t:
Hence, we have that
1 lim sup E t t!1 Let m ¼ min
# 2 Z t" K 2 c r2 c r2 K ðl r21 Þ SðuÞ þ l þ e þ 2 I2 ðuÞ þ l 3 R2 ðuÞ du 6 r21 : l 2 2 2 2 l 0
n
2
2
ð5:10Þ
o
l r21 ; l þ e þ 2c r22 ; l 2c r23 , it follows from (5.9) that
1 lim sup E t t!1
# Z t " K 2 r 2 K2 2 2 SðuÞ þ I ðuÞ þ R ðuÞ du 6 1 2 : l ml 0
The proof of Theorem 5.1 is thus completed. h Notice that when r1 ¼ 0; E0 is also the disease-free equilibrium of system (1.4) and (5.9) is reduced to the following form
K 2 c r2 c r2 LVðtÞ 6 l SðtÞ l þ e þ 2 I2 ðtÞ l 3 R2 ðtÞ; l 2 2 2 2 which means that LVðtÞ is negative definite provided that r22 < 2ðl þ eÞ þ c and r23 < 2l c. Therefore, we can have the following result [14]. Corollary 5.1. Assume that R0 1 and r1 ¼ 0. Then the disease-free equilibrium E0 of system (1.4) is stochastically asymptotically stable in the large provided r22 < 2ðl þ eÞ þ c and r23 < 2l c. 6. Discussions and numerical simulations In order to make the readers understand our results better, we will perform some numerical simulations to illustrate our theoretical results. The numerical simulations are given by the following Milstein scheme [21]. Consider the discretization of model (1.4) for t ¼ 0; Mt; 2Mt; ; nMt: Group(a)
S(t)
Relative frequency density
2
100
1
50
0
0
100
200
300
400
0 −1
t I(t)
Group(b)
100
0.5
50 0
100
200
300
400
0 −1
t R(t)
Group(c)
100
1
50 0
100
200
t
0 1 I(t) at time 150
2
Relative frequency density
2
0
2
Relative frequency density
1
0
0 1 S(t) at time 150
300
400
0 −1
0 1 R(t) at time 150
Fig. 1. The solution of stochastic system (1.4) and its histogram with Sð0Þ ¼ 0:62; Ið0Þ ¼ 0:54; Rð0Þ ¼ 0:54; K ¼ 0:4; b ¼ 0:8; r1 ¼ 0:1; r2 ¼ 0:1; r3 ¼ 0:1; r4 ¼ 0:1.
2
l ¼ 0:2; e ¼ 0:1; c ¼ 0:2;
129
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
Group(a)
S(t)
Relative frequency density
2
4
1
2
0
0
100
200
300
0 −1
400
4
1
2 0
100
200
300
0 −1
400
0 1 I(t) at time 150
t R(t)
Group(c)
2
Relative frequency density
2
0
0 1 S(t) at time 150
t I(t)
Group(b)
2
Relative frequency density
4
5
2 0
0
100
200
300
0 −1
400
0 1 R(t) at time 150
t Fig. 2. The solution of stochastic system (1.4) and e ¼ 0:1; c ¼ 0:2; r1 ¼ 0:2; r2 ¼ 0:2; r3 ¼ 0:2; r4 ¼ 0:2.
histogram
with
Sð0Þ ¼ 0:62; Ið0Þ ¼ 0:54; Rð0Þ ¼ 0:54; K ¼ 0:4; b ¼ 0:8;
pffiffiffiffiffiffiffi 1 1 2 2
pffiffiffiffiffiffiffi 1 2 2 pffiffiffiffiffiffi 1 2 2 Ikþ1 ¼ Ik þ ½bSk Ik ðc þ l þ ÞIk Mt þ Ik r2 n2;k Mt þ r2 ðn21;k 1ÞMt þ Sk Ik r4 n4;k 4t þ r4 ðn4;k 1ÞMt ; 2 2
pffiffiffiffiffiffi 1 2 2 Rkþ1 ¼ Rk þ ðcIk lRk ÞMt þ Rk r3 n3;k Mt þ r3 ðn3;k 1ÞMt ; 2
Skþ1 ¼ Sk þ ðK bSk Ik lSk ÞMt þ Sk
its
2
l ¼ 0:2;
pffiffiffiffiffiffi
r1 n1;k Mt þ r21 ðn21;k 1ÞMt þ Sk Ik r4 n4;k 4t þ r24 ðn24;k 1ÞMt ;
ð6:1Þ
where Mt is time increment n1k ; n2k ; n3k and n4k ðk ¼ 1; 2; . . . ; nÞ are Nð0; 1Þ-distributed independent random variables. In the figures, the blue lines and the red lines represent solutions of deterministic system (1.1) and stochastic system (1.4) respectively. In Figs. 1–3, we choose the parameters K ¼ 0:4; b ¼ 0:8; l ¼ 0:2; e ¼ 0:1; c ¼ 0:2 and the initial point Sð0Þ ¼ 0:62; Ið0Þ ¼ 0:54; Rð0Þ ¼ 0:54, for numerical simulations of stochastic system (1.4). We use different values of r1 ; r2 ; r3 and r4 in order to understand their role on the dynamics.
3
1.4
1.2
2.5
1 2
I(t)
S(t)
0.8 1.5
0.6 1 0.4 0.5
0
0.2
0
50
100
t
150
0
0
50
100
150
t
Fig. 3. Deterministic and stochastic trajectories of system for initial condition and parameter values: Sð0Þ ¼ 0:62; Ið0Þ ¼ 0:54; Rð0Þ ¼ 0:54; K ¼ 0:4; b ¼ 0:8; l ¼ 0:2; e ¼ 0:1; c ¼ 0:2; r1 ¼ 0:1; r2 ¼ 0:7; r3 ¼ 0:1; r4 ¼ 0:6.
130
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
Group(a) 1.5 0.2
I(t)
R(t)
0.2
S(t)
1 0.5
0
50
100
0
150
0.1
0
50
t
Group(b) 1.5
100
0.3
0.3
0.2
0.2
0.1
50
50
100
0
150
100
150
100
150
t
I(t)
S(t) 0.5
0
0
t
1
0
0
150
R(t)
0
0.1
0.1
0
50
t
100
0
150
0
50
t
t
Fig. 4. Deterministic and stochastic trajectories of system for initial condition and parameter values: Sð0Þ ¼ 0:7; Ið0Þ ¼ 0:2; Rð0Þ ¼ 0:1; K ¼ 0:2; b ¼ 0:4; l ¼ 0:2; e ¼ 0:1; c ¼ 0:2; Group (a): r1 ¼ 0:1; r2 ¼ 0:1; r3 ¼ 0:1, r4 ¼ 0:1; Group (b): r1 ¼ 0:05; r2 ¼ 0:05; r3 ¼ 0:05; r4 ¼ 0:05.
For a clear understanding of this idea we can look at the stationary distribution for parameter values as mentioned in Section 3 and for two different sets of forcing intensities r1 ¼ 0:2; r2 ¼ 0:2; r3 ¼ 0:2; r4 ¼ 0:2 and r1 ¼ 0:1; r2 ¼ 0:1; r3 ¼ 0:1; r4 ¼ 0:1. Parameter values chosen above are consistent with the conditions required for the existence of stationary distribution (see Theorem 3.1). The distributions of the population obtained through numerical simulation for two different sets of forcing intensities are reported in Figs. 1 and 2. These figures report the results of one simulation run and stationary distribution about the populations of susceptible individuals, infected individuals and removed individuals. After some initial transients the population densities fluctuate around the deterministic steady-state values S ¼ 0:625; I ¼ 0:55 and R ¼ 0:55. From the Figs. 1 and 2, with intensities increasing, the amplitude of fluctuation is larger. In Fig. 3, we choose r2 ¼ 0:7; r4 ¼ 0:6, other parameters take the same values as Fig. 1. As a result, the population of infected individuals goes to extinction. And, we can verify that the conditions in Theorem 4 are satisfied. That is to say, large noises can lead the disease to extinction, which is a phenomenon different from its corresponding deterministic system (1.1). Group(a) 1.5 0.2
I(t)
R(t)
0.2
S(t)
1 0.5
0
50
100
0
150
0.1
0
50
t
Group(b) 1.5
100
0.3
0.3
0.2
0.2
0.1
50
100
t
50
150
100
150
100
150
t
I(t)
S(t) 0.5
0
0
t
1
0
0
150
R(t)
0
0.1
0
0.1
0
50
100
t
150
0
0
50
t
Fig. 5. Deterministic and stochastic trajectories of system for initial condition and parameter values: Sð0Þ ¼ 0:7; Ið0Þ ¼ 0:2; Rð0Þ ¼ 0:1; K ¼ 0:2; b ¼ 0:4; l ¼ 0:2; e ¼ 0:1; c ¼ 0:2; Group (a): r1 ¼ 0:1; r2 ¼ 0:2; r3 ¼ 0:3; r4 ¼ 0:2; Group (b): r1 ¼ 0; r2 ¼ 0:1; r3 ¼ 0:1; r4 ¼ 0:1.
Y. Zhou et al. / Applied Mathematics and Computation 244 (2014) 118–131
131
In the following, we consider the long time behavior of system (1.4) in the case of R0 < 1. In Figs. 4 and 5, we always choose initial value Sð0Þ ¼ 0:7; Ið0Þ ¼ 0:2, Rð0Þ ¼ 0:1 and parameters K ¼ 0:2; b ¼ 0:4; l ¼ 0:2; e ¼ 0:1; c ¼ 0:2 with different intensities of white noise which satisfy the conditions in Theorem 5.1. In Fig. 4(a), we choose r1 ¼ 0:1; r2 ¼ 0:1; r3 ¼ 0:1; r4 ¼ 0:1 and In Fig. 4(b), we choose r1 ¼ 0:05; r2 ¼ 0:05; r3 ¼ 0:05; r4 ¼ 0:05. As we can see, the curves of system (1.4) always oscillate around the curves of system (1.1) in time, and moreover, the larger the intensities of the white noises are, the larger the fluctuations of the solutions will be. Furthermore, if we take r1 the same value as in Fig. 5(a), but r2 ; r3 and r4 take larger values r2 ¼ 0:2; r3 ¼ 0:3; r4 ¼ 0:2. We can see from Fig. 5(a) that the fluctuation of the solution will be very smaller. That is to say, the intensities of r2 ; r3 and r4 have little effects on the solution. When r1 ¼ 0, then E0 becomes the disease-free equilibrium of system (1.4). We can see from Fig. 5(b) that solution of system (1.4) is stochastically asymptotically stable in the large, which supports the conclusion of Corollary 5.1. Acknowledgments The authors thank the editors and reviewers for their valuable comments. This research is supported by National Natural Science Foundation of China (11271260, 11071164) and the Innovation Program of Shanghai Municipal Education Commission (13ZZ116). Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/ j.amc.2014.06.100. References [1] W.O. Kermack, A.G. McKendrick, Contributions to the mathematical theory of epidemics (Part I), Proc. Soc. London Ser. A 115 (1927) 700–721. [2] J.M. Tchuenche, A. Nwagwo, R. Levins, Global behaviour of an SIR epidemic model with time delay, Math. Methods Appl. Sci. 30 (2007) 733–749. [3] T.L. Zhang, Z.D. Teng, Permanence and extinction for a nonautonomous SIRS epidemic model with time delay, Appl. Math. Model. 33 (2009) 1058– 1071. [4] F.P. Zhang, Z.Z. Li, F. Zhang, Global stability of an SIR epidemic model with constant infectious period, Appl. Math. Comput. 199 (2008) 285–291. [5] A. Korobeinikov, Lyapunov functions and global stability for SIR and SIRS epidemiological models with non-linear transmission, Bull. Math. Biol. 68 (2006) 615–626. [6] R.M. Anderson, R.M. May, Population biology of infectious diseases. Part I, Nature 280 (1979) 361–367. [7] R.M. May, Stability and Complexity in Model Ecosystems, Princeton University Press, New Jersey, 1973. [8] C.Y. Ji, D.Q. Jiang, N.Z. Shi, The behavior of an SIR epidemic model with stochastic perturbation, Stoch. Anal. Appl. 30 (2012) 755–773. [9] A. Gray, D. Greenhalgh, L. Hu, X.R. Mao, J.F. Pan, A stochastic differential equation SIS epidemic model, SIAM J. Appl. Math. 71 (2011) 876–902. [10] D.Q. Jiang, J.J. Yu, C.Y. Ji, N.Z. Shi, Asymptotic behavior of global positive solution to a stochastic SIR model, Math. Comput. Model. 54 (2011) 221–232. [11] J.R. Artalejo, A. Economou, M.J. Lopez-Herrero, On the number of recovered individuals in the SIS and SIR stochastic epidemic models, Math. Biosci. 228 (2010) 45–55. [12] Q. Lu, Stability of SIRS system with random perturbations, Physica A 388 (2009) 3677–3686. [13] Y.G. Lin, D.Q. Jiang, Long-time behavior of perturbed SIR model by white noise, Discrete Contin. Dyn. Syst. Ser. B 18 (2013) 1873–1887. [14] X.R. Mao, Stochastic Differential Equations and Applications, 2nd Edition., Horwood, Chichester, UK, 2007. [15] C.Y. Ji, D.Q. Jiang, Dynamics of a stochastic density dependent predator-prey system with Beddington–DeAngelis functional response, J. Math. Anal. Appl 381 (2011) 441–453. [16] R.Z. Hasminskii, Stochastic Stability of Differential Equations, Sijthoff and Noordhoff, 1980. [17] T.C. Gard, Introduction to Stochastic Differential Equations, Marcel Dekker, New York, 1987. [18] G. Strang, Linear Algebra and Its Applications, Thomson Learning Inc, 1988. [19] C. Zhu, G. Yin, Asymptotic properties of hybrid diffusion systems, SIAM J. Control Optim. 46 (2007) 1155–1179. [20] X.R. Mao, Almost sure asymptotic bounds for a class of stochastic differential equations, Stoch. Stoch. Rep. 41 (1992) 57–69. [21] D.J. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Rev. 43 (2001) 525–546.