# The Complex Dynamics of a Stochastic Predator-Prey Model

Authors: Xixi Wang; Huilin Huang; Yongli Cai; Weiming Wang
Journal: Abstract and Applied Analysis
 (2012)

---

## 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.


---

## Body

## 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)dxdt=ax(1-bx)-p(xy)y,dydt=ey(1-fyx),
where x=x(t), y=y(t) stand for the population (the density) of the prey and the predator at time t, respectively, and p(x/y) is the predator functional response to prey. And we assumed that the prey grows logistically with growth rate a and carrying capacity 1/b in the absence of predation. The predator consumes the prey according to the functional response p(x/y) and grow logistically with growth rate e and carrying capacity x/f proportional to the population size of prey (or prey abundance). The parameter f is a measure of the food quality that the prey provides for conversion into predator birth. The term fy/x of this equation is called the Leslie-Gower term.On the other hand, the predatory 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 x, is not in abundance. In this situation, a positive constant m can be added to the denominator, m measures the extent to which the environment provides protection to the predator [8, 9], and the second equation of model (1.1) becomes
(1.2)dydt=ey(1-fym+x).Based on the above discussions, in the paper, we will focus on the following ratio-dependent Leslie-Gower model:(1.3)dxdt=ax(1-bx-cyx+ny)≜f(x,y),dydt=ey(1-fym+x)≜g(x,y).The rest of the paper is organized as follows. In Section2, 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
(2.1)limsupt→∞x(t)≤1b.
Hence, there exists a T>0 such that, for t>T, x≤1/b.From the second equation of model (1.3), we see that, for t>T, we have dy/dt≤ey(1-bfy/(mb+1)). A standard comparison argument shows that
(2.2)limsupt→∞y(t)≤(1+bm)bf.Thus, we have the following conclusion.Lemma 2.1.
Model (1.3) is dissipative.Lemma 2.2.
Ifc<n, then model (1.3) is permanent.Proof.
Ifc<n, from the first equation of model (1.3), we have dx/dt≥ax(1-bx-c/n). Therefore, by standard comparison argument, we have
(2.3)liminft→∞x(t)≥n-cbn.
Hence, for anyɛ>0 and large t, x(t)>(n-c)/bn-ɛ, and
(2.4)dydt≥ey(bmn+nc-bnɛ-bfnybmn+nc-bnɛ).
From the arbitrariness ofɛ>0, we can get that
(2.5)liminft→∞y(t)≥bmn+n-cbfn.


