The Complex Dynamics of a Stochastic Predator-Prey Model

Abstract

A modified stochastic ratio-dependent Leslie-Gower predator-prey model is formulated and analyzed. For the deterministic model, we focus on the existence of equilibria, local, and global stability; for the stochastic model, by applying Itô formula and constructing Lyapunov functions, some qualitative properties are given, such as the existence of global positive solutions, stochastic boundedness, and the global asymptotic stability. Based on these results, we perform a series of numerical simulations and make a comparative analysis of the stability of the model system within deterministic and stochastic environments.

1. Introduction

The dynamic relationship between predators and their prey has long been and will continue to be one of the dominant themes in both ecology and mathematical ecology due to its universal existence and importance [1]. In recent years, one of important predator-prey models is Leslie-Gower model [2, 3], which has been extensively studied [4, 5]. And, more and more obvious evidences of biology and physiology show that in many conditions, especially when the predators have to search for food (consequently, have to share or compete for food), a more realistic and general predator-prey system should rely on the theory of ratio-dependence, this is confirmed by lots of experimental results [4, 6]. A ratio-dependent Leslie-Gower predator-prey model [7], takes the form 𝑑𝑥()𝑥𝑑𝑡=𝑎𝑥1−𝑏𝑥−𝑝𝑦𝑦,𝑑𝑦𝑑𝑡=𝑒𝑦1−𝑓𝑦𝑥, (1.1) where 𝑥=𝑥(𝑡), 𝑦=𝑦(𝑡) stand for the population (the density) of the prey and the predator at time 𝑡, respectively, and 𝑝(𝑥/𝑦) is the predator functional response to prey. And we assumed that the prey grows logistically with growth rate 𝑎 and carrying capacity 1/𝑏 in the absence of predation. The predator consumes the prey according to the functional response 𝑝(𝑥/𝑦) and grow logistically with growth rate 𝑒 and carrying capacity 𝑥/𝑓 proportional to the population size of prey (or prey abundance). The parameter 𝑓 is a measure of the food quality that the prey provides for conversion into predator birth. The term 𝑓𝑦/𝑥 of this equation is called the Leslie-Gower term.

On the other hand, the predator 𝑦 can switch over to other population when the prey population is severely scarce, but its growth will be limited, because we cannot forget the fact that its most favorite food, the prey 𝑥, is not in abundance. In this situation, a positive constant 𝑚 can be added to the denominator, 𝑚 measures the extent to which the environment provides protection to the predator [8, 9], and the second equation of model (1.1) becomes 𝑑𝑦𝑑𝑡=𝑒𝑦1−𝑓𝑦.𝑚+𝑥 (1.2)

Based on the above discussions, in the paper, we will focus on the following ratio-dependent Leslie-Gower model: 𝑑𝑥𝑑𝑡=𝑎𝑥1−𝑏𝑥−𝑐𝑦,𝑥+𝑛𝑦≜𝑓𝑥,𝑦𝑑𝑦𝑑𝑡=𝑒𝑦1−𝑓𝑦.𝑚+𝑥≜𝑔𝑥,𝑦 (1.3)

The rest of the paper is organized as follows. In Section 2, we give some theorems about the stability property of the equilibria of model (1.3). In Section 3, we establish a stochastic model based on model (1.3) and focus on the existence of global positive solutions, stochastic boundedness, and the global asymptotic stability of the stochastic model. In Section 4, we give some numerical examples and make a comparative analysis of the stability of the model system within deterministic and stochastic environments.

2. The Dynamics of Model (1.3)

2.1. Dissipativeness

By standard simple arguments, one can show that the solution of model (1.3) always exists and stays positive. In fact, from the first equation of model (1.3), we can get limsup_𝑡→∞𝑥(𝑡)≤1𝑏. (2.1) Hence, there exists a 𝑇>0 such that, for 𝑡>𝑇, 𝑥≤1/𝑏.

From the second equation of model (1.3), we see that, for 𝑡>𝑇, we have 𝑑𝑦/𝑑𝑡≤𝑒𝑦(1−𝑏𝑓𝑦/(𝑚𝑏+1)). A standard comparison argument shows that limsup_𝑡→∞𝑦(𝑡)≤(1+𝑏𝑚).𝑏𝑓 (2.2)

Thus, we have the following conclusion.

Lemma 2.1. Model (1.3) is dissipative.

Lemma 2.2. If 𝑐<𝑛, then model (1.3) is permanent.

Proof. If 𝑐<𝑛, from the first equation of model (1.3), we have 𝑑𝑥/𝑑𝑡≥𝑎𝑥(1−𝑏𝑥−𝑐/𝑛). Therefore, by standard comparison argument, we have liminf_𝑡→∞𝑥(𝑡)≥𝑛−𝑐.𝑏𝑛 (2.3)
Hence, for any 𝜀>0 and large 𝑡, 𝑥(𝑡)>(𝑛−𝑐)/𝑏𝑛−𝜀, and 𝑑𝑦𝑑𝑡≥𝑒𝑦𝑏𝑚𝑛+𝑛𝑐−𝑏𝑛𝜀−𝑏𝑓𝑛𝑦.𝑏𝑚𝑛+𝑛𝑐−𝑏𝑛𝜀 (2.4)
From the arbitrariness of 𝜀>0, we can get that liminf_𝑡→∞𝑦(𝑡)≥𝑏𝑚𝑛+𝑛−𝑐.𝑏𝑓𝑛 (2.5)

2.2. Stability Analysis of the Equilibria

In this section, we will focus on the existence of equilibria and their stabilities of model (1.3).

It is easy to find that model (1.3) always has three boundary equilibria 𝐸₀=(0,0), 𝐸₁=(1/𝑏,0), 𝐸₂=(0,𝑚/𝑓). And the positive equilibria (𝑥,𝑦) satisfies the equations 𝑎𝑥1−𝑏𝑥−𝑐𝑦𝑥+𝑛𝑦=0,𝑒𝑦1−𝑓𝑦𝑚+𝑥=0, (2.6) which yields 𝑏𝑓+𝑛𝑥²−𝑓−𝑐−𝑏𝑚𝑛+𝑛𝑥+(𝑐−𝑛)𝑚=0. (2.7)

For simplicity, we define 𝐴₁=𝑓−𝑐−𝑏𝑚𝑛+𝑛𝑏,𝐵𝑓+𝑛₁=(𝑐−𝑛)𝑚𝑏,𝐶𝑓+𝑛₁=𝑓−𝑐−𝑏𝑚𝑛+𝑛²𝑏²𝑓+𝑛²−4(𝑐−𝑛)𝑚𝑏=𝑓+𝑛𝐴²₁−4𝐵₁, (2.8) then (2.7) can be rewritten as 𝑥²−𝐴₁𝑥+𝐵₁=0. (2.9)

