Dynamics of a Multigroup SIR Epidemic Model with Nonlinear Incidence and Stochastic Perturbation

Abstract

We introduce stochasticity into a multigroup SIR model with nonlinear incidence. We prove that when the intensity of white noise is small, the solution of stochastic system converges weakly to a singular measure (i.e., a distribution) if ℛ₀≤1 and there exists an invariant distribution which is ergodic if ℛ₀>1. This is the same situation as the corresponding deterministic case. When the intensity of white noise is large, white noise controls this system. This means that the disease will extinct exponentially regardless of the magnitude of ℛ₀.

1. Introduction

Considering different contact patterns, distinct number of sexual partners, or different geography, and so forth, individual hosts are often divided into groups in modeling epidemic diseases. Multigroup models have been proposed in the literature to describe the transmission dynamics of infectious diseases in heterogeneous host populations. One of the earliest work on multigroup models is the seminal paper by Lajmanovich and Yorke [1] on a class of SIS multigroup models for the transmission dynamics of Gonorrhea. From then on, much research has been done on various forms of multi-group models, see, for example, [2–6]. It is well known that the global stability of the endemic equilibrium of multigroup models is a very challenging problem. Recently, Guo et al. [7, 8] and Li and Shuai [9] proposed a graph-theoretic approach to the method of global Lyapunov functions and completely solved this problem for some multigroup epidemic models. Subsequently, this approach or ideas of [7–9] were applied to the investigation into the dynamics of several classes of multigroup epidemic models (e.g., [10–13]).

Now we consider a deterministic multigroup SIR (susceptible, infective, and recovered) epidemic model with nonlinear incidence:𝑆̇_𝑘=Λ_𝑘−^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝜙_𝑗𝐼_𝑗−𝑑_𝑘𝑆_𝑘,𝐼̇_𝑘=^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝜙_𝑗𝐼_𝑗−𝜖_𝑘+𝛾_𝑘𝐼_𝑘,𝑅̇_𝑘=𝛾_𝑘𝐼_𝑘−𝛿_𝑘𝑅_𝑘,𝑘=1,2,…,𝑛. (1) Here 𝑆_𝑘(𝑡), 𝐼_𝑘(𝑡), and 𝑅_𝑘(𝑡) denote the susceptible, infective, and recovered population at time 𝑡 in the 𝑘th group, respectively, 𝑘=1,2,…,𝑛. Suppose the death rates of 𝑆_𝑘, 𝐼_𝑘, 𝑅_𝑘 in the 𝑘th group are different. The parameters in the model (1) are assumed to be positive and summarized in the following list: Λ_𝑘: influx of individuals into the 𝑘th group; 𝛽_𝑘𝑗: transmission coefficient between compartments 𝑆_𝑘 and 𝐼_𝑗; 𝑑_𝑘: death rate of 𝑆 compartment in the 𝑘th group; 𝜖_𝑘: death rate of 𝐼 compartment in the 𝑘th group; 𝛿_𝑘: death rate of 𝑅 compartment in the 𝑘th group; 𝛾_𝑘: recovery rate of infectious individuals in the 𝑘th group.

Considering that the death rates of infective and recovered population are usually no less than the susceptible’s, we assume 𝑑_𝑘≤min{𝜖_𝑘,𝛿_𝑘} for all 𝑘.

Throughout this paper, we consider the following basic assumptions on functions 𝜙_𝑗(𝑥),  𝑥∈(0,+∞), 𝑗=1,2,…,𝑛: (𝐶₁)𝜙_𝑗(𝑥) is a 𝐶¹-function on (0,+∞), 𝜙_𝑗(𝑥)>0 for 𝑥>0 and 𝜙_𝑗(0)=0; (𝐶₂) there exist positive constants 𝑐_𝑗 such that 𝜙_𝑗(𝑥)≤𝑐_𝑗𝑥, 𝑥∈(0,+∞); (𝐶₃)lim_𝑥→0⁺(𝜙_𝑗(𝑥)/𝑥)=𝑐_𝑗; (𝐶₄)[𝜙_𝑗(𝑥)−𝜙_𝑗(𝑦)][𝜙_𝑗(𝑥)/𝑥−𝜙_𝑗(𝑦)/𝑦]≤0, for any 𝑥>0 and 𝑦>0.

It is easy to check that classes of 𝜙_𝑗(𝑥) satisfying (𝐶₁)--(𝐶₃) include common incidence functions such as 𝜙_𝑗(𝐼_𝑗)=𝐼_𝑗 [7], 𝜙_𝑗(𝐼_𝑗)=𝐼_𝑗/(1+𝛼_𝑗𝐼_𝑗) [14], 𝜙_𝑗(𝐼_𝑗)=𝐼_𝑗/(1+𝛼𝐼²_𝑗) [15]. Furthermore, 𝜙_𝑗(𝐼_𝑗)=𝐼_𝑗, 𝜙_𝑗(𝐼_𝑗)=𝐼_𝑗/(1+𝛼_𝑗𝐼_𝑗) satisfy (𝐶₄).

In Sun [13], the author considered a general multigroup SIR models with nonlinear incidence by using the same method as in [9]. Noting that model (1) is a special case of multigroup SIR models appeared in [13], according to the results of [13], the following results hold for system (1).

Proposition 1. Assume that 𝐵=(𝛽_𝑘𝑗) is irreducible and 𝜙_𝑗(𝑥) satisfies (𝐶₁)--(𝐶₃).(1)If ℛ₀≤1, then for system (1), 𝑃₀=(𝑆⁰₁,0,0,𝑆⁰₂,0,0,…,𝑆⁰_𝑛,0,0) is the unique equilibrium and it is globally asymptotically stable in Γ;(2)if ℛ₀>1 and (𝐶₄) is satisfied, then there exists an endemic equilibrium 𝑃^∗=(𝑆^∗₁,𝐼^∗₁,𝑅^∗₁,…,𝑆^∗_𝑛,𝐼^∗_𝑛,𝑅^∗_𝑛) for system (1) and 𝑃^∗ is globally asymptotically stable in Γ,
where 𝑆⁰_𝑘=Λ_𝑘/𝑑_𝑘,  𝑘=1,2,…,𝑛, ℛ₀=𝜌(𝑀₀) (the spectral radius of 𝑀₀), 𝑀₀ = 𝑀(𝑆⁰₁,𝑆⁰₂,…,𝑆⁰_𝑛) = (𝛽_𝑘𝑗𝑆⁰_𝑘𝑐_𝑗/(𝜖_𝑘+𝛾_𝑘))_𝑛×𝑛,  and Γ = {(𝑆₁,𝐼₁,𝑅₁,…,𝑆_𝑛,𝐼_𝑛,𝑅_𝑛)∈ℝ₊^3𝑛∶𝑆_𝑘≤Λ_𝑘/𝑑_𝑘,𝑆_𝑘+𝐼_𝑘≤Λ_𝑘/𝑑_𝑘,𝑘= 1,2,…,𝑛}.

The aim of this paper is to evaluate the effect of stochastic parameter perturbation on system (1). Considering the existence of environmental noise, we introduce randomness into the model (1) by replacing the parameters 𝑑_𝑘, 𝜖_𝑘, and 𝛿_𝑘 by 𝑑_𝑘⟶𝑑_𝑘+𝛼_𝑘𝐵̇_1𝑘(𝑡),𝜖_𝑘⟶𝜖_𝑘+𝛽_𝑘𝐵̇_2𝑘(𝑡),𝛿_𝑘⟶𝛿_𝑘+𝜎_𝑘𝐵̇_3𝑘(𝑡), (2) where 𝐵_1𝑘(𝑡), 𝐵_2𝑘(𝑡), 𝐵_3𝑘(𝑡), 𝑘=1,2,…,𝑛 are mutual independent standard Brownian motions with 𝐵_1𝑘(0)=0,  𝐵_2𝑘(0)=0,  𝐵_3𝑘(0)=0, and the intensities of white noises 𝛼²_𝑘>0, 𝛽²_𝑘>0,  𝜎²_𝑘>0, respectively. Then the stochastic version corresponding to the deterministic model (1) takes the following form: 𝑑𝑆_𝑘(𝑡)=Λ_𝑘−^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘(𝑡)𝜙_𝑗𝐼_𝑗(𝑡)−𝑑_𝑘𝑆_𝑘(𝑡)𝑑𝑡−𝛼_𝑘𝑆_𝑘(𝑡)𝑑𝐵_1𝑘(𝑡),𝑑𝐼_𝑘(𝑡)=^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘(𝑡)𝜙_𝑗𝐼_𝑗(𝑡)−𝜖_𝑘+𝛾_𝑘𝐼_𝑘(𝑡)𝑑𝑡−𝛽_𝑘𝐼_𝑘(𝑡)𝑑𝐵_2𝑘(𝑡),𝑑𝑅_𝑘(𝑡)=𝛾_𝑘𝐼_𝑘(𝑡)−𝛿_𝑘𝑅_𝑘(𝑡)𝑑𝑡−𝜎_𝑘𝑅_𝑘(𝑡)𝑑𝐵_3𝑘(𝑡),𝑘=1,2,…,𝑛. (3) System (3) with bilinear incidence has been researched by Ji et al. [16]. They proved that if ℛ₀≤1, the solution of the model is fluctuating around a steady state, which is the disease-free equilibrium of the corresponding deterministic model; while if ℛ₀>1, there is a stationary distribution, which means the disease will prevail.

