Extinction and Ergodic Property of Stochastic SIS Epidemic Model with Nonlinear Incidence Rate

Abstract

We investigate a stochastic SIS model with nonlinear incidence rate. We show that there exists a unique nonnegative solution to the system, and condition for the infectious individuals 𝐼(𝑡) to be extinct is given. Moreover, we prove that the system has ergodic property. Finally, computer simulations are carried out to verify our results.

1. Introduction

More attention has been paid to the epidemics models in order to monitor and curb the spread of some human diseases. A classical model is proposed by Kermack and McKendrick in 1927 [1]. They divided the population into three classes denoted by 𝑆(𝑡), 𝐼(𝑡), and 𝑅(𝑡), which expressed the number of susceptible individuals, infective individuals, and removed individuals at time 𝑡, respectively. The model is called susceptible-infected-removed (SIR) model, and SIR models were investigated by many researchers [2–4].

However some diseases, such as some sexually transmitted and bacterial diseases, do not have permanent immunity. In [5], they introduced a SIS model to describe the spread of the disease, which takes the following form: (𝑡)=(𝑡)𝐼(𝑡)(𝑡)(𝑡(𝑡)=(𝑡)𝐼(𝑡)−𝐼(𝑡𝑑𝑆𝜇𝑁−𝛽𝑆+𝛾𝐼−𝜇𝑆)𝑑𝑡,𝑑𝐼𝛽𝑆𝜇+𝛾)𝑑𝑡 (1) with initial values 𝑆₀+𝐼₀=𝑁. 𝑁 is the total size of the population. Where 𝑆(𝑡) and 𝐼(𝑡) express the number of susceptible individuals and infective individuals at time 𝑡, respectively, 𝜇 is the per capita birth (and death) rate, 𝛾 is the rate at which infected individuals become cured, and 𝛽 is the per capita contact rate. In model (1), they assumed that a rate of contacts by an infective individual with a susceptible individual is proportional to population size, and model (1) has been well studied [6].

In fact, population dynamics is inevitably affected by environmental white noise, which is always present. Since the parameters in the deterministic models are constant, they have some limitations when we describe the epidemics systems. Some researchers have paid their attention to the stochastic epidemics model [7–9]. Especially, in [10], Gray et al. consider that the parameter 𝛽 in (1) is perturbed with 𝛽⟶𝛽+𝜎𝐵̇(𝑡), (2) where 𝐵(𝑡) is Brownian motions and 𝜎 represents the intensities of the white noise. Corresponding to the deterministic model system (1), the stochastic system takes the following form: (𝑡)=(𝑡)𝐼(𝑡)(𝑡)(𝑡(𝑡)𝐼(𝑡)(𝑡),(𝑡)=(𝑡)𝐼(𝑡)−𝐼(𝑡(𝑡)𝐼(𝑡)(𝑡).𝑑𝑆𝜇𝑁−𝛽𝑆+𝛾𝐼−𝜇𝑆)𝑑𝑡−𝜎𝑆𝑑𝐵𝑑𝐼𝛽𝑆𝜇+𝛾)𝑑𝑡+𝜎𝑆𝑑𝐵 (3) Since 𝑆(𝑡)+𝐼(𝑡)=𝑁, then (3) is reduced to (𝑡)=𝛽((𝑡𝐼(𝑡)−𝐼(𝑡((𝑡𝐼(𝑡)(𝑡).𝑑𝐼𝑁−𝐼))𝜇+𝛾)𝑑𝑡+𝜎𝑁−𝐼))𝑑𝐵 (4) For model (4), they pointed out that(i)if (𝛽𝑁/(𝜇+𝛾)−𝜎²𝑁²/2(𝜇+𝛾))<1 and 𝜎²≤𝛽/𝑁, the disease 𝐼(𝑡) will die out with probability one;(ii)if (𝛽𝑁/(𝜇+𝛾)−𝜎²𝑁²/2(𝜇+𝛾))>1, then model (4) has a unique stationary distribution.

