# Dynamics of a Stochastic SIS Epidemic Model with Saturated Incidence

Authors: Can Chen; Yanmei Kang
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We introduce stochasticity into the SIS model with saturated incidence. The existence and uniqueness of the positive solution are proved by employing the Lyapunov analysis method. Then, we carry out a detailed analysis on both its almost sure exponential stability and itspth moment exponential stability, which indicates that the pth moment exponential stability implies the almost sure exponential stability. Additionally, the results show that the conditions for the disease to become extinct are much weaker than those in the corresponding deterministic model. The conditions for the persistence in the mean and the existence of a stationary distribution are also established. Finally, we derive the expressions for the mean and variance of the stationary distribution. Compared with the corresponding deterministic model, the threshold value for the disease to die out is affected by the half saturation constant. That is, increasing the saturation effect can reduce the disease transmission. Computer simulations are presented to illustrate our theoretical results.


---

## Body

## 1. Introduction

Epidemiology has been investigated by mathematicians through establishing mathematical models for a long time; see, for instance, [1–4]. Mathematical models can help people to better understand the spread of the disease and thus take effective measures to reduce its transmission as much as possible. Particularly, the classical SIS epidemic model [5, 6] is often used to model the dynamics of the diseases with no protective immunity such as gonorrhea.Usually, the incidence function describes the number of new infections per unit time and it plays a vital role in the deterministic models. As it is known, in the most existing literature, the bilinear incidence rateβSI is frequently used. However, this kind of incidence rate has some limitations, since it does not consider the behavioral change of susceptible individuals or control measures taken by the government when the number of infective individuals gets large. For example, people would wash their hands frequently, wear masks, and limit their time of going out, and the government would also take measures such as quarantine and isolation when the A/H1N1 influenza became serious in 2009 [7]. Therefore, it is more reasonable to adopt saturated incidence rate βSI/(1+hI) [8] in the mathematical models. When the number of infective individuals gets large, βI/(1+hI) tends to a saturation effect, where βI is the infection force of the disease and 1/(1+hI) is the inhibition effect from the behavioral change or control measures. Here, h represents the half saturation constant [9]. With the saturated incidence taken into consideration, the classical SIS epidemic model is transformed into
(1)dS(t)dt=μN-βS(t)I(t)1+hI(t)+γI(t)-μS(t),dI(t)dt=βS(t)I(t)1+hI(t)-(μ+γ)I(t),
where S(t) and I(t) are the number of susceptible individuals and infective individuals at time t, respectively. The total population N=S(t)+I(t) is a constant. Consider the following parameters:  β is the transmission rate, μ is the birth and death rate, and γ is the recovery rate. Assume all these parameters are positive; then, the classical stability analysis has shown that the system (1) always has a disease free equilibrium:
(2)E0=(S0,I0)=(N,0).
The basic reproduction number, R0d=βN/(μ+γ), is a threshold value, which determines the extinction and persistence of the disease. When R0d≤1, the disease free equilibrium E0 is globally stable. When R0d>1, the disease free equilibrium E0 is unstable, and model (1) has a unique endemic equilibrium:
(3)E*=(S*,I*)=((1+hN)(μ+γ)β+h(μ+γ),βN-(μ+γ)β+h(μ+γ)),
which is globally stable.Since the noise exists almost everywhere, epidemic models are inevitably affected by it. Hence, it is more reasonable to introduce random perturbations into mathematical models. At the same time, we note that taking the effect of environment noise on epidemic models into consideration has been a popular trend in disease spread modeling [10–13]. Gray et al. presented a stochastic SIS epidemic model in [14]. They established conditions for extinction and persistence of the disease and the existence of a stationary distribution. Stochastic SIR and SIRS models with and without distributed time delay were investigated in [15, 16], respectively, where the asymptotic behavior was discussed. In [17], stochastically perturbed SIR and SEIR epidemic models with saturated incidence were studied and the results of extinction and ergodicity were concluded. Ji et al. investigated a stochastic multigroup SIR model with fluctuations around the transmission coefficient and the death rate of each subpopulation in [18, 19] separately. Zhao et al. [20] discussed the extinction and persistence of the stochastic SIS epidemic model with vaccination. However, among all these researches, there are few literatures considering the parameter perturbation in SIS epidemic model with saturated incidence.Motivated by the above consideration, in this paper, we introduce random effect into the SIS model (1) by replacing the transmission rate β with β+σdW(t), where W(t) is standard Brownian motion and σ is the intensity of white noise. The stochastic version corresponding to the deterministic model (1) can be represented as
(4)dS(t)=(μN-βS(t)I(t)1+hI(t)+γI(t)-μS(t))dtdS(t)=-σS(t)I(t)1+hI(t)dW(t),dI(t)=(βS(t)I(t)1+hI(t)-(μ+γ)I(t))dt+σS(t)I(t)1+hI(t)dW(t).
Given that S(t)+I(t)=N, it is sufficient to study the following equation:
(5)dI(t)=(β(N-I(t))I(t)1+hI(t)-(μ+γ)I(t))dt+σ(N-I(t))I(t)1+hI(t)dW(t),
with initial value I(0)=I0∈(0,N). In the following sections, we will restrict ourselves to (5).The paper is organized as follows. In Section2, the existence of unique positive solution is shown. Both its almost sure exponential stability and its pth moment exponential stability are then investigated in Section 3. The conditions for persistence in the mean of the disease are established in Section 4. The existence of a stationary distribution and the expressions for the mean and variance are presented in Section 5. In Section 6, we give a brief conclusion. Besides, the computer simulations which support our results are given in each section. Finally, we give an appendix containing some theory used in the previous sections.