In this paper, we will prove that when the intensity of white noise is small, the solution of system (3) converges weakly to a singular measure (i.e., a distribution) if ℛ₀≤1 and there exists an invariant distribution which is ergodic if ℛ₀>1. This is the same situation as in Proposition 1. When the intensity of white noise is large, the disease will extinct exponentially regardless of the magnitude of ℛ₀. Compared to the results of Ji et al. [16] our results can provide a deep insight into the dynamics of corresponding multigroup model.

It is worth mentioning that by now there are many excellent works about stochastic single-group SIR model. Beretta et al. [17] considered a stochastic SIR model with time delays and obtained asymptotic mean square stability conditions for positive equilibrium. In this paper, the authors assumed that stochastic perturbations are of white noise type, which are directly proportional to distances 𝑆(𝑡),  𝐼(𝑡),  𝑅(𝑡) from values of 𝑆^∗,  𝐼^∗,  𝑅^∗, influence on the 𝑆(̇𝑡),  𝐼(̇𝑡),  𝑅̇(𝑡), respectively. Tornatore et al. [18] discuss a single-group case of (3). They proved that 0<𝛽<min{𝜆+𝜇−𝜎²/2,2𝜇} is a sufficient condition for the asymptotic stability of the disease-free equilibrium. And only by computer simulations they showed that if min{𝜆+𝜇−𝜎²/2,2𝜇}<𝛽<𝜆+𝜇+𝜎²/2, the disease-free equilibrium is stable and the disease does not occur; if (𝜆+𝜇)(𝑅₀−1)>𝜎²/2, the disease-free equilibrium is unstable. Ji et al. [19] considered the same model as in [18]. They deduce the globally asymptotical stability and exponential mean square stability of the disease-free equilibrium under some conditions and investigate the asymptotic behavior of the solution around the endemic equilibrium of the deterministic model.

The rest of this paper is organized as follows. In Section 2, we show there is a unique nonnegative solution of system (3). In Section 3, if ℛ₀≤1 combined with small or large enough intensity of white noise, we show that the solution converges weakly to a singular measure. Section 4 focuses on the persistence of the disease. By choosing appropriate Lyapunov function, we show that there is a stationary distribution for system (3) and it is ergodic, if ℛ₀>1. Finally, for the self-contained, we present an Appendix which contains some results used in the previous sections.

Throughout this paper, let (Ω,ℱ,{ℱ_𝑡}_𝑡≥0,𝑃) be a complete probability space with a filtration {ℱ_𝑡}_𝑡≥0 satisfying the usual conditions (i.e., it is right continuous and ℱ₀ contains all 𝑃-null sets). Denote ℝ₊^3𝑛=𝑥∈ℝ^3𝑛∶𝑥_𝑘>0∀1≤𝑘≤3𝑛,ℝ₊^3𝑛=𝑥∈ℝ^3𝑛∶𝑥_𝑘≥0∀1≤𝑘≤3𝑛. (4) In general, consider a 𝑑-dimensional stochastic differential equation 𝑑𝑥(𝑡)=𝑓(𝑥(𝑡),𝑡)𝑑𝑡+𝑔(𝑥(𝑡),𝑡)𝑑𝐵(𝑡)on𝑡≥𝑡₀ (5) with initial value 𝑥(𝑡₀)=𝑥₀∈ℝ^𝑑. 𝐵(𝑡) denotes 𝑑-dimensional standard Brownian motion defined on the above probability space. Define the differential operator 𝐿 associated with (5) by 𝜕𝐿=+𝜕𝑡^𝑑_𝑘=1𝑓_𝑘()𝜕𝑥,𝑡𝜕𝑥_𝑘+12^𝑑_𝑘,𝑗=1𝑔^𝑇(𝑥,𝑡)𝑔(𝑥,𝑡)_𝑘𝑗𝜕²𝜕𝑥_𝑘𝜕𝑥_𝑗. (6) If 𝐿 acts on a function 𝑉∈𝐶^2,1(𝑅^𝑑×ℝ₊;ℝ₊), then 𝐿𝑉(𝑥,𝑡)=𝑉_𝑡(𝑥,𝑡)+𝑉_𝑥()𝑓()+1𝑥,𝑡𝑥,𝑡2trace𝑔^𝑇(𝑥,𝑡)𝑉_𝑥𝑥(𝑥,𝑡)𝑔(𝑥,𝑡), (7) where 𝑉_𝑡=𝜕𝑉/𝜕𝑡,  𝑉_𝑥=(𝜕𝑉/𝜕𝑥₁,…,𝜕𝑉/𝜕𝑥_𝑑), 𝑉_𝑥𝑥=(𝜕²𝑉/𝜕𝑥_𝑘𝜕𝑥_𝑗)_𝑑×𝑑. By Itô's formula, if 𝑥(𝑡)∈𝑅^𝑑, then 𝑑𝑉(𝑥(𝑡),𝑡)=𝐿𝑉(𝑥(𝑡),𝑡)𝑑𝑡+𝑉_𝑥(𝑥(𝑡),𝑡)𝑔(𝑥(𝑡),𝑡)𝑑𝐵(𝑡). (8)

Consider (5), assume 𝑓(0,𝑡)=0 and 𝑔(0,𝑡)=0 for all 𝑡≥𝑡₀. So 𝑥(𝑡)≡0 is a solution of (5), called the trivial solution or equilibrium position.

2. Existence and Uniqueness of the Nonnegative Solution

For a population model, one is interested in whether the solution is nonnegative and global. Hence in this section we show that the solution of system (3) is global and nonnegative. By making the change of variables and Lyapunov analysis method (see [20]), we will show global existence and uniqueness of the positive solution. From now on, we denote the solution (𝑆₁(𝑡),𝐼₁(𝑡),𝑅₁(𝑡),…,𝑆_𝑛(𝑡),𝐼_𝑛(𝑡),𝑅_𝑛(𝑡)) of system (3) as 𝑌(𝑡).

Theorem 2. If assumption (𝐶₁) holds and 𝐵=(𝛽_𝑘𝑗)_𝑛×𝑛 is irreducible, then for any initial value 𝑌(0)∈𝑅₊^3𝑛, there is a unique solution 𝑌(𝑡) of system (3) on 𝑡≥0, and the solution will remain in 𝑅₊^3𝑛 with probability 1.

Proof. First consider system 𝑑𝑢_𝑘(𝑡)=Λ_𝑘𝑒^{−𝑢𝑘(𝑡)}−^𝑛_𝑗=1𝛽_𝑘𝑗𝜙_𝑗𝑒^{𝑣𝑗(𝑡)}−𝑑_𝑘−𝛼²_𝑘2𝑑𝑡−𝛼_𝑘𝑑𝐵_1𝑘(𝑡),𝑑𝑣_𝑘(𝑡)=^𝑛_𝑗=1𝛽_𝑘𝑗𝑒^{𝑢𝑘(𝑡)−𝑣𝑘(𝑡)}𝜙_𝑗𝑒^{𝑣𝑗(𝑡)}−𝜖_𝑘+𝛾_𝑘−𝛽²_𝑘2𝑑𝑡−𝛽_𝑘𝑑𝐵_2𝑘(𝑡),𝑑𝑤_𝑘(𝑡)=𝛾_𝑘𝑒^{𝑣𝑘(𝑡)}−𝛿_𝑘𝑤_𝑘(𝑡)𝑑𝑡−𝜎_𝑘𝑤_𝑘(𝑡)𝑑𝐵_3𝑘(𝑡),𝑘=1,2,…,𝑛. (9) With initial value 𝑢_𝑘(0)=ln𝑆_𝑘(0), 𝑣_𝑘(0)=ln𝐼_𝑘(0), 𝑤_𝑘(0)=𝑅_𝑘(0), 1≤𝑘≤𝑛. Noting that (𝐶₁) holds, the coefficients of system (9) satisfy the local Lipschitz condition; then there is a unique local solution (𝑢₁(𝑡),𝑣₁(𝑡),𝑤₁(𝑡),…,𝑢_𝑛(𝑡),𝑣_𝑛(𝑡),𝑤_𝑛(𝑡)) on 𝑡∈[0,𝜏_𝑒), where 𝜏_𝑒 is the explosion time. Therefore, by Itô’s formula, it is easy to see (𝑒^𝑢1(𝑡),𝑒^𝑣1(𝑡),𝑤₁(𝑡),…,𝑒^{𝑢𝑛(𝑡)},𝑒^{𝑣𝑛(𝑡)},𝑤_𝑛(𝑡)) is the unique positive local solution to system (3) with initial value 𝑌(0)∈𝑅₊^3𝑛.
Next, we will prove that this solution is global. The following proof is almost the same as in the proof of Theorem 3.1 in [16]. We do not alter any words except replacing 𝐼_𝑗 by 𝜙_𝑗(𝐼_𝑗) and therefore we can obtain that 𝐿𝑉≤𝑀. Here we omit the details.

3. Exponential Stability of Infectious Disease

It is clear that 𝑃₀=(Λ₁/𝑑₁,0,0,Λ₂/𝑑₂,0,0,…,Λ_𝑛/𝑑_𝑛,0,0) is the disease-free equilibrium of system (1) but not (3). For system (1), which has been mentioned in Introduction, 𝑃₀ is globally stable if ℛ₀≤1, 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. Although there is none of disease-free equilibrium of system (3), in this section, we can still present sufficient conditions for the disease to extinct exponentially.