The incidence rate in (4) is bilinear, and several authors pointed out that the disease transmission process may have a nonlinear incidence rate [11, 12]. In [13], Xiao and Ruan propose an incidence rate 𝑔(𝐼)𝑆=𝛽𝐼𝑆1+𝛼𝐼², (5) where 𝛼 is the parameter that measures the psychological or inhibitory effect, and 1/(1+𝛼𝐼²) describes the psychological or inhibitory effect from the behavioral change of the susceptible individuals when the number of infective individuals is very large. It can be used to explain some phenomena; for example, the outbreak of severe acute respiratory syndrome (SARS) had such psychological effects on the general public [14]: for a very large number of infective individuals, the infection force may decrease as the number of infective individuals increases. Some control measures and policies, such as border screening, mask wearing, quarantine, isolation, and so forth, can decrease the infection rate although the number of infective individuals was getting relatively larger. Equation (1) with nonlinear incidence rate (5) and the disturbed parameter 𝛽 (𝛽→𝛽+𝜎𝐵̇(𝑡)) can be written as follows: (𝑡)=(𝑡)𝐼(𝑡)𝑑𝑆𝜇𝑁−𝛽𝑆1+𝛼𝐼²(𝑡)(𝑡)(𝑡−(𝑡)𝐼(𝑡)+𝛾𝐼−𝜇𝑆)𝑑𝑡𝜎𝑆1+𝛼𝐼²(𝑡)(𝑡),(𝑡)=(𝑡)𝐼(𝑡)𝑑𝐵𝑑𝐼𝛽𝑆1+𝛼𝐼²(𝑡)−𝐼(𝑡(𝑡)𝐼(𝑡)𝜇+𝛾)𝑑𝑡+𝜎𝑆1+𝛼𝐼²(𝑡)𝑑𝐵(𝑡). (6) The parameters appearing in (6) have the same meaning as those above. Given that 𝑆(𝑡)+𝐼(𝑡)=𝑁, it is sufficient to study the SDE for 𝐼(𝑡), (𝑡)=𝛽((𝑡𝐼(𝑡)𝑑𝐼𝑁−𝐼))1+𝛼𝐼²(𝑡)−𝐼(𝑡+𝜎((𝑡𝐼(𝑡)𝜇+𝛾)𝑑𝑡𝑁−𝐼))1+𝛼𝐼²(𝑡)𝑑𝐵(𝑡), (7) with initial value 𝐼(0)=𝐼₀∈(0,𝑁). Notice that when 𝛼=0, system (7) becomes (4). In this paper, we will analyze the dynamical behaviors of (7).

The organization of this paper is as follows. In the next section, we show that there exists a unique positive solution to (7). In Section 3, we carry out a qualitative analysis of the model (7) and extinction conditions for 𝐼(𝑡) is derived. We prove that the system has ergodic property under some condition in Section 4. In Section 5, we present some numerical simulations to illustrate our mathematical findings. A brief conclusion is given in Section 6.

Throughout this paper, let (Ω,ℱ,𝑃) be a complete probability space with a filtration {ℱ_𝑡}_𝑡≥0 satisfying the usual conditions (i.e., it is increasing and right continuous while ℱ₀ contains all 𝑃-null sets) and let 𝐵(𝑡) be a scalar Brownian motion defined on the probability space.

2. Existence and Uniqueness of the Global Positive Solution

In order for the model to make sense, we need to show the solution is global and nonnegative. However, theorem of existence and uniqueness (cf. Arnold [15] and Mao [16]) is not satisfied in (7). By using tools established by Mao et al. [17], we will show existence and uniqueness of the global positive solution of (7).

Theorem 1. For any given initial data 𝐼(0)=𝐼₀∈(0,𝑁), there exists a unique solution 𝐼(𝑡)∈(0,𝑁) for all 𝑡≥0 with probability 1.