## 2. Existence of Unique Positive Solution

To investigate the dynamical behavior of the epidemic model, we need to show that the model has a unique global solution and the solution will remain within(0,N) whenever it starts there. Hence, in this section, employing the Lyapunov analysis method (mentioned in [18, 21]), we establish Theorem 1.Theorem 1.
There is a unique global solutionI(t)∈(0,N) of system (5) on t≥0 for any given initial value I(0)=I0∈(0,N) with probability 1; namely,
(6)P{I(t)∈(0,N),∀t≥0}=1.Proof.
Since the coefficients of (5) are locally Lipschitz continuous for any given initial value I(0)=I0∈(0,N), there is a unique local solution I(t) on t∈[0,τe), where τe is the explosion time [22]. To show that the solution is global, we need to show that τe=∞ a.s. Let m0>0 be sufficiently large such that 1/m0<I0<N-1/m0. For each integer m≥m0, define the stopping time
(7)τm=inf⁡{t∈[0,τe):I(t)∉(1m,N-1m)},
where throughout this paper we set inf⁡∅=∞ (as usual ∅ denotes the empty set). Clearly, τm is increasing as m→∞. Set τ∞=lim⁡m→∞τm, whence τ∞≤τe a.s. If we can show that τ∞=∞ a.s., then τe=∞ a.s. and I(t)∈(0,N) a.s. for all t≥0. In other words, to complete the proof, all what we need to show is that τ∞=∞ a.s. If this statement is not true, then there exists a pair of constants T>0 and ϵ∈(0,1) such that
(8)P{τ∞≤T}>ϵ.
Therefore, there is an integer m1≥m0 such that
(9)P{τm≤T}≥ϵ,∀m≥m1.
Define a function V:(0,N)→R+ by
(10)V(x)=1x+1N-x.
Using the Itô formula, we get
(11)dV(x)=(-1x2+1(N-x)2)dx+(1x3+1(N-x)3)(dx)2=LV(x)dt+(-σ(N-x)x(1+hx)+σx(N-x)(1+hx))dW(t),
where LV:(0,N)→R is defined by
(12)LV=(-1x2+1(N-x)2)(β(N-x)x1+hx-(μ+γ)x)+(1x3+1(N-x)3)σ2(N-x)2x2(1+hx)2=-β(N-x)x(1+hx)+μ+γx+βx(N-x)(1+hx)+σ2(N-x)2x(1+hx)2+σ2x2(N-x)(1+hx)2-(μ+γ)x(N-x)2≤μ+γx+βx(N-x)(1+hx)+σ2(N-x)2x(1+hx)2+σ2x2(N-x)(1+hx)2≤1x(μ+γ+σ2N2)+1N-x(βh+σ2h2)≤CV(x),
where C=(μ+γ+σ2N2)∨(β/h+σ2/h2). Therefore,
(13)∫0τm∧t‍dV(I(s))=∫0τm∧t‍LV(I(s))ds+∫0τm∧t‍(-σ(N-I(s))I(s)(1+hI(s))00000000000+σI(s)(N-I(s))(1+hI(s)))dW(s)≤∫0τm∧t‍CV(I(s))ds+∫0τm∧t‍(-σ(N-I(s))I(s)(1+hI(s))00000000000+σI(s)(N-I(s))(1+hI(s)))dW(s),
which implies that
(14)EV(I(τm∧t))≤V(I0)+CE∫0τm∧t‍V(I(s))ds≤V(I0)+C∫0t‍EV(I(τm∧s))ds.
The Gronwall inequality yields that
(15)EV(I(τm∧T))≤V(I0)eCT.
Set Ωm={τm≤T} for all m≥m1. Then, by (9), P(Ωm)≥ϵ. For every ω∈Ωm, I(τm,ω) equals either 1/m or N-1/m, and therefore
(16)V(I(τm,ω))≥m.
It then follows from (9) and (15) that
(17)V(I0)eCT≥E[1Ωm(ω)V(I(τm,ω))]≥ϵm,
where 1Ωm(ω) is the indicator function of Ωm. Letting m→∞ leads to the contradiction that
(18)∞>V(I0)eCT=∞.
So, we must have τ∞=∞ a.s., whence the proof is complete.


## 3. Extinction

In discussing the extinction of system (5), we focus on the two kinds of exponential stabilities, almost sure exponential stability and pth moment exponential stability.


### 3.1. Almost Sure Exponential Stability

