# Dynamics of Stochastically Perturbed SIS Epidemic Model  with Vaccination

Authors: Yanan Zhao; Daqing Jiang
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

We introduce stochasticity into an SIS epidemic model with vaccination. The stochasticity in the model is a standard technique in stochastic population modeling. In the deterministic models, the basic reproduction numberR0 is a threshold which  determines the persistence or extinction of the disease. When the perturbation and the disease-related death rate are small, we carry out a detailed analysis on the dynamical behavior of the stochastic model, also regarding of the value of R0. If R0≤1, the solution of the model is oscillating around a steady state, which is the disease-free equilibrium of the corresponding deterministic model, whereas, if R0>1, there is a stationary distribution, which means that the disease will prevail. The results are illustrated by computer simulations.


---

## Body

## 1. Introduction

Epidemiology is the study of the spread of diseases with the objective to trace factors that are responsible for or contribute to their occurrence. Significant progress has been made in the theory and application of epidemiology modeling by mathematical research. Controlling infectious diseases has been an increasingly complex issue in recent years. Vaccination is an important strategy for the elimination of infectious diseases [1–3]. The vaccination enables the vaccinated to acquire a permanent or temporary immunity. When the immunity is temporary, the immunity can be lost after a period of time. This temporary immunity is used in many references [4–7] that assume the process of losing immunity is in the exponential form.In [7], Li and Ma discussed an SIS model with vaccination. The system has following form:
(1)S˙(t)=(1-q)A-βS(t)I(t)-(μ+p)S(t)hh+γI(t)+εV(t),I˙(t)=βS(t)I(t)-(μ+γ+a)I(t),V˙(t)=qA+pS(t)-(μ+ε)V(t).
Here S(t) denotes the number of members of a population who are susceptible to an infection at time t. I(t) denotes the number of members who are infective at time t. V(t) denotes the number of members who are immune to an infection at time t as the result of vaccination. The derivatives dS(t)/dt, dI(t)/dt, and dV(t)/dt are denoted by S˙(t), I˙(t), and V˙(t), respectively.The parameters in the model are summarized in the following list:A: a constant input of new members into the population per unit time;q: a fraction of vaccinated for newborns;β: transmission coefficient between compartments S and I;μ: natural death rate of S, I, and V compartments;p: the proportional coefficient of vaccinated for the susceptible;γ: recovery rate of infectious individuals;ε: the rate of losing their immunity for vaccinated individuals;α: disease-caused death rate of infectious individuals.All parameter values are assumed to be nonnegative andμ,A>0.In Li and Ma [7], system (1) always has the disease-free equilibrium P0=(S0,I0,V0)=((A/μ)((μ(1-q)+ε)/(μ+ε+p)),0,(A/μ)((μq+p)/(μ+ε+p))). If R0⩽1, then P0 is the unique equilibrium of (1), and it is globally stable in the invariant set Γ, where Γ={(S,I,V):S>0,I≥0,V≥0,S+I+V≤A/μ}. If R0>1, then P0 is unstable and there is an endemic equilibrium P*=(S*,I*,V*)=((μ+γ+α)/β,(μ(μ+γ+α)(μ+ε+p)/β(μ+α)(μ+ε))(R0-1),(qA+(p/β)(μ+γ+α))/(μ+ε)) of (1), which is globally asymptotically stable under the sufficient condition in the invariant set Γ. The basic reproduction number of system (1) is as follows:
(2)R0=Aβ(μ(1-q)+ε)μ(μ+γ+α)(μ+ε+p)=βS0μ+γ+α.Considering real life is full of randomness and stochasticity, using stochastic models can gain more real benefits. Consequently, there are many authors who studied stochastic biological systems and stochastic epidemic models; see [8–24]. There may be different approaches which result in different effects on the population system to include random perturbations in the models. In [20–23], the situation of the parameter perturbation was considered. Gray et al. [23] discuss the following stochastic SIS model:
(3)dS(t)=[μN-βS(t)I(t)+γI(t)-μS(t)]dt-σS(t)I(t)dB(t),dI(t)=[βS(t)I(t)-(μ+γ)I(t)]dt+σS(t)I(t)dB(t).
Here S and I denote the susceptible and infected numbers of individuals, respectively. With S(t)+I(t)≡N, they simplify system (3) into a single equation and show that if R0S<1 and σ2≤β/N, then the disease will die out with probability one. If R0S>1, then the disease will be persistent, where R0S=R0D-(σ2N2/2(μ+γ)) and R0D=βN/(μ+γ) is the threshold of the corresponding deterministic system. R0S can be considered as the threshold of system (3), which is less than the value of R0D.On the other hand, white noise stochastic perturbations around the positive endemic equilibrium of epidemic models were considered in [19, 25]. Beretta et al. [25] considered the following system differential equation:
(4)dS(t)=-βS(t)∫0h‍f(s)I(t-s)ds-μ1S(t)+b+σ1(S(t)-S*)dB1(t),dI(t)=βS(t)∫0h‍f(s)I(t-s)ds-(μ2+λ)I(t)+σ2(I(t)-I*)dB2(t),dR(t)=λI(t)-μ3R(t)+σ3(R(t)-R*)dB3(t),
where S*, I*, and R*are the positive points of equilibrium for the corresponding deterministic system and σi are constants, Bi(t) are independent from other standard Wiener processes (i=1,2,3). In their paper, they proved the above stochastic system was stable in probability by using the general method of Lyapunov functionals construction and obtained the stability condition immediately in terms of the parameters of the system under consideration.Another different approach to include stochastic perturbations in a biological model was considered by Imhof and Walcher in [11]. They introduced and analyzed a variant of the deterministic single-substrate chemostat model. Also they found all species were persistent in the multiple species by using a comparison principle. Subsequently, they turned to modeling the influence of random fluctuations on the deterministic chemostat model. They set up the following stochastic chemostat model:
(5)dX0=(r-δX0-a1(X0,X1)-a2(X0,X2))dt+σ0X0dB0(t),dX1=(a1(X0,X1)-s1(X1,X2))dt+σ1X1dB1(t),dX2=(a2(X0,X2)-s2(X1,X2))dt+σ2X2dB1(t),
where σ0, σ1, σ2≥0 and Bi(t)(i=1,2,3) are independent Brownian motions. By analyzing this stochastic differential equation model, they proved that the stochastic model led to extinction even though the deterministic counterpart predicts persistence.In this paper, our approach to include stochastic perturbation is analogous to that of Imhof and Walcher [11]. We are devoted to studying the following stochastic system:
(6)dS(t)=[(1-q)A-βSI-(μ+p)S+γI+εV]dt+σ1SdB1(t),dI(t)=[βS-(μ+γ+α)]Idt+σ2IdB2(t),dV(t)=[qA+pS-(μ+ε)V]dt+σ3VdB3(t),
where Bi(t)(i=1,2,3) are independent Brownian motions and σi(i=1,2,3) are their intensities.This paper is organized as follows. In Section2, we show there is a unique positive solution of system (6) by the way mentioned in [20, 23]. In Section 3, when R0≤1, we derive that the solution of the system (6) oscillates around the disease-free equilibrium P0 of system (1). In Section 4, when R0>1, although the solution of system (6) does not converge to the endemic proportion equilibrium P* of system (1), we conclude there is a unique stationary distribution for system (6), and it has ergodic property under some conditions. In Section 5, we make simulations to confirm our analytical results. Finally, in order to be self-contained, we have the appendix.


