Dynamics of a Stochastic SIS Epidemic Model with Saturated Incidence

Abstract

We introduce stochasticity into the SIS model with saturated incidence. The existence and uniqueness of the positive solution are proved by employing the Lyapunov analysis method. Then, we carry out a detailed analysis on both its almost sure exponential stability and its pth moment exponential stability, which indicates that the pth moment exponential stability implies the almost sure exponential stability. Additionally, the results show that the conditions for the disease to become extinct are much weaker than those in the corresponding deterministic model. The conditions for the persistence in the mean and the existence of a stationary distribution are also established. Finally, we derive the expressions for the mean and variance of the stationary distribution. Compared with the corresponding deterministic model, the threshold value for the disease to die out is affected by the half saturation constant. That is, increasing the saturation effect can reduce the disease transmission. Computer simulations are presented to illustrate our theoretical results.

1. Introduction

Epidemiology has been investigated by mathematicians through establishing mathematical models for a long time; see, for instance, [1–4]. Mathematical models can help people to better understand the spread of the disease and thus take effective measures to reduce its transmission as much as possible. Particularly, the classical SIS epidemic model [5, 6] is often used to model the dynamics of the diseases with no protective immunity such as gonorrhea.

Usually, the incidence function describes the number of new infections per unit time and it plays a vital role in the deterministic models. As it is known, in the most existing literature, the bilinear incidence rate 𝛽𝑆𝐼 is frequently used. However, this kind of incidence rate has some limitations, since it does not consider the behavioral change of susceptible individuals or control measures taken by the government when the number of infective individuals gets large. For example, people would wash their hands frequently, wear masks, and limit their time of going out, and the government would also take measures such as quarantine and isolation when the A/H1N1 influenza became serious in 2009 [7]. Therefore, it is more reasonable to adopt saturated incidence rate 𝛽𝑆𝐼/(1+ℎ𝐼) [8] in the mathematical models. When the number of infective individuals gets large, 𝛽𝐼/(1+ℎ𝐼) tends to a saturation effect, where 𝛽𝐼 is the infection force of the disease and 1/(1+ℎ𝐼) is the inhibition effect from the behavioral change or control measures. Here, ℎ represents the half saturation constant [9]. With the saturated incidence taken into consideration, the classical SIS epidemic model is transformed into 𝑑𝑆(𝑡)(𝑡)𝐼(𝑡)𝑑𝑡=𝜇𝑁−𝛽𝑆(𝑡)(𝑡)(𝑡),(𝑡)1+ℎ𝐼+𝛾𝐼−𝜇𝑆𝑑𝐼=(𝑡)𝐼(𝑡)𝑑𝑡𝛽𝑆(𝑡)−𝐼(𝑡),1+ℎ𝐼𝜇+𝛾 (1) where 𝑆(𝑡) and 𝐼(𝑡) are the number of susceptible individuals and infective individuals at time 𝑡, respectively. The total population 𝑁=𝑆(𝑡)+𝐼(𝑡) is a constant. Consider the following parameters:  𝛽 is the transmission rate, 𝜇 is the birth and death rate, and 𝛾 is the recovery rate. Assume all these parameters are positive; then, the classical stability analysis has shown that the system (1) always has a disease free equilibrium: 𝐸₀=𝑆₀,𝐼₀=(𝑁,0). (2) The basic reproduction number, 𝑅^𝑑₀=𝛽𝑁/(𝜇+𝛾), is a threshold value, which determines the extinction and persistence of the disease. When 𝑅^𝑑₀≤1, the disease free equilibrium 𝐸₀ is globally stable. When 𝑅^𝑑₀>1, the disease free equilibrium 𝐸₀ is unstable, and model (1) has a unique endemic equilibrium: 𝐸^∗=𝑆^∗,𝐼^∗=(1+ℎ𝑁)𝜇+𝛾,𝛽+ℎ𝜇+𝛾𝛽𝑁−𝜇+𝛾,𝛽+ℎ𝜇+𝛾 (3) which is globally stable.

Since the noise exists almost everywhere, epidemic models are inevitably affected by it. Hence, it is more reasonable to introduce random perturbations into mathematical models. At the same time, we note that taking the effect of environment noise on epidemic models into consideration has been a popular trend in disease spread modeling [10–13]. Gray et al. presented a stochastic SIS epidemic model in [14]. They established conditions for extinction and persistence of the disease and the existence of a stationary distribution. Stochastic SIR and SIRS models with and without distributed time delay were investigated in [15, 16], respectively, where the asymptotic behavior was discussed. In [17], stochastically perturbed SIR and SEIR epidemic models with saturated incidence were studied and the results of extinction and ergodicity were concluded. Ji et al. investigated a stochastic multigroup SIR model with fluctuations around the transmission coefficient and the death rate of each subpopulation in [18, 19] separately. Zhao et al. [20] discussed the extinction and persistence of the stochastic SIS epidemic model with vaccination. However, among all these researches, there are few literatures considering the parameter perturbation in SIS epidemic model with saturated incidence.

Motivated by the above consideration, in this paper, we introduce random effect into the SIS model (1) by replacing the transmission rate 𝛽 with 𝛽+𝜎𝑑𝑊(𝑡), where 𝑊(𝑡) is standard Brownian motion and 𝜎 is the intensity of white noise. The stochastic version corresponding to the deterministic model (1) can be represented as (𝑡)=(𝑡)𝐼(𝑡)𝑑𝑆𝜇𝑁−𝛽𝑆(𝑡)(𝑡)(𝑡)−(𝑡)𝐼(𝑡)1+ℎ𝐼+𝛾𝐼−𝜇𝑆𝑑𝑡𝜎𝑆(𝑡)(𝑡),(𝑡)=(𝑡)𝐼(𝑡)1+ℎ𝐼𝑑𝑊𝑑𝐼𝛽𝑆(𝑡)−𝐼(𝑡)(𝑡)𝐼(𝑡)1+ℎ𝐼𝜇+𝛾𝑑𝑡+𝜎𝑆(𝑡)(𝑡).1+ℎ𝐼𝑑𝑊 (4) Given that 𝑆(𝑡)+𝐼(𝑡)=𝑁, it is sufficient to study the following equation: (𝑡)=𝛽((𝑡𝐼(𝑡)𝑑𝐼𝑁−𝐼))(𝑡)−𝐼(𝑡)+𝜎((𝑡𝐼(𝑡)1+ℎ𝐼𝜇+𝛾𝑑𝑡𝑁−𝐼))(𝑡)(𝑡),1+ℎ𝐼𝑑𝑊 (5) with initial value 𝐼(0)=𝐼₀∈(0,𝑁). In the following sections, we will restrict ourselves to (5).

