# On a Stochastic SEIS Model with Treatment Rate of Latent Population

Authors: Shujing Gao; Yanfei Dai; Yan Zhang; Yujiang Liu
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

The asymptotic dynamics of a stochastic SEIS epidemic model with treatment rate of latent population is investigated. First, we show that the system provides a unique positive global solution starting from the positive initial value. Then, the long-term asymptotic behavior of the model is studied: ifR
0, which is called the basic reproduction number of the corresponding deterministic model, is not more than unity, the solution of the model is oscillating around the disease-free equilibrium of the corresponding deterministic system, whereas if R
0 is larger than unity, we show how the solution spirals around the endemic equilibrium of deterministic system under certain parametric restrictions. Finally, numerical simulations are carried out to support our theoretical findings.


---

## Body

## 1. Introduction

Mathematical models describing the population dynamics of infectious diseases have made a significant progress in better understanding of disease transmissions and behavior of epidemics. There have been a large number of works on the dynamics of epidemic models described by ordinary differential equations [1–9] and the references cited therein. In most of the literatures, the incubation period is often to be ignored usually. However, for some diseases, such as tuberculosis, schistosomiasis, measles, and AIDS, once a susceptible individual adequate contact with an infective, it becomes exposed, that is, infected but not infective. This individual remains in the exposed class for a certain latent period before becoming infective [10–13]. Particularly, Fan and Li in [12] established a class of SEIS epidemic model that incorporates constant recruitment, disease-caused death, and disease latency as follows:
(1)
S
′
=
A
-
λ
I
S
-
d
S
+
γ
I
,
E
′
=
λ
I
S
-
(
d
+
ϵ
)
E
,
I
′
=
ϵ
E
-
(
d
+
γ
+
α
)
I
,
where the meaning of the parameters can be found in the literature [12]. The author obtained that the global dynamics is completely determined by the basic reproduction number R
0
=
A
λ
ϵ
/
d
(
d
+
ϵ
)
(
d
+
γ
+
α
). If R
0
≤
1, then the disease-free equilibrium P
0
=
(
A
/
d
,
0,0
) is the only equilibrium and it is globally asymptotically stable, implying that the disease dies out. If R
0
>
1, then P
0 becomes unstable and there exists a unique endemic equilibrium P
*
=
(
S
*
,
E
*
,
I
*
) with S
*
=
(
d
+
ϵ
)
(
d
+
γ
+
α
)
/
λ
ϵ, I
*
=
(
A
-
d
S
*
)
/
(
λ
S
*
-
γ
), and E
*
=
(
(
d
+
γ
+
α
)
/
ϵ
)
I
*, and P
* is globally asymptotically stable in the interior of the feasible region, meaning that the disease persists at the endemic equilibrium. Furthermore, Castillo-Chavez and Feng in [14] considered an SEIS model which described the transmission of tuberculosis with standard incidence ratio and treatment rates of latent individuals. They show that there is a global stability switch from the disease-free equilibrium to the positive endemic equilibrium when the basic reproductive number passes through the critical value 1. Xu in [15] studied an SEIS epidemiological model with a saturation incidence rate and a time delay representing the latent period of the disease. He obtained the basic reproduction number R
0 and the global asymptotical stability of disease-free equilibrium and endemic equilibrium.In this paper, we consider a class of deterministic SEIS model that incorporates bilinear incidence rate, disease-caused death rate, disease latency, and treatment rates of latent and infectious individuals as follows:(2)
S
˙
(
t
)
=
Λ
-
β
S
(
t
)
I
(
t
)
-
μ
S
(
t
)
+
δ
I
(
t
)
+
γ
E
(
t
)
,
E
˙
(
t
)
=
β
S
(
t
)
I
(
t
)
-
(
μ
+
ϵ
+
γ
)
E
(
t
)
,
I
˙
(
t
)
=
ϵ
E
(
t
)
-
(
μ
+
δ
+
d
)
I
(
t
)
,
where Λ is the recruitment rate of individuals (including newborns and immigrants) in the susceptible population, the natural death rate is assumed to be the same constant μ for all hosts, infectious hosts suffer an extra disease-related death with constant rate d, γ and δ are the per capita treatment rates of latent and infectious individuals, respectively, ϵ is the rate at which the exposed individuals become infective so that 1
/
ϵ is the mean latent period, and the incidence term is taken as the bilinear mass-action form β
S
(
t
)
I
(
t
). Assume that all parameters of (2) are positive constants. Using a similar argument as in [12], we can easily get the following results.(i)
Model (2) always has a disease-free equilibrium F
0
=
(
Λ
/
μ
,
0,0
). And the basic reproduction number R
0
=
β
Λ
ϵ
/
μ
(
μ
+
ϵ
+
γ
)
(
μ
+
δ
+
d
) is the threshold of the system for an epidemic to occur.
(ii)
WhenR
0
≤
1, then the disease-free equilibrium F
0 is the only equilibrium and globally asymptotically stable, which implies that the disease will disappear eventually.
(iii)
WhenR
0
>
1, then F
0 becomes unstable and there exists a unique endemic equilibrium F
*
=
(
S
*
,
E
*
,
I
*
), where S
*
=
(
μ
+
ϵ
+
γ
)
(
μ
+
δ
+
d
)
/
β
ϵ, I
*
=
μ
(
μ
+
ϵ
+
γ
)
(
μ
+
δ
+
d
)
(
R
0
-
1
)
/
β
[
μ
(
μ
+
δ
+
d
)
+
ϵ
(
μ
+
d
)
], and E
*
=
(
(
μ
+
δ
+
d
)
/
ϵ
)
I
*. Furthermore, F
* is globally asymptotically stable in the interior of the feasible region, which means that the disease will always prevail and persist in the population.However, in the real world, epidemic dynamics is inevitably affected by the environmental noise which is an important component in the epidemic systems. As a matter of fact, the epidemic models are often subject to environmental noise; that is, due to environmental fluctuations, parameters involved in epidemic models are not absolute constants, and they may fluctuate around some average values. So, inclusion of random perturbations in such models makes them more realistic in comparison to their deterministic model. In recent years, epidemic models under environmental noise described by stochastic different equations have been studied by many researchers [16–23]. To the best of our knowledge, there are few papers to deal with the stochastic epidemic model with latent individuals [22], as it is difficult to choose appropriate Lyapunov functions. Yang et al. in [22] include stochastic perturbations into SIR and SEIR epidemic models with saturated incidence and investigate their dynamics according to the basic reproduction number R
0.In addition, from a biological and mathematical perspective, there are different possible approaches which result in different effects on the population system to include random effects in the model. There are usually four approaches to be mentioned [17]. The first one is through time Markov chain model to consider environment noise [18]. The second is with parameters perturbation [19]. The third one is the environmental noise that is proportional to the variables [20], and the last one is to robust the positive equilibria of deterministic models [21]. In this paper, we will consider a stochastic counterpart of model (2) by the third approach. That is, the stochastic perturbation is assumed to be of a white noise type which is directly proportional to S
(
t
), E
(
t
), and I
(
t
) and influenced by the S
˙
(
t
), E
˙
(
t
), and I
˙
(
t
) in model (2). By this way, a reasonable stochastic analogue of system (2) is given by
(3)
d
S
(
t
)
=
(
Λ
-
β
S
(
t
)
I
(
t
)
-
μ
S
(
t
)
+
δ
I
(
t
)
+
γ
E
(
t
)
)
d
t
+
σ
1
S
(
t
)
d
B
1
(
t
)
,
d
E
(
t
)
=
(
β
S
(
t
)
I
(
t
)
-
(
μ
+
ϵ
+
γ
)
E
(
t
)
)
d
t
+
σ
2
E
(
t
)
d
B
2
(
t
)
,
d
I
(
t
)
=
(
ϵ
E
(
t
)
-
(
μ
+
δ
+
d
)
I
(
t
)
)
d
t
+
σ
3
I
(
t
)
d
B
3
(
t
)
,
where B
i
(
t
) are mutually independent Brownian motions and σ
i represent the noise intensities of B
i
(
t
), i
=
1,2
,
3.Since system (3) is constructed by adding stochastic perturbation in a deterministic system (2), it seems reasonable to investigate whether there are similar properties as in system (2). But there is neither a disease-free equilibrium F
0 nor an endemic equilibrium F
* for system (3). Hence, in order to show the stability to some extent, we discuss the behavior around F
0 and F
*, respectively. Our main purpose is to investigate how the solution of system (3) spirals around the disease-free equilibrium and endemic equilibrium of the corresponding deterministic system under some conditions. We first show the existence and uniqueness of a global positive solution of model (3). Then the main results will be seen in Sections 3 and 4. Finally, numerical simulations are present in Section 5 to illustrate our results.Throughout this paper, unless otherwise specified, we let(
Ω
,
F
,
{
F
t
}
t
≥
0
,
P
) be a complete probability space with a filtration satisfying the usual conditions (i.e., it is right continuous and F
0 contains allP-null sets). Let B
i
(
t
) (i
=
1,2
,
3) denote the independent standard Brownian motions defined on this probability space. We also denote R
+
3
=
{
x
∈
R
3
:
x
i
>
0 for all 1
≤
i
≤
3
} and x
(
t
)
=
(
S
(
t
)
,
E
(
t
)
,
I
(
t
)
)
T.Here we show the following auxiliary statements which are introduced in [24].Consider them-dimensional stochastic differential equation
(4)
d
x
(
t
)
=
f
(
x
(
t
)
,
t
)
d
t
+
g
(
x
(
t
)
,
t
)
d
B
(
t
)
,
for
t
≥
t
0
.
Denote by C
2,1
(
R
m
×
[
t
0
,
∞
]
;
R
) the family of all nonnegative functions V
(
x
,
t
) defined on R
m
×
[
t
0
,
∞
], such that they are continuously twice differentiable in x and once in t. Define the differential operator L associated with (4) by
(5)
L
=
∂
∂
t
+
∑
i
=
1
m
‍
f
i
(
x
,
t
)
∂
∂
x
i
+
1
2
∑
i
=
1
m
‍
[
g
T
(
x
,
t
)
g
(
x
,
t
)
]
i
j
∂
2
∂
x
i
∂
x
i
.
If L acts on a function V
∈
C
2,1
(
R
m
×
[
t
0
,
∞
]
;
R
), then
(6)
L
V
(
x
,
t
)
=
V
t
(
x
,
t
)
+
V
x
(
x
,
t
)
f
(
x
,
t
)
+
1
2
trace
[
g
T
(
x
,
t
)
V
x
x
(
x
,
t
)
g
(
x
,
t
)
]
,
where V
t
=
∂
V
/
∂
t
,
V
x
=
(
∂
V
/
∂
x
1
,
…
,
∂
V
/
∂
x
m
) and V
x
x
=
(
∂
2
V
/
∂
x
i
∂
x
j
)
m
×
m. By Itô’s formula, then
(7)
d
V
(
x
(
t
)
,
t
)
=
L
V
(
x
(
t
)
,
t
)
d
t
+
V
x
(
x
(
t
)
,
t
)
g
(
x
(
t
)
,
t
)
d
B
(
t
)
.


