The Stability of SI Epidemic Model in Complex Networks with Stochastic Perturbation

Abstract

We investigate a stochastic SI epidemic model in the complex networks. We show that this model has a unique global positive solution. Then we consider the asymptotic behavior of the model around the disease-free equilibrium and show that the solution will oscillate around the disease-free equilibrium of deterministic system when 𝑅₀≤1. Furthermore, we derive that the disease will be persistent when 𝑅₀>1. Finally, a series of numerical simulations are presented to illustrate our mathematical findings. A new result is given such that, when 𝑅₀≤1, with the increase of noise intensity the solution of stochastic system converging to the disease-free equilibrium is faster than that of the deterministic system.

1. Introduction

Epidemiology is the science to study the distribution of disease and influencing factors, so as to explore the etiology, clarify the popular rule of the disease, and formulate the countermeasures and measures for preventing, controlling, and eliminating the disease. Many mathematical models of diseases spreading help us to understand the propagation of diseases [1, 2]. The transmission of diseases can be influenced by many factors, such as the age and social structure of the population, the contact network among individuals, and the metapopulation characteristics. So it is difficult to establish an accurate epidemic model which is completely consistent with the real world. In recent years, a lot of compartmental epidemic models have been studied by many researchers [3–5], and complex networks also have been used to study the spread of diseases [6–17].

In this paper we consider an SI model with the birth and death in complex networks. As mentioned in the paper [6, 13], the birth and death do not affect the degree of nodes. Suppose 𝑆_𝑘(𝑡),𝐼_𝑘(𝑡) are the number of the healthy and infected nodes with the degree 𝑘 at time 𝑡; the mean-field equations can be written as 𝑑𝑆_𝑘𝑑𝑡=𝑏_𝑘−𝜆𝑘𝑆_𝑘𝜃−𝑑𝑆_𝑘,𝑑𝐼_𝑘𝑑𝑡=𝜆𝑘𝑆_𝑘𝜃−(𝑑+𝜖)𝐼_𝑘, (1) where 𝜃=(1/⟨𝑘⟩)∑^𝑛_𝑘=1𝜆𝑘𝑃(𝑘). For system (1), it can be written as the following form: 𝑑𝑆_𝑘(𝑡)=𝑏_𝑘−1⟨𝑘⟩^𝑛_𝑗=1𝜆𝑘𝑗𝑃𝑗𝑆_𝑘(𝑡)𝐼_𝑗(𝑡)−𝑑𝑆_𝑘(𝑡)𝑑𝑡,𝑑𝐼_𝑘(𝑡)=1⟨𝑘⟩^𝑛_𝑗=1𝜆𝑘𝑗𝑃𝑗𝑆_𝑘(𝑡)𝐼_𝑗(𝑡)−(𝑑+𝜖)𝐼_𝑘(𝑡)𝑑𝑡. (2) We denote 𝛽_𝑘𝑗=(1/⟨𝑘⟩)𝜆𝑘𝑗𝑃(𝑗), so we obtain 𝑑𝑆_𝑘(𝑡)=𝑏_𝑘−^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘(𝑡)𝐼_𝑗(𝑡)−𝑑𝑆_𝑘(𝑡)𝑑𝑡,𝑑𝐼_𝑘(𝑡)=^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘(𝑡)𝐼_𝑗(𝑡)−(𝑑+𝜖)𝐼_𝑘(𝑡)𝑑𝑡. (3) It always has the disease-free equilibrium 𝐸₀=(𝑆⁰₁,0,…,𝑆⁰_𝑛,0), where 𝑆⁰_𝑘=𝑏_𝑘/𝑑, 𝑘=1,2,…,𝑛. If 𝐴=(𝛽_𝑘𝑗)_𝑛×𝑛 is irreducible and 𝑅₀≤1, then 𝐸₀ is globally stable in 𝐷, while if 𝑅₀>1, 𝐸₀ is unstable and there is an endemic equilibrium 𝐸^∗=(𝑆^∗₁,𝐼^∗₁,…,𝑆^∗_𝑛,𝐼^∗_𝑛) belonging to 𝐷 which is globally asymptotically stable in 𝐷; here 𝐷=𝑆₁,𝐼₁,…,𝑆_𝑛,𝐼_𝑛∈ℝ₊^2𝑛∶𝑆_𝑘,𝐼_𝑘≤𝑏_𝑘𝑑,𝑆_𝑘+𝐼_𝑘≤𝑏_𝑘𝑑,,𝑀𝑘=1,2,…,𝑛₀=𝑀𝑆⁰₁,…,𝑆⁰_𝑛=𝛽_𝑘𝑗𝑆⁰_𝑘𝑑+𝜖_𝑛×𝑛,𝑅₀=𝜌𝑀₀, (4) and 𝜌(𝑀₀) denotes the spectral radius of 𝑀₀.

The deterministic models have some limitations in describing the spread of disease. The accident in the process of disease transmission can not be reflected by the deterministic models. This is because of the fact that the deterministic models ignore the effect of the environmental noise. In an ecosystem, the environmental noise is inevitably in the real world; thus stochastic models are more realistic. In the research of stochastic epidemic models, many researchers make a lot of contributions [17–26].

In this paper, we consider the following stochastic system: 𝑑𝑆_𝑘(𝑡)=𝑏_𝑘−^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘(𝑡)𝐼_𝑗(𝑡)−𝑑𝑆_𝑘(𝑡)𝑑𝑡+𝜎_𝑘1𝑆_𝑘(𝑡)𝑑𝐵_𝑘1(𝑡),𝑑𝐼_𝑘(𝑡)=^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘(𝑡)𝐼_𝑗(𝑡)−(𝑑+𝜖)𝐼_𝑘(𝑡)𝑑𝑡+𝜎_𝑘2𝐼_𝑘(𝑡)𝑑𝐵_𝑘2(𝑡), (5) where 𝐵_𝑘𝑖(𝑡), 𝑘=1,2,…,𝑛, 𝑖=1,2, are independent standard Brownian motions with 𝐵_𝑘𝑖(0)=0, and 𝜎²_𝑘𝑖≥0, 𝑘=1,2,…,𝑛, 𝑖=1,2, represent the intensities of 𝐵_𝑘𝑖(𝑡).

The remaining parts of this paper are as follows. In the next section we show the existence and uniqueness of a global positive solution of model (5). In Section 3, we analyze the asymptotic behavior around the disease-free equilibrium. In Section 4, we study the dynamic of system (5) around the endemic of the deterministic model. In Section 5, numerical simulations and conclusions are carried out.

2. Global Positive Solution

When we study a dynamical behavior, a global solution is important for the system. In this section we show that the solution of system (5) is global and nonnegative. As we know, for a stochastic differential equation, the coefficients of the equation are generally required to satisfy the linear growth condition and the local Lipschitz condition. It is a sufficient condition for a stochastic differential equation has a unique global (i.e., no explosion in a finite time) solution for any given initial value [27, 28]. Although the coefficients of system (5) satisfy locally Lipschitz continuous, they are not satisfied with the linear growth condition, so the solution of system (5) may explode at a finite time. In this section, Lyapunov analysis method (mentioned in [29]) is used to show that the solution of system (5) is positive and global.