## 2. Existence and Uniqueness of Positive Solution

To investigate the dynamical behavior, the first concern is whether the solution has a global existence. Moreover, for a population dynamics model, whether the value is nonnegative is also considered. Hence in this section, we first show that the solution of system (6) is global and nonnegative.Theorem 1.
There is a unique solutionY(t)=(S(t),I(t),V(t)) of system (6) on t⩾0 for any initial value Y(0)=(S(0),I(0),V(0))∈ℝ+3, and the solution will remain in ℝ+3 with probability 1, namely, Y(t)∈ℝ+3  for all t⩾0 almost surely.Proof.
Our proof is motivated by the works of Mao et al. [17]. Since the coefficients of (6) are locally Lipschitz continuous for any given initial value (S(0),I(0),V(0))∈ℝ+3, there is a unique local solution (S(t),I(t),V(t)) on t∈[0,τe), where τe is the explosion time (see [20]). To show that this solution is global, we need to show that τe=∞ a.s. Let k0≥0 be sufficiently large so that S(0), I(0), and V(0) all lie within the interval [1/k0,k0]. For each integer k⩾k0, define the stopping time:
(7)τk=inf{1kt∈[0,τe):min{S(t),I(t),V(t)}≤1kormax{S(t),I(t),V(t)}≥k},
where throughout this paper, we set inf∅=∞ (as usual ∅ denotes the empty set). According to the definition, τk is increasing as k→∞.  Set τ∞=limk→∞τk, whence τ∞⩽τe a.s. If we can show that τ∞=∞ a.s., then τe=∞ and (S(t),I(t),V(t))∈ℝ+3 a.s. for all t⩾0. In other words, to complete the proof, all we need to show is that τ∞=∞ a.s. If this statement is false, then there exist a pair of constants T>0 and ε∈(0,1) such that
(8)P{τ∞⩽T}>ε.
Hence there is an integer k1⩾k0 such that
(9)P{τk⩽T}⩾ε∀k⩾k1.
Define a C2-function W:ℝ+3→ℝ- (see the appendix) by
(10)W=(S-a-alogSa)+(I-1-logI)+(V-1-logV).
The nonnegativity of this function can be seen from u-1-logu⩾0, for all u>0. Let k⩾k0 and T>0 be arbitrary. Applying the Itô formula, we obtain
(11)dW=(1-aS)[(1-q)A-βSI-(μ+p)S+γI+εV]dt+(1-aS)σ1SdB1(t)+12aσ12dt+(1-1I)[βS-(μ+γ+γ)]Idt+(1-1I)σ2IdB2(t)+12σ22dt+(1-1V)[qA+pS-(μ+ε)V]dt+(1-1V)σ3VdB3(t)+12σ32dt:=LWdt+σ1(S-a)dB1(t)+σ2(I-1)dB2(t)+σ3(V-1)dB3(t),
where LW:ℝ+3→ℝ+ is defined by
(12)LW(S,I,V)=A-μS-(1-q)AaS+aβI+a(μ+p)-aγIS-aεVS+12aσ12-(μ+α)I-βS+(μ+α+γ)+12σ22-μV-qAV-pSV+(μ+ε)+12σ32⩽A+[aβ-(μ+α)]I+a(μ+p)+12aσ12+12σ22+12σ32+(μ+α+γ)+(μ+ε).
Let a=(μ+α)/β, s.t. aβ-(μ+α)=0. Substituting this into (12), we get
(13)LW(S,I,V)⩽A+(μ+α)(μ+p)β+μ+α2βσ12h.+(μ+γ+α)+12σ22h.+(μ+ε)+12σ32:=K.
The remainder of the proof follows that in Mao et al. [17].