## 2. Existence and Uniqueness of the Positive Solution

In order to investigate the dynamical behavior, the first problem considered is the global existence of the solution. SinceS
(
t
), E
(
t
), and I
(
t
) in model (2) are the sizes of the susceptible individuals, latent individuals, and infected individuals at time t, respectively, they should be nonnegative. Therefore, we are only interested in positive solutions. In order to obtain 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 [24]. However, the coefficients of model (3) do not satisfy the linear growth condition, though they are locally Lipschitz continuous. Hence, the solution of model (3) may explode at a finite time. In what follows, we will prove that the solution of model (3) is positive and global.Theorem 1.
There is a unique positive solution (S
(
t
), E
(
t
), I
(
t
)) of system (3) on t
≥
0 for any given initial value (
S
(
0
)
,
E
(
0
)
,
I
(
0
)
)
∈
R
+
3, and the solution will remain in R
+
3 with probability 1; that is, (
S
(
t
)
,
E
(
t
)
,
I
(
t
)
)
∈
R
+
3 for all t
≥
0 almost surely.Proof.
Since the coefficients of system (3) are locally Lipschitz continuous, then for any initial value (
S
(
0
)
,
E
(
0
)
,
I
(
0
)
)
∈
R
+
3, there exists a unique local solution (S
(
t
), E
(
t
), I
(
t
)) on t
∈
[
0
,
τ
e
), where τ
e is the explosion time [23]. Hence it suffices to prove that the unique local solution of system (3) is global and positive. To show this solution is global, we need to show that τ
e
=
+
∞ a.s. Let m
0
>
0 be sufficiently large for S
(
0
)
∈
[
1
/
m
0
,
m
0
], E
(
0
)
∈
[
1
/
m
0
,
m
0
], and I
(
0
)
∈
[
1
/
m
0
,
m
0
]. For each integer m
≥
m
0, define the stopping time
(8)
τ
m
=
inf
⁡
{
t
∈
[
0
,
τ
e
)
:
S
(
t
)
∉
[
1
m
,
m
]
or
E
(
t
)
∉
[
1
m
,
m
]
or
I
(
t
)
∉
[
1
m
,
m
]
}
,
where throughout this paper we set inf
⁡
⌀
=
∞ (as usual ⌀ denotes the empty set). Clearly, τ
m is increasing as m
→
∞. Set τ
∞
=
lim
⁡
m
→
∞
τ
m; hence, τ
∞
≤
τ
e a.s. If it holds that τ
∞
=
∞ a.s., then τ
e
=
∞ a.s. and (
S
(
t
)
,
E
(
t
)
,
I
(
t
)
)
∈
R
+
3 a.s. for t
≥
0. In other words, to complete the proof we need to show that τ
∞
=
∞ a.s. If this statement is false, then there exist a pair of constants T
>
0 and ε
∈
(
0,1
) such that
(9)
P
{
τ
∞
≤
T
}
>
ε
.
Hence, there is an integer m
1
≥
m
0 such that
(10)
P
{
τ
m
≤
T
}
≥
ε
∀
m
≥
m
1
.
Define a C
2-function V
:
R
+
3
→
R
+ by
(11)
V
(
S
,
E
,
I
)
=
(
S
-
a
-
a
log
⁡
S
a
)
+
(
E
-
1
-
log
⁡
E
)
+
(
I
-
1
-
log
⁡
I
)
,
where a is a positive constant to be defined later. The nonnegativity of this function can be derived from u
-
1
-
log
⁡
u
≥
0, for all u
>
0. Let m
≥
m
0 and T
>
0 be arbitrary. Using Itô’s formula, we obtain
(12)
d
V
(
S
,
E
,
I
)
=
(
1
-
a
S
)
d
S
+
a
(
d
S
)
2
2
S
2
+
(
1
-
1
E
)
d
E
+
(
d
E
)
2
2
E
2
+
(
1
-
1
I
)
d
I
+
(
d
I
)
2
2
I
2
=
(
1
-
a
S
)
[
(
Λ
-
β
S
I
-
μ
S
+
δ
I
+
γ
E
)
d
t
+
σ
1
S
d
B
1
(
t
)
]
+
1
2
a
σ
1
2
d
t
+
(
1
-
1
E
)
×
[
(
β
S
I
-
(
μ
+
ϵ
+
γ
)
E
)
d
t
+
σ
2
E
d
B
2
(
t
)
]
+
1
2
σ
2
2
d
t
+
(
1
-
1
I
)
[
(
ϵ
E
-
I
×
(
μ
+
δ
+
d
)
)
d
t
+
σ
3
I
d
B
3
(
t
)
]
+
1
2
σ
3
2
d
t
=
[
(
1
-
a
S
)
(
Λ
-
β
S
I
-
μ
S
+
δ
I
+
γ
E
)
+
(
1
-
1
E
)
(
β
S
I
-
(
μ
+
ϵ
+
γ
)
E
)
+
(
1
-
1
I
)
(
ϵ
E
-
(
μ
+
δ
+
d
)
I
)
]
d
t
+
[
σ
1
(
S
-
a
)
d
B
1
(
t
)
+
σ
2
(
E
-
1
)
d
B
2
(
t
)
+
σ
3
(
I
-
1
)
d
B
3
(
t
)
]
=
L
V
(
S
,
E
,
I
)
d
t
+
[
σ
1
(
S
-
a
)
d
B
1
(
t
)
+
σ
2
(
E
-
1
)
d
B
2
(
t
)
+
σ
3
(
I
-
1
)
d
B
3
(
t
)
]
,
where L
V
:
R
+
3
→
R
+ is defined by
(13)
L
V
(
S
,
E
,
I
)
=
(
1
-
a
S
)
(
Λ
-
β
S
I
-
μ
S
+
δ
I
+
γ
E
)
+
(
1
-
1
E
)
(
β
S
I
-
(
μ
+
ϵ
+
γ
)
E
)
+
(
1
-
1
I
)
(
ϵ
E
-
(
μ
+
δ
+
d
)
I
)
=
Λ
+
a
μ
+
2
μ
+
ϵ
+
γ
+
δ
+
d
+
(
a
β
-
μ
-
d
)
I
+
1
2
a
σ
1
2
+
1
2
σ
2
2
+
1
2
σ
3
2
-
μ
E
-
a
Λ
S
-
a
δ
I
S
-
a
γ
E
S
-
β
S
I
E
-
ϵ
E
I
.
Choosing a
=
(
μ
+
d
)
/
β such that (
a
β
-
μ
-
d
)
I
=
0, then
(14)
L
V
(
S
,
E
,
I
)
=
Λ
+
a
μ
+
2
μ
+
ϵ
+
γ
+
δ
+
d
+
1
2
a
σ
1
2
+
1
2
σ
2
2
+
1
2
σ
3
2
-
μ
E
-
a
Λ
S
-
a
δ
I
S
-
a
γ
E
S
-
β
S
I
E
-
ϵ
E
I
≤
Λ
+
a
μ
+
2
μ
+
ϵ
+
γ
+
δ
+
d
+
1
2
a
σ
1
2
+
1
2
σ
2
2
+
1
2
σ
3
2
=
:
K
.
The remainder of the proof follows that Theorem 2.1 in Ji et al. [25].


## 3. Asymptotic Behavior around the Disease-Free Equilibrium of the Deterministic Model