Proof. It is obvious that the coefficients of the SDE (7) are locally continuous. For any given initial data 𝐼(0)=𝐼₀∈(0,𝑁), there exists a unique maximal local solution 𝐼(𝑡) on 𝑡∈[0,𝜏_𝑒), where 𝜏_𝑒 is the explosion time. In order to show that the solution is global, it is sufficient to show 𝜏_𝑒=∞ a.s. Let 𝑚₀>0 be sufficiently large so that 𝐼₀ lies within the interval [1/𝑚₀,𝑁−1/𝑚₀]. For each integer 𝑚≥𝑚₀, we define the stopping time 𝜏_𝑚=inf𝑡∈0,𝜏_𝑒(𝑡)∉1∶𝐼𝑚1,𝑁−𝑚, (8) where, throughout this paper, we set inf∅=∞ (as usual ∅ denotes the empty set). It is clear that 𝜏_𝑚 is increasing as 𝑚→∞. Let 𝜏_∞=lim_𝑚→∞𝜏_𝑚, whence 𝜏_∞≤𝜏_𝑒. It is easy to show that 𝜏_∞=∞ a.s. implies 𝜏_𝑒=∞ a.s. and 𝐼(𝑡)∈(0,𝑁) a.s. for all 𝑡≥0. Therefore, to complete this proof, it is enough to show that 𝜏_∞=∞ a.s.
Define a function 𝑉∶(0,𝑁)→𝑅₊ as follows: 𝑉(𝑥)=1𝑥+1.𝑁−𝑥 (9) By Itô’s formula, we get (𝑥)=𝑥−1𝑑𝑉𝑥²+1(𝑁−𝑥)²𝛽(𝑁−𝑥)1+𝛼𝑥²+𝜎−𝜇−𝛾₂𝑥²(𝑁−𝑥)²1+𝛼𝑥²²1𝑥³+1(𝑁−𝑥)³+−1𝑑𝑡𝑥²+1(𝑁−𝑥)²𝜎𝑥(𝑁−𝑥)1+𝛼𝑥²(𝑡)(𝑥)−1𝑑𝐵∶=𝐿𝑉𝑑𝑡+𝑥²+1(𝑁−𝑥)²𝜎𝑥(𝑁−𝑥)1+𝛼𝑥²𝑑𝐵(𝑡), (10) where (𝑥)−1𝐿𝑉=𝑥𝑥²+1(𝑁−𝑥)²𝛽(𝑁−𝑥)1+𝛼𝑥²+𝜎−𝜇−𝛾₂𝑥²(𝑁−𝑥)²1+𝛼𝑥²²1𝑥³+1(𝑁−𝑥)³≤𝜇+𝛾𝑥+𝛽𝑁𝑁−𝑥+𝜎²𝑁²1𝑥+1≤∨𝑁−𝑥𝜇+𝛾𝛽𝑁+𝜎²𝑁²𝑉(𝑥). (11) By almost the same method in the proof of [10], the desired result will be obtained.

3. Extinction

In this section, we will point out the condition for 𝐼(𝑡) to be extinct. We firstly do some preparation work.

Consider the following stochastic equation: 𝑑𝑋(𝑡)=𝑏(𝑋(𝑡))𝑑𝑡+𝜎(𝑋(𝑡))𝑑𝐵(𝑡), (12) and assume that the coefficients 𝜎∶𝐽→𝑅, 𝑏∶𝐽→𝑅 satisfy(1)𝜎²(𝑥)>0, ∀𝑥∈𝐽,(2)∀𝑥∈𝐼, ∃𝜖>0 such that ∫^𝑥+𝜖_𝑥−𝜖((1+|𝑏(𝑦)|)/𝜎²(𝑦))𝑑𝑦<∞,

where 𝐽=(𝑙,𝑟); −∞≤𝑙<𝑟≤∞.

Lemma 2 (see [18]). Assume that (1) and (2) hold, and let 𝑋(𝑡) be a weak solution of (12) in 𝐽, with nonrandom initial condition 𝑋₀=𝑥∈𝐽. Let 𝑝 be given by 𝑝(𝑥)=^𝑥_𝑐𝑒^{−∫𝑣𝑐(2𝑏(𝑦)/𝜎2(𝑦))𝑑𝑦}𝑑𝑣,𝑐∈𝐽. (13) If 𝑝(𝑙+)>−∞, 𝑝(𝑟−)=∞, then 𝑃lim_𝑡→∞𝑋(𝑡)=𝑙=𝑃sup_𝑡≥0𝑋(𝑡)<𝑟=1. (14)

Theorem 3. If 𝑅^𝑠₀∶=2(𝛽𝑁−𝜇−𝛾)/𝜎²𝑁²<1, then for any initial data 𝐼(0)=𝐼₀∈(0,𝑁), the solution of SDE (7) has the following property: 𝑃lim_𝑡→∞𝐼(𝑡)=0=1; (15) that is, the disease dies out with probability one.