The paper is organized as follows. In Section 2, the existence of unique positive solution is shown. Both its almost sure exponential stability and its 𝑝th moment exponential stability are then investigated in Section 3. The conditions for persistence in the mean of the disease are established in Section 4. The existence of a stationary distribution and the expressions for the mean and variance are presented in Section 5. In Section 6, we give a brief conclusion. Besides, the computer simulations which support our results are given in each section. Finally, we give an appendix containing some theory used in the previous sections.

2. Existence of Unique Positive Solution

To investigate the dynamical behavior of the epidemic model, we need to show that the model has a unique global solution and the solution will remain within (0,𝑁) whenever it starts there. Hence, in this section, employing the Lyapunov analysis method (mentioned in [18, 21]), we establish Theorem 1.

Theorem 1. There is a unique global solution 𝐼(𝑡)∈(0,𝑁) of system (5) on 𝑡≥0 for any given initial value 𝐼(0)=𝐼₀∈(0,𝑁) with probability 1; namely, 𝑃{𝐼(𝑡)∈(0,𝑁),∀𝑡≥0}=1. (6)

Proof. Since the coefficients of (5) are locally Lipschitz continuous for any given initial value 𝐼(0)=𝐼₀∈(0,𝑁), there is a unique local solution 𝐼(𝑡) on 𝑡∈[0,𝜏_𝑒), where 𝜏_𝑒 is the explosion time [22]. To show that the solution is global, we need to show that 𝜏_𝑒=∞ a.s. Let 𝑚₀>0 be sufficiently large such that 1/𝑚₀<𝐼₀<𝑁−1/𝑚₀. For each integer 𝑚≥𝑚₀, define the stopping time 𝜏_𝑚=inf𝑡∈0,𝜏_𝑒(𝑡)∉1∶𝐼𝑚1,𝑁−𝑚, (7) where throughout this paper we set inf∅=∞ (as usual ∅ denotes the empty set). Clearly, 𝜏_𝑚 is increasing as 𝑚→∞. Set 𝜏_∞=lim_𝑚→∞𝜏_𝑚, whence 𝜏_∞≤𝜏_𝑒 a.s. If we can show that 𝜏_∞=∞ a.s., then 𝜏_𝑒=∞ a.s. and 𝐼(𝑡)∈(0,𝑁) a.s. for all 𝑡≥0. In other words, to complete the proof, all what we need to show is that 𝜏_∞=∞ a.s. If this statement is not true, then there exists a pair of constants 𝑇>0 and 𝜖∈(0,1) such that 𝑃𝜏_∞≤𝑇>𝜖. (8) Therefore, there is an integer 𝑚₁≥𝑚₀ such that 𝑃𝜏_𝑚≤𝑇≥𝜖,∀𝑚≥𝑚₁. (9) Define a function 𝑉∶(0,𝑁)→𝑅₊ by 𝑉(𝑥)=1𝑥+1.𝑁−𝑥 (10) Using the Itô formula, we get (𝑥)=−1𝑑𝑉𝑥²+1(𝑁−𝑥)²1𝑑𝑥+𝑥³+1(𝑁−𝑥)³(𝑑𝑥)²(𝑥)+−𝜎()=𝐿𝑉𝑑𝑡𝑁−𝑥𝑥()+1+ℎ𝑥𝜎𝑥()(𝑡),𝑁−𝑥)(1+ℎ𝑥𝑑𝑊 (11) where 𝐿𝑉∶(0,𝑁)→𝑅 is defined by −1𝐿𝑉=𝑥²+1(𝑁−𝑥)²𝛽(𝑁−𝑥)𝑥−𝑥+11+ℎ𝑥𝜇+𝛾𝑥³+1(𝑁−𝑥)³𝜎²(𝑁−𝑥)²𝑥²(1+ℎ𝑥)²𝛽()=−𝑁−𝑥𝑥()+1+ℎ𝑥𝜇+𝛾𝑥+𝛽𝑥()+𝜎𝑁−𝑥)(1+ℎ𝑥₂(𝑁−𝑥)²𝑥(1+ℎ𝑥)²+𝜎²𝑥²(𝑁−𝑥)(1+ℎ𝑥)²−𝜇+𝛾𝑥(𝑁−𝑥)²≤𝜇+𝛾𝑥+𝛽𝑥()+𝜎𝑁−𝑥)(1+ℎ𝑥²(𝑁−𝑥)²𝑥(1+ℎ𝑥)²+𝜎²𝑥²(𝑁−𝑥)(1+ℎ𝑥)²≤1𝑥𝜇+𝛾+𝜎²𝑁²+1𝛽𝑁−𝑥ℎ+𝜎²ℎ²≤𝐶𝑉(𝑥), (12) where 𝐶=(𝜇+𝛾+𝜎²𝑁²)∨(𝛽/ℎ+𝜎²/ℎ²). Therefore, ^𝜏𝑚₀^∧𝑡(𝐼(𝑠=𝑑𝑉))_{𝜏𝑚0∧𝑡}(𝐼(𝑠+𝐿𝑉))𝑑𝑠_{𝜏𝑚0∧𝑡}−𝜎(𝑁−𝐼(𝑠))𝐼(𝑠(𝑠+(𝑠))(1+ℎ𝐼))𝜎𝐼((𝑠(𝑠(𝑠)≤𝑁−𝐼))(1+ℎ𝐼))𝑑𝑊_{𝜏𝑚0∧𝑡}(𝐼(𝑠+𝐶𝑉))𝑑𝑠_{𝜏𝑚0∧𝑡}−𝜎(𝑁−𝐼(𝑠))𝐼(𝑠(𝑠+(𝑠))(1+ℎ𝐼))𝜎𝐼((𝑠(𝑠(𝑠),𝑁−𝐼))(1+ℎ𝐼))𝑑𝑊 (13) which implies that 𝐸𝑉𝐼𝜏_𝑚∧𝑡≤𝑉𝐼₀+𝐶𝐸^𝜏𝑚₀^∧𝑡𝑉(𝐼(𝑠𝐼))𝑑𝑠≤𝑉₀+𝐶^𝑡₀𝐸𝑉𝐼𝜏_𝑚∧𝑠𝑑𝑠. (14) The Gronwall inequality yields that 𝐸𝑉𝐼𝜏_𝑚∧𝑇≤𝑉𝐼₀𝑒^𝐶𝑇. (15) Set Ω_𝑚={𝜏_𝑚≤𝑇} for all 𝑚≥𝑚₁. Then, by (9), 𝑃(Ω_𝑚)≥𝜖. For every 𝜔∈Ω_𝑚, 𝐼(𝜏_𝑚,𝜔) equals either 1/𝑚 or 𝑁−1/𝑚, and therefore 𝑉𝐼𝜏_𝑚,𝜔≥𝑚. (16) It then follows from (9) and (15) that 𝑉𝐼₀𝑒^𝐶𝑇≥𝐸1_Ω𝑚(𝜔)𝑉𝐼𝜏_𝑚,𝜔≥𝜖𝑚, (17) where 1_Ω𝑚(𝜔) is the indicator function of Ω_𝑚. Letting 𝑚→∞ leads to the contradiction that ∞>𝑉𝐼₀𝑒^𝐶𝑇=∞. (18) So, we must have 𝜏_∞=∞ a.s., whence the proof is complete.

3. Extinction

In discussing the extinction of system (5), we focus on the two kinds of exponential stabilities, almost sure exponential stability and 𝑝th moment exponential stability.

3.1. Almost Sure Exponential Stability

Theorem 2. Let 𝐼(𝑡) be the solution of system (5) with initial value 𝐼(0)=𝐼₀∈(0,𝑁). If (a)𝑅^𝑠₀₁∶=𝛽𝑁/(𝜇+𝛾)−𝜎²𝑁²/(2(𝜇+𝛾)(1+ℎ𝑁)²)=𝑅^𝑑₀−𝜎²𝑁²/(2(𝜇+𝛾)(1+ℎ𝑁)²)<1 and 𝜎²≤𝛽(1+ℎ𝑁)²/𝑁,(b)𝜎²>𝛽(1+ℎ𝑁)²/𝑁∨𝛽²(1+ℎ𝑁)²/2(𝜇+𝛾).then, limsup_𝑡→∞(log(𝐼(𝑡))/𝑡)≤𝛽𝑁−𝜇−𝛾−𝜎²𝑁²/2(1+ℎ𝑁)²<0   𝑎.𝑠. if (a) holds; limsup_𝑡→∞(log(𝐼(𝑡))/𝑡)≤𝛽²(1+ℎ𝑁)²/2𝜎²−𝜇−𝛾<0  𝑎.𝑠. if (b) holds;that is, 𝐼(𝑡) tends to zero exponentially 𝑎.𝑠.; namely, the disease will die out with probability one.

Proof. Applying Itô formula to system (5) leads to (𝐼(𝑡(𝐼(𝑡𝜎((𝑡𝑑log))=𝐹))𝑑𝑡+𝑁−𝐼))(𝑡)(𝑡),1+ℎ𝐼𝑑𝑊 (19) where 𝐹(𝐼(𝑡))=𝛽(𝑁−𝐼(𝑡))/(1+ℎ𝐼(𝑡))−𝜇−𝛾−𝜎²(𝑁−𝐼(𝑡))²/2(1+ℎ𝐼(𝑡))². Integrating both sides of (19) from 0 to 𝑡, we have log(𝐼(𝑡))=log𝐼₀+^𝑡₀𝐹(𝐼(𝑠))𝑑𝑠+𝑀(𝑡), (20) where 𝑀(𝑡)=∫^𝑡₀(𝜎(𝑁−𝐼(𝑠))/(1+ℎ𝐼(𝑠)))𝑑𝑊(𝑠). Obviously, 𝐹(𝐼(𝑠(𝑠)𝜎))≤𝛽𝑁−𝛽𝐼−𝜇−𝛾−²(𝑁−𝐼(𝑠))²2(1+ℎ𝑁)²∶=𝐺(𝐼(𝑠)), (21) where 𝐺(𝐼(𝑠))=𝛽𝑁−𝜇−𝛾−𝜎²𝑁²/2(1+ℎ𝑁)²−(𝛽−𝜎²𝑁/(1+ℎ𝑁)²)𝐼(𝑠)−(𝜎²/2(1+ℎ𝑁)²)𝐼(𝑠)², for 𝐼(𝑠)∈(0,𝑁).
If condition (a) is satisfied, it then follows from (20) that log(𝐼(𝑡))≤log𝐼₀+𝜎𝛽𝑁−𝜇−𝛾−²𝑁²2(1+ℎ𝑁)²𝑡+𝑀(𝑡). (22) This implies that limsup_𝑡→∞log(𝐼(𝑡))𝑡𝜎≤𝛽𝑁−𝜇−𝛾−²𝑁²2(1+ℎ𝑁)²+limsup_𝑡→∞𝑀(𝑡)𝑡a.s. (23) Clearly, 𝑀(𝑡)=∫^𝑡₀(𝜎(𝑁−𝐼(𝑠))/(1+ℎ𝐼(𝑠)))𝑑𝑊(𝑠) is a local martingale and vanishes at 𝑡=0. Moreover, limsup_𝑡→∞⟨𝑀,𝑀⟩_𝑡𝑡≤𝜎²𝑁²<∞a.s. (24) By strong law of large numbers [22], we get lim_𝑡→∞𝑀(𝑡)𝑡=0a.s. (25) Then, limsup_𝑡→∞log(𝐼(𝑡))𝑡𝜎≤𝛽𝑁−𝜇−𝛾−²𝑁²2(1+ℎ𝑁)²<0a.s. (26)
If condition (b) is satisfied, the function 𝐺(𝑥) takes its maximum value 𝐺(𝑥) at 𝑥=𝑁−𝛽(1+ℎ𝑁)²/𝜎²∈(0,𝑁)); that is, 𝐺(𝑥)=𝛽²(1+ℎ𝑁)²/2𝜎²−(𝜇+𝛾). It then follows from (20) that log(𝐼(𝑡))≤log𝐼₀+𝛽²(1+ℎ𝑁)²2𝜎²−𝜇+𝛾𝑡+𝑀(𝑡). (27) Then, limsup_𝑡→∞log(𝐼(𝑡))𝑡≤𝛽²(1+ℎ𝑁)²2𝜎²−𝜇−𝛾<0a.s. (28) This finishes the proof.