As mentioned in Section1, for the deterministic SEIS system (2), there always exists a disease-free equilibrium F
0
=
(
Λ
/
μ
,
0,0
). And if R
0
=
β
Λ
ϵ
/
μ
(
μ
+
ϵ
+
γ
)
(
μ
+
δ
+
d
)
≤
1, then F
0 is globally asymptotically stable, which means that the disease will vanish after some period of time. However, there is no disease-free equilibrium in stochastic model (3); it is natural to ask how we can consider the disease will go to extinction. In this section we mainly use the method of estimating the oscillation around F
0 to reflect how the solution of model (3) spirals closely around F
0. We have the following results.Theorem 2.
IfR
0
=
β
Λ
ϵ
/
μ
(
μ
+
ϵ
+
γ
)
(
μ
+
δ
+
d
)
≤
1 and the following condition is satisfied
(15)
σ
1
2
<
μ
,
1
2
σ
2
2
<
μ
,
1
2
σ
3
2
<
A
,
where A
=
(
2
μ
(
μ
+
δ
+
d
)
(
2
μ
+
ϵ
+
δ
+
d
) + ϵ
(
2
μ
+
ϵ
)
(
μ
+
d
) + ϵ
d
(
2
μ
+
ϵ
+
d
)
)
/
(
2
μ
(
2
μ
+
ϵ
+
δ
+
d
) + ϵ
(
2
μ
+
ϵ
+
d
)
)
>
0, then for any given initial value (
S
(
0
)
,
E
(
0
)
,
I
(
0
)
)
∈
R
+
3, the solution of model (3) has the property
(16)
lim
⁡
sup
⁡
t
→
∞
1
t
E
∫
0
t
‍
[
(
S
(
s
)
-
Λ
μ
)
2
+
E
2
(
s
)
+
I
2
(
s
)
]
d
s
≤
ϵ
σ
1
2
Λ
2
μ
2
K
1
,
where K
1
=
min
⁡
{
ϵ
(
μ
-
σ
1
2
)
,
ϵ
(
μ
-
(
1
/
2
)
σ
2
2
)
,
a
2
(
A
-
(
1
/
2
)
σ
3
2
)
} and a
2
=
(
2
μ
(
2
μ
+
ϵ
+
δ
+
d
)
+
ϵ
(
2
μ
+
ϵ
+
d
)
)
/
ϵ.Proof.
Define a functionV
:
R
+
3
→
R
+ by
(17)
V
(
x
)
=
ϵ
2
(
S
-
Λ
μ
+
E
)
2
+
(
2
μ
+
ϵ
)
(
S
-
Λ
μ
+
E
)
I
+
a
1
E
+
a
2
2
I
2
+
a
3
I
,
where a
1, a
2, and a
3 are positive constants to be chosen later. For simplicity, we divide (17) into two functions: V
(
x
)
=
V
1
(
x
)
+
V
2
(
x
), where
(18)
V
1
(
x
)
=
ϵ
2
(
S
-
Λ
μ
+
E
)
2
+
(
2
μ
+
ϵ
)
(
S
-
Λ
μ
+
E
)
I
,
V
2
(
x
)
=
a
1
E
+
a
2
2
I
2
+
a
3
I
.
From Itô’s formula, we compute
(19)
d
V
1
(
x
)
=
L
V
1
d
t
+
[
ϵ
(
S
-
Λ
μ
+
E
)
+
(
2
μ
+
ϵ
)
I
]
×
(
σ
1
S
d
B
1
(
t
)
+
σ
2
E
d
B
2
(
t
)
)
+
σ
3
(
2
μ
+
ϵ
)
(
S
-
Λ
μ
+
E
)
I
d
B
3
(
t
)
,
d
V
2
(
x
)
=
L
V
2
d
t
+
σ
2
a
1
E
d
B
2
(
t
)
+
σ
3
(
a
2
I
+
a
3
)
I
d
B
3
(
t
)
.
In detail,
(20)
L
V
1
(
x
)
=
[
ϵ
(
S
-
Λ
μ
+
E
)
+
(
2
μ
+
ϵ
)
I
]
[
Λ
-
μ
S
+
δ
I
-
(
μ
+
ϵ
)
E
]
+
(
2
μ
+
ϵ
)
(
S
-
Λ
μ
+
E
)
×
(
ϵ
E
-
(
μ
+
δ
+
d
)
I
)
+
1
2
ϵ
σ
1
2
S
2
+
1
2
ϵ
σ
2
2
E
2
=
-
μ
ϵ
(
S
-
Λ
μ
)
2
+
[
ϵ
δ
-
(
2
μ
+
ϵ
)
(
2
μ
+
δ
+
d
)
]
×
(
S
-
Λ
μ
)
I
-
[
2
μ
(
2
μ
+
ϵ
+
δ
+
d
)
+
ϵ
(
2
μ
+
ϵ
+
d
)
]
E
I
-
μ
ϵ
E
2
+
δ
(
2
μ
+
ϵ
)
I
2
+
1
2
ϵ
σ
1
2
(
S
-
Λ
μ
+
Λ
μ
)
2
+
1
2
ϵ
σ
2
2
E
2
≤
-
ϵ
(
μ
-
σ
1
2
)
(
S
-
Λ
μ
)
2
-
ϵ
(
μ
-
1
2
σ
2
2
)
E
2
+
δ
(
2
μ
+
ϵ
)
I
2
-
[
2
μ
(
2
μ
+
δ
+
d
)
+
ϵ
(
2
μ
+
d
)
]
(
S
-
Λ
μ
)
I
-
[
2
μ
(
2
μ
+
ϵ
+
δ
+
d
)
+
ϵ
(
2
μ
+
ϵ
+
d
)
]
E
I
+
ϵ
Λ
2
μ
2
σ
1
2
,
L
V
2
(
x
)
=
a
1
(
β
S
I
-
(
μ
+
ϵ
+
γ
)
I
)
+
(
a
2
I
+
a
3
)
(
ϵ
E
-
(
μ
+
δ
+
d
)
I
)
+
1
2
a
2
σ
3
2
I
2
=
β
a
1
(
S
-
Λ
μ
)
I
-
a
2
(
μ
+
δ
+
d
)
I
2
+
ϵ
a
2
E
I
+
[
ϵ
a
3
-
(
μ
+
ϵ
+
γ
)
a
1
]
E
+
[
β
Λ
μ
a
1
-
(
μ
+
δ
+
d
)
a
3
]
I
+
1
2
a
2
σ
3
2
I
2
.
Taking (20) together, we get
(21)
L
V
(
x
)
=
L
V
1
(
x
)
+
L
V
2
(
x
)
≤
-
ϵ
(
μ
-
σ
1
2
)
(
S
-
Λ
μ
)
2
-
ϵ
(
μ
-
1
2
σ
2
2
)
E
2
-
[
a
2
(
μ
+
δ
+
d
)
-
δ
(
2
μ
+
ϵ
)
-
1
2
a
2
σ
3
2
]
I
2
-
[
2
μ
(
2
μ
+
δ
+
d
+
ϵ
)
+
ϵ
d
-
β
a
1
]
(
S
-
Λ
μ
)
I
-
[
2
μ
(
2
μ
+
ϵ
+
δ
+
d
)
+
ϵ
(
2
μ
+
ϵ
+
d
)
-
ϵ
a
2
]
E
I
+
[
ϵ
a
3
-
(
μ
+
ϵ
+
γ
)
a
1
]
E
+
[
β
Λ
μ
a
1
-
(
μ
+
δ
+
d
)
a
3
]
I
+
ϵ
Λ
2
μ
2
σ
1
2
.
Choose a
1
=
(
2
μ
(
2
μ
+
δ
+
d
+
ϵ
)
+
ϵ
d
)
/
β, a
2
=
(
2
μ
(
2
μ
+
ϵ
+
δ
+
d
)
+
ϵ
(
2
μ
+
ϵ
+
d
)
)
/
ϵ, and a
3
=
(
μ
+
ϵ
+
γ
)
a
1
/
ϵ; then 2
μ
(
2
μ
+
δ
+
d
+
ϵ
)
+
ϵ
d
-
β
a
1
=
0, 2
μ
(
2
μ
+
ϵ
+
δ
+
d
)
+
ϵ
(
2
μ
+
ϵ
+
d
)
-
ϵ
a
2
=
0, and ϵ
a
3
-
(
μ
+
ϵ
+
γ
)
a
1
=
0. Besides, noting that R
0
≤
1, thus (21) takes the following:
(22)
L
V
(
x
)
≤
-
ϵ
(
μ
-
σ
1
2
)
(
S
-
Λ
μ
)
2
-
ϵ
(
μ
-
1
2
σ
2
2
)
E
2
-
[
a
2
(
μ
+
δ
+
d
)
-
δ
(
2
μ
+
ϵ
)
-
1
2
a
2
σ
3
2
]
I
2
+
(
μ
+
ϵ
+
γ
)
(
μ
+
δ
+
d
)
a
1
ϵ
(
R
0
-
1
)
I
+
ϵ
Λ
2
μ
2
σ
1
2
≤
-
ϵ
(
μ
-
σ
1
2
)
(
S
-
Λ
μ
)
2
-
ϵ
(
μ
-
1
2
σ
2
2
)
E
2
-
a
2
(
A
-
1
2
σ
3
2
)
I
2
+
ϵ
Λ
2
μ
2
σ
1
2
,
where A
=
(
μ
+
δ
+
d
)
-
(
δ
(
2
μ
+
ϵ
)
/
a
2
) = (
2
μ
(
μ
+
δ
+
d
)
(
2
μ
+
ϵ
+
δ
+
d
) + ϵ
(
2
μ
+
ϵ
)
(
μ
+
d
)
+
ϵ
d
(
2
μ
+
ϵ
+
d
)
)
/
(
2
μ
(
2
μ
+
ϵ
+
δ
+
d
) + ϵ
(
2
μ
+
ϵ
+
d
)
)
>
0. Therefore,
(23)
d
V
(
x
)
≤
[
-
ϵ
(
μ
-
σ
1
2
)
(
S
-
Λ
μ
)
2
-
ϵ
(
μ
-
1
2
σ
2
2
)
E
2
-
a
2
(
A
-
1
2
σ
3
2
)
I
2
+
ϵ
Λ
2
μ
2
σ
1
2
]
d
t
+
σ
1
S
×
[
ϵ
(
S
-
Λ
μ
+
E
)
+
(
2
μ
+
ϵ
)
I
]
d
B
1
(
t
)
+
[
ϵ
(
S
-
Λ
μ
+
E
)
+
(
2
μ
+
ϵ
)
I
+
a
1
]
σ
2
E
×
d
B
2
(
t
)
+
σ
3
(
a
2
I
+
a
3
)
I
d
B
3
(
t
)
.
Since a
2
=
(
2
μ
(
2
μ
+
ϵ
+
δ
+
d
)
+
ϵ
(
2
μ
+
ϵ
+
d
)
)
/
ϵ
>
(
2
μ
+
ϵ
)
2
/
ϵ, the positive definiteness of the quadratic polynomial of the Lyapunov function V is satisfied; thus, the Lyapunov function V is nonnegative. Integrating both sides of (23) from 0 to t and then taking expectation yield
(24)
0
≤
E
V
(
x
(
t
)
)
≤
V
(
x
(
0
)
)
+
E
∫
0
t
[
-
ϵ
(
μ
-
σ
1
2
)
(
S
(
s
)
-
Λ
μ
)
2
-
ϵ
(
μ
-
1
2
σ
2
2
)
E
2
(
s
)
-
a
2
(
A
-
1
2
σ
3
2
)
I
2
(
s
)
+
ϵ
Λ
2
μ
2
σ
1
2
]
d
s
,
which implies
(25)
E
∫
0
t
[
ϵ
(
μ
-
σ
1
2
)
(
S
(
s
)
-
Λ
μ
)
2
+
ϵ
(
μ
-
1
2
σ
2
2
)
E
2
(
s
)
+
a
2
(
A
-
1
2
σ
3
2
)
I
2
(
s
)
]
d
s
≤
V
(
x
(
0
)
)
+
ϵ
σ
1
2
Λ
2
μ
2
t
.
Therefore,
(26)
lim
⁡
sup
⁡
t
→
∞
1
t
E
∫
0
t
[
ϵ
(
μ
-
σ
1
2
)
(
S
(
s
)
-
Λ
μ
)
2
+
ϵ
(
μ
-
1
2
σ
2
2
)
E
2
(
s
)
+
a
2
(
A
-
1
2
σ
3
2
)
I
2
(
s
)
(
S
(
s
)
-
Λ
μ
)
2
]
d
s
≤
ϵ
σ
1
2
Λ
2
μ
2
.
If we let K
1
=
min
⁡
{
ϵ
(
μ
-
σ
1
2
)
,
ϵ
(
μ
-
(
1
/
2
)
σ
2
2
)
,
a
2
(
A
-
(
1
/
2
)
σ
3
2
)
}, then
(27)
limsup
⁡
t
→
∞
1
t
E
∫
0
t
‍
[
(
S
(
s
)
-
Λ
μ
)
2
+
E
2
(
s
)
+
I
2
(
s
)
]
d
s
≤
ϵ
σ
1
2
Λ
2
μ
2
K
1
.
The proof of Theorem 2 is thus completed.Remark 3.
Theorem2 reveals that the solution will oscillate around the disease-free equilibrium, and the intensity is relevant to the values of σ
1, σ
2, and σ
3. The weaker the values are, the smaller the fluctuation is. In other words, if the stochastic perturbations become small, the solution of system (3) will be close to the disease-free equilibrium of system (2). In addition, when σ
1
=
0, then (22) can be reduced to the following form:
(28)
L
V
(
x
)
≤
-
ϵ
μ
(
S
-
Λ
μ
)
2
-
ϵ
(
μ
-
1
2
σ
2
2
)
E
2
-
a
2
(
A
-
1
2
σ
3
2
)
I
2
,
which means that L
V
(
x
) is negative definite provided that (
1
/
2
)
σ
2
2
<
μ and (
1
/
2
)
σ
3
2
<
A. Therefore, F
0 is stochastically asymptotically stable in the large [24].


