# Dynamics of a Stochastic Multigroup SEIR Epidemic Model

Authors: Xiaojing Zhong; Feiqi Deng
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

To be more precise about the real world activity, a stochastic multigroup SEIR epidemic model is formulated. we define the basic reproduction numberR
0
S and show that it is a sharp threshold for the dynamics of SDE model. If R
0
S
<
1, the disease-free equilibrium is asymptotically stable; and if R
0
S
>
1, the disease persists and there exists a globally asymptotically stable stationary distribution. Numerical simulation examples are carried out to substantiate the analytical results.


---

## Body

## 1. Introduction

For the past decades, many epidemic models have been proposed for modeling the spread process of infectious diseases, and in the meantime considerable attention has been paid to study the dynamical properties of these various models. Most models descend from the classical SIR model of Kermack and McKendrick [1]; it is the earliest triumphs in mathematical epidemiology. After that, many researchers worked on epidemic models and established different type of epidemic models [2–10]. In particular, multigroup models have been proposed to describe the transmission dynamics of infectious diseases in heterogeneous host populations, such as meals, mumps, gonorrhea, HIV/AIDS, West-Nile virus, and vector borne diseases such as Malaria. One of the earliest works on multigroup models is the seminal paper by Lajmanovich and Yorke [11] on a class of SIS multigroup models for the transmission dynamics of Gonorrhea; they established a complete analysis of the global dynamics. The global stability of the unique equilibrium is proved by using a complete analysis of the global Lyapunov function. Recently, a group-theoretic approach to the method of global Lyapunov function was proposed by Li and Shuai [12]; they studied the following SEIR model:
(1)
S
i
′
=
A
i
-
d
i
S
S
i
-
∑
j
=
1
n
‍
β
i
j
S
i
I
j
,
E
i
′
=
∑
j
=
1
n
‍
β
i
j
S
i
I
j
-
(
d
i
E
+
ϵ
i
)
E
i
,
i
=
1,2
,
…
,
n
,
I
i
′
=
ϵ
i
E
i
-
(
d
i
I
+
γ
i
)
I
i
.The model describes the spread of an infectious disease in a heterogeneous population, which is partitioned inton homogeneous group. Each group i is further compartmentalized into S
i, E
i, and I
i; here S
i, E
i, and I
i denote the subpopulation that are susceptible to the disease, infected but noninfectious, and infectious, respectively. All parameters in the above model are nonnegative constants and summarized in the following list:β
i
j: transmission coefficient between compartments S
i and I
j;
d
i
S
;
d
i
E
;
d
i
I: nature death rates of S, E, and I compartments in the ith group, respectively;
A
i: influx of individuals into the ith group;
ɛ
i: rate of becoming infectious after latent period in the ith group who are immunized;
γ
i: recovery rate of infectious individuals in the ith group.There are two equilibria of this model as disease-free equilibriumP
0
=
(
S
1
0
,
0,0
,
…
,
S
n
0
,
0,0
), where S
i
0
=
A
i
/
d
i
S, and endemic equilibrium P
*
=
(
S
1
*
,
E
1
*
,
I
1
*
,
…
,
S
n
*
,
E
n
*
,
I
n
*
). A threshold R
0 is defined which decides the epidemic will prevalent or not; here,
(2)
R
0
=
ρ
(
M
0
)
denote the spectral radius of the matrix
(3)
M
0
=
M
(
S
1
0
,
S
2
0
,
…
,
S
n
0
)
=
(
β
i
j
ϵ
i
S
i
0
(
d
i
E
+
ϵ
i
)
(
d
i
I
+
γ
i
)
)
1
≤
i
,
j
≤
n
.For more details of it, ifR
0
≤
1, then P
0 is the unique equilibrium and it is globally asymptotically stable in Γ, where Γ is the limit set of system (1):
(4)
Γ
=
{
A
k
min
⁡
⁡
{
d
k
S
,
d
k
E
,
d
k
I
}
(
S
1
,
E
1
,
I
1
,
…
,
S
n
,
E
n
,
I
n
)
∣
S
k
mm
≤
S
k
0
,
S
k
+
E
k
+
I
k
≤
A
k
min
⁡
⁡
{
d
k
S
,
d
k
E
,
d
k
I
}
,
1
≤
k
≤
n
}
.IfR
0
>
1, then P
0 is unstable and it is uniformly persistent. Furthermore, there exists an endemic equilibrium P
* and it is globally asymptotically stable in Γ. In the whole proof, it used a very important group theorem [12].Given a nonnegative matrixA
=
(
β
i
j
), the directed graph G
(
A
) associated with A
=
(
β
i
j
) has vertices 1,2
,
…
,
n with a directed arc (i, j) from i to j if and only if β
i
j
≠
0. It is strongly connected if any two distinct vertices are joined by an oriented path. The matrix A is irreducible if and only if G
(
A
) is strongly connected. A tree is a connected graph with no cycles. A subtree T of a graph G is said to be spanning if T contains all the vertices of G. A directed tree is a tree in which each edge has been replaced by an arc directed one way to the other. A directed tree is said to be rooted at a vertex, called the root, if every arc is oriented in the direction towards the root. An oriented cycle in a directed graph is a simple closed oriented path. A unicyclic graph is a directed graph consisting of a collection of disjoint rooted directed trees whose roots are on an oriented cycle. For a given nonnegative matrix A
=
(
β
i
j
), let
(5)
L
=
[
∑
l
≠
1
‍
β
1
l
-
β
21
⋯
-
β
n
1
-
β
12
∑
l
≠
2
‍
β
2
l
⋯
β
n
2
⋮
⋮
⋯
⋮
-
β
1
n
-
β
2
n
⋯
∑
l
≠
n
‍
β
n
l
]
,
be the Laplacian matrix of the directed graph G
(
A
) and C
i
j denote the cofactor of the (
i
,
j
) entry of L. In light of these results, complete determination of the global dynamic of these models is essential for their application and further development.Whereas the statement above, the large-scale biological system's parameters are assumed as constants, but in the real situation, parameters involved with the model always fluctuate around some average value due to continuous fluctuation in the environment. In order to study the dynamics of interacting population under realistic situation, we need to analyse the associated stochastic model. Stochastic epidemic models have been studied by many authors [13–29], Tornatore et al. [23], Yu et al. [24], Ji et al. [25], Liu et al. [26], and Ji et al. [27] using Lyanupov methods to find out sufficient conditions of the stability of the steady-state based on the deterministic threshold R
0. Gray et al. [28] established a stochastic SIS model and found out a sufficient and necessary condition of the disease-free equilibrium. Hasminskii [29] work on the stochastic persistence of epidemic model and give many stochastic persistence definitions about epidemic model.In the present paper, we introduce white noise into system (1) by perturbing model parameters d
i
S, d
i
E, d
i
I to arrive at the following system of stochastic differential equations:
(6)
d
S
i
=
(
A
i
-
d
i
s
S
i
-
∑
j
=
1
n
‍
β
i
j
S
i
I
j
)
d
t
-
σ
S
i
d
B
(
t
)
,
d
E
i
=
[
∑
j
=
1
n
‍
β
i
j
S
i
I
j
-
(
d
i
E
+
ϵ
i
)
E
i
]
d
t
-
σ
E
i
d
B
(
t
)
,
m
m
m
m
m
m
m
m
m
m
m
m
d
m
m
i
=
1,2
,
…
,
n
,
d
I
i
=
[
ϵ
i
E
i
-
(
d
i
I
+
γ
i
)
I
i
]
d
t
-
σ
I
i
d
B
(
t
)
,
where B
(
t
) is standard brownian motions. Our main objective is to derive a sharp threshold for the extinction and persistence of the disease. We proved that the dynamics of our model is determined by a noise modified basic reproduction number:
(7)
R
0
S
=
ρ
(
(
β
i
j
ϵ
i
S
i
0
(
d
i
E
+
ϵ
i
+
σ
2
/
2
)
(
d
i
I
+
γ
i
+
σ
2
/
2
)
)
1
≤
i
,
j
≤
n
)
.
More specifically, if R
0
S
<
1, the disease-free equilibrium P
0 is asymptotically stable and the disease dies out. If R
0
S
>
1, then P
0 is unstable, system (6) is stochastically persistent, and there exists a stationary distribution. Our definition of R
0
S includes as a special case of a basic reproduction number for a stochastic SIS model in [22]. From (7), we see that R
0
S
<
R
0 if σ are nonzero. This implies that the presence of noise lowers the threshold for the extinction of disease and hence results in a larger parameter region for disease to die out. This agrees with an earlier result on stochastic SIS models in [22] and findings on stochastic logistic equations that the presence of noise increases the parameter region in which the species becomes extinct. Unlike the standard approach of using Lyapunov functions in the literature of SDE epidemic models, our stability analysis of P
0 applied the method of linearization. And we use the recurrence condition to prove the existence of stationary distribution.In this paper, we establish the global existence of positive solutions in Section2. Stability analysis of the disease-free equilibrium is carried out in Section 3. In Section 4, we prove the existence of a globally stable stationary distribution when R
0
S
>
1. Numerical simulations are provided at the end of Sections 3 and 4 to illustrate our analytical results.