Lemma 2.3 (existence of equilibria). (a) Suppose 𝐴₁>0 and 𝐵₁>0, in (2.9), then one has(a1) If 𝐶²₁>0, then it has two positive roots given by 𝑥_𝑒1=𝐴₁+𝐶₁2,𝑥_𝑒2=𝐴₁−𝐶₁2. (2.10) Therefore, model (1.3) has two positive equilibria 𝐸₃=(𝑥_𝑒1,𝑦_𝑒1)=((𝐴₁+𝐶₁)/2,(2𝑚 + 𝐴₁+𝐶₁)/2𝑓) and 𝐸₄=(𝑥_𝑒2,𝑦_𝑒2)=((𝐴₁−𝐶₁)/2,(2𝑚+𝐴₁−𝐶₁)/2𝑓).(a2) If 𝐶²₁=0, (2.9) has a unique positive root of multiplicity-2 given by 𝑥_𝑒=𝑚(𝑐−𝑛)/𝑏(𝑓+𝑛). Thus, model (1.3) has a unique positive equilibrium 𝐸_𝑒=𝑚(𝑐−𝑛)𝑏,+𝑓+𝑛𝑏𝑚𝑓+𝑛𝑏𝑚𝑓+𝑛(𝑐−𝑛)𝑏𝑓𝑓+𝑛. (2.11)(a3) If 𝐶²₁<0, there are no positive roots that exist and then, model (1.3) has no positive equilibrium.
(b)𝐵₁=0,𝐴₁>0, (2.9) has a unique positive root 𝑥=(𝑓−𝑏𝑚𝑛)/𝑏(𝑓+𝑛), model (1.3) has a unique positive equilibrium 𝐸=(𝑥,𝑦)=((𝑓−𝑏𝑚𝑛)/𝑏(𝑓+𝑛),(𝑏𝑚+1)/𝑏(𝑓+𝑛)).
(c)𝐵₁<0, (2.9) has a unique positive root 𝑥^∗=(𝐴₁+𝐶₁)/2, model (1.3) has a unique positive equilibrium 𝐸^∗=(𝑥^∗,𝑦^∗)=((𝐴₁+𝐶₁)/2,(2𝑚+𝐴₁+𝐶₁)/2𝑓).

Next, we discuss the local stabilities of these equilibria. Easy to obtain the following results:

Lemma 2.4. (i)  𝐸₁=(1/𝑏,0) is a saddle point;(ii) if 𝑐>𝑛, 𝐸₂=(0,𝑚/𝑓) is a stable node point. If 𝑐≤𝑛, 𝐸₂=(0,𝑚/𝑓) is a saddle point.

Following, let 𝐸=(𝑥,𝑦) be an arbitrary positive equilibrium. The Jacobian matrix for 𝐸=(𝑥,𝑦) is given by 𝐽(𝐸)=𝑎𝑥−𝑏+𝑐𝑦𝑥+𝑛𝑦²−𝑐𝑎𝑥+𝑥+𝑛𝑦𝑐𝑛𝑦𝑥+𝑛𝑦²𝑒𝑓−𝑒. (2.12)

Then we can get the following: tr(𝐽(𝐸))=𝑎𝑥−𝑏+𝑐𝑦𝑥+𝑛𝑦²(𝐽(𝐸=𝑏−𝑒,det))𝑎𝑒𝑥𝑓+𝑛²𝑥²+2𝑏𝑚𝑛𝑓+𝑛𝑥−𝑐𝑓𝑚+𝑏𝑚²𝑛²𝑓𝑥+𝑚𝑛+𝑛𝑥². (2.13)

So the sign of det(𝐽(𝐸)) is determined by 𝐹(𝑥)≜𝑏𝑓+𝑛²𝑥²+2𝑏𝑚𝑛𝑓+𝑛𝑥−𝑐𝑓𝑚+𝑏𝑚²𝑛². (2.14)

Theorem 2.5. For model (1.3), the stabilities of 𝐸₃ and 𝐸₄ are as follows.(a) The positive equilibrium 𝐸₃=(𝑥_𝑒1,(𝑚+𝑥_𝑒1)/𝑓)=((𝐴₁+𝐶₁)/2,(2𝑚+𝐴₁+𝐶₁)/2𝑓) is(a1) stable if and only if 𝑒>𝑎𝑥_𝑒1−𝑏+𝑐𝑓𝑚+𝑥_𝑒1𝑓𝑥_𝑒1+𝑛𝑚+𝑥_𝑒1²; (2.15)(a2) unstable if and only if 𝑒<𝑎𝑥_𝑒1−𝑏+𝑐𝑓𝑚+𝑥_𝑒1𝑓𝑥_𝑒1+𝑛(𝑚+𝑥_𝑒1)². (2.16)(b) The positive equilibrium 𝐸₄=(𝑥_𝑒2,(𝑚+𝑥_𝑒2)/𝑓)=((𝐴₁−𝐶₁)/2,(2𝑚+𝐴₁−𝐶₁)/2𝑓) is a saddle point.

Proof. (a) At the point 𝐸₃, we have the following: 𝐹𝑥_𝑒1=𝑏𝑓+𝑛²𝐴₁+𝐶₁²4+𝑏𝑚𝑛𝑓+𝑛𝐴₁+𝐶₁+𝑏𝑚²𝑛²=𝑏−𝑐𝑓𝑚𝑓+𝑛²𝐶²₁2+𝐶₁+𝑏𝑏𝑚𝑛𝑓+𝑛𝑓+𝑛²𝐴₁2>0, (2.17) thus, det(𝐽(𝐸₃))>0, and the nature of singularity 𝐸₃ depends on the trace given by tr𝐽𝐸₃=𝑎𝑥_𝑒1−𝑏+𝑐𝑓𝑚+𝑥_𝑒1𝑓𝑥_𝑒1+𝑛(𝑚+𝑥_𝑒1)²−𝑒. (2.18) Clearly, (a1) If 𝑒>𝑎𝑥_𝑒1(−𝑏+𝑐𝑓(𝑚+𝑥_𝑒1)/(𝑓𝑥_𝑒1+𝑛(𝑚+𝑥_𝑒1))²), then tr(𝐽(𝐸₃))<0, and 𝐸₃=(𝑥_𝑒1,(𝑚+𝑥_𝑒1)/𝑓) is stable.(a2) If 𝑒<𝑎𝑥_𝑒1(−𝑏+𝑐𝑓(𝑚+𝑥_𝑒1)/(𝑓𝑥_𝑒1+𝑛(𝑚+𝑥_𝑒1))²), then tr(𝐽(𝐸₃))>0, and 𝐸₃=(𝑥_𝑒1,(𝑚+𝑥_𝑒1)/𝑓) is unstable.(b) At the point 𝐸₄, we have the following: 𝐹𝑥_𝑒2=𝑏𝑓+𝑛²𝐴₁−𝐶₁²4+𝑏𝑚𝑛𝑓+𝑛𝐴₁−𝐶₁+𝑏𝑚²𝑛²=𝑏−𝑐𝑓𝑚𝑓+𝑛²2𝐶₁𝐴₁−𝐶₁−𝑏𝑚𝑛𝑓+𝑛𝐶₁−2𝑓𝑚𝑓+𝑛𝐵₁<0, (2.19) then det(𝐽(𝐸₄))<0, and 𝐸₄ is a saddle point.

Figure 1 shows the dynamics of model (1.3). In this case, 𝐸₀=(0,0) is a nodal source, 𝐸₁=(4,0) is a saddle point, 𝐸₂=(0,071429) is a nodal sink, 𝐸₄=(0.31772,1.5846) is a saddle point, and 𝐸₃=(0.96454,2.3932) is locally asymptotically stable. There exists a separatrix curve determined by the stable manifold of 𝐸₄, which divides the behavior of trajectories, that is, the stable manifold of saddle 𝐸₄ split the feasible region into two parts such that orbits initiating inside tend to the positive equilibrium 𝐸₃, while orbits initiating outside tend to 𝐸₃ except for the stable manifolds of 𝐸₄.

Figure 1 

The phase portraits of model (1.3). The parameters are taken as 𝑎=2,𝑏=0.25,𝑐=0.875,𝑒=1,𝑚=0.95,and𝑛=0.75,𝑓=0.8. In this case, 𝐸0=(0,0) is a source, 𝐸1=(4,0) and 𝐸4=(0.31772,1.5846) are saddle points, and 𝐸2=(0,071429) is a nodal sink, which is stable. The positive equilibrium 𝐸3=(0.96454,2.3932) is locally asymptotically stable. The dashed curve is the 𝑥-nullcline 𝑓(𝑥,𝑦)=0, and the dotted curve is the 𝑦-nullcline 𝑔(𝑥,𝑦)=0.