## 4. Asymptotic Behavior around the Endemic Equilibrium of the Deterministic Model

In this section, we assumeR
0
>
1. Then there is an endemic equilibrium F
* for system (2) but not the endemic equilibrium F
* for system (3), as there is no endemic equilibrium in system (3). Similarly, we also expect to find out whether or not the solution goes around F
*. The following result gives a positive answer under certain parametric restrictions.Theorem 4.
Assume thatR
0
=
β
Λ
ϵ
/
μ
(
μ
+
ϵ
+
γ
)
(
μ
+
δ
+
d
)
>
1; then one has the following two cases to discuss.
(i) If(29)
0
<
B
2
<
4
μ
γ
(
2
μ
+
ϵ
+
γ
)
(
μ
+
δ
+
d
)
×
[
(
μ
+
ϵ
)
(
μ
+
δ
)
+
μ
δ
]
,
and parameters σ
1, σ
2, and σ
3 satisfy conditions
(30)
σ
1
2
<
2
ρ
1
2
μ
+
ϵ
+
γ
,
σ
2
2
<
2
ρ
2
γ
,
σ
3
2
<
4
a
ρ
3
b
,
where
(31)
B
=
(
2
μ
+
ϵ
)
[
(
μ
+
ϵ
)
(
μ
+
d
)
+
μ
δ
]
-
ϵ
δ
γ
,
ρ
1
=
2
μ
ϵ
(
2
μ
+
ϵ
+
γ
)
C
-
|
B
|
2
C
,
ρ
2
=
4
a
γ
(
μ
+
ϵ
)
D
-
(
2
a
γ
δ
+
b
ϵ
)
4
a
D
,
ρ
3
=
2
ϵ
b
(
μ
+
δ
+
d
)
-
2
a
C
|
B
|
-
ϵ
D
(
2
a
γ
δ
+
b
ϵ
)
4
a
ϵ
,
a
=
μ
ϵ
4
(
2
μ
+
ϵ
+
γ
)
,
b
=
2
μ
γ
ϵ
2
(
2
μ
+
ϵ
+
γ
)
[
2
(
μ
+
ϵ
)
(
μ
+
d
)
+
2
μ
δ
+
ϵ
δ
]
,
C
=
|
B
|
3
μ
ϵ
(
2
μ
+
ϵ
+
γ
)
+
8
a
b
ϵ
γ
(
μ
+
δ
+
d
)
(
μ
+
ϵ
)
-
ϵ
(
2
a
γ
δ
+
b
ϵ
)
2
24
γ
a
2
|
B
|
(
μ
+
ϵ
)
,
D
=
2
a
γ
δ
+
b
ϵ
6
a
γ
(
μ
+
ϵ
)
+
2
μ
ϵ
2
b
(
μ
+
δ
+
d
)
(
2
μ
+
ϵ
+
γ
)
-
a
B
2
3
μ
ϵ
2
(
2
a
γ
δ
+
b
ϵ
)
(
2
μ
+
ϵ
+
γ
)
,
then for any given initial value (
S
(
0
)
,
E
(
0
)
,
I
(
0
)
)
∈
R
+
3, the solution of model (3) has the property
(32)
lim
⁡
sup
⁡
t
→
∞
1
t
E
∫
0
t
[
(
S
(
r
)
-
2
ρ
1
2
ρ
1
-
(
2
μ
+
ϵ
+
γ
)
σ
1
2
S
*
)
2
+
(
E
(
r
)
-
2
ρ
2
2
ρ
2
-
γ
σ
2
2
E
*
)
2
+
(
I
(
r
)
-
4
a
ρ
3
4
a
ρ
3
-
b
σ
3
2
I
*
)
2
]
d
r
≤
K
σ
M
,
where
(33)
K
σ
=
ρ
1
(
2
μ
+
ϵ
+
γ
)
σ
1
2
2
ρ
1
-
(
2
μ
+
ϵ
+
γ
)
σ
1
2
S
*
2
+
ρ
2
γ
σ
2
2
2
ρ
2
-
γ
σ
2
2
E
*
2
+
ρ
3
b
σ
3
2
4
a
ρ
3
-
b
σ
3
2
I
*
2
,
M
=
min
⁡
{
ρ
1
-
2
μ
+
ϵ
+
γ
2
σ
1
2
,
ρ
2
-
γ
2
σ
2
2
,
ρ
3
-
b
4
a
σ
3
2
}
.
(ii) IfB
=
0 and parameters σ
1, σ
2, and σ
3 satisfy conditions
(34)
σ
1
2
<
2
μ
,
σ
2
2
<
2
μ
,
σ
3
2
<
2
(
μ
+
d
)
,
then for any given initial value (
S
(
0
)
,
E
(
0
)
,
I
(
0
)
)
∈
R
+
3, the solution of model (3) has the property
(35)
lim
⁡
sup
⁡
t
→
∞
1
t
E
∫
0
t
[
(
S
(
r
)
-
2
μ
2
μ
-
σ
1
2
S
*
)
2
+
(
E
(
r
)
-
2
μ
2
μ
-
σ
2
2
E
*
)
2
+
(
I
(
r
)
-
2
(
μ
+
d
)
2
(
μ
+
d
)
-
σ
3
2
I
*
)
2
]
d
r
≤
K
σ
′
M
′
,
where
(36)
K
σ
′
=
μ
ϵ
(
2
μ
+
ϵ
+
γ
)
σ
1
2
γ
(
2
μ
-
σ
1
2
)
S
*
2
+
μ
ϵ
σ
2
2
2
μ
-
σ
2
2
E
*
2
+
δ
(
μ
+
d
)
σ
3
2
2
(
μ
+
d
)
-
σ
3
2
I
*
2
,
M
′
=
min
⁡
{
ϵ
(
2
μ
+
ϵ
+
γ
)
(
2
μ
-
σ
1
2
)
2
γ
,
ϵ
(
2
μ
-
σ
2
2
)
2
,
δ
[
2
(
μ
+
d
)
-
σ
3
2
]
2
}
.Proof.
Define aC
2-function V
:
R
+
3
→
R
+ by
(37)
V
(
x
)
=
b
1
2
(
S
-
S
*
+
E
-
E
*
)
2
+
b
2
2
(
S
-
S
*
)
2
+
b
3
2
(
I
-
I
*
)
2
,
where b
1, b
2, and b
3 are positive constants to be chosen later. For simplicity, we divide (37) into two functions: V
(
x
)
=
V
1
(
x
)
+
V
2
(
x
), where
(38)
V
1
(
x
)
=
b
1
2
(
S
-
S
*
+
E
-
E
*
)
2
+
b
2
2
(
S
-
S
*
)
2
,
V
2
(
x
)
=
b
3
2
(
I
-
I
*
)
2
.
From Itô’s formula, we have
(39)
d
V
1
(
x
)
=
L
V
1
d
t
+
[
b
1
(
S
-
S
*
+
E
-
E
*
)
+
b
2
(
S
-
S
*
)
]
σ
1
S
d
B
1
(
t
)
+
b
1
σ
1
E
(
S
-
S
*
+
E
-
E
*
)
d
B
2
(
t
)
,
d
V
2
(
x
)
=
L
V
2
d
t
+
σ
3
I
[
b
3
(
I
-
I
*
)
+
(
E
-
E
*
)
]
d
B
3
(
t
)
.
In detail,
(40)
L
V
1
(
x
)
=
b
1
(
S
-
S
*
+
E
-
E
*
)
[
Λ
-
μ
S
+
δ
I
-
(
μ
+
ϵ
)
E
]
+
(
Λ
-
β
S
I
-
μ
S
+
δ
I
+
γ
E
)
×
b
2
(
S
-
S
*
)
+
1
2
b
1
σ
1
2
S
2
+
1
2
b
1
σ
2
2
E
2
+
1
2
b
2
σ
1
2
S
2
=
b
1
(
S
-
S
*
+
E
-
E
*
)
×
[
-
μ
(
S
-
S
*
)
+
δ
(
I
-
I
*
)
-
(
μ
+
ϵ
)
(
E
-
E
*
)
]
+
b
2
(
S
-
S
*
)
×
[
-
β
S
*
(
I
-
I
*
)
-
β
I
(
S
-
S
*
)
-
μ
(
S
-
S
*
)
+
δ
(
I
-
I
*
)
+
γ
(
E
-
E
*
)
]
+
b
1
+
b
2
2
σ
1
2
S
2
+
1
2
b
1
σ
2
2
E
2
=
-
μ
b
1
(
S
-
S
*
)
2
+
b
1
δ
(
S
-
S
*
)
(
I
-
I
*
)
-
b
1
(
μ
+
ϵ
)
(
E
-
E
*
)
2
+
b
1
δ
(
E
-
E
*
)
(
I
-
I
*
)
-
b
1
(
2
μ
+
ϵ
)
(
S
-
S
*
)
(
E
-
E
*
)
-
b
2
β
S
*
(
S
-
S
*
)
(
I
-
I
*
)
-
β
b
2
I
(
S
-
S
*
)
2
-
μ
b
2
(
S
-
S
*
)
2
+
b
2
δ
(
S
-
S
*
)
(
I
-
I
*
)
+
b
2
γ
(
S
-
S
*
)
(
E
-
E
*
)
+
b
1
+
b
2
2
σ
1
2
S
2
+
1
2
b
1
σ
2
2
E
2
≤
-
μ
(
b
1
+
b
2
)
(
S
-
S
*
)
2
-
b
1
(
μ
+
ϵ
)
(
E
-
E
*
)
2
+
[
-
b
1
(
2
μ
+
ϵ
)
+
b
2
γ
]
(
S
-
S
*
)
(
E
-
E
*
)
+
(
b
1
δ
+
b
2
δ
-
b
2
β
S
*
)
(
S
-
S
*
)
(
I
-
I
*
)
+
b
1
δ
(
E
-
E
*
)
(
I
-
I
*
)
+
b
1
+
b
2
2
σ
1
2
S
2
+
1
2
b
1
σ
2
2
E
2
,
L
V
2
(
x
)
=
b
3
(
I
-
I
*
)
(
ϵ
E
-
(
μ
+
δ
+
d
)
I
)
+
1
2
b
3
σ
3
2
I
2
=
b
3
(
I
-
I
*
)
[
ϵ
(
E
-
E
*
)
-
(
μ
+
δ
+
d
)
(
I
-
I
*
)
]
+
1
2
b
3
σ
3
2
I
2
=
b
3
ϵ
(
E
-
E
*
)
(
I
-
I
*
)
-
b
3
(
μ
+
δ
+
d
)
(
I
-
I
*
)
2
+
1
2
b
3
σ
3
2
I
2
.
Taking (40) together, noting that S
*
=
(
μ
+
ϵ
)
(
μ
+
δ
+
d
)
/
β
ϵ, we get
(41)
L
V
(
x
)
=
L
V
1
(
x
)
+
L
V
2
(
x
)
≤
-
μ
(
b
1
+
b
2
)
(
S
-
S
*
)
2
-
b
1
(
μ
+
ϵ
)
(
E
-
E
*
)
2
-
b
3
(
μ
+
δ
+
d
)
(
I
-
I
*
)
2
-
[
b
1
(
2
μ
+
ϵ
)
-
b
2
γ
]
(
S
-
S
*
)
(
E
-
E
*
)
+
[
b
1
δ
-
(
μ
+
ϵ
)
(
μ
+
d
)
+
μ
δ
ϵ
b
2
]
(
S
-
S
*
)
(
I
-
I
*
)
+
(
b
1
δ
+
b
3
ϵ
)
(
E
-
E
*
)
(
I
-
I
*
)
+
b
1
+
b
2
2
σ
1
2
S
2
+
1
2
b
1
σ
2
2
E
2
+
1
2
b
3
σ
3
2
I
2
.
Choose b
2
=
(
(
2
μ
+
ϵ
)
/
γ
)
b
1, from the inequality ±
x
y
≤
(
1
/
2
c
)
x
2
+
2
c
y
2, where c
>
0; then (41) takes the following:
(42)
L
V
(
x
)
≤
-
μ
(
2
μ
+
ϵ
+
γ
)
γ
b
1
(
S
-
S
*
)
2
-
(
μ
+
ϵ
)
b
1
(
E
-
E
*
)
2
-
b
3
(
μ
+
δ
+
d
)
(
I
-
I
*
)
2
-
B
b
1
ϵ
γ
(
S
-
S
*
)
(
I
-
I
*
)
+
(
b
1
δ
+
b
3
ϵ
)
(
E
-
E
*
)
(
I
-
I
*
)
+
2
μ
+
ϵ
+
γ
2
γ
b
1
σ
1
2
S
2
+
1
2
b
1
σ
2
2
E
2
+
1
2
b
3
σ
3
2
I
2
≤
-
μ
(
2
μ
+
ϵ
+
γ
)
γ
b
1
(
S
-
S
*
)
2
-
(
μ
+
ϵ
)
b
1
(
E
-
E
*
)
2
-
b
3
(
μ
+
δ
+
d
)
(
I
-
I
*
)
2
+
|
B
|
b
1
2
ϵ
γ
C
(
S
-
S
*
)
2
+
C
|
B
|
b
1
2
ϵ
γ
(
I
-
I
*
)
2
+
b
1
δ
+
b
3
ϵ
2
D
(
E
-
E
*
)
2
+
D
(
b
1
δ
+
b
3
ϵ
)
2
×
(
I
-
I
*
)
2
+
2
μ
+
ϵ
+
γ
2
γ
b
1
σ
1
2
S
2
+
1
2
b
1
σ
2
2
E
2
+
1
2
b
3
σ
3
2
I
2
=
-
[
μ
(
2
μ
+
ϵ
+
γ
)
-
|
B
|
2
ϵ
C
]
b
1
γ
(
S
-
S
*
)
2
-
[
(
μ
+
ϵ
)
b
1
-
b
1
δ
+
b
3
ϵ
2
D
]
(
E
-
E
*
)
2
-
[
b
3
(
μ
+
δ
+
d
)
-
C
|
B
|
b
1
2
ϵ
γ
-
D
(
b
1
δ
+
b
3
ϵ
)
2
]
×
(
I
-
I
*
)
2
+
2
μ
+
ϵ
+
γ
2
γ
b
1
σ
1
2
S
2
+
1
2
b
1
σ
2
2
E
2
+
1
2
b
3
σ
3
2
I
2
,
where B
=
(
2
μ
+
ϵ
)
[
(
μ
+
ϵ
)
(
μ
+
d
)
+
μ
δ
]
-
ϵ
δ
γ and C, D are positive constants to be chosen later. We have to consider the following two cases: (i) B
≠
0 and (ii) B
=
0.
First we consider case (i); that is,B
2
>
0; we choose b
1
=
γ, b
3
=
b
/
2
a, C
=
(
|
B
|
/
3
μ
ϵ
(
2
μ
+
ϵ
+
γ
)
)
+
(
(
8
a
b
ϵ
γ
(
μ
+
δ
+
d
)
(
μ
+
ϵ
)
-
ϵ
(
2
a
γ
δ
+
b
ϵ
)
2
)
/
24
γ
a
2
|
B
|
(
μ
+
ϵ
)
), and D
=
(
(
2
a
γ
δ
+
b
ϵ
)
/
6
a
γ
(
μ
+
ϵ
)
)
+
(
(
2
μ
ϵ
2
b
(
μ
+
δ
+
d
)
(
2
μ
+
ϵ
+
γ
)
-
a
B
2
)
/
3
μ
ϵ
2
(
2
a
γ
δ
+
b
ϵ
)
(
2
μ
+
ϵ
+
γ
)
), where a
=
μ
ϵ
4
(
2
μ
+
ϵ
+
γ
) and b
=
2
μ
γ
ϵ
2
(
2
μ
+
ϵ
+
γ
)
[
2
(
μ
+
ϵ
)
(
μ
+
d
)
+
2
μ
δ
+
ϵ
δ
] (The choice of parameters
C
and
D
can be seen in Remark5.). Then the quantities in the brackets of (42) are positive. For convenience sake, we introduce ρ
1
=
(
2
μ
ϵ
(
2
μ
+
ϵ
+
γ
)
C
-
|
B
|
)
/
2
C
>
0, ρ
2
=
(
4
a
γ
(
μ
+
ϵ
)
D
-
(
2
a
γ
δ
+
b
ϵ
)
)
/
4
a
D
>
0, and ρ
3
=
(
2
ϵ
b
(
μ
+
δ
+
d
)
-
2
a
C
|
B
|
-
ϵ
D
(
2
a
γ
δ
+
b
ϵ
)
)
/
4
a
ϵ
>
0. Then, from (42), we have
(43)
L
V
(
x
)
≤
-
ρ
1
(
S
-
S
*
)
2
-
ρ
2
(
E
-
E
*
)
2
-
ρ
3
(
I
-
I
*
)
2
+
2
μ
+
ϵ
+
γ
2
σ
1
2
S
2
+
γ
2
σ
2
2
E
2
+
b
4
a
σ
3
2
I
2
=
-
(
ρ
1
-
2
μ
+
ϵ
+
γ
2
σ
1
2
)
S
2
+
2
ρ
1
S
*
S
-
ρ
1
S
*
2
-
(
ρ
2
-
γ
2
σ
2
2
)
E
2
+
2
ρ
2
E
*
E
-
ρ
2
E
*
2
+
(
ρ
3
-
b
4
a
σ
3
2
)
I
2
+
2
ρ
3
I
*
I
-
ρ
3
I
*
2
=
-
(
ρ
1
-
2
μ
+
ϵ
+
γ
2
σ
1
2
)
×
(
S
-
2
ρ
1
2
ρ
1
-
(
2
μ
+
ϵ
+
γ
)
σ
1
2
S
*
)
2
-
(
ρ
2
-
γ
2
σ
2
2
)
(
E
-
2
ρ
2
2
ρ
2
-
γ
σ
2
2
E
*
)
2
-
(
ρ
3
-
b
4
a
σ
3
2
)
(
I
-
4
a
ρ
3
4
a
ρ
3
-
b
σ
3
2
I
*
)
2
+
ρ
1
(
2
μ
+
ϵ
+
γ
)
σ
1
2
2
ρ
1
-
(
2
μ
+
ϵ
+
γ
)
σ
1
2
S
*
2
+
ρ
2
γ
σ
2
2
2
ρ
2
-
γ
σ
2
2
E
*
2
+
ρ
3
b
σ
3
2
4
a
ρ
3
-
b
σ
3
2
I
*
2
.
Thus
(44)
d
V
(
x
)
=
L
V
d
t
+
[
γ
(
S
-
S
*
+
E
-
E
*
)
+
(
2
μ
+
ϵ
)
(
S
-
S
*
)
]
×
σ
1
S
d
B
1
(
t
)
+
γ
σ
1
E
(
S
-
S
*
+
E
-
E
*
)
d
B
2
(
t
)
+
σ
3
I
[
b
1
(
I
-
I
*
)
+
(
E
-
E
*
)
]
d
B
3
(
t
)
.
Integrating both sides of (44) from 0 to t and then taking expectation yield
(45)
0
≤
E
V
(
x
(
t
)
)
≤
V
(
x
(
0
)
)
-
E
∫
0
t
[
(
I
(
r
)
-
4
a
ρ
3
4
a
ρ
3
-
b
σ
3
2
I
*
)
2
(
ρ
1
-
2
μ
+
ϵ
+
γ
2
σ
1
2
)
×
(
S
(
r
)
-
2
ρ
1
2
ρ
1
-
(
2
μ
+
ϵ
+
γ
)
σ
1
2
S
*
)
2
+
(
ρ
2
-
γ
2
σ
2
2
)
×
(
E
(
r
)
-
2
ρ
2
2
ρ
2
-
γ
σ
2
2
E
*
)
2
+
(
ρ
3
-
b
4
a
σ
3
2
)
(

I
(
r
)
-
4
a
ρ
3
4
a
ρ
3
-
b
σ
3
2
I
*

)
2
]