## 2. Existence and Uniqueness of the Positive Global Solution

In this section, we prove the positive global existence of our stochastic system’s (6) solution. As a stochastic differential equation, the functions involved with stochastic system are generally required to satisfy the lipschitz condition and linear growth condition. Obviously, the function of system (6) does not satisfy the linear growth condition, so the solution may explode at a finite time, only if we prove that the explosion time is infinite. We use the lyapunov analysis method to confirm our assumption that the solution of our system is global existence and positive.Theorem 1.
IfB
=
(
β
i
j
)
n
×
n is irreducible, then, for any initial value (
S
1
(
0
)
,
E
1
(
0
)
,
I
1
(
0
)
,
…
,
S
n
(
0
)
,
E
n
(
0
)
,
I
n
(
0
)
)
∈
R
+
3
n of system (6), there exists a unique solution (
S
1
(
t
)
,
E
1
(
t
)
,
I
1
(
t
)
,
…
,
S
n
(
t
)
,
E
n
(
t
)
,
I
n
(
t
)
)
∈
R
+
3
n, and it satisfies
(8)
P
(
(
S
1
(
t
)
,
E
1
(
t
)
,
I
1
(
t
)
,
…
,
S
n
(
t
)
,
E
n
(
t
)
,
I
n
(
t
)
)
∣
(
S
1

(
t
)
,
E
1

(
t
)
,
I
1

(
t
)
,
…
,
S
n
(
t
)
,
E
n
(
t
)
,
I
n
(
t
)
)
∈
R
+
3
n
)
=
1
which means (
S
1
(
t
)
,
E
1
(
t
)
,
I
1
(
t
)
,
…
,
S
n
(
t
)
,
E
n
(
t
)
,
I
n
(
t
)
)
∈
R
+
3
n.Proof.
Since the coefficients of the equation are locally Lipschitz continuous, there is a unique local solution ont
∈
[
0
,
τ
e
), where τ
e is the explosion time [30]. We assume the solution (
S
1
(
t
)
,
E
1
(
t
)
I
1
(
t
)
,
…
,
S
n
(
t
)
,
E
n
(
t
)
,
I
n
(
t
)
)
=
Y
(
t
); now, we need to prove Y
(
t
) is global. Let m
0 be sufficiently large so that S
i
(
0
), E
i
(
0
), I
i
(
0
) all lie within the interval [
1
/
m
0
,
m
0
]. For each integer m
≥
m
0, define the stopping time τ
m
=
inf
⁡
{
t
∈
[
0
,
τ
e
)
:
min
⁡
{
S
i
(
t
)
,
E
i
(
t
)
,
I
i
(
t
)
}
≤
1
/
m
or
max
⁡
{
S
i
(
t
)
,
E
i
(
t
)
,
I
i
(
t
)
}
≥
m
}. To complete the proof, we need to show that lim
⁡
m
→
∞
τ
m
=
∞. If this statement is false, then there are constants T
>
0, ɛ
∈
(
0,1
), and m
1
≥
m
0, such that P
(
τ
m
≤
T
)
≥
ɛ for all m
≥
m
1.
Letc
i denote the cofactor of the ith diagonal entry of L
B, which is the Laplacian matrix of (
G
,
B
) and g
(
x
)
=
x
-
1
-
ln
⁡
x. We define
(9)
V
(
Y
(
t
)
)
=
∑
i
=
1
n
‍
[
a
c
i
g
(
S
i
a
c
i
)
+
g
(
E
i
)
+
g
(
I
i
)
]
,
where a is a positive constant. Obviously, V
(
Y
(
t
)
) is positive. Using Itô’s formula, we obtain
(10)
d
V
(
Y
(
t
)
)
=
∑
i
=
1
n
‍
[
+
a
c
i
2
S
i
2
(
d
S
i
)
2
+
1
2
E
i
2
(
d
e
i
)
2
+
1
2
I
i
2
(
d
I
i
)
2
(
1
-
a
c
i
S
i
)
d
S
i
+
(
1
-
1
E
i
)
d
E
i
+
(
1
-
1
I
i
)
d
I
i
m
m
i
m
i
m
m
+
a
c
i
2
S
i
2
(
d
S
i
)
2
+
1
2
E
i
2
(
d
e
i
)
2
+
1
2
I
i
2
(
d
I
i
)
2
]
=
L
V
-
∑
i
=
1
n
‍
[
σ
(
1
-
a
c
i
S
i
)
S
i
d
B
+
σ
(
1
-
1
E
i
)
E
i
d
B
mmmimimimm
+
σ
(
1
-
1
I
i
)
I
i
d
B
]
,
where
(11)
L
V
(
Y
(
t
)
)
=
∑
i
=
1
n
‍
[
K
i
-
(
d
i
I
+
γ
i
)
I
i
+
∑
j
=
1
n
a
c
i
β
i
j
I
j
]
-
∑
i
=
1
n
‍
(
d
i
S
S
i
+
d
i
E
E
i
+
a
c
i
A
i
S
i
+
∑
j
=
1
n
β
i
j
S
i
I
j
E
i
+
ϵ
i
E
i
I
i
)
≤
∑
i
=
1
n
‍
[
K
i
-
(
d
i
I
+
γ
i
)
I
i
+
∑
j
=
1
n
a
c
i
β
i
j
I
j
]
,
where K
i
=
A
i
+
a
c
i
d
i
S
+
d
i
E
+
d
i
I
+
ϵ
i
+
γ
i
+
(
a
c
i
/
2
)
σ
2
+
(
1
/
2
)
σ
2
+
(
1
/
2
)
σ
2. As B
=
(
β
i
j
)
n
×
n is irreducible [12], we know that ∑
i
=
1
n
∑
j
=
1
n
β
i
j
a
c
i
I
j
=
∑
i
=
1
n
∑
j
=
1
n
β
i
j
a
c
i
I
j, which implies that L
V
≤
[
K
i
-
(
d
i
I
+
γ
i
-
∑
j
=
1
n
β
i
j
a
c
i
)
I
i
]. We define a
=
min
⁡
{
(
d
i
I
+
γ
i
)
/
(
∑
j
=
1
n
β
i
j
c
i
)
,
i
=
1,2
,
…
,
n
}; then, we obtain
(12)
L
V
≤
∑
i
=
1
n
‍
K
i
.
Therefore,
(13)
E
(
V
(
Y
(
τ
m
∧
T
)
)
)
=
E
(
∫
0
τ
m
∧
T
d
V
(
Y
(
t
)
)
)
V
(
Y
(
0
)
)
+
∫
0
τ
m
∧
T
d
V
(
Y
(
t
)
)
∫
0
τ
m
∧
T
d
V
(
Y
(
t
)
)
)
)
≤
V
(
Y
(
0
)
)
+
E
∫
0
τ
m
∧
T
M
~
d
t
≤
V
(
Y
(
0
)
)
+
∑
i
=
1
n
‍
K
i
T
.
Set Ω
m
=
{
τ
m
∣
τ
m
≤
T
f
o
r
m
≥
m
1
}. Then, P
(
Ω
m
)
≥
ɛ. Note that, for every ω
∈
Ω
m, there is at least one of S
i
(
τ
m
,
ω
)
,
E
i
(
τ
m
,
ω
)
,
I
i
(
τ
m
,
ω
) that equals either m or 1
/
m. Then,
(14)
V
(
Y
(
τ
m
)
)
≥
min
⁡
0
≤
k
≤
n
(
m
-
a
c
i
-
a
c
i
ln
⁡
m
a
c
i
)
∧
min
⁡
0
≤
k
≤
n
(
1
m
-
a
c
i
-
a
c
i
ln
⁡
1
a
c
i
m
)
,
where we define c
0 such that a
c
0
=
1. Then, we obtain
(15)
V
(
Y
(
0
)
)
+
M
~
T
≥
E
[
1
Ω
m
(
ω
)
V
(
Y
(
τ
m
)
)
]
≥
ɛ
[
min
⁡
0
≤
k
≤
n
(
m
-
a
c
i
-
a
c
i
ln
⁡
m
a
c
i
)
mimm
∧
min
⁡
0
≤
k
≤
n
1
m
-
a
c
i
-
a
c
i
ln
⁡
1
a
c
i
m
]
.
Letting m
→
∞, then ∞
>
V
(
Y
(
0
)
)
+
∑
i
=
1
n
‍
K
i
T
=
∞. So, the assumption is wrong, and we obtain τ
∞
=
∞. The proof is complete.