Theorem 2.6. The singularity 𝐸=(𝑥,(𝑚+𝑥)/𝑓)=((𝑓−𝑏𝑚𝑛)/𝑏(𝑓+𝑛),(𝑏𝑚+1)/𝑏(𝑓+𝑛)) is(a) stable if and only if 𝑒>𝑎𝑥(−𝑏+𝑐𝑓𝑚+𝑥)𝑓𝑥+𝑛(𝑚+𝑥)²; (2.20)(b) unstable if and only if 𝑒<𝑎𝑥(−𝑏+𝑐𝑓𝑚+𝑥)𝑓𝑥+𝑛(𝑚+𝑥)². (2.21)

Proof. At the point 𝐸=(𝑥,(𝑚+𝑥)/𝑓), we have 𝐹(𝑥)=𝑏𝑓+𝑛²𝐴²₁+2𝑏𝑚𝑛𝑓+𝑛𝐴₁+𝑏𝑚²𝑛²−𝑐𝑓𝑚=𝑏𝑓+𝑛²𝐴²₁+𝑏𝑚𝑛𝑛+𝑓𝐴₁>0, (2.22) as 𝐹(𝑥)>0, then det(𝐽(𝐸))>0, and the nature of singularity 𝐸=(𝑥,(𝑚+𝑥)/𝑓) depends on the trace given by tr𝐽𝐸=𝑎𝑥(−𝑏+𝑐𝑓𝑚+𝑥)𝑓𝑥+𝑛(𝑚+𝑥)²−𝑒. (2.23)
Clearly,(a) If 𝑒>𝑎𝑥(−𝑏+𝑐𝑓(𝑚+𝑥)/(𝑓𝑥+𝑛(𝑚+𝑥))²), then tr(𝐽(𝐸))<0, and 𝐸=(𝑥,(𝑚+𝑥)/𝑓) is stable.(b) If 𝑒<𝑎𝑥(−𝑏+𝑐𝑓(𝑚+𝑥)/(𝑓𝑥+𝑛(𝑚+𝑥))²), then tr(𝐽(𝐸))>0, and 𝐸=(𝑥,(𝑚+𝑥)/𝑓) is unstable.

Theorem 2.7. Singularity 𝐸^∗=(𝑥^∗,(𝑚+𝑥^∗)/𝑓) is(a) stable if and only if 𝑒>𝑎𝑥^∗−𝑏+𝑐𝑓𝑚+𝑥^∗𝑓𝑥^∗+𝑛(𝑚+𝑥^∗)²; (2.24)(b) unstable if and only if 𝑒<𝑎𝑥^∗−𝑏+𝑐𝑓𝑚+𝑥^∗𝑓𝑥^∗+𝑛(𝑚+𝑥^∗)². (2.25)

Proof. At the point 𝐸^∗=(𝑥^∗,(𝑚+𝑥^∗)/𝑓), as 𝐵₁<0,𝐶²₁=𝐴²₁−4𝐵, so 𝐹𝑥^∗=𝑏𝑓+𝑛²𝐴₁+𝐶₁²4+𝑏𝑚𝑛𝑓+𝑛𝐴₁+𝐶₁+𝑏𝑚²𝑛²=𝑏−𝑐𝑓𝑚𝑓+𝑛²𝐴₁+𝐶₁²4+𝑏𝑚𝑛𝑓+𝑛𝐶₁+𝑚𝑛+𝑓(𝑛−𝑐)>0, (2.26) as 𝐹(𝑥^∗)>0, then det(𝐽(𝐸^∗))>0, and the nature of singularity 𝐸^∗=(𝑥^∗,(𝑚+𝑥^∗)/𝑓) depends on the trace given by tr𝐽𝐸^∗=𝑎𝑥^∗−𝑏+𝑐𝑓𝑚+𝑥^∗𝑓𝑥^∗+𝑛(𝑚+𝑥^∗)²−𝑒. (2.27) Clearly,(a) if 𝑒>𝑎𝑥^∗(−𝑏+𝑐𝑓(𝑚+𝑥^∗)/(𝑓𝑥^∗+𝑛(𝑚+𝑥^∗))²), then tr(𝐽(𝐸^∗))<0, 𝐸^∗=(𝑥^∗,(𝑚+𝑥^∗)/𝑓) is stable;(b) if 𝑒<𝑎𝑥^∗(−𝑏+𝑐𝑓(𝑚+𝑥^∗)/(𝑓𝑥^∗+𝑛(𝑚+𝑥^∗))²), then tr(𝐽(𝐸^∗))>0, 𝐸^∗=(𝑥^∗,(𝑚+𝑥^∗)/𝑓) is unstable.

Theorem 2.8. Singularity 𝐸_𝑒=(𝑥_𝑒,(𝑚+𝑥_𝑒)/𝑓) is(a) a non-hyperbolic attractor node if and only if 𝑒>𝑎𝑥_𝑒−𝑏+𝑐𝑓𝑚+𝑥_𝑒𝑓𝑥_𝑒+𝑛(𝑚+𝑥_𝑒)²; (2.28)(b) a non-hyperbolic repellor node if and only if 𝑒<𝑎𝑥_𝑒−𝑏+𝑐𝑓𝑚+𝑥_𝑒𝑓𝑥_𝑒+𝑛(𝑚+𝑥_𝑒)². (2.29)

Proof. At the point (𝑥_𝑒,(𝑚+𝑥_𝑒)/𝑓), we have 𝐹𝑥_𝑒=𝑏𝑓+𝑛²𝐴²₁4+𝑏𝑚𝑛𝑓+𝑛𝐴₁+𝑏𝑚²𝑛²−𝑐𝑓𝑚=0, (2.30) then, det(𝐽(𝐸_𝑒))=0, and tr(𝐽(𝐸_𝑒))=𝑎𝑥_𝑒(−𝑏+𝑐𝑓(𝑚+𝑥_𝑒)/(𝑓𝑥_𝑒+𝑛(𝑚+𝑥_𝑒))²)−𝑒.
Hence we can conclude (2.28) and (2.29).

Theorem 2.9. If 𝑛<𝑐<𝑛(𝑏𝑚+1) and 𝑓𝑛+(𝑛−𝑐)(𝑏𝑚𝑛+𝑛−𝑐)<0 hold the boundary equilibria 𝐸₂=(0,𝑚/𝑓) of model (1.3) is globally asymptotically stable.

Proof. Since limsup_𝑡→∞𝑥(𝑡)≤1/𝑏 and liminf_𝑡→∞𝑦(𝑡)≥(𝑏𝑚𝑛+𝑛−𝑐)/𝑏𝑓𝑛, from the first equation of model (1.3), for any 𝜇>0, there exists a 𝑇₁>0, for all 𝑡≥𝑇₁, we have 𝑑𝑥𝑐𝑑𝑡=𝑎𝑥1−𝑏𝑥−𝑐𝑥/𝑦+𝑛<𝑎𝑥1−𝑏𝑥−𝑏𝑚𝑛+𝑛−𝑐−𝑏𝑓𝜇.𝑓𝑛1+𝑏𝜇+𝑛𝑏𝑚𝑛+𝑛−𝑐−𝑏𝑓𝜇 (2.31) From the arbitrariness of 𝜇>0, we can get that 𝑑𝑥𝑐()𝑑𝑡<𝑎𝑥1−𝑏𝑥−𝑏𝑚𝑛+𝑛−𝑐𝑛𝑏𝑚𝑛+𝑓+𝑛−𝑐. (2.32)
As 1−𝑐(𝑏𝑚𝑛+𝑛−𝑐)/𝑛(𝑏𝑚𝑛+𝑓+𝑛−𝑐)<0, by standard comparison arguments, it follows that limsup_𝑡→∞𝑥(𝑡)≤0, (2.33) thus, lim_𝑡→∞𝑥(𝑡)=0. (2.34)
As a result, using the second equation of model (1.3), one can easily know that lim_𝑡→∞𝑦(𝑡)=𝑚/𝑓. The proof is complete.