d
r
+
K
σ
t
,
where
(46)
K
σ
=
ρ
1
(
2
μ
+
ϵ
+
γ
)
σ
1
2
2
ρ
1
-
(
2
μ
+
ϵ
+
γ
)
σ
1
2
S
*
2
+
ρ
2
γ
σ
2
2
2
ρ
2
-
γ
σ
2
2
E
*
2
+
ρ
3
b
σ
3
2
4
a
ρ
3
-
b
σ
3
2
I
*
2
,
which implies
(47)
E
∫
0
t
[
(
I
(
r
)
-
4
a
ρ
3
4
a
ρ
3
-
b
σ
3
2
I
*
)
2
(
ρ
1
-
2
μ
+
ϵ
+
γ
2
σ
1
2
)
×
(
S
(
r
)
-
2
ρ
1
2
ρ
1
-
(
2
μ
+
ϵ
+
γ
)
σ
1
2
S
*
)
2
+
(
ρ
2
-
γ
2
σ
2
2
)
×
(
E
(
r
)
-
2
ρ
2
2
ρ
2
-
γ
σ
2
2
E
*
)
2
+
(
ρ
3
-
b
4
a
σ
3
2
)
(
I
(
r
)
-
4
a
ρ
3
4
a
ρ
3
-
b
σ
3
2
I
*
)
2
]
d
r
≤
V
(
x
(
0
)
)
+
K
σ
t
.
Dividing both sides by t and letting t
→
∞, we get
(48)
lim
⁡
sup
⁡
t
→
∞
1
t
E
∫
0
t
[
(
I
(
r
)
-
4
a
ρ
3
4
a
ρ
3
-
b
σ
3
2
I
*
)
2
(
ρ
1
-
2
μ
+
ϵ
+
γ
2
σ
1
2
)
×
(
S
(
r
)
-
2
ρ
1
2
ρ
1
-
(
2
μ
+
ϵ
+
γ
)
σ
1
2
S
*
)
2
+
(
ρ
2
-
γ
2
σ
2
2
)
×
(
E
(
r
)
-
2
ρ
2
2
ρ
2
-
γ
σ
2
2
E
*
)
2
+
(

ρ
3
-
b
4
a
σ
3
2

)
(

I
(
r
)
-
4
a
ρ
3
4
a
ρ
3
-
b
σ
3
2
I
*

)
2

]
d
r
≤
K
σ
.
Introduce M
=
min
⁡
{
ρ
1
-
(
(
2
μ
+
ϵ
+
γ
)
/
2
)
σ
1
2
,
ρ
2
-
(
γ
/
2
)
σ
2
2
,
ρ
3
-
(
b
/
4
a
)
σ
3
2
}; then
(49)
lim
⁡
sup
⁡
t
→
∞
1
t
E
∫
0
t
[
(
S
(
r
)
-
2
ρ
1
2
ρ
1
-
(
2
μ
+
ϵ
+
γ
)
σ
1
2
S
*
)
2
+
(
E
(
r
)
-
2
ρ
2
2
ρ
2
-
γ
σ
2
2
E
*
)
2
+
(
I
(
r
)
-
4
a
ρ
3
4
a
ρ
3
-
b
σ
3
2
I
*
)
2
]
d
r
≤
K
σ
M
.
For case (ii), we can choose b
1
=
ϵ, b
3
=
δ, and D
=
δ
/
ϵ; then from (42), we can obtain
(50)
L
V
(
x
)
≤
-
μ
ϵ
(
2
μ
+
ϵ
+
γ
)
γ
(
S
-
S
*
)
2
-
ϵ
(
μ
+
ϵ
)
(
E
-
E
*
)
2
-
δ
(
μ
+
δ
+
d
)
(
I
-
I
*
)
2
+
2
ϵ
δ
(
E
-
E
*
)
(
I
-
I
*
)
+
ϵ
(
2
μ
+
ϵ
+
γ
)
2
γ
σ
1
2
S
2
+
ϵ
2
σ
2
2
E
2
+
δ
2
σ
3
2
I
2
≤
-
μ
ϵ
(
2
μ
+
ϵ
+
γ
)
γ
(
S
-
S
*
)
2
-
μ
ϵ
(
E
-
E
*
)
2
-
δ
(
μ
+
d
)
(
I
-
I
*
)
2
+
ϵ
(
2
μ
+
ϵ
+
γ
)
2
γ
σ
1
2
S
2
+
ϵ
2
σ
2
2
E
2
+
δ
2
σ
3
2
I
2
=
-
ϵ
(
2
μ
+
ϵ
+
γ
)
(
2
μ
-
σ
1
2
)
2
γ
(
S
(
s
)
-
2
μ
2
μ
-
σ
1
2
S
*
)
2
-
ϵ
(
2
μ
-
σ
2
2
)
2
(
E
(
s
)
-
2
μ
2
μ
-
σ
2
2
E
*
)
2
-
δ
[
2
(
μ
+
d
)
-
σ
3
2
]
2
(
I
(
s
)
-
2
(
μ
+
d
)
2
(
μ
+
d
)
-
σ
3
2
I
*
)
2
+
K
σ
′
,
where
(51)
K
σ
′
=
μ
ϵ
(
2
μ
+
ϵ
+
γ
)
σ
1
2
γ
(
2
μ
-
σ
1
2
)
S
*
2
+
μ
ϵ
σ
2
2
2
μ
-
σ
2
2
E
*
2
+
δ
(
μ
+
d
)
σ
3
2
2
(
μ
+
d
)
-
σ
3
2
I
*
2
.
In a similar way, we can get
(52)
lim
⁡
sup
⁡
t
→
∞
1
t
E
∫
0
t
[
ϵ
(
2
μ
+
ϵ
+
γ
)
(
2
μ
-
σ
1
2
)
2
γ
×
(
S
(
r
)
-
2
μ
2
μ
-
σ
1
2
S
*
)
2
+
ϵ
(
2
μ
-
σ
2
2
)
2
(
E
(
r
)
×
-
2
μ
2
μ
-
σ
2
2
E
*
)
2
+
δ
[
2
(
μ
+
d
)
-
σ
3
2
]
2
×
(
I
(
r
)
-
2
(
μ
+
d
)
2
(
μ
+
d
)
-
σ
3
2
I
*
)
2
]
d
r
≤
K
σ
′
.
Set M
′
=
min
⁡
{
ϵ
(
2
μ
+
ϵ
+
γ
)
(
2
μ
-
σ
1
2
)
/
2
γ
,
ϵ
(
2
μ
-
σ
2
2
)
/
2
,
δ
[
2
(
μ
+
d
)
-
σ
3
2
]
/
2
}; then
(53)
lim
⁡
sup
⁡
t
→
∞
1
t
E
∫
0
t
[
(
S
(
r
)
-
2
μ
2
μ
-
σ
1
2
S
*
)
2
+
(
E
(
r
)
-
2
μ
2
μ
-
σ
2
2
E
*
)
2
+
(
I
(
r
)
-
2
(
μ
+
d
)
2
(
μ
+
d
)
-
σ
3
2
I
*
)
2
]
d
r
≤
K
σ
′
M
′
.
This completes the theorem.Remark 5 (the choice of parameters
C
and
D).
In order to ensure that the quantities in the brackets of (42) are positive, C and D must satisfy the following three inequalities:
(54)
C
>
|
B
|
2
μ
ϵ
(
2
μ
+
ϵ
+
γ
)
,
D
>
γ
δ
+
b
3
ϵ
2
γ
(
μ
+
ϵ
)
,
C
|
B
|
2
ϵ
+
D
(
γ
δ
+
b
3
ϵ
)
2
<
b
3
(
μ
+
δ
+
d
)
.
In order to guarantee the existence of positive constants C and D, from (54), the following inequality must be satisfied:
(55)
B
2
4
μ
ϵ
2
(
2
μ
+
ϵ
+
γ
)
+
(
γ
δ
+
b
3
ϵ
)
2
4
γ
(
μ
+
ϵ
)
<
b
3
(
μ
+
δ
+
d
)
.
If we introduce y
=
b
3, then the above inequality is equivalent to the following form:
(56)
h
(
y
)
:
=
a
y
2
-
b
y
+
c
<
0
,
where a
=
μ
ϵ
4
(
2
μ
+
ϵ
+
γ
)
>
0, b
=
2
μ
γ
ϵ
2
(
2
μ
+
ϵ
+
γ
)
[
2
(
μ
+
ϵ
)
(
μ
+
d
)
+
2
μ
δ
+
ϵ
δ
]
>
0, and c
=
γ
(
μ
+
ϵ
)
B
2
+
μ
ϵ
2
γ
2
δ
2
(
2
μ
+
ϵ
+
γ
)
>
0. Furthermore, h
(
0
)
=
c
>
0. To ensure that inequality (56) holds, there exists a positive solution; the constants a
,
b, and c must satisfy the following condition:
(57)
b
2
-
4
a
c
=
4
μ
γ
ϵ
4
(
2
μ
+
ϵ
+
γ
)
(
μ
+
ϵ
)
×
{
[
(
μ
+
ϵ
)
B
2
(
μ
+
d
)
+
μ
δ
]
4
μ
γ
(
2
μ
+
ϵ
+
γ
)
(
μ
+
δ
+
d
)
×
[
(
μ
+
ϵ
)
(
μ
+
d
)
+
μ
δ
]
-
B
2
}
>
0
.
By condition (29), inequality (57) holds; then we choose y
=
b
/
2
a; that is b
3
=
b
/
2
a; therefore, (54) can be written in the following form;
(58)
C
>
|
B
|
2
μ
ϵ
(
2
μ
+
ϵ
+
γ
)
,
D
>
2
a
γ
δ
+
b
ϵ
4
a
γ
(
μ
+
ϵ
)
,
2
a
|
B
|
C
+
(
2
a
γ
δ
+
b
ϵ
)
ϵ
D
<
2
ϵ
b
(
μ
+
δ
+
d
)
.
The corresponding equalities of the above inequalities give three straight lines. The feasible region of C and D is a triangle area surrounded by these three straight lines. The third straight line can be written as follows:
(59)
l
:
2
a
|
B
|
C
+
(
2
a
γ
δ
+
b
ϵ
)
ϵ
D
=
2
ϵ
b
(
μ
+
δ
+
d
)
.
For convenience sake, we denote C
*
=
|
B
|
/
2
μ
ϵ
(
2
μ
+
ϵ
+
γ
) and D
*
=
(
2
a
γ
δ
+
b
ϵ
)
/
4
a
γ
(
μ
+
ϵ
). Substituting C
=
C
* into (59) gets D
=
(
2
ϵ
b
(
μ
+
δ
+
d
)
-
2
a
|
B
|
C
*
)
/
ϵ
(
2
a
γ
δ
+
b
ϵ
)
:
=
D
*, and by substituting D
=
D
* into it, we obtain C
=
(
2
ϵ
b
(
μ
+
δ
+
d
)
-
(
2
a
γ
δ
+
b
ϵ
)
ϵ
D
*
)
/
2
a
|
B
|
:
=
C
*. We choose C = C
* + (
1
/
3
)
(
C
*
-
C
*
) = (
|
B
|
/
3
μ
ϵ
(
2
μ
+
ϵ
+
γ
)
) + (
(
8
a
b
ϵ
γ
(
μ
+
δ
+
d
)
(
μ
+
ϵ
)
-
ϵ
(
2
a
γ
δ
+
b
ϵ
)
2
)
/
24
γ
a
2
|
B
|
(
μ
+
ϵ
)
) and D = D
* + (
1
/
3
)
(
D
*
-
D
*
) = (
(
2
a
γ
δ
+
b
ϵ
)
/
6
a
γ
(
μ
+
ϵ
)
) + (
(
2
μ
ϵ
2
b
(
μ
+
δ
+
d
)
(
2
μ
+
ϵ
+
γ
)
-
a
B
2
)
/
3
μ
ϵ
2
(
2
a
γ
δ
+
b
ϵ
)
(
2
μ
+
ϵ
+
γ
)
), which are easily satisfied the inequalities (54); that is, the quantities in the brackets of (42) are positive.Remark 6.
Theorem4 shows that if 0
<
B
2
<
4
μ
γ
(
2
μ
+
ϵ
+
γ
)
(
μ
+
δ
+
d
)
[
(
μ
+
ϵ
)
(
μ
+
δ
)
+
μ
δ
], the solution of model (3) fluctuates around a certain level which is relevant to P
* and σ
i, where P
*
=
(
(
2
ρ
1
/
(
2
ρ
1
-
(
2
μ
+
ϵ
+
γ
)
σ
1
2
)
)
S
*
,
(
2
ρ
2
/
(
2
ρ
2
-
γ
σ
2
2
)
)
E
*
,
(
4
a
ρ
3
/
(
4
a
ρ
3
-
b
σ
3
2
)
)
I
*
), i
=
1,2
,
3. With the value of σ
i decreasing, P
* will be close to F
* and the difference between X and P
* will also decrease, where X denotes the solution of system (3). If B
=
0, we have similar results.Remark 7.
As the limitation of method, whether or not the solutions of model (3) have similar results is unknown to us provided that B
2
≥
4
μ
γ
(
2
μ
+
ϵ
+
γ
)
(
μ
+
δ
+
d
)
[
(
μ
+
ϵ
)
(
μ
+
δ
)
+
μ
δ
]. This is our further research work.