Theorem 2.
LetI(t) be the solution of system (5) with initial value I(0)=I0∈(0,N). If (a)
R01s:=βN/(μ+γ)-σ2N2/(2(μ+γ)(1+hN)2)=R0d-σ2N2/(2(μ+γ)(1+hN)2)<1 and σ2≤β(1+hN)2/N,(b)
σ2>β(1+hN)2/N∨β2(1+hN)2/2(μ+γ).then,limsup⁡t→∞(log⁡(I(t))/t)≤βN-μ-γ-σ2N2/2(1+hN)2<0a.s. if (a) holds;limsup⁡t→∞(log⁡(I(t))/t)≤β2(1+hN)2/2σ2-μ-γ<0a.s. if (b) holds;that is, I(t) tends to zero exponentially a.s.; namely, the disease will die out with probability one.Proof.
Applying Itô formula to system (5) leads to
(19)dlog⁡(I(t))=F(I(t))dt+σ(N-I(t))1+hI(t)dW(t),
where F(I(t))=β(N-I(t))/(1+hI(t))-μ-γ-σ2(N-I(t))2/2(1+hI(t))2. Integrating both sides of (19) from 0 to t, we have
(20)log⁡(I(t))=log⁡(I0)+∫0t‍F(I(s))ds+M(t),
where M(t)=∫0t‍(σ(N-I(s))/(1+hI(s)))dW(s). Obviously,
(21)F(I(s))≤βN-βI(s)-μ-γ-σ2(N-I(s))22(1+hN)2:=G(I(s)),
where G(I(s))=βN-μ-γ-σ2N2/2(1+hN)2-(β-σ2N/(1+hN)2)I(s)-(σ2/2(1+hN)2)I(s)2, for I(s)∈(0,N).
If condition (a) is satisfied, it then follows from (20) that
(22)log⁡(I(t))≤log⁡(I0)+(βN-μ-γ-σ2N22(1+hN)2)t+M(t).
This implies that
(23)limsup⁡t→∞log⁡⁡(I(t))t≤βN-μ-γ-σ2N22(1+hN)2+limsup⁡t→∞M(t)ta.s.
Clearly, M(t)=∫0t‍(σ(N-I(s))/(1+hI(s)))dW(s) is a local martingale and vanishes at t=0. Moreover,
(24)limsup⁡t→∞〈M,M〉tt≤σ2N2<∞a.s.
By strong law of large numbers [22], we get
(25)lim⁡t→∞M(t)t=0a.s.
Then,
(26)limsup⁡t→∞log⁡(I(t))t≤βN-μ-γ-σ2N22(1+hN)2<0a.s.
If condition (b) is satisfied, the functionG(x) takes its maximum value G(x~) at x~=N-β(1+hN)2/σ2∈(0,N)); that is, G(x~)=β2(1+hN)2/2σ2-(μ+γ). It then follows from (20) that
(27)log⁡(I(t))≤log⁡(I0)+(β2(1+hN)22σ2-(μ+γ))t+M(t).
Then,
(28)limsup⁡t→∞log⁡(I(t))t≤β2(1+hN)22σ2-μ-γ<0a.s.
This finishes the proof.From Theorem2, we know that when R01s=R0d-σ2N2/2(μ+γ)(1+hN)2<1 and the perturbation is small, the disease will die out a.s., whereas the disease may still persist for the corresponding deterministic model (1). If the noise intensity σ2 is larger than β(1+hN)2/N∨β2(1+hN)2/2(μ+γ), then the disease will also die out a.s. For h=0, system (5) becomes the model discussed in [14] and the conditions for the disease to become extinct we obtain here are in agreement with those derived in [14]. If there exists no environment noise, that is, σ=0, then R01s=R0d is the basic reproduction number of the deterministic model (1). In addition, the half saturation constant also influences R01s, which shows that increasing the saturation effect can reduce the value of R01s and thus can reduce the disease spread. The following simulations indicate these results.Example 3.
We assume that the unit time is one day. The parameters are given by(29)N=100,β=0.013,μ=0.4,γ=0.6,h=0.05,σ=0.06.
Note that
(30)R01s=βNμ+γ-σ2N22(μ+γ)(1+hN)2=0.8<1,σ2=0.0036≤β(1+hN)2N=0.0047.
That is, condition (a) is satisfied and the solution I(t) satisfies
(31)limsup⁡t→∞log⁡(I(t))t≤βN-μ-γ-σ2N22(1+hN)2=-0.2<0a.s.
Hence, I(t) tends to zero exponentially a.s.; that is, the disease will die out a.s.
On the other hand, for the deterministic model (1), the basic reproduction number is
(32)R0d=βNμ+γ=1.3>1.
The solution I(t) obeys
(33)lim⁡t→∞I(t)=βN-(μ+γ)β+(μ+γ)h=4.7619.
Employing the Euler-Maruyama (EM) method in [23], the computer simulations are given in Figure 1, which support our results.Computer simulations ofI(t) for model (5) and the deterministic model (1), using the EM method with step size 0.001 and with initial values (a) I(0)=80 and (b) I(0)=10.
(a)
(b)Example 4.
We keep the parameters the same as Example3 but increase σ to 0.09. Note that
(34)σ2=0.0081>max⁡⁡{β(1+hN)2N=0.00468,β2(1+hN)22(μ+γ)=0.00304}=0.00468.
That is, condition (b) is satisfied and
(35)limsup⁡t→∞log⁡(I(t))t≤β2(1+hN)22σ2-μ-γ=-0.624<0a.s.
Therefore, I(t) tends to zero exponentially a.s.; namely, the disease will die out a.s. The simulations are shown in Figure 2 to support our results completely.Computer simulations ofI(t) for model (5) and the deterministic model (1), using the EM method with step size 0.001 and with initial values (a) I(0)=80 and (b) I(0)=10.
(a)
(b)