From Theorem 2, we know that when 𝑅^𝑠₀₁=𝑅^𝑑₀−𝜎²𝑁²/2(𝜇+𝛾)(1+ℎ𝑁)²<1 and the perturbation is small, the disease will die out a.s., whereas the disease may still persist for the corresponding deterministic model (1). If the noise intensity 𝜎² is larger than 𝛽(1+ℎ𝑁)²/𝑁∨𝛽²(1+ℎ𝑁)²/2(𝜇+𝛾), then the disease will also die out a.s. For ℎ=0, system (5) becomes the model discussed in [14] and the conditions for the disease to become extinct we obtain here are in agreement with those derived in [14]. If there exists no environment noise, that is, 𝜎=0, then 𝑅^𝑠₀₁=𝑅^𝑑₀ is the basic reproduction number of the deterministic model (1). In addition, the half saturation constant also influences 𝑅^𝑠₀₁, which shows that increasing the saturation effect can reduce the value of 𝑅^𝑠₀₁ and thus can reduce the disease spread. The following simulations indicate these results.

Example 3. We assume that the unit time is one day. The parameters are given by 𝑁=100,𝛽=0.013,𝜇=0.4,𝛾=0.6,ℎ=0.05,𝜎=0.06. (29) Note that 𝑅^𝑠₀₁=𝛽𝑁−𝜎𝜇+𝛾²𝑁²2𝜇+𝛾(1+ℎ𝑁)²𝜎=0.8<1,₂𝛽()=0.0036≤1+ℎ𝑁²𝑁=0.0047. (30) That is, condition (a) is satisfied and the solution 𝐼(𝑡) satisfies limsup_𝑡→∞log(𝐼(𝑡))𝑡𝜎≤𝛽𝑁−𝜇−𝛾−²𝑁²2(1+ℎ𝑁)²=−0.2<0a.s. (31) Hence, 𝐼(𝑡) tends to zero exponentially a.s.; that is, the disease will die out a.s.
On the other hand, for the deterministic model (1), the basic reproduction number is 𝑅^𝑑₀=𝛽𝑁𝜇+𝛾=1.3>1. (32) The solution 𝐼(𝑡) obeys lim_𝑡→∞𝐼(𝑡)=𝛽𝑁−𝜇+𝛾𝛽+𝜇+𝛾ℎ=4.7619. (33) Employing the Euler-Maruyama (EM) method in [23], the computer simulations are given in Figure 1, which support our results.