Figure 2 shows the dynamics of model (1.3). In this case, 𝐸₀=(0,0) is a nodal source, 𝐸₁=(4,0) is a saddle point, and 𝐸₂=(0,1.1875) is globally asymptotically stable, that is, all orbits tend to the equilibrium 𝐸₂ for any initial values.

Figure 2 

The phase portraits of model (1.3). The parameters are taken as 𝑎=2,𝑏=0.25,𝑐=0.875,𝑒=1,𝑚=0.95,and𝑛=0.7,𝑓=0.8. In this case, 𝐸0=(0,0) is a nodal source, 𝐸1=(4,0) is a saddle point. The boundary equilibria 𝐸2=(0,1.1875) is globally asymptotically stable. The dashed curve is the 𝑥-nullcline 𝑓(𝑥,𝑦)=0, and the dotted curve is the 𝑦-nullcline 𝑔(𝑥,𝑦)=0.

Theorem 2.10. Assume 𝐵₁<0, 𝐸^∗=(𝑥^∗,𝑦^∗) is globally asymptotically stable, if the following conditions hold(i)1−𝑏𝑥^∗<𝑏(𝜁+𝑛𝛿); (ii)(𝑒(𝜂+𝑛𝜍)−(𝑥^∗/𝑦^∗)(1−𝑏𝑥^∗)(𝜁+𝑚))²<4𝑎𝑒𝑓(𝑏(𝜁+𝑛𝛿)−(1−𝑏𝑥^∗))(𝜂+𝑚)(𝜂+𝑛𝜍), where 𝜂=1/𝑏,𝜍=(1+𝑏𝑚)/𝑏𝑓,𝜁=(𝑛−𝑐)/𝑏𝑛,𝛿=(𝑏𝑚𝑛+𝑛−𝑐)/𝑏𝑛𝑓.

Proof. Define a Lyapunov function: 𝑉𝑥,𝑦=^𝑥_𝑥∗𝜉−𝑥^∗𝜉𝑑𝜉+^𝑦_𝑦∗𝜂−𝑦^∗𝜂𝑑𝜂, (2.35) so, 𝑑𝑉=𝑑𝑡𝑥−𝑥^∗𝑥𝑑𝑥+𝑑𝑡𝑦−𝑦^∗𝑦𝑑𝑦𝑑𝑡=𝑎𝑥−𝑥^∗1−𝑏𝑥−𝑐𝑦𝑥+𝑛𝑦+𝑒𝑦−𝑦^∗1−𝑓𝑦𝑥+𝑚=𝑎−𝑏+𝑐𝑦^∗𝑥^∗+𝑛𝑦^∗𝑥+𝑛𝑦𝑥−𝑥^∗2+−𝑎𝑐𝑥^∗𝑥^∗+𝑛𝑦^∗+𝑥+𝑛𝑦𝑒𝑓𝑦^∗(𝑥^∗+𝑚)(𝑥+𝑚)𝑥−𝑥^∗𝑦−𝑦^∗−𝑒𝑓𝑦−𝑦^∗2𝑥+𝑚=−𝐴₂𝑥−𝑥^∗2+𝐵₂𝑥−𝑥^∗𝑦−𝑦^∗−𝐶₂𝑦−𝑦^∗2, (2.36) let 𝐴₂=𝑎𝑏−𝑐𝑦^∗𝑥^∗+𝑛𝑦^∗𝑥+𝑛𝑦=𝑎𝑏−1−𝑏𝑥^∗,𝐵𝑥+𝑛𝑦₂=−𝑎𝑐𝑥^∗𝑥^∗+𝑛𝑦^∗+𝑥+𝑛𝑦𝑒𝑓𝑦^∗(𝑥^∗+𝑚)(𝑥+𝑚)=−𝑎𝑥^∗1−𝑏𝑥^∗𝑦^∗+𝑒𝑥+𝑛𝑦,𝐶𝑥+𝑚₂=𝑒𝑓.𝑥+𝑚 (2.37) If 𝐴₂>0,thatis,1−𝑏𝑥^∗(),𝐵<𝑏𝑥+𝑛𝑦<𝑏𝜁+𝑛𝛿₂₂<4𝐴₂𝐶₂,thatis,−𝑎𝑥^∗(1−𝑏𝑥^∗)𝑦^∗+𝑒(𝑥+𝑛𝑦)𝑥+𝑚²<4𝑎𝑒𝑓𝑥+𝑚𝑏−1−𝑏𝑥^∗𝑥+𝑛𝑦hold,the𝑑𝑉/𝑑𝑡<0. (2.38) Hence, 𝑒(𝑥+𝑛𝑦)−𝑎𝑥^∗𝑦^∗(1−𝑏𝑥^∗)(𝑥+𝑚)²<4𝑎𝑒𝑓𝑏𝑥+𝑛𝑦−1−𝑏𝑥^∗(𝑥+𝑚)𝑥+𝑛𝑦, (2.39) according to Lemma 2.2, we can obtain that 𝑒(𝑥+𝑛𝑦)−𝑎𝑥^∗𝑦^∗(1−𝑏𝑥^∗)(𝑥+𝑚)²<𝑒(𝜂+𝑛𝜍)−𝑎𝑥^∗𝑦^∗1−𝑏𝑥^∗(𝜁+𝑚)²,4𝑎𝑒𝑓𝑏(𝜁+𝑛𝛿)−1−𝑏𝑥^∗𝜂+𝑚𝜂+𝑛𝜍<4𝑎𝑒𝑓𝑏𝑥+𝑛𝑦−1−𝑏𝑥^∗(𝑥+𝑚)𝑥+𝑛𝑦, (2.40) so, 𝑥𝑒(𝜂+𝑛𝜍)−^∗𝑦^∗(1−𝑏𝑥^∗)(𝜁+𝑚)²<4𝑎𝑒𝑓𝑏(𝜁+𝑛𝛿)−1−𝑏𝑥^∗𝜂+𝑚𝜂+𝑛𝜍. (2.41)
Considering 𝐴₂>0,𝐵²₂<4𝐴₂𝐶₂, we obtain 𝑑𝑉/𝑑𝑡<0. This ends the proof.

Figure 3 shows the dynamics of model (1.3) for the case of 𝐵1<0. In this case, 𝐸₀=(0,0) is a nodal source, 𝐸₁=(4,0) is a saddle point, 𝐸₂=(0,1.1875) is a saddle point, and 𝐸^∗=(1.4761,3.0326) is globally asymptotically stable, that is, all orbits tend to the equilibrium 𝐸^∗ for any initial values.

Figure 3 

The phase portraits of model (1.3). The parameters are taken as 𝑎=2,𝑏=0.25,𝑐=0.875,𝑒=1,  𝑚=0.95,𝑛=0.9,and𝑓=0.8. In this case, 𝐸1=(4,0) and 𝐸2=(0,1.1875) are saddle points, the positive interior equilibrium point 𝐸∗=(1.4761,3.0326) is globally asymptotically stable. The dashed curve is the 𝑥-nullcline 𝑓(𝑥,𝑦)=0, and the dotted curve is the 𝑦-nullcline 𝑔(𝑥,𝑦)=0.