Theorem 3. Assume (𝐶₁),(𝐶₂) hold and 𝐵=(𝛽_𝑘𝑗)_𝑛×𝑛 is irreducible. If 2𝑑_𝑘>𝛼²_𝑘,  𝑘=1,2,…,𝑛, then the solution 𝑌(𝑡) of system (3) with initial value 𝑌(0)∈𝑅₊^3𝑛 has the property lim_𝑡→+∞(log(∑^𝑛_𝑘=1𝑎_𝑘𝐼_𝑘)/𝑡)≤𝑐,𝑎.𝑠. If 𝑐<0, the disease will extinct exponentially almost surely. Here 𝑀₀=𝑀𝑆⁰=𝛽_𝑘𝑗𝑆⁰_𝑘𝑐_𝑗𝜖_𝑘+𝛾_{𝑘𝑛×𝑛},𝑆⁰_𝑘=Λ_𝑘𝑑_𝑘,𝑎_𝑘=𝜔_𝑘𝜖_𝑘+𝛾_𝑘𝜔,𝑘=1,2,…,𝑛,₁,𝜔₂,…,𝜔_𝑛𝑖𝑠𝑎𝑙𝑒𝑓𝑡𝑒𝑖𝑔𝑒𝑛𝑣𝑒𝑐𝑡𝑜𝑟𝑜𝑓𝑀₀𝑐𝑜𝑟𝑟𝑒𝑠𝑝𝑜𝑛𝑑𝑖𝑛𝑔𝑡𝑜𝜌𝑀₀,𝑐∶=max_{0≤𝑘,𝑗≤𝑛}𝛽_𝑘𝑗𝑐_𝑗𝑎_𝑗^𝑛_𝑘=1𝑎_𝑘𝛼_𝑘𝑆⁰_𝑘2𝑑_𝑘−𝛼²_𝑘+min_{1≤𝑘≤𝑛}𝜖_𝑘+𝛾_𝑘ℛ₀−1𝟏_ℛ0≤1+max_{1≤𝑘≤𝑛}𝜖_𝑘+𝛾_𝑘ℛ₀−1𝟏_ℛ0≥1−12∑^𝑛_𝑘=11/𝛽²_𝑘. (10)

Proof. Let 𝑋_𝑘, 1≤𝑘≤𝑛 be the solution of the equation: 𝑑𝑋_𝑘(𝑡)=Λ_𝑘−𝑑_𝑘𝑋_𝑘(𝑡)𝑑𝑡−𝛼_𝑘𝑋_𝑘(𝑡)𝑑𝐵_1𝑘(𝑡),𝑋_𝑘(0)=𝑆_𝑘(0). (11) By comparison theorem for stochastic equations, we have 𝑆_𝑘(𝑡)≤𝑋_𝑘(𝑡),a.s. (12) It is easy to know that 𝑋_𝑘(𝑡) has the following property: lim_𝑡→∞∫^𝑡₀𝑋_𝑘𝑑𝑡𝑡=𝑆⁰_𝑘a.s.,lim_𝑡→∞∫^𝑡₀𝑋_𝑘−𝑆⁰_𝑘𝑑𝑡𝑡≤𝛼_𝑘𝑆⁰_𝑘2𝑑_𝑘−𝛼²_𝑘,if2𝑑_𝑘>𝛼²_𝑘,𝑘=1,2,…,𝑛. (13)(14) As for the proof of the above two properties, the reader may refer to the proof of Theorems 3.1 and 4.1 in [21] for details.
Since 𝐵=(𝛽_𝑘𝑗)_𝑛×𝑛 is irreducible, 𝛽_𝑘𝑗≥0,  𝑘,𝑗=1,2,…,𝑛 and 𝑆⁰_𝑘>0,  𝜖_𝑘+𝛾_𝑘>0,  𝑐_𝑘>0,   𝑘=1,2,…,𝑛, then 𝑀₀ is also nonnegative and irreducible. Hence by Lemma A.1, there is a left eigenvector of 𝑀₀ corresponding to 𝜌(𝑀₀), which is denoted as (𝜔₁,𝜔₂,…,𝜔_𝑛) and 𝜔_𝑘>0,  𝑘=1,2,…,𝑛; that is, 𝜔₁,𝜔₂,…,𝜔_𝑛𝜌𝑀₀=𝜔₁,𝜔₂,…,𝜔_𝑛𝑀₀. (15) Let 𝑎_𝑘=𝜔_𝑘/(𝜖_𝑘+𝛾_𝑘), 𝑘=1,2,…,𝑛. Define 𝐶²-function 𝑉∶𝑅^𝑛₊→𝑅₊ as follows: 𝑉𝐼₁,𝐼₂,…,𝐼_𝑛=^𝑛_𝑘=1𝑎_𝑘𝐼_𝑘. (16) By Itô’s formula, we compute 𝑑1log𝑉=𝐿log𝑉𝑑𝑡−𝑉^𝑛_𝑘=1𝑎_𝑘𝛽_𝑘𝐼_𝑘𝑑𝐵_2𝑘, (17) where 𝐿=1log𝑉𝑉^𝑛_𝑘=1𝑎_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝜙_𝑗𝐼_𝑗−𝜖_𝑘+𝛾_𝑘𝐼_𝑘−∑^𝑛_𝑘=1𝑎²_𝑘𝛽²_𝑘𝐼²_𝑘2𝑉²≤1𝑉^𝑛_𝑘=1𝑎_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑋_𝑘𝑐_𝑗𝐼_𝑗−𝜖_𝑘+𝛾_𝑘𝐼_𝑘−∑^𝑛_𝑘=1𝑎²_𝑘𝛽²_𝑘𝐼²_𝑘2𝑉²=1𝑉^𝑛_𝑘=1𝑎_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑐_𝑗𝑋_𝑘−𝑆⁰_𝑘𝐼_𝑗+^𝑛_𝑗=1𝛽_𝑘𝑗𝑐_𝑗𝑆⁰_𝑘𝐼_𝑗−𝜖_𝑘+𝛾_𝑘𝐼_𝑘−∑^𝑛_𝑘=1𝑎²_𝑘𝛽²_𝑘𝐼²_𝑘2𝑉² (18) according to (𝐶₂) and (12). Let 𝐻₁=1𝑉^𝑛_𝑘=1𝑎_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑐_𝑗𝑋_𝑘−𝑆⁰_𝑘𝐼_𝑗,𝐻₂=1𝑉^𝑛_𝑘=1𝑎_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑐_𝑗𝑆⁰_𝑘𝐼_𝑗−𝜖_𝑘+𝛾_𝑘𝐼_𝑘,𝐻₃∑=−^𝑛𝑘=1𝑎²_𝑘𝛽²_𝑘𝐼²_𝑘2𝑉². (19) Then 𝐻₁≤1𝑉^𝑛_𝑘=1𝑎_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑐_𝑗𝑋_𝑘−𝑆⁰_𝑘𝐼_𝑗=1𝑉^𝑛_𝑘=1𝑎_𝑘𝑋_𝑘−𝑆⁰_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑐_𝑗𝐼_𝑗≤max_{0≤𝑘,𝑗≤𝑛}𝛽_𝑘𝑗𝑐_𝑗𝑎_𝑗1𝑉^𝑛_𝑘=1𝑎_𝑘𝑋_𝑘−𝑆⁰_𝑘^𝑛_𝑗=1𝑎_𝑗𝐼_𝑗=max_{0≤𝑘,𝑗≤𝑛}𝛽_𝑘𝑗𝑐_𝑗𝑎_𝑗^𝑛_𝑘=1𝑎_𝑘𝑋_𝑘−𝑆⁰_𝑘. (20) In view of the definition of 𝑎_𝑘 and (15), we have 𝐻₂=1𝑉^𝑛𝑛_𝑘=1𝑗=1𝑎_𝑘𝛽_𝑘𝑗𝑐_𝑗𝑆⁰_𝑘𝐼_𝑗−^𝑛_𝑘=1𝑎_𝑘𝜖_𝑘+𝛾_𝑘𝐼_𝑘=1𝑉^𝑛𝑛_𝑘=1𝑗=1𝜔_𝑘𝛽_𝑘𝑗𝑐_𝑗𝑆⁰_𝑘𝜖_𝑘+𝛾_𝑘𝐼_𝑗−^𝑛_𝑘=1𝜔_𝑘𝐼_𝑘=1𝑉𝜔₁,𝜔₂,…,𝜔_𝑛×𝑀₀𝐼₁,𝐼₂,…,𝐼_𝑛^𝑇−𝐼₁,𝐼₂,…,𝐼_𝑛^𝑇=1𝑉ℛ₀−1^𝑛_𝑘=1𝜔_𝑘𝐼_𝑘=1𝑉ℛ₀−1^𝑛_𝑘=1𝜖_𝑘+𝛾_𝑘𝑎_𝑘𝐼_𝑘≤min_{1≤𝑘≤𝑛}𝜖_𝑘+𝛾_𝑘ℛ₀−1𝟏_ℛ0≤1+max_{1≤𝑘≤𝑛}𝜖_𝑘+𝛾_𝑘ℛ₀−1𝟏_ℛ0≥1∶=𝑚₁,𝐻₃∑=−^𝑛𝑘=1𝑎²_𝑘𝛽²_𝑘𝐼²_𝑘2𝑉²∑=−^𝑛𝑘=1𝑎²_𝑘𝛽²_𝑘𝐼²_𝑘2∑^𝑛_𝑘=1𝑎_𝑘/𝛽_𝑘𝛽_𝑘𝐼_𝑘²∑≤−^𝑛𝑘=1𝑎²_𝑘𝛽²_𝑘𝐼²_𝑘2∑^𝑛_𝑘=1𝑎²_𝑘𝛽²_𝑘𝐼²_𝑘∑^𝑛_𝑘=11/𝛽²_𝑘1=−2∑^𝑛_𝑘=11/𝛽²_𝑘∶=𝑚₂. (21) Hence, 𝐿log𝑉≤max_{0≤𝑘,𝑗≤𝑛}𝛽_𝑘𝑗𝑐_𝑗𝑎_𝑗×^𝑛_𝑘=1𝑎_𝑘𝑋_𝑘−𝑆⁰_𝑘+𝑚₁+𝑚₂, (22) which combined with (17) yields log𝑉(𝑡)𝑡≤log𝑉(0)𝑡+max_{0≤𝑘,𝑗≤𝑛}𝛽_𝑘𝑗𝑐_𝑗𝑎_𝑗^𝑛_𝑘=1𝑎_𝑘∫^𝑡₀𝑋_𝑘−𝑆⁰_𝑘𝑑𝑡𝑡+𝑚₁+𝑚₂−^𝑛_𝑘=1𝑎_𝑘𝛽_𝑘∫^𝑡₀𝐼_𝑘/𝑉𝑑𝐵_2𝑘𝑡. (23) By [22, Lemma 2.6], we get lim_𝑡→+∞∫^𝑡₀𝐼_𝑘/𝑉𝑑𝐵_2𝑘𝑡=0,a.s. (24) This together with (14) implies lim_𝑡→∞log𝑉𝑡≤max_{0≤𝑘,𝑗≤𝑛}𝛽_𝑘𝑗𝑐_𝑗𝑎_𝑗^𝑛_𝑘=1𝑎_𝑘𝛼_𝑘𝑆⁰_𝑘2𝑑_𝑘−𝛼²_𝑘+𝑚₁+𝑚₂,a.s. (25) Thus, the proof of Theorem 3 is completed.