(a)

(b)

Figure 1 

Computer simulations of 𝐼(𝑡) for model (5) and the deterministic model (1), using the EM method with step size 0.001 and with initial values (a) 𝐼(0)=80 and (b) 𝐼(0)=10.

Example 4. We keep the parameters the same as Example 3 but increase 𝜎 to 0.09. Note that 𝜎²𝛽()=0.0081>max1+ℎ𝑁²𝑁𝛽=0.00468,²(1+ℎ𝑁)²2𝜇+𝛾=0.00304=0.00468. (34) That is, condition (b) is satisfied and limsup_𝑡→∞log(𝐼(𝑡))𝑡≤𝛽²(1+ℎ𝑁)²2𝜎²−𝜇−𝛾=−0.624<0a.s. (35) Therefore, 𝐼(𝑡) tends to zero exponentially a.s.; namely, the disease will die out a.s. The simulations are shown in Figure 2 to support our results completely.

(a)

(b)

Figure 2 

Computer simulations of 𝐼(𝑡) for model (5) and the deterministic model (1), using the EM method with step size 0.001 and with initial values (a) 𝐼(0)=80 and (b) 𝐼(0)=10.

3.2. 𝑝th Moment Exponential Stability

Following (20), we obtain 𝐼(𝑡)=𝐼₀exp^𝑡₀𝐹(𝐼(𝑠))𝑑𝑠+𝑀(𝑡). (36) As 𝑊(𝑡)∼𝑁(0,𝑡), then [24] 𝐸𝑒^𝑀(𝑡)=𝐸exp^𝑡₀𝜎(𝑁−𝐼(𝑠))(𝑠)(𝑠)1+ℎ𝐼𝑑𝑊≤𝐸𝑒^{𝜎𝑁𝑊(𝑡)}=𝑒^{(𝜎2𝑁2/2)𝑡}. (37) Using the fact that 𝐹(𝐼(𝑠))≤𝛽𝑁−𝜇−𝛾, (38) we find 𝐸𝐼(𝑡)=𝐸𝐼₀exp^𝑡₀𝐹(𝐼(𝑠))𝑑𝑠+𝑀(𝑡)≤𝐼₀𝜎exp𝛽𝑁−𝜇−𝛾+²𝑁²2𝑡. (39) Since log(𝐼(𝑡))^𝑝(𝐼(𝑡𝐼=𝑝log))≤𝑝log₀+𝛽𝑁−𝜇−𝛾𝑡+𝜎𝑁𝑊(𝑡), (40) then 𝐸(𝐼(𝑡))^𝑝≤𝐼₀^𝑝⋅𝑒^{𝑝(𝛽𝑁−𝜇−𝛾)𝑡}⋅𝑒^{(𝑝2𝜎2𝑁2/2)𝑡}=𝐼₀^𝑝exp𝑝𝛽𝑁−𝜇−𝛾+𝑝𝜎²𝑁²2𝑡. (41)

Theorem 5. Assume that 𝛽𝑁−𝜇−𝛾+(𝑝𝜎²𝑁²/2)<0; the solution 𝐼(𝑡) with initial value 𝐼(0)=𝐼₀∈(0,𝑁) is 𝑝th moment exponentially stable.

Next, we will make a comparison between the almost sure exponential stability and the 𝑝th moment exponential stability, since 𝑥𝛽(𝑁−𝑥)𝑥−𝑥=𝛽()1+ℎ𝑥𝜇+𝛾𝑁−𝑥−𝑥1+ℎ𝑥𝜇+𝛾²≤𝛽𝑁−𝜇+𝛾𝑥². (42) From Theorem 5, we know that 𝛽𝑁−(𝜇+𝛾)<0; that is, 𝑥𝛽(𝑁−𝑥)𝑥−𝑥≤𝑥1+ℎ𝑥𝜇+𝛾𝛽𝑁−𝜇+𝛾²≤0. (43) Moreover, 𝜎²(𝑁−𝑥)²𝑥²(1+ℎ𝑥)²≤𝜎²(𝑁−𝑥)²𝑥²≤𝜎²𝑁²𝑥². (44) Let 𝐾=𝜎²𝑁² be positive; then 𝑥𝛽(𝑁−𝑥)𝑥−𝑥∨𝜎1+ℎ𝑥𝜇+𝛾²(𝑁−𝑥)²𝑥²(1+ℎ𝑥)²≤𝐾𝑥². (45) Therefore, the condition of Lemma A.4 (in the appendix) [22] is satisfied, which suggests that the 𝑝th moment exponential stability of system (5) implies the almost sure exponential stability.

