Stochastic Dynamics of an SIRS Epidemic Model with Ratio-Dependent Incidence Rate

Abstract

We investigate the complex dynamics of an epidemic model with nonlinear incidence rate of saturated mass action which depends on the ratio of the number of infectious individuals to that of susceptible individuals. We first deal with the boundedness, dissipation, persistence, and the stability of the disease-free and endemic points of the deterministic model. And then we prove the existence and uniqueness of the global positive solutions, stochastic boundedness, and permanence for the stochastic epidemic model. Furthermore, we perform some numerical examples to validate the analytical findings. Needless to say, both deterministic and stochastic epidemic models have their important roles.

1. Introduction

Since the pioneer work of Kermack and McKendrick [1], mathematical models are used extensively in analyzing the spread, and control of infectious diseases qualitatively and quantitatively. The research results are helpful for predicting the developing tendencies of the infectious disease, for determining the key factors of the disease spreading, and for seeking the optimum strategies for preventing and controlling the spread of infectious diseases [2]. And in modeling communicable diseases, the incidence function has been considered to play a key role in ensuring that the models indeed give reasonable qualitative description of the transmission dynamics of the diseases [3–7].

Let 𝑆(𝑡) be the number of susceptible individuals, 𝐼(𝑡) the number of infective individuals, and 𝑅(𝑡) the number of removed individuals at time 𝑡, respectively. We consider the general SIRS epidemic model: 𝑑𝑆()𝑑𝑡=𝑏−𝑑𝑆−𝐻𝐼,𝑆+𝛾𝑅,𝑑𝐼()−𝑑𝑡=𝐻𝐼,𝑆𝑑+𝜇+𝛿𝐼,𝑑𝑅𝑑𝑡=𝜇𝐼−𝑑+𝛾𝑅, (1) where 𝑏 is the recruitment rate of the population, 𝑑 is the natural death rate of the population, 𝜇 is the natural recovery rate of the infective individuals, 𝛾 is the rate at which recovered individuals lose immunity and return to the susceptible class, and 𝛿 is the disease-induced death rate. And the transmission of the infection is governed by an incidence rate 𝐻(𝐼,𝑆).

In [8], Liu et al. proposed the general saturated nonlinear incidence rate: 𝐻(𝐼,𝑆)=𝑆𝑔(𝐼),𝑔(𝐼)=𝑘𝐼^𝑙1+𝛼𝐼^ℎ, (2) where the parameters 𝑙 and ℎ are positive constants, 𝑘 the proportionality constant, and 𝛼 is a nonnegative constant, which measures the psychological or inhibitory effect. 𝑘𝐼^𝑙 measures the infection force of the disease, and 1/(1+𝛼𝐼^ℎ) measures the inhibition effect from the behavioral change of the susceptible individuals when their number increases or from the crowding effect of the infective individuals. And the other nonlinear incidence rates are considered in [6, 9–19].

Note that the infectious force 𝑔(𝐼) of classical disease transmission models typically is only a function of infective individuals. But in the transmission of communicable diseases, it involves both infective individuals and susceptible individuals. Thus, Yuan et al. [18, 19] studied the infections force function with a ratio-dependent nonlinear incident rate which takes the following form: 𝑔𝐼𝑆=𝑘(𝐼/𝑆)^𝑙1+𝛼(𝐼/𝑆)^ℎ. (3) And in [19], Li et al. focus on an epidemic disease of SIRS type, in which they assume that the infectious force takes the form of (3) with 𝑙=1 and ℎ=1, and the model is as follows: 𝑑𝑆𝑑𝑡=𝑏−𝑑𝑆−𝑘𝐼𝑆𝑆+𝛼𝐼+𝛾𝑅,𝑑𝐼=𝑑𝑡𝑘𝐼𝑆−𝑆+𝛼𝐼𝑑+𝜇𝐼,𝑑𝑅𝑑𝑡=𝜇𝐼−𝑑+𝛾𝑅, (4) where all the parameters are nonnegative and have the same definitions as in model (1).

From the standpoint of epidemiology, we are only interested in the dynamics of model (4) in the closed first quadrant ℝ³₊={(𝑆,𝐼,𝑅)∶𝑆≥0,𝐼≥0,𝑅≥0}. Thus, we consider only the epidemiological meaningful initial conditions 𝑆(0)>0, 𝐼(0)>0, 𝑅(0)>0. Straightforward computation shows that model (4) is continuous and Lipschizian in ℝ³₊ if we redefine that when (𝑆,𝐼,𝑅)=(0,0,0), 𝑑𝑆/𝑑𝑡=𝑏, 𝑑𝐼/𝑑𝑡=0, 𝑑𝑅/𝑑𝑡=0. Hence, the solution of model (4) with positive initial conditions exists and is unique.

It is clear that the limit set of model (4) is on the plane 𝑆+𝐼+𝑅=𝑏/𝑑, and the model can be reduced to the following: 𝑑𝑆=𝑑𝑡𝑏+𝛾𝑏𝑑−𝑑+𝛾𝑆−𝛾𝐼−𝑘𝐼𝑆,𝑆+𝛼𝐼𝑑𝐼=𝑑𝑡𝑘𝐼𝑆−𝑆+𝛼𝐼𝑑+𝜇𝐼, (5) when (𝑆,𝐼)=(0,0), 𝑑𝑆/𝑑𝑡=𝑏+(𝛾𝑏/𝑑), 𝑑𝐼/𝑑𝑡=0.

For mathematical simplicity, let us nondimensionalize model (5) as in [19] with the following scaling: 𝑑𝑥=𝑑+𝜇𝑏𝑑+𝛾𝑆,𝑦=𝑑𝛾𝑏𝑑+𝛾𝐼,𝜏=𝑑+𝜇𝑡. (6) We still use variable 𝑡 instead of 𝜏, and model (5) takes the following form: 𝑑𝑥𝑑𝑡=1−𝑞𝑥−𝑦−𝑎𝑥𝑦,𝑥+𝑝𝑦𝑑𝑦=𝑅𝑑𝑡⁰𝑥𝑥+𝑝𝑦−1𝑦, (7) where 𝑞=(𝑑+𝛾)/(𝑑+𝜇), 𝑝=𝛼(𝑑+𝜇)/𝛾, 𝑎=𝑘/𝛾 are positive constants. 𝑅₀=𝑘/(𝑑+𝜇) is the basic reproduction number. And when (𝑆,𝐼)=(0,0), 𝑑𝑥/𝑑𝑡=1, 𝑑𝑦/𝑑𝑡=0.