Proof. Applying Lemma 2 with 𝑏(𝑥)=𝛽(𝑁−𝑥)𝑥/(1+𝛼𝑥²)−(𝜇+𝛾)𝑥, 𝜎(𝑥)=𝜎(𝑁−𝑥)𝑥/(1+𝛼𝑥²) and 𝑐∈𝐽=(0,𝑁), we can compute 2𝑏(^𝑥_𝑐𝜏)𝜎²(𝜏)=2𝑑𝜏𝜎^2𝑥_𝑐𝛽(𝑁−𝜏)1+𝛼𝜏²−𝜇+𝛾1+𝛼𝜏²²𝜏(𝑁−𝜏)²=2𝑑𝜏𝜎²𝛽𝑁−𝜇+𝛾𝑁²−ln𝑥𝛽𝑁−𝜇+𝛾𝑁²𝑁+𝛼𝛽𝑁+2𝜇+𝛾+3𝛼𝜇+𝛾²()−𝑁ln𝑁−𝑥𝜇+𝛾/𝑁²𝛼+1²−𝛼𝑁−𝑥𝜇+𝛾²2𝑥²−2𝜇+𝛾𝛼²𝑁+𝛼𝛽𝑥+𝐶₀. (16) Clearly, conditions (1) and (2) are satisfied. By calculation, we have 𝑝(𝑥)=^𝑥_𝑐exp−^𝑠_𝑐2𝑏(𝜏)𝜎²(𝜏)𝑑𝜏𝑑𝑠=𝑒^−𝑐0𝑥_𝑐𝑠^{−2(𝛽𝑁−(𝜇+𝛾))/𝜎2𝑁2}×(𝑁−𝑠)^{2[(𝛽𝑁−(𝜇+𝛾))/𝑁2+𝛼(𝛽𝑁+2(𝜇+𝛾)+3𝛼(𝜇+𝛾)𝑁2)]/𝜎2}×𝑒^{2(𝜇+𝛾)(𝑁2𝛼+1)2/𝑁𝜎2(𝑁−𝑠)}×𝑒^{((𝜇+𝛾)𝛼2𝑠2+(4(𝜇+𝛾)𝛼2𝑁+2𝛼𝛽)𝑠)/𝜎2}𝑑𝑠. (17)
We see that 𝑝(𝑁−)=𝑒^−𝑐0𝑁_𝑐𝑠^{−2(𝛽𝑁−(𝜇+𝛾))/𝜎2𝑁2}×(𝑁−𝑠)^{2[(𝛽𝑁−(𝜇+𝛾))/𝑁2+𝛼(𝛽𝑁+2(𝜇+𝛾)+3𝛼(𝜇+𝛾)𝑁2)]/𝜎2}×𝑒^{2(𝜇+𝛾)(𝑁2𝛼+1)2/𝑁𝜎2(𝑁−𝑠)}×𝑒^{((𝜇+𝛾)𝛼2𝑠2+(4(𝜇+𝛾)𝛼2𝑁+2𝛼𝛽)𝑠)/𝜎2}𝑑𝑠. (18) Let 1/(𝑁−𝑠)=𝑡; we obtain that 𝑝(𝑁−)=𝑒^−𝑐0∞_{1/(𝑁−𝑐)}(𝑁𝑡−1)^{−2(𝛽𝑁−(𝜇+𝛾))/𝜎2𝑁2}𝑡^{2(𝛽𝑁−(𝜇+𝛾))/𝜎2𝑁2}×𝑡^{−2[(𝛽𝑁−(𝜇+𝛾))/𝑁2+𝛼(𝛽𝑁+2(𝜇+𝛾)+3𝛼(𝜇+𝛾)𝑁2)]/𝜎2}×𝑒^{2(𝜇+𝛾)(𝑁2𝛼+1)2𝑡/𝑁𝜎2}×𝑒^{((𝜇+𝛾)𝛼2(𝑁−1/𝑡)2+(4(𝜇+𝛾)𝛼2𝑁+2𝛼𝛽)(𝑁−1/𝑡))/𝜎2}𝑡⁻²𝑑𝑡=𝑒^−𝑐0∞_{1/(𝑁−𝑐)}(𝑁𝑡−1)^{−2(𝛽𝑁−(𝜇+𝛾))/𝜎2𝑁2}×𝑡^{−(2𝛼(𝛽𝑁+2(𝜇+𝛾)+3𝛼(𝜇+𝛾)𝑁2)/𝜎2)−2}×𝑒^{2(𝜇+𝛾)(𝑁2𝛼+1)2𝑡/𝑁𝜎2}×𝑒^{((𝜇+𝛾)𝛼2(𝑁−1/𝑡)2+(4(𝜇+𝛾)𝛼2𝑁+2𝛼𝛽)(𝑁−1/𝑡))/𝜎2}𝑑𝑡. (19) Clearly, 𝑝(𝑁−)=∞. (20) When 𝑅^𝑠₀<1, we have −𝑝(0+)=𝑒^−𝑐0𝑐₀𝑠^{−2(𝛽𝑁−(𝜇+𝛾))/𝜎2𝑁2}×(𝑁−𝑠)^{2[(𝛽𝑁−(𝜇+𝛾))/𝑁2+𝛼(𝛽𝑁+2(𝜇+𝛾)+3𝛼(𝜇+𝛾)𝑁2)]/𝜎2}×𝑒^{2(𝜇+𝛾)(𝑁2𝛼+1)2/𝑁𝜎2(𝑁−𝑠)}×𝑒^{((𝜇+𝛾)𝛼2𝑠2+(4(𝜇+𝛾)𝛼2𝑁+2𝛼𝛽)𝑠)/𝜎2}𝑑𝑠<∞; (21) that is, 𝑝(0+)>−∞. (22) It can be straightly shown by Lemma 2 that 𝑃lim_𝑡→∞𝐼(𝑡)=0=1. (23) The proof is therefore completed.