4. Persistence

We will investigate persistence in the mean for system (5) in this section. The definition of the persistence in the mean was initially proposed for the deterministic system [25] and was also used for the stationary stochastic system with ergodicity [19, 20]. According to the ergodicity for the stationary process, the time average in the long time limit is equal to the ensemble expectation over phase space, so we introduce the notation ⟨𝐼(𝑡=1)⟩𝑡^𝑡₀𝐼(𝑠)𝑑𝑠. (46)

Definition 6. System (5) is said to be persistent in the mean, if liminf_𝑡→∞1𝑡^𝑡₀𝐼(𝑠)𝑑𝑠>0a.s. (47)

Theorem 7. If 𝑅^𝑠₀₂∶=𝛽𝑁()−𝜎𝜇+𝛾1+ℎ𝑁²𝑁²2𝜎𝜇+𝛾>1,₂<𝛽(1+ℎ𝑁)²𝑁, (48) then, for any initial value 𝐼(0)=𝐼₀∈(0,𝑁), the solution 𝐼(𝑡) of system (5) has the following property: 𝐼≤liminf_𝑡→∞⟨𝐼(𝑡)⟩≤limsup_𝑡→∞⟨𝐼(𝑡)⟩≤𝐼a.s., (49) where 𝜎𝐼=𝛽𝑁−𝜇−𝛾−²𝑁²/2(1+ℎ𝑁)²𝛽−𝜎²𝑁/(1+ℎ𝑁)²,𝐼=𝛽𝑁/(1+ℎ𝑁)−𝜇−𝛾−𝜎²𝑁²/2().𝛽/1+ℎ𝑁 (50)

Proof. By (20), we obtain log(𝐼(𝑡))−log𝐼₀𝑡≤1𝑡^𝑡₀(𝑠)𝜎𝛽𝑁−𝛽𝐼−𝜇−𝛾−²(𝑁−𝐼(𝑠))²2(1+ℎ𝑁)²+𝑀(𝑡)𝑑𝑠𝑡𝜎=𝛽𝑁−𝜇−𝛾−²𝑁²2(1+ℎ𝑁)²⟨𝐼(𝑡+𝜎−𝛽)⟩₂𝑁(1+ℎ𝑁)²⟨𝐼(𝑡−𝜎)⟩²2(1+ℎ𝑁)²𝐼(𝑡)²+𝑀(𝑡)𝑡𝜎≤𝛽𝑁−𝜇−𝛾−²𝑁²2(1+ℎ𝑁)²−𝜎𝛽−²𝑁(1+ℎ𝑁)²⟨𝐼(𝑡+𝑀(𝑡))⟩𝑡. (51) The inequality can be written as log(𝐼(𝑡))𝑡𝜎≤𝛽𝑁−𝜇−𝛾−²𝑁²2(1+ℎ𝑁)²−𝜎𝛽−²𝑁(1+ℎ𝑁)²⟨𝐼(𝑡+𝑀(𝑡))⟩𝑡+log𝐼₀𝑡. (52) Since 𝑅^𝑠₀₁≥𝑅^𝑠₀₂, if condition (48) is satisfied, then 𝜎𝛽𝑁−𝜇−𝛾−²𝑁²2(1+ℎ𝑁)²𝜎>0,𝛽−²𝑁(1+ℎ𝑁)²>0. (53) Together with Lemma A.2. in [20], we have lim_𝑡→∞⟨𝐼(𝑡≤𝜎sup)⟩𝛽𝑁−𝜇−𝛾−²𝑁²/2(1+ℎ𝑁)²𝛽−𝜎²𝑁/(1+ℎ𝑁)²=𝐼a.s. (54) On the other hand, log(𝐼(𝑡))−log𝐼₀𝑡≥1𝑡^𝑡₀𝛽(𝑁−𝐼(𝑠))11+ℎ𝑁−𝜇−𝛾−2𝜎²(𝑁−𝐼(𝑠))²+𝑀(𝑡)𝑑𝑠𝑡≥𝛽𝑁11+ℎ𝑁−𝜇−𝛾−2𝜎²𝑁²−𝛽(𝑡)𝑀(𝑡)1+ℎ𝑁⟨𝐼⟩+𝑡. (55) That is, ⟨𝐼(𝑡≥1)⟩𝛽/(1+ℎ𝑁)𝛽𝑁11+ℎ𝑁−𝜇−𝛾−2𝜎²𝑁²+𝑀(𝑡)𝑡−log(𝐼(𝑡))−log𝐼₀𝑡. (56) Since −∞<log(𝐼(𝑡))<log(𝑁), then liminf_𝑡→∞⟨𝐼(𝑡≥()𝜎)⟩𝛽𝑁/1+ℎ𝑁−𝜇−𝛾−²𝑁²/2𝛽/(1+ℎ𝑁)=𝐼∈(0,𝑁)a.s. (57) Moreover, 𝜎𝐼=𝛽𝑁−𝜇−𝛾−²𝑁²/2(1+ℎ𝑁)²𝛽−𝜎²𝑁/(1+ℎ𝑁)²≥𝐼=𝛽𝑁/(1+ℎ𝑁)−𝜇−𝛾−𝜎²𝑁²/2().𝛽/1+ℎ𝑁 (58) Therefore, we complete the proof.

Example 8. The parameters are given by 𝑁=100,𝛽=0.2,𝜇=0.4,𝛾=0.6,ℎ=0.05,𝜎=0.02. (59) Note that 𝑅^𝑠₀₂∶=𝛽𝑁()−𝜎𝜇+𝛾1+ℎ𝑁²𝑁²2𝜎𝜇+𝛾=1.33>1,₂𝛽()=0.0004<1+ℎ𝑁²𝑁=0.072. (60) Then the solution 𝐼(𝑡) of system (5) satisfies 10≤liminf_𝑡→∞⟨𝐼(𝑡)⟩≤limsup_𝑡→∞⟨𝐼(𝑡)⟩≤95.2938a.s. (61) For the corresponding deterministic model (1), 𝑅^𝑑₀=𝛽𝑁𝜇+𝛾=20>1,lim_𝑡→∞𝐼(𝑡)=𝛽𝑁−𝜇+𝛾𝛽+ℎ𝜇+𝛾=76. (62) The simulation results are shown in Figure 3.