### 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 E0=(0,0), E1=(1/b,0), E2=(0,m/f). And the positive equilibria (x,y) satisfies the equations
(2.6)ax(1-bx-cyx+ny)=0,ey(1-fym+x)=0,
which yields
(2.7)b(f+n)x2-(f-c-bmn+n)x+(c-n)m=0.For simplicity, we define(2.8)A1=f-c-bmn+nb(f+n),B1=(c-n)mb(f+n),C1=(f-c-bmn+n)2b2(f+n)2-4(c-n)mb(f+n)=A12-4B1,
then (2.7) can be rewritten as
(2.9)x2-A1x+B1=0.Lemma 2.3 (existence of equilibria).
(a) Suppose A1>0 and B1>0, in (2.9), then one has(a1)
IfC12>0, then it has two positive roots given by
(2.10)xe1=A1+C12,xe2=A1-C12.Therefore, model (1.3) has two positive equilibria E3=(xe1,ye1)=((A1+C1)/2,(2m + A1+C1)/2f) and E4=(xe2,ye2)=((A1-C1)/2,(2m+A1-C1)/2f).(a2)
IfC12=0, (2.9) has a unique positive root of multiplicity-2 given by xe=m(c-n)/b(f+n). Thus, model (1.3) has a unique positive equilibrium
(2.11)Ee=(m(c-n)b(f+n),bm(f+n)+bm(f+n)(c-n)bf(f+n)).(a3)
IfC12<0, there are no positive roots that exist and then, model (1.3) has no positive equilibrium.
(b)
B
1
=
0
,
A
1
>
0, (2.9) has a unique positive root x¯=(f-bmn)/b(f+n), model (1.3) has a unique positive equilibrium E¯=(x¯,y¯)=((f-bmn)/b(f+n),(bm+1)/b(f+n)).
(c)
B
1
<
0, (2.9) has a unique positive root x*=(A1+C1)/2, model (1.3) has a unique positive equilibrium E*=(x*,y*)=((A1+C1)/2,(2m+A1+C1)/2f).Next, we discuss the local stabilities of these equilibria. Easy to obtain the following results:Lemma 2.4.
(
i
)  E1=(1/b,0) is a saddle point;(ii)
ifc>n, E2=(0,m/f) is a stable node point. If c≤n, E2=(0,m/f) is a saddle point.Following, letE=(x,y) be an arbitrary positive equilibrium. The Jacobian matrix for E=(x,y) is given by
(2.12)J(E)=(ax(-b+cy(x+ny)2)ax(-cx+ny+cny(x+ny)2)ef-e).Then we can get the following:(2.13)tr(J(E))=ax(-b+cy(x+ny)2)-e,det(J(E))=aex[b(f+n)2x2+2bmn(f+n)x-cfm+bm2n2](fx+mn+nx)2.So the sign ofdet(J(E)) is determined by
(2.14)F(x)≜b(f+n)2x2+2bmn(f+n)x-cfm+bm2n2.Theorem 2.5.
For model (1.3), the stabilities of E3 and E4 are as follows.(a)
The positive equilibriumE3=(xe1,(m+xe1)/f)=((A1+C1)/2,(2m+A1+C1)/2f) is
(a1)
stable if and only if(2.15)e>axe1(-b+cf(m+xe1)(fxe1+n(m+xe1))2);(a2)
unstable if and only if(2.16)e<axe1(-b+cf(m+xe1)(fxe1+n(m+xe1))2).(b)
The positive equilibriumE4=(xe2,(m+xe2)/f)=((A1-C1)/2,(2m+A1-C1)/2f) is a saddle point.Proof.
(a) At the pointE3, we have the following:
(2.17)F(xe1)=b(f+n)2(A1+C1)24+bmn(f+n)(A1+C1)+bm2n2-cfm=b(f+n)2C122+C1(bmn(f+n)+b(f+n)2A12)>0,thus,det(J(E3))>0, and the nature of singularity E3 depends on the trace given by
(2.18)tr(J(E3))=axe1(-b+cf(m+xe1)(fxe1+n(m+xe1))2)-e.Clearly,
(a1)
Ife>axe1(-b+cf(m+xe1)/(fxe1+n(m+xe1))2), then tr(J(E3))<0, and E3=(xe1,(m+xe1)/f) is stable.(a2)
Ife<axe1(-b+cf(m+xe1)/(fxe1+n(m+xe1))2), then tr(J(E3))>0, and E3=(xe1,(m+xe1)/f) is unstable.(b)
At the pointE4, we have the following:
(2.19)F(xe2)=b(f+n)2(A1-C1)24+bmn(f+n)(A1-C1)+bm2n2-cfm=b(f+n)22C1(A1-C1)-bmn(f+n)C1-2fm(f+n)B1<0,thendet(J(E4))<0, and E4 is a saddle point.Figure1 shows the dynamics of model (1.3). In this case, E0=(0,0) is a nodal source, E1=(4,0) is a saddle point, E2=(0,071429) is a nodal sink, E4=(0.31772,1.5846) is a saddle point, and E3=(0.96454,2.3932) is locally asymptotically stable. There exists a separatrix curve determined by the stable manifold of E4, which divides the behavior of trajectories, that is, the stable manifold of saddle E4 split the feasible region into two parts such that orbits initiating inside tend to the positive equilibrium E3, while orbits initiating outside tend to E3 except for the stable manifolds of E4.Figure 1
The phase portraits of model (1.3). The parameters are taken as a=2,b=0.25,c=0.875,e=1,m=0.95,andn=0.75,f=0.8. In this case, E0=(0,0) is a source, E1=(4,0) and E4=(0.31772,1.5846) are saddle points, and E2=(0,071429) is a nodal sink, which is stable. The positive equilibrium E3=(0.96454,2.3932) is locally asymptotically stable. The dashed curve is the x-nullcline f(x,y)=0, and the dotted curve is the y-nullcline g(x,y)=0.Theorem 2.6.
The singularityE¯=(x¯,(m+x¯)/f)=((f-bmn)/b(f+n),(bm+1)/b(f+n)) is(a)
stable if and only if(2.20)e>ax¯(-b+cf(m+x¯)(fx¯+n(m+x¯))2);(b)
unstable if and only if(2.21)e<ax¯(-b+cf(m+x¯)(fx¯+n(m+x¯))2).Proof.
At the pointE¯=(x¯,(m+x¯)/f), we have
(2.22)F(x¯)=b(f+n)2A12+2bmn(f+n)A1+bm2n2-cfm=b(f+n)2A12+bmn(n+f)A1>0,
as F(x¯)>0, then det(J(E¯))>0, and the nature of singularity E¯=(x¯,(m+x¯)/f) depends on the trace given by
(2.23)tr(J(E¯))=ax¯(-b+cf(m+x¯)(fx¯+n(m+x¯))2)-e.
Clearly,(a)
Ife>ax¯(-b+cf(m+x¯)/(fx¯+n(m+x¯))2), then tr(J(E¯))<0, and E¯=(x¯,(m+x¯)/f) is stable.(b)
Ife<ax¯(-b+cf(m+x¯)/(fx¯+n(m+x¯))2), then tr(J(E¯))>0, and E¯=(x¯,(m+x¯)/f) is unstable.Theorem 2.7.
SingularityE*=(x*,(m+x*)/f) is(a)
stable if and only if(2.24)e>ax*(-b+cf(m+x*)(fx*+n(m+x*))2);(b)
unstable if and only if(2.25)e<ax*(-b+cf(m+x*)(fx*+n(m+x*))2).Proof.
At the pointE*=(x*,(m+x*)/f), as B1<0,C12=A12-4B, so
(2.26)F(x*)=b(f+n)2(A1+C1)24+bmn(f+n)(A1+C1)+bm2n2-cfm=b(f+n)2(A1+C1)24+bmn(f+n)C1+m(n+f)(n-c)>0,
as F(x*)>0, then det(J(E*))>0, and the nature of singularity E*=(x*,(m+x*)/f) depends on the trace given by
(2.27)tr(J(E*))=ax*(-b+cf(m+x*)(fx*+n(m+x*))2)-e.
Clearly,(a)
ife>ax*(-b+cf(m+x*)/(fx*+n(m+x*))2), then tr(J(E*))<0, E*=(x*,(m+x*)/f) is stable;(b)
ife<ax*(-b+cf(m+x*)/(fx*+n(m+x*))2), then tr(J(E*))>0, E*=(x*,(m+x*)/f) is unstable.Theorem 2.8.
SingularityEe=(xe,(m+xe)/f) is(a)
a non-hyperbolic attractor node if and only if(2.28)e>axe(-b+cf(m+xe)(fxe+n(m+xe))2);(b)
a non-hyperbolic repellor node if and only if(2.29)e<axe(-b+cf(m+xe)(fxe+n(m+xe))2).Proof.
At the point(xe,(m+xe)/f), we have
(2.30)F(xe)=b(f+n)2A124+bmn(f+n)A1+bm2n2-cfm=0,
then, det(J(Ee))=0, and tr(J(Ee))=axe(-b+cf(m+xe)/(fxe+n(m+xe))2)-e.
Hence we can conclude (2.28) and (2.29).Theorem 2.9.
Ifn<c<n(bm+1) and fn+(n-c)(bmn+n-c)<0 hold the boundary equilibria E2=(0,m/f) of model (1.3) is globally asymptotically stable.Proof.
Sincelimsupt→∞x(t)≤1/b and liminft→∞y(t)≥(bmn+n-c)/bfn, from the first equation of model (1.3), for any μ>0, there exists a T1>0, for all t≥T1, we have
(2.31)dxdt=ax(1-bx-cx/y+n)<ax(1-bx-c(bmn+n-c-bfμ)fn(1+bμ)+n(bmn+n-c-bfμ)).
From the arbitrariness of μ>0, we can get that
(2.32)dxdt<ax(1-bx-c(bmn+n-c)n(bmn+f+n-c)).
As1-c(bmn+n-c)/n(bmn+f+n-c)<0, by standard comparison arguments, it follows that
(2.33)limsupt→∞x(t)≤0,
thus,
(2.34)limt→∞x(t)=0.
As a result, using the second equation of model (1.3), one can easily know that limt→∞y(t)=m/f. The proof is complete.Figure2 shows the dynamics of model (1.3). In this case, E0=(0,0) is a nodal source, E1=(4,0) is a saddle point, and E2=(0,1.1875) is globally asymptotically stable, that is, all orbits tend to the equilibrium E2 for any initial values.Figure 2
The phase portraits of model (1.3). The parameters are taken as a=2,b=0.25,c=0.875,e=1,m=0.95,andn=0.7,f=0.8. In this case, E0=(0,0) is a nodal source, E1=(4,0) is a saddle point. The boundary equilibria E2=(0,1.1875) is globally asymptotically stable. The dashed curve is the x-nullcline f(x,y)=0, and the dotted curve is the y-nullcline g(x,y)=0.Theorem 2.10.
AssumeB1<0, E*=(x*,y*) is globally asymptotically stable, if the following conditions hold(i)
1-bx*<b(ζ+nδ);(ii)
(e(η+nς)-(x*/y*)(1-bx*)(ζ+m))2<4aef(b(ζ+nδ)-(1-bx*))(η+m)(η+nς), where η=1/b,ς=(1+bm)/bf,ζ=(n-c)/bn,δ=(bmn+n-c)/bnf.Proof.
Define a Lyapunov function:(2.35)V(x,y)=∫x*x‍ξ-x*ξdξ+∫y*y‍η-y*ηdη,
so,
(2.36)dVdt=x-x*xdxdt+(y-y*)ydydt=a(x-x*)(1-bx-cyx+ny)+e(y-y*)(1-fyx+m)=a(-b+cy*(x*+ny*)(x+ny))(x-x*)2+(-acx*(x*+ny*)(x+ny)+efy*(x*+m)(x+m))(x-x*)(y-y*)-ef(y-y*)2x+m=-A2(x-x*)2+B2(x-x*)(y-y*)-C2(y-y*)2,
let
(2.37)A2=a(b-cy*(x*+ny*)(x+ny))=a(b-1-bx*x+ny),B2=-acx*(x*+ny*)(x+ny)+efy*(x*+m)(x+m)=-ax*(1-bx*)y*(x+ny)+ex+m,C2=efx+m.
If
(2.38)A2>0,thatis,1-bx*<b(x+ny)<b(ζ+nδ),B22<4A2C2,thatis,(-ax*(1-bx*)y*(x+ny)+ex+m)2<4aefx+m(b-1-bx*x+ny)hold,thedV/dt<0.
Hence,
(2.39)(e(x+ny)-ax*y*(1-bx*)(x+m))2<4aef(b(x+ny)-(1-bx*))(x+m)(x+ny),
according to Lemma 2.2, we can obtain that
(2.40)(e(x+ny)-ax*y*(1-bx*)(x+m))2<(e(η+nς)-ax*y*(1-bx*)(ζ+m))2,4aef(b(ζ+nδ)-(1-bx*))(η+m)(η+nς)<4aef(b(x+ny)-(1-bx*))(x+m)(x+ny),
so,
(2.41)(e(η+nς)-x*y*(1-bx*)(ζ+m))2<4aef(b(ζ+nδ)-(1-bx*))(η+m)(η+nς).
ConsideringA2>0,B22<4A2C2, we obtain dV/dt<0. This ends the proof.Figure3 shows the dynamics of model (1.3) for the case of B1<0. In this case, E0=(0,0) is a nodal source, E1=(4,0) is a saddle point, E2=(0,1.1875) is a saddle point, and E*=(1.4761,3.0326) is globally asymptotically stable, that is, all orbits tend to the equilibrium E* for any initial values.Figure 3
The phase portraits of model (1.3). The parameters are taken as a=2,b=0.25,c=0.875,e=1,  m=0.95,n=0.9,andf=0.8. In this case, E1=(4,0) and E2=(0,1.1875) are saddle points, the positive interior equilibrium point E*=(1.4761,3.0326) is globally asymptotically stable. The dashed curve is the x-nullcline f(x,y)=0, and the dotted curve is the y-nullcline g(x,y)=0.