On the other hand, if the environment is randomly varying, the population is subject to a continuous spectrum of disturbances [20, 21]. That is to say, population systems are often subject to environmental noise; that is, due to environmental fluctuations, parameters involved in epidemic models are not absolute constants, and they may fluctuate around some average values. Based on these factors, more and more people began to be concerned about stochastic epidemic models describing the randomness and stochasticity [22–34], and the stochastic epidemic models can provide an additional degree of realism if compared to their deterministic counterparts [10, 35–47]. In Particular, Mao et al. [26] obtained the interesting and surprising conclusion: even a sufficiently small noise can suppress explosions in population dynamics. Beretta et al. [35] obtained the stability of epidemic model with stochastic time delays influenced by probability under certain conditions. Carletti [36] studied the stable properties of a stochastic model for phage-bacteria interaction in open marine environment analytically and numerically. In [37], establishing some stochastic models and studying of several endemic infections with demography, Nåsell found that some deterministic models are unacceptable approximations of the stochastic models for a large range of realistic parameter values. Dalal et al. [39, 40] showed that stochastic models had nonnegative solutions and carried out analysis on the asymptotic stability of models. In [41], Yu et al. presented stochastic asymptotic stability of the epidemic point of the two-group SIR model with random perturbation. It is shown in [45] that the SIR model has a unique global positive and asymptotic solution. But to our knowledge, the research on the stochastic dynamics of the epidemic model with ratio-dependent nonlinear incidence rate seems rare.

There are different possible approaches to including random effects in the model, both from a biological and from a mathematical perspectives [48]. Our basic approach is analogous to that of Beddington and May [20], which is pursued in [48], and also, for example, in [45, 47] to epidemic models, in which they considered that the environmental noise was proportional to the variables. Following them, in this paper, we assume that stochastic perturbations are of a white noise type which is directly proportional to 𝑥(𝑡), 𝑦(𝑡), influenced on the 𝑑𝑥(𝑡)/𝑑𝑡 and 𝑑𝑦(𝑡)/𝑑𝑡 in model (4). In this way, we introduce stochastic perturbation terms into the growth equations of susceptible and infected individuals to incorporate the effect of randomly fluctuating environment, and the following stochastic differential equation is corresponding to model (7): 𝑑𝑥=1−𝑞𝑥−𝑦−𝑎𝑥𝑦𝑥+𝑝𝑦𝑑𝑡+𝜎₁𝑥𝑑𝐵₁(𝑡),𝑅𝑑𝑦=⁰𝑥𝑦𝑥+𝑝𝑦−𝑦𝑑𝑡+𝜎₂𝑦𝑑𝐵₂(𝑡), (8) where 𝜎₁, 𝜎₂ are real constants and known as the intensity of environmental fluctuations, and 𝐵₁(𝑡), 𝐵₂(𝑡) are independent standard Brownian motions.

The aim of this paper is to consider the dynamics of the epidemic models (7) and (8). The paper is organized as follows. In Section 2, we give some properties about deterministic model (7). In Section 3, we carry out the analysis of the dynamical properties of stochastic model (8). And in Section 4, we give some numerical examples and make a comparative analysis of the stability of the model with deterministic and stochastic environments and have some discussions.

2. Dynamics of the Deterministic Model

Let us begin to determine the location and number of the equilibria of model (7). It is easy to see that if 𝑅₀<1, the disease-free point 𝐸₀=(1/𝑞,0) is the unique equilibrium, corresponding to the extinction of the disease; if 𝑅₀>1, in addition to the disease-free point 𝐸₀, there is a unique endemic point 𝐸^∗=(𝑥^∗,𝑦^∗), corresponding to the survival of the disease, described by the following expressions: 𝑥^∗=𝑝𝑅₀𝑝𝑞𝑅₀+𝑅₀+𝑎𝑅₀−1,𝑦^∗=𝑅₀−1𝑝𝑥^∗. (9)

The Jacobian matrix of model (7) at 𝐸₀ is as follows: −𝑞−1−𝑎0𝑅₀−1. (10) It follows that 𝐸₀ is asymptotically stable if 𝑅₀<1 and unstable if 𝑅₀>1.

The Jacobian matrix of model (7) at 𝐸^∗ is as follows: 𝐽^∗=𝐽₁₁𝐽₁₂𝐽₂₁𝐽₂₂, (11) where 𝐽₁₁=−𝑝𝑞𝑅²₀+𝑎𝑅₀−1²𝑝𝑅²₀,𝐽₁₂𝑅=−²⁰+𝑎𝑅²₀,𝐽₂₁=𝑅₀−1²𝑝𝑅₀,𝐽₂₂𝑅=−⁰−1𝑅₀. (12) It is easy, by simple computations, to see that tr𝐽^∗=𝐽₁₁+𝐽₂₂𝐽<0,det^∗=𝑝𝑞𝑅₀+𝑎+𝑅₀𝑅₀−1𝑝𝑅²₀>0. (13)

Summarizing the above, we have the following results on the dynamics of model (7).

Theorem 1. (i) If 𝑅₀<1, then model (7) has a unique disease-free equilibrium 𝐸₀ which is asymptotically stable.
(ii) If 𝑅₀>1, then model (7) has two equilibria, a disease-free equilibrium 𝐸₀ which is an unstable saddle and an endemic equilibrium 𝐸^∗ which is asymptotically stable.

As a matter of fact, we can prove that the endemic point 𝐸^∗=(𝑥^∗,𝑦^∗) is also global asymptotically stable. For more details, see [19].

In Figure 1, we show the dynamics of the deterministic model (7) with the following parameters: 𝑎=0.3,𝑝=0.5,𝑞=2,𝑅₀=4.5. (14) In this case, 𝐸₀=(0.5,0) is a saddle point. 𝐸^∗=(0.10563,0.73944) is globally asymptotically stable.

Figure 1 

The dynamics of model (7). The parameters are taken as (14). 𝐸0=(0.5,0) is a saddle point. 𝐸∗=(0.11,0.74) is globally asymptotically stable.