4. Ergodic Property

In this section, we show that the small perturbation forces the infective individuals to be ergodic.

Theorem 4. If 𝑅^𝑠₀∶=2(𝛽𝑁−𝜇−𝛾)/𝜎²𝑁²>1, then, for any initial data 𝐼(0)=𝐼₀∈(0,𝑁), the solution of SDE (7) is ergodic.

Proof. If 𝑅^𝑠₀>1, we get exp−^𝑐₀^𝑠_𝑐2𝑏(𝜏)𝜎²(𝜏)𝑑𝜏𝑑𝑠=𝑒^−𝑐0𝑐₀𝑠^{−2(𝛽𝑁−(𝜇+𝛾))/𝜎2𝑁2}×(𝑁−𝑠)^{2[(𝛽𝑁−(𝜇+𝛾))/𝑁2+𝛼(𝛽𝑁+2(𝜇+𝛾)+3𝛼(𝜇+𝛾)𝑁2)]/𝜎2}×𝑒^{2(𝜇+𝛾)(𝑁2𝛼+1)2/𝑁𝜎2(𝑁−𝑠)}×𝑒^{((𝜇+𝛾)𝛼2𝑠2+(4(𝜇+𝛾)𝛼2𝑁+2𝛼𝛽)𝑠)/𝜎2}𝑑𝑠=∞. (24) It follows from (20) that exp−^𝑁_𝑐^𝑠_𝑐2𝑏(𝜏)𝜎²(𝜏)𝑑𝜏𝑑𝑠=𝑝(𝑁−)=∞. (25) Besides, 𝑚=^𝑁₀1𝜎²(𝑠)exp^𝑠_𝑐2𝑏(𝜏)𝜎²(𝜏)=𝑒𝑑𝜏𝑑𝑠_𝑐0𝜎^2𝑁₀1+𝛼𝑠²²×𝑠^{(2(𝛽𝑁−(𝜇+𝛾))/𝜎2𝑁2)−2}×(𝑁−𝑠)^{−(2[(𝛽𝑁−(𝜇+𝛾))/𝑁2+𝛼(𝛽𝑁+2(𝜇+𝛾)+3𝛼𝜇𝑁2)]/𝜎2)−2}×𝑒^{−2(𝜇+𝛾)(𝑁2𝛼+1)2/𝑁𝜎2(𝑁−𝑠)}×𝑒^{−((𝜇+𝛾)𝛼2𝑠2+(4(𝜇+𝛾)𝛼2𝑁+2𝛼𝛽)𝑠)/𝜎2}𝑑𝑠. (26) Let 𝑡=𝑁/(𝑁−𝑥)−1; then 𝑒𝑚=^𝑐0𝜎^2∞₀1+𝛼𝑁𝑡𝑡+1²²×𝑁𝑡𝑡+1^{(2(𝛽𝑁−(𝜇+𝛾))/𝜎2𝑁2)−2}×(𝑡+1)^{(2[(𝛽𝑁−(𝜇+𝛾))/𝑁2+𝛼(𝛽𝑁+2(𝜇+𝛾)+3𝛼(𝜇+𝛾)𝑁2)]/𝜎2)+2}×𝑒^{−2(𝜇+𝛾)(𝑁2𝛼+1)2(𝑡+1)/𝑁2𝜎2}×𝑒^{−((𝜇+𝛾)𝛼2(𝑁𝑡/(𝑡+1))2+(4(𝜇+𝛾)𝛼2+2𝛼𝛽)(𝑁𝑡/(𝑡+1)))/𝜎2}×(𝑡+1)⁻²=𝑒𝑑𝑡_𝑐0𝜎²𝑁^{(2(𝛽𝑁−(𝜇+𝛾))/𝜎2𝑁2)−2}×𝑒^{−2(𝜇+𝛾)(𝑁2𝛼+1)2/𝑁2𝜎2}×^∞₀𝛼𝑁²𝑡²+(𝑡+1)²²𝑡^{(2(𝛽𝑁−(𝜇+𝛾))/𝜎2𝑁2)−2}×(𝑡+1)^{(2𝛼(𝛽𝑁+2(𝜇+𝛾)+3𝛼(𝜇+𝛾)𝑁2)/𝜎2)−2}×𝑒^{−2(𝜇+𝛾)(𝑁2𝛼+1)2𝑡/𝑁2𝜎2}×𝑒^{−(((𝜇+𝛾)𝛼2(𝑁𝑡/(𝑡+1))2+(4(𝜇+𝛾)𝛼2+2𝛼𝛽)(𝑁𝑡/(𝑡+1)))/𝜎2)}𝑑𝑡. (27) When 𝑅^𝑠₀>1, that is, 2−(2(𝛽𝑁−(𝜇+𝛾))/𝜎²𝑁²)<1, we have 𝑚=^𝑁₀1𝜎²(𝑠)exp^𝑠_𝑐2𝑏(𝜏)𝜎²(𝜏)𝑑𝜏𝑑𝑠<∞. (28) The conditions of Theorem 1.16 in Kutoyants [19] follow from (24), (25), and (28). Therefore 𝐼(𝑡) is ergodic, and the invariant density is given by 𝜋(𝑥)=𝐶1+𝛼𝑥²²×𝑥^{(2(𝛽𝑁−(𝜇+𝛾))/𝜎2𝑁2)−2}×(𝑁−𝑥)^{−(2[(𝛽𝑁−(𝜇+𝛾))/𝑁2+𝛼(𝛽𝑁+2(𝜇+𝛾)+3𝛼(𝜇+𝛾)𝑁2)]/𝜎2)−2}×𝑒^{−2(𝜇+𝛾)(𝑁2𝛼+1)2/𝑁𝜎2(𝑁−𝑥)}×𝑒^{−((𝜇+𝛾)𝛼2𝑥2+(4(𝜇+𝛾)𝛼2𝑁+2𝛼𝛽)𝑥)/𝜎2},𝑥∈(0,𝑁), (29) where 𝐶 is a constant.