## 3. Extinction of the Epidemic

In the study of population systems, extinction and persistence are two of the most important issues. For the deterministic model (1), extinction is implied by showing that the disease-free equilibrium (
S
1
0
,
0,0
,
…
,
S
n
0
,
0,0
) is asymptotically stable. For our stochastic model (6), there does not exist the disease-free equilibrium because of the first perturbation term σ
S
i. If E
i
(
t
)
-
I
i
(
t
)
=
0, the first equation in (6) changes to
(16)
d
S
¯
i
=
(
A
i
-
d
i
s
S
¯
i
)
d
t
-
σ
S
¯
i
d
B
(
t
)
.
For this kind of equations, Gray et al. [21] have shown that the solution satisfies that
(17)
lim
⁡
T
→
∞
1
T
∫
0
T
S
¯
i
(
t
)
d
t
=
A
i
d
i
S
.We make the change of variablesu
i
(
t
)
=
S
i
(
t
)
-
S
¯
i
,
v
i
(
t
)
=
E
i
(
t
) and w
(
t
)
=
I
i
(
t
) so that the origin will represent (
S
¯
(
t
)
,
0,0
); then, we consider the linearized system:
(18)
d
u
i
=
(
-
d
i
s
u
i
-
∑
j
=
1
n
‍
β
i
j
S
¯
i
w
j
)
d
t
-
σ
u
i
d
B
(
t
)
,
d
v
i
=
[
∑
j
=
1
n
‍
β
i
j
S
¯
i
w
j
-
(
d
i
E
+
ϵ
i
)
v
i
]
d
t
-
σ
v
i
d
B
(
t
)
,
M
M
M
M
M
M
M
I
I
I
M
M
I
I
I
I
M
M
i
=
1,2
,
…
,
n
,
d
w
i
=
[
ϵ
i
v
i
-
(
d
i
I
+
γ
i
)
w
i
]
d
t
-
σ
w
i
d
B
(
t
)
.To be simplified, we rewrite the second and third equations in (18) as
(19)
d
x
(
t
)
=
F
x
(
t
)
d
t
+
G
x
(
t
)
d
B
(
t
)
,
where
(20)
x
(
t
)
=
[
v
1
(
t
)
,
w
1
(
t
)
,
…
,
v
n
(
t
)
,
w
n
(
t
)
]
T
,
F
=
[
-
(
d
1
E
+
ϵ
1
)
⋯
β
1
n
S
¯
1
⋮
⋱
⋮
0
⋯
-
(
d
n
I
+
γ
n
)
]
,
G
=
[
σ
σ
⋱
σ
]
.
F and G are commute, and the explicit solution of the linearized system (18) is
(21)
x
(
t
)
=
x
(
0
)
exp
⁡
[
(
F
-
1
2
G
2
)
t
+
G
B
(
t
)
]
,
where
(22)
F
-
1
2
G
=
[
-
(
d
1
E
+
ϵ
1
+
1
2
σ
2
)
⋯
β
1
n
S
1
0
⋮
⋱
⋮
0
⋯
-
(
d
n
I
+
γ
n
+
1
2
σ
2
)
]
.
Let R
0
S
=
ρ
(
M
0
S
) denote the spectral radius of the matrix
(23)
M
0
S
=
(
β
i
j
ϵ
i
S
i
0
(
d
i
E
+
ϵ
i
+
σ
2
/
2
)
(
d
i
I
+
γ
i
+
σ
2
/
2
)
)
1
≤
i
,
j
≤
n
if R
0
S
<
1, which means all the eigenvalue of F
-
(
1
/
2
)
G
2 have negative real parts. Then, there is a pair of positive constants C and λ such that
(24)
∥
exp
⁡
[
(
F
-
1
2
G
2
)
t
]
∥
≤
C
e
-
λ
t
.
It then follows that
(25)
|
x
(
t
)
|
≤
C
|
x
(
0
)
|
exp
⁡
⁡
[
-
λ
t
+
∥
G
∥
|
B
(
t
)
|
]
.
Using the strong law of large numbers states that lim
⁡
t
→
∞
(
B
(
t
)
/
t
)
=
0 a.s.; we obtain
(26)
limsup
⁡
t
→
∞
1
t
log
⁡
|
x
(
t
)
|
≤
-
λ
a
.
s
.
In other words, the solution of (18) is almost surely exponentially stable. Next, we give estimate for u
i
(
t
), using Ito
˙’s formula; we derive that
(27)
u
i
(
t
)
=
e
-
(
d
i
S
+
σ
2
/
2
)
t
-
σ
B
(
t
)
×
[
u
i
(
0
)
+
∫
0
t
∑
j
=
1
n
‍
β
i
j
S
¯
i
w
j
(
s
)
e
(
d
i
S
+
σ
2
/
2
)
s
+
σ
B
(
s
)
d
s
]
.
According to (19), there exist T
>
0, w
j
(
t
)
≤
e
-
λ
t. Substituting it to (26), we get
(28)
u
i
(
t
)
=
e
-
(
d
i
S
+
σ
2
/
2
)
t
-
σ
B
(
t
)
×
[
u
i
(
0
)
+
∫
0
T
∑
j
=
1
n
‍
β
i
j
S
¯
i
w
j
(
s
)
e
(
d
i
S
+
σ
2
/
2
)
s
+
σ
B
(
s
)
d
s
m
m
m
+
∫
T
t
∑
j
=
1
n
‍
β
i
j
S
¯
i
w
j
(
s
)
e
(
d
i
S
+
σ
2
/
2
)
s
+
σ
B
(
s
)
d
s
]
≤
e
-
d
i
S
t
[
u
i
(
0
)
+
∫
0
T
∑
j
=
1
n
‍
β
i
j
S
¯
i
w
j
(
s
)
e
d
i
S
s
d
s
]
+
∑
j
=
1
n
‍
β
i
j
S
¯
i
d
i
S
-
λ
[
e
-
λ
t
+
t
l
n
l
n
t
-
e
(
d
i
S
-
λ
)
T
-
d
i
S
t
+
t
l
n
l
n
t
]
.Therefore,(29)
limsup
⁡
t
→
∞
1
t
log
⁡
|
u
i
(
t
)
|
=
-
d
i
S
∨
-
λ
<
0
.
In this way, we proved that the solution of (7) is exponentially stable. According to the Oseledec multiplicative ergodic theorem [30], the necessary and sufficient condition for the almost sure asymptotic stability of the trivial solution of the system is the largest lyapunov exponent of the linearized system is negative. Therefore, we have the following results.Theorem 2.
Assume thatB
=
(
β
i
j
)
n
×
n
) is irreducible.
(1)  IfR
0
S
<
1, then the disease-free equilibrium P
0 is almost sure asymptotically stable, which means the disease will die out almost surely.
(2)  IfR
0
S
>
1, then the disease-free equilibrium P
0 is unstable.Remark 3.
It is useful to observe that for either the classical deterministic model or the stochastic model, there is a threshold which reflects the prevalent or extinction of the epidemic, but the thresholds are different between them; the stochastic thresholdR
0
S is smaller then the deterministic one. In other words, the conditions for I
(
t
) to become extinct in the SDE epidemic model are weaker than in the classical deterministic epidemic model. We give the following example that illustrates this result more explicitely.Example 4.
For simplicity, letk
=
2 and we choose the following system parameters:
(30)
A
1
=
100
;
A
2
=
300
;
d
1
S
=
3
;
d
2
S
=
3
;
d
1
E
=
3
;
d
2
E
=
5
;
d
1
I
=
3
;
d
2
I
=
5
;
β
11
=
0.1
;
β
12
=
0.2
;
β
21
=
0.3
;
β
22
=
0.4
;
ϵ
1
=
1
;
ϵ
2
=
1
;
γ
1
=
1
;
γ
2
=
1
;
so the stochastic multigroup SEIR model (6) becomes
(31)
d
S
1
=
(
100
-
2
S
1
-
0.1
S
1
I
1
-
0.2
S
1
I
2
)
d
t
-
σ
S
1
d
B
(
t
)
,
d
E
1
=
[
0.1
S
1
I
1
+
0.2
S
1
I
2
-
4
E
1
]
d
t
-
σ
E
1
d
B
(
t
)
,
d
I
1
=
(
E
1
-
4
I
1
)
d
t
-
σ
I
1
d
B
(
t
)
,
d
S
2
=
(
300
-
3
S
1
-
0.3
S
2
I
1
-
0.4
S
2
I
2
)
d
t
-
σ
S
2
d
B
(
t
)
,
d
E
2
=
[
0.3
S
2
I
1
+
0.4
S
2
I
2
-
6
E
2
]
dt
-
σ
E
2
d
B
(
t
)
,
d
I
2
=
(
E
2
-
6
I
2
)
d
t
-
σ
I
2
d
B
(
t
)
.
Clearly, if σ
=
0, system (27) becomes the related deterministic multigroup SEIR model. We start our numerical simulation with σ
=
2, and the initial values are I
1
(
0
)
=
10
,
I
2
(
0
)
=
20. Note that R
0
S
=
ρ
(
M
0
S
)
=
(
40
/
63
)
<
1. By Theorem 2, I
1
(
t
)
,
I
2
(
t
) will tend to zero exponentially. If we consider the corresponding deterministic model, R
0
=
(
4
/
3
)
>
1, I
1
(
t
), I
2
(
t
) will tend to their endemic equilibrium. The computer simulation in Figure 1 illustrates extinction of the disease.Computer simulation of pathsI
1
(
t
), I
2
(
t
) for the system (18) and its corresponding deterministic model, using the EM method with step size 0.001, with initial values I
1
(
0
)
=
10, I
2
(
0
)
=
20. The full line expresses stochastic model's simulation, and the dotted line expresses the related deterministic model.
(a)
(b)Next, we keep the parameter value and start our computer simulation at the initial valueI
1
(
0
)
=
I
2
(
0
)
=
1; we gain the same results in Figure 2.Computer simulation of pathsI
1
(
t
), I
2
(
t
) for the system (18) and its corresponding deterministic model, using the EM method with step size 0.001, with initial values I
1
(
0
)
=
1
=
I
2
(
0
)
=
1. The full line expresses stochastic model's simulation, and the dotted line expresses the related deterministic model.
(a)
(b)If we decrease the environment intensity toσ
=
0.5 and starting from I
1
(
0
)
=
I
2
(
0
)
=
1, this means that R
0
S
>
1. From Theorem 2, the disease-free equilibrium will be unstable; results of one simulation run in Figure 3 proved our results.Computer simulation of pathsI
1
(
t
), I
2
(
t
) for the system (18) and its corresponding deterministic model, using the EM method with step size 0.001, with initial values I
1
(
0
)
=
1, I
2
(
0
)
=
1. The full line expresses stochastic model's simulation, and the dotted line expresses the related deterministic model.
(a)
(b)