## 3. Asymptotic Behavior aroundP0

It is clear thatP0=(S0,0,V0) is the disease-free equilibrium of system (1). If R0≤1, then P0 is globally stable, which means the disease will die out after some period of time. Hence, it is interesting to study the disease-free equilibrium for controlling infectious disease. But there is no disease-free equilibrium in system (6). In this section, we show the average oscillation around P0 in time to exhibit whether the disease will die out.Theorem 2.
IfR0⩽1,  α2<4μ(μ+α)(1+2μ/ε), and the following condition is satisfied
(14)μ(2(μ+p)ε+1)>(c2+1)σ12,(μ+α)(1+2με)-α4ρ>(c2+1)σ22,(μ-ρα)>σ32,
then the solution Y(t) of system (6) with initial value Y(0)∈ℝ+3 has the property
(15)limsupt→∞1tE∫0t‍{[μ(2(μ+p)ε+1)-(c2+1)σ12]hhhhhhhhhhhhh×(S(r)-S0)2hhhhhhhhhhhhh+[α4ρ-(c2+1)σ22(μ+α)(1+2με)hhhhhhhhhhhhhhhh-α4ρ-(c2+1)σ22(μ+α)(1+2με)]I(r)2hhhhhhhhhhhhh+[(μ-ρα)-σ32](V(r)-V0)22(μ+p)ε}dshhh⩽(c2+1)σ12S02+σ32V02,
where P0=(S0,0,V0) is the disease-free equilibrium of system (1), c2=2μ/ε, and ρ is positive constant that satisfied α/4(μ+α)(1+2μ/ε)<ρ<μ/α.Proof.
Letx=S-S0, y=I, and z=V-V0. Then
(16)dx=[-β(x+S0)y-(μ+p)x+γy+εz]dt+σ1(x+S0)dB1(t),dy={βxy+[βS0-(μ+γ+α)]y}dt+σ2ydB2(t),dz=[px-(μ+ε)z]dt+σ3(z+V0)dB3(t).
Define a C2-function W:ℝ+3→ℝ+ by
(17)W(x,y,z)=c1y+c22(x+y)2+12(x+y+z)2:=c1W1(x,y,z)+c2W2(x,y,z)+W3(x,y,z),
where c1=1/β[2μ(2μ+p+α)/ε+(2μ+α)]. By Itô's formula, we compute
(18)dW1={βxy+[βS0-(μ+γ+α)]y}dt+σ2ydB2(t)=[βxy+(μ+γ+α)(R0-1)y]dt+σ2ydB2(t)=LW1dt+σ2ydB2(t),dW2=(x+y){[12-(μ+p)x-(μ+α)yHHHHHHH.+εz+12σ12(x+S0)2+12σ22y2]dtHHHHHH..+σ1(x+S0)dB1(t)+σ2ydB2(t)12}=LW2dt+(x+y)[σ1(x+S0)dB1(t)hhhhhhhhhhhhhhh.+σ2ydB2(t)],dW3=(x+y+z){[12-μx-(μ+α)y-μzhhhhhhhhhhl+12σ12(x+S0)2+12σ22y2hhhhhhhhhhl+12σ32(z+V0)2]dthhhhhhhhhl+σ1(x+S0)dB1(t)+σ2yB2(t)hhhhhhhhhl+σ3(z+V0)dB3(t)12}=LW3dt+(x+y+z)×[σ1(x+S0)dB1(t)+σ2yB2(t)hhh+σ3(z+V0)dB3(t)].
If R0⩽1, then
(19)LW=c1LW1+c2LW2+LW3⩽[c1β-c2(2μ+p+α)-(2μ+α)]xy-[c2(μ+p)+μ]x2-(c2+1)(μ+α)y2-μz2+(c2ε-2μ)xz+[c2ε-(2μ+α)]yz+12(c2+1)σ12(x+S0)2+12(c2+1)σ22y2+12σ32(z+V0)2.
Note that
(20)c1β-c2(2μ+p+α)-(2μ+α)=0,c2ε-2μ=0.
Therefore, we can obtain
(21)LW⩽-μ[2(μ+p)ε+1]x2-(1+2με)(μ+α)y2-μz2-αyz+12(c2+1)σ12(x+S0)2+12(c2+1)σ22y2+12σ32(z+V0)2⩽-{μ[2(μ+p)ε+1]-(c2+1)σ12}x2-[(μ+α)(1+2με)-α4ρ-(c2+1)σ22(μ+α)(1+2με)-α4ρ]y2-[(μ-ρα)-σ32]z2+(c2+1)σ12S02+σ32V02,
where ρ is positive constant to be specified later and Young’s inequality is used. If
(22)α2<4μ(μ+α)(1+2με),
we choose ρ such that α/4(μ+α)(1+2μ/ε)<ρ<μ/α. This implies
(23)(μ+α)(1+2με)-α4ρ>0,μ-ρα>0.
Then we have
(24)dW⩽-{[μ(2(μ+p)ε+1)-(c2+1)σ122(μ+p)ε]x2+[(μ+α)(1+2με)-α4ρ-(c2+1)σ22]y2+[(μ-ρα)-σ32]z2-(c2+1)σ12S02-σ32V02(2(μ+p)ε+1)}dt+σ1(x+S0)[(c2+1)(x+y)+z]dB1(t)+σ2y[c1+(c2+1)(x+y)+z]dB2(t)+σ3(x+y+z)(z+V0)dB3(t).
Integrating this from 0 to t and taking the expectation, we have
(25)EW(t)-W(0)⩽-E∫0t‍{[μ(2(μ+p)ε+1)-(c2+1)σ12]x2(r)hhhhhhhh+[(μ+α)(1+2με)-α4ρhhhhhhhhhHHHHHHH-(c2+1)σ22(μ+α)(1+2με)]y2(r)hhhhhhhh+[(μ-ρα)-σ32]z2(r)hhhhhhhh-(c2+1)σ12S02-σ32V022(μ+p)ε}dr.
Hence,
(26)limsupt→∞1tE∫0t‍{[μ(2(μ+p)ε+1)hhhhhhhhhhhhh-(c2+1)σ12(2(μ+p)ε+1)]x2(r)hhhhhhhhhhhh+[(μ+α)(1+2με)-α4ρhhhhhhhhhhhhhhh-(c2+1)σ22(μ+α)(1+2με)-α4ρ]y2(r)hhhhhhhhhhhh+[(μ-ρα)-σ32]z2(r)2(μ+p)ε}dr⩽(c2+1)σ12S02+σ32V02.
Consequently,
(27)limsupt→∞1tE∫0t‍{[μ(2(μ+p)ε+1)hhhhhhhhhhhhh-(c2+1)σ122(μ+p)ε]hhhhhhhhhhhh×(S(r)-S0)2hhhhhhhhhhhh+[α4ρ-(c2+1)σ22(μ+α)(1+2με)hhhhhhhhhhhhhhh-α4ρ-(c2+1)σ22(1+2με)]I(r)2hhhhhhhhhhhh+[(μ-ρα)-σ32]hhhhhhhhhhhh×(V(r)-V0)22(μ+p)ε}dr⩽(c2+1)σ12S02+σ32V02.
This finishes the proof of Theorem 2.Remark 3.
From Theorem2, we show the solution of system (6) will oscillate around the disease-free equilibrium of system (1) under some conditions and the disturbance intensity is proportional to the intensity of the white noise. In a biological view, as the intensity of stochastic perturbations is small, we consider the disease will die out.
Besides, ifσ1=0 and σ3=0, then P0 is also the disease-free equilibrium of system (6). From the proof of Theorem 2, we get
(28)LW⩽-{μ[2(μ+p)ε+1]-(c2+1)σ122(μ+p)ε}x2-[(μ+α)(1+2με)-α4ρ-(c2+1)σ22]y2-[(μ-ρα)-σ32]z2.
Thus, the solution of system (6) is stochastically asymptotically stable in the large if (μ+2)(1+2μ/ε)-α/4ρ>(c2+1)σ22.