Remark 5. If 𝛼=0, then 𝜋(𝑥)=𝐶𝑥^{(2[𝛽𝑁−(𝜇+𝛾)]/𝜎2𝑁2)−2}×(𝑁−𝑥)^{(−2[𝛽𝑁−(𝜇+𝛾)]/𝜎2𝑁2)−2}×𝑒^{−2(𝜇+𝛾)/𝑁𝜎2(𝑁−𝑥)},𝑥∈(0,𝑁). (30) By the property of density function, we have 1=^𝑁₀𝐶𝑥^{(2[𝛽𝑁−(𝜇+𝛾)]/𝜎2𝑁2)−2}×(𝑁−𝑥)^{(−2[𝛽𝑁−(𝜇+𝛾)]/𝜎2𝑁2)−2}×𝑒^{−2(𝜇+𝛾)/𝑁𝜎2(𝑁−𝑥)}𝑑𝑥. (31) Let 𝑥=𝑁𝑡/(𝑡+1); then ^𝑁₀𝐶𝑥^{(2[𝛽𝑁−(𝜇+𝛾)]/(𝜎2𝑁2))−2}×(𝑁−𝑥)^{(−2[𝛽𝑁−(𝜇+𝛾)]/𝜎2𝑁2)−2}𝑒^{−2(𝜇+𝛾)/𝑁𝜎2(𝑁−𝑥)}𝑑𝑥=𝐶𝑁⁻³𝑒^{−2(𝜇+𝛾)/𝜎2𝑁2}×^∞₀𝑡^{2[𝛽𝑁−(𝜇+𝛾)]/𝜎2𝑁2}(𝑡+1)²𝑒^{−2(𝜇+𝛾)𝑡/𝜎2𝑁2}𝑑𝑡=𝐶𝑁⁻³𝑒^{−2(𝜇+𝛾)/𝜎2𝑁2}𝑁²𝜎²2𝜇+𝛾^{2[𝛽𝑁−(𝜇+𝛾)]/𝜎2𝑁2}×𝛽𝑁−𝜇−𝛾2𝛽𝑁−𝜇−𝛾−𝜎²𝑁²𝜇+𝛾𝜎²𝑁²+4𝛽𝑁−𝜇+𝛾−2𝜎²𝑁²𝜎²𝑁²+2𝜇+𝛾𝜎²𝑁²2×Γ𝛽𝑁−𝜇−𝛾𝜎²𝑁²−1, (32) which implies 1𝐶=𝑁⁻³𝑒^{−2(𝜇+𝛾)/𝜎2𝑁2}×𝑁²𝜎²2𝜇+𝛾^{2[𝛽𝑁−(𝜇+𝛾)]/𝜎2𝑁2}×𝛽𝑁−𝜇−𝛾2𝛽𝑁−𝜇−𝛾−𝜎²𝑁²𝜇+𝛾𝜎²𝑁²+4𝛽𝑁−𝜇−𝛾−2𝜎²𝑁²𝜎²𝑁²+2𝜇+𝛾𝜎²𝑁²2×Γ𝛽𝑁−𝜇−𝛾𝜎²𝑁²−1. (33) Next, we compute 𝐸(𝑋): 𝐸(𝑋)=^𝑁₀𝐶𝑥^{(2[𝛽𝑁−(𝜇+𝛾)]/𝜎2𝑁2)−1}×(𝑁−𝑥)^{(−2[𝛽𝑁−(𝜇+𝛾)]/𝜎2𝑁2)−2}𝑒^{−2(𝜇+𝛾)/𝑁𝜎2(𝑁−𝑥)}𝑑𝑥=𝐶𝑁⁻²𝑒^{−2(𝜇+𝛾)/𝜎2𝑁2}𝑁²𝜎²2𝜇+𝛾^{2[𝛽𝑁−(𝜇+𝛾)]/𝜎2𝑁2}×2𝛽𝑁−𝜇−𝛾−𝜎²𝑁²𝜎²𝑁²𝛽𝑁2𝜇+𝛾×Γ𝛽𝑁−𝜇−𝛾𝜎²𝑁²−1. (34) This together with (33) shows that 𝐸(𝑋)=𝛽2𝛽𝑁−𝜇−𝛾−𝜎²𝑁²2𝛽²−𝜇+𝛾+𝛽𝑁𝜎². (35) Similarly, we have 𝐸𝑋²=𝛽𝑁−𝜇−𝛾2𝛽𝑁−𝜇−𝛾−𝜎²𝑁²2𝛽²−𝜇+𝛾+𝛽𝑁𝜎²=𝛽𝑁−𝜇−𝛾𝐸(𝑋)𝛽. (36) Consequently, (𝑋)=𝐸(𝑋)Var𝛽𝑁−𝜇−𝛾𝛽−[𝐸(𝑋)]². (37) When 𝛼=0, model (7) becomes (4), and the mean and variance of the stationary distribution of model (4) are the same as the results in [10].