3. The Stochastic Model

Those important and useful works on deterministic models provide a great insight into the effect of the pollution. In the real world, population dynamics is inevitably subjected to environmental noise (see e.g., [10, 11]), which is an important component in an ecosystem. May [12] pointed out the fact that due to environmental noise, the birth rates, carrying capacity, competition coefficients, and other parameters involved in the system exhibit random fluctuation to a greater or lesser extent.

In this part, we focus on the stochastic stability analysis of model (1.3).

Taking into account the effect of randomly fluctuating environment, we incorporate white noise in each equations of model (1.3). We assume that fluctuations in the environment will manifest themselves mainly as fluctuations in the growth rates of the prey population and the predator population, set 𝑎→𝑎+𝛼𝐵̇₁(𝑡),𝑒→𝑒+𝛽𝐵̇₂(𝑡), (3.1) then the stochastic version of model (1.3) is given by the following Itô type 𝑑𝑥=𝑥1−𝑏𝑥−𝑐𝑦𝑥+𝑛𝑦𝑎𝑑𝑡+𝛼𝑑𝐵₁(𝑡),𝑑𝑦=𝑦1−𝑓𝑦𝑥+𝑚𝑒𝑑𝑡+𝛽𝑑𝐵₂(𝑡), (3.2) where 𝐵_𝑖(𝑡),𝑖=(1,2) are the 1-dimensional standard Brownian motion defined on a complete probability space (Ω,ℱ,𝒫) with a filtration {ℱ_𝑡}_𝑡∈𝑅+ satisfying the usual conditions (i.e., it is right continuous and increasing while 𝐹₀ contains all 𝒫-null sets) and 𝐵̇₁(𝑡), 𝐵̇₂(𝑡) are, respectively, white noises with possible intensity 𝛼², 𝛽².

3.1. Existence of Global Positive Solutions

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

Proof. Consider the equations 𝑑𝑢=𝑎1−𝑏𝑒^𝑢−𝑐𝑒^𝑣𝑛𝑒^𝑣+𝑒^𝑢−12𝛼²1−𝑏𝑒^𝑢−𝑐𝑒^𝑣𝑛𝑒^𝑣+𝑒^𝑢2𝑑𝑡+𝛼1−𝑏𝑒^𝑢−𝑐𝑒^𝑣𝑛𝑒^𝑣+𝑒^𝑢𝑑𝐵₁(𝑡),𝑑𝑣=𝑒1−𝑓𝑒^𝑣𝑚+𝑒^𝑢−12𝛽²1−𝑓𝑒^𝑣𝑚+𝑒^𝑢2𝑑𝑡+𝛽1−𝑓𝑒^𝑣𝑚+𝑒^𝑢𝑑𝐵₂(𝑡), (3.3) on 𝑡≥0 with initial value 𝑢(0)=ln𝑥₀, 𝑣(0)=ln𝑦₀.
It is easy to see that the coefficients of model (3.3) satisfy the local Lipschitz condition, then there is a unique local solution 𝑢(𝑡),𝑣(𝑡) on [0,𝜏_𝑒) [13]. Therefore, by Itô formula, 𝑥(𝑡)=𝑒^𝑢(𝑡),𝑦(𝑡)=𝑒^𝑣(𝑡) are the unique positive local solutions to model (3.3) with initial value 𝑥₀>0,𝑦₀>0.

Lemma 3.1 only tells us that there has a unique local positive solution to model (3.2). Next, we show this solution is global, that is, 𝜏_𝑒=∞.

Theorem 3.2. Consider model (3.2), for any given initial value (𝑥₀,𝑦₀)∈𝑅²₊, there is a unique solution (𝑥(𝑡),𝑦(𝑡)) on 𝑡≥0 and the solution will remain in 𝑅²₊ with probability 1, where 𝑅²₊={𝑥∈𝑅²∣𝑥_𝑖>0,𝑖=1,2}.

Proof. For convenience of statement, we introduce some notations. Define 𝐹𝑥,𝑦=𝑐𝑦=𝑛𝑦+𝑥,𝐺𝑥,𝑦𝑓𝑦.𝑚+𝑥 (3.4) Let 𝑛₀>0 be sufficiently large for 𝑥₀ and 𝑦₀ lying within the interval [1/𝑛₀,𝑛₀]. For each integer 𝑛>𝑛₀, define the stopping times 𝜏_𝑛=inf𝑡∈0,𝜏_𝑒(𝑡)∉1∶𝑥𝑛(𝑡)∉1,𝑛or𝑦𝑛,𝑛. (3.5)
Assume inf∅=∞ (as usual ∅= the empty set). Clearly, 𝜏_𝑛 is increasing as 𝑛→∞. Let 𝜏_∞=lim_𝑛→∞𝜏_𝑛, then 𝜏_∞≤𝜏_𝑒 a.s. If we can show that 𝜏_∞=∞ a.s., then 𝜏_𝑒=∞ a.s., and 𝑁(𝑡)∈𝑅₊ a.s., for all 𝑡≥0. In other words, to complete the proof all we need to prove that 𝜏_∞=∞ a.s. If this statement is false, then there is a pair of constants 𝑇>0 and 𝜀∈(0,1) such that 𝒫{𝜏_∞≤𝑇}>𝜀. There is an integer 𝑛₁≥𝑛₀, such that 𝒫𝜏_𝑛≤𝑇≥𝜀,𝑛≥𝑛₁. (3.6) Define a 𝐶² function 𝑉∶𝑅²₊→𝑅₊ by 𝑉𝑥,𝑦=√𝑥−1−ln𝑥2+√𝑦−1−ln𝑦2. (3.7) Obviously, 𝑉(𝑥,𝑦) is nonnegative. And 𝑎√𝑑𝑉=𝑥−12𝛼1−𝑏𝑥−𝐹𝑥,𝑦𝑑𝑡+²2−√𝑥81−𝑏𝑥−𝐹𝑥,𝑦²+𝛼√𝑑𝑡𝑥−121−𝑏𝑥−𝐹𝑥,𝑦𝑑𝐵₁𝑒√𝑡+𝑦−12+𝛽1−𝐺𝑥,𝑦𝑑𝑡₂2−√𝑦81−𝐺𝑥,𝑦²𝛽√𝑑𝑡+𝑦−121−𝐺𝑥,𝑦𝑑𝐵₂𝑡≤𝑎2√𝑒𝑥+𝑏𝑥+𝐹𝑥,𝑦𝑑𝑡+2√+𝛼𝑦+𝐺𝑥,𝑦𝑑𝑡₂82𝑏𝑥^3/2−𝑏²𝑥^5/2+2√𝑥𝐹𝑥,𝑦+2+2𝑏²𝑥²+4𝑏𝑥𝐹𝑥,𝑦+2𝐹²+𝛽𝑥,𝑦𝑑𝑡₂82√𝑦𝐺𝑥,𝑦+2+2𝐺𝑥,𝑦²𝛼𝑑𝑡+2√𝑥−11−𝑏𝑥𝑎−𝐹𝑥,𝑦𝑑𝐵₁(𝑡)+𝛽2√𝑦−11−𝐺𝑥,𝑦𝑑𝐵₂(𝑡)≤𝑘₁+𝑘₂+𝛼2√𝑥−11−𝑏𝑥−𝐹𝑥,𝑦𝑑𝐵₁(𝑡)+𝛽2√𝑦−11−𝐺𝑥,𝑦𝑑𝐵₂(𝑡), (3.8) where 𝑘₁,𝑘₂ are positive numbers. Integrating both sides of the above inequality from 0 to 𝜏_𝑛∧𝑇 and then taking expectations, yields 𝐸𝑉𝑥𝜏_𝑛∧𝑇,𝑦𝜏_𝑛∧𝑇≤𝑉𝑥₀,𝑦₀+𝐾₁+𝐾₂𝑇. (3.9) Set Ω_𝑛={𝜏_𝑛≤𝑇} for 𝑛≥𝑛₁ and by (3.6), we have 𝒫(Ω_𝑛)≥𝜀. Note that for every 𝜔∈Ω_𝑛, there is some 𝑖 such that 𝑥_𝑖(𝜏_𝑛,𝜔) equals either 𝑛 or 1/𝑛 for 𝑖=1,2, hence 𝑉(𝑥(𝜏_𝑛,𝜔),𝑦(𝜏_𝑛,𝜔)) is no less than min{(√𝑛−1−ln𝑛/2),(1/𝑛−1−ln(1/𝑛)/2)}. It then follows from (3.9) that 𝑉𝑥₀,𝑦₀+𝐾₁+𝐾₂𝑇≥𝐸𝐼_Ω𝑛(𝜔)𝑉𝑥𝜏_𝑛,𝑦𝜏_𝑛≥𝜀min√𝑛−1−ln𝑛2,1𝑛−1−ln𝑛2, (3.10) where 𝐼_Ω𝑛 is the indicator function of Ω_𝑛. Let 𝑛→∞, then ∞>𝑉𝑥₀,𝑦₀+𝐾₁+𝐾₂𝑇=∞a.s. (3.11)
This completes the proof.