In the following, we will focus on the boundedness, dissipation, and persistence of mode (7).

Theorem 2. All the solutions of model (7) with the positive initial condition (𝑥(0),𝑦(0)) are uniformly bounded within a region Γ, where Γ=𝑥,𝑦∈ℝ²₊𝑎∶𝑥+𝑅₀1𝑦≤min𝑞,𝑅₀𝑅₀.+𝑎 (15)

Proof. Define function 𝑁(𝑡)(𝑡)+𝑎=𝑥𝑅₀𝑦(𝑡). (16)
Differentiating 𝑁(𝑡) with respect to time 𝑡 along the solutions of model (7), we can get the following: 𝑑𝑁(𝑡)=𝑑𝑡𝑑𝑥+𝑎𝑑𝑡𝑅₀𝑑𝑦𝑎𝑑𝑡=1−𝑞𝑥−1+𝑅₀𝑦. (17)
Thus, we obtain the following: 𝑑𝑁(𝑡)(𝑡)𝑎𝑑𝑡+𝜂𝑁=1−𝑞−𝜂𝑥−1+𝑅₀−𝜂𝑦<1, (18) where 𝜂<min{𝑞,1+(𝑎/𝑅₀)}. And we obtain the following: ≤10<𝑁𝑥,𝑦𝜂+𝑁𝑥(0),𝑦(0)𝑒^−𝜂𝑡. (19) As 𝑡→∞, 0<𝑁≤1/𝜂. Therefore, all solutions of model (7) enter into the region Γ. This completes the proof.

Theorem 3. If 𝑅₀>1, model (7) is dissipative.

Proof. Since all solutions of model (7) are positive, by the first equation of (7), we have the following: 𝑑𝑥𝑑𝑡≤1−𝑞𝑥. (20) A standard comparison theorem shows that limsup_𝑡→∞𝑥(𝑡)≤1𝑞. (21)
Hence, for any 0<𝜀≪1 and large 𝑡, 𝑥≤(1/𝑞)+𝜀. It then follows that 𝑦 satisfies the following: 𝑑𝑦≤𝑦𝑅𝑑𝑡⁰−1/𝑞+𝜀𝑅₀−1−𝑝𝑦1/𝑞+𝜀+𝑝𝑦. (22) The arbitrariness of 𝜀 then implies that limsup_𝑡→∞𝑦(𝑡)≤𝑅₀−1.𝑝𝑞 (23)

Theorem 4. If 𝑅₀>1 and 𝑝𝑞<(1+𝑎)(𝑅₀−1), then model (7) is permanent; that is, there exists 𝜀>0 (independent of initial conditions), such that liminf_𝑡→∞𝑥(𝑡)>𝜀, liminf_𝑡→∞𝑦(𝑡)>𝜀.

Proof. By the first equation in (7), we have the following: 𝑑𝑥()𝑑𝑡=1−𝑞𝑥−1+𝑎𝑦+𝑎𝑝𝑦²()𝑥+𝑝𝑦>1−𝑞𝑥−1+𝑎𝑦. (24) If 𝑅₀>1 and 𝑝𝑞<(1+𝑎)(𝑅₀−1), from the proof of Theorem 3, we see that limsup_𝑡→∞𝑦(𝑡)≤(𝑅₀−1)/𝑝𝑞. Thus, for any 0<𝜀≤(𝑅₀−1)/𝑝𝑞 and large 𝑡, 𝑦(𝑡)>((𝑅₀−1)/𝑝𝑞)−𝜀. As a result, we have the following: 𝑑𝑥()𝑅𝑑𝑡>1−1+𝑎⁰−1𝑝𝑞−𝜀−𝑞𝑥. (25) With the comparison principle, the arbitrariness of 𝜀 implies that liminf_𝑡→∞𝑥(𝑡)≥𝑝𝑞−(1+𝑎)𝑅₀−1𝑝𝑞²≜𝑥. (26) Hence, for any 0<𝜀<(𝑝𝑞−(1+𝑎)(𝑅₀−1))/𝑝𝑞² and large 𝑡, 𝑥(𝑡)>𝑥−𝜀.
And for large 𝑡, we have the following: 𝑑𝑦>𝑦𝑥𝑑𝑡−𝜀𝑅₀−1−𝑝𝑦𝑥.−𝜀+𝑝𝑦 (27) Therefore, liminf_𝑡→∞𝑥(𝑡)≥𝑥−𝜀𝑅₀−1𝑝. (28) The arbitrariness of 𝜀 then implies that liminf_𝑡→∞𝑥(𝑡)≥𝑥𝑅₀−1𝑝≜𝑦. (29) Choosing a positive number 𝜖 such that 𝜖<min{𝑥/2,𝑦/2}, we see that liminf_𝑡→∞𝑥(𝑡)>𝜖,liminf_𝑡→∞𝑦(𝑡)>𝜀. (30) This ends the proof.

Noting that if the parameters of model (7) are fixed as (14), we can obtain the following: 𝑅₀>1,𝑝𝑞=0.6<(1+𝑎)𝑅₀−1=4.55, (31) and from Theorems 3 and 4, we can conclude that model (7) is dissipation and persistence.

3. Dynamics of the Stochastic Model