Remark 4. From Theorem 3 we know that if ℛ₀<1 and 𝛼_𝑘, 1≤𝑘≤𝑛 are small, the disease will extinct exponentially; that is, 𝐼_𝑘→0 as 𝑡→∞. By using the same arguments as in [20, 21], we get 𝑆_𝑘(𝑡)−→^𝑤𝜈_𝑘 and 𝑅_𝑘(𝑡)→0 (𝑡→+∞), where −→^𝑤 means the weak convergence and 𝜈_𝑘 is a distribution in 𝑅₊ such that ∫^∞₀𝑥𝜈_𝑘(𝑑𝑥)=𝑆⁰_𝑘 and its density is (𝑀_𝑘𝑥²𝑝(𝑥))⁻¹, where 𝑀_𝑘 is a normal constant and 𝑝(𝑥)=𝑥^{2𝑑𝑘/𝛼2𝑘}exp(2𝑑_𝑘𝑆⁰_𝑘/𝛼²_𝑘𝑥), 𝑥>0. That is, the solution of system (3) converges weakly to a singular measure (i.e., a distribution) when ℛ₀<1 and 𝛼_𝑘, 1≤𝑘≤𝑛 are small.

Theorem 5. Assume (𝐶₁), (𝐶₂) hold and 𝐵=(𝛽_𝑘𝑗)_𝑛×𝑛 is irreducible. Then the solution 𝑌(𝑡) of system (3) with initial value 𝑌(0)∈𝑅₊^3𝑛 has the property lim_𝑡→+∞(log(∑^𝑛_𝑘=1𝐼_𝑘)/𝑡)≤𝑐,𝑎.𝑠. If 𝑐<0, the disease will extinct exponentially almost surely. Here 𝑐∶=max_{1≤𝑘,𝑗≤𝑛}{𝛽_𝑘𝑗𝑐_𝑗}∑^𝑛_𝑘=1𝑆⁰_𝑘 − min_{1≤𝑘≤𝑛}{𝜖_𝑘+𝛾_𝑘}−1/2∑^𝑛_𝑘=1(1/𝛽²_𝑘).

Proof. Define 𝐶²-function 𝑉∶𝑅^𝑛₊→𝑅₊ by 𝑉(𝐼₁,𝐼₂,…,𝐼_𝑛)=∑^𝑛_𝑘=1𝐼_𝑘. By Itô's formula, we compute 𝑑1log𝑉=𝐿log𝑉𝑑𝑡−𝑉^𝑛_𝑘=1𝛽_𝑘𝐼_𝑘𝑑𝐵_2𝑘, (26) where 𝐿=1log𝑉𝑉^𝑛_𝑘=1^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝜙_𝑗𝐼_𝑗−𝜖_𝑘+𝛾_𝑘𝐼_𝑘−∑^𝑛_𝑘=1𝛽²_𝑘𝐼²_𝑘2𝑉²≤1𝑉^𝑛_𝑘=1^𝑛_𝑗=1𝛽_𝑘𝑗𝑋_𝑘𝑐_𝑗𝐼_𝑗−𝜖_𝑘+𝛾_𝑘𝐼_𝑘−∑^𝑛_𝑘=1𝛽²_𝑘𝐼²_𝑘2∑^𝑛_𝑘=11/𝛽_𝑘𝛽_𝑘𝐼_𝑘²≤max_{1≤𝑘,𝑗≤𝑛}𝛽_𝑘𝑗𝑐_𝑗^𝑛_𝑘=1𝑋_𝑘−min_{1≤𝑘≤𝑛}𝜖_𝑘+𝛾_𝑘−12∑^𝑛_𝑘=11/𝛽²_𝑘. (27) Then log𝑉(𝑡)𝑡≤log𝑉(0)𝑡+max_{1≤𝑘,𝑗≤𝑛}𝛽_𝑘𝑗𝑐_𝑗^𝑛_𝑘=1∫^𝑡₀𝑋_𝑘𝑑𝑡𝑡−min_{1≤𝑘≤𝑛}𝜖_𝑘+𝛾_𝑘−12∑^𝑛_𝑘=11/𝛽²_𝑘−^𝑛_𝑘=1𝛽_𝑘∫^𝑡₀𝐼_𝑘/𝑉𝑑𝐵_2𝑘𝑡. (28) Taking 𝑡→∞, in view of (24) and (13), we get lim_𝑡→+∞log𝑉(𝑡)𝑡≤max_{1≤𝑘,𝑗≤𝑛}𝛽_𝑘𝑗𝑐_𝑗^𝑛_𝑘=1𝑆⁰_𝑘−min_{1≤𝑘≤𝑛}𝜖_𝑘+𝛾_𝑘−12∑^𝑛_𝑘=11/𝛽²_𝑘,a.s. (29)

Remark 6. Theorem 5 tells us that the large perturbation forces the infective to expire regardless of the magnitude of ℛ₀.

4. Ergodicity of System (3)

When studying epidemic dynamical systems, we are also interested in when the disease will prevail and persist in a population. For a deterministic model, this is usually solved by showing that the endemic equilibrium is a global attractor or is globally asymptotically stable. But, for the stochastic system (3), there is no endemic equilibrium. In this section, we explore a weak stability. We show that there is a stationary distribution based on the theory of [23] (see Appendix), which reveals the disease will prevail.

From the proof of Theorem 2, we obtain 𝐿𝑉≤𝑀. Let 𝑉=𝑉+𝑀. Then 𝐿𝑉≤𝑉 and it is clear that inf_𝑌∈𝑅+^3𝑛_⧵𝐷𝑘𝑉(𝑌)→∞ as 𝑘→∞, where 𝐷_𝑘=(1/𝑘,𝑘)×(1/𝑘,𝑘)×⋯×(1/𝑘,𝑘). Hence, by Remark 2 of Theorem 4.1 of Hasminskii ([23], p. 86), we obtain that the solution 𝑌(𝑡) is a homogeneous Markov process in 𝑅₊^3𝑛.

Theorem 7. Assume (𝐶₁)--(𝐶₄) hold, 𝐵=(𝛽_𝑘𝑗)_𝑛×𝑛 is irreducible, and ℛ₀=𝜌(𝑀₀)>1. If 0<𝛼²_𝑘<𝑑_𝑘,  0<𝛽²_𝑘<𝛾_𝑘+𝜖_𝑘,  0<𝜎²_𝑘<𝛿_𝑘,  𝑘=1,2,…,𝑛 such that 𝛿<min𝑐_𝑘2𝑑_𝑘−𝛼²_𝑘𝑆^∗_𝑘,𝑏_𝑘4𝜖_𝑘+𝛾_𝑘−𝛽²_𝑘𝐼^∗_𝑘²,𝑚_𝑘2𝛿_𝑘−𝜎²_𝑘𝑅^∗_𝑘²,𝑘=1,2,…,𝑛. (30) Then, for any initial value 𝑌(0)∈𝑅₊^3𝑛, there is a stationary distribution 𝜇(⋅) for system (3) and it has ergodic property, where 𝑀₀=𝑀(𝑆⁰)=(𝛽_𝑘𝑗𝑆⁰_𝑘𝑐_𝑗/(𝜖_𝑘+𝛾_𝑘))_𝑛×𝑛, 𝑆⁰_𝑘=Λ_𝑘/𝑑_𝑘, 𝛿=^𝑛_𝑘=1𝑎+22𝑐_𝑘+𝑏_𝑘𝑆^∗_𝑘𝑆^∗_𝑘𝛼²_𝑘+𝑎+12𝑐_𝑘+𝑏_𝑘𝐼^∗_𝑘𝐼^∗_𝑘𝛽²_𝑘+𝑚_𝑘𝑅^∗_𝑘²𝜎²_𝑘, (31)𝑎,  𝑏_𝑘,  𝑚_𝑘,  𝑘=1,2,…,𝑛 are defined as in the proof, 𝑃^∗=(𝑆^∗₁,𝐼^∗₁,𝑅^∗₁,…,𝑆^∗_𝑛,𝐼^∗_𝑛,𝑅^∗_𝑛) is the endemic equilibrium of system (1), and 𝑐_𝑘,  𝑘=1,2,…,𝑛 denote the cofactor of the 𝑘th diagonal element of 𝐿_𝐵, 𝐵=(𝛽_𝑘𝑗)_𝑛×𝑛=(𝛽_𝑘𝑗𝑆^∗_𝑘𝜙_𝑗(𝐼^∗_𝑗))_𝑛×𝑛.

