# Dynamics Analysis of a Stochastic SIR Epidemic Model

Authors: Feng Rao
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We investigate an SIR epidemic model with stochastic perturbations. We assume that stochastic perturbations are of a white noise type which is directly proportional to the distances of three variables from the steady-state values, respectively. By constructing suitable Lyapunov functions and applying Itô’s formula, some qualitative properties are obtained, such as the existence of global positive solutions, stochastic boundedness, and permanence. A series of numerical simulations to illustrate these mathematical findings are presented.


---

## Body

## 1. Introduction

Almost all mathematical models for the transmission of infectious diseases descend from the classical susceptible-infective-removed (SIR) model of Kermack and McKendrick [1]. The dynamic behavior of different epidemic models and a lot of their extensions is well investigated by a number of scholars; see [2–11]. The basic and important research subjects for recent studies are the existence of the threshold values which distinguish whether the disease dies out, the stability of the disease-free and the endemic equilibria, permanence, and extinction [12]. During the last few decades, a number of realistic transmission functions have become the focus of considerable attention, and many authors are interested in the formulation of nonlinear incidence rate (see [13–17]). A nonlinear incidence rate can arise from saturation effects that if the proportion of the infection in a population is very high, so that exposure to the disease agent is virtually certain, then the transmission rate may respond more slowly than linear to the increase in the number of infection [18]. For example, Capasso and Serio [19] introduced a saturated transmission rate f(S,I)=kSI/(1+αI), where kI measures the infection force of the disease and 1/(1+αI) measures the inhibition effect from the behavioral change of the susceptible individuals when their number increases or from the crowding effect of the infective individuals. To be biologically feasible, the function of the incidence rate f(S,I)=kSI/(1+αI) is a positive continuous and differentiable and satisfies the conditions
(1)f(S,0)=0=f(0,I),∂f(S,I)∂S=kI1+αI>0,∂f(S,I)∂I=kS(1+αI)2>0,
for all S,I>0. It is easy to know that the function f(S,I) is concave with respect to the variable I; that is,
(2)∂2f(S,I)∂I2=-2αkS(1+αI)3<0forS,I>0,
which implies that when the number of infections is very high that the exposure to the disease agent is virtually certain, the incidence rate will respond more slowly than linearly to the disease in I.In the real world, population dynamics is inevitably subjected to environmental noise, which is an important component in an ecosystem. Most natural phenomena do not follow strictly deterministic laws but rather oscillate randomly about some average values, so that the population density never attains a fixed value with the advancement of time [20, 21]. Recent advances in stochastic differential equations enable a lot of authors to introduce randomness into deterministic model of physical phenomena to reveal the effect of environmental variability, whether it is a random noise in the system of differential equations or environmental fluctuations in parameters; see [12, 13, 22–30]. Of them, Tuckwell and Williams [28] investigated the properties of a simple discrete time stochastic epidemic model. A classical model of an SIRS epidemic in an open population was considered by El Maroufy et al. [12]. They established the global stability of disease-free and endemic equilibrium points for both the deterministic and stochastic models. Based on the theory of stochastic differential equation, Cai et al. [13] studied the dynamics of an SIRS epidemic model with a ratio-dependent incidence rate. In [29], the authors extended the classical SIRS epidemic model incorporating media coverage from a deterministic framework to a stochastic differential equation and focused on how environmental fluctuations of the contact coefficient affect the extinction of the disease.To the best of our knowledge, a small amount of work has been done with stochastic perturbation on an SIR epidemic model with a saturated transmission ratekSI/(1+αI). The purpose of this paper is to study that the stochastic factor has a significant effect on the dynamics of SIR epidemic model with a saturated incidence rate. The organization of this paper is as follows. In the next section, we present the formulation of mathematical model with environmental noise. We give some properties about deterministic model (4) and carry out the analysis of the dynamical properties of stochastic model (3), respectively. Finally, we give a concluding section.


## 2. Model and Dynamics Analysis

LetS(t) be the number of susceptible individuals, I(t) the number of infective individuals, and R(t) the number of removed individuals at time t, respectively. Motivated by [31], we assume that stochastic perturbations are of white noise type, which are directly proportional to distances S(t),I(t),R(t) from the steady-state values of S*,I*,R*and influence on dS(t)/dt,dI(t)/dt, dR(t)/dt, respectively. In this way, an SIR epidemic model with a saturated transmission rate and stochastic fluctuations will be reduced to the following form:
(3)dS=(b-dS-kSI1+αI+γR)dt+σ1(S-S*)dB(t),dI=(kSI1+αI-(d+μ)I)dt+σ2(I-I*)dB(t),dR=(μI-(d+γ)R)dt+σ3(R-R*)dB(t).
All parameters are positive constants, b is the recruitment rate of the population, d is the natural death rate of the population, k is the proportionality constant, α is the parameter that measures the psychological or inhibitory effect, γ is the rate at which recovered individuals lose immunity and return to the susceptible class, and μ is the natural recovery rate of the infective individuals. Note that σ1, σ2, and σ3 are real constants and known as the intensity of the stochastic environment and B(t) is standard Brownian motion.


### 2.1. Dynamics of the Deterministic Model