### 3.2.pth Moment Exponential Stability

Following (20), we obtain
(36)I(t)=I0exp⁡(∫0t‍F(I(s))ds+M(t)).
As W(t)~N(0,t), then [24]
(37)EeM(t)=Eexp⁡(∫0t‍σ(N-I(s))1+hI(s)dW(s))≤EeσNW(t)=e(σ2N2/2)t.
Using the fact that
(38)F(I(s))≤βN-μ-γ,
we find
(39)EI(t)=EI0exp⁡(∫0t‍F(I(s))ds+M(t))≤I0exp⁡((βN-μ-γ+σ2N22)t).
Since
(40)log⁡(I(t))p=plog⁡(I(t))≤p(log⁡(I0)+(βN-μ-γ)t+σNW(t)),
then
(41)E(I(t))p≤(I0)p·ep(βN-μ-γ)t·e(p2σ2N2/2)t=(I0)pexp⁡(p(βN-μ-γ+pσ2N22)t).Theorem 5.
Assume thatβN-μ-γ+(pσ2N2/2)<0; the solution I(t) with initial value I(0)=I0∈(0,N) is pth moment exponentially stable.Next, we will make a comparison between the almost sure exponential stability and thepth moment exponential stability, since
(42)x(β(N-x)x1+hx-(μ+γ)x)=(β(N-x)1+hx-(μ+γ))x2≤(βN-(μ+γ))x2.
From Theorem 5, we know that βN-(μ+γ)<0; that is,
(43)x(β(N-x)x1+hx-(μ+γ)x)≤(βN-(μ+γ))x2≤0.
Moreover,
(44)σ2(N-x)2x2(1+hx)2≤σ2(N-x)2x2≤σ2N2x2.
Let K=σ2N2 be positive; then
(45)x(β(N-x)x1+hx-(μ+γ)x)∨σ2(N-x)2x2(1+hx)2≤Kx2.
Therefore, the condition of Lemma A.4 (in the appendix) [22] is satisfied, which suggests that the pth moment exponential stability of system (5) implies the almost sure exponential stability.


## 3.1. Almost Sure Exponential Stability