(a)

(b)

Figure 3 

Computer simulations of 𝐼(𝑡) for model (5) (red line) and the deterministic model (1) (blue line), using the EM method with step size 0.001 and with initial values (a) 𝐼(0)=80 and (b) 𝐼(0)=10. 𝐼=10 and 𝐼=95.2938.

5. Stationary Distribution

Applying Lemma A.5 (in the appendix) [26], if we can show that conditions (B.1) and (B.2) are satisfied, then system (5) has a stationary distribution. Obviously, the square of the diffusion coefficient of system (5), that is, 𝜎²(𝑁−𝐼)²𝐼²/(1+ℎ𝐼)², is bounded away from zero for 𝐼∈(𝑎,𝑏)⊂(0,𝑁)); then, condition (B.1) is satisfied. Next, we will prove that (B.2) is also valid.

Theorem 9. Assume that 𝑅^𝑠₀₂∶=𝛽𝑁/((𝜇+𝛾)(1+ℎ𝑁))−𝜎²𝑁²/(2(𝜇+𝛾))>1 and 𝜎²<(2𝛽/(𝑁²(1+ℎ𝐼^∗)𝐼^∗))min{(𝐼^∗)²,(𝑁−𝐼^∗)²}; then, there is a unique stationary distribution of system (5). Here, 𝐼^∗ is the unique endemic equilibrium of system (1); that is, 𝐼^∗=(𝛽𝑁−(𝜇+𝛾))/(𝛽+ℎ(𝜇+𝛾)). Moreover, lim_𝑡→∞1𝑡𝐸^𝑡₀𝛽1+ℎ𝐼^∗𝐼(𝑠)−𝐼^∗21𝑑𝑠≤2𝐼^∗𝜎²𝑁². (63)

Proof. If 𝑅^𝑠₀₂∶=𝛽𝑁/((𝜇+𝛾)(1+ℎ𝑁))−𝜎²𝑁²/2(𝜇+𝛾)>1, it is clear that 𝑅^𝑑₀=𝛽𝑁/(𝜇+𝛾)>1; then, model (1) has a unique endemic equilibrium 𝐼^∗; that is, 𝛽𝑁−𝐼^∗1+ℎ𝐼^∗=𝜇+𝛾,𝐼^∗=𝛽𝑁−𝜇+𝛾.𝛽+ℎ𝜇+𝛾 (64) Define 𝑉(𝐼)=𝐼−𝐼^∗−𝐼^∗𝐼log𝐼^∗. (65) Then 𝑉 is positive definite. By Itô formula, we yield 𝐼𝑑𝑉=1−𝐼^∗𝐼𝑑𝐼+^∗2𝐼²(𝑑𝐼)²𝜎()=𝐿𝑉𝑑𝑡+𝑁−𝐼𝐼−𝐼^∗(𝑡),1+ℎ𝐼𝑑𝑊 (66) where 𝛽()𝐿𝑉=𝑁−𝐼−1+ℎ𝐼𝜇+𝛾𝐼−𝐼^∗+𝐼^∗𝜎²(𝑁−𝐼)²2(1+ℎ𝐼)²=𝛽(𝑁−𝐼)−𝛽1+ℎ𝐼𝑁−𝐼^∗1+ℎ𝐼^∗𝐼−𝐼^∗+𝐼^∗𝜎²(𝑁−𝐼)²2(1+ℎ𝐼)²𝛽(=−1+ℎ𝑁)𝐼−𝐼^∗2(1+ℎ𝐼)(1+ℎ𝐼^∗)+𝐼^∗𝜎²(𝑁−𝐼)²2(1+ℎ𝐼)²𝛽≤−1+ℎ𝐼^∗𝐼−𝐼^∗2+12𝐼^∗𝜎²𝑁². (67) Note that if 𝜎²<(2𝛽/𝑁²(1+ℎ𝐼^∗)𝐼^∗)min{(𝐼^∗)²,(𝑁−𝐼^∗)²}, then the neighborhood 𝑈∶−(𝛽/(1+ℎ𝐼^∗))(𝐼−𝐼^∗)²+(1/2)𝐼^∗𝜎²𝑁²=0 lies entirely in (0,𝑁). Hence, for 𝐼∈(0,𝑁)⧵𝑈, 𝐿𝑉≤−𝑘 (𝑘 is a positive constant), which implies that condition (B.2) is satisfied; for more details, we refer to [27]. Therefore, system (5) has a stationary distribution.
Moreover, 0≤𝐸^𝑡₀𝑑𝑉(𝐼(𝑠))=𝐸^𝑡₀𝐿𝑉(𝐼(𝑠))𝑑𝑠≤−𝐸^𝑡₀𝛽1+ℎ𝐼^∗𝐼(𝑠)−𝐼^∗21𝑑𝑠+2𝐼^∗𝜎²𝑁²𝑡. (68) Then lim_𝑡→∞1𝑡𝐸^𝑡₀𝛽1+ℎ𝐼^∗𝐼(𝑠)−𝐼^∗21𝑑𝑠≤2𝐼^∗𝜎²𝑁². (69)

Next, we will give the expressions for the mean and variance of the stationary distribution.

Theorem 10. If 𝑅^𝑠₀₂∶=𝛽𝑁/((𝜇+𝛾)(1+ℎ𝑁))−𝜎²𝑁²/(2(𝜇+𝛾)) >1 and 𝜎²<(2𝛽/(𝑁²(1+ℎ𝐼^∗)𝐼^∗))min{(𝐼^∗)²,(𝑁−𝐼^∗)²}. Let 𝑚 and 𝑣 denote the mean and variance of the stationary distribution of model (5). Then, 1𝑚=𝛽𝑁−𝜇−𝛾−2𝜎²𝑁²×𝛽ℎ+𝜇+𝛾ℎ²+12𝜎²𝛽𝑁−𝜇−𝛾ℎ−𝛽()𝛽+𝜇+𝛾ℎ𝑁−1−2ℎ𝜇+𝛾+𝜎²𝑁⁻¹,𝑣=𝛽𝑁−𝜇−𝛾𝛽+𝜇+𝛾ℎ𝑚−𝑚². (70)