3.2. Stochastic Boundedness

Let us now recall the definition of stochastically ultimate boundedness [14–16].

Definition 3.3. The solution of model (3.2) is said to be stochastically ultimately bounded if for any 𝜀∈(0,1), there is a positive constant 𝜒=𝜒(𝜔), such that for any initial data (𝑥₀,𝑦₀)∈𝑅²₊, the solution of model (3.2) has the property that limsup_𝑡→∞𝑃𝑥(𝑡),𝑦(𝑡)>𝜒≤𝜀. (3.12)

Lemma 3.4. For any 𝜃∈(0,1), there is a positive constant 𝐻=𝐻(𝜃)>0, which is independent of the initial data (𝑥₀,𝑦₀)∈𝑅²₊, such that the solution of model (3.2) has the property that limsup_𝑡→∞𝐸(𝑥(𝑡),𝑦(𝑡))^𝜃≤𝐻. (3.13)

Proof. Define 𝑉𝑥,𝑦=𝑥^𝜃+𝑦^𝜃,𝑥,𝑦∈𝑅²₊. (3.14) For the sake of discussion, we rewrite the above as 𝑉𝑥,𝑦=𝑉₁+𝑉₂,𝑉₁=𝑥^𝜃,𝑉₂=𝑦^𝜃. (3.15)
By using the Itô formula, we have 𝑑𝑉𝑥,𝑦=𝐿𝑉𝑥,𝑦𝑑𝑡+𝜃𝛼𝑥^𝜃1−𝑏𝑥−𝐹𝑥,𝑦𝑑𝐵₁(𝑡)+𝜃𝛽𝑦^𝜃1−𝐺𝑥,𝑦𝑑𝐵₂(𝑡),𝐿𝑉𝑥,𝑦=𝑎𝜃𝑥^𝜃+𝜃()𝛼1−𝑏𝑥−𝐹𝑥,𝑦𝜃−1²𝑥^𝜃21−𝑏𝑥−𝐹(𝑥,𝑦)²+𝑒𝜃𝑦^𝜃+𝜃()𝛽1−𝐺𝑥,𝑦𝜃−1²𝑦^𝜃21−𝐺(𝑥,𝑦)², (3.16) as 0<𝜃<1, we have 𝐿𝑉𝑥,𝑦≤𝑎𝜃𝑥^𝜃−𝜃(1−𝜃)𝛼²𝑥^𝜃21−𝑏𝑥−𝐹(𝑥,𝑦)²+𝑒𝜃𝑦^𝜃−𝜃()𝛽1−𝐺𝑥,𝑦1−𝜃²𝑦^𝜃21−𝐺(𝑥,𝑦)²≤𝑎𝜃𝑥^𝜃−𝜃(1−𝜃)𝛼²𝑥^𝜃2𝑏²𝑥²−2𝑏𝑥−2𝑐𝑛+𝑒𝜃𝑦^𝜃−𝜃(1−𝜃)𝛽²𝑦^𝜃2−2𝑓𝑦𝑚=𝑎𝜃𝑥^𝜃−𝜃(1−𝜃)𝛼²𝑥^𝜃2𝑏²𝑥²−2𝑏𝑥−2𝑐𝑛+𝑥^𝜃+𝑒𝜃𝑦^𝜃−𝜃(1−𝜃)𝛽²𝑦^𝜃2−2𝑓𝑦𝑚+𝑦^𝜃−𝑉𝑥,𝑦≤𝐻₁−𝑉𝑥,𝑦, (3.17) here 𝐻₁ is an integer. So 𝑑𝑉𝑥,𝑦≤𝐻₁−𝑉𝑥,𝑦𝑑𝑡+𝜃𝛼𝑥^𝜃1−𝑏𝑥−𝐹𝑥,𝑦𝑑𝐵₁(𝑡). (3.18)
Now, by using the Itô formula again, we have 𝑑𝑒^𝑡𝑉𝑥,𝑦=𝑒^𝑡𝑉𝑥,𝑦𝑑𝑡+𝑑𝑉𝑥,𝑦≤𝑒^𝑡𝐻₁𝑑𝑡+𝑒^𝑡𝜃𝛼𝑥^𝜃1−𝑏𝑥−𝐹𝑥,𝑦𝑑𝐵₁(𝑡)+𝑒^𝑡𝜃𝛽𝑦^𝜃1−𝐺𝑥,𝑦𝑑𝐵₂(𝑡). (3.19)
According to the above, we can easily get 𝑒^𝑡𝐸𝑉𝑥,𝑦≤𝑉𝑥₀,𝑦₀+𝐻₁𝑒^𝑡−1,limsup_𝑡→∞𝐸𝑉𝑥(𝑡),𝑦(𝑡)≤𝐻. (3.20)
On the other hand, (𝑥,𝑦)^𝜃≤2^𝜃/2max𝑥^𝜃,𝑦^𝜃≤2^𝜃/2𝑉𝑥,𝑦, (3.21) so, we can obtain limsup_𝑡→∞𝐸𝑥(𝑡),𝑦(𝑡)^𝜃≤2^𝜃/2,limsup_𝑡→∞𝐸𝑉𝑥(𝑡),𝑦(𝑡)≤2^𝜃/2𝐻. (3.22)
Set 𝐾(𝜃)=2^𝜃/2𝐻, this completes the proof.

Theorem 3.5. The solutions of model (3.2) is stochastically ultimately bounded.

Proof. By Lemma 3.4, let 𝜃=1/2, there is a 𝐾₁>0 such that limsup_𝑡→∞𝐸𝑥(𝑡),𝑦(𝑡)^1/2≤𝐾₁. (3.23)
Now, for any 𝜀>0, let 𝜒=𝐾²₁/𝜀². Then by Chebyshevs inequality 𝑃𝑥(𝑡)(𝑡≤𝐸𝑥(𝑡)(𝑡),𝑦)>𝜒,𝑦^1/2√𝜒. (3.24)
Hence, limsup_𝑡→∞𝑃𝑥(𝑡),𝑦(𝑡)>𝜒≤𝜀. (3.25)
This ends the proof.

Next, we study the asymptotic properties of the moment solutions of model (3.2).