In this subsection, we investigate the dynamical behavior of the stochastic model (8). Throughout this paper, let (Ω,ℱ,𝒫) be a complete probability space with a filtration {ℱ_𝑡}_𝑡∈ℝ+ satisfying the usual conditions (i.e., it is right continuous and increasing while ℱ₀ contains all 𝒫-null sets). 𝐵₁(𝑡), 𝐵₂(𝑡) are the Brownian motions defined on this probability space. We denote by 𝑋(𝑡)=((𝑥(𝑡),𝑦(𝑡)) and |𝑋(𝑡)|=(𝑥²(𝑡)+𝑦²(𝑡))^1/2. Denote Λ={(𝑥,𝑦)∈ℝ²₊∶𝑥≥𝑎/𝑅₀,𝑦>0}.

Denote by 𝐶^2,1(ℝ^𝑑×(0,∞);ℝ₊) the family of all nonnegative functions 𝑉(𝑥,𝑡) defined on ℝ^𝑑×(0,∞) such that they are continuously twice differentiable in 𝑥 and once in 𝑡. Define the differential operator 𝐿 associated with 𝑑-dimensional stochastic differential equation: 𝑑𝑥(𝑡)=𝑓(𝑥(𝑡),𝑡)𝑑𝑡+ℎ(𝑥(𝑡),𝑡)𝑑𝐵(𝑡) (32) by 𝜕𝐿=+𝜕𝑡^𝑑_𝑖=1𝑓_𝑖()𝜕𝑥,𝑡𝜕𝑥_𝑖+12^𝑑_𝑖,𝑗=1ℎ^𝑇(𝑥,𝑡)ℎ(𝑥,𝑡)_𝑖𝑗𝜕²𝜕𝑥_𝑖𝜕𝑥_𝑗. (33) If 𝐿 acts in a function 𝑉∈𝐶^2,1(ℝ^𝑑×(0,∞);ℝ₊), then 𝐿𝑉(𝑥,𝑡)=𝑉_𝑡(𝑥,𝑡)+𝑉_𝑥()𝑓()+1𝑥,𝑡𝑥,𝑡2traceℎ^𝑇(𝑥,𝑡)𝑉_𝑥𝑥(𝑥,𝑡)ℎ(𝑥,𝑡), (34) where 𝑇 means transposition.

3.1. Existence and Uniqueness of Global Positive Solutions

To investigate the dynamical behavior of model (8), the first thing concerned is whether the solution is global existent. In this section, using the Lyapunov analysis method (mentioned in [24]), we will show the solution of model (8) is global and nonnegative.

Lemma 5. There is a unique local positive solution (𝑥(𝑡),𝑦(𝑡)) for 𝑡∈[0,𝜏_𝑒) to model (8) almost surely (a.s.) for the initial value (𝑥(0),𝑦(0))∈Λ, where 𝜏_𝑒 is the explosion time.

Proof. Set 𝑢(𝑡)=ln𝑥(𝑡),𝑣(𝑡)=ln𝑦(𝑡), (35) by Itô formula, we have the following: 1𝑑𝑢=𝑒^𝑢𝑒−𝑞−^𝑣𝑒^𝑢−𝑎𝑒^𝑣𝑒^𝑢+𝑝𝑒^𝑣−𝜎²₁2𝑑𝑡+𝜎₁𝑑𝐵₁(𝑡),𝑅𝑑𝑣=⁰𝑒^𝑢𝑒^𝑢+𝑝𝑒^𝑣𝜎−1−²²2𝑑𝑡+𝜎₂𝑑𝐵₂(𝑡), (36) at 𝑡≥0 with initial value 𝑢(0)=ln𝑥(0), 𝑣(0)=ln𝑦(0).
It is easy to see that the coefficients of model (36) satisfy the local Lipschitz condition, and there is a unique local solution 𝑢(𝑡), 𝑣(𝑡) on [0,𝜏_𝑒) [24]. Therefore, 𝑥(𝑡)=𝑒^𝑢(𝑡), 𝑦(𝑡)=𝑒^𝑣(𝑡) are the unique positive local solutions to model (36) with the initial value (𝑥(0),𝑦(0))∈Λ.

Lemma 5 only tells us that there exists a unique local positive solution to model (8). In the following, we show this solution is global; that is, 𝜏_𝑒=∞, which is motived by the work of Luo and Mao [29].

Theorem 6. Consider model (8), for any given initial value (𝑥(0),𝑦(0))∈Λ, there is a unique solution (𝑥(𝑡),𝑦(𝑡)) on 𝑡≥0 and the solution will remain in Λ with probability 1.

Proof. Let 𝑛₀>0 be sufficiently large for 𝑥(0) and 𝑦(0) lying within the interval [1/𝑛₀,𝑛₀]. For each integer 𝑛>𝑛₀, define the stopping times: 𝜏_𝑛=inf𝑡∈0,𝜏_𝑒(𝑡)∉1∶𝑥𝑛(𝑡)∉1,𝑛or𝑦𝑛,𝑛. (37) We set inf∅=∞  (∅ represents the empty set) in this paper. 𝜏_𝑛 is increasing as 𝑛→∞. Let 𝜏_∞=lim_𝑛→∞𝜏_𝑛; then 𝜏_∞≤𝜏_𝑒 a.s..
In the following, we need to show 𝜏_∞=∞ a.s. If this statement is violated, there exist constants 𝑇>0 and 𝜀∈(0,1) such that 𝒫{𝜏_∞≤𝑇}>𝜀. As a consequence, there exists an integer 𝑛₁≥𝑛₀ such that 𝒫𝜏_𝑛≤𝑇≥𝜀,𝑛≥𝑛₁. (38)
Define a function 𝑉₁∶Λ→ℝ₊ by the following: 𝑉₁=𝑅𝑥,𝑦⁰𝑎𝑅𝑥−1−ln⁰𝑎𝑥+𝑦−1−ln𝑦, (39) which is a non-negativity function.
If (𝑥(𝑡),𝑦(𝑡))∈Λ, by the Itô formula, we compute the following: 𝑑𝑉₁=𝑅₀𝑎−1𝑥1−𝑞𝑥−𝑦−𝑎𝑥𝑦+𝜎𝑥+𝑝𝑦2²¹𝑑𝑡+𝜎₁𝑅₀𝑎𝑥−1𝑑𝐵₁(𝑡)+11−𝑦𝑅₀𝑥𝜎𝑥+𝑝𝑦−1𝑦+²²2𝑑𝑡+𝜎₂𝑦−1𝑑𝐵₂(𝑡)=𝐿𝑉₁𝑑𝑡+𝜎₁𝑅₀𝑎𝑥−1𝑑𝐵₁(𝑡)+𝜎₂𝑦−1𝑑𝐵₂(𝑡), (40) where 𝐿𝑉₁=𝑞+𝑎𝑦−𝑅₀𝑥−𝑅𝑥+𝑝𝑦⁰𝑞𝑎𝑥+𝑅₀𝑎−1𝑥𝜎1−𝑦+1−𝑦+²¹2+𝜎²₂2𝑎≤𝑞+𝑝+𝑅₀𝑎−1𝑥𝜎1−𝑦+1−𝑦+²¹2+𝜎²₂2. (41)Case  1 (assume 𝑎≥𝑅₀). In this case, we have 𝑥≥1. It follows that 𝑅₀𝑎−1𝑥≤𝑅1−𝑦+1−𝑦₀𝑎1+𝑦𝑥−1−1𝑥𝑅+1≤1+⁰𝑎. (42)Case  2 (assume 𝑎<𝑅₀). If 𝑎/𝑅₀≤𝑥<1, one has the following:(a)((𝑅₀/𝑎)−(1/𝑥))(1−𝑦)+1−𝑦≤(𝑅₀/𝑎)+((𝑦−1)/𝑥)+1−𝑦≤1+(𝑅₀/𝑎) provided that 0<𝑦≤1;(b)((𝑅₀/𝑎)−(1/𝑥))(1−𝑦)+1−𝑦≤1 provided that 𝑦>1.
Hence, there exists a positive number 𝑀 independent on 𝑥, 𝑦 and 𝑡 such that 𝐿𝑉₁≤𝑀. Substituting this inequality into (40), we can get the following: 𝑑𝑉₁≤𝑀𝑑𝑡+𝜎₁𝑅₀𝑎𝑥−1𝑑𝐵₁(𝑡)+𝜎₂𝑦−1𝑑𝐵₂(𝑡). (43)
Integrating both sides of the above inequality from 0 to 𝜏_𝑛∧𝑇 and taking expectations leads to the following: 𝐸𝑉₁𝑥𝜏_𝑛∧𝑇,𝑦𝜏_𝑛∧𝑇≤𝑉₁𝑥(0),𝑦(0)+𝑀𝑇. (44)
Set Ω_𝑛={𝜏_𝑛≤𝑇}, for 𝑛≥𝑛₁ and consider inequality (38), we can get 𝒫(Ω_𝑛)≥𝜀. Note that for every 𝜔∈Ω_𝑛, there exists some 𝑖 such that 𝑥_𝑖(𝜏_𝑛,𝜔) equals either 𝑛 or 1/𝑛 for 𝑖=1,2; hence, 𝑉₁𝑥𝜏_𝑛,𝜔,𝑦𝜏_𝑛),1,𝜔≥min(𝑛−1−ln𝑛𝑛1−1−ln𝑛. (45) It then follows from (44) that 𝑉₁𝑥(0)(0)𝐼,𝑦+𝑀𝑇≥𝐸_Ω𝑛(𝜔)𝑉₁𝑥𝜏_𝑛,𝑦𝜏_𝑛),1≥𝜖min(𝑛−1−ln𝑛𝑛1−1−ln𝑛, (46) where 𝐼_Ω𝑛 is the indicator function of Ω_𝑛.
As 𝑛→∞ we have the following: ∞>𝑉₁𝑥(0),𝑦(0)+𝑀𝑇=∞a.s., (47) which leads to the contradiction. This completes the proof.

3.2. Stochastic Boundedness and Permanence

Theorem 6 shows that the solutions to model (8) will remain in Λ. Generally speaking, the nonexplosion property, the existence, and the uniqueness of the solution are not enough but the property of boundedness and permanence are more desirable since they mean the long-time survival in the population dynamics. Now, we present the definition of stochastic ultimate boundedness and stochastic permanence [31].

Definition 7. The solutions 𝑋(𝑡)=(𝑥(𝑡),𝑦(𝑡)) of model (8) are said to be stochastically ultimately bounded, if for any 𝜀∈(0,1), there is a positive constant 𝛿=𝛿(𝜀), such that for any initial value (𝑥(0),𝑦(0))∈Λ, the solution 𝑋(𝑡) of model (8) has the property that limsup_𝑡→∞𝑃{|𝑋(𝑡)|>𝛿}<𝜀. (48)

Definition 8. The solutions 𝑋(𝑡)=(𝑥(𝑡),𝑦(𝑡)) of model (8) are said to be stochastically permanent if for any 𝜀∈(0,1), there exists a pair of positive constants 𝛿=𝛿(𝜀) and 𝜒=𝜒(𝜀), such that for any initial value (𝑥(0),𝑦(0))∈Λ, the solution 𝑋(𝑡) of model (8) has the property that liminf_𝑡→∞𝑃{|𝑋(𝑡)|≥𝛿}≥1−𝜀,liminf_𝑡→∞𝑃|𝑋(𝑡)|≤𝜒≥1−𝜀. (49)

Theorem 9. The solutions of model (8) are stochastically ultimately bounded for any initial value (𝑥(0),𝑦(0))∈Λ.

Proof. Denote functions 𝑉₂=𝑒^𝑡𝑥^𝜃,𝑉₃=𝑒^𝑡𝑦^𝜃 (50) for (𝑥,𝑦)∈Λ and 0<𝜃<1.
Applying the Itô formula leads to the following: 𝑑𝑉₂=𝐿𝑉₂𝑑𝑡+𝜎₁𝜃𝑒^𝑡𝑥^𝜃𝑑𝐵₁(𝑡),𝑑𝑉₃=𝐿𝑉₃𝑑𝑡+𝜎₂𝜃𝑒^𝑡𝑦^𝜃𝑑𝐵₂(𝑡), (51) where 𝐿𝑉₂=𝑒^𝑡𝑥^𝜃11+𝜃𝑥𝑦−𝑞−𝑥−𝑎𝑦+𝜎𝑥+𝑝𝑦²¹𝜃(𝜃−1)2,𝐿𝑉₃=𝑒^𝑡𝑦^𝜃𝑅1+𝜃⁰𝑥+𝜎𝑥+𝑝𝑦−1²²𝜃(𝜃−1)2. (52) Thus, there exists the positive constants 𝑀₁ and 𝑀₂ such that we have 𝐿𝑉₂<𝑀₁𝑒^𝑡 and 𝐿𝑉₃<𝑀₂𝑒^𝑡. It follows that 𝑒^𝑡𝐸𝑥^𝜃−𝐸𝑥(0)^𝜃≤𝑀₁𝑒^𝑡 and 𝑒^𝑡𝐸𝑦^𝜃−𝐸𝑦(0)^𝜃≤𝑀₂𝑒^𝑡. Then we get the following: limsup_𝑡→∞𝐸𝑥^𝜃≤𝑀₁<+∞,limsup_𝑡→∞𝐸𝑦^𝜃≤𝑀₂<+∞. (53)
Note that |𝑋(𝑡)|^𝜃=𝑥²(𝑡)+𝑦²(𝑡)^𝜃/2≤2^𝜃/2max𝑥^𝜃(𝑡),𝑦^𝜃(𝑡)≤2^𝜃/2𝑥^𝜃+𝑦^𝜃. (54)
Therefore, we obtain the following: limsup_𝑡→∞𝐸|𝑋(𝑡)|^𝜃≤2^𝜃/2𝑀₁+𝑀₂<+∞. (55)
As a result, there exists a positive constant 𝛿₁ such that limsup_𝑡→∞𝐸|𝑋(𝑡)|<𝛿₁. (56) Now, for any 𝜀>0, let 𝛿=𝛿²₁/𝜀²; then by Chebyshev’s inequality, 𝒫𝑋(𝑡}≤𝐸{|)|>𝛿|𝑋(𝑡)|√𝛿. (57) Hence, limsup_𝑡→∞𝒫𝑋(𝑡}≤𝛿{|)|>𝛿¹√𝛿=𝜀, (58) which yields the required assertion.

We are now in the position to show the stochastic permanence. Let us present some hypothesis and a useful lemma.

Lemma 10. Assume 𝑅₀>𝑎+max{4,2𝑝𝑞}. For any initial value (𝑥(0),𝑦(0))∈Λ, the solution (𝑥(𝑡),𝑦(𝑡)) satisfies that limsup_𝑡→∞𝐸1|𝑋(𝑡)|^𝜌≤𝐻, (59) where 𝜌 is an arbitrary positive constant satisfying 𝜌+12max𝜎₁,𝜎₂²𝑅<1+min⁰−𝑎𝑅2𝑝−𝑞,⁰−𝑎2,2−2𝐻=^𝜌𝐶₂+4𝑘𝐶₁4𝑘𝐶₁×max1,2𝐶₁+𝐶₂+𝐶²₂+4𝐶₁𝐶₂2𝐶₁^𝜌−2 (60)(61) in which 𝑘 is an arbitrary positive constant satisfying 𝜌𝜌+12max𝜎₁,𝜎₂²𝑅+𝑘<𝜌+𝜌min⁰−𝑎𝑅2𝑝−𝑞,⁰−𝑎2−2 (62) with 𝐶₁𝑅=𝜌+𝜌min⁰−𝑎𝑅2𝑝−𝑞,⁰−𝑎2−𝜌−2𝜌+12max𝜎₁,𝜎₂²𝐶−𝑘>0,₂=𝜌max𝑞,2+𝑎+𝜌𝑅₀𝑅₀−1max1,𝑝²𝜎2𝑎𝑝+𝜌max₁,𝜎₂²+2𝑘>0. (63)

Proof. Set 𝑈(𝑥,𝑦)=1/(𝑥+𝑦) for (𝑥(𝑡),𝑦(𝑦))∈Λ, by the Itô formula, we have the following: 𝑑𝑈=−𝑈²1−𝑞𝑥−𝑦−𝑎𝑥𝑦+𝑅𝑥+𝑝𝑦⁰𝑥𝑦𝑥+𝑝𝑦−𝑦𝑑𝑡+𝑈³𝜎²₁𝑥²+𝜎²₂𝑦²𝑑𝑡−𝑈²𝜎₁𝑥𝑑𝐵₁(𝑡)+𝜎₂𝑦𝑑𝐵₂(𝑡)=𝐿𝑈𝑑𝑡−𝑈²𝜎₁𝑥𝑑𝐵₁(𝑡)+𝜎₂𝑦𝑑𝐵₂(𝑡), (64) where 𝐿𝑈=−𝑈²𝑅1−𝑞𝑥−2𝑦+⁰−𝑎𝑥𝑦𝑥+𝑝𝑦+𝑈³𝜎²₁𝑥²+𝜎²₂𝑦². (65)
Choose a positive constant 𝜌 such that it satisfies (60). Applying the Itô formula again, we can get the following: 𝐿(1+𝑈)^𝜌=𝜌(1+𝑈)^𝜌−1+𝜌𝐿𝑈𝜌−12𝑈⁴(1+𝑈)^𝜌−2𝜎²₁𝑥²+𝜎²₂𝑦²=(1+𝑈)^𝜌−2Φ, (66) where Φ=−𝜌𝑈²𝑅1−𝑞𝑥−2𝑦+⁰−𝑎𝑥𝑦𝑥+𝑝𝑦−𝜌𝑈³𝑅1−𝑞𝑥−2𝑦+⁰−𝑎𝑥𝑦𝑥+𝑝𝑦+𝜌𝑈³𝜎²₁𝑥²+𝜎²₂𝑦²+𝜌1+𝜌𝑈⁴2𝜎²₁𝑥²+𝜎²₂𝑦²≤−𝜌𝑈²+𝜌𝑈²𝑞𝑥+(2+𝑎)𝑦−𝜌𝑈³𝑅₀−𝑎𝑅2𝑝−𝑞𝑥+⁰−𝑎2−2𝑦+𝜌𝑈³𝑅₀−𝑎𝑥²+𝑝²𝑦²2𝑝𝑥+𝑝𝑦+𝜌𝑈³𝜎²₁𝑥²+𝜎²₂𝑦²+𝜌1+𝜌𝑈⁴2𝜎²₁𝑥²+𝜎²₂𝑦². (67)
Using the facts that 𝑈³𝜎²₁𝑥²+𝜎²₂𝑦²<max𝜎₁,𝜎₂²𝑈𝑈,₄𝜎²₁𝑥²+𝜎²₂𝑦²<max𝜎₁,𝜎₂²𝑈², (68) so, Φ≤−𝑈²𝑅𝜌+𝜌min⁰−𝑎𝑅2𝑝−𝑞,⁰−𝑎2−𝜌−2𝜌+12max𝜎₁,𝜎₂²+𝑈𝜌max𝑞,2+𝑎+𝜌𝑅₀𝑅₀−1max1,𝑝²𝜎2𝑎𝑝+𝜌max₁,𝜎₂². (69) Now, let 𝑘>0 sufficiently small such that it satisfies (62), by the Itô formula; then 𝐿𝑒^𝑘𝑡(1+𝑈)^𝜌=𝑘𝑒^𝑘𝑡(1+𝑈)^𝜌+𝑒^𝑘𝑡𝐿(1+𝑈)^𝜌=𝑒^𝑘𝑡(1+𝑈)^𝜌−2𝑘(1+𝑈)²+Φ≤𝑒^𝑘𝑡(1+𝑈)^𝜌−2−𝐶₁𝑈²+𝐶₂𝑢+𝑘≤𝐻₁𝑒^𝑘𝑡, (70) where 𝐻₁=((𝐶₂+4𝑘𝐶₁)/4𝐶₁)max{1,((2𝐶₁+𝐶₂+𝐶²₂+4𝐶₁𝐶₂)/2𝐶₁)^𝜌−2} and 𝐶₁, 𝐶₂ have been defined in the statement of the theorem. Thus, 𝐸𝑒^𝑘𝑡(1+𝑈)^𝜌≤(1+𝑈(0))^𝜌+𝐻₁𝑘𝑒^𝑘𝑡. (71) So we can have the following: limsup_𝑡→∞𝐸𝑈(𝑡)^𝜌≤limsup_𝑡→∞𝐸(1+𝑈)^𝜌≤𝐻₁𝑘. (72)
In addition, we know that (𝑥+𝑦)^𝜌≤2^𝜌(𝑥²+𝑦²)^𝜌/2=2^𝜌|𝑋(𝑡)|^𝜌; consequently, limsup_𝑡→∞𝐸1|𝑋(𝑡)|^𝜌≤2^𝜌limsup_𝑡→∞𝐸𝑈(𝑡)^𝜌≤2^𝜌𝐻₁𝑘=𝐻, (73) which complets the proof.

Consider Chebyshev inequality, Theorem 9, and Lemma 10 together, we immediately obtain the following result.

Theorem 11. If the following conditions are satisfied(i)𝑎+max{2𝑝𝑞,4}<𝑅₀;(ii)(1/2)(max{𝜎₁,𝜎₂})²<1+min{((𝑅₀−𝑎)/2𝑝)−𝑞,((𝑅₀−𝑎)/2)−2},
then the solutions of model (8) is stochastically permanent.

4. Conclusions and Discussions

In this paper, by using the theory of stochastic differential equation, we investigate the dynamics of an SIRS epidemic model with a ratio-dependent incidence rate. The value of this study lies in two aspects. First, it presents some relevant properties of the deterministic model (7), including boundedness, dissipation, persistence, and the stability of the disease-free and endemic points. Second, it verifies the existence of global positive solutions, stochastic boundedness, and permanence for the stochastic model (8).

As an example, we give some numerical examples to illustrate the dynamical behavior of stochastic model (8) by using the Milstein method mentioned in Higham [49]. In this way, model (8) can be rewritten as the following discretization equations: 𝑥_𝑘+1=𝑥_𝑘+1−𝑞𝑥_𝑘−𝑦_𝑘−𝑎𝑥_𝑘𝑦_𝑘𝑥_𝑘+𝑝𝑦_𝑘Δ𝑡+𝜎₁𝑥_𝑘√Δ𝑡𝜉_𝑘+𝜎²₁2𝑥²_𝑘𝜉²_𝑘𝑦−1Δ𝑡,_𝑘+1=𝑦_𝑘+𝑅₀𝑥_𝑘𝑦_𝑘𝑥_𝑘+𝑝𝑦_𝑘−𝑦_𝑘Δ𝑡+𝜎₂𝑦_𝑘√Δ𝑡𝜂_𝑘+𝜎²₂2𝑦²_𝑘𝜂²_𝑘−1Δ𝑡, (74) where 𝜉_𝑘 and 𝜂_𝑘, 𝑘=1,2,…,𝑛, are the Gaussian random variables 𝑁(0,1).

The parameters of model (8) are fixed as (14). In this case, model (7) has the endemic point 𝐸^∗=(0.11,0.74). And model (8) becomes as follows: 𝑑𝑥=1−2𝑥−𝑦−0.3𝑥𝑦𝑥+0.5𝑦𝑑𝑡+𝜎₁𝑥𝑑𝐵₁(𝑡),𝑑𝑦=4.5𝑥𝑦𝑥+0.5𝑦−𝑦𝑑𝑡+𝜎₂𝑦𝑑𝐵₂(𝑡). (75) Simple computations show that 𝑎+max2𝑝𝑞,4=4.3<4.5=𝑅₀,0.03²2=12max𝜎₁,𝜎₂²𝑅<1+min⁰−𝑎𝑅2𝑝−𝑞,⁰−𝑎2𝜎−2=1.1,if₁,𝜎₂=(0.03,0.01),0.5²2=12max𝜎₁,𝜎₂²𝑅<1+min⁰−𝑎𝑅2𝑝−𝑞,⁰−𝑎2𝜎−2=1.1,if₁,𝜎₂=(0.5,0.3). (76) It is easy to see that, all the conditions of Theorem 11 are satisfied, and we can therefore conclude that, with (𝜎₁,𝜎₂)=(0.03,0.01) and (𝜎₁,𝜎₂)=(0.5,0.3), the solutions of model (8) is stochastically permanent. The numerical examples shown in Figures 2 and 3 clearly support these results. In Figure 2, with (𝜎₁,𝜎₂)=(0.03,0.01), the solutions of model (8) will be oscillating slightly around the endemic point 𝐸^∗=(0.11,0.74) of model (7). And in Figure 3, with (𝜎₁,𝜎₂)=(0.5,0.3), the solutions of model (8) will be oscillating strongly around the endemic point 𝐸^∗=(0.11,0.74) of model (7).

Figure 2 

The solution of the stochastic model (8) with initial values 𝑥(0)=0.2, 𝑦(0)=0.15. The parameters are taken as (14), 𝜎1=0.03, 𝜎2=0.01.

Figure 3 

The solution of the stochastic model (8) with initial values 𝑥(0)=0.2, 𝑦(0)=0.15. The parameters are taken as (14), 𝜎1=0.5, 𝜎2=0.3.

It is worthy to note that, throughout this paper, the parameters for model (7), also for model (8), are fixed as the set (14). The reason is that with this parameter set, the conditions of our theoretical results hold. Of course, one can adopt other parameters set to show the numerical results.

From the theoretical and numerical results, we can know that, when the noise density is not large, the stochastic model (8) preserves the property of the stability of the deterministic model (7). To a great extent, we can ignore the noise and use the deterministic model (7) to describe the population dynamics. However, when the noise is sufficiently large, it can force the population to become largely fluctuating. In this case, we cannot use deterministic model (7) but stochastic model (8) to describe the population dynamics. Needless to say, both deterministic and stochastic epidemic models have their important roles.

Furthermore, from the numerical results in Figure 2, one can see that model (8) is stochastically stable. But we cannot prove the stochastic stability because of the complexity of model (8). This can be further investigated.

On the other hand, we know that there are different possible approaches to including random effects in the epidemic models affected by environmental white noise, here we consider another method to introduce random effects in the epidemic model (7). The martingale approach was initiated by Beretta et al. [35] and applied in [27, 30, 45, 47]. They introduced stochastic perturbation terms into the growth equations to incorporate the effect of a randomly fluctuating environment. In detail, assume that the stochastic perturbations of the state variables around their steady-state 𝐸^∗ are of a white noise type which is proportional to the distances of 𝑥, 𝑦 from their steady-state values 𝑥^∗ and 𝑦^∗, respectively. In this way, model (7) will be reduced to the following form: 𝑑𝑥=1−𝑞𝑥−𝑦−𝑎𝑥𝑦𝑥+𝑝𝑦𝑑𝑡+𝜎₁𝑥−𝑥^∗𝑑𝐵₁(𝑡),𝑅𝑑𝑦=⁰𝑥𝑦𝑥+𝑝𝑦−𝑦𝑑𝑡+𝜎₂𝑦−𝑦^∗𝑑𝐵₂(𝑡), (77) where the definitions of 𝜎₁, 𝜎₂ and 𝐵₁(𝑡), 𝐵₂(𝑡) are the same as in (8).

If 𝑅₀>1, stochastic model (77) can center at its endemic point 𝐸^∗, with the change of variables 𝑢=𝑥−𝑥^∗, 𝑣=𝑦−𝑦^∗. The linearized version of model (77) is as follows: 𝑑𝑧(𝑡)=𝑓₁(𝑧(𝑡))𝑑𝑡+𝑓₂(𝑧(𝑡))𝑑𝐵(𝑡), (78) where 𝑧(𝑡)=𝑢(𝑡)𝑣(𝑡),𝑓₁=𝐽₁₁𝑢(𝑡)+𝐽₁₂𝑣(𝑡)𝐽₂₁𝑢(𝑡)+𝐽₂₂𝑣(𝑡),𝑓₂=𝜎₁𝑢(𝑡)00𝜎₂𝑣(𝑡), (79) where 𝐽₁₁, 𝐽₁₂, 𝐽₂₁, 𝐽₂₂ are defined as (12).

It is easy to see that the stability of the endemic point 𝐸^∗ of model (77) is equivalent to the stability of zero solution of model (78).

Before proving the stochastic stability of the zero solution of model (78), we put forward a lemma in [50].

Lemma 12. Suppose there exists a function 𝑉(𝑧,𝑡)∈𝐶²(Ω) satisfying the following inequalities: 𝐾₁|𝑧|^𝜔≤𝑉(𝑧,𝑡)≤𝐾₂|𝑧|^𝜔,𝐿𝑉(𝑧,𝑡)≤−𝐾₃|𝑧|^𝜔, (80)(81) where 𝜔>0 and 𝐾_𝑖(𝑖=1,2,3) is positive constant. Then the zero solution of mode (78) is exponentially 𝜔-stable for all time 𝑡≥0.

From the lemma above, note that if 𝜔=2 in (80) and (81), then the zero solution of model (78) is stochastically asymptotically stable in probability. Thus, we obtain the following theorem.

Theorem 13. Assume that 𝜎²₁<2(𝑝𝑞𝑅²₀+𝑎(𝑅₀−1)²)/𝑝𝑅²₀,𝜎²₂<2(𝑅₀−1)/𝑅₀ hold; then the zero solution of model (78) is asymptotically mean square stable. And the endemic point 𝐸^∗ of model (77) is asymptotically mean square stable.

The details of the proof are shown in the Appendix.

We should point out that the results obtained in this paper are only for the simple case when 𝑙=ℎ=1 of the incidence rate (3). The dynamical behaviors of the stochastic epidemic model with general ratio-dependent incidence rate (3) are desirable in future studies.

Appendix

The proof of Theorem 13

Proof. Let us consider the Lyapunov function 𝑉₅(𝑧(𝑡=1))2𝑢²+𝜅𝑣², (A.1) where 𝜅=(𝑅²₀+𝑎)/𝑅₀(𝑅₀−1)².
It is easy to check that inequality (80) holds with 𝜔=2. Moreover, 𝐿𝑉₅(𝑧(𝑡))=𝑢𝐽₁₁𝑢+𝐽₁₂𝑣+𝜅𝑣𝐽₂₁𝑢+𝐽₂₂𝑣+12𝜎²₁𝑢²+𝜅𝜎²₂𝑣²=𝐽₁₁+𝜎²₁2𝑢²+𝜅𝐽₂₂+𝜎²₂2𝑣²=−𝑧^𝑇𝑄𝑧, (A.2) where 𝐽𝑄=¹¹+𝜎²₁200𝜅𝐽₂₂+𝜎²₂2. (A.3)
When 𝜎²₁<2(𝑝𝑞𝑅²₀+𝑎(𝑅₀−1)²)/𝑝𝑅²₀, 𝜎²₂<2(𝑅₀−1)/𝑅₀, the two eigenvalues 𝜆₁, 𝜆₂ of the matrix 𝑄 will be positive. Set 𝜆_min=min{𝜆₁,𝜆₂}, it follows from (A.2) immediately that 𝐿𝑉₅(𝑧(𝑡))≤−𝜆_min|𝑧(𝑡)|². (A.4) We therefore have the assertion.

Acknowledgments

The authors would like to thank the editors and referees for their helpful comments and suggestions. They also thank Natural Science Foundation of Zhejiang Province (LY12A01014) and the National Basic Research Program of China (2012CB426510).