Proof. Multiply both sides of (5) by 1+ℎ𝐼(𝑡); that is, ((𝑡(𝑡)=𝛽((𝑡𝐼(𝑡)−𝐼(𝑡(𝑡((𝑡𝐼(𝑡)(𝑡).1+ℎ𝐼))𝑑𝐼𝑁−𝐼))𝜇+𝛾)(1+ℎ𝐼))𝑑𝑡+𝜎𝑁−𝐼))𝑑𝑊 (71) Integrating (71) from 0 to 𝑡 and dividing by 𝑡 on both sides, then 1𝑡^𝑡₀((𝑠(𝑠)=11+ℎ𝐼))𝑑𝐼𝑡^𝑡₀𝛽((𝑠𝐼(𝑠)−𝐼(𝑠(𝑠+1𝑁−𝐼))𝜇+𝛾)(1+ℎ𝐼))𝑑𝑠𝑡^𝑡₀𝜎(𝑁−𝐼(𝑠))𝐼(𝑠)𝑑𝑊(𝑠). (72) Since 𝐼(𝑡)−𝐼₀𝑡≤1𝑡^𝑡₀((𝑠(𝑠)≤()𝐼(𝑡)1+ℎ𝐼))𝑑𝐼1+ℎ𝑁−𝐼₀𝑡, (73) letting 𝑡→∞, we have lim_𝑡→∞1𝑡^𝑡₀(1+ℎ𝐼(𝑠))𝑑𝐼(𝑠)=0a.s. (74) Then, applying the ergodic property of the stationary distribution and strong law of large numbers, we get 𝛽𝑁−𝜇−𝛾𝑚−𝛽+𝜇+𝛾ℎ𝐼(𝑡)²=0. (75) Here, ⟨𝐼(𝑡)²⟩ denotes the second moment of the stationary distribution [14].
In a similar way, multiplying both sides of (19) by (1+ℎ𝐼(𝑡))², we obtain that (1+ℎ𝐼(𝑡))²𝑑log(𝐼(𝑡))=(1+ℎ𝐼(𝑡))²𝐹(𝐼(𝑡((𝑡(𝑡(𝑡).))𝑑𝑡+𝜎𝑁−𝐼))(1+ℎ𝐼))𝑑𝑊 (76) Then, integrating (76) from 0 to 𝑡 and dividing by 𝑡 on both sides, 1𝑡^𝑡₀(1+ℎ𝐼(𝑠))²(𝐼(𝑠=1𝑑log))𝑡^𝑡₀(1+ℎ𝐼(𝑠))²𝐹(𝐼(𝑠+1))𝑑𝑠𝑡^𝑡₀𝜎(𝑁−𝐼(𝑠))(1+ℎ𝐼(𝑠))𝑑𝑊(𝑠). (77) Since log𝐼(𝑡)−log𝐼₀𝑡≤1𝑡^𝑡₀(1+ℎ𝐼(𝑠))²(𝐼(𝑠≤()𝑑log))1+ℎ𝑁²log𝐼(𝑡)−log𝐼₀𝑡, (78) then lim_𝑡→∞1𝑡^𝑡₀(1+ℎ𝐼(𝑠))²𝑑log(𝐼(𝑠))=0a.s. (79) Namely, lim_𝑡→∞1𝑡^𝑡₀(1+ℎ𝐼(𝑠))²𝐹(𝐼(𝑠))𝑑𝑠+lim_𝑡→∞1𝑡^𝑡₀𝜎(𝑁−𝐼(𝑠))(1+ℎ𝐼(𝑠))𝑑𝑊(𝑠)=0a.s. (80) Therefore, 1𝛽𝑁−𝜇−𝛾−2𝜎²𝑁²+𝛽(ℎ𝑁−1)−2ℎ𝜇+𝛾+𝜎²𝑁𝑚−𝛽ℎ+𝜇+𝛾ℎ²+12𝜎²𝐼(𝑡)²=0. (81) From (75) and (81), we can get 1𝑚=𝛽𝑁−𝜇−𝛾−2𝜎²𝑁²×𝛽ℎ+𝜇+𝛾ℎ²+12𝜎²𝛽𝑁−𝜇−𝛾ℎ−𝛽()𝛽+𝜇+𝛾ℎ𝑁−1−2ℎ𝜇+𝛾+𝜎²𝑁⁻¹,𝑣=𝛽𝑁−𝜇−𝛾𝛽+𝜇+𝛾ℎ𝑚−𝑚². (82)

Example 11. We keep the parameters the same as Example 8 but reduce 𝜎 to 0.005, 0.001, and 0.0005. Note that 𝑅^𝑠₀₂ is 3.2083, 3.3283, and 3.3321, respectively, and all of them are bigger than 1. Furthermore, 𝐼^∗=𝛽𝑁−𝜇+𝛾𝛽+ℎ𝜇+𝛾=76,2𝛽𝑁²(1+ℎ𝐼^∗)𝐼^∗min𝐼^∗2,𝑁−𝐼^∗2=0.000063. (83) The values of 𝜎 are satisfied for the conditions of Theorem 10. The histograms of 𝐼(𝑡) are shown in Figure 4, respectively, for different values of 𝜎=0.005, 𝜎=0.001, and 𝜎=0.0005. The simulations were run for 100,000 iterations with step size 0.001, and the first 90,000 iterations were discarded to allow 𝐼(𝑡) to reach its recurrent level. The distribution is slightly negative skewed for 𝜎=0.005 and positive skewed for 𝜎=0.001 and 𝜎=0.0005. Additionally, the corresponding sample skewness coefficients are −0.00067, 0.0016, and 0.0017, separately. These simulation results illustrate that the distribution of 𝐼(𝑡) has reached the stationary distribution. When 𝜎=0.005, from (70), the mean and variance of the stationary distribution are 75.9710 and 2.2016, respectively, while the sample mean and variance are 76.3754 and 0.3209. For 𝜎=0.001, the mean and variance of the stationary distribution are 75.9988 and 0.0876, compared to the sample mean and variance which are 75.9285 and 0.0157, separately. When 𝜎 reduces to 0.0005, the mean and variance of the stationary distribution are 75.9997 and 0.0219 and the corresponding sample mean and variance are 76.0248 and 0.0050. As 𝜎 decreases, 𝐼(𝑡) becomes more symmetric about 𝐼^∗=76 and the perturbations become much weaker. What is more, the normal quantile-quantile plots in Figure 5 suggest that these data are not far from being normally distributed.