## 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
(2.1)limsupt→∞x(t)≤1b.
Hence, there exists a T>0 such that, for t>T, x≤1/b.From the second equation of model (1.3), we see that, for t>T, we have dy/dt≤ey(1-bfy/(mb+1)). A standard comparison argument shows that
(2.2)limsupt→∞y(t)≤(1+bm)bf.Thus, we have the following conclusion.Lemma 2.1.
Model (1.3) is dissipative.Lemma 2.2.
Ifc<n, then model (1.3) is permanent.Proof.
Ifc<n, from the first equation of model (1.3), we have dx/dt≥ax(1-bx-c/n). Therefore, by standard comparison argument, we have
(2.3)liminft→∞x(t)≥n-cbn.
Hence, for anyɛ>0 and large t, x(t)>(n-c)/bn-ɛ, and
(2.4)dydt≥ey(bmn+nc-bnɛ-bfnybmn+nc-bnɛ).
From the arbitrariness ofɛ>0, we can get that
(2.5)liminft→∞y(t)≥bmn+n-cbfn.


## 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 E0=(0,0), E1=(1/b,0), E2=(0,m/f). And the positive equilibria (x,y) satisfies the equations
(2.6)ax(1-bx-cyx+ny)=0,ey(1-fym+x)=0,
which yields
(2.7)b(f+n)x2-(f-c-bmn+n)x+(c-n)m=0.For simplicity, we define(2.8)A1=f-c-bmn+nb(f+n),B1=(c-n)mb(f+n),C1=(f-c-bmn+n)2b2(f+n)2-4(c-n)mb(f+n)=A12-4B1,
then (2.7) can be rewritten as
(2.9)x2-A1x+B1=0.Lemma 2.3 (existence of equilibria).
(a) Suppose A1>0 and B1>0, in (2.9), then one has(a1)
IfC12>0, then it has two positive roots given by
(2.10)xe1=A1+C12,xe2=A1-C12.Therefore, model (1.3) has two positive equilibria E3=(xe1,ye1)=((A1+C1)/2,(2m + A1+C1)/2f) and E4=(xe2,ye2)=((A1-C1)/2,(2m+A1-C1)/2f).(a2)
IfC12=0, (2.9) has a unique positive root of multiplicity-2 given by xe=m(c-n)/b(f+n). Thus, model (1.3) has a unique positive equilibrium
(2.11)Ee=(m(c-n)b(f+n),bm(f+n)+bm(f+n)(c-n)bf(f+n)).(a3)
IfC12<0, there are no positive roots that exist and then, model (1.3) has no positive equilibrium.
(b)
B
1
=
0
,
A
1
>
0, (2.9) has a unique positive root x¯=(f-bmn)/b(f+n), model (1.3) has a unique positive equilibrium E¯=(x¯,y¯)=((f-bmn)/b(f+n),(bm+1)/b(f+n)).
(c)
B
1
<
0, (2.9) has a unique positive root x*=(A1+C1)/2, model (1.3) has a unique positive equilibrium E*=(x*,y*)=((A1+C1)/2,(2m+A1+C1)/2f).Next, we discuss the local stabilities of these equilibria. Easy to obtain the following results:Lemma 2.4.
(
i
)  E1=(1/b,0) is a saddle point;(ii)
ifc>n, E2=(0,m/f) is a stable node point. If c≤n, E2=(0,m/f) is a saddle point.Following, letE=(x,y) be an arbitrary positive equilibrium. The Jacobian matrix for E=(x,y) is given by
(2.12)J(E)=(ax(-b+cy(x+ny)2)ax(-cx+ny+cny(x+ny)2)ef-e).Then we can get the following:(2.13)tr(J(E))=ax(-b+cy(x+ny)2)-e,det(J(E))=aex[b(f+n)2x2+2bmn(f+n)x-cfm+bm2n2](fx+mn+nx)2.So the sign ofdet(J(E)) is determined by
(2.14)F(x)≜b(f+n)2x2+2bmn(f+n)x-cfm+bm2n2.Theorem 2.5.
For model (1.3), the stabilities of E3 and E4 are as follows.(a)
The positive equilibriumE3=(xe1,(m+xe1)/f)=((A1+C1)/2,(2m+A1+C1)/2f) is
(a1)
stable if and only if(2.15)e>axe1(-b+cf(m+xe1)(fxe1+n(m+xe1))2);(a2)
unstable if and only if(2.16)e<axe1(-b+cf(m+xe1)(fxe1+n(m+xe1))2).(b)
The positive equilibriumE4=(xe2,(m+xe2)/f)=((A1-C1)/2,(2m+A1-C1)/2f) is a saddle point.Proof.
(a) At the pointE3, we have the following:
(2.17)F(xe1)=b(f+n)2(A1+C1)24+bmn(f+n)(A1+C1)+bm2n2-cfm=b(f+n)2C122+C1(bmn(f+n)+b(f+n)2A12)>0,thus,det(J(E3))>0, and the nature of singularity E3 depends on the trace given by
(2.18)tr(J(E3))=axe1(-b+cf(m+xe1)(fxe1+n(m+xe1))2)-e.Clearly,
(a1)
Ife>axe1(-b+cf(m+xe1)/(fxe1+n(m+xe1))2), then tr(J(E3))<0, and E3=(xe1,(m+xe1)/f) is stable.(a2)
Ife<axe1(-b+cf(m+xe1)/(fxe1+n(m+xe1))2), then tr(J(E3))>0, and E3=(xe1,(m+xe1)/f) is unstable.(b)
At the pointE4, we have the following:
(2.19)F(xe2)=b(f+n)2(A1-C1)24+bmn(f+n)(A1-C1)+bm2n2-cfm=b(f+n)22C1(A1-C1)-bmn(f+n)C1-2fm(f+n)B1<0,thendet(J(E4))<0, and E4 is a saddle point.Figure1 shows the dynamics of model (1.3). In this case, E0=(0,0) is a nodal source, E1=(4,0) is a saddle point, E2=(0,071429) is a nodal sink, E4=(0.31772,1.5846) is a saddle point, and E3=(0.96454,2.3932) is locally asymptotically stable. There exists a separatrix curve determined by the stable manifold of E4, which divides the behavior of trajectories, that is, the stable manifold of saddle E4 split the feasible region into two parts such that orbits initiating inside tend to the positive equilibrium E3, while orbits initiating outside tend to E3 except for the stable manifolds of E4.Figure 1
The phase portraits of model (1.3). The parameters are taken as a=2,b=0.25,c=0.875,e=1,m=0.95,andn=0.75,f=0.8. In this case, E0=(0,0) is a source, E1=(4,0) and E4=(0.31772,1.5846) are saddle points, and E2=(0,071429) is a nodal sink, which is stable. The positive equilibrium E3=(0.96454,2.3932) is locally asymptotically stable. The dashed curve is the x-nullcline f(x,y)=0, and the dotted curve is the y-nullcline g(x,y)=0.Theorem 2.6.
The singularityE¯=(x¯,(m+x¯)/f)=((f-bmn)/b(f+n),(bm+1)/b(f+n)) is(a)
stable if and only if(2.20)e>ax¯(-b+cf(m+x¯)(fx¯+n(m+x¯))2);(b)
unstable if and only if(2.21)e<ax¯(-b+cf(m+x¯)(fx¯+n(m+x¯))2).Proof.
At the pointE¯=(x¯,(m+x¯)/f), we have
(2.22)F(x¯)=b(f+n)2A12+2bmn(f+n)A1+bm2n2-cfm=b(f+n)2A12+bmn(n+f)A1>0,
as F(x¯)>0, then det(J(E¯))>0, and the nature of singularity E¯=(x¯,(m+x¯)/f) depends on the trace given by
(2.23)tr(J(E¯))=ax¯(-b+cf(m+x¯)(fx¯+n(m+x¯))2)-e.
Clearly,(a)
Ife>ax¯(-b+cf(m+x¯)/(fx¯+n(m+x¯))2), then tr(J(E¯))<0, and E¯=(x¯,(m+x¯)/f) is stable.(b)
Ife<ax¯(-b+cf(m+x¯)/(fx¯+n(m+x¯))2), then tr(J(E¯))>0, and E¯=(x¯,(m+x¯)/f) is unstable.Theorem 2.7.
SingularityE*=(x*,(m+x*)/f) is(a)
stable if and only if(2.24)e>ax*(-b+cf(m+x*)(fx*+n(m+x*))2);(b)
unstable if and only if(2.25)e<ax*(-b+cf(m+x*)(fx*+n(m+x*))2).Proof.
At the pointE*=(x*,(m+x*)/f), as B1<0,C12=A12-4B, so
(2.26)F(x*)=b(f+n)2(A1+C1)24+bmn(f+n)(A1+C1)+bm2n2-cfm=b(f+n)2(A1+C1)24+bmn(f+n)C1+m(n+f)(n-c)>0,
as F(x*)>0, then det(J(E*))>0, and the nature of singularity E*=(x*,(m+x*)/f) depends on the trace given by
(2.27)tr(J(E*))=ax*(-b+cf(m+x*)(fx*+n(m+x*))2)-e.
Clearly,(a)
ife>ax*(-b+cf(m+x*)/(fx*+n(m+x*))2), then tr(J(E*))<0, E*=(x*,(m+x*)/f) is stable;(b)
ife<ax*(-b+cf(m+x*)/(fx*+n(m+x*))2), then tr(J(E*))>0, E*=(x*,(m+x*)/f) is unstable.Theorem 2.8.
SingularityEe=(xe,(m+xe)/f) is(a)
a non-hyperbolic attractor node if and only if(2.28)e>axe(-b+cf(m+xe)(fxe+n(m+xe))2);(b)
a non-hyperbolic repellor node if and only if(2.29)e<axe(-b+cf(m+xe)(fxe+n(m+xe))2).Proof.
At the point(xe,(m+xe)/f), we have
(2.30)F(xe)=b(f+n)2A124+bmn(f+n)A1+bm2n2-cfm=0,
then, det(J(Ee))=0, and tr(J(Ee))=axe(-b+cf(m+xe)/(fxe+n(m+xe))2)-e.
Hence we can conclude (2.28) and (2.29).Theorem 2.9.
Ifn<c<n(bm+1) and fn+(n-c)(bmn+n-c)<0 hold the boundary equilibria E2=(0,m/f) of model (1.3) is globally asymptotically stable.Proof.
Sincelimsupt→∞x(t)≤1/b and liminft→∞y(t)≥(bmn+n-c)/bfn, from the first equation of model (1.3), for any μ>0, there exists a T1>0, for all t≥T1, we have
(2.31)dxdt=ax(1-bx-cx/y+n)<ax(1-bx-c(bmn+n-c-bfμ)fn(1+bμ)+n(bmn+n-c-bfμ)).
From the arbitrariness of μ>0, we can get that
(2.32)dxdt<ax(1-bx-c(bmn+n-c)n(bmn+f+n-c)).
As1-c(bmn+n-c)/n(bmn+f+n-c)<0, by standard comparison arguments, it follows that
(2.33)limsupt→∞x(t)≤0,
thus,
(2.34)limt→∞x(t)=0.
As a result, using the second equation of model (1.3), one can easily know that limt→∞y(t)=m/f. The proof is complete.Figure2 shows the dynamics of model (1.3). In this case, E0=(0,0) is a nodal source, E1=(4,0) is a saddle point, and E2=(0,1.1875) is globally asymptotically stable, that is, all orbits tend to the equilibrium E2 for any initial values.Figure 2
The phase portraits of model (1.3). The parameters are taken as a=2,b=0.25,c=0.875,e=1,m=0.95,andn=0.7,f=0.8. In this case, E0=(0,0) is a nodal source, E1=(4,0) is a saddle point. The boundary equilibria E2=(0,1.1875) is globally asymptotically stable. The dashed curve is the x-nullcline f(x,y)=0, and the dotted curve is the y-nullcline g(x,y)=0.Theorem 2.10.
AssumeB1<0, E*=(x*,y*) is globally asymptotically stable, if the following conditions hold(i)
1-bx*<b(ζ+nδ);(ii)
(e(η+nς)-(x*/y*)(1-bx*)(ζ+m))2<4aef(b(ζ+nδ)-(1-bx*))(η+m)(η+nς), where η=1/b,ς=(1+bm)/bf,ζ=(n-c)/bn,δ=(bmn+n-c)/bnf.Proof.
Define a Lyapunov function:(2.35)V(x,y)=∫x*x‍ξ-x*ξdξ+∫y*y‍η-y*ηdη,
so,
(2.36)dVdt=x-x*xdxdt+(y-y*)ydydt=a(x-x*)(1-bx-cyx+ny)+e(y-y*)(1-fyx+m)=a(-b+cy*(x*+ny*)(x+ny))(x-x*)2+(-acx*(x*+ny*)(x+ny)+efy*(x*+m)(x+m))(x-x*)(y-y*)-ef(y-y*)2x+m=-A2(x-x*)2+B2(x-x*)(y-y*)-C2(y-y*)2,
let
(2.37)A2=a(b-cy*(x*+ny*)(x+ny))=a(b-1-bx*x+ny),B2=-acx*(x*+ny*)(x+ny)+efy*(x*+m)(x+m)=-ax*(1-bx*)y*(x+ny)+ex+m,C2=efx+m.
If
(2.38)A2>0,thatis,1-bx*<b(x+ny)<b(ζ+nδ),B22<4A2C2,thatis,(-ax*(1-bx*)y*(x+ny)+ex+m)2<4aefx+m(b-1-bx*x+ny)hold,thedV/dt<0.
Hence,
(2.39)(e(x+ny)-ax*y*(1-bx*)(x+m))2<4aef(b(x+ny)-(1-bx*))(x+m)(x+ny),
according to Lemma 2.2, we can obtain that
(2.40)(e(x+ny)-ax*y*(1-bx*)(x+m))2<(e(η+nς)-ax*y*(1-bx*)(ζ+m))2,4aef(b(ζ+nδ)-(1-bx*))(η+m)(η+nς)<4aef(b(x+ny)-(1-bx*))(x+m)(x+ny),
so,
(2.41)(e(η+nς)-x*y*(1-bx*)(ζ+m))2<4aef(b(ζ+nδ)-(1-bx*))(η+m)(η+nς).
ConsideringA2>0,B22<4A2C2, we obtain dV/dt<0. This ends the proof.Figure3 shows the dynamics of model (1.3) for the case of B1<0. In this case, E0=(0,0) is a nodal source, E1=(4,0) is a saddle point, E2=(0,1.1875) is a saddle point, and E*=(1.4761,3.0326) is globally asymptotically stable, that is, all orbits tend to the equilibrium E* for any initial values.Figure 3
The phase portraits of model (1.3). The parameters are taken as a=2,b=0.25,c=0.875,e=1,  m=0.95,n=0.9,andf=0.8. In this case, E1=(4,0) and E2=(0,1.1875) are saddle points, the positive interior equilibrium point E*=(1.4761,3.0326) is globally asymptotically stable. The dashed curve is the x-nullcline f(x,y)=0, and the dotted curve is the y-nullcline g(x,y)=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)a→a+αB˙1(t),e→e+βB˙2(t),
then the stochastic version of model (1.3) is given by the following Itô type
(3.2)dx=x(1-bx-cyx+ny)(adt+αdB1(t)),dy=y(1-fyx+m)(edt+βdB2(t)),
where Bi(t),i=(1,2) are the 1-dimensional standard Brownian motion defined on a complete probability space (Ω,ℱ,𝒫) with a filtration {ℱt}t∈R+ satisfying the usual conditions (i.e., it is right continuous and increasing while F0 contains all 𝒫-null sets) and B˙1(t), B˙2(t) are, respectively, white noises with possible intensity α2, β2.