5. Simulations

We utilize the method developed in [20] to illustrate our findings. We have the corresponding discretization equation: 𝐼_𝑘+1=𝐼_𝑘+𝐼_𝑘𝛽𝑁−𝐼_𝑘1+𝛼𝐼²_𝑘−𝜇+𝛾𝐼_𝑘Δ𝑡+𝜎𝐼_𝑘𝑁−𝐼_𝑘1+𝛼𝐼²_𝑘𝜀_𝑘√Δ𝑡, (38) where 𝜀_𝑘 are the Gaussian random variables 𝑁(0,1).

Using the discretized equation and with the help of Matlab software, choosing the appropriate parameters 𝑁=2, 𝛽=0.5, 𝜇=0.2, 𝛾=0.3, and 𝛼=0.5, we get simulations of system (7) and the corresponding deterministic system.

In Figure 1, we choose initial values 𝐼(0)=1.8 and 0.4, respectively, 𝜎=0.6, noting that 2(𝛽𝑁−𝜇−𝛾)/𝜎²𝑁²=0.694<1, as the case in Theorem 3 expected, for any initial value 𝐼₀∈(0,2), the large white noise leads to the extinction of 𝐼(𝑡) and the solution of system (7) tends to zero; that is, the disease dies out. While the solution of corresponding deterministic system does not tend to zero.