Theorem 3.6. For any given 𝜃∈(0,1), there is a 𝐾=𝐾(𝜃)>0, let (𝑥(𝑡),𝑦(𝑡)) be the solution of model (3.2) with any initial value (𝑥₀,𝑦₀), then limsup_𝑡→∞1𝑡^𝑡₀𝐸𝑥^2+𝜃(𝑠)+𝑦^2+𝜃(𝑠)𝑑𝑠≤𝐾. (3.26)

Proof. Set 𝑉(𝑥,𝑦)∶𝑅²₊→𝑅₊, from Lemma 3.4, we have 𝐿𝑉𝑥,𝑦≤𝑎𝜃𝑥^𝜃−𝜃(1−𝜃)𝛼²𝑥^𝜃2𝑏²𝑥²−2𝑏𝑥−2𝑐𝑛+𝑒𝜃𝑦^𝜃−𝜃(1−𝜃)𝛽²𝑦^𝜃2−2𝑓𝑦𝑚, (3.27) so, we get 𝜃()𝛼𝐿𝑉+1−𝜃²𝑏²𝑥^𝜃+24+𝜃(1−𝜃)𝛽²𝑓𝑦^𝜃+24𝑚≤𝑎𝜃𝑥^𝜃−𝜃(1−𝜃)𝛼²𝑏²𝑥^𝜃+24+𝜃(1−𝜃)𝛼²𝑥^𝜃𝑐𝑏𝑥+𝑛+𝑒𝜃𝑦^𝜃+𝜃(1−𝜃)𝛽²𝑓𝑦^𝜃+1𝑚+𝜃(1−𝜃)𝛽²𝑓𝑦^𝜃+24𝑚≤𝑀, (3.28) here 𝑀 is a positive number, so 𝜃()𝛾𝐿𝑉+1−𝜃²𝜎²𝑥^𝜃+2+𝑦^𝜃+24≤𝜃(1−𝜃)𝛼²𝑏²𝑥^𝜃+24+𝜃(1−𝜃)𝛽²𝑓𝑦^𝜃+24𝑚≤𝑀. (3.29)
Let 𝛾²𝜎²=min{𝛼²𝑏²,𝛽²𝑓/𝑚}, we can get +𝜃()𝛾𝐸𝑉𝑥,𝑦1−𝜃²𝜎²4^𝑡₀𝐸𝑥^2+𝜃(𝑠)+𝑦^2+𝜃(𝑠𝑥)𝑑𝑠≤𝑉₀+𝑦₀+𝑀𝑡. (3.30)
That is, ^𝑡₀𝐸𝑥^2+𝜃(𝑠)+𝑦^2+𝜃(𝑠4𝑉𝑥)𝑑𝑠≤⁰+𝑦₀+𝑀𝑡𝜃(1−𝜃)𝛾²𝜎². (3.31)
Set 𝐾=4𝑀/𝜃(1−𝜃)𝛾²𝜎², so limsup_𝑡→∞1𝑡^𝑡₀𝐸𝑥^2+𝜃(𝑠)+𝑦^2+𝜃(𝑠)𝑑𝑠≤𝐾. (3.32)

3.3. Stochastic Asymptotic Stability

Note that a solution of model (1.3) is also a solution of model (3.2), so, in the following, we will focus on stochastic asymptotic stability of the positive equilibria of model (3.2). As an example, we only give the proof of the unique positive equilibrium 𝐸^∗=(𝑥^∗,𝑦^∗) of model (3.2).