## 5. Conclusion and Numerical Simulations

As mentioned in the Introduction, on one hand, deterministic system (2) always has a disease-free equilibrium F
0 and it is globally asymptotically stable if R
0
=
β
Λ
ϵ
/
μ
(
μ
+
ϵ
+
γ
)
(
μ
+
δ
+
d
)
≤
1. For the stochastic system (3), the expectations of S
(
t
)
,
I
(
t
)
,
R
(
t
) are bounded in time average when condition (15) holds. And the smaller the value of σ
1 is, the less the boundedness is. In addition, if σ
1 is reducing to zero, then F
0 is stochastically asymptotically stable in the large provided that (
1
/
2
)
σ
2
2
<
μ and (
1
/
2
)
σ
3
2
<
A. On the other hand, if R
0
>
1, then the disease-free equilibrium of deterministic system (2) becomes unstable and there exists a unique endemic equilibrium F
*, which is globally asymptotically stable in the interior of the feasible region. For the stochastic system (3), the solution of system (3) spirals around the endemic equilibrium of the corresponding deterministic system under some conditions. Furthermore, the fluctuation becomes weaker with intensities decreasing.In order to conform the results above, we numerically simulate the solution of system (3). Using Milsteins higher order method [26], we obtain the following discretization equation:
(60)
S
k
+
1
=
S
k
+
(
Λ
-
β
S
k
I
k
-
μ
S
k
+
δ
I
k
+
γ
E
k
)
Δ
t
+
σ
1
S
k
Δ
t
ξ
1
,
k
+
σ
1
2
2
S
k
Δ
t
(
ξ
1
,
k
2
-
1
)
,
E
k
+
1
=
E
k
+
(
β
S
k
I
k
-
(
μ
+
ϵ
+
γ
)
E
k
)
Δ
t
+
σ
2
E
k
Δ
t
ξ
2
,
k
+
σ
2
2
2
E
k
Δ
t
(
ξ
2
,
k
2
-
1
)
,
I
k
+
1
=
I
k
+
(
μ
E
k
-
(
μ
+
δ
+
d
)
I
k
)
Δ
t
+
σ
3
I
k
Δ
t
ξ
3
,
k
+
σ
3
2
2
I
k
Δ
t
(
ξ
3
,
k
2
-
1
)
,
where time increment Δ
t
>
0 and ξ
1
,
k, ξ
2
,
k, and ξ
3
,
k (k
=
1,2
,
…
,
n) are N
(
0,1
)-distributed mutually independent random variables which can be generated numerically by pseudorandom number generators.In the following figures, the blue lines and the red lines represent solutions of the deterministic system (2) and the stochastic system (3), respectively, and initial values are always S
(
0
)
=
0.7, E
(
0
)
=
0.2, and I
(
0
)
=
0.1. In Figures 1 and 2, we fix Λ
=
0.6, β
=
0.4, μ
=
0.2, δ
=
0.1, γ
=
0.2, ϵ
=
0.15, and d
=
0.1 with different intensities of white noise which satisfy condition (15) in Theorem 2. We can easily compute that R
0
=
0.8182
<
1, F
0
=
(
Λ
/
μ
,
0,0
)
=
(
3,0
,
0
). By MATLAB software, we simulate the solutions of model (3) with different values of σ
i, i
=
1,2
,
3, and the solution of model (2). In Figure 1 Group (a), we choose σ
1
=
0.04, σ
2
=
0.03, and σ
3
=
0.03 and in Figure 1 Group (b), we fix σ
1
=
0.01, σ
2
=
0.02, and σ
3
=
0.02. We can see from Figure 1 that the solution of model (3) will oscillate around the disease-free equilibrium in time, and comparisons of Group (a) and Group (b) suggest that the fluctuations reduce as the noise level decreases.S
(
0
)
=
0.7, E
(
0
)
=
0.2, I
(
0
)
=
0.1, Λ
=
0.6, β
=
0.4, μ
=
0.2, δ
=
0.1, γ
=
0.2, ϵ
=
0.15, d
=
0.1, and Δ
t
=
0.002. Group (a) σ
1
=
0.04, σ
2
=
0.03, and σ
3
=
0.03. Group (b) σ
1
=
0.01, σ
2
=
0.02, and σ
3
=
0.02.
(a)
(b)S
(
0
)
=
0.7, E
(
0
)
=
0.2, I
(
0
)
=
0.1, Λ
=
0.6, β
=
0.4, μ
=
0.2, δ
=
0.1, γ
=
0.2, ϵ
=
0.15, d
=
0.1, and Δ
t
=
0.002. Group (a) σ
1
=
0.04, σ
2
=
0.1, and σ
3
=
0.2. Group (b) σ
1
=
0, σ
2
=
0.1, and σ
3
=
0.2.
(a)
(b)Furthermore, in Figure2 Group (a), the intensities are much larger than those in Figure 1 Group (a), but σ
1 keeps the same. We can conclude from the comparison of Figure 2 Group (a) and Figure 1 Group (a) that the intensities of σ
2 and σ
3 have little effect on the fluctuation. For σ
1
=
0, Figure 2 Group (b) shows that the fluctuation for S
(
t
) almost cannot be seen. Thus we can believe that the solution is stochastically asymptotically stable in the large.In the following, we consider the long behavior of model (3) in the case of R
0
>
1. Let β
=
0.6 and fix the other parameter values as those in Figure 1; we can obtain that R
0
=
1.2273
>
1, F
*
=
(
2.4444,0.3556,0.1333
), and 0.0043
=
B
2
<
4
μ
γ
(
2
μ
+
ϵ
+
γ
)
(
μ
+
δ
+
d
)
[
(
μ
+
ϵ
)
(
μ
+
δ
)
+
μ
δ
]
=
0.006, implying that condition (29) in Theorem 4 holds. In Figure 3 Group (a), we fix σ
1
=
0.04, σ
2
=
0.03, and σ
3
=
0.02. And it is easy to verify that condition (30) holds. Figure 3 Group (a) shows that the solution of system (3) fluctuates around the solution of system (2), which supports the conclusion of Theorem 4. Furthermore, with intensities of white noise decreasing, the fluctuation becomes weaker, which is shown in Figure 3 Group (b).S
(
0
)
=
0.7, E
(
0
)
=
0.2, I
(
0
)
=
0.1, Λ
=
0.6, β
=
0.6, μ
=
0.2, δ
=
0.1, γ
=
0.2, ϵ
=
0.15, d
=
0.1, and Δ
t
=
0.002. Group (a) σ
1
=
0.04, σ
2
=
0.03, and σ
3
=
0.02. Group (b) σ
1
=
0.01, σ
2
=
0.02, and σ
3
=
0.01.
(a)
(b)On the other hand, takeΛ
=
0.4, β
=
0.8, μ
=
0.1, δ
=
0.3, γ
=
0.6, ϵ
=
0.2, and d
=
0.1. By calculating, we have R
0
=
1.4222
>
1, F
*
=
(
2.8125,0.6597,0.2639
), and B
=
0, which satisfy condition (34) in Theorem 4. Choosing different intensities of white noise, we have similar results as above. The specific time sequence diagram to compare the fluctuation is shown in Figure 4 Groups (a) and (b).S
(
0
)
=
0.7, E
(
0
)
=
0.2, I
(
0
)
=
0.1, Λ
=
0.4, β
=
0.8, μ
=
0.1, δ
=
0.3, γ
=
0.6, ϵ
=
0.2, d
=
0.1, and Δ
t
=
0.002. Group (a) σ
1
=
0.04, σ
2
=
0.03, and σ
3
=
0.02. Group (b) σ
1
=
0.01, σ
2
=
0.02, and σ
3
=
0.01.
(a)
(b)Finally, we chooseΛ
=
0.8, β
=
0.6, μ
=
0.2, δ
=
0.2, γ
=
0.3, ϵ
=
0.3, and d
=
0.4. By calculating, we have R
0
=
1.125
>
1, F
*
=
(
3.5556,0.2092,0.0784
), and 0.0484
=
B
2
>
4
μ
γ
(
2
μ
+
ϵ
+
γ
)
(
μ
+
δ
+
d
)
[
(
μ
+
ϵ
)
(
μ
+
δ
)
+
μ
δ
]
=
0.0461; that is, the conditions in Theorem 4 do not hold. From the numerical simulations (see Figure 5 Groups (a) and (b)), it is clear that the solutions of system (3) still fluctuate around the solution of system (2) and the fluctuation becomes weaker with intensities reducing. Therefore, we guess: the solutions of system (3) may fluctutate around the solution of system (2) by removing the condition B
2
≤
4
μ
γ
(
2
μ
+
ϵ
+
γ
)
(
μ
+
δ
+
d
)
[
(
μ
+
ϵ
)
(
μ
+
δ
)
+
μ
δ
]. This will be an open problem to be studied in the future work.S
(
0
)
=
0.7, E
(
0
)
=
0.2, I
(
0
)
=
0.1, Λ
=
0.8, β
=
0.6, μ
=
0.2, δ
=
0.2, γ
=
0.3, ϵ
=
0.3, and d
=
0.4. Group (a) σ
1
=
0.04, σ
2
=
0.03, and σ
3
=
0.02. Group (b) σ
1
=
0.01, σ
2
=
0.02, and σ
3
=
0.01.
(a)
(b)



---

*Source: 894842-2014-10-14.xml*