(a)

(b)

Figure 1 

The solution 𝐼(𝑡) of SIS model (7) compared to the deterministic system with different initial value 𝐼(0)=1.8 and 0.4, 𝜎=0.6. The blue line represents the solution of stochastic system, and the black line represents the solution of the corresponding deterministic system.

In Figure 2, we choose 𝜎=0.2 and 0.05, respectively, with initial values 𝐼(0)=1.5, which satisfy the cases in Theorem 4; that is, 2(𝛽𝑁−𝜇−𝛾)/𝜎²𝑁²>1. From the left picture in Figure 2, we can see that the solution of system (7) is fluctuating in a small neighborhood, and there is a stationary distribution (see the histogram in Figures 2(b) and 2(d)); the disease becomes endemic.

(a)

(b)

(c)

(d)

Figure 2 

The solution 𝐼(𝑡) of SIS model (7) and its histogram with differing value of 𝜎=0.2,0.05.

6. Conclusion

In this paper, we analyze the dynamic behaviors of a stochastic SIS model with nonlinear incidence rate, under the assumption that the population lives in an environment subjected to random fluctuations which mainly affect the disease transmission term. First of all, we show that there exists a unique positive solution in system (7). Moreover, we obtain the threshold between prevalence and extinction of the disease; that is, if 𝑅^𝑠₀=2(𝛽𝑁−𝜇−𝛾)/𝜎²𝑁²<1, the disease will die out with probability one; if 𝑅^𝑠₀=2(𝛽𝑁−𝜇−𝛾)/𝜎²𝑁²>1, 𝐼(𝑡) is ergodic, which means the disease will become endemic. Finally, we illustrate our results with computer simulations.

We see that if 𝛼=0, (7) becomes (4), which is studied in [10]; they show that if 𝛽𝑁/(𝜇+𝛾)−𝜎²𝑁²/2(𝜇+𝛾)<1 and 𝜎²≤𝛽/𝑁, the disease 𝐼(𝑡) will die out with probability one. Obviously, 𝛽𝑁/(𝜇+𝛾)−𝜎²𝑁²/2(𝜇+𝛾)<1 and 𝑅^𝑠₀=2(𝛽𝑁−𝜇−𝛾)/𝜎²𝑁²<1 are equivalent. The condition 𝜎²≤𝛽/𝑁, however, is not necessary in our investigation. In addition, if 𝑅^𝑠₀=2(𝛽𝑁−𝜇−𝛾)/𝜎²𝑁²>1, the invariant density is obtained, which is not mentioned in [10].

An extension of our work is to consider a stochastic SIS model with the general incidence rate 𝛽𝐼^𝑝𝑆/(1+𝛼𝐼^𝑞), and it is currently a work in progress.

Acknowledgments

The work was supported by the Program for Changjiang Scholars and Innovative Research Team in University, NSF of China (nos. 11371085 and 11101072); the Ph.D. Programs Foundation of Ministry of China (no. 200918); and the Nature Science Foundation of Changchun Normal University (no. 2010007).