## 4. Ergodicity of System (6)

In studying epidemic dynamical system, we are also interested in when the disease will prevail and persist in a population. In the deterministic models, the problem is solved by showing that the endemic equilibrium of corresponding model is a global attractor or is globally asymptotically stable. But there is no endemic equilibrium in system (6). In this section, based on the theory of Has’minskiĭ [26] (see the appendix), we show that there is a stationary distribution subjected to some conditions on R0 and the parameters of the model which reveals that the disease will prevail also.System (6) can be written as a form of system (A.8) (see the appendix).Consider(29)d(SIV)=((1-q)A-βS(t)I(t)-(μ+p)S(t)+γI(t)+εV(t)βS(t)I(t)-(μ+γ+a)I(t)qA+pS(t)-(μ+ε)V(t))dt+(σ1S00)dB1(t)+(0σ2I0)dB2(t)+(00σ3V)dB3(t),
and the diffusion matrix is A=diag(σ12S2,σ22I2,σ32V2).Remark 4.
Theorem1 shows that there exists a unique positive solution Y(t)=(S(t),I(t),V(t)) of system (6). Besides, from the proof of Theorem 1, we obtain LV⩽K. Define V~=V+K, then LV~≤V~, and it is clear that V~k=infX∈R+3∖CkV~(Y)→∞ as k→∞, where Ck=(1/k,k)×(1/k,k)×(1/k,k). Hence, by Remark 2 of Theorem 4.1 of Has’minskiĭ (1980, page 86) [26], we obtain that the solution Y(t) is a homogeneous Markov process in R+3.Theorem 5.
IfR0>1, α2<4μ(μ+α)(1+2μ/ε), and
(30)0<σ12<μc2+1(2(μ+p)ε+1),0<σ22<1c2+1[(μ+α)(1+2με)-α4ρ],0<σ32<μ-ρα,
such that
(31)δ<min{[μ(2(μ+p)ε+1)hhhhhhhh-(c2+1)σ122(μ+p)ε](S*)2,hhh1hhh[(c2+1)σ22(μ+α)(1+2με)-α4ρhhhhhhhhh-(c2+1)σ22(1+2με)](I*)2,hhh1hhh(μ-ρα-σ32)(V*)22(μ+p)ε},
then, for any initial value Y(0)∈R+3, there is a stationary distribution μ(·) for system (6) and it has an ergodic property, where δ=σ22I*/2+(c2+1)σ12(S*)2+(c2+1)σ22(I*)2+σ32(V*)2, P*=(S*,I*,V*) is the endemic proportion equilibrium of system (1), and c2 and ρ are the same constants as in Theorem 2.Proof.
IfR0>1, there is an endemic proportion equilibrium P*=(S*,I*,V*) of system (1); then
(32)(1-q)A=βS*I*+(μ+p)S*-γI*-εV*,βS*I*=(μ+γ+α)I*,qA+pS*=(μ+ε)V*.
Define aC2-function W:ℝ+3→ℝ+ by
(33)W(S,I,V)=c1(I-I*-I*logII*)+c22[(S-S*)+(I-I*)]2+12[(S-S*)+(I-I*)+(V-V*)]2:=c1W1+c2W2+W3,
where c1=1/β[2μ(2μ+p+α)/ε+(2μ+α)]. The nonnegativity of the function can be observed from u-1-logu≥0 on u>0. Let L be the generating operator of system (6). Then we get
(34)LW1=β(S-S*)(I-I*)+12σ22I*,LW2=-(μ+p)(S-S*)2-(μ+α)(I-I*)2-[(μ+p)+(μ+α)](S-S*)(I-I*)+ε(S-S*)(V-V*)+ε(I-I*)(V-V*)+12(σ12S2+σ22I2),LW3=-μ(S-S*)2-(μ+α)(I-I*)2hl-μ(V-V*)2-(2μ+α)(S-S*)(I-I*)hl-2μ(S-S*)(V-V*)hl-(2μ+α)(I-I*)(V-V*)hl+12(σ12S2+σ22I2+σ32V2).
Then
(35)LW=c1LW1+c2LW2+LW3=[c1β-c2(2μ+p+α)-(2μ+α)]×(S-S*)(I-I*)-[c2(μ+p)+μ](S-S*)2-(c2+1)(μ+α)(I-I*)2-μ(V-V*)2+(c2ε-2μ)(S-S*)(V-V*)+[c2ε-(2μ+α)](I-I*)(V-V*)+12σ22I*+c22(σ12S2+σ22I2)+12(σ12S2+σ22I2+σ32V2).
Note that
(36)c1β-c2(2μ+p+α)-(2μ+α)=0,c2ε-2μ=0.
Then
(37)LW=-μ[2(μ+p)ε+1](S-S*)2-(1+2με)(μ+α)(I-I*)2-μ(V-V*)2-α(I-I*)(V-V*)+12σ22I*+c2+12σ12S2+c2+12σ22I2+σ32V2⩽-[μ(2(μ+p)ε+1)-(c2+1)σ12](S-S*)2-[α4ρ-(c2+1)σ22(μ+α)(1+2με)hhhh-α4ρ-(c2+1)σ22(1+2με)](I-I*)2-(μ-ρα-σ32)(V-V*)2+12σ22I*+(c2+1)σ12(S*)2+(c2+1)σ22(I*)2+σ32(V*)2∶=-[μ(2(μ+p)ε+1)-(c2+1)σ12](S-S*)2-[α4ρ-(c2+1)σ22(μ+α)(1+2με)hhhh-α4ρ-(c2+1)σ22](I-I*)2-(μ-ρα-σ32)(V-V*)2+δ,
where Young’s inequality is used and
(38)δ=12σ22I*+(c2+1)σ12(S*)2+(c2+1)σ22(I*)2+σ32(V*)2.
If α2<4μ(μ+α)(1+2μ/ε), choose ρ such that
(39)(μ+α)(1+2με)-α4ρ>0,μ-ρα>0.
Note that
(40)δ<min{[μ(2(μ+p)ε+1)hhhhhhh-(c2+1)σ12(2(μ+p)ε+1)](S*)2,hhhhhh[(μ+α)(1+2με)-α4ρhhhhhhh-(c2+1)σ22(1+2με)-α4ρ](I*)2,hhhhhh[(μ-ρα-σ32)](V*)2(2(μ+p)ε+1)}.
Thus the ellipsoid
(41)-[μ(2(μ+p)ε+1)-(c2+1)σ12](S-S*)2-[α4ρ-(c2+1)σ22(μ+α)(1+2με)hhh-α4ρ-(c2+1)σ22(1+2με)](I-I*)2-(μ-ρα-σ32)(V-V*)2+δ=0
lies entirely in R+3. We can take U to be a neighborhood of the ellipsoid with U-⊆El=R+3, so, for x∈R+3∖U, LV≤-M (M is a positive constant), which implies that condition (B.2) in Lemma A.2 is satisfied. Hence, the solution Y(t)=(S(t),I(t),V(t)) is recurrent in the domain U, which, together with Lemma A.4 and Remark 4, implies that Y(t)=(S(t),I(t),V(t)) is recurrent in any bounded domain D⊂R+3. Besides, for all D, there is
(42)M=min{σ12S2,σ22I2,σ32V2}>0,
such that
(43)∑i,j=13‍aijξiξj=σ12S2ξ12+σ22I2ξ22+σ32V2ξ32≥M∥ξ∥2,
for all Y∈D-, ξ∈R3, which implies that condition (B.1) is also satisfied. Therefore, the stochastic system (6) has a stationary distribution μ(·), and it is ergodic.