Figure 4 

Histograms of the solution 𝐼(𝑡) for model (5) for 𝐼(0)=10 and differing values of 𝜎=0.005, 𝜎=0.001, and 𝜎=0.0005.

Figure 5 

Normal quantile-quantile plots of the solution 𝐼(𝑡) for model (5) for 𝐼(0)=10 and differing values of 𝜎=0.005, 𝜎=0.001, and 𝜎=0.0005.

From (70), we know that when 𝜎=0, the mean value 𝑚 equals the value of the endemic equilibrium 𝐼^∗ and the variance 𝑣 is zero. Furthermore, if 𝜎 increases, the mean value 𝑚 decreases and the variance 𝑣 increases as shown in Figure 6. Obviously, the greater the intensity of the noise is, the larger the variance is, that is, the stronger the fluctuations of the distribution are. And this is in accordance with the actual case.

(a)

(b)

Figure 6 

Plot the mean 𝑚 and variance 𝑣 against 𝜎. The endemic equilibrium 𝐼∗ of the deterministic model (1) is 76.

6. Conclusion

In this paper, we established a stochastic SIS epidemic model with saturated incidence to investigate the effect of environment noise. After proving the existence and uniqueness of the positive solution, we considered two kinds of stabilities: almost sure exponential stability and 𝑝th moment exponential stability. Our deduction shows that the 𝑝th moment exponential stability implies the almost sure exponential stability. As for the extinction of the disease, our results indicate that when the noise 𝜎² is larger than 𝛽(1+ℎ𝑁)²/𝑁∨𝛽²(1+ℎ𝑁)²/2(𝜇+𝛾), the disease will die out a.s. If the noise is small, the disease will also die out a.s. when 𝑅^𝑠₀₁<1, instead, the disease will be persistent in the mean when 𝑅^𝑠₀₂>1. Our investigation shows that the conditions for the disease to become extinct are much weaker than those in the corresponding deterministic model; that is, even if the disease dies out for the stochastic model, it may still persist for the corresponding deterministic model. Additionally, the half saturation constant influences the threshold value 𝑅^𝑠₀₁, which is not obtained in the deterministic model. Increasing the saturation effect by changing the behavior of the susceptible individuals and taking effective control measures of the government can reduce 𝑅^𝑠₀₁ and thus can reduce the spread of the disease. Finally, we proved the existence of a stationary distribution and derived the expressions for the mean and variance when the noise is small and 𝑅^𝑠₀₂>1. Our theoretical results were further verified by computer simulations.

Appendix

First, we give some theory in stochastic differential equations.

Let (Ω,ℱ,{ℱ_𝑡}_𝑡≥0,𝑃) be a complete probability space with {ℱ_𝑡}_𝑡≥0 satisfying the usual conditions (i.e., it is increasing and right continuous while {ℱ₀} contains all 𝑃-null sets). We use 𝑎∨𝑏 to denote max(𝑎,𝑏), 𝑎∧𝑏 to denote min(𝑎,𝑏), and a.s. to denote almost surely.

Definition A.1 (see [22]). Consider the 𝑑-dimensional stochastic differential equation 𝑑𝑥(𝑡)=𝑓(𝑥(𝑡),𝑡)𝑑𝑡+𝑔(𝑥(𝑡),𝑡)𝑑𝐵(𝑡),for𝑡≥𝑡₀, (A.1) with initial value 𝑥(𝑡₀)=𝑥₀∈𝑅^𝑑 and 𝐵(𝑡) being the 𝑚-dimensional standard Brownian motion. The solution is represented as 𝑥(𝑡;𝑡₀,𝑥₀). Assume that, for any 𝑡≥𝑡₀, there is 𝑓(0,𝑡)=0,𝑔(0,𝑡)=0. (A.2) So (A.1) has the solution 𝑥(𝑡)≡0. This solution is called the trivial solution.

Definition A.2 (see [22]). The trivial solution of (A.1) is said to be almost surely exponentially stable if limsup_𝑡→∞1𝑡log𝑥𝑡;𝑡₀,𝑥₀<0a.s., (A.3) for all 𝑥₀∈𝑅^𝑑.

Definition A.3 (see [22]). The trivial solution of (A.1) is said to be 𝑝th moment exponentially stable if there is a pair of positive constants 𝜆 and 𝐶 such that 𝐸𝑥𝑡;𝑡₀,𝑥₀^𝑝≤𝐶𝑥₀^𝑝𝑒^{−𝜆(𝑡−𝑡0)}on𝑡≥𝑡₀, (A.4) for all 𝑥₀∈𝑅^𝑑. When 𝑝=2, it is said to be exponentially stable in mean square.

Lemma A.4 (see [22]). Assume that there is a positive constant 𝐾 such that 𝑥^𝑇𝑓(𝑥,𝑡)∨𝑔(𝑥,𝑡)²≤𝐾|𝑥|²∀(𝑥,𝑡)∈𝑅^𝑑×𝑡₀,∞. (A.5) Then, the 𝑝th moment exponential stability of the trivial solution of (A.1) implies the almost sure exponential stability.

Next, we give some theory about stationary distributions [26].

Let 𝑋(𝑡) be a homogeneous Markov process in 𝐸_𝑙 (𝐸_𝑙 denotes Euclidean 𝑙-space) described by 𝑑𝑋(𝑡)=𝑓(𝑋)𝑑𝑡+^𝑚_𝑠=1𝑔_𝑠(𝑋)𝑑𝐵_𝑠(𝑡). (A.6) The diffusion matrix is 𝐴(𝑥)=(𝑎_𝑖𝑗(𝑥)), 𝑎_𝑖𝑗(𝑥)=∑^𝑚_𝑠=1𝑔^𝑖_𝑠(𝑥)𝑔^𝑗_𝑠(𝑥).

Assumption B. There exists a bounded domain 𝑈⊂𝐸_𝑙 with regular boundary Γ, having the following properties.(B.1)In the domain 𝑈 and some neighborhood 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 [26]). If (B.1) and (B.2) hold, then the Markov process 𝑋(𝑡) has a stationary distribution.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

The work is financially supported by the National Natural Science Foundation of China (nos. 11072182 and 11372233).