Theorem 2.
LetI(t) be the solution of system (5) with initial value I(0)=I0∈(0,N). If (a)
R01s:=βN/(μ+γ)-σ2N2/(2(μ+γ)(1+hN)2)=R0d-σ2N2/(2(μ+γ)(1+hN)2)<1 and σ2≤β(1+hN)2/N,(b)
σ2>β(1+hN)2/N∨β2(1+hN)2/2(μ+γ).then,limsup⁡t→∞(log⁡(I(t))/t)≤βN-μ-γ-σ2N2/2(1+hN)2<0a.s. if (a) holds;limsup⁡t→∞(log⁡(I(t))/t)≤β2(1+hN)2/2σ2-μ-γ<0a.s. if (b) holds;that is, I(t) tends to zero exponentially a.s.; namely, the disease will die out with probability one.Proof.
Applying Itô formula to system (5) leads to
(19)dlog⁡(I(t))=F(I(t))dt+σ(N-I(t))1+hI(t)dW(t),
where F(I(t))=β(N-I(t))/(1+hI(t))-μ-γ-σ2(N-I(t))2/2(1+hI(t))2. Integrating both sides of (19) from 0 to t, we have
(20)log⁡(I(t))=log⁡(I0)+∫0t‍F(I(s))ds+M(t),
where M(t)=∫0t‍(σ(N-I(s))/(1+hI(s)))dW(s). Obviously,
(21)F(I(s))≤βN-βI(s)-μ-γ-σ2(N-I(s))22(1+hN)2:=G(I(s)),
where G(I(s))=βN-μ-γ-σ2N2/2(1+hN)2-(β-σ2N/(1+hN)2)I(s)-(σ2/2(1+hN)2)I(s)2, for I(s)∈(0,N).
If condition (a) is satisfied, it then follows from (20) that
(22)log⁡(I(t))≤log⁡(I0)+(βN-μ-γ-σ2N22(1+hN)2)t+M(t).
This implies that
(23)limsup⁡t→∞log⁡⁡(I(t))t≤βN-μ-γ-σ2N22(1+hN)2+limsup⁡t→∞M(t)ta.s.
Clearly, M(t)=∫0t‍(σ(N-I(s))/(1+hI(s)))dW(s) is a local martingale and vanishes at t=0. Moreover,
(24)limsup⁡t→∞〈M,M〉tt≤σ2N2<∞a.s.
By strong law of large numbers [22], we get
(25)lim⁡t→∞M(t)t=0a.s.
Then,
(26)limsup⁡t→∞log⁡(I(t))t≤βN-μ-γ-σ2N22(1+hN)2<0a.s.
If condition (b) is satisfied, the functionG(x) takes its maximum value G(x~) at x~=N-β(1+hN)2/σ2∈(0,N)); that is, G(x~)=β2(1+hN)2/2σ2-(μ+γ). It then follows from (20) that
(27)log⁡(I(t))≤log⁡(I0)+(β2(1+hN)22σ2-(μ+γ))t+M(t).
Then,
(28)limsup⁡t→∞log⁡(I(t))t≤β2(1+hN)22σ2-μ-γ<0a.s.
This finishes the proof.From Theorem2, we know that when R01s=R0d-σ2N2/2(μ+γ)(1+hN)2<1 and the perturbation is small, the disease will die out a.s., whereas the disease may still persist for the corresponding deterministic model (1). If the noise intensity σ2 is larger than β(1+hN)2/N∨β2(1+hN)2/2(μ+γ), then the disease will also die out a.s. For h=0, system (5) becomes the model discussed in [14] and the conditions for the disease to become extinct we obtain here are in agreement with those derived in [14]. If there exists no environment noise, that is, σ=0, then R01s=R0d is the basic reproduction number of the deterministic model (1). In addition, the half saturation constant also influences R01s, which shows that increasing the saturation effect can reduce the value of R01s and thus can reduce the disease spread. The following simulations indicate these results.Example 3.
We assume that the unit time is one day. The parameters are given by(29)N=100,β=0.013,μ=0.4,γ=0.6,h=0.05,σ=0.06.
Note that
(30)R01s=βNμ+γ-σ2N22(μ+γ)(1+hN)2=0.8<1,σ2=0.0036≤β(1+hN)2N=0.0047.
That is, condition (a) is satisfied and the solution I(t) satisfies
(31)limsup⁡t→∞log⁡(I(t))t≤βN-μ-γ-σ2N22(1+hN)2=-0.2<0a.s.
Hence, I(t) tends to zero exponentially a.s.; that is, the disease will die out a.s.
On the other hand, for the deterministic model (1), the basic reproduction number is
(32)R0d=βNμ+γ=1.3>1.
The solution I(t) obeys
(33)lim⁡t→∞I(t)=βN-(μ+γ)β+(μ+γ)h=4.7619.
Employing the Euler-Maruyama (EM) method in [23], the computer simulations are given in Figure 1, which support our results.Computer simulations ofI(t) for model (5) and the deterministic model (1), using the EM method with step size 0.001 and with initial values (a) I(0)=80 and (b) I(0)=10.
(a)
(b)Example 4.
We keep the parameters the same as Example3 but increase σ to 0.09. Note that
(34)σ2=0.0081>max⁡⁡{β(1+hN)2N=0.00468,β2(1+hN)22(μ+γ)=0.00304}=0.00468.
That is, condition (b) is satisfied and
(35)limsup⁡t→∞log⁡(I(t))t≤β2(1+hN)22σ2-μ-γ=-0.624<0a.s.
Therefore, I(t) tends to zero exponentially a.s.; namely, the disease will die out a.s. The simulations are shown in Figure 2 to support our results completely.Computer simulations ofI(t) for model (5) and the deterministic model (1), using the EM method with step size 0.001 and with initial values (a) I(0)=80 and (b) I(0)=10.
(a)
(b)


## 3.2.pth Moment Exponential Stability

Following (20), we obtain
(36)I(t)=I0exp⁡(∫0t‍F(I(s))ds+M(t)).
As W(t)~N(0,t), then [24]
(37)EeM(t)=Eexp⁡(∫0t‍σ(N-I(s))1+hI(s)dW(s))≤EeσNW(t)=e(σ2N2/2)t.
Using the fact that
(38)F(I(s))≤βN-μ-γ,
we find
(39)EI(t)=EI0exp⁡(∫0t‍F(I(s))ds+M(t))≤I0exp⁡((βN-μ-γ+σ2N22)t).
Since
(40)log⁡(I(t))p=plog⁡(I(t))≤p(log⁡(I0)+(βN-μ-γ)t+σNW(t)),
then
(41)E(I(t))p≤(I0)p·ep(βN-μ-γ)t·e(p2σ2N2/2)t=(I0)pexp⁡(p(βN-μ-γ+pσ2N22)t).Theorem 5.
Assume thatβN-μ-γ+(pσ2N2/2)<0; the solution I(t) with initial value I(0)=I0∈(0,N) is pth moment exponentially stable.Next, we will make a comparison between the almost sure exponential stability and thepth moment exponential stability, since
(42)x(β(N-x)x1+hx-(μ+γ)x)=(β(N-x)1+hx-(μ+γ))x2≤(βN-(μ+γ))x2.
From Theorem 5, we know that βN-(μ+γ)<0; that is,
(43)x(β(N-x)x1+hx-(μ+γ)x)≤(βN-(μ+γ))x2≤0.
Moreover,
(44)σ2(N-x)2x2(1+hx)2≤σ2(N-x)2x2≤σ2N2x2.
Let K=σ2N2 be positive; then
(45)x(β(N-x)x1+hx-(μ+γ)x)∨σ2(N-x)2x2(1+hx)2≤Kx2.
Therefore, the condition of Lemma A.4 (in the appendix) [22] is satisfied, which suggests that the pth moment exponential stability of system (5) implies the almost sure exponential stability.