In this subsection, whenσ1=σ2=σ3=0, we consider the deterministic SIR epidemic model:
(4)dSdt=b-dS-kSI1+αI+γR,dIdt=kSI1+αI-(d+μ)I,dRdt=μI-(d+γ)R.Because of the biological meaning of the components(S(t),I(t),R(t)), we focus on the model in the first quadrant ℝ+3={(S,I,R)∈ℝ3:S≥0,I≥0,R≥0}. Model (4) always has a disease-free equilibrium E0=(b/d,0,0), which corresponds to the extinction of the disease.Define the basic reproduction number as(5)R0=bkd(d+μ),
which denotes the number of individuals infected by a single infected individual placed in a totally susceptible population.Theorem 1.
From model (4), it follows that(i)
ifR0≤1, there is no positive equilibrium;(ii)
ifR0>1, there is a unique endemic equilibrium E*=(S*,I*,R*), which corresponds to the coexistence of S, I, and R and is given by
(6)S*=(d+μ)(d(d+μ+γ)+αb(d+γ))αd(d+μ)(d+γ)+dk(d+μ+γ),I*=(d+γ)(bk-d(d+μ))(d+μ)(d(d+μ+γ)+αb(d+γ))S*,R*=μd+γI*.In other words, when R0>1, the disease can invade a totally susceptible population and the number of cases will increase, whereas when R0≤1, the disease will always fail to spread.Lemma 2.
The planeS+I+R=b/d is a manifold of model (4), which is attracting in the first octant.Proof.
Summing up the three equations in (4) and denoting N(t)=S(t)+I(t)+R(t), we obtain
(7)dNdt=b-dN.
It is clear thatN(t)=b/d is a solution of (7) and for any N(t0)≥0, the general solution of (7) is
(8)N(t)=1d(b-(b-dN(t0))e-d(t-t0)).
Hence,
(9)lim⁡t→∞N(t)=bd,
which implies the conclusion.Theorem 3.
The endemic equilibrium pointE*=(S*,I*,R*) is globally asymptotically stable in ℝ+3.Proof.
The Jacobian matrixJ(S,I,R) at equilibrium point E*=(S*,I*,R*) is given by
(10)J(E*)=(J11J12J13J21J22J23J31J32J33),
where
(11)J11=-d-kI*1+αI*,J12=-kS*(1+αI*)2,J13=γ,J21=kI*1+αI*,J22=kS*(1+αI*)2-d-μ,J23=0,J31=0,J32=μ,J33=-d-γ.
The characteristic equation at the interior equilibrium pointE* is
(12)λ3+Q1λ2+Q2λ+Q3=0,
where
(13)Q1=-(J11+J22+J33),Q2=J11J22-J12J21+(J11+J22)J33,Q3=J12J21J33-J21J13J32-J11J22J33.
It is clear that
(14)J11<0,J12<0,J13>0,J21>0,J22<0,J32>0,J33<0.
Here Q1>0,  Q2>0 and Q3>0.
NowQ1Q2-Q3=(J11+J22)J12J21+J21J13J32-J112(J22+J33)-J222(J11+J33)-J332(J11+J22)-2J11J22J33>0. Therefore, model (4) is globally stable at the equilibrium E*=(S*,I*,R*).


### 2.2. Dynamics of the Stochastic Model