Proof. Since ℛ₀>1 and (𝐶₁)--(𝐶₄) hold, from Proposition 1 there is an endemic equilibrium 𝑃^∗ of system (1). Then ^𝑛_𝑗=1𝛽_𝑘𝑗𝑆^∗_𝑘𝜙_𝑗𝐼^∗_𝑗+𝑑_𝑘𝑆^∗_𝑘=Λ_𝑘,_𝑛𝑗=1𝛽_𝑘𝑗𝑆^∗_𝑘𝜙_𝑗𝐼^∗_𝑗=𝜖_𝑘+𝛾_𝑘𝐼^∗_𝑘. (32) Define 𝑉(𝑌)=𝑎^𝑛_𝑘=1𝑐_𝑘𝑆_𝑘−𝑆^∗_𝑘−𝑆^∗_𝑘𝑆log^𝑘𝑆^∗_𝑘+𝐼_𝑘−𝐼^∗_𝑘−𝐼^∗_𝑘𝐼log^𝑘𝐼^∗_𝑘+^𝑛_𝑘=1𝑐_𝑘𝐼_𝑘−𝐼^∗_𝑘𝐼−log^𝑘𝐼^∗_𝑘+12^𝑛_𝑘=1𝑏_𝑘𝑆_𝑘−𝑆^∗_𝑘+𝐼_𝑘−𝐼^∗_𝑘²+12^𝑛_𝑘=1𝑙_𝑘𝑆_𝑘−𝑆^∗_𝑘²+12^𝑛_𝑘=1𝑚_𝑘𝑅_𝑘−𝑅^∗_𝑘²∶=𝑎𝑉₁+𝑉₂+𝑉₃+𝑉₄+𝑉₅, (33) where 𝑎,  𝑏_𝑘,  𝑚_𝑘,  𝑙_𝑘,  𝑘=1,2,…,𝑛 are positive constants to be determined later, and 𝑐_𝑘>0,  𝑘=1,2,…,𝑛 according to Lemma A.1. Then 𝑉 is positive definite. In view of (32) and the inequality (𝑥+𝑦)²≤2(𝑥²+𝑦²), by direct calculation, we get 𝐿𝑉_𝑖  (𝑖=1,2,…,5) as follows: 𝐿𝑉₁=^𝑛_𝑘=1𝑐_𝑘𝑆1−^∗𝑘𝑆_𝑘Λ_𝑘−^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝜙_𝑗𝐼_𝑗−𝑑_𝑘𝑆_𝑘+^𝑛_𝑘=1𝑐_𝑘𝑆^∗_𝑘𝛼²_𝑘2+^𝑛_𝑘=1𝑐_𝑘𝐼1−^∗𝑘𝐼_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝜙_𝑗𝐼_𝑗−𝜖_𝑘+𝛾_𝑘𝐼_𝑘+^𝑛_𝑘=1𝑐_𝑘𝐼^∗_𝑘𝛽²_𝑘2=^𝑛_𝑘=1𝑐_𝑘𝑆1−^∗𝑘𝑆_𝑘×^𝑛_𝑗=1𝛽_𝑘𝑗𝑆1−^𝑘𝜙_𝑗𝐼_𝑗𝑆^∗_𝑘𝜙_𝑗𝐼^∗_𝑗−𝑑_𝑘𝑆_𝑘−𝑆^∗_𝑘+^𝑛_𝑘=1𝑐_𝑘𝑆^∗_𝑘𝛼²_𝑘2+^𝑛_𝑘=1𝑐_𝑘𝐼1−^∗𝑘𝐼_𝑘×^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝜙_𝑗𝐼_𝑗𝑆^∗_𝑘𝜙_𝑗𝐼^∗_𝑗−^𝑛_𝑗=1𝛽_𝑘𝑗𝐼_𝑘𝐼^∗_𝑘+^𝑛_𝑘=1𝑐_𝑘𝐼^∗_𝑘𝛽²_𝑘2=^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆1−^∗𝑘𝑆_𝑘+𝜙_𝑗𝐼_𝑗𝜙_𝑗𝐼^∗_𝑗−𝑆_𝑘𝜙_𝑗𝐼_𝑗𝑆^∗_𝑘𝜙_𝑗𝐼^∗_𝑗−𝑑_𝑘𝑆_𝑘−𝑆^∗_𝑘²𝑆_𝑘+^𝑛_𝑘=1𝑐_𝑘𝑆^∗_𝑘𝛼²_𝑘2+^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝐼1−^𝑘𝐼^∗_𝑘+𝑆_𝑘𝜙_𝑗𝐼_𝑗𝑆^∗_𝑘𝜙_𝑗𝐼^∗_𝑗−𝑆_𝑘𝜙_𝑗𝐼_𝑗𝐼^∗_𝑘𝑆^∗_𝑘𝜙_𝑗𝐼^∗_𝑗𝐼_𝑘+^𝑛_𝑘=1𝑐_𝑘𝐼^∗_𝑘𝛽²_𝑘2=^𝑛_𝑘=1𝑐_𝑘−𝑑_𝑘𝑆_𝑘−𝑆^∗_𝑘²𝑆_𝑘+^𝑛_𝑗=1𝛽_𝑘𝑗𝑆2−^∗𝑘𝑆_𝑘−𝐼_𝑘𝐼^∗_𝑘+𝜙_𝑗𝐼_𝑗𝜙_𝑗𝐼^∗_𝑗−𝑆_𝑘𝜙_𝑗𝐼_𝑗𝐼^∗_𝑘𝑆^∗_𝑘𝜙_𝑗𝐼^∗_𝑗𝐼_𝑘+𝑆^∗_𝑘𝛼²_𝑘2+𝐼^∗_𝑘𝛽²_𝑘2. (34) From the calculation of 𝐿𝑉₁, we directly get 𝐿𝑉₂=^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝐼1−^𝑘𝐼^∗_𝑘+𝑆_𝑘𝜙_𝑗𝐼_𝑗𝑆^∗_𝑘𝜙_𝑗𝐼^∗_𝑗−𝑆_𝑘𝜙_𝑗𝐼_𝑗𝐼^∗_𝑘𝑆^∗_𝑘𝜙_𝑗𝐼^∗_𝑗𝐼_𝑘+^𝑛_𝑘=1𝑐_𝑘𝐼^∗_𝑘𝛽²_𝑘2𝐿𝑉₃=^𝑛_𝑘=1𝑏_𝑘𝑆_𝑘−𝑆^∗_𝑘+𝐼_𝑘−𝐼^∗_𝑘×Λ_𝑘−𝑑_𝑘𝑆_𝑘−𝜖_𝑘+𝛾_𝑘𝐼_𝑘+𝑆²_𝑘𝛼²_𝑘2+𝐼²_𝑘𝛽²_𝑘2=^𝑛_𝑘=1𝑏_𝑘𝑆_𝑘−𝑆^∗_𝑘+𝐼_𝑘−𝐼^∗_𝑘×−𝑑_𝑘𝑆_𝑘−𝑆^∗_𝑘−𝜖_𝑘+𝛾_𝑘𝐼_𝑘−𝐼^∗_𝑘+𝑆²_𝑘𝛼²_𝑘2+𝐼²_𝑘𝛽²_𝑘2≤−^𝑛_𝑘=1𝑏_𝑘𝑑_𝑘𝑆_𝑘−𝑆^∗_𝑘²+𝜖_𝑘+𝛾_𝑘𝐼_𝑘−𝐼^∗_𝑘²−𝛾_𝑘+𝑑_𝑘+𝜖_𝑘𝑆_𝑘−𝑆^∗_𝑘𝐼_𝑘−𝐼^∗_𝑘+^𝑛_𝑘=1𝑏_𝑘𝛼²_𝑘𝑆_𝑘−𝑆^∗_𝑘²+𝑆^∗_𝑘²+^𝑛_𝑘=1𝑏_𝑘𝛽²_𝑘𝐼_𝑘−𝐼^∗_𝑘²+𝐼^∗_𝑘²=−^𝑛_𝑘=1𝑏_𝑘𝑑_𝑘−𝛼²_𝑘𝑆_𝑘−𝑆^∗_𝑘²+𝜖_𝑘+𝛾_𝑘−𝛽²_𝑘𝐼_𝑘−𝐼^∗_𝑘²−^𝑛_𝑘=1𝑏_𝑘𝛾_𝑘+𝑑_𝑘+𝜖_𝑘𝑆_𝑘−𝑆^∗_𝑘𝐼_𝑘−𝐼^∗_𝑘+^𝑛_𝑘=1𝑏_𝑘𝛼²_𝑘𝑆^∗_𝑘²+𝛽²_𝑘𝐼^∗_𝑘²≤−^𝑛_𝑘=1𝑏_𝑘𝑑_𝑘−𝛼²_𝑘−𝑑_𝑘+𝛾_𝑘+𝜖_𝑘²2𝛾_𝑘+𝑑_𝑘+𝜖_𝑘−𝛽²_𝑘×𝑆_𝑘−𝑆^∗_𝑘²+12𝜖_𝑘+𝛾_𝑘−𝛽²_𝑘𝐼_𝑘−𝐼^∗_𝑘²+^𝑛_𝑘=1𝑏_𝑘𝛼²_𝑘𝑆^∗_𝑘²+𝛽²_𝑘𝐼^∗_𝑘²,𝐿𝑉₄=^𝑛_𝑘=1𝑙_𝑘𝑆_𝑘−𝑆^∗_𝑘×Λ_𝑘−^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝜙_𝑗𝐼_𝑗−𝑑_𝑘𝑆_𝑘+^𝑛_𝑘=1𝑙_𝑘𝑆²_𝑘𝛼²_𝑘2=^𝑛_𝑘=1𝑙_𝑘𝑆_𝑘−𝑆^∗_𝑘×^𝑛_𝑗=1𝛽_𝑘𝑗𝑆^∗_𝑘𝜙_𝑗𝐼^∗_𝑗−𝑆_𝑘𝜙_𝑗𝐼_𝑗−𝑑_𝑘𝑆_𝑘−𝑆^∗_𝑘+^𝑛_𝑘=1𝑙_𝑘𝑆²_𝑘𝛼²_𝑘2=−^𝑛𝑛_𝑘=1𝑗=1𝑙_𝑘𝛽_𝑘𝑗𝑆^∗_𝑘×𝑆_𝑘−𝑆^∗_𝑘𝜙_𝑗𝐼_𝑗−𝜙_𝑗𝐼^∗_𝑗−^𝑛𝑛_𝑘=1𝑗=1𝑙_𝑘𝛽_𝑘𝑗𝑆_𝑘−𝑆^∗_𝑘²𝜙_𝑗𝐼_𝑗−^𝑛_𝑘=1𝑙_𝑘𝑑_𝑘𝑆_𝑘−𝑆^∗_𝑘²+^𝑛_𝑘=1𝑙_𝑘𝑆²_𝑘𝛼²_𝑘2≤−^𝑛𝑛_𝑘=1𝑗=1𝑙_𝑘𝛽_𝑘𝑗𝑆^∗_𝑘𝑆_𝑘−𝑆^∗_𝑘𝜙_𝑗𝐼_𝑗−𝜙_𝑗𝐼^∗_𝑗−^𝑛_𝑘=1𝑙_𝑘𝑑_𝑘−𝛼²_𝑘𝑆_𝑘−𝑆^∗_𝑘²+^𝑛_𝑘=1𝑙_𝑘𝑆^∗_𝑘²𝛼²_𝑘,𝐿𝑉₅=^𝑛_𝑘=1𝑚_𝑘𝑅_𝑘−𝑅^∗_𝑘𝛾_𝑘𝐼_𝑘−𝛿_𝑘𝑅_𝑘+^𝑛_𝑘=1𝑚_𝑘𝑅²_𝑘𝜎²_𝑘2=^𝑛_𝑘=1𝑚_𝑘𝑅_𝑘−𝑅^∗_𝑘×𝛾_𝑘𝐼_𝑘−𝐼^∗_𝑘−𝛿_𝑘𝑅_𝑘−𝑅^∗_𝑘+^𝑛_𝑘=1𝑚_𝑘𝑅²_𝑘𝜎²_𝑘2≤^𝑛_𝑘=1𝑚_𝑘𝛾_𝑘𝐼_𝑘−𝐼^∗_𝑘𝑅_𝑘−𝑅^∗_𝑘−^𝑛_𝑘=1𝑚_𝑘𝛿_𝑘−𝜎²_𝑘𝑅_𝑘−𝑅^∗_𝑘²+^𝑛_𝑘=1𝑚_𝑘𝑅^∗_𝑘²𝜎²_𝑘≤^𝑛_𝑘=1𝑚_𝑘𝛾²_𝑘2𝛿_𝑘−𝜎²_𝑘𝐼_𝑘−𝐼^∗_𝑘²−12^𝑛_𝑘=1𝑚_𝑘𝛿_𝑘−𝜎²_𝑘𝑅_𝑘−𝑅^∗_𝑘²+^𝑛_𝑘=1𝑚_𝑘𝑅^∗_𝑘²𝜎²_𝑘. (35)(36) By Lemma A.3, we know ^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝜙_𝑗𝐼_𝑗𝜙_𝑗𝐼^∗_𝑗−^𝑛_𝑗=1𝛽_𝑘𝑗𝜙_𝑘𝐼_𝑘𝜙_𝑘𝐼^∗_𝑘=0,_𝑛𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝜙log^𝑗𝐼_𝑗𝜙_𝑗𝐼^∗_𝑗−^𝑛_𝑗=1𝛽_𝑘𝑗𝜙log^𝑘𝐼_𝑘𝜙_𝑘𝐼^∗_𝑘=0. (37)(38) Besides note that 𝑎−1−log𝑎≥0 for 𝑎>0; then 𝑆1−^∗𝑘𝑆_𝑘𝑆≤−log^∗𝑘𝑆_𝑘,𝑆1−^𝑘𝜙_𝑗𝐼_𝑗𝐼^∗_𝑘𝑆^∗_𝑘𝜙_𝑗𝐼^∗_𝑗𝐼_𝑘𝑆≤−log^𝑘𝜙_𝑗𝐼_𝑗𝐼^∗_𝑘𝑆^∗_𝑘𝜙_𝑗𝐼^∗_𝑗𝐼_𝑘𝑆=−log^𝑘𝑆^∗_𝑘𝜙−log^𝑗𝐼_𝑗𝜙_𝑗𝐼^∗_𝑗𝐼−log^∗𝑘𝐼_𝑘,𝜙−log^𝑘𝐼_𝑘𝜙_𝑘𝐼^∗_𝑘𝐼−log^∗𝑘𝐼_𝑘𝜙=−log^𝑘𝐼_𝑘𝐼^∗_𝑘𝜙_𝑘𝐼^∗_𝑘𝐼_𝑘≤𝜙_𝑘𝐼^∗_𝑘𝐼_𝑘𝜙_𝑘𝐼_𝑘𝐼^∗_𝑘−1. (39)(40)(41) Successively substituting (39), (40), (37), (38), and (41) into (34) yields 𝐿𝑉₁≤−^𝑛_𝑘=1𝑐_𝑘𝑑_𝑘𝑆_𝑘−𝑆^∗_𝑘²𝑆_𝑘+^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆−log^∗𝑘𝑆_𝑘𝑆−log^𝑘𝑆^∗_𝑘𝜙−log^𝑗𝐼_𝑗𝜙_𝑗𝐼^∗_𝑗𝐼−log^∗𝑘𝐼_𝑘−𝐼_𝑘𝐼^∗_𝑘+𝜙_𝑗𝐼_𝑗𝜙_𝑗𝐼^∗_𝑗+^𝑛_𝑘=1𝑐_𝑘𝑆^∗_𝑘𝛼²_𝑘2+𝐼^∗_𝑘𝛽²_𝑘2=−^𝑛_𝑘=1𝑐_𝑘𝑑_𝑘𝑆_𝑘−𝑆^∗_𝑘²𝑆_𝑘+^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝜙−log^𝑘𝐼_𝑘𝜙_𝑘𝐼^∗_𝑘𝐼−log^∗𝑘𝐼_𝑘−𝐼_𝑘𝐼^∗_𝑘+𝜙_𝑘𝐼_𝑘𝜙_𝑘𝐼^∗_𝑘+^𝑛_𝑘=1𝑐_𝑘𝑆^∗_𝑘𝛼²_𝑘2+𝐼^∗_𝑘𝛽²_𝑘2≤−^𝑛_𝑘=1𝑐_𝑘𝑑_𝑘𝑆_𝑘−𝑆^∗_𝑘²𝑆_𝑘+^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝜙_𝑘𝐼^∗_𝑘𝐼_𝑘𝜙_𝑘𝐼_𝑘𝐼^∗_𝑘𝐼−1−^𝑘𝐼^∗_𝑘+𝜙_𝑘𝐼_𝑘𝜙_𝑘𝐼^∗_𝑘+^𝑛_𝑘=1𝑐_𝑘𝑆^∗_𝑘𝛼²_𝑘2+𝐼^∗_𝑘𝛽²_𝑘2. (42) In view of (𝐶₄), we obtain 𝜙_𝑘𝐼^∗_𝑘𝐼_𝑘𝜙_𝑘𝐼_𝑘𝐼^∗_𝑘𝐼−1−^𝑘𝐼^∗_𝑘+𝜙_𝑘𝐼_𝑘𝜙_𝑘𝐼^∗_𝑘=1𝜙_𝑘𝐼^∗_𝑘𝐼_𝑘𝜙_𝑘𝐼_𝑘𝜙_𝑘𝐼^∗_𝑘−𝜙_𝑘𝐼_𝑘×𝜙_𝑘𝐼^∗_𝑘𝐼^∗_𝑘−𝜙_𝑘𝐼_𝑘𝐼_𝑘≤0. (43) Then 𝐿𝑉₁≤−^𝑛_𝑘=1𝑐_𝑘𝑑_𝑘𝑆_𝑘−𝑆^∗_𝑘²𝑆_𝑘+^𝑛_𝑘=1𝑐_𝑘𝑆^∗_𝑘𝛼²_𝑘2+𝐼^∗_𝑘𝛽²_𝑘2. (44) By using inequality 𝑎−1−log𝑎≥0 for 𝑎>0 again, it follows that 𝑆−log^𝑘𝑆^∗_𝑘≤𝑆^∗_𝑘𝑆_𝑘−1. (45) Successively substituting (40), (37), (45), (41), and (43) into (35), we obtain 𝐿𝑉₂≤^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗−𝐼_𝑘𝐼^∗_𝑘+𝑆_𝑘𝜙_𝑗𝐼_𝑗𝑆^∗_𝑘𝜙_𝑗𝐼^∗_𝑗𝑆−log^𝑘𝑆^∗_𝑘𝜙−log^𝑗𝐼_𝑗𝜙_𝑗𝐼^∗_𝑗𝐼−log^∗𝑘𝐼_𝑘+^𝑛_𝑘=1𝑐_𝑘𝐼^∗_𝑘𝛽²_𝑘2=^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗−𝐼_𝑘𝐼^∗_𝑘+𝑆_𝑘𝜙_𝑗𝐼_𝑗𝑆^∗_𝑘𝜙_𝑗𝐼^∗_𝑗𝑆−log^𝑘𝑆^∗_𝑘𝜙−log^𝑘𝐼_𝑘𝜙_𝑘𝐼^∗_𝑘𝐼−log^∗𝑘𝐼_𝑘+^𝑛_𝑘=1𝑐_𝑘𝐼^∗_𝑘𝛽²_𝑘2≤^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗−𝐼_𝑘𝐼^∗_𝑘+𝑆_𝑘𝜙_𝑗𝐼_𝑗𝑆^∗_𝑘𝜙_𝑗𝐼^∗_𝑗+𝑆^∗_𝑘𝑆_𝑘𝐼−1+^𝑘𝜙_𝑘𝐼^∗_𝑘𝐼^∗_𝑘𝜙_𝑘𝐼_𝑘+−1_𝑛𝑘=1𝑐_𝑘𝐼^∗_𝑘𝛽²_𝑘2=^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝑆^∗_𝑘𝜙−1^𝑗𝐼_𝑗𝜙_𝑗𝐼^∗_𝑗+𝑆−1_𝑘𝑆^∗_𝑘+𝑆^∗_𝑘𝑆_𝑘+𝜙−2_𝑘𝐼^∗_𝑘𝐼_𝑘𝜙_𝑘𝐼_𝑘𝐼^∗_𝑘𝐼−1−^𝑘𝐼^∗_𝑘+𝜙_𝑘𝐼_𝑘𝜙_𝑘𝐼^∗_𝑘+^𝑛_𝑘=1𝑐_𝑘𝐼^∗_𝑘𝛽²_𝑘2≤^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘𝑆^∗_𝑘𝜙−1^𝑗𝐼_𝑗𝜙_𝑗𝐼^∗_𝑗+𝑆−1_𝑘𝑆^∗_𝑘+𝑆^∗_𝑘𝑆_𝑘+−2_𝑛𝑘=1𝑐_𝑘𝐼^∗_𝑘𝛽²_𝑘2=^𝑛_𝑘=1𝑐_𝑘^𝑛_𝑗=1𝛽_𝑘𝑗𝑆_𝑘−𝑆^∗_𝑘𝜙_𝑗𝐼_𝑗−𝜙_𝑗𝐼^∗_𝑗+^𝑛_𝑗=1𝛽_𝑘𝑗𝜙_𝑗𝐼^∗_𝑗𝑆_𝑘−𝑆^∗_𝑘²𝑆_𝑘+^𝑛_𝑘=1𝑐_𝑘𝐼^∗_𝑘𝛽²_𝑘2. (46) Hence 𝑎𝐿𝑉₁+𝐿𝑉₂+𝐿𝑉₄≤−^𝑛_𝑘=1^𝑛_𝑗=1𝛽_𝑘𝑗𝑙_𝑘𝑆^∗_𝑘−𝑐_𝑘𝑆_𝑘−𝑆^∗_𝑘𝜙_𝑗𝐼_𝑗−𝜙_𝑗𝐼^∗_𝑗+𝑐_𝑘𝑎𝑑_𝑘−^𝑛_𝑗=1𝛽_𝑘𝑗𝜙_𝑗𝐼^∗_𝑗𝑆_𝑘−𝑆^∗_𝑘²𝑆_𝑘−^𝑛_𝑘=1𝑙_𝑘𝑑_𝑘−𝛼²_𝑘𝑆_𝑘−𝑆^∗_𝑘²+^𝑛_𝑘=1𝑐_𝑘𝛽²_𝑘2(𝑎+1)𝐼^∗_𝑘+𝛼²_𝑘2(𝑎+2)𝑆^∗_𝑘. (47) Choose 𝑙_𝑘=𝑐_𝑘/𝑆^∗_𝑘 and 𝑎=max{∑^𝑛_𝑗=1𝛽_𝑘𝑗𝐼^∗_𝑗/𝑑_𝑘,𝑘=1,2,…,𝑛} such that 𝑙_𝑘𝑆^∗_𝑘−𝑐_𝑘=0 and 𝑎𝑑_𝑘 − ∑^𝑛_𝑗=1𝛽_𝑘𝑗𝐼^∗_𝑗 ≥ 0,  𝑘=1,2,…,𝑛. Then 𝑎𝐿𝑉₁+𝐿𝑉₂+𝐿𝑉₄≤−^𝑛_𝑘=1𝑐_𝑘𝑆^∗_𝑘𝑑_𝑘−𝛼²_𝑘𝑆_𝑘−𝑆^∗_𝑘²+^𝑛_𝑘=1𝑐_𝑘𝛽²_𝑘2(𝑎+1)𝐼^∗_𝑘+𝛼²_𝑘2(𝑎+2)𝑆^∗_𝑘. (48) Furthermore, 𝑎𝐿𝑉₁+𝐿𝑉₂+𝐿𝑉₃+𝐿𝑉₄+𝐿𝑉₅≤−^𝑛_𝑘=1𝑏_𝑘𝑑_𝑘−𝛼²_𝑘−𝑑_𝑘+𝛾_𝑘+𝜖_𝑘²2𝛾_𝑘+𝑑_𝑘+𝜖_𝑘−𝛽²_𝑘+𝑐_𝑘𝑆^∗_𝑘𝑑_𝑘−𝛼²_𝑘𝑆_𝑘−𝑆^∗_𝑘²−12^𝑛_𝑘=1𝑏_𝑘𝜖_𝑘+𝛾_𝑘−𝛽²_𝑘−𝑚_𝑘𝛾²_𝑘𝛿_𝑘−𝜎²_𝑘𝐼_𝑘−𝐼^∗_𝑘²−12^𝑛_𝑘=1𝑚_𝑘𝛿_𝑘−𝜎²_𝑘𝑅_𝑘−𝑅^∗_𝑘²+^𝑛_𝑘=1𝑎+22𝑐_𝑘+𝑏_𝑘𝑆^∗_𝑘𝑆^∗_𝑘𝛼²_𝑘+𝑎+12𝑐_𝑘+𝑏_𝑘𝐼^∗_𝑘×𝐼^∗_𝑘𝛽²_𝑘+𝑚_𝑘𝑅^∗_𝑘²𝜎²_𝑘. (49) Choose 𝑏_𝑘>0 such that 𝑏_𝑘(𝑑_𝑘−𝛼²_𝑘−(𝑑_𝑘+𝛾_𝑘+𝜖_𝑘)²/2(𝛾_𝑘+𝑑_𝑘+𝜖_𝑘−𝛽²_𝑘)) + (𝑐_𝑘/𝑆^∗_𝑘)(𝑑_𝑘−𝛼²_𝑘) ≥ (𝑐_𝑘/2𝑆^∗_𝑘)(𝑑_𝑘−𝛼²_𝑘) and choose 𝑚_𝑘 = (𝑏_𝑘/2)(𝜖_𝑘+𝛾_𝑘−𝛽²_𝑘)((𝛿_𝑘−𝜎²_𝑘)/𝛾²_𝑘) such that 𝑏_𝑘(𝜖_𝑘+𝛾_𝑘−𝛽²_𝑘) − 𝑚_𝑘𝛾²_𝑘/(𝛿_𝑘−𝜎²_𝑘) = (𝑏_𝑘/2)(𝜖_𝑘+𝛾_𝑘−𝛽²_𝑘); then 𝑎𝐿𝑉₁+𝐿𝑉₂+𝐿𝑉₃+𝐿𝑉₄+𝐿𝑉₅≤−^𝑛_𝑘=1𝑐_𝑘2𝑆^∗_𝑘𝑑_𝑘−𝛼²_𝑘𝑆_𝑘−𝑆^∗_𝑘²−14^𝑛_𝑘=1𝑏_𝑘𝜖_𝑘+𝛾_𝑘−𝛽²_𝑘𝐼_𝑘−𝐼^∗_𝑘²−12^𝑛_𝑘=1𝑚_𝑘𝛿_𝑘−𝜎²_𝑘𝑅_𝑘−𝑅^∗_𝑘²+^𝑛_𝑘=1𝑎+22𝑐_𝑘+𝑏_𝑘𝑆^∗_𝑘𝑆^∗_𝑘𝛼²_𝑘+𝑎+12𝑐_𝑘+𝑏_𝑘𝐼^∗_𝑘𝐼^∗_𝑘𝛽²_𝑘+𝑚_𝑘𝑅^∗_𝑘²𝜎²_𝑘∶=−^𝑛_𝑘=1𝑐_𝑘2𝑆^∗_𝑘𝑑_𝑘−𝛼²_𝑘𝑆_𝑘−𝑆^∗_𝑘²−^𝑛_𝑘=1𝑏_𝑘4𝜖_𝑘+𝛾_𝑘−𝛽²_𝑘𝐼_𝑘−𝐼^∗_𝑘²−^𝑛_𝑘=1𝑚_𝑘2𝛿_𝑘−𝜎²_𝑘𝑅_𝑘−𝑅^∗_𝑘²+𝛿. (50) Note that 𝛿<min𝑐_𝑘2𝑑_𝑘−𝛼²_𝑘𝑆^∗_𝑘,𝑏_𝑘4𝜖_𝑘+𝛾_𝑘−𝛽²_𝑘𝐼^∗_𝑘²,𝑚_𝑘2𝛿_𝑘−𝜎²_𝑘𝑅^∗_𝑘²,𝑘=1,2,…,𝑛. (51) Then the ellipsoid −^𝑛_𝑘=1𝑐_𝑘2𝑆^∗_𝑘𝑑_𝑘−𝛼²_𝑘𝑆_𝑘−𝑆^∗_𝑘²−^𝑛_𝑘=1𝑏_𝑘4𝜖_𝑘+𝛾_𝑘−𝛽²_𝑘𝐼_𝑘−𝐼^∗_𝑘²−^𝑛_𝑘=1𝑚_𝑘2𝛿_𝑘−𝜎²_𝑘𝑅_𝑘−𝑅^∗_𝑘²+𝛿=0 (52) lies entirely in 𝑅₊^3𝑛. We can take 𝑈 to be a neighborhood of the ellipsoid with 𝑈⊆𝐸_3𝑛=𝑅₊^3𝑛, so for 𝑥∈𝐸_3𝑛⧵𝑈,  𝐿𝑉≤−𝐾(𝐾 is a positive constant), which implies condition (B.2) in Lemma A.5 is satisfied.
Note that the diffusion matrix of system (3) is 𝐴=𝑎_{𝑖𝑗𝑛×𝑛}=diag𝛼²₁𝑆²₁,𝛽²₁𝐼²₁,𝜎²₁𝑅²₁,…,𝛼²_𝑛𝑆²_𝑛,𝛽²_𝑛𝐼²_𝑛,𝜎²_𝑛𝑅²_𝑛. (53) For any bounded domain 𝐷⊂𝑅₊^3𝑛 there is 𝑀=min𝛼²_𝑖𝑆²_𝑖,𝛽²_𝑖𝐼²_𝑖,𝜎²_𝑖𝑅²_𝑖,𝑖=1,2,…,𝑛,𝑌∈𝐷>0 (54) such that ^3𝑛_𝑖,𝑗=1𝑎_𝑖𝑗𝜉_𝑖𝜉_𝑗=^𝑛_𝑖=1𝛼²_𝑖𝑆²_𝑖𝜉²_3𝑖−2+^𝑛_𝑖=1𝛽²_𝑖𝐼²_𝑖𝜉²_3𝑖−1+^𝑛_𝑖=1𝜎²_𝑖𝑅²_𝑖𝜉²_3𝑖≥𝑀𝜉² (55) for all 𝑌∈𝐷,  𝜉∈𝑅^3𝑛, which implies condition (B.1) is also satisfied. Therefore, by Lemma A.5, the stochastic system (3) has a stationary distribution 𝜇(⋅) and it is ergodic.