### 3.1. Existence of Global Positive Solutions

Lemma 3.1.
There is a unique local positive solution(x(t),y(t)) for t∈[0,τe) to model (3.2) almost surely (a.s.) for the initial value x0>0, y0>0, where τe is the explosion time.Proof.
Consider the equations(3.3)du={a[1-beu-cevnev+eu]-12α2[1-beu-cevnev+eu]2}dt+α[1-beu-cevnev+eu]dB1(t),dv={e[1-fevm+eu]-12β2[1-fevm+eu]2}dt+β[1-fevm+eu]dB2(t),
on t≥0 with initial value u(0)=lnx0, v(0)=lny0.
It is easy to see that the coefficients of model (3.3) satisfy the local Lipschitz condition, then there is a unique local solution u(t),v(t) on [0,τe) [13]. Therefore, by Itô formula, x(t)=eu(t),y(t)=ev(t) are the unique positive local solutions to model (3.3) with initial value x0>0,y0>0.Lemma3.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, τe=∞.Theorem 3.2.
Consider model (3.2), for any given initial value (x0,y0)∈R+2, there is a unique solution (x(t),y(t)) on t≥0 and the solution will remain in R+2 with probability 1, where R+2={x∈R2∣xi>0,i=1,2}.Proof.
For convenience of statement, we introduce some notations. Define(3.4)F(x,y)=cyny+x,G(x,y)=fym+x.
Let n0>0 be sufficiently large for x0 and y0 lying within the interval [1/n0,n0]. For each integer n>n0, define the stopping times
(3.5)τn=inf{t∈[0,τe]:x(t)∉(1n,n)ory(t)∉(1n,n)}.
Assumeinf∅=∞ (as usual ∅= the empty set). Clearly, τn is increasing as n→∞. Let τ∞=limn→∞τn, then τ∞≤τe a.s. If we can show that τ∞=∞ a.s., then τe=∞ a.s., and N(t)∈R+ a.s., for all t≥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 T>0 and ɛ∈(0,1) such that 𝒫{τ∞≤T}>ɛ. There is an integer n1≥n0, such that
(3.6)𝒫{τn≤T}≥ɛ,n≥n1.
Define a C2 function V:R+2→R+ by
(3.7)V(x,y)=(x-1-lnx2)+(y-1-lny2).
Obviously, V(x,y) is nonnegative. And
(3.8)dV=a(x-1)2(1-bx-F(x,y))dt+α2(2-x)8[1-bx-F(x,y)]2dt+α(x-1)2[1-bx-F(x,y)]dB1t+e(y-1)2[1-G(x,y)]dt+β2(2-y)8[1-G(x,y)]2dt+β(y-1)2[1-G(x,y)]dB2t≤a2[x+bx+F(x,y)]dt+e2[y+G(x,y)]dt+α28[2bx3/2-b2x5/2+2xF(x,y)+2+2b2x2+4bxF(x,y)+2F2(x,y)]dt+β28[2yG(x,y)+2+2G(x,y)2]dt+α2(x-1)[1-bxa-F(x,y)]dB1(t)+β2(y-1)[1-G(x,y)]dB2(t)≤k1+k2+α2(x-1)[1-bx-F(x,y)]dB1(t)+β2(y-1)[1-G(x,y)]dB2(t),
where k1,k2 are positive numbers. Integrating both sides of the above inequality from 0 to τn∧T and then taking expectations, yields
(3.9)EV(x(τn∧T),y(τn∧T))≤V(x0,y0)+(K1+K2)T.
Set Ωn={τn≤T} for n≥n1 and by (3.6), we have 𝒫(Ωn)≥ɛ. Note that for every ω∈Ωn, there is some i such that xi(τn,ω) equals either n or 1/n for i=1,2, hence V(x(τn,ω),y(τn,ω)) is no less than min{(n-1-lnn/2),(1/n-1-ln(1/n)/2)}. It then follows from (3.9) that
(3.10)V(x0,y0)+(K1+K2)T≥E[IΩn(ω)V(x(τn),y(τn))]≥εmin{(n-1-lnn2),(1n-1-lnn2)},
where IΩn is the indicator function of Ωn. Let n→∞, then
(3.11)∞>V(x0,y0)+(K1+K2)T=∞a.s.
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 (x0,y0)∈R+2, the solution of model (3.2) has the property that
(3.12)limsupt→∞P{|x(t),y(t)|>χ}≤ɛ.Lemma 3.4.
For anyθ∈(0,1), there is a positive constant H=H(θ)>0, which is independent of the initial data (x0,y0)∈R+2, such that the solution of model (3.2) has the property that
(3.13)limsupt→∞(E|(x(t),y(t))|θ)≤H.Proof.
Define(3.14)V(x,y)=xθ+yθ,(x,y)∈R+2.
For the sake of discussion, we rewrite the above as
(3.15)V(x,y)=V1+V2,V1=xθ,V2=yθ.
By using the Itô formula, we have(3.16)dV(x,y)=LV(x,y)dt+θαxθ[1-bx-F(x,y)]dB1(t)+θβyθ[1-G(x,y)]dB2(t),LV(x,y)=aθxθ[1-bx-F(x,y)]+θ(θ-1)α2xθ2[1-bx-F(x,y)]2+eθyθ[1-G(x,y)]+θ(θ-1)β2yθ2[1-G(x,y)]2,
as 0<θ<1, we have
(3.17)LV(x,y)≤aθxθ-θ(1-θ)α2xθ2[1-bx-F(x,y)]2+eθyθ[1-G(x,y)]-θ(1-θ)β2yθ2[1-G(x,y)]2≤aθxθ-θ(1-θ)α2xθ2[b2x2-2bx-2cn]+eθyθ-θ(1-θ)β2yθ2[-2fym]=aθxθ-θ(1-θ)α2xθ2[b2x2-2bx-2cn]+xθ+eθyθ-θ(1-θ)β2yθ2[-2fym]+yθ-V(x,y)≤H1-V(x,y),
here H1 is an integer. So
(3.18)dV(x,y)≤(H1-V(x,y))dt+θαxθ[1-bx-F(x,y)]dB1(t).
Now, by using the Itô formula again, we have(3.19)d(etV(x,y))=et[V(x,y)dt+dV(x,y)]≤etH1dt+etθαxθ[1-bx-F(x,y)]dB1(t)+etθβyθ[1-G(x,y)]dB2(t).
According to the above, we can easily get(3.20)etE(V(x,y))≤V(x0,y0)+H1(et-1),limsupt→∞EV(x(t),y(t))≤H.
On the other hand,(3.21)|(x,y)|θ≤2θ/2max{xθ,yθ}≤2θ/2V(x,y),
so, we can obtain
(3.22)limsupt→∞E|(x(t),y(t))|θ≤2θ/2,limsupt→∞EV(x(t),y(t))≤2θ/2H.
SetK(θ)=2θ/2H, this completes the proof.Theorem 3.5.
The solutions of model (3.2) is stochastically ultimately bounded.Proof.
By Lemma3.4, let θ=1/2, there is a K1>0 such that
(3.23)limsupt→∞E|(x(t),y(t))|1/2≤K1.
Now, for anyɛ>0, let χ=K12/ɛ2. Then by Chebyshevs inequality
(3.24)P{|x(t),y(t)|>χ}≤E|(x(t),y(t))|1/2χ.
Hence,(3.25)limsupt→∞P{|x(t),y(t)|>χ}≤ɛ.
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 K=K(θ)>0, let (x(t),y(t)) be the solution of model (3.2) with any initial value (x0,y0), then
(3.26)limsupt→∞1t∫0t‍E[x2+θ(s)+y2+θ(s)]ds≤K.Proof.
SetV(x,y):R+2→R+, from Lemma 3.4, we have
(3.27)LV(x,y)≤aθxθ-θ(1-θ)α2xθ2[b2x2-2bx-2cn]+eθyθ-θ(1-θ)β2yθ2(-2fym),
so, we get
(3.28)LV+θ(1-θ)α2b2xθ+24+θ(1-θ)β2fyθ+24m≤aθxθ-θ(1-θ)α2b2xθ+24+θ(1-θ)α2xθ(bx+cn)+eθyθ+θ(1-θ)β2fyθ+1m+θ(1-θ)β2fyθ+24m≤M,
here M is a positive number, so
(3.29)LV+θ(1-θ)γ2σ2[xθ+2+yθ+2]4≤θ(1-θ)α2b2xθ+24+θ(1-θ)β2fyθ+24m≤M.
Letγ2σ2=min{α2b2,β2f/m}, we can get
(3.30)EV(x,y)+θ(1-θ)γ2σ24∫0t‍E[x2+θ(s)+y2+θ(s)]ds≤V(x0+y0)+Mt.
That is,(3.31)∫0t‍E[x2+θ(s)+y2+θ(s)]ds≤4[V(x0+y0)+Mt]θ(1-θ)γ2σ2.
SetK=4M/θ(1-θ)γ2σ2, so
(3.32)limsupt→∞1t∫0t‍E[x2+θ(s)+y2+θ(s)]ds≤K.