Let(Ω,ℱ,P) be a complete probability space with a filtration {ℱt}t∈ℝ+ satisfying the usual conditions; that is, it is right continuous and increasing while ℱ0 contains all P-null sets. Denote
(15)X(t)=(S(t),I(t),R(t))≜(x1(t),x2(t),x3(t)),
and the norm |X(t)|=S2(t)+I2(t)+R2(t). And denote C2,1(ℝ3×(0,∞);ℝ+) as the family of all nonnegative functions V(X,t) defined on ℝ3×(0,∞) such that they are continuously twice differentiable in X and once in t.We define the differential operatorL associated with three-dimensional stochastic differential equation
(16)dX(t)=f(X,t)dt+g(X,t)dB(t),
as
(17)L=∂∂t+∑i=13‍fi(X,t)∂∂xi+12∑i,j=13‍(gT(X,t)g(X,t))ij∂2∂xi∂xj,
where
(18)f=(b-dS-kSI1+αI+γRkSI1+αI-(d+μ)IμI-(d+γ)R),g=diag⁡(σ1(S-S*),σ2(I-I*),σ3(R-R*)).IfL acts on a function V∈C2,1(ℝ3×(0,∞);ℝ+), then we denote
(19)LV(X,t)=Vt(X,t)+VX(X,t)f(X,t)+12trace(gT(X,t)VXX(X,t)g(X,t)),
where T means transposition.In this subsection, we first show the existence of a unique positive global solution of the stochastic model (3).Theorem 4.
For model (3) and any given initial value (S(0),I(0),R(0))∈ℝ+3, there is a unique solution (S(t),I(t),R(t)) on t≥0 and will remain in ℝ+3 with probability one.Proof.
Since the coefficients of model (3) satisfy the local Lipschitz condition, there is a unique local solution on [0,τe), where τe is the explosion time. Therefore, by Itô’s formula, the unique local solution of model (3) is positive. Next, let us show that this solution is global; that is, τe=∞ a.s.
Letn0>0 be sufficiently large for S(0), I(0), and R(0) lying with the interval [1/n0,n0]. For each integer n≥n0, define a sequence of stopping times by
(20)τn=inf⁡{t∈[0,τe]:S(t)∉(1n,n)orI(t)∉(1n,n)orR(t)∉(1n,n)},
where we set inf⁡∅=∞ (∅ represents the empty set) in this paper. Since τn is nondecreasing as n→∞, there exists the limit
(21)τ∞=lim⁡n→∞τn.
Then τ∞≤τe a.s. Now, we need to show τ∞=∞ a.s. If this statement is violated, then there exist T>0 and ɛ∈(0,1) such that
(22)P{τ∞≤T}>ɛ.
Thus, there is an integer n1≥n0 such that
(23)P{τn≤T}≥ɛ∀n≥n1.
Define aC3-function V:ℝ+3→ℝ+ by
(24)V(S,I,R)=(S-1-ln⁡S)+(I-1-ln⁡I)+(R-1-ln⁡R),
which is a nonnegative function. If (S(t),I(t),R(t))∈ℝ+3, by using Itô’s formula, we compute
(25)dV(S,I,R)=((1-R*R)2(1-1S)(b-dS-kSI1+αI+γR)+(1-1I)(kSI1+αI-(d+μ)I)+(1-1R)(μI-(d+γ)R)+σ122(1-S*S)2+σ222(1-I*I)2+σ322(1-R*R)2)dt+(σ1(1-1S)(S-S*)+σ2(1-1I)(I-I*)+σ3(1-1R)(R-R*))dB(t)=(b+3d+μ+γ+σ122(1-S*S)2+σ222(1-I*I)2+σ322(1-R*R)2+kI1+αI-kS1+αI-d(S+I+R)-bS-γRS-μIR(1-S*S)2)dt+(σ1(1-1S)(S-S*)+σ2(1-1I)(I-I*)+σ3(1-1R)(R-R*))dB(t)≤(b+3d+μ+γ+σ122(1-S*S)2+σ222(1-I*I)2+σ322(1-R*R)2+kα)dt+(σ1(1-1S)(S-S*)+σ2(1-1I)(I-I*)+σ3(1-1R)(R-R*))dB(t)≤Mdt+(σ1(1-1S)(S-S*)+σ2(1-1I)(I-I*)+σ3(1-1R)(R-R*))dB(t),
where M is a positive constant. Integrating both sides of the above inequality from 0 to τn∧T, we get
(26)∫0τn∧T‍dV(S(s),I(s),R(s))≤∫0τn∧T‍Mds+∫0τn∧T‍(σ1(1-1S)(S-S*)+σ2(1-1I)(I-I*)+σ3(1-1R)(R-R*))dB(s),
where τn∧T=min⁡{τn,T}. Then taking the expectations leads to
(27)EV(S(τn∧T),I(τn∧T),R(τn∧T))≤V(S(0),I(0),R(0))+MT.
SetΩn={τn≤T} for n≥n1 and from (23), we have P(Ωn)≥ɛ. For every ν∈Ωn, there are some i such that xi(τn,ν) equals either n or 1/n for i=1,2,3; hence V(S(τn,ν),I(τn,ν),R(τn,ν)) is no less than min⁡{n-1-ln⁡n,1/n-1-ln⁡(1/n)}. Then we obtain
(28)V(S(0),I(0),R(0))+MT≥E(1Ωn(ν)V(S(τn),I(τn),R(τn)))≥ɛmin⁡{n-1-ln⁡n,1n-1-ln⁡1n},
where 1Ωn(ν) is the indicator function of Ωn. Letting n→∞ leads to the contradiction ∞=V(S(0),I(0),R(0))+MT<∞. This completes the proof.Theorem4 shows that the solution to model (3) will remain in ℝ+3. The property makes us continue to discuss how the solution varies in ℝ+3 in more detail. Here, we present that the definition of stochastic ultimate boundedness [32] is one of the important topics in population dynamics and is defined as follows.Definition 5.
The solutionsX(t)=(S(t),I(t),R(t)) of model (3) are said to be stochastically ultimately bounded, if for any ɛ∈(0,1), there is a positive constant δ=δ(ɛ), such that for any initial value (S(0),I(0),R(0))∈ℝ+3, the solution X(t) to model (3) has the property that
(29)limsup⁡t→∞P{|X(t)|>δ}<ɛ.Theorem 6.
The solutions of model (3) are stochastically ultimately bounded for any initial value (S(0),I(0),R(0))∈ℝ+3.Proof.
From Theorem4, the solution will remain in ℝ+3 for all t≥0 almost surely. Define a function
(30)V(S,I,R)=et(Sθ+Iθ+Rθ),
for (S,I,R)∈ℝ+3 and θ>1. By Itô’s formula we obtain
(31)dV(S,I,R)=et(+σ32Rθ(1-R*R)2)θSθ-1(b-dS-kSI1+αI+γR)+θIθ-1(kSI1+αI-(d+μ)I)+θRθ-1(μI-(d+γ)R)+θ(θ-1)2(σ12Sθ(1-S*S)2+σ22Iθ(1-I*I)2+σ32Rθ(1-R*R)2))dt+etθ(σ1Sθ(1-S*S)+σ2Iθ(1-I*I)+σ3Rθ(1-R*R))dB(t)≤Cetdt+etθ(σ1Sθ(1-S*S)+σ2Iθ(1-I*I)+σ3Rθ(1-R*R))dB(t),
where C>0 is a suitable constant.
Based on Theorem4 and from (31), we have
(32)E(et∧τnV(S(t∧τn),I(t∧τn),R(t∧τn)))≤V(S(0),I(0),R(0))+CE∫0t∧τn‍esds.
Letting n→∞ yields
(33)etEV(S(t),I(t),R(t))≤V(S(0),I(0),R(0))+C(et-1),
which implies
(34)EV(S(t),I(t),R(t))≤e-tV(S(0),I(0),R(0))+C.
Note that(35)|X(t)|θ=(S2(t)+I2(t)+R2(t))θ/2≤3θ/2max⁡⁡{Sθ(t),Iθ(t),Rθ(t)}≤3θ/2(Sθ+Iθ+Rθ).
Then we get
(36)E|X(t)|θ≤3θ/2(e-tV(S(0),I(0),R(0))+C),
which means
(37)limsup⁡t→∞E|X(t)|θ≤3θ/2C<∞.
Therefore, there exists a positive constantδ1 such that
(38)limsup⁡t→∞E|X(t)|<δ1.
For any ɛ>0, set δ=δ12/ɛ2, then by Chebyshev’s inequality,
(39)P{|X(t)>δ|}≤E|X(t)|δ.
Thus, we obtain
(40)limsup⁡t→∞P{|X(t)>δ|}≤δ1δ=ɛ,
which yields the required assertion.Generally speaking, the nonexplosion property, the existence, and the uniqueness of the solution are not enough but the property of permanence is more desirable since it means the long time survival in a population dynamics. Now, the definition of stochastic permanence [33] will be given below.Definition 7.
The solutionsX(t)=(S(t),I(t),R(t)) of model (3) are said to be stochastically permanent, if for any ɛ∈(0,1), there exists a pair of positive constants δ=δ(ɛ) and χ=χ(ɛ) such that for any initial value (S(0),I(0),R(0))∈ℝ+3, the solution X(t) to model (3) has the properties
(41)liminf⁡t→∞P{|X(t)|≤δ}≥1-ɛ,liminf⁡t→∞P{|X(t)|≥χ}≥1-ɛ.Theorem 8.
Assumed<b and for any initial value (S(0),I(0),R(0))∈ℝ+3, the solution (S(t),I(t),R(t)) satisfies
(42)limsup⁡t→∞E(|X(t)|-ϑ)≤Q,
where ϑ is an arbitrary positive constant satisfying
(43)ϑ+12max⁡{σ12,σ22,σ32}<b-d,(44)Q=3ϑ(4ωC1+C2)4ωC1×max⁡{1,(2C1+C2+C22+4C1C22C1)ϑ-2},
in which ω is an arbitrary positive constant satisfying
(45)ω<b-d-ϑ+12max⁡{σ12,σ22,σ32},(46)C1=b-d-ϑ+12max⁡{σ12,σ22,σ32}-ω,C2=d+max⁡{σ12,σ22,σ32}+2ω.Proof.
Define a function(47)V(S,I,R)=1S+I+R,
for (S(t),I(t),R(t))∈ℝ+3; using Itô’s formula, we get
(48)dV(S,I,R)=-V2(b-d(S+I+R))dt+V3(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2)dt-V2(σ1(S-S*)+σ2(I-I*)+σ3(R-R*))dB(t).
Choosing a positive constantϑ that satisfies (43) and applying Itô’s formula, we obtain
(49)L(1+V)ϑ=ϑ(1+V)ϑ-1×((σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2)-V2(b-d(S+I+R))+V3(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2))+ϑ(ϑ-1)2V4(1+V)ϑ-2×(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2)=ϑ(1+V)ϑ-2×(ϑ+12-V2(b-d(S+I+R))-V3(b-d(S+I+R))+V3(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2)+ϑ+12V4(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2))=ϑ(1+V)ϑ-2W,
where
(50)W=-V2(b-d(S+I+R))-V3(b-d(S+I+R))+V3(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2)+ϑ+12V4(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2)≤dV-(b-d)V2+V3(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2)+ϑ+12V4(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2).
Using the facts that(51)V3(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2)<max⁡⁡{σ12,σ22,σ32}V,V4(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2)<max⁡⁡{σ12,σ22,σ32}V2,
then,
(52)W≤(d+max⁡{σ12,σ22,σ32})V-(b-d-ϑ+12max⁡{σ12,σ22,σ32})V2.
Letω>0 be sufficiently small such that it satisfies (45), by Itô’s formula; then
(53)L(eωt(1+V)ϑ)=ωeωt(1+V)ϑ+eωtL(1+V)ϑ=eωt(1+V)ϑ-2(ω(1+V)2+W)≤eωt(1+V)ϑ-2(ω-C1V2+C2V)≤Q0eωt,
where Q0=(4ωC1+C2)/(4C1)max⁡{1,((2C1+C2+C22+4C1C2)/(2C1))ϑ-2} and C1,C2 have been defined in the statement of the theorem. Thus,
(54)E(eωt(1+V)ϑ)≤(1+V(0))ϑ+Q0ωeωt.
Therefore we obtain
(55)limsup⁡t→∞E(V(t)ϑ)≤limsup⁡t→∞E(1+V)ϑ≤Q0ω.
For(S,I,R)∈ℝ+3, we know that (S+I+R)ϑ≤3ϑ(S2+I2+R2)ϑ/2≤3ϑ|X(t)|ϑ; consequently,
(56)limsup⁡t→∞E(1|X(t)|ϑ)≤3ϑlimsup⁡t→∞E(V(t)ϑ)≤3ϑQ0ω=Q,
which completes the proof.Considering Chebyshev inequality, Theorems6 and 8, we immediately obtain the following result.Theorem 9.
Assumemax⁡{σ12,σ22,σ32}<2(b-d); then the solutions of model (3) are stochastically permanent.Proof.
From Theorem6, we have P{|X(t)|>δ}≤ɛ which implies
(57)P{|X(t)|≤δ}≥1-ɛ.
This follows that
(58)liminf⁡t→∞P{|X(t)|≤δ}≥1-ɛ.
By Theorem8, we get
(59)limsup⁡t→∞E(1|X(t)|ϑ)≤Q.
For any ɛ>0, let χ=ɛϑ/Qϑ; then
(60)P{|X(t)|<χ}=P{1|X(t)|>1χ}≤χ1/ϑE(|X(t)|-ϑ).
Hence,
(61)limsup⁡t→∞P{|X(t)|<χ}≤χ1/ϑQ=ɛ,
which follows that
(62)liminf⁡t→∞P{|X(t)|≥χ}≥1-ɛ.
The proof is complete.