Appendix

For the completeness of the paper, in this section, we list some theory used in the upper sections.

(1) Some graph theory (see [9, 24]): a directed graph or a digraph 𝒢=(𝑉,𝐸) contains a set 𝑉={1,2,…,𝑛} of vertices and a set 𝐸 of arcs (𝑘,𝑗) leading from initial vertex 𝑘 to terminal vertex 𝑗. A digraph 𝒢 is weighted if each arc (𝑗,𝑘) is assigned a positive weight 𝑎_𝑘𝑗. Given a weighted digraph 𝒢 with 𝑛 vertices, define the weight matrix 𝐴=(𝑎_𝑘𝑗)_𝑛×𝑛 whose entry 𝑎_𝑘𝑗 equals the weight of arc (𝑗,𝑘) if it exists, and 0 otherwise. A weighted digraph is denoted as (𝒢,𝐴). A digraph 𝒢 is strongly connected if for any pair of distinct vertices, there exists a directed path from one to the other. A weighted digraph (𝒢,𝐴) is strongly connected if and only if the weight matrix 𝐴 is irreducible [25].

Lemma A.1. If 𝐴 is nonnegative and irreducible, then the spectral radius 𝜌(𝐴) of 𝐴 is a simple eigenvalue, and 𝐴 has a positive eigenvector 𝜔=(𝜔₁,𝜔₂,…,𝜔_𝑛) corresponding to 𝜌(𝐴).