Theorem 1. For any given initial value (𝑆₁(0),𝐼₁(0),…,𝑆_𝑛(0),𝐼_𝑛(0))∈ℝ₊^2𝑛, there is a unique positive solution (𝑆₁(𝑡),𝐼₁(𝑡),…,𝑆_𝑛(𝑡),𝐼_𝑛(𝑡)) of model (5) on 𝑡≥0 and the solution will remain in ℝ₊^2𝑛 with probability 1, namely, (𝑆₁(𝑡),𝐼₁(𝑡),…,𝑆_𝑛(𝑡),𝐼_𝑛(𝑡))∈ℝ₊^2𝑛 for all 𝑡≥0 a.s.

Proof. Due to the fact that the coefficients of the system (5) are locally Lipschitz continuous, for any given initial value (𝑆₁(0),𝐼₁(0),…,𝑆_𝑛(0),𝐼_𝑛(0))∈ℝ₊^2𝑛, it has a unique local solution (𝑆₁(𝑡),𝐼₁(𝑡),…,𝑆_𝑛(𝑡),𝐼_𝑛(𝑡)) on 𝑡∈[0,𝜏_𝑒), where 𝜏_𝑒 is the explosion time [30]. If we show that 𝜏_𝑒=∞ a.s., it suggests that this solution is global. Let 𝑙₀>0 be sufficiently large so that 𝑆_𝑘(0), 𝐼_𝑘(0)   (𝑘=1,2,…,𝑛) all lie within the interval [1/𝑙₀,𝑙₀]. For each integer 𝑙≥𝑙₀, defining the stopping time 𝜏_𝑙=inf𝑡∈0,𝜏_𝑒∶min𝑆_𝑘(𝑡),𝐼_𝑘(𝑡)≤1,𝑘=1,…,𝑛𝑙ormax𝑆_𝑘(𝑡),𝐼_𝑘(𝑡),𝑘=1,…,𝑛≥𝑙, (6) we set inf∅=∞ (as usual ∅ denotes the empty set). Obviously, 𝜏_𝑙 is increasing as 𝑙→∞. Set 𝜏_∞=lim_𝑙→∞𝜏_𝑙; therefore 𝜏_∞≤𝜏_𝑒 a.s. If 𝜏_∞=∞ a.s. is true, then 𝜏_𝑒=∞ a.s. and (𝑆₁(𝑡),𝐼₁(𝑡),…,𝑆_𝑛(𝑡),𝐼_𝑛(𝑡))∈ℝ₊^2𝑛 a.s. for 𝑡≥0. In other words, to complete the proof it is required to show that 𝜏_∞=∞ a.s. If this statement is false, then there is a pair of constants 𝑇>0 and 𝜀∈(0,1) such that 𝑃{𝜏_∞≤𝑇}>𝜀. Thus there is an integer 𝑙₁≥𝑙₀, such that 𝑃𝜏_𝑙≤𝑇≥𝜀,∀𝑙≥𝑙₁. (7) Define a 𝐶²-function 𝑉∶ℝ₊^2𝑛→ℝ₊ as follows: 𝑉𝑆₁,𝐼₁,…,𝑆_𝑛,𝐼_𝑛=^𝑛_𝑘=1𝑆_𝑘−1−ln𝑆_𝑘+𝐼_𝑘−1−ln𝐼_𝑘. (8) Applying Itô’s formula, we obtain 𝑑𝑆₁,𝐼₁,…,𝑆_𝑛,𝐼_𝑛=^𝑛_𝑘=111−𝑆_𝑘𝑑𝑆_𝑘+121𝑆²_𝑘𝑑𝑆_𝑘²+11−𝐼_𝑘𝑑𝐼_𝑘+121𝐼²_𝑘𝑑𝐼_𝑘²=^𝑛_𝑘=111−𝑆_𝑘𝑏_𝑘−^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝐼_𝑗−𝑑𝑆_𝑘𝑑𝑡+𝜎_𝑘1𝑆_𝑘(𝑡)𝑑𝐵_𝑘1(𝑡)+12𝜎²_𝑘11𝑑𝑡+1−𝐼_𝑘×^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝐼_𝑗−(𝑑+𝜖)𝐼_𝑘𝑑𝑡+𝜎_𝑘2𝐼_𝑘𝑑𝐵_𝑘2(𝑡)+12𝜎²_𝑘2=𝑑𝑡_𝑛𝑘=1𝑏_𝑘−∑^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝐼_𝑗−𝑑𝑆_𝑘−𝑏_𝑘𝑆_𝑘+∑^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝐼_𝑗+1+𝑑2𝜎²_𝑘1+^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝐼_𝑗−(𝑑+𝜖)𝐼_𝑘−∑^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝐼_𝑗𝐼_𝑘+()+1𝑑+𝜖2𝜎²_𝑘2+𝑑𝑡_𝑛𝑘=1𝜎_𝑘1𝑆_𝑘−1𝑑𝐵_𝑘1(𝑡)+𝜎_𝑘2𝐼_𝑘𝑑𝐵_𝑘2(𝑡≤_𝑛𝑘=1)𝑏_𝑘+2𝑑+𝜖+^𝑛_𝑗=1𝛽_𝑘𝑗𝐼_𝑗+12𝜎²_𝑘1+𝜎²_𝑘2+𝑑𝑡_𝑛𝑘=1𝜎_𝑘1𝑆_𝑘−1𝑑𝐵_𝑘1(𝑡)+𝜎_𝑘2𝐼_𝑘−1𝑑𝐵_𝑘2(𝑡𝜎)=∶𝐾𝑑𝑡+_𝑘1𝑆_𝑘−1𝑑𝐵_𝑘1(𝑡)+𝜎_𝑘2𝐼_𝑘𝑑𝐵_𝑘2(𝑡). (9) We can now integrate both sides of (9) from 0 to 𝜏_𝑙∧𝑇 and then take the expectations 𝐸𝑉𝑆₁𝜏_𝑙∧𝑇,𝐼₁𝜏_𝑙∧𝑇,…,𝑆_𝑛𝜏_𝑙∧𝑇,𝐼_𝑛𝜏_𝑙𝑆∧𝑇≤𝑉₁(0),𝐼₁(0),…,𝑆_𝑛(0),𝐼_𝑛(0)+𝐸^𝜏𝑙₀^∧𝑇𝑆𝐾𝑑𝑡≤𝑉₁(0),𝐼₁(0),…,𝑆_𝑛(0),𝐼_𝑛(0)+𝐾𝑇. (10) Let Ω_𝑙={𝜏_𝑙≤𝑇} for 𝑙≥𝑙₁ and, by (7), 𝑃(Ω_𝑙)≥𝜀. Note that, for every 𝜔∈Ω_𝑙, there is at least one of 𝑆_𝑘(𝜏_𝑙,𝜔) and 𝐼_𝑘(𝜏_𝑙,𝜔), 𝑘=1,2,…,𝑛, that equals either 𝑙 or 1/𝑙, and therefore 𝑉(𝑆₁(𝜏_𝑙,𝜔),𝐼₁(𝜏_𝑙,𝜔),…,𝑆_𝑛(𝜏_𝑙,𝜔),𝐼_𝑛(𝜏_𝑙,𝜔)) is not less than either 1𝑙−1−ln𝑙or𝑙1−1−ln𝑙=1𝑙−1+ln𝑙. (11) Hence, 𝑉𝑆₁𝜏_𝑙,𝜔,𝐼₁𝜏_𝑙,𝜔,…,𝑆_𝑛𝜏_𝑙,𝜔,𝐼_𝑛𝜏_𝑙≥()∧1,𝜔𝑙−1−ln𝑙𝑙−1+ln𝑙. (12) It then follows from (7) and (10) that 𝑉𝑆₁(0),𝐼₁(0),…,𝑆_𝑛(0),𝐼_𝑛(0)1+𝐾𝑇≥𝐸_Ω𝑙(𝜔)𝑉𝑆₁𝜏_𝑙,𝜔,𝐼₁𝜏_𝑙,𝜔,…,𝑆_𝑛𝜏_𝑙,𝜔,𝐼_𝑛𝜏_𝑙)∧1,𝜔≥𝜀(𝑙−1−ln𝑙𝑙−1+ln𝑙, (13) where 1_Ω𝑙(𝜔) is the indicator function of Ω_𝑙. Letting 𝑙→∞, we have that ∞>𝑉𝑆₁(0),𝐼₁(0),…,𝑆_𝑛(0),𝐼_𝑛(0)+𝐾𝑇≥∞ (14) is a contradiction. So we must have 𝜏_∞=∞. Therefore, it implies 𝑆_𝑘(𝑡), 𝐼_𝑘(𝑡), 𝑘=1,2,…,𝑛, will not explode in a finite time with probability one.

3. Asymptotic Behavior around the Disease-Free Equilibrium

As mentioned in the Introduction, 𝐸₀=(𝑏₁/𝑑,0,…,𝑏_𝑛/𝑑,0) is the disease-free equilibrium of system (3), and when 𝑅₀≤1, 𝐸₀ is globally stable, which means that the disease will be extinct in the limited time. In this section, we will study the asymptotic behavior around 𝐸₀ of system (5).

Lemma 2. If 𝐴 is nonnegative and irreducible, then the spectral radius 𝜌(𝐴) of 𝐴 is a simple eigenvalue, and 𝐴 has a positive eigenvector 𝜔=(𝜔₁,𝜔₂,…,𝜔_𝑛) corresponding to 𝜌(𝐴). Besides, if 0≤𝐴≤𝐵, then 𝜌(𝐴)≤𝜌(𝐵). (This lemma can be found in [20].)

Theorem 3. Assume 𝐴=(𝛽_𝑘𝑗)_𝑛×𝑛 is irreducible. If 𝑅₀≤1 and the following condition is satisfied: 𝜎²_𝑘1≤43𝑑,𝜎²_𝑘2≤2(𝑑+𝜖), (15) then for any given initial value (𝑆₁(0),𝐼₁(0),…,𝑆_𝑛(0),𝐼_𝑛(0))∈ℝ₊^2𝑛, the solution of system (5) has the property limsup_𝑡→∞1𝑡𝐸^𝑡₀𝑆_𝑘(𝑟)−𝑏_𝑘𝑑²+𝐼²_𝑘(𝑟≤3)𝑑𝑟2^𝑛_𝑘∑₌₁𝑎_𝑘𝜎²_𝑘1𝑏²_𝑘𝑑²𝐾₁, (16) where 𝐾₁𝜔=min^𝑘𝛽_𝑘𝑘()32𝑑+𝜖)(𝑑+𝜖𝑑−4𝜎²_𝑘11,𝑑+𝜖−2𝜎²_𝑘2. (17)

Proof. First change the variables 𝑠_𝑘=𝑆_𝑘−𝑏_𝑘/𝑑, 𝑖_𝑘=𝐼_𝑘; then −𝑏_𝑘/𝑑≤𝑠_𝑘≤0, 𝑖_𝑘≥0 and system (5) can be written as 𝑑𝑠_𝑘=−^𝑛_𝑗=1𝛽_𝑘𝑗𝑠_𝑘+𝑏_𝑘𝑑𝑖_𝑗−𝑑𝑠_𝑘𝑑𝑡+𝜎_𝑘1𝑠_𝑘+𝑏_𝑘𝑑𝑑𝐵_𝑘1(𝑡),𝑑𝑖_𝑘=^𝑛_𝑗=1𝛽_𝑘𝑗𝑠_𝑘+𝑏_𝑘𝑑𝑖_𝑗−(𝑑+𝜖)𝑖_𝑘𝑑𝑡+𝜎_𝑘2𝑖_𝑘𝑑𝐵_𝑘2(𝑡). (18) Let 𝑆⁰=(𝑆⁰₁,𝑆⁰₂,…,𝑆⁰_𝑛), where 𝑆⁰_𝑘=𝑏_𝑘/𝑑, 𝑘=1,2,…,𝑛. Define 𝑀(𝑠)=𝛽₁₁𝑆⁰₁𝛽𝑑+𝜖¹²𝑆⁰₁+𝑠₁⋯𝛽𝑑+𝜖^1𝑛𝑆⁰₁+𝑠₁𝛽𝑑+𝜖₂₁𝑆⁰₂+𝑠₂𝛽𝑑+𝜖²²𝑆⁰₂⋯𝛽𝑑+𝜖^2𝑛𝑆⁰₂+𝑠₂𝛽𝑑+𝜖⋮⋮⋱⋮_𝑛1𝑆⁰_𝑛+𝑠_𝑛𝛽𝑑+𝜖^𝑛2𝑆⁰_𝑛+𝑠_𝑛⋯𝛽𝑑+𝜖^𝑛𝑛𝑆⁰_𝑛𝑑+𝜖; (19) then it is nonnegative and irreducible. By Lemma 2, there is a positive eigenvector 𝜔=(𝜔₁,𝜔₂,…,𝜔_𝑛) of 𝑀(𝑠) corresponding to 𝜌(𝑀(𝑠)), such that 𝜔₁,𝜔₂,…,𝜔_𝑛𝑀(𝑠)=𝜔₁,𝜔₂,…,𝜔_𝑛𝜌(𝑀(𝑠)). (20) Define a 𝐶²-function 𝑉∶ℝ₊^2𝑛→ℝ₊ by 𝑉𝑠₁,𝑖₁,…,𝑠_𝑛,𝑖_𝑛=12^𝑛_𝑘=1𝑎_𝑘𝑠_𝑘+𝑖_𝑘²+^𝑛_𝑘=1𝜔_𝑘𝑖𝑑+𝜖^𝑘, (21) where 𝑎_𝑘, 𝑘=1,2,…,𝑛, are positive constants. Then the function 𝑉 is positive definite, and 𝑑𝑉=𝐿𝑉𝑑𝑡+^𝑛_𝑘=1𝑎_𝑘𝑠_𝑘+𝑖_𝑘×𝜎_𝑘1𝑠_𝑘+𝑏_𝑘𝑑𝑑𝐵_𝑘1(𝑡)+𝜎_𝑘2𝑖_𝑘𝑑𝐵_𝑘2(𝑡)+^𝑛_𝑘=1𝜔_𝑘𝜎𝑑+𝜖^𝑘2𝑖_𝑘𝑑𝐵_𝑘2(𝑡), (22) where 𝐿𝑉=^𝑛_𝑘=1𝑎_𝑘𝑠_𝑘+𝑖_𝑘−𝑑𝑠_𝑘−(𝑑+𝜖)𝑖_𝑘+12^𝑛_𝑘=1𝑎_𝑘𝜎²_𝑘1𝑠_𝑘+𝑏_𝑘𝑑²+𝜎²_𝑘2𝑖²_𝑘+^𝑛_𝑘=1𝜔_𝑘𝑑+𝜖^𝑛_𝑗=1𝛽_𝑘𝑗𝑠_𝑘+𝑏_𝑘𝑑𝑖_𝑗−(𝑑+𝜖)𝑖_𝑘=−^𝑛_𝑘=1𝑎_𝑘1𝑑−2𝜎²_𝑘1𝑠²_𝑘+1𝑑+𝜖−2𝜎²_𝑘2𝑖²_𝑘+^𝑛_{𝑘=1𝑗≠𝑘}𝜔_𝑘𝛽𝑑+𝜖^𝑘𝑗𝑠_𝑘𝑖_𝑘−^𝑛_𝑘=1𝑎_𝑘()−𝜔2𝑑+𝜖^𝑘𝛽𝑑+𝜖^𝑘𝑘𝑠_𝑘𝑖_𝑘+^𝑛𝑛_𝑘=1𝑗=1𝜔_𝑘𝛽𝑑+𝜖^𝑘𝑗𝑏_𝑘𝑑𝑖_𝑘−^𝑛_𝑘=1𝜔_𝑘𝑖_𝑘+12^𝑛_𝑘=1𝑎_𝑘𝜎²_𝑘12𝑠_𝑘𝑏_𝑘𝑑+𝑏²_𝑘𝑑². (23) Choose  𝑎_𝑘=𝜔_𝑘𝛽_𝑘𝑘/(2𝑑+𝜖)(𝑑+𝜖),  𝑘=1,2,…,𝑛; then 𝑎_𝑘(2𝑑+𝜖)−(𝜔_𝑘/(𝑑+𝜖))𝛽_𝑘𝑘=0. And we note that ^𝑛_{𝑘=1𝑗≠𝑘}𝜔_𝑘𝛽𝑑+𝜖^𝑘𝑗𝑠_𝑘𝑖_𝑘+^𝑛𝑛_𝑘=1𝑗=1𝜔_𝑘𝛽𝑑+𝜖^𝑘𝑗𝑏_𝑘𝑑𝑖_𝑘=𝜔₁,𝜔₂,…,𝜔_𝑛×𝛽₁₁𝑆⁰₁𝛽𝑑+𝜖¹²𝑆⁰₁+𝑠₁⋯𝛽𝑑+𝜖^1𝑛𝑆⁰₁+𝑠₁𝛽𝑑+𝜖₂₁𝑆⁰₂+𝑠₂𝛽𝑑+𝜖²²𝑆⁰₂⋯𝛽𝑑+𝜖^2𝑛𝑆⁰₂+𝑠₂𝛽𝑑+𝜖⋮⋮⋱⋮_𝑛1𝑆⁰_𝑛+𝑠_𝑛𝛽𝑑+𝜖^𝑛2𝑆⁰_𝑛+𝑠_𝑛⋯𝛽𝑑+𝜖^𝑛𝑛𝑆⁰_𝑛×𝑖𝑑+𝜖₁𝑖₂⋮𝑖_𝑛=𝜔𝑀(𝑠)𝑖. (24) Then 𝐿𝑉=−^𝑛_𝑘=1𝑎_𝑘1𝑑−2𝜎²_𝑘1𝑠²_𝑘+1𝑑+𝜖−2𝜎²_𝑘2𝑖²_𝑘(𝑠)1−𝜔𝑖+𝜔𝑀𝑖+2^𝑛_𝑘=1𝑎_𝑘𝜎²_𝑘12𝑠_𝑘𝑏_𝑘𝑑+𝑏²_𝑘𝑑²≤−^𝑛_𝑘=1𝑎_𝑘3𝑑−4𝜎²_𝑘1𝑠²_𝑘+1𝑑+𝜖−2𝜎²_𝑘2𝑖²_𝑘(𝑀(𝑠))3+𝜔−1𝑖+2^𝑛_𝑘=1𝑎_𝑘𝜎²_𝑘1𝑏²_𝑘𝑑²=−^𝑛_𝑘=1𝑎_𝑘3𝑑−4𝜎²_𝑘1𝑠²_𝑘+1𝑑+𝜖−2𝜎²_𝑘2𝑖²_𝑘𝜌(𝑀(𝑠3+𝜔))−1𝑖+2^𝑛_𝑘=1𝑎_𝑘𝜎²_𝑘1𝑏²_𝑘𝑑², (25) where the last equality is derived from (20). Since −𝑏_𝑘/𝑑≤𝑠_𝑘≤0, then 0≤𝑀(𝑠)≤𝑀(𝑆⁰)=(𝑆⁰_𝑘𝛽_𝑘𝑗/(𝑑+𝜖))_𝑛×𝑛=𝑀₀, and so 𝜌(𝑀(𝑠))≤𝜌(𝑀₀) according to Lemma 2. Besides, 𝑅₀≤1, and then 𝜌(𝑀(𝑠))≤1. Therefore 𝑑𝑉≤−^𝑛_𝑘=1𝑎_𝑘3𝑑−4𝜎²_𝑘1𝑠²_𝑘+1𝑑+𝜖−2𝜎²_𝑘2𝑖²_𝑘+32^𝑛_𝑘=1𝑎_𝑘𝜎²_𝑘1𝑏²_𝑘𝑑²+𝑑𝑡_𝑛𝑘=1𝑎_𝑘𝑠_𝑘+𝑖_𝑘𝜎_𝑘1𝑠_𝑘+𝑏_𝑘𝑑𝑑𝐵_𝑘1(𝑡)+^𝑛_𝑘=1𝑎_𝑘𝑠_𝑘+𝑖_𝑘𝜎_𝑘2+𝜔_𝑘𝜎𝑑+𝜖^𝑘2𝑖_𝑘𝑑𝐵_𝑘2(𝑡). (26) Integrating both sides of (26) from 0 to 𝑡, and taking expectation, yields 0≤𝐸𝑉𝑠₁(𝑡),𝑖₁(𝑡),…,𝑠_𝑛(𝑡),𝑖_𝑛(𝑡)≤𝐸𝑉𝑠₁(0),𝑖₁(0),…,𝑠_𝑛(0),𝑖_𝑛(0)+𝐸^𝑡₀−^𝑛_𝑘=1𝑎_𝑘3𝑑−4𝜎²_𝑘1𝑠²_𝑘(𝑟)+1𝑑+𝜖−2𝜎²_𝑘2𝑖²_𝑘(𝑟+3)2^𝑛_𝑘=1𝑎_𝑘𝜎²_𝑘1𝑏²_𝑘𝑑²𝑑𝑟, (27) which implies 𝐸^𝑡₀^𝑛_𝑘=1𝑎_𝑘3𝑑−4𝜎²_𝑘1𝑠²_𝑘(𝑟)+1𝑑+𝜖−2𝜎²_𝑘2𝑖²_𝑘(𝑟𝑉𝑠)𝑑𝑟≤𝐸₁(0),𝑖₁(0),…,𝑠_𝑛(0),𝑖_𝑛(0)+32^𝑛_𝑘=1𝑎_𝑘𝜎²_𝑘1𝑏²_𝑘𝑑²𝑡. (28) Therefore limsup_𝑡→∞1𝑡𝐸^𝑡₀^𝑛_𝑘=1𝑎_𝑘3𝑑−4𝜎²_𝑘1𝑠²_𝑘(𝑟)+1𝑑+𝜖−2𝜎²_𝑘2𝑖²_𝑘(𝑟≤3)𝑑𝑟2^𝑛_𝑘∑₌₁𝑎_𝑘𝜎²_𝑘1𝑏²_𝑘𝑑²; (29)
that is, limsup_𝑡→∞1𝑡𝐸^𝑡₀^𝑛_𝑘=1𝑎_𝑘3𝑑−4𝜎²_𝑘1𝑆_𝑘(𝑟)−𝑏_𝑘𝑑²+1𝑑+𝜖−2𝜎²_𝑘2𝐼²_𝑘(𝑟≤3)𝑑𝑟2^𝑛_𝑘=1𝑎_𝑘𝜎²_𝑘1𝑏²_𝑘𝑑². (30) If we let 𝐾₁=min𝑎_𝑘3𝑑−4𝜎²_𝑘11,𝑑+𝜖−2𝜎²_𝑘2, (31) then limsup_𝑡→∞1𝑡𝐸^𝑡₀𝑆_𝑘(𝑟)−𝑏_𝑘𝑑²+𝐼²_𝑘(𝑟≤3)𝑑𝑟2^𝑛_𝑘∑₌₁𝑎_𝑘𝜎²_𝑘1𝑏²_𝑘𝑑²𝐾₁, (32) as the theorem is proved.

Remark 4. From Theorem 3, we can get the conclusion that the solution of the stochastic system will oscillate around the disease-free equilibrium of the deterministic model; the values of 𝜎_𝑘1 and 𝜎_𝑘2 have bearing on the intensity of turbulence. If the stochastic perturbations become small, the solution of system (5) will be close to the disease-free equilibrium of system (3).
Besides, if 𝜎_𝑘1=0, then 𝐸₀ is also the disease-free equilibrium of system (5). From the proof of Theorem 3, we can obtain 𝐿𝑉≤−^𝑛_𝑘=1𝑎_𝑘3𝑑−4𝜎²_𝑘1𝑠²_𝑘+1𝑑+𝜖−2𝜎²_𝑘2𝑖²_𝑘≤0. (33) Therefore, 𝐸₀ is globally asymptotically stable.

4. The Dynamic of System (5) around the Endemic of System (3)

In the deterministic model, if 𝑅₀>1, there exists the endemic equilibrium 𝐸^∗. But 𝐸^∗ is not the endemic equilibrium of stochastic system (5), because there is no endemic equilibrium for the stochastic system (5). In fact, we still want to find the relation between the solution of stochastic system and 𝐸^∗.

Given a weighted digraph (𝒢,𝐴) with 𝑛 vertices, where 𝐴=(𝑎_𝑘𝑗)_𝑛𝑛 is the weight matrix, whose entry 𝑎_𝑘𝑗 equals the weight of arc (𝑗,𝑘) if it exists, and 0 otherwise, the Laplacian matrix of 𝐴 is defined as 𝐿_𝐴=_𝑘∑_≠1𝑎_1𝑘−𝑎₁₂⋯−𝑎_1𝑛−𝑎₂₁∑_𝑘≠2𝑎_2𝑘⋯−𝑎_2𝑛⋮⋮⋱⋮−𝑎_𝑛1−𝑎_𝑛2⋯∑_𝑘≠𝑛𝑎_𝑛𝑘. (34) Let 𝑐_𝑘 denote the cofactor of the 𝑘th diagonal element of 𝐿_𝐴, and we have the following results.

Theorem 5. Assume 𝐴=(𝛽_𝑘𝑗)_𝑛×𝑛 is irreducible and 𝑅₀>1. For any given initial value (𝑆₁(0),𝐼₁(0),…,𝑆_𝑛(0),𝐼_𝑛(0))∈ℝ₊^2𝑛, the solution of system (5) has the property limsup_𝑡→∞1𝑡^𝑛_𝑘=1^𝑡₀𝑝_𝑘𝑑𝑆_𝑘−𝑆^∗_𝑘²+𝑚_𝑘(𝑑+𝜖)𝐼_𝑘−𝐼^∗_𝑘²≤𝑑𝑠_𝑛𝑘=1𝑎𝑐_𝑘𝑏_𝑘𝑑+𝑚_𝑘+𝑝_𝑘𝑏²_𝑘𝑑²𝜎²_𝑘1+(𝑎+1)𝑐_𝑘𝑏_𝑘𝑑+𝑚_𝑘𝑏²_𝑘𝑑²𝜎²_𝑘2,𝑎.𝑠., (35) where 𝐸^∗=(𝑆^∗₁,𝐼^∗₁,…,𝑆^∗_𝑛,𝐼^∗_𝑛) is the endemic equilibrium of system (3) and 𝑐_𝑘, 𝑘=1,2,…,𝑛, denote the cofactor of the 𝑘th diagonal element of 𝐿_𝐴(𝐴=(𝛽_𝑘𝑗)_𝑛×𝑛=(𝛽_𝑘𝑗𝑆^∗_𝑘𝐼^∗_𝑗)_𝑛×𝑛), and 𝑎, 𝑚_𝑘, 𝑝_𝑘, 𝑘=1,2,…,𝑛, are positive constants defined as in the proof.

Proof. Since 𝐸^∗ is the endemic equilibrium of system (3), we have ^𝑛_𝑗=1𝛽_𝑘𝑗𝑆^∗_𝑘𝐼^∗_𝑗+𝑑𝑆^∗_𝑘=𝑏_𝑘,^𝑛_𝑗=1𝛽_𝑘𝑗𝑆^∗_𝑘𝐼^∗_𝑗=(𝑑+𝜖)𝐼^∗_𝑘. (36) Define 𝑉𝑆₁,𝐼₁,…,𝑆_𝑛,𝐼_𝑛=𝑎^𝑛_𝑘=1𝑐_𝑘𝑆_𝑘−𝑆^∗_𝑘−𝑆^∗_𝑘𝑆ln^𝑘𝑆^∗_𝑘+𝐼_𝑘−𝐼^∗_𝑘−𝐼^∗_𝑘𝐼ln^𝑘𝐼^∗_𝑘+^𝑛_𝑘=1𝑐_𝑘𝐼_𝑘−𝐼^∗_𝑘−𝐼^∗_𝑘𝐼ln^𝑘𝐼^∗_𝑘+12^𝑛_𝑘=1𝑚_𝑘𝑆_𝑘−𝑆^∗_𝑘+𝐼_𝑘−𝐼^∗_𝑘²+12^𝑛_𝑘=1𝑝_𝑘𝑆_𝑘−𝑆^∗_𝑘²∶=𝑎𝑉₁+𝑉₂+𝑉₃+𝑉₄, (37) where 𝑎, 𝑚_𝑘, 𝑝_𝑘, 𝑘=1,2,…,𝑛, are positive constants to be determined later. From the property (1) of Lemma A.2 (see [20]), we know 𝑐_𝑘>0, 𝑘=1,2,…,𝑛. Hence 𝑉 is positive definite. Let 𝐿 be the generating operator of system (5). Then we get 𝐿𝑉₁=^𝑛_𝑘=1𝑐_𝑘𝑆1−^∗𝑘𝑆_𝑘𝑏_𝑘−^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝐼_𝑗−𝑑𝑆_𝑘+^𝑛_𝑘=1𝑐_𝑘𝑆^∗_𝑘𝜎²_𝑘12+^𝑛_𝑘=1𝑐_𝑘𝐼1−^∗𝑘𝐼_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝐼_𝑗−(𝑑+𝜖)𝐼_𝑘+^𝑛_𝑘=1𝑐_𝑘𝐼^∗_𝑘𝜎²_𝑘22=^𝑛_𝑘=1𝑐_𝑘𝑏_𝑘−𝑑𝑆_𝑘−𝑏_𝑘𝑆^∗_𝑘𝑆_𝑘+^𝑛_𝑗=1𝛽_𝑘𝑗𝑆^∗_𝑘𝐼_𝑗+𝑑𝑆^∗_𝑘−(𝑑+𝜖)𝐼_𝑘+(𝑑+𝜖)𝐼^∗_𝑘−^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝐼_𝑗𝐼^∗_𝑘𝐼_𝑘+12𝑆^∗_𝑘𝜎²_𝑘1+𝐼^∗_𝑘𝜎²_𝑘2=^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗+𝑑𝑆^∗_𝑘−𝑑𝑆_𝑘−^𝑛_𝑗=1𝛽_𝑘𝑗𝑆^∗_𝑘𝑆_𝑘−𝑑𝑆^∗_𝑘𝑆^∗_𝑘𝑆_𝑘+^𝑛_𝑗=1𝛽_𝑘𝑗𝐼_𝑗𝐼^∗_𝑗+𝑑𝑆^∗_𝑘−^𝑛_𝑗=1𝛽_𝑘𝑗𝐼_𝑘𝐼^∗_𝑘−^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝐼_𝑗𝐼^∗_𝑘𝐼_𝑘𝑆^∗_𝑘𝐼^∗_𝑗+^𝑛_𝑗=1𝛽_𝑘𝑗+12𝑆^∗_𝑘𝜎²_𝑘1+𝐼^∗_𝑘𝜎²_𝑘2=^𝑛_𝑘=1𝑐_𝑘−𝑑𝑆^∗_𝑘𝑆^∗_𝑘𝑆_𝑘+𝑆_𝑘𝑆^∗_𝑘−2+^𝑛_𝑗=1𝛽_𝑘𝑗𝐼_𝑗𝐼^∗_𝑗−𝐼_𝑘𝐼^∗_𝑘+^𝑛_𝑗=1𝛽_𝑘𝑗𝑆2−^∗𝑘𝑆_𝑘−𝑆_𝑘𝐼_𝑗𝐼^∗_𝑘𝐼_𝑘𝑆^∗_𝑘𝐼^∗_𝑗+12𝑆^∗_𝑘𝜎²_𝑘1+𝐼^∗_𝑘𝜎²_𝑘2.𝐿𝑉₂=^𝑛_𝑘=1𝑐_𝑘𝐼1−^∗𝑘𝐼_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝐼_𝑗−(𝑑+𝜖)𝐼_𝑘+^𝑛_𝑘=1𝑐_𝑘𝐼^∗_𝑘𝜎²_𝑘22=^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝐼_𝑗−(𝑑+𝜖)𝐼_𝑘−𝐼^∗_𝑘𝐼_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝐼_𝑗+(𝑑+𝜖)𝐼^∗_𝑘+𝐼^∗_𝑘𝜎²_𝑘22=^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝐼_𝑗−^𝑛_𝑗=1𝛽_𝑘𝑗𝐼_𝑘𝐼^∗_𝑘−^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝐼_𝑗𝐼^∗_𝑘𝐼_𝑘𝑆^∗_𝑘𝐼^∗_𝑗+^𝑛_𝑗=1𝛽_𝑘𝑗+𝐼^∗_𝑘𝜎²_𝑘22.𝐿𝑉₃=^𝑛_𝑘=1𝑚_𝑘𝑆_𝑘−𝑆^∗_𝑘+𝐼_𝑘−𝐼^∗_𝑘𝑏_𝑘−𝑑𝑆_𝑘−(𝑑+𝜖)𝐼_𝑘+^𝑛_𝑘=1𝑚_𝑘𝜎²_𝑘1𝑆²_𝑘+𝜎²_𝑘2𝐼²_𝑘2=^𝑛_𝑘=1𝑚_𝑘𝑆_𝑘−𝑆^∗_𝑘+𝐼_𝑘−𝐼^∗_𝑘×−𝑑𝑆_𝑘−𝑆^∗_𝑘−(𝑑+𝜖)𝐼_𝑘−𝐼^∗_𝑘+^𝑛_𝑘=1𝑚_𝑘𝜎²_𝑘1𝑆²_𝑘+𝜎²_𝑘2𝐼²_𝑘2=^𝑛_𝑘=1𝑚_𝑘−𝑑𝑆_𝑘−𝑆^∗_𝑘²−(𝑑+𝜖)𝐼_𝑘−𝐼^∗_𝑘²−(2𝑑+𝜖)𝑆_𝑘−𝑆^∗_𝑘𝐼_𝑘−𝐼^∗_𝑘+𝜎²_𝑘1𝑆²_𝑘+𝜎²_𝑘2𝐼²_𝑘2≤^𝑛_𝑘=1𝑚_𝑘−()𝑑−2𝑑+𝜖²2()𝑆𝑑+𝜖^𝑘−𝑆^∗_𝑘²−(𝑑+𝜖)2𝐼_𝑘−𝐼^∗_𝑘²+𝜎²_𝑘1𝑆²_𝑘+𝜎²_𝑘2𝐼²_𝑘2=^𝑛_𝑘=1𝑚_𝑘(𝑑+𝜖)²+𝑑²2()𝑆𝑑+𝜖^𝑘−𝑆^∗_𝑘²−(𝑑+𝜖)2𝐼_𝑘−𝐼^∗_𝑘²+𝜎²_𝑘1𝑆²_𝑘+𝜎²_𝑘2𝐼²_𝑘2.𝐿𝑉₄=^𝑛_𝑘=1𝑝_𝑘𝑆_𝑘−𝑆^∗_𝑘𝑏_𝑘−^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝐼_𝑗−𝑑𝑆_𝑘+^𝑛_𝑘=1𝑝_𝑘𝜎²_𝑘1𝑆²_𝑘2=^𝑛_𝑘=1𝑝_𝑘𝑆_𝑘−𝑆^∗_𝑘×^𝑛_𝑗=1𝛽_𝑘𝑗𝑆^∗_𝑘𝐼^∗_𝑗−𝑆_𝑘𝐼_𝑗−𝑑𝑆_𝑘−𝑆^∗_𝑘+^𝑛_𝑘=1𝑝_𝑘𝜎²_𝑘1𝑆²_𝑘2=−^𝑛_𝑘=1𝑝_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆^∗_𝑘𝑆_𝑘−𝑆^∗_𝑘𝐼_𝑗−𝐼^∗_𝑗+^𝑛_𝑗=1𝛽_𝑘𝑗𝐼_𝑗𝑆_𝑘−𝑆^∗_𝑘²+𝑑𝑆_𝑘−𝑆^∗_𝑘²−𝜎²_𝑘1𝑆²_𝑘2. (38)(39)(40) By property (2) of Lemma A.2 (see [20]), we know ^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝐼_𝑗𝐼^∗_𝑗−^𝑛_𝑗=1𝛽_𝑘𝑗𝐼_𝑘𝐼^∗_𝑘=0,_𝑛𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝐼ln^𝑗𝐼^∗_𝑗−^𝑛_𝑗=1𝛽_𝑘𝑗𝐼ln^𝑘𝐼^∗_𝑘=0. (41) Besides, note that 𝑥−1−ln𝑥≥0 for 𝑥>0; then 𝑆^∗_𝑘𝑆_𝑘𝑆≥1+ln^∗𝑘𝑆_𝑘,𝑆_𝑘𝐼_𝑗𝐼^∗_𝑘𝐼_𝑘𝑆^∗_𝑘𝐼^∗_𝑗𝑆≥1+ln^𝑘𝐼_𝑗𝐼^∗_𝑘𝐼_𝑘𝑆^∗_𝑘𝐼^∗_𝑗. (42) According to (41) and (42), we get ^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆2−^∗𝑘𝑆_𝑘−𝑆_𝑘𝐼_𝑗𝐼^∗_𝑘𝐼_𝑘𝑆^∗_𝑘𝐼^∗_𝑗≤^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆2−1+ln^∗𝑘𝑆_𝑘−𝑆1+ln^𝑘𝐼_𝑗𝐼^∗_𝑘𝐼_𝑘𝑆^∗_𝑘𝐼^∗_𝑗=^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝐼ln^𝑘𝐼^∗_𝑘𝐼−ln^𝑗𝐼^∗_𝑗=0,_𝑛𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝐼_𝑗𝐼^∗_𝑘𝐼_𝑘𝑆^∗_𝑘𝐼^∗_𝑗≥^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆1+ln^𝑘𝐼_𝑗𝐼^∗_𝑘𝐼_𝑘𝑆^∗_𝑘𝐼^∗_𝑗=^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆1−ln^∗𝑘𝑆_𝑘+^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝐼ln^𝑗𝐼^∗_𝑗𝐼−ln^𝑘𝐼^∗_𝑘≥^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆2−^∗𝑘𝑆_𝑘. (43)(44) Substituting (41) and (43) into (38), we get 𝐿𝑉₁≤^𝑛_𝑘=1𝑐_𝑘−𝑑𝑆^∗_𝑘𝑆^∗_𝑘𝑆_𝑘+𝑆_𝑘𝑆^∗_𝑘+1−22𝑆^∗_𝑘𝜎²_𝑘1+𝐼^∗_𝑘𝜎²_𝑘2=−^𝑛_𝑘=1𝑐_𝑘𝑑𝑆_𝑘−𝑆^∗_𝑘²𝑆_𝑘+^𝑛_𝑘=1𝑐_𝑘2𝑆^∗_𝑘𝜎²_𝑘1+𝐼^∗_𝑘𝜎²_𝑘2. (45) Substituting (44) into (39), we get 𝐿𝑉₂≤^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝐼_𝑗−^𝑛_𝑗=1𝛽_𝑘𝑗𝐼_𝑘𝐼^∗_𝑘−^𝑛_𝑗=1𝛽_𝑘𝑗𝑆2−^∗𝑘𝑆_𝑘+^𝑛_𝑗=1𝛽_𝑘𝑗+𝐼^∗_𝑘𝜎²_𝑘22=^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘−𝑆^∗_𝑘𝐼_𝑗−𝐼^∗_𝑗−^𝑛_𝑗=1𝛽_𝑘𝑗𝐼_𝑘𝐼^∗_𝑘−𝐼_𝑗𝐼^∗_𝑗−^𝑛_𝑗=1𝛽_𝑘𝑗𝑆2−^∗𝑘𝑆_𝑘−𝑆_𝑘𝑆^∗_𝑘+𝐼^∗_𝑘𝜎²_𝑘22=^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘−𝑆^∗_𝑘𝐼_𝑗−𝐼^∗_𝑗+^𝑛_𝑗=1𝛽_𝑘𝑗𝐼^∗_𝑗𝑆_𝑘−𝑆^∗_𝑘²𝑆_𝑘+𝐼^∗_𝑘𝜎²_𝑘22. (46) Therefore, 𝐿𝑉=𝑎𝐿𝑉₁+𝐿𝑉₂+𝐿𝑉₃+𝐿𝑉₄≤−^𝑛_𝑘=1𝑝_𝑘𝑑−𝑚_𝑘(𝑑+𝜖)²+𝑑²2()𝑆𝑑+𝜖^𝑘−𝑆^∗_𝑘²−12^𝑛_𝑘=1𝑚_𝑘(𝑑+𝜖)𝐼_𝑘−𝐼^∗_𝑘²−^𝑛_𝑘=1𝑐_𝑘𝑎𝑑−^𝑛_𝑗=1𝛽_𝑘𝑗𝐼^∗_𝑗𝑆_𝑘−𝑆^∗_𝑘²𝑆_𝑘+^𝑛𝑛_𝑘=1𝑗=1𝛽_𝑘𝑗𝑐_𝑘−𝑝_𝑘𝑆^∗_𝑘𝑆_𝑘−𝑆^∗_𝑘𝐼_𝑗−𝐼^∗_𝑗−^𝑛𝑛_𝑘=1𝑗=1𝑝_𝑘𝛽_𝑘𝑗𝐼_𝑗𝑆_𝑘−𝑆^∗_𝑘²+^𝑛_𝑘=1𝑎𝑐_𝑘𝑆^∗_𝑘+𝑚_𝑘+𝑝_𝑘𝑆²_𝑘2𝜎²_𝑘1+^𝑛_𝑘=1(𝑎+1)𝑐_𝑘𝐼^∗_𝑘+𝑚_𝑘𝐼²_𝑘2𝜎²_𝑘2≤−^𝑛_𝑘=1𝑝_𝑘𝑑−𝑚_𝑘(𝑑+𝜖)²+𝑑²2()𝑆𝑑+𝜖^𝑘−𝑆^∗_𝑘²−12^𝑛_𝑘=1𝑚_𝑘(𝑑+𝜖)𝐼_𝑘−𝐼^∗_𝑘²−^𝑛_𝑘=1𝑐_𝑘𝑎𝑑−^𝑛_𝑗=1𝛽_𝑘𝑗𝐼^∗_𝑗𝑆_𝑘−𝑆^∗_𝑘²𝑆_𝑘+^𝑛𝑛_𝑘=1𝑗=1𝛽_𝑘𝑗𝑐_𝑘−𝑝_𝑘𝑆^∗_𝑘𝑆_𝑘−𝑆^∗_𝑘𝐼_𝑗−𝐼^∗_𝑗+^𝑛_𝑘=1𝑎𝑐_𝑘𝑆^∗_𝑘+𝑚_𝑘+𝑝_𝑘𝑆²_𝑘2𝜎²_𝑘1+^𝑛_𝑘=1(𝑎+1)𝑐_𝑘𝐼^∗_𝑘+𝑚_𝑘𝐼²_𝑘2𝜎²_𝑘2. (47) Choose  𝑎=max{(∑^𝑛_𝑗=1𝛽_𝑘𝑗𝐼^∗_𝑗)/𝑑,𝑘=1,2,…,𝑛}, 𝑚_𝑘=(𝑑+𝜖)𝑝_𝑘𝑑/((𝑑+𝜖)²+𝑑²), 𝑝_𝑘=𝑐_𝑘/𝑆^∗_𝑘, 𝑘=1,2,…,𝑛; then 𝐿𝑉≤−^𝑛_𝑘=1𝑝_𝑘𝑑2𝑆_𝑘−𝑆^∗_𝑘²−12^𝑛_𝑘=1𝑚_𝑘(𝑑+𝜖)𝐼_𝑘−𝐼^∗_𝑘²+12^𝑛_𝑘=1𝑎𝑐_𝑘𝑏𝑑+𝑚_𝑘+𝑝_𝑘𝑏²𝑑²𝜎²_𝑘1+(𝑎+1)𝑐_𝑘𝑏𝑑+𝑚_𝑘𝑏²𝑑²𝜎²_𝑘2∶=𝐹(𝑡). (48) Therefore, (𝑡)+𝑑𝑉≤𝐹𝑑𝑡_𝑛𝑘=1𝑎𝑐_𝑘𝜎_𝑘1𝑆_𝑘−𝑆^∗_𝑘+𝑚_𝑘+𝑝_𝑘𝜎_𝑘1𝑆_𝑘−𝑆^∗_𝑘𝑆_𝑘×𝑑𝐵_𝑘1(𝑡)+𝑐_𝑘(𝑎+1)𝜎_𝑘2𝐼_𝑘−𝐼^∗_𝑘+𝑚_𝑘𝜎_𝑘2𝐼_𝑘−𝐼^∗_𝑘𝐼_𝑘×𝑑𝐵_𝑘2(𝑡). (49) Integrating both sides of (49) from 0 to 𝑡 yields 𝑉(𝑡)(0)≤−𝑉_𝑡0𝐹(𝑠)+𝑑𝑠_{𝑡0𝑛𝑘=1}𝑎𝑐_𝑘𝜎_𝑘1𝑆_𝑘−𝑆^∗_𝑘+𝑚_𝑘+𝑝_𝑘𝜎_𝑘1𝑆_𝑘−𝑆^∗_𝑘𝑆_𝑘×𝑑𝐵_𝑘1(𝑠)+^𝑡₀^𝑛_𝑘=1𝑐_𝑘(𝑎+1)𝜎_𝑘2𝐼_𝑘−𝐼^∗_𝑘+𝑚_𝑘𝜎_𝑘2𝐼_𝑘−𝐼^∗_𝑘𝐼_𝑘×𝑑𝐵_𝑘2(𝑠). (50) Let 𝑀₁(𝑡)∶=∫^𝑡₀∑^𝑛_𝑘=1[𝑎𝑐_𝑘𝜎_𝑘1(𝑆_𝑘−𝑆^∗_𝑘)+(𝑚_𝑘+𝑝_𝑘)𝜎_𝑘1(𝑆_𝑘−𝑆^∗_𝑘)𝑆_𝑘]𝑑𝐵_𝑘1(𝑠), 𝑀₂(𝑡)∶=∫^𝑡₀∑^𝑛_𝑘=1[𝑐_𝑘(𝑎+1)𝜎_𝑘2(𝐼_𝑘−𝐼^∗_𝑘)+𝑚_𝑘𝜎_𝑘2(𝐼_𝑘−𝐼^∗_𝑘)𝐼_𝑘]𝑑𝐵_𝑘2(𝑠), which are local continuous martingale, and 𝑀₁(0)=𝑀₂(0)=0. Moreover limsup_𝑡→∞⟨𝑀₁,𝑀₁⟩_𝑡𝑡≤8^𝑛_𝑘=1𝜎²_𝑘1𝑎²𝑐²_𝑘+(𝑚_𝑘+𝑝_𝑘)²𝑏²_𝑘𝑑²𝑏²_𝑘𝑑²<∞,limsup_𝑡→∞⟨𝑀₂,𝑀₂⟩_𝑡𝑡≤8^𝑛_𝑘=1𝜎²_𝑘2𝑐²_𝑘(𝑎+1)²+𝑚²_𝑘𝑏²_𝑘𝑑²𝑏²_𝑘𝑑²<∞. (51) By Lemma A.4 (see [20]), we obtain lim_𝑡→∞𝑀₁(𝑡)𝑡=0,lim_𝑡→∞𝑀₂(𝑡)𝑡=0a.s., (52) which together with (50) implies liminf_𝑡→∞∫^𝑡₀𝐹(𝑠)𝑑𝑠𝑡≥0a.s. (53) Consequently, limsup_𝑡→∞1𝑡^𝑛_𝑘=1^𝑡₀𝑝_𝑘𝑑𝑆_𝑘−𝑆^∗_𝑘²+𝑚_𝑘(𝑑+𝜖)𝐼_𝑘−𝐼^∗_𝑘²≤𝑑𝑠_𝑛𝑘=1𝑎𝑐_𝑘𝑏_𝑘𝑑+𝑚_𝑘+𝑝_𝑘𝑏²_𝑘𝑑²𝜎²_𝑘1+(𝑎+1)𝑐_𝑘𝑏_𝑘𝑑+𝑚_𝑘𝑏²_𝑘𝑑²𝜎²_𝑘2,a.s. (54) Thus Theorem 5 is proved.

Remark 6. Theorem 5 shows that the solution of system (5) fluctuates around the certain level which is relevant to 𝐸^∗ of system (3) and 𝜎²_𝑘1, 𝜎²_𝑘2, 𝑘=1,2,…,𝑛. The distance between the solution 𝑋(𝑡)=(𝑆₁(𝑡),𝐼₁(𝑡),…,𝑆_𝑛(𝑡),𝐼_𝑛(𝑡)) and 𝐸^∗ of system (3) has the following form: limsup_𝑡→∞1𝑡^𝑡₀𝑋(𝑠)−𝐸^∗2𝑑𝑠≤𝐶‖𝜎‖², (55) where 𝐶 is a positive constant and ‖𝜎‖²=∑^𝑛_𝑘=1(𝜎²_𝑘1+𝜎²_𝑘2). Although the solution of system (5) does not have stability as the deterministic system, we can draw a conclusion that system (5) is persistent on the basis of the result of Theorem 5, which also accounts for the fact that the disease is prevalent.

5. Simulations and Conclusions

5.1. Numerical Simulations

In order to confirm the results above, we numerically simulate the solution of system (5) with 𝑛=50, 𝑏_𝑘=0.25, 𝑑=0.3, 𝜖=0.01, and initial value 𝑆_𝑘(0)=0.5, 𝐼_𝑘(0)=0.1, 𝑘=1,2,…,50. Using Milstein’s Higer Order Method [31], we get the discretization equation: 𝑆_𝑘,𝑖+1=𝑆_𝑘,𝑖+Δ𝑡𝑏_𝑘−²_𝑗=1𝛽_𝑘𝑗𝑆_𝑘,𝑖𝐼_𝑗,𝑖−𝑑𝑆_𝑘,𝑖+𝜎_𝑘1𝑆_𝑘,𝑖√Δ𝑡𝜉_𝑘1,𝑖+𝜎²_𝑘12𝑆_𝑘,𝑖Δ𝑡𝜉²_𝑘1,𝑖,𝐼−1_𝑘,𝑖+1=𝐼_𝑘,𝑖+Δ𝑡²_𝑗=1𝛽_𝑘𝑗𝑆_𝑘,𝑖𝐼_𝑗,𝑖−(𝑑+𝜖)𝐼_𝑘,𝑖+𝜎_𝑘2𝐼_𝑘,𝑖√Δ𝑡𝜉_𝑘2,𝑖+𝜎²_𝑘22𝐼_𝑘,𝑖Δ𝑡𝜉²_𝑘2,𝑖−1, (56) where 𝑘=1,2,…,𝑛 and 𝜉_𝑘1,𝑖, 𝜉_𝑘2,𝑖, 𝑖=1,2,…,𝑁, are the independent Gaussian random variables 𝑁(0,1).

From Theorem 3 and Remark 4, it is shown that the expectations of 𝑆_𝑘(𝑡), 𝐼_𝑘(𝑡), 𝑘=1,2,…,𝑛, are converging under some conditions, and the solution of system (5) will oscillate around the disease-free equilibrium of system (3). In Figure 1, we choose parameters 𝜆=0.025,