## 4. Persistence

We will investigate persistence in the mean for system (5) in this section. The definition of the persistence in the mean was initially proposed for the deterministic system [25] and was also used for the stationary stochastic system with ergodicity [19, 20]. According to the ergodicity for the stationary process, the time average in the long time limit is equal to the ensemble expectation over phase space, so we introduce the notation
(46)〈I(t)〉=1t∫0t‍I(s)ds.Definition 6.
System (5) is said to be persistent in the mean, if
(47)liminf⁡t→∞1t∫0t‍I(s)ds>0a.s.Theorem 7.
If(48)R02s≔βN(μ+γ)(1+hN)-σ2N22(μ+γ)>1,σ2<β(1+hN)2N,
then, for any initial value I(0)=I0∈(0,N), the solution I(t) of system (5) has the following property:
(49)I_≤liminf⁡t→∞〈I(t)〉≤limsup⁡t→∞〈I(t)〉≤I¯a.s.,
where
(50)I¯=βN-μ-γ-(σ2N2/2(1+hN)2)β-(σ2N/(1+hN)2),I_=βN/(1+hN)-μ-γ-(σ2N2/2)β/(1+hN).Proof.
By (20), we obtain
(51)log⁡⁡(I(t))-log⁡⁡(I0)t≤1t∫0t‍(βN-βI(s)-μ-γ-σ2(N-I(s))22(1+hN)2)ds+M(t)t=βN-μ-γ-σ2N22(1+hN)2-β〈I(t)〉+σ2N(1+hN)2〈I(t)〉-σ22(1+hN)2〈I(t)2〉+M(t)t≤βN-μ-γ-σ2N22(1+hN)2-(β-σ2N(1+hN)2)〈I(t)〉+M(t)t.
The inequality can be written as
(52)log⁡⁡(I(t))t≤βN-μ-γ-σ2N22(1+hN)2-(β-σ2N(1+hN)2)〈I(t)〉+M(t)t+log⁡⁡(I0)t.
Since R01s≥R02s, if condition (48) is satisfied, then
(53)βN-μ-γ-σ2N22(1+hN)2>0,β-σ2N(1+hN)2>0.
Together with Lemma A.2. in [20], we have
(54)lim⁡t→∞⁡sup⁡〈I(t)〉≤βN-μ-γ-(σ2N2/2(1+hN)2)β-(σ2N/(1+hN)2)=I¯a.s.
On the other hand,
(55)log⁡⁡(I(t))-log⁡⁡(I0)t≥1t∫0t‍(β(N-I(s))1+hN-μ-γ-12σ2(N-I(s))2)ds+M(t)t≥βN1+hN-μ-γ-12σ2N2-β1+hN〈I(t)〉+M(t)t.
That is,
(56)〈I(t)〉≥1β/(1+hN)(log⁡⁡(I(t))-log⁡⁡(I0)tβN1+hN-μ-γ-12σ2N200000000000000+M(t)t-log⁡⁡(I(t))-log⁡⁡(I0)t).
Since -∞<log⁡(I(t))<log⁡(N), then
(57)liminf⁡t→∞〈I(t)〉≥βN/(1+hN)-μ-γ-(σ2N2/2)β/(1+hN)=I_∈(0,N)a.s.
Moreover,
(58)I¯=βN-μ-γ-(σ2N2/2(1+hN)2)β-(σ2N/(1+hN)2)≥I_=βN/(1+hN)-μ-γ-(σ2N2/2)β/(1+hN).
Therefore, we complete the proof.Example 8.
The parameters are given by(59)N=100,β=0.2,μ=0.4,γ=0.6,h=0.05,σ=0.02.
Note that
(60)R02s≔βN(μ+γ)(1+hN)-σ2N22(μ+γ)=1.33>1,σ2=0.0004<β(1+hN)2N=0.072.
Then the solution I(t) of system (5) satisfies
(61)10≤liminf⁡t→∞〈I(t)〉≤limsup⁡t→∞〈I(t)〉≤95.2938a.s.
For the corresponding deterministic model (1),
(62)R0d=βN(μ+γ)=20>1,lim⁡t→∞I(t)=βN-(μ+γ)β+h(μ+γ)=76.
The simulation results are shown in Figure 3.Computer simulations ofI(t) for model (5) (red line) and the deterministic model (1) (blue line), using the EM method with step size 0.001 and with initial values (a) I(0)=80 and (b) I(0)=10. I_=10 and I¯=95.2938.
(a)
(b)