## 5. Numerical Simulations

In order to confirm the results above, we numerically simulate the solution of system (6) with given initial value and the parameters. Using Milstein’s higher order method mentioned in [27], we get the discretization equation of system (6):
(44)Sk+1=Sk+((1-q)A-βSkIkhhhhhhh-(μ+p)Sk+γIk+εVk(1-q))Δt+σ1SkΔtξ1,k+σ122SkΔt(ξ1,k2-1),Ik+1=Ik+(βSkIk-(μ+γ+α)Ik)Δt+σ2IkΔtξ2,k+σ222IkΔt(ξ2,k2-1),Vk+1=Vk+(qA+pSk-(μ+ε)Vk)Δt+σ3VkΔtξ3,k+σ322VkΔt(ξ3,k2-1),
where ξ1,k, ξ2,k, and ξ3,k, k=1,2,…,n, are the independent Gaussian random variables N(0,1). Choosing suitable parameters in the system, by Matlab we get the simulation figures with initial value (S(0),I(0),V(0))=(0.8,0.4,0.5) and time step Δt=0.001.Example 6.
Throughout the paper we will assume that the unit of time is one day, and the population sizes are measured in units of 1 million. Choose the parameters in system (6) as follows: A=0.2, q=0.5, β=0.5, μ=0.2, p=0.3, γ=0.2, ε=0.1, α=0.3, σ1=0.05, σ2=0.1, and σ3=0.05. Note that
(45)R0=Aβ(μ(1-q)+ε)μ(μ+γ+α)(μ+ε+p)=0.4<1,α2=0.09<4μ(μ+α)(1+2με)=2,
and σ1, σ2, and σ3 satisfy condition (14), then by Theorem 2, we show the solution of system (6) will oscillate around the disease-free equilibrium P0(S0,I0,V0) of system (1) in time which is globally asymptotically stable. Besides Theorem 2 tells us that the difference between the perturbed solution and P0 is only related with white noises σ1 and σ3. Numerical simulations are shown in Figure 1.Computer simulation of the pathsS(t), I(t), and V(t) for the SDE SIS model (6) with vaccination and its corresponding deterministic model (1) for parameter values A=0.2, q=0.5, β=0.5, μ=0.2, p=0.3, γ=0.2, ε=0.1, α=0.3, σ1=0.05, σ2=0.1, and σ3=0.05.
(a)
(b)Example 7.
In Figures2 and 3, we choose the parameters in system (6) as follows: A=0.9, q=0.1, β=0.8, μ=0.4, p=0.1, γ=0.1, ε=0.2, α=0.1, σ1=0.1, σ2=0.1, and σ3=0.1. Note that R0=2.486>1 and α2<4μ(μ+α)(1+2μ/ε)=4; then the equilibrium P*(S*,I*,V*) of system (1) is globally asymptotically stable. But the white noise may make system (1) appear as different phenomena. In detail, the conditions (30) and (31) are also satisfied. Therefore, as Theorem 5 stated, there is a stationary distribution (see the histograms in the middle in Figure 2). In addition, the left pictures in Figure 2 show that the solution of system (6) has fluctuating in a small neighborhood. Moreover, from Figure 3, we find that 95% or more of the population distribution lies within a neighborhood, which can be imagined a circular or elliptic region centered at P*(S*,I*,V*) (see the red point in Figure 3). All of these imply system (6) has stochastic stability.
Testing these data for normality, all tests used were highly significant, conclusively rejecting normality in all cases. This is not surprising in view of the very large sample sizes (10,000), as even moderate deviations from the tested distribution will be significant; however the normal quantile-quantile plots (see theQQ plots in the right in Figure 2) suggest that these data are not far from being normally distributed for smaller values of σi (i=1,2,3).The solution of the stochastic system (6), its histogram, and normal quantile-quantile plots for parameter values A=0.9, q=0.1, β=0.8, μ=0.4, p=0.1, γ=0.1, ε=0.2, α=0.1, σ1=0.1, σ2=0.1, and σ3=0.1.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 3
Population distribution around pointP*(S*,I*,V*)=(0.75,0.98,0.275) (red point) corresponding to Figure 2.


## 6. Conclusions

As most real world problems are not deterministic, incorporating stochastic effects into the model give us a more realistic way of modeling epidemic models. In this paper, we have considered stochastic SIS epidemic models with vaccination. We first proved the positivity of the solutions. When the perturbation and the disease-related death rate are small, we illustrated the dynamical behavior of the SDE model according toR0≤1 or R0>1. When R0≤1, we proved that the asymptotic behavior around the disease-free equilibrium of the deterministic model. When R0>1, we also proved that the SDE model has the ergodic property as the fluctuation is small, where the positive solution converges weakly to the unique stationary distribution. Simulations are also carried out to verify our analytical results.



---

*Source: 517439-2013-09-03.xml*