## 2.1. Dynamics of the Deterministic Model

In this subsection, whenσ1=σ2=σ3=0, we consider the deterministic SIR epidemic model:
(4)dSdt=b-dS-kSI1+αI+γR,dIdt=kSI1+αI-(d+μ)I,dRdt=μI-(d+γ)R.Because of the biological meaning of the components(S(t),I(t),R(t)), we focus on the model in the first quadrant ℝ+3={(S,I,R)∈ℝ3:S≥0,I≥0,R≥0}. Model (4) always has a disease-free equilibrium E0=(b/d,0,0), which corresponds to the extinction of the disease.Define the basic reproduction number as(5)R0=bkd(d+μ),
which denotes the number of individuals infected by a single infected individual placed in a totally susceptible population.Theorem 1.
From model (4), it follows that(i)
ifR0≤1, there is no positive equilibrium;(ii)
ifR0>1, there is a unique endemic equilibrium E*=(S*,I*,R*), which corresponds to the coexistence of S, I, and R and is given by
(6)S*=(d+μ)(d(d+μ+γ)+αb(d+γ))αd(d+μ)(d+γ)+dk(d+μ+γ),I*=(d+γ)(bk-d(d+μ))(d+μ)(d(d+μ+γ)+αb(d+γ))S*,R*=μd+γI*.In other words, when R0>1, the disease can invade a totally susceptible population and the number of cases will increase, whereas when R0≤1, the disease will always fail to spread.Lemma 2.
The planeS+I+R=b/d is a manifold of model (4), which is attracting in the first octant.Proof.
Summing up the three equations in (4) and denoting N(t)=S(t)+I(t)+R(t), we obtain
(7)dNdt=b-dN.
It is clear thatN(t)=b/d is a solution of (7) and for any N(t0)≥0, the general solution of (7) is
(8)N(t)=1d(b-(b-dN(t0))e-d(t-t0)).
Hence,
(9)lim⁡t→∞N(t)=bd,
which implies the conclusion.Theorem 3.
The endemic equilibrium pointE*=(S*,I*,R*) is globally asymptotically stable in ℝ+3.Proof.
The Jacobian matrixJ(S,I,R) at equilibrium point E*=(S*,I*,R*) is given by
(10)J(E*)=(J11J12J13J21J22J23J31J32J33),
where
(11)J11=-d-kI*1+αI*,J12=-kS*(1+αI*)2,J13=γ,J21=kI*1+αI*,J22=kS*(1+αI*)2-d-μ,J23=0,J31=0,J32=μ,J33=-d-γ.
The characteristic equation at the interior equilibrium pointE* is
(12)λ3+Q1λ2+Q2λ+Q3=0,
where
(13)Q1=-(J11+J22+J33),Q2=J11J22-J12J21+(J11+J22)J33,Q3=J12J21J33-J21J13J32-J11J22J33.
It is clear that
(14)J11<0,J12<0,J13>0,J21>0,J22<0,J32>0,J33<0.
Here Q1>0,  Q2>0 and Q3>0.
NowQ1Q2-Q3=(J11+J22)J12J21+J21J13J32-J112(J22+J33)-J222(J11+J33)-J332(J11+J22)-2J11J22J33>0. Therefore, model (4) is globally stable at the equilibrium E*=(S*,I*,R*).