## 5. Stationary Distribution

Applying LemmaA.5 (in the appendix) [26], if we can show that conditions (B.1) and (B.2) are satisfied, then system (5) has a stationary distribution. Obviously, the square of the diffusion coefficient of system (5), that is, σ2(N-I)2I2/(1+hI)2, is bounded away from zero for I∈(a,b)⊂(0,N)); then, condition (B.1) is satisfied. Next, we will prove that (B.2) is also valid.Theorem 9.
Assume thatR02s:=βN/((μ+γ)(1+hN))-σ2N2/(2(μ+γ))>1 and σ2<(2β/(N2(1+hI*)I*))min⁡{(I*)2,(N-I*)2}; then, there is a unique stationary distribution of system (5). Here, I* is the unique endemic equilibrium of system (1); that is, I*=(βN-(μ+γ))/(β+h(μ+γ)). Moreover,
(63)lim⁡t→∞1tE∫0t‍β1+hI*(I(s)-I*)2ds≤12I*σ2N2.Proof.
IfR02s:=βN/((μ+γ)(1+hN))-σ2N2/2(μ+γ)>1, it is clear that R0d=βN/(μ+γ)>1; then, model (1) has a unique endemic equilibrium I*; that is,
(64)β(N-I*)1+hI*=μ+γ,I*=βN-(μ+γ)β+h(μ+γ).
Define
(65)V(I)=I-I*-I*log⁡(II*).
Then V is positive definite. By Itô formula, we yield
(66)dV=(1-II*)dI+I*2I2(dI)2=LVdt+σ(N-I)(I-I*)1+hIdW(t),
where
(67)LV=(β(N-I)1+hI-(μ+γ))(I-I*)+I*σ2(N-I)22(1+hI)2=(β(N-I)1+hI-β(N-I*)1+hI*)(I-I*)+I*σ2(N-I)22(1+hI)2=-β(1+hN)(I-I*)2(1+hI)(1+hI*)+I*σ2(N-I)22(1+hI)2≤-β1+hI*(I-I*)2+12I*σ2N2.
Note that if σ2<(2β/N2(1+hI*)I*)min⁡{(I*)2,(N-I*)2}, then the neighborhood U:-(β/(1+hI*))(I-I*)2+(1/2)I*σ2N2=0 lies entirely in (0,N). Hence, for I∈(0,N)∖U, LV≤-k (k is a positive constant), which implies that condition (B.2) is satisfied; for more details, we refer to [27]. Therefore, system (5) has a stationary distribution.
Moreover,(68)0≤E∫0t‍dV(I(s))=E∫0t‍LV(I(s))ds≤-E∫0t‍β1+hI*(I(s)-I*)2ds+12I*σ2N2t.
Then
(69)lim⁡t→∞1tE∫0t‍β1+hI*(I(s)-I*)2ds≤12I*σ2N2.Next, we will give the expressions for the mean and variance of the stationary distribution.Theorem 10.
IfR02s:=βN/((μ+γ)(1+hN))-σ2N2/(2(μ+γ))>1 and σ2<(2β/(N2(1+hI*)I*))min⁡{(I*)2,(N-I*)2}. Let m and v denote the mean and variance of the stationary distribution of model (5). Then,
(70)m=(βN-μ-γ-12σ2N2)×((βh+(μ+γ)h2+12σ2)βN-μ-γβ+(μ+γ)h-(β(hN-1)-2h(μ+γ)+σ2N)(βh+(μ+γ)h2βN-μ-γβ+(μ+γ)h+12σ2))-1,v=βN-μ-γβ+(μ+γ)hm-m2.Proof.
Multiply both sides of (5) by 1+hI(t); that is,
(71)(1+hI(t))dI(t)=(β(N-I(t))I(t)-(μ+γ)I(t)(1+hI(t)))dt+σ(N-I(t))I(t)dW(t).
Integrating (71) from 0 to t and dividing by t on both sides, then
(72)1t∫0t‍(1+hI(s))dI(s)=1t∫0t‍(β(N-I(s))I(s)-(μ+γ)I(s)(1+hI(s)))ds+1t∫0t‍σ(N-I(s))I(s)dW(s).
Since
(73)I(t)-I0t≤1t∫0t‍(1+hI(s))dI(s)≤(1+hN)I(t)-I0t,
letting t→∞, we have
(74)lim⁡t→∞1t∫0t‍(1+hI(s))dI(s)=0a.s.
Then, applying the ergodic property of the stationary distribution and strong law of large numbers, we get
(75)(βN-μ-γ)m-(β+(μ+γ)h)〈I(t)2〉=0.
Here, 〈I(t)2〉 denotes the second moment of the stationary distribution [14].
In a similar way, multiplying both sides of (19) by (1+hI(t))2, we obtain that
(76)(1+hI(t))2dlog⁡(I(t))=(1+hI(t))2F(I(t))dt+σ(N-I(t))(1+hI(t))dW(t).
Then, integrating (76) from 0 to t and dividing by t on both sides,
(77)1t∫0t‍(1+hI(s))2dlog⁡(I(s))=1t∫0t‍(1+hI(s))2F(I(s))ds+1t∫0t‍σ(N-I(s))(1+hI(s))dW(s).
Since
(78)log⁡I(t)-log⁡I0t≤1t∫0t‍(1+hI(s))2dlog⁡(I(s))≤(1+hN)2log⁡I(t)-log⁡I0t,
then
(79)lim⁡t→∞1t∫0t‍(1+hI(s))2dlog⁡(I(s))=0a.s.
Namely,
(80)lim⁡t→∞1t∫0t‍(1+hI(s))2F(I(s))ds+lim⁡t→∞1t∫0t‍σ(N-I(s))(1+hI(s))dW(s)=0a.s.
Therefore,
(81)βN-μ-γ-12σ2N2+(β(hN-1)-2h(μ+γ)+σ2N)m-(βh+(μ+γ)h2+12σ2)〈I(t)2〉=0.
From (75) and (81), we can get
(82)m=(βN-μ-γ-12σ2N2)×((βh+(μ+γ)h2+12σ2)βN-μ-γβ+(μ+γ)h-(β(hN-1)-2h(μ+γ)+σ2N)(βh+(μ+γ)h2+12σ2)βN-μ-γβ+(μ+γ)h)-1,v=βN-μ-γβ+(μ+γ)hm-m2.Example 11.
We keep the parameters the same as Example8 but reduce σ to 0.005, 0.001, and 0.0005. Note that R02s is 3.2083, 3.3283, and 3.3321, respectively, and all of them are bigger than 1. Furthermore,
(83)I*=βN-(μ+γ)β+h(μ+γ)=76,2βN2(1+hI*)I*min⁡⁡{(I*)2,(N-I*)2}=0.000063.
The values of σ are satisfied for the conditions of Theorem 10. The histograms of I(t) are shown in Figure 4, respectively, for different values of σ=0.005, σ=0.001, and σ=0.0005. The simulations were run for 100,000 iterations with step size 0.001, and the first 90,000 iterations were discarded to allow I(t) to reach its recurrent level. The distribution is slightly negative skewed for σ=0.005 and positive skewed for σ=0.001 and σ=0.0005. Additionally, the corresponding sample skewness coefficients are −0.00067, 0.0016, and 0.0017, separately. These simulation results illustrate that the distribution of I(t) has reached the stationary distribution. When σ=0.005, from (70), the mean and variance of the stationary distribution are 75.9710 and 2.2016, respectively, while the sample mean and variance are 76.3754 and 0.3209. For σ=0.001, the mean and variance of the stationary distribution are 75.9988 and 0.0876, compared to the sample mean and variance which are 75.9285 and 0.0157, separately. When σ reduces to 0.0005, the mean and variance of the stationary distribution are 75.9997 and 0.0219 and the corresponding sample mean and variance are 76.0248 and 0.0050. As σ decreases, I(t) becomes more symmetric about I*=76 and the perturbations become much weaker. What is more, the normal quantile-quantile plots in Figure 5 suggest that these data are not far from being normally distributed.Figure 4
Histograms of the solutionI(t) for model (5) for I(0)=10 and differing values of σ=0.005, σ=0.001, and σ=0.0005.Figure 5
Normal quantile-quantile plots of the solutionI(t) for model (5) for I(0)=10 and differing values of σ=0.005, σ=0.001, and σ=0.0005.From (70), we know that when σ=0, the mean value m equals the value of the endemic equilibrium I* and the variance v is zero. Furthermore, if σ increases, the mean value m decreases and the variance v increases as shown in Figure 6. Obviously, the greater the intensity of the noise is, the larger the variance is, that is, the stronger the fluctuations of the distribution are. And this is in accordance with the actual case.Plot the meanm and variance v against σ. The endemic equilibrium I* of the deterministic model (1) is 76.
(a)
(b)


## 6. Conclusion

In this paper, we established a stochastic SIS epidemic model with saturated incidence to investigate the effect of environment noise. After proving the existence and uniqueness of the positive solution, we considered two kinds of stabilities: almost sure exponential stability andpth moment exponential stability. Our deduction shows that the pth moment exponential stability implies the almost sure exponential stability. As for the extinction of the disease, our results indicate that when the noise σ2 is larger than β(1+hN)2/N∨β2(1+hN)2/2(μ+γ), the disease will die out a.s. If the noise is small, the disease will also die out a.s. when R01s<1, instead, the disease will be persistent in the mean when R02s>1. Our investigation shows that the conditions for the disease to become extinct are much weaker than those in the corresponding deterministic model; that is, even if the disease dies out for the stochastic model, it may still persist for the corresponding deterministic model. Additionally, the half saturation constant influences the threshold value R01s, which is not obtained in the deterministic model. Increasing the saturation effect by changing the behavior of the susceptible individuals and taking effective control measures of the government can reduce R01s and thus can reduce the spread of the disease. Finally, we proved the existence of a stationary distribution and derived the expressions for the mean and variance when the noise is small and R02s>1. Our theoretical results were further verified by computer simulations.



---

*Source: 723825-2014-06-15.xml*