Remark A.2. The result of this lemma is given in P. 27, [25, Theorem 1.4].

The Laplacian matrix of (𝒢,𝐴) is defined as 𝐿_𝐴=_𝑘∑_≠1𝑎_1𝑘−𝑎₁₂⋯−𝑎_1𝑛−𝑎₂₁∑_𝑘≠2𝑎_2𝑘⋯−𝑎_2𝑛⋮⋮⋱⋮−𝑎_𝑛1−𝑎_𝑛2⋯∑_𝑘≠𝑛𝑎_𝑛𝑘. (A.1) Let 𝑐_𝑘 denote the cofactor of the 𝑘th diagonal element of 𝐿_𝐴, which has following property.

Lemma A.3 (see [9]). Assume 𝑛≥2. (1)If (𝒢,𝐴) is strongly connected, then 𝑐_𝑘>0 for 1≤𝑘≤𝑛. (2)The following identity holds: ^𝑛_𝑘,𝑗=1𝑐_𝑘𝑎_𝑘𝑗𝐺_𝑘𝑥_𝑘=^𝑛_𝑘,𝑗=1𝑐_𝑘𝑎_𝑘𝑗𝐺_𝑗𝑥_𝑗, (A.2)where 𝐺_𝑘(𝑥_𝑘),  1≤𝑘≤𝑛 are arbitrary functions.