### 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 E*=(x*,y*) of model (3.2).Theorem 3.7.
Let(3.33)A3=-ab+acy*(ny*+x*)+α2b2x*2+α2c2x*(y*)22(ny*+x*)2+β2f2(y*)32(m+x*)2,B3=α2c2(x*)32(ny*+x*)2+β2(fm+x*)2y*2(m+x*)2,C3=α2bc(x*)2(ny*+x*)+α2c2y*(x*)2(ny*+x*)2+fy*(m+x*)+β2(fm+x*)f(y*)2(m+x*)2,
if
(3.34)(bm+1)n2+[(1-b)f-c]n-cf≥0,[(m-1)b+1]n-c≥0,A3<0,C32-4A3B3<0,
when B1<0, then the equilibrium position E*=(x*,y*) of model (3.2) is stochastically asymptotically stable in the large, that is, for any initial data (x(0),y(0)), the solution of model (3.2) has the property that
(3.35)limt→∞(x(t))=x*,limt→∞(y(t))=y*,a.s.Proof.
From the theory of stability of stochastic differential equations [13], we only need to establish a Lyapunov function V(z) satisfying LV(z)≤0 and the identity holds if and only if z=z*, where z=z(t) is the solution of the stochastic differential equation
(3.36)dz=f(z(t),t)dt+g(z(t),t)dB(t),z* is the equilibrium position of model (3.36), and
(3.37)LV(z,t)=Vz(z,t)f(z,t)+12tr[gT(z,t)Vzz(z,t)g(z,t)].
Define Lyapunov functions
(3.38)V1(x)=∫x*x‍ξ-x*ξdξ,V2(x)=∫y*y‍η-y*ηdη,
the nonnegativity of this function can be observed from u-1-lnu≥0 on u>0. Applying Itô formula leads to
(3.39)LV1(x)=a(x-x*)[1-bx-F(x,y)]+α2x*2[1-bx-F(x,y)]2=a(x-x*)[bx*-bx+F(x*,y*)-F(x,y)]+α2x*2[bx*-bx+F(x*,y*)-F(x,y)]2={-ab+acy*(ny*+x*)(ny+x)+α2b2x*2-bcα2x*y*(ny*+x*)(ny+x)+c2α2x*(y*)22(ny*+x*)2(ny+x)2}(x-x*)2+c2α2(x*)32(ny*+x*)2(ny+x)2(y-y*)2-acx*(ny*+x*)(ny+x)(x-x*)(y-y*)+bcα2(x*)2(ny*+x*)(ny+x)(x-x*)(y-y*)-c2α2y*(x*)2(ny*+x*)2(ny+x)2(x-x*)(y-y*),
when ny+x≥1, according to Lemma 2.2, we have (bm+1)n2+[(1-b)f-c]n-cf≥0, so,
(3.40)LV1(x)≤{-ab+acy*(ny*+x*)+α2b2x*2+c2α2x*(y*)22(ny*+x*)2}(x-x*)2+c2α2(x*)32(ny*+x*)2(y-y*)2+{bcα2(x*)2(ny*+x*)+c2α2y*(x*)2(ny*+x*)2}|(x-x*)||(y-y*)|.
Similarly, by using the Itô formula, we can obtain(3.41)LV2(y)=e(y-y*)[1-G(x,y)]+β2y*2[1-G(x,y)]2,
when m+x≥1, according to Lemma 2.2, we have [(m-1)b+1]n-c≥0, so
(3.42)LV2(y)≤β2(fm+x*)2y*2(m+x*)2(y-y*)2+β2f2(y*)32(m+x*)2(x-x*)2+{fy*(m+x*)+β2(fm+x*)f(y*)2(m+x*)2}|(x-x*)||(y-y*)|,
hence,
(3.43)LV(x,y)=LV1(x)+LV2(y)≤{-ab+acy*(ny*+x*)+α2b2x*2+c2α2x*(y*)22(ny*+x*)2+β2f2(y*)32(m+x*)2}(x-x*)2+{c2α2(x*)32(ny*+x*)2+β2(fm+x*)2y*2(m+x*)2}(y-y*)2+{bcα2(x*)2a2(ny*+x*)+c2α2y*(x*)2(ny*+x*)2+fy*(m+x*)+β2(fm+x*)f(y*)2(m+x*)2}|(x-x*)||(y-y*)|.
So,(3.44)LV(x,y)≤A3(x-x*)2+B3(y-y*)2+C3|(x-x*)||(y-y*)|.
Let(3.45)|z-z*|=(|(x-x*)|,|(y-y*)|)T,
then we get
(3.46)LV(x,y)≤12(|z-z*|)T(2A3C3C32B3)|z-z*|=A3|(x-x*)|2+C3|(x-x*)||(y-y*)|+B3|(y-y*)|2.
Clearly, if (3.34) hold, then the above inequality implies LV(x,y)<0 along all trajectories in the first quadrant except (x*,y*). Then the desired assertion (3.35) follows immediately. This completes the proof.