## 2.2. Dynamics of the Stochastic Model

Let(Ω,ℱ,P) be a complete probability space with a filtration {ℱt}t∈ℝ+ satisfying the usual conditions; that is, it is right continuous and increasing while ℱ0 contains all P-null sets. Denote
(15)X(t)=(S(t),I(t),R(t))≜(x1(t),x2(t),x3(t)),
and the norm |X(t)|=S2(t)+I2(t)+R2(t). And denote C2,1(ℝ3×(0,∞);ℝ+) as the family of all nonnegative functions V(X,t) defined on ℝ3×(0,∞) such that they are continuously twice differentiable in X and once in t.We define the differential operatorL associated with three-dimensional stochastic differential equation
(16)dX(t)=f(X,t)dt+g(X,t)dB(t),
as
(17)L=∂∂t+∑i=13‍fi(X,t)∂∂xi+12∑i,j=13‍(gT(X,t)g(X,t))ij∂2∂xi∂xj,
where
(18)f=(b-dS-kSI1+αI+γRkSI1+αI-(d+μ)IμI-(d+γ)R),g=diag⁡(σ1(S-S*),σ2(I-I*),σ3(R-R*)).IfL acts on a function V∈C2,1(ℝ3×(0,∞);ℝ+), then we denote
(19)LV(X,t)=Vt(X,t)+VX(X,t)f(X,t)+12trace(gT(X,t)VXX(X,t)g(X,t)),
where T means transposition.In this subsection, we first show the existence of a unique positive global solution of the stochastic model (3).Theorem 4.
For model (3) and any given initial value (S(0),I(0),R(0))∈ℝ+3, there is a unique solution (S(t),I(t),R(t)) on t≥0 and will remain in ℝ+3 with probability one.Proof.
Since the coefficients of model (3) satisfy the local Lipschitz condition, there is a unique local solution on [0,τe), where τe is the explosion time. Therefore, by Itô’s formula, the unique local solution of model (3) is positive. Next, let us show that this solution is global; that is, τe=∞ a.s.
Letn0>0 be sufficiently large for S(0), I(0), and R(0) lying with the interval [1/n0,n0]. For each integer n≥n0, define a sequence of stopping times by
(20)τn=inf⁡{t∈[0,τe]:S(t)∉(1n,n)orI(t)∉(1n,n)orR(t)∉(1n,n)},
where we set inf⁡∅=∞ (∅ represents the empty set) in this paper. Since τn is nondecreasing as n→∞, there exists the limit
(21)τ∞=lim⁡n→∞τn.
Then τ∞≤τe a.s. Now, we need to show τ∞=∞ a.s. If this statement is violated, then there exist T>0 and ɛ∈(0,1) such that
(22)P{τ∞≤T}>ɛ.
Thus, there is an integer n1≥n0 such that
(23)P{τn≤T}≥ɛ∀n≥n1.
Define aC3-function V:ℝ+3→ℝ+ by
(24)V(S,I,R)=(S-1-ln⁡S)+(I-1-ln⁡I)+(R-1-ln⁡R),
which is a nonnegative function. If (S(t),I(t),R(t))∈ℝ+3, by using Itô’s formula, we compute
(25)dV(S,I,R)=((1-R*R)2(1-1S)(b-dS-kSI1+αI+γR)+(1-1I)(kSI1+αI-(d+μ)I)+(1-1R)(μI-(d+γ)R)+σ122(1-S*S)2+σ222(1-I*I)2+σ322(1-R*R)2)dt+(σ1(1-1S)(S-S*)+σ2(1-1I)(I-I*)+σ3(1-1R)(R-R*))dB(t)=(b+3d+μ+γ+σ122(1-S*S)2+σ222(1-I*I)2+σ322(1-R*R)2+kI1+αI-kS1+αI-d(S+I+R)-bS-γRS-μIR(1-S*S)2)dt+(σ1(1-1S)(S-S*)+σ2(1-1I)(I-I*)+σ3(1-1R)(R-R*))dB(t)≤(b+3d+μ+γ+σ122(1-S*S)2+σ222(1-I*I)2+σ322(1-R*R)2+kα)dt+(σ1(1-1S)(S-S*)+σ2(1-1I)(I-I*)+σ3(1-1R)(R-R*))dB(t)≤Mdt+(σ1(1-1S)(S-S*)+σ2(1-1I)(I-I*)+σ3(1-1R)(R-R*))dB(t),
where M is a positive constant. Integrating both sides of the above inequality from 0 to τn∧T, we get
(26)∫0τn∧T‍dV(S(s),I(s),R(s))≤∫0τn∧T‍Mds+∫0τn∧T‍(σ1(1-1S)(S-S*)+σ2(1-1I)(I-I*)+σ3(1-1R)(R-R*))dB(s),
where τn∧T=min⁡{τn,T}. Then taking the expectations leads to
(27)EV(S(τn∧T),I(τn∧T),R(τn∧T))≤V(S(0),I(0),R(0))+MT.
SetΩn={τn≤T} for n≥n1 and from (23), we have P(Ωn)≥ɛ. For every ν∈Ωn, there are some i such that xi(τn,ν) equals either n or 1/n for i=1,2,3; hence V(S(τn,ν),I(τn,ν),R(τn,ν)) is no less than min⁡{n-1-ln⁡n,1/n-1-ln⁡(1/n)}. Then we obtain
(28)V(S(0),I(0),R(0))+MT≥E(1Ωn(ν)V(S(τn),I(τn),R(τn)))≥ɛmin⁡{n-1-ln⁡n,1n-1-ln⁡1n},
where 1Ωn(ν) is the indicator function of Ωn. Letting n→∞ leads to the contradiction ∞=V(S(0),I(0),R(0))+MT<∞. This completes the proof.Theorem4 shows that the solution to model (3) will remain in ℝ+3. The property makes us continue to discuss how the solution varies in ℝ+3 in more detail. Here, we present that the definition of stochastic ultimate boundedness [32] is one of the important topics in population dynamics and is defined as follows.Definition 5.
The solutionsX(t)=(S(t),I(t),R(t)) of model (3) are said to be stochastically ultimately bounded, if for any ɛ∈(0,1), there is a positive constant δ=δ(ɛ), such that for any initial value (S(0),I(0),R(0))∈ℝ+3, the solution X(t) to model (3) has the property that
(29)limsup⁡t→∞P{|X(t)|>δ}<ɛ.Theorem 6.
The solutions of model (3) are stochastically ultimately bounded for any initial value (S(0),I(0),R(0))∈ℝ+3.Proof.
From Theorem4, the solution will remain in ℝ+3 for all t≥0 almost surely. Define a function
(30)V(S,I,R)=et(Sθ+Iθ+Rθ),
for (S,I,R)∈ℝ+3 and θ>1. By Itô’s formula we obtain
(31)dV(S,I,R)=et(+σ32Rθ(1-R*R)2)θSθ-1(b-dS-kSI1+αI+γR)+θIθ-1(kSI1+αI-(d+μ)I)+θRθ-1(μI-(d+γ)R)+θ(θ-1)2(σ12Sθ(1-S*S)2+σ22Iθ(1-I*I)2+σ32Rθ(1-R*R)2))dt+etθ(σ1Sθ(1-S*S)+σ2Iθ(1-I*I)+σ3Rθ(1-R*R))dB(t)≤Cetdt+etθ(σ1Sθ(1-S*S)+σ2Iθ(1-I*I)+σ3Rθ(1-R*R))dB(t),
where C>0 is a suitable constant.
Based on Theorem4 and from (31), we have
(32)E(et∧τnV(S(t∧τn),I(t∧τn),R(t∧τn)))≤V(S(0),I(0),R(0))+CE∫0t∧τn‍esds.
Letting n→∞ yields
(33)etEV(S(t),I(t),R(t))≤V(S(0),I(0),R(0))+C(et-1),
which implies
(34)EV(S(t),I(t),R(t))≤e-tV(S(0),I(0),R(0))+C.
Note that(35)|X(t)|θ=(S2(t)+I2(t)+R2(t))θ/2≤3θ/2max⁡⁡{Sθ(t),Iθ(t),Rθ(t)}≤3θ/2(Sθ+Iθ+Rθ).
Then we get
(36)E|X(t)|θ≤3θ/2(e-tV(S(0),I(0),R(0))+C),
which means
(37)limsup⁡t→∞E|X(t)|θ≤3θ/2C<∞.
Therefore, there exists a positive constantδ1 such that
(38)limsup⁡t→∞E|X(t)|<δ1.
For any ɛ>0, set δ=δ12/ɛ2, then by Chebyshev’s inequality,
(39)P{|X(t)>δ|}≤E|X(t)|δ.
Thus, we obtain
(40)limsup⁡t→∞P{|X(t)>δ|}≤δ1δ=ɛ,
which yields the required assertion.Generally speaking, the nonexplosion property, the existence, and the uniqueness of the solution are not enough but the property of permanence is more desirable since it means the long time survival in a population dynamics. Now, the definition of stochastic permanence [33] will be given below.Definition 7.
The solutionsX(t)=(S(t),I(t),R(t)) of model (3) are said to be stochastically permanent, if for any ɛ∈(0,1), there exists a pair of positive constants δ=δ(ɛ) and χ=χ(ɛ) such that for any initial value (S(0),I(0),R(0))∈ℝ+3, the solution X(t) to model (3) has the properties
(41)liminf⁡t→∞P{|X(t)|≤δ}≥1-ɛ,liminf⁡t→∞P{|X(t)|≥χ}≥1-ɛ.Theorem 8.
Assumed<b and for any initial value (S(0),I(0),R(0))∈ℝ+3, the solution (S(t),I(t),R(t)) satisfies
(42)limsup⁡t→∞E(|X(t)|-ϑ)≤Q,
where ϑ is an arbitrary positive constant satisfying
(43)ϑ+12max⁡{σ12,σ22,σ32}<b-d,(44)Q=3ϑ(4ωC1+C2)4ωC1×max⁡{1,(2C1+C2+C22+4C1C22C1)ϑ-2},
in which ω is an arbitrary positive constant satisfying
(45)ω<b-d-ϑ+12max⁡{σ12,σ22,σ32},(46)C1=b-d-ϑ+12max⁡{σ12,σ22,σ32}-ω,C2=d+max⁡{σ12,σ22,σ32}+2ω.Proof.
Define a function(47)V(S,I,R)=1S+I+R,
for (S(t),I(t),R(t))∈ℝ+3; using Itô’s formula, we get
(48)dV(S,I,R)=-V2(b-d(S+I+R))dt+V3(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2)dt-V2(σ1(S-S*)+σ2(I-I*)+σ3(R-R*))dB(t).
Choosing a positive constantϑ that satisfies (43) and applying Itô’s formula, we obtain
(49)L(1+V)ϑ=ϑ(1+V)ϑ-1×((σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2)-V2(b-d(S+I+R))+V3(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2))+ϑ(ϑ-1)2V4(1+V)ϑ-2×(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2)=ϑ(1+V)ϑ-2×(ϑ+12-V2(b-d(S+I+R))-V3(b-d(S+I+R))+V3(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2)+ϑ+12V4(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2))=ϑ(1+V)ϑ-2W,
where
(50)W=-V2(b-d(S+I+R))-V3(b-d(S+I+R))+V3(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2)+ϑ+12V4(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2)≤dV-(b-d)V2+V3(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2)+ϑ+12V4(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2).
Using the facts that(51)V3(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2)<max⁡⁡{σ12,σ22,σ32}V,V4(σ12(S-S*)2+σ22(I-I*)2+σ32(R-R*)2)<max⁡⁡{σ12,σ22,σ32}V2,
then,
(52)W≤(d+max⁡{σ12,σ22,σ32})V-(b-d-ϑ+12max⁡{σ12,σ22,σ32})V2.
Letω>0 be sufficiently small such that it satisfies (45), by Itô’s formula; then
(53)L(eωt(1+V)ϑ)=ωeωt(1+V)ϑ+eωtL(1+V)ϑ=eωt(1+V)ϑ-2(ω(1+V)2+W)≤eωt(1+V)ϑ-2(ω-C1V2+C2V)≤Q0eωt,
where Q0=(4ωC1+C2)/(4C1)max⁡{1,((2C1+C2+C22+4C1C2)/(2C1))ϑ-2} and C1,C2 have been defined in the statement of the theorem. Thus,
(54)E(eωt(1+V)ϑ)≤(1+V(0))ϑ+Q0ωeωt.
Therefore we obtain
(55)limsup⁡t→∞E(V(t)ϑ)≤limsup⁡t→∞E(1+V)ϑ≤Q0ω.
For(S,I,R)∈ℝ+3, we know that (S+I+R)ϑ≤3ϑ(S2+I2+R2)ϑ/2≤3ϑ|X(t)|ϑ; consequently,
(56)limsup⁡t→∞E(1|X(t)|ϑ)≤3ϑlimsup⁡t→∞E(V(t)ϑ)≤3ϑQ0ω=Q,
which completes the proof.Considering Chebyshev inequality, Theorems6 and 8, we immediately obtain the following result.Theorem 9.
Assumemax⁡{σ12,σ22,σ32}<2(b-d); then the solutions of model (3) are stochastically permanent.Proof.
From Theorem6, we have P{|X(t)|>δ}≤ɛ which implies
(57)P{|X(t)|≤δ}≥1-ɛ.
This follows that
(58)liminf⁡t→∞P{|X(t)|≤δ}≥1-ɛ.
By Theorem8, we get
(59)limsup⁡t→∞E(1|X(t)|ϑ)≤Q.
For any ɛ>0, let χ=ɛϑ/Qϑ; then
(60)P{|X(t)|<χ}=P{1|X(t)|>1χ}≤χ1/ϑE(|X(t)|-ϑ).
Hence,
(61)limsup⁡t→∞P{|X(t)|<χ}≤χ1/ϑQ=ɛ,
which follows that
(62)liminf⁡t→∞P{|X(t)|≥χ}≥1-ɛ.
The proof is complete.