## 4. Stationary Distribution

As we know stochastic persistence means if solution trajectories start from a positive initial condition, then they will remain within the positive interior and bounded at all future times. If we prove the existence of stationary distribution of our stochastic multigroup SEIR model, it means the disease will persist. Before proving the main theorem, we reference to the book by Hasminskii [29]. Let X
(
t
) be a regular time-homogeneous Markov process described by the SDE
(32)
d
X
(
t
)
=
b
(
X
)
d
t
+
∑
r
=
1
k
‍
σ
r
(
X
)
d
B
r
(
t
)
.The diffusion matrix is defined as follows:(33)
A
(
x
)
=
(
a
i
j
(
x
)
)
,
a
i
j
(
x
)
=
∑
r
=
1
k
‍
σ
r
i
(
x
)
σ
r
j
(
x
)
.Lemma 5.
The Markov processX
(
t
) has a unique stationary distribution μ if there exists an open bounded domain U
⊂
R
l, and the conditions are satisfied.
(P1). In the domainU and some neighborhood thereof, the smallest eigenvalue of the diffusion matrix A
(
x
) is bounded away from zero.
(P2). Ifx
∈
R
l
∖
U, the mean time τ at which a path issuing from x reaches the set U is finite, and sup
⁡
x
∈
K
E
x
τ
<
∞ for every compact subset K
⊂
R
l. Let f
(
·
) be a function integrable with respect to the measure μ. Then,
(34)
P
(
lim
⁡
T
→
∞
⁡
1
T
∫
0
T
f
(
X
x
(
t
)
)
d
t
=
∫
R
l
f
(
x
)
μ
(
(
d
)
x
)
)
=
1
,
for all x
∈
R
l.Remark 6.
The proof of the above lemma is given in Hasminskii [29]. The existence of a stationary distribution with density is given in Theorem 4.1 at page 119 and Lemma 9.4 at page 138.
To validate (P1), it sufficient to prove thatF is uniformly elliptical in U, where F
u
=
b
(
x
)
u
x
+
(
1
/
2
)
(
tr
⁡
(
A
(
x
)
u
x
x
)
), which means there is a positive number M such that ∑
i
j
=
1
k
a
i
j
(
x
)
ξ
i
ξ
j
≥
M
|
ξ
|
2 for any x
∈
U. To validate (P2), it is enough to show that there exist some neighborhood U and a nonnegative C
2-function V such that, for any x
∈
R
l
∖
U, L
V is negative definite function (for details, see page 1163 in [31]).Theorem 7.
Assume thatB
=
(
β
i
j
)
n
×
n is irreducible and R
0
S
>
1. There is a stationary distribution for system (7) and it has an ergodic property.Proof.
Since1
<
R
0
S
<
R
0, there is an endemic equilibrium P
*
=
(
S
1
*
,
E
1
*
,
I
1
*
,
…
,
S
n
*
,
E
n
*
,
I
n
*
) for the deterministic system of (7). We obtain the following equation:
(35)
A
i
=
d
i
S
S
i
*
+
∑
j
=
1
n
‍
β
i
j
S
i
*
I
j
*
;
∑
j
=
1
n
‍
β
i
j
S
i
*
I
j
*
=
(
d
i
E
+
ϵ
i
)
E
i
;
ϵ
i
E
i
=
(
d
i
I
+
γ
i
)
I
i
.
Using the same method in the proof of Theorem 1.1 [7], we choose β
¯
i
j
=
β
i
j
S
i
*
I
j
*, 1
≤
i, j
≤
n, B
¯
=
(
β
¯
i
j
), {
v
1
,
…
,
v
n
}, v
i
>
0 such that B
¯
v
=
0. Set
(36)
V
=
∑
i
=
1
n
‍
v
i
[
d
i
E
+
ϵ
i
ϵ
i
(
I
i
-
I
i
*
ln
⁡
I
i
)
(
S
i
-
S
i
*
ln
⁡
S
i
)
+
(
E
i
-
E
i
*
ln
⁡
E
i
)
m
i
m
m
+
d
i
E
+
ϵ
i
ϵ
i
(
I
i
-
I
i
*
ln
⁡
I
i
)
]
.
Applying Itô’s formula, we can calculate
(37)
d
V
=
∑
i
=
1
n
‍
v
i
[
d
i
E
+
ϵ
i
ϵ
i
(
1
-
I
i
*
I
i
)
d
I
i
+
(
d
i
E
+
ϵ
i
)
I
i
*
2
ϵ
i
I
i
2
(
d
I
i
)
2
(
1
-
S
i
*
S
i
)
d
S
i
+
S
i
*
2
S
i
2
(
d
S
i
)
2
mmmmmm
+
(
1
-
E
i
*
E
i
)
d
E
i
+
E
i
*
2
E
i
2
(
d
E
i
)
2
mmmmmm
+
d
i
E
+
ϵ
i
ϵ
i
(
1
-
I
i
*
I
i
)
d
I
i
+
(
d
i
E
+
ϵ
i
)
I
i
*
2
ϵ
i
I
i
2
(
d
I
i
)
2
]
.
Substituting (6) to it and using (31), we obtain
(38)
L
V
=
∑
i
=
1
n
‍
v
i
[
A
i
-
d
i
S
S
i
-
∑
j
=
1
n
‍
β
i
j
S
i
I
j
-
A
i
S
i
*
S
i
-
d
i
S
S
i
*
mmmmim
+
∑
j
=
1
n
‍
β
i
j
S
i
*
I
j
+
S
i
*
σ
2
2
mmmmim
+
∑
j
=
1
n
‍
β
i
j
S
i
I
j
-
(
d
i
E
+
ϵ
i
)
E
i
mmmmim
-
∑
j
=
1
n
‍
β
i
j
E
i
*
S
i
I
j
E
i
+
(
d
i
E
+
ϵ
i
)
E
i
*
+
E
i
*
σ
2
2
mmmmim
+
(
d
i
E
+
ϵ
i
)
E
i
-
(
d
i
E
+
ϵ
i
)
(
d
i
I
+
γ
i
)
ϵ
i
I
i
mmmmim
-
(
d
i
E
+
ϵ
i
)
I
i
*
E
i
I
i
+
(
d
i
E
+
ϵ
i
)
(
d
i
I
+
γ
i
)
ϵ
i
I
i
*
mmmmim
+
(
d
i
E
+
ϵ
i
)
I
i
*
σ
2
2
ϵ
i
A
i
-
d
i
S
S
i
-
∑
j
=
1
n
‍
β
i
j
S
i
I
j
-
A
i
S
i
*
S
i
-
d
i
S
S
i
*
]
=
∑
i
=
1
n
‍
v
i
[
S
i
*
σ
2
+
E
i
*
σ
2
2
+
(
d
i
E
+
ϵ
i
)
I
i
*
2
ϵ
i
σ
2
d
i
S
S
i
*
(
2
-
S
i
*
S
i
-
S
i
S
i
*
)
mmiimim
+
(
∑
j
=
1
n
‍
β
i
j
S
i
*
I
j
-
(
d
i
E
+
ϵ
i
)
(
d
i
I
+
γ
i
)
ϵ
i
I
i
)
mmiimim
+
(
3
∑
j
=
1
n
‍
β
i
j
S
i
*
I
j
*
-
∑
j
=
1
n
‍
β
i
j
I
j
*
(
S
i
*
)
2
S
i
mmiimmmmim
-
∑
j
=
1
n
‍
β
i
j
S
i
I
j
E
i
*
E
i
-
(
d
i
E
+
ϵ
i
)
E
i
I
i
*
I
i
)
mmiimim
+
S
i
*
σ
2
+
E
i
*
σ
2
2
+
(
d
i
E
+
ϵ
i
)
I
i
*
2
ϵ
i
σ
2
]
≤
∑
i
=
1
n
‍
v
i
[
(
∑
j
=
1
n
‍
β
i
j
S
i
*
I
j
-
(
d
i
E
+
ϵ
i
)
(
d
i
I
+
γ
i
)
ϵ
i
I
i
)
mmiimim
+
S
i
*
σ
2
+
E
i
*
σ
2
2
+
(
d
i
E
+
ϵ
i
)
I
i
*
2
ϵ
i
σ
2
mmiimim
+
(
∑
j
=
1
n
‍
β
i
j
S
i
I
j
E
i
*
E
i
-
(
d
i
E
+
ϵ
i
)
E
i
I
i
*
I
i
3
∑
j
=
1
n
‍
β
i
j
S
i
*
I
j
*
-
∑
j
=
1
n
‍
β
i
j
I
j
*
(
S
i
*
)
2
S
i
mmiimmmmim
-
∑
j
=
1
n
‍
β
i
j
S
i
I
j
E
i
*
E
i
-
(
d
i
E
+
ϵ
i
)
E
i
I
i
*
I
i
)
]
.
Since B
¯
v
=
0, it follows that ∑
j
=
1
n
β
¯
j
i
v
j
=
∑
k
=
1
n
β
¯
i
k
v
i; using (31), we obtain
(39)
∑
j
=
1
n
‍
β
j
i
S
j
*
S
j
*
I
i
*
v
j
=
∑
k
=
1
n
‍
β
i
k
S
i
*
I
k
*
v
i
=
∑
j
=
1
n
‍
β
i
j
S
i
*
I
j
*
v
i
=
(
d
i
E
+
ϵ
i
)
(
d
i
I
+
γ
i
)
ϵ
i
I
i
*
v
i
which means
(40)
L
V
≤
∑
i
=
1
n
‍
v
i
[
∑
j
=
1
n
‍
β
i
j
S
i
I
j
E
i
*
E
i
-
(
d
i
E
+
ϵ
i
)
E
i
I
i
*
I
i
S
i
*
σ
2
+
E
i
*
σ
2
2
+
(
d
i
E
+
ϵ
i
)
I
i
*
2
ϵ
i
σ
2
mmmmmm
+
3
∑
j
=
1
n
‍
β
i
j
S
i
*
I
j
*
-
∑
j
=
1
n
‍
β
i
j
I
j
*
(
S
i
*
)
2
S
i
m
m
m
m
m
m
-
∑
j
=
1
n
‍
β
i
j
S
i
I
j
E
i
*
E
i
-
(
d
i
E
+
ϵ
i
)
E
i
I
i
*
I
i
]
=
H
(
S
1
,
E
1
,
I
1
,
…
,
S
n
,
E
n
,
I
n
)
.
Note that
(41)
lim
⁡
S
i
,
E
i
,
I
i
→
0
H
(
S
1
,
E
1
,
I
1
,
…
,
S
n
,
E
n
,
I
n
)
=
-
∞
,
lim
⁡
S
i
,
E
i
,
I
i
→
∞
H
(
S
1
,
E
1
,
I
1
,
…
,
S
n
,
E
n
,
I
n
)
=
-
∞
.
So, there exists a domain U lying entirely in R
+
3
n. For (
S
1
,
E
1
,
I
1
,
…
,
S
n
,
E
n
,
I
n
)
∈
U
∖
R
+
3
n, L
V
<
-
M, where M is a positive constant. It implies that condition (P2) is satisfied. Besides, there is a K
=
min
⁡
{
σ
S
i
2
,
σ
E
i
2
,
σ
I
i
2
,
i
=
1,2
,
…
,
n
}
>
0 such that ∑
i
,
j
=
1
3
n
a
i
j
ξ
i
ξ
j
=
∑
i
=
1
n
σ
2
S
i
2
ξ
3
i
-
2
2
+
∑
i
=
1
n
σ
2
I
i
2
ξ
3
i
2
+
∑
i
=
1
n
σ
2
E
i
2
ξ
3
i
-
1
2
≥
K
∥
ξ
2
∥, which implies that condition (P2) is satisfied. Therefore, according to Lemma 5, our stochastic SEIR model (6) has a stationary distribution and it is ergodic. The proof is complete.Example 8.
To substantiate the analytic findings above, we provide numerical simulation results for the stochastic model (27). We also use the same parameters in Example 4, and let σ
=
0.5. We have shown in Figure 3 that I
1
(
t
), I
2
(
t
) will not tend to 0. Theorem 4.3 tells us that there is a stationary distribution. Figure 4 shows histograms of the approximate stationary distribution of the infective classes.Frequency histograms of pathsI
1
(
t
), I
2
(
t
) for the system (31) based on 10000 stochastic sumulations for each population at time t
=
100, using the EM method with step size 0.001, with initial values I
1
(
0
)
=
1, I
2
(
0
)
=
1.
(a)
(b)



---

*Source: 258915-2014-07-24.xml*