## 3.1. Existence of Global Positive Solutions

Lemma 3.1.
There is a unique local positive solution(x(t),y(t)) for t∈[0,τe) to model (3.2) almost surely (a.s.) for the initial value x0>0, y0>0, where τe is the explosion time.Proof.
Consider the equations(3.3)du={a[1-beu-cevnev+eu]-12α2[1-beu-cevnev+eu]2}dt+α[1-beu-cevnev+eu]dB1(t),dv={e[1-fevm+eu]-12β2[1-fevm+eu]2}dt+β[1-fevm+eu]dB2(t),
on t≥0 with initial value u(0)=lnx0, v(0)=lny0.
It is easy to see that the coefficients of model (3.3) satisfy the local Lipschitz condition, then there is a unique local solution u(t),v(t) on [0,τe) [13]. Therefore, by Itô formula, x(t)=eu(t),y(t)=ev(t) are the unique positive local solutions to model (3.3) with initial value x0>0,y0>0.Lemma3.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, τe=∞.Theorem 3.2.
Consider model (3.2), for any given initial value (x0,y0)∈R+2, there is a unique solution (x(t),y(t)) on t≥0 and the solution will remain in R+2 with probability 1, where R+2={x∈R2∣xi>0,i=1,2}.Proof.
For convenience of statement, we introduce some notations. Define(3.4)F(x,y)=cyny+x,G(x,y)=fym+x.
Let n0>0 be sufficiently large for x0 and y0 lying within the interval [1/n0,n0]. For each integer n>n0, define the stopping times
(3.5)τn=inf{t∈[0,τe]:x(t)∉(1n,n)ory(t)∉(1n,n)}.
Assumeinf∅=∞ (as usual ∅= the empty set). Clearly, τn is increasing as n→∞. Let τ∞=limn→∞τn, then τ∞≤τe a.s. If we can show that τ∞=∞ a.s., then τe=∞ a.s., and N(t)∈R+ a.s., for all t≥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 T>0 and ɛ∈(0,1) such that 𝒫{τ∞≤T}>ɛ. There is an integer n1≥n0, such that
(3.6)𝒫{τn≤T}≥ɛ,n≥n1.
Define a C2 function V:R+2→R+ by
(3.7)V(x,y)=(x-1-lnx2)+(y-1-lny2).
Obviously, V(x,y) is nonnegative. And
(3.8)dV=a(x-1)2(1-bx-F(x,y))dt+α2(2-x)8[1-bx-F(x,y)]2dt+α(x-1)2[1-bx-F(x,y)]dB1t+e(y-1)2[1-G(x,y)]dt+β2(2-y)8[1-G(x,y)]2dt+β(y-1)2[1-G(x,y)]dB2t≤a2[x+bx+F(x,y)]dt+e2[y+G(x,y)]dt+α28[2bx3/2-b2x5/2+2xF(x,y)+2+2b2x2+4bxF(x,y)+2F2(x,y)]dt+β28[2yG(x,y)+2+2G(x,y)2]dt+α2(x-1)[1-bxa-F(x,y)]dB1(t)+β2(y-1)[1-G(x,y)]dB2(t)≤k1+k2+α2(x-1)[1-bx-F(x,y)]dB1(t)+β2(y-1)[1-G(x,y)]dB2(t),
where k1,k2 are positive numbers. Integrating both sides of the above inequality from 0 to τn∧T and then taking expectations, yields
(3.9)EV(x(τn∧T),y(τn∧T))≤V(x0,y0)+(K1+K2)T.
Set Ωn={τn≤T} for n≥n1 and by (3.6), we have 𝒫(Ωn)≥ɛ. Note that for every ω∈Ωn, there is some i such that xi(τn,ω) equals either n or 1/n for i=1,2, hence V(x(τn,ω),y(τn,ω)) is no less than min{(n-1-lnn/2),(1/n-1-ln(1/n)/2)}. It then follows from (3.9) that
(3.10)V(x0,y0)+(K1+K2)T≥E[IΩn(ω)V(x(τn),y(τn))]≥εmin{(n-1-lnn2),(1n-1-lnn2)},
where IΩn is the indicator function of Ωn. Let n→∞, then
(3.11)∞>V(x0,y0)+(K1+K2)T=∞a.s.
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 (x0,y0)∈R+2, the solution of model (3.2) has the property that
(3.12)limsupt→∞P{|x(t),y(t)|>χ}≤ɛ.Lemma 3.4.
For anyθ∈(0,1), there is a positive constant H=H(θ)>0, which is independent of the initial data (x0,y0)∈R+2, such that the solution of model (3.2) has the property that
(3.13)limsupt→∞(E|(x(t),y(t))|θ)≤H.Proof.
Define(3.14)V(x,y)=xθ+yθ,(x,y)∈R+2.
For the sake of discussion, we rewrite the above as
(3.15)V(x,y)=V1+V2,V1=xθ,V2=yθ.
By using the Itô formula, we have(3.16)dV(x,y)=LV(x,y)dt+θαxθ[1-bx-F(x,y)]dB1(t)+θβyθ[1-G(x,y)]dB2(t),LV(x,y)=aθxθ[1-bx-F(x,y)]+θ(θ-1)α2xθ2[1-bx-F(x,y)]2+eθyθ[1-G(x,y)]+θ(θ-1)β2yθ2[1-G(x,y)]2,
as 0<θ<1, we have
(3.17)LV(x,y)≤aθxθ-θ(1-θ)α2xθ2[1-bx-F(x,y)]2+eθyθ[1-G(x,y)]-θ(1-θ)β2yθ2[1-G(x,y)]2≤aθxθ-θ(1-θ)α2xθ2[b2x2-2bx-2cn]+eθyθ-θ(1-θ)β2yθ2[-2fym]=aθxθ-θ(1-θ)α2xθ2[b2x2-2bx-2cn]+xθ+eθyθ-θ(1-θ)β2yθ2[-2fym]+yθ-V(x,y)≤H1-V(x,y),
here H1 is an integer. So
(3.18)dV(x,y)≤(H1-V(x,y))dt+θαxθ[1-bx-F(x,y)]dB1(t).
Now, by using the Itô formula again, we have(3.19)d(etV(x,y))=et[V(x,y)dt+dV(x,y)]≤etH1dt+etθαxθ[1-bx-F(x,y)]dB1(t)+etθβyθ[1-G(x,y)]dB2(t).
According to the above, we can easily get(3.20)etE(V(x,y))≤V(x0,y0)+H1(et-1),limsupt→∞EV(x(t),y(t))≤H.
On the other hand,(3.21)|(x,y)|θ≤2θ/2max{xθ,yθ}≤2θ/2V(x,y),
so, we can obtain
(3.22)limsupt→∞E|(x(t),y(t))|θ≤2θ/2,limsupt→∞EV(x(t),y(t))≤2θ/2H.
SetK(θ)=2θ/2H, this completes the proof.Theorem 3.5.
The solutions of model (3.2) is stochastically ultimately bounded.Proof.
By Lemma3.4, let θ=1/2, there is a K1>0 such that
(3.23)limsupt→∞E|(x(t),y(t))|1/2≤K1.
Now, for anyɛ>0, let χ=K12/ɛ2. Then by Chebyshevs inequality
(3.24)P{|x(t),y(t)|>χ}≤E|(x(t),y(t))|1/2χ.
Hence,(3.25)limsupt→∞P{|x(t),y(t)|>χ}≤ɛ.
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 K=K(θ)>0, let (x(t),y(t)) be the solution of model (3.2) with any initial value (x0,y0), then
(3.26)limsupt→∞1t∫0t‍E[x2+θ(s)+y2+θ(s)]ds≤K.Proof.
SetV(x,y):R+2→R+, from Lemma 3.4, we have
(3.27)LV(x,y)≤aθxθ-θ(1-θ)α2xθ2[b2x2-2bx-2cn]+eθyθ-θ(1-θ)β2yθ2(-2fym),
so, we get
(3.28)LV+θ(1-θ)α2b2xθ+24+θ(1-θ)β2fyθ+24m≤aθxθ-θ(1-θ)α2b2xθ+24+θ(1-θ)α2xθ(bx+cn)+eθyθ+θ(1-θ)β2fyθ+1m+θ(1-θ)β2fyθ+24m≤M,
here M is a positive number, so
(3.29)LV+θ(1-θ)γ2σ2[xθ+2+yθ+2]4≤θ(1-θ)α2b2xθ+24+θ(1-θ)β2fyθ+24m≤M.
Letγ2σ2=min{α2b2,β2f/m}, we can get
(3.30)EV(x,y)+θ(1-θ)γ2σ24∫0t‍E[x2+θ(s)+y2+θ(s)]ds≤V(x0+y0)+Mt.
That is,(3.31)∫0t‍E[x2+θ(s)+y2+θ(s)]ds≤4[V(x0+y0)+Mt]θ(1-θ)γ2σ2.
SetK=4M/θ(1-θ)γ2σ2, so
(3.32)limsupt→∞1t∫0t‍E[x2+θ(s)+y2+θ(s)]ds≤K.