## 3. Conclusions

In this paper, we propose an SIR epidemic model with a nonlinear incidence rate of the formkSI/(1+αI). We extend to consider and analyze the epidemic model with stochastic perturbations. The value of this study lies in two aspects. First, it presents existence and global stability analysis of the endemic equilibrium for the deterministic model (4). Second, it verifies some relevant properties of the corresponding stochastic model (3) and reveals the effect of environmental noise on the epidemic model.To study the effect of environmental noise on the deterministic model (4), we stochastically perturb model (4) with respect to white noise around its endemic equilibrium. By constructing suitable Lyapunov functions and applying Itô’s formula, we obtain that there is a unique positive solution to model (3) for any positive initial value and derive that the solution is stochastically bounded and permanent under some conditions. These conditions depend on the intensities of noise σ1, σ2, and σ3. When the intensities of noise satisfy some conditions and are not sufficiently large, the population of the stochastic model may be stochastically permanent.As an example, we perform some numerical simulations to illustrate the analytical results of stochastic model (3) by referring to the method mentioned in Higham [34]. Then model (3) can be rewritten as the following discretization equations:
(63)Si+1=Si+(b-dSi-kSiIi1+αIi+γRi)Δt+σ1(Si-S*)Δtξi+σ122(Si-S*)2(ξi2-1)Δt,Ii+1=Ii+(kSiIi1+αIi-(d+μ)Ii)+σ2(Ii-I*)Δtξi+σ222(Ii-I*)2(ξi2-1)Δt,Ri+1=Ri+(μIi-(d+γ)Ri)Δt+σ3(Ri-R*)Δtξi+σ322(Ri-R*)2(ξi2-1)Δt,
where ξi(i=1,2,…,n) is the Gaussian random variables N(0,1).Figure1 shows time-series plots for model (3) with and without stochastic perturbations. The parameters are taken as b=1, d=0.2, k=1, α=0.5, γ=0.25, and μ=0.3 and initial value (S(0),I(0),R(0))=(1.35,0.9,0.45). In this case, model (4) has the endemic point E*=(1.087,2.3478,1.5652). The only difference between conditions of Figures 1(a) and 1(b) is that the values of environmental noise intensities σ1, σ2, and σ3 are different. In Figure 1(a), with σ1=0.05, σ2=0.01, and σ3=0.03 and in Figure 1(b), with σ1=0.1, σ2=0.06,  and   σ3=0.12, the condition of Theorem 9 is satisfied. That is, the solutions of model (3) are stochastically permanent. From Figures 1(a) and 1(b), one can see that with increasing the noise intensities, the solutions of model (3) will be oscillating strongly around the endemic point E* of model (4).Solutions of model (3) with different noise. Other parameters and initial condition are given in text. (a) σ1=0.05,σ2=0.01, and σ3=0.03 and (b) σ1=0.1, σ2=0.06, andσ3=0.12.
(a)
(b)To study the effect of noise in model (3) further, in Figure 2(a), we choose σ1=σ3=0, σ2=0.06 which satisfies the condition of Theorem 9, while in Figure 2(b), σ1=1.7, σ2=σ3=0 that does not satisfy the condition of Theorem 9. From Figure 2(a), one can see that the infective population I will be oscillating slightly around I*=2.3478, and both the susceptible S and the removed R population will be affected by the noise but the effect is very small. From Figure 2(b), when the condition of Theorem 9 is not satisfied, the noise can force the population to become largely fluctuating. In this case, the solution of model (3) is not stochastically permanent.Solutions of model (3) with different noise. Other parameters and initial condition are the same as Figure 1. (a) σ1=0, σ2=0.06, and σ3=0 and (b) σ1=1.7 and σ2=σ3=0.
(a)
(b)



---

*Source: 356013-2014-03-02.xml*