Remark A.4. Lemma 2.1 of [8] gives the result (1) of Lemma A.3, and the result (2) is referred to [9, Theorem 2.3].

(2) Some theory about stationary distribution (see [23]): let 𝑋(𝑡) be a homogeneous Markov process in 𝐸_𝑙 (𝐸_𝑙 denotes Euclidean 𝑙-space) described by the stochastic equation 𝑑𝑋(𝑡)=𝑏(𝑋)𝑑𝑡+^𝑘_𝑟=1𝑔_𝑟(𝑋)𝑑𝐵_𝑟(𝑡). (A.3) The diffusion matrix is 𝐴(𝑥)=𝑎_𝑖𝑗(𝑥),𝑎_𝑖𝑗(𝑥)=^𝑘_𝑟=1𝑔^𝑖_𝑟(𝑥)𝑔^𝑗_𝑟(𝑥). (A.4)

Assumption B. There exists a bounded domain 𝑈⊂𝐸_𝑙 with regular boundary Γ, having the following properties.(B.1) In the domain 𝑈 and some neighbourhood thereof, the smallest eigenvalue of the diffusion matrix 𝐴(𝑥) is bounded away from zero.(B.2) If 𝑥∈𝐸_𝑙⧵𝑈, the mean time 𝜏 at which a path issuing from 𝑥 reaches the set 𝑈 is finite, and sup_𝑥∈𝐾𝐸_𝑥𝜏<∞ for every compact subset 𝐾⊂𝐸_𝑙.

Lemma A.5 (see [23]). If (B) holds, then the Markov process 𝑋(𝑡) has a stationary distribution 𝜇(⋅). Let 𝑓(⋅) be a function integrable with respect to the measure 𝜇. Then 𝑃_𝑥lim_𝑇→∞1𝑇^𝑇₀𝑓(𝑋(𝑡))𝑑𝑡=_𝐸𝑙𝑓(𝑥)𝜇(𝑑𝑥)=1 (A.5) for all 𝑥∈𝐸_𝑙.

Remark A.6. The proof is given in [23]. Exactly, the existence of stationary distribution with density is referred to Theorem 4.1, P. 119 and Lemma 9.4, P. 138. The weak convergence and the ergodicity are obtained in Theorem 5.1, P. 121 and Theorem 7.1, P. 130.
To validate (B.1), it suffices to prove 𝐹 is uniformly elliptical in 𝑈, where 𝐹𝑢=𝑏(𝑥)⋅𝑢_𝑥+(1/2)tr(𝐴(𝑥)𝑢_𝑥𝑥); that is, there is a positive number 𝑀 such that ^𝑘_𝑖,𝑗=1𝑎_𝑖𝑗(𝑥)𝜉_𝑖𝜉_𝑗≥𝑀𝜉²,𝑥∈𝑈,𝜉∈𝑅^𝑘 (A.6) (see Chapter 3, P. 103 of [26] and Rayleigh's principle in [27] Chapter 6, P. 349). To verify (B.2), it is sufficient to show that there exists some neighborhood 𝑈 and a nonnegative 𝐶²-function such that and for any 𝐸_𝑙⧵𝑈, 𝐿𝑉 is negative (for details refer to [28], P. 1163).

Acknowledgments

The work was supported by the Ministry of Education of China (no. 109051), NSFC of China (no. 10971021 and no. 11201008), and Program for Changjiang Scholars and Innovative Research Team in University.