Theorem 3.7. Let 𝐴₃=−𝑎𝑏+𝑎𝑐𝑦^∗𝑛𝑦^∗+𝑥^∗+𝛼²𝑏²𝑥^∗2+𝛼²𝑐²𝑥^∗𝑦^∗22𝑛𝑦^∗+𝑥^∗2+𝛽²𝑓²𝑦^∗32(𝑚+𝑥^∗)²,𝐵₃=𝛼²𝑐²𝑥^∗32𝑛𝑦^∗+𝑥^∗2+𝛽²𝑓𝑚+𝑥^∗2𝑦^∗2(𝑚+𝑥^∗)²,𝐶₃=𝛼²𝑏𝑐𝑥^∗2𝑛𝑦^∗+𝑥^∗+𝛼²𝑐²𝑦^∗𝑥^∗2𝑛𝑦^∗+𝑥^∗2+𝑓𝑦^∗(𝑚+𝑥^∗)+𝛽²𝑓𝑚+𝑥^∗𝑓𝑦^∗2(𝑚+𝑥^∗)², (3.33) if (𝑏𝑚+1)𝑛²+))]𝐴(1−𝑏𝑓−𝑐𝑛−𝑐𝑓≥0,[(𝑚−1𝑏+1𝑛−𝑐≥0,₃<0,𝐶²₃−4𝐴₃𝐵₃<0, (3.34) when 𝐵₁<0, then the equilibrium position 𝐸^∗=(𝑥^∗,𝑦^∗) of model (3.2) is stochastically asymptotically stable in the large, that is, for any initial data (𝑥(0),𝑦(0)), the solution of model (3.2) has the property that lim_𝑡→∞(𝑥(𝑡))=𝑥^∗,lim_𝑡→∞𝑦(𝑡)=𝑦^∗,𝑎.𝑠. (3.35)

Proof. From the theory of stability of stochastic differential equations [13], we only need to establish a Lyapunov function 𝑉(𝑧) satisfying 𝐿𝑉(𝑧)≤0 and the identity holds if and only if 𝑧=𝑧^∗, where 𝑧=𝑧(𝑡) is the solution of the stochastic differential equation 𝑑𝑧=𝑓(𝑧(𝑡),𝑡)𝑑𝑡+𝑔(𝑧(𝑡),𝑡)𝑑𝐵(𝑡), (3.36)𝑧^∗ is the equilibrium position of model (3.36), and 𝐿𝑉(𝑧,𝑡)=𝑉_𝑧()𝑓()+1𝑧,𝑡𝑧,𝑡2tr𝑔^𝑇(𝑧,𝑡)𝑉_𝑧𝑧(𝑧,𝑡)𝑔(𝑧,𝑡). (3.37) Define Lyapunov functions 𝑉₁(𝑥)=^𝑥_𝑥∗𝜉−𝑥^∗𝜉𝑑𝜉,𝑉₂(𝑥)=^𝑦_𝑦∗𝜂−𝑦^∗𝜂𝑑𝜂, (3.38) the nonnegativity of this function can be observed from 𝑢−1−ln𝑢≥0 on 𝑢>0. Applying Itô formula leads to 𝐿𝑉₁(𝑥)=𝑎𝑥−𝑥^∗+𝛼1−𝑏𝑥−𝐹𝑥,𝑦²𝑥^∗21−𝑏𝑥−𝐹(𝑥,𝑦)²=𝑎𝑥−𝑥^∗𝑏𝑥^∗−𝑏𝑥+𝐹𝑥^∗,𝑦^∗+𝛼−𝐹𝑥,𝑦²𝑥^∗2𝑏𝑥^∗−𝑏𝑥+𝐹𝑥^∗,𝑦^∗−𝐹𝑥,𝑦²=−𝑎𝑏+𝑎𝑐𝑦^∗𝑛𝑦^∗+𝑥^∗+𝛼𝑛𝑦+𝑥²𝑏²𝑥^∗2−𝑏𝑐𝛼²𝑥^∗𝑦^∗𝑛𝑦^∗+𝑥^∗+𝑐𝑛𝑦+𝑥₂𝛼²𝑥^∗𝑦^∗22𝑛𝑦^∗+𝑥^∗2𝑛𝑦+𝑥²𝑥−𝑥^∗2+𝑐²𝛼²𝑥^∗32𝑛𝑦^∗+𝑥^∗2𝑛𝑦+𝑥²𝑦−𝑦^∗2−𝑎𝑐𝑥^∗𝑛𝑦^∗+𝑥^∗𝑛𝑦+𝑥𝑥−𝑥^∗𝑦−𝑦^∗+𝑏𝑐𝛼²𝑥^∗2𝑛𝑦^∗+𝑥^∗𝑛𝑦+𝑥𝑥−𝑥^∗𝑦−𝑦^∗−𝑐²𝛼²𝑦^∗𝑥^∗2𝑛𝑦^∗+𝑥^∗2𝑛𝑦+𝑥²𝑥−𝑥^∗𝑦−𝑦^∗, (3.39) when 𝑛𝑦+𝑥≥1, according to Lemma 2.2, we have (𝑏𝑚+1)𝑛²+[(1−𝑏)𝑓−𝑐]𝑛−𝑐𝑓≥0, so, 𝐿𝑉₁(𝑥)≤−𝑎𝑏+𝑎𝑐𝑦^∗𝑛𝑦^∗+𝑥^∗+𝛼²𝑏²𝑥^∗2+𝑐²𝛼²𝑥^∗𝑦^∗22𝑛𝑦^∗+𝑥^∗2𝑥−𝑥^∗2+𝑐²𝛼²𝑥^∗32𝑛𝑦^∗+𝑥^∗2𝑦−𝑦^∗2+𝑏𝑐𝛼²𝑥^∗2𝑛𝑦^∗+𝑥^∗+𝑐²𝛼²𝑦^∗𝑥^∗2𝑛𝑦^∗+𝑥^∗2𝑥−𝑥^∗𝑦−𝑦^∗. (3.40)
Similarly, by using the Itô formula, we can obtain 𝐿𝑉₂𝑦=𝑒𝑦−𝑦^∗+𝛽1−𝐺𝑥,𝑦²𝑦^∗21−𝐺𝑥,𝑦², (3.41) when 𝑚+𝑥≥1, according to Lemma 2.2, we have [(𝑚−1)𝑏+1]𝑛−𝑐≥0, so 𝐿𝑉₂𝑦≤𝛽²𝑓𝑚+𝑥^∗2𝑦^∗2(𝑚+𝑥^∗)²𝑦−𝑦^∗2+𝛽²𝑓²𝑦^∗32(𝑚+𝑥^∗)²𝑥−𝑥^∗2+𝑓𝑦^∗(𝑚+𝑥^∗)+𝛽²𝑓𝑚+𝑥^∗𝑓𝑦^∗2(𝑚+𝑥^∗)²𝑥−𝑥^∗𝑦−𝑦^∗, (3.42) hence, 𝐿𝑉𝑥,𝑦=𝐿𝑉₁(𝑥)+𝐿𝑉₂𝑦≤−𝑎𝑏+𝑎𝑐𝑦^∗𝑛𝑦^∗+𝑥^∗+𝛼²𝑏²𝑥^∗2+𝑐²𝛼²𝑥^∗𝑦^∗22𝑛𝑦^∗+𝑥^∗2+𝛽²𝑓²𝑦^∗32(𝑚+𝑥^∗)²𝑥−𝑥^∗2+𝑐²𝛼²𝑥^∗32𝑛𝑦^∗+𝑥^∗2+𝛽²𝑓𝑚+𝑥^∗2𝑦^∗2(𝑚+𝑥^∗)²𝑦−𝑦^∗2+𝑏𝑐𝛼²𝑥^∗2𝑎²𝑛𝑦^∗+𝑥^∗+𝑐²𝛼²𝑦^∗𝑥^∗2𝑛𝑦^∗+𝑥^∗2+𝑓𝑦^∗(𝑚+𝑥^∗)+𝛽²𝑓𝑚+𝑥^∗𝑓𝑦^∗2(𝑚+𝑥^∗)²𝑥−𝑥^∗𝑦−𝑦^∗. (3.43)
So, 𝐿𝑉𝑥,𝑦≤𝐴₃𝑥−𝑥^∗2+𝐵₃𝑦−𝑦^∗2+𝐶₃𝑥−𝑥^∗𝑦−𝑦^∗. (3.44)
Let 𝑧−𝑧^∗=𝑥−𝑥^∗,𝑦−𝑦^∗𝑇, (3.45) then we get ≤1𝐿𝑉𝑥,𝑦2𝑧−𝑧^∗𝑇2𝐴₃𝐶₃𝐶₃2𝐵₃𝑧−𝑧^∗=𝐴₃𝑥−𝑥^∗2+𝐶₃𝑥−𝑥^∗𝑦−𝑦^∗+𝐵₃𝑦−𝑦^∗2. (3.46)
Clearly, if (3.34) hold, then the above inequality implies 𝐿𝑉(𝑥,𝑦)<0 along all trajectories in the first quadrant except (𝑥^∗,𝑦^∗). Then the desired assertion (3.35) follows immediately. This completes the proof.

4. Conclusions and Remarks

In this paper, we consider a modified stochastic Leslie-Gower predator-prey model. The value of this study lies in two aspects. First, it presents the analysis of stability for the equilibria of model (1.3). Second, it verifies some relevant properties of the stochastic model (3.2) with white noise, which shows that the existence of global positive solutions, stochastic boundedness, and stochastic asymptotic stability.

Next, we give some numerical examples to illustrate the dynamical behaviors of model (3.2) by using the method mentioned in [17].

In Figure 4(a), we choose 𝛼=𝛽=0, that is, without noise, we observe that the positive equilibrium 𝐸^∗=(1.4761,3.0326) is globally stable. In Figure 4(b), with noise densities 𝛼=0.8,𝛽=0.4, starting with a homogeneous state 𝐸^∗=(1.4761,3.0326), the random white noise leads to a slight oscillations, and the later random noise makes the oscillations decay, ending with the time-independent stability. Comparing Figure 4(a) with Figure 4(b), one can realize that, if the white noise is not strong, the stochastic perturbation does not cause sharp changes of the dynamics of the system. However, in Figure 5, we choose that 𝛼=2 and 𝛽=4, which violates condition (3.34), we find that the stochastic model (3.2) is not permanent. This shows that strong white noise might make a permanent system be nonpersistent.

(a)

(b)

Figure 4 

The solution of the stochastic model (3.2). The parameters are taken as 𝑎=2,𝑏=0.25,𝑐=0.875,𝑒=1,𝑚=0.95,𝑛=0.9,and𝑓=0.8. (a) 𝛼=𝛽=0; (b) 𝛼=0.8and𝛽=0.4.

Figure 5 

The solution of the stochastic model (3.2). The parameters are taken as 𝛼=2,𝛽=4, 𝑎=2,𝑏=0.25,𝑐=0.875,𝑒=1,𝑚=0.95,𝑛=0.9,and𝑓=0.8.

We note that, when the noise is not large, the stochastic model preserves the property of the global stability, that is to say, when the noise is not sufficiently large, the populations may be stochastic permanence and stochastic persistent in mean. In this case, we can ignore the noise and use the deterministic model to describe the population dynamics. But, when the noise is sufficiently large, the noise can force the population to become extinct. In this case, we cannot ignore the effect of the noise. That’s to say, in the case of sufficiently large noise, we cannot use deterministic model but stochastic model to describe the population dynamics. Our complete analysis of the model will give new suggestions to the studies of the population dynamics.

Acknowledgments

The authors would like to thank the anonymous referee for very helpful suggestions and comments which led to improvements of our original manuscript. This research was supported by the National Science Foundation of China (11201344 & 60970065), the National Basic Research Program of China (2012CB426510), and Zhejiang Provincial Natural Science Foundation (R1110261 & LY12A01014).