## 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 E*=(x*,y*) of model (3.2).Theorem 3.7.
Let(3.33)A3=-ab+acy*(ny*+x*)+α2b2x*2+α2c2x*(y*)22(ny*+x*)2+β2f2(y*)32(m+x*)2,B3=α2c2(x*)32(ny*+x*)2+β2(fm+x*)2y*2(m+x*)2,C3=α2bc(x*)2(ny*+x*)+α2c2y*(x*)2(ny*+x*)2+fy*(m+x*)+β2(fm+x*)f(y*)2(m+x*)2,
if
(3.34)(bm+1)n2+[(1-b)f-c]n-cf≥0,[(m-1)b+1]n-c≥0,A3<0,C32-4A3B3<0,
when B1<0, then the equilibrium position E*=(x*,y*) of model (3.2) is stochastically asymptotically stable in the large, that is, for any initial data (x(0),y(0)), the solution of model (3.2) has the property that
(3.35)limt→∞(x(t))=x*,limt→∞(y(t))=y*,a.s.Proof.
From the theory of stability of stochastic differential equations [13], we only need to establish a Lyapunov function V(z) satisfying LV(z)≤0 and the identity holds if and only if z=z*, where z=z(t) is the solution of the stochastic differential equation
(3.36)dz=f(z(t),t)dt+g(z(t),t)dB(t),z* is the equilibrium position of model (3.36), and
(3.37)LV(z,t)=Vz(z,t)f(z,t)+12tr[gT(z,t)Vzz(z,t)g(z,t)].
Define Lyapunov functions
(3.38)V1(x)=∫x*x‍ξ-x*ξdξ,V2(x)=∫y*y‍η-y*ηdη,
the nonnegativity of this function can be observed from u-1-lnu≥0 on u>0. Applying Itô formula leads to
(3.39)LV1(x)=a(x-x*)[1-bx-F(x,y)]+α2x*2[1-bx-F(x,y)]2=a(x-x*)[bx*-bx+F(x*,y*)-F(x,y)]+α2x*2[bx*-bx+F(x*,y*)-F(x,y)]2={-ab+acy*(ny*+x*)(ny+x)+α2b2x*2-bcα2x*y*(ny*+x*)(ny+x)+c2α2x*(y*)22(ny*+x*)2(ny+x)2}(x-x*)2+c2α2(x*)32(ny*+x*)2(ny+x)2(y-y*)2-acx*(ny*+x*)(ny+x)(x-x*)(y-y*)+bcα2(x*)2(ny*+x*)(ny+x)(x-x*)(y-y*)-c2α2y*(x*)2(ny*+x*)2(ny+x)2(x-x*)(y-y*),
when ny+x≥1, according to Lemma 2.2, we have (bm+1)n2+[(1-b)f-c]n-cf≥0, so,
(3.40)LV1(x)≤{-ab+acy*(ny*+x*)+α2b2x*2+c2α2x*(y*)22(ny*+x*)2}(x-x*)2+c2α2(x*)32(ny*+x*)2(y-y*)2+{bcα2(x*)2(ny*+x*)+c2α2y*(x*)2(ny*+x*)2}|(x-x*)||(y-y*)|.
Similarly, by using the Itô formula, we can obtain(3.41)LV2(y)=e(y-y*)[1-G(x,y)]+β2y*2[1-G(x,y)]2,
when m+x≥1, according to Lemma 2.2, we have [(m-1)b+1]n-c≥0, so
(3.42)LV2(y)≤β2(fm+x*)2y*2(m+x*)2(y-y*)2+β2f2(y*)32(m+x*)2(x-x*)2+{fy*(m+x*)+β2(fm+x*)f(y*)2(m+x*)2}|(x-x*)||(y-y*)|,
hence,
(3.43)LV(x,y)=LV1(x)+LV2(y)≤{-ab+acy*(ny*+x*)+α2b2x*2+c2α2x*(y*)22(ny*+x*)2+β2f2(y*)32(m+x*)2}(x-x*)2+{c2α2(x*)32(ny*+x*)2+β2(fm+x*)2y*2(m+x*)2}(y-y*)2+{bcα2(x*)2a2(ny*+x*)+c2α2y*(x*)2(ny*+x*)2+fy*(m+x*)+β2(fm+x*)f(y*)2(m+x*)2}|(x-x*)||(y-y*)|.
So,(3.44)LV(x,y)≤A3(x-x*)2+B3(y-y*)2+C3|(x-x*)||(y-y*)|.
Let(3.45)|z-z*|=(|(x-x*)|,|(y-y*)|)T,
then we get
(3.46)LV(x,y)≤12(|z-z*|)T(2A3C3C32B3)|z-z*|=A3|(x-x*)|2+C3|(x-x*)||(y-y*)|+B3|(y-y*)|2.
Clearly, if (3.34) hold, then the above inequality implies LV(x,y)<0 along all trajectories in the first quadrant except (x*,y*). 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 Figure4(a), we choose α=β=0, that is, without noise, we observe that the positive equilibrium E*=(1.4761,3.0326) is globally stable. In Figure 4(b), with noise densities α=0.8,β=0.4, starting with a homogeneous state E*=(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.The solution of the stochastic model (3.2). The parameters are taken as a=2,b=0.25,c=0.875,e=1,m=0.95,n=0.9,andf=0.8. (a) α=β=0; (b) α=0.8andβ=0.4.
(a)
(b)Figure 5
The solution of the stochastic model (3.2). The parameters are taken as α=2,β=4, a=2,b=0.25,c=0.875,e=1,m=0.95,n=0.9,andf=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.



---

*Source: 401031-2012-09-18.xml*