# Dynamics of a Stochastic Cooperative Predator-Prey System with Beddington-DeAngelis Functional Response

Authors: Dongwei Huang; Yu Li; Yongfeng Guo
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

Stochastic cooperative predator-prey system with Beddington-DeAngelis functional response is studied. It presents an investigation of dynamic properties of the system. Our results show that there exists a unique positive solution to the system for any positive initial value, and the positive solution is stochastically bounded. Moreover, under some conditions, we analyze global asymptotic stability of the positive solutions. With small environmental noises, the stochastic system is getting more similar to the corresponding deterministic system. Neither of the species in the system will die out. Finally, simulations are carried out to conform to our result.


---

## Body

## 1. Introduction

As we all know, in mathematical biology, predator-prey system, competitive system, and cooperative system are the three rudimentary and important ecological systems. The dynamic relationship between species has long been and will continue to be a dominant theme in ecology due to its universal existence and importance. It is well-known that predator-prey systems are very important and extensive in the nature fields. One significant component of the predator-prey relationship is predator’s functional response, that is, the rate of prey consumption by an average predator. There are many significant functional responses in order to model various different situations. In fact, most of the functional responses are prey-dependent; however, some biologists have argued that in many cases, especially when predators have to search for food and therefore have to share or compete for food, the traditional predator-prey systems with prey-dependent functional response fail to model the interference among predators, the functional response should be predator-dependent. In [1, 2], Beddington and DeAngelis proposed the following predator-prey model with Beddington-DeAngelis functional response:
(1)dxdt=x(r-θx)-Exya+bx+cy,dydt=-dy+βxya+bx+cy.
Skalski and Gilliam [3] compared statistical evidence from 19 predator-prey systems, and then they claimed that three predator-dependent functional responses (Hassell-Varley, Beddington-DeAngelis, and Growley-Martin) can provide better description of predator feeding over a range of predator-prey abundances. And the Beddington-DeAngelis type functional response was even suitable in some cases.But most of this work is restricted to predator-prey systems, little has been done for cooperative systems [4, 5]. May [6] suggested the following set of equations:
(2)dxdt=r1x(1-xa1+b1y-c1x),dydt=r2y(1-ya2+b2x-c2y),
to describe a pair of mutualists.However, there is often the interaction among multiple species in nature, whose relationship is more complex than those in two species. Therefore, it is more realistic to consider the multiple-species predator-prey systems. In order to continue studying such models, in this paper, we consider a cooperative predator-prey system with Beddington-DeAngelis functional responses at first:(3)dx=x(a1-b1x-c2y1+α2x+β2y-c3z1+α3x+β3z)dt,dy=y(a2-b2y-h2yf2+g2z+d2x1+α2x+β2y)dt,dz=z(a3-b3z-h3zf3+g3y+d3x1+α3x+β3z)dt,
where species x is the prey of y and z, y, and z are cooperative species. All the parameters in system (3) are positive constants.In fact, population dynamics is inevitably affected by environmental white noise which is an important component in an ecosystem. But the model (3) is deterministic and does not incorporate the effect of environmental noise. May [6] also pointed out the fact that due to environmental fluctuation, the birth rates, carrying capacity, competition coefficients, and other parameters involved in system exhibit random fluctuation to a greater or a lesser extent. Therefore many scholars rewrote the deterministic models as stochastic ones subjected to stochastic noises, for studying the effect of environmental variability on the population dynamics [7–9].The parameters in the real ecosystems are often subject to lots of environmental noises, since they relate to climate, geographical distribution, geological features, human disaster, human intervention, and other environmental factors. Therefore, the logistics and energy flow, in which they are determined by groups, are fluctuating. The oscillation in population biomass is directly manifested as birth and death rates of random perturbation. Currently, one of the main ways considered in the literature to model the effect of the environmental fluctuations in population dynamics is to assume that the most sensitive parameter is the intrinsic growth rate. Thus, in this paper we introduce some stochastic perturbation into the intrinsic growth rate. Therefore, the intrinsic growth rate can be written as an average growth rate adding some small random perturbed terms. In general, by the well-known central limit theorem, the small terms follow some normal distributions, so we can use standard Brownian motions to represent the environmental fluctuations.In this paper, taking into account the effect of randomly fluctuating environment, we introduce stochastic perturbation into growth ratesa1, a2, and a3 to become a1+σ1B˙1(t), a2+σ2B˙2(t), and a3+σ3B˙3(t) in system, where σi2 represents the intensity of the noise and B˙i(t) is a standard white noise, namely, is Bi(t) a standard Brownian motion defined on a complete probability space (Ω,F,P). Then the stochastic system takes the following form:
(4)dx=x(a1-b1x-c2y1+α2x+β2y-c3z1+α3x+β3z)dt+σ1xdB1(t),dy=y(a2-b2y-h2yf2+g2z+d2x1+α2x+β2y)dt+σ2ydB2(t),dz=z(a3-b3z-h3zf3+g3y+d3x1+α3x+β3z)dt+σ3zdB3(t).
Considering system (4), the initial conditions x(0)>0, y(0)>0, and z(0)>0 will be referred to.


## 2. Global Positive Solutions

Lemma 1.
For any initial value(x0,y0,z0)∈R+3, where R+3={(x,y,z)∣x>0,y>0,z>0}, system (4) has a unique positive local solution (x(t),y(t),z(t)) for t∈[0,τe] almost surely, where τe is the explosion time.Proof.
Consider the following equations:(5)du=(a1-12σ12-b1eu-c2ev1+α2ev+β2eu-c3ew1+α3ev+β3ew)dt+σ1dB1(t),dv=(a2-12σ22-b2ev-h2evf2+g2ew+d2eu1+α2ev+β2eu)dt+σ2dB2(t),dw=(a3-12σ32-b3ew-h3evf3+g3ew+d3eu1+α3ew+β3eu)dt+σ3dB3(t),
on t>0 with initial value u(0)=ln⁡x(0), v(0)=ln⁡y(0), and w(0)=ln⁡z(0). The coefficients of (5) satisfy the local Lipschitz condition, thus there is a unique local solution (u(t),v(t),w(t)) on [0,τe). Then x(t)=eu(t), y(t)=ev(t), and z(t)=ew(t) are the unique positive local solutions with initial value x0>0, y0>0, and z0>0 by Itô’s formula.Theorem 2.
For any initial value(x0,y0,z0)∈R+3, there is a unique solution (x(t),y(t),z(t)) of system (4) on t⩾0, and the solution will remain in R+3 with probability 1.Proof.
According to Lemma1, we only need to show that τe=∞. Let n0>0 be sufficiently large for x0, y0, and z0 lying within the interval [1/n0,n0]. For each integer n>n0, we define the stopping times
(6)τn=inf⁡{t∈[0,τe):x(t)∉(1n,n),y(t)∉(1n,n)orz(t)∉(1n,n)}.
Obviously, τn is increasing as n→∞. Set τ∞=lim⁡n→∞τn; hence, τ∞⩽τe almost surely. Now, we only need to show that τ∞=∞. If this statement is false, there is a pair of constants T>0 and ϵ∈(0,1) such that P{τ∞⩽T}>ϵ. Thus there exists an integer n1>n0 such that
(7)P{τn⩽T}⩾ϵ,n⩾n1.
Define a C3 function V: R+3→R+ by
(8)V(x,y,z)=(x-1-ln⁡x)+(y-1-ln⁡y)+(z-1-ln⁡z).
The nonnegativity of this function can be seen from u-1-ln⁡u⩾0 and ∀u>0. If (x0,y0,z0)∈R+3, we have
(9)dV=(x-1)(a1-b1x-c2y1+α2x+β2y-c3z1+α3x+β3z)dt+σ122dt+(y-1)(a2-b2y-h2yf2+g2z+d2x1+α2x+β2y)dt+σ222dt+(z-1)(a3-b3z-h3zf3+g3y+d3x1+α3x+β3z)dt+σ322dt+(x-1)σ1dB1(t)+(y-1)σ2dB2(t)+(z-1)σ3dB3(t),
then
(10)LV=(x-1)(a1-b1x-c2y1+α2x+β2y-c3z1+α3x+β3z)+σ122+(y-1)(a2-b2y-h2yf2+g2z+d2x1+α2x+β2y)+σ222+(z-1)(a3-b3z-h3zf3+g3y+d3x1+α3x+β3z)+σ322⩽(x-1)(a1-b1x)+(y-1)(a2-b2y+d2α2)+(z-1)(a3-b3z+d3α3)+σ12+σ22+σ322⩽x(a1-b1x)+y(a2-b2y+d2α2)+z(a3-b3z+d3α3)+σ12+σ22+σ322⩽K,
where K is a positive number. Substituting this inequality into (9), we see that
(11)dV(x,y,z)⩽Kdt+(x-1)σ1dB1(t)+(y-1)σ2dB2(t)+(z-1)σ3dB3(t).
Integrating both sides of the above inequality from 0 to τk∧T and then taking the expectations leads to
(12)EV(x(τn∧T),y(τn∧T),z(τn∧T))⩽V(x0,y0,z0)+KT.
Set Ωn={τn⩽T}, then we get P(Ωn)⩾ϵ by inequality (7). Obviously, for every ω∈Ωn, there are at least x(τn,ω), y(τn,ω), and z(τn,ω) which equal either k or 1/k, then V(x(τn,ω),y(τn,ω),z(τn,ω)) is no less than
(13)min⁡{(n-1-ln⁡n),(1n-1-ln⁡1n)}.
It then follows from (9) that
(14)V(x0,y0,z0)+KT⩾E[1Ωn(ω)V(x(τn,ω),y(τn,ω),z(τn,ω))]⩾ϵmin⁡{(n-1-ln⁡n),(1n-1-ln⁡1n)},
where 1Ωn(ω) is the indicator function of Ωn, letting n→∞, we have that
(15)∞>V(x0,y0,z0)+KT=∞
This completes the proof.


## 3. Stochastic Boundedness

Definition 3.
The solution(x(t),y(t),z(t)) of system (4) is said to be stochastically ultimately bounded, if for any ϵ∈(0,1), there is a positive constant δ=δ(ϵ), such that for any initial value (x0,y0,z0)∈R+3, the solution (x(t),y(t),z(t)) of system (4) has the property that
(16)limsup⁡n→∞P{|x(t),y(t),z(t)|x2(t)+y2(t)+z2(t)=x2(t)+y2(t)+z2(t)>δ}<ϵ.Assumption A.
For any initial value(x0,y0,z0)∈R+3, there exists p⩾1 such that
(17)x(0)<a1+(1/2)(p-1)σ12b1,y(0)<a2+(d2/α2)+(1/2)(p-1)σ22b2,z(0)<a3+(d3/α3)+(1/2)(p-1)σ32b3.Lemma 4.
Assume that Assumption A holds. Let(x(t),y(t),z(t)) be a positive solution of (4) with any initial value (x0,y0,z0)∈R+3, for all p>1, then
(18)E[xp(t)]⩽K1(p),E[yp(t)]⩽K2(p),E[zp(t)]⩽K3(p),
where
(19)K1(p):=(a1+(1/2)(p-1)σ12b1)p,K2(p):=(a2+(d2/α2)+(1/2)(p-1)σ22b2)p,K3(p):=(a3+(d3/α3)+(1/2)(p-1)σ32b3)p.Proof.
Define the functionV1=xp, for (x(t),y(t),z(t))∈R+3 and p>0. By Itô’s formula we get
(20)dV1=px1p-1dx+12p(p-1)x1p-2(dx)2=pxp[(a1-b1x-c2y1+α2x+β2y(a1-b1x-c2y1+α2x+β2y-c3z1+α3x+β3z)-c3z1+α3x+β3z)dt+σ1dB1(t)]+12p(p-1)σ12dt=pxp[(a1-b1x-c2y1+α2x+β2y-c3z1+α3x+β3z-c3z1+α3x+β3z+12p(p-1)σ12)dt+σ1dB1(t)c2y1+α2x+β2y].
Integrating from 0 to t and taking expectations yields
(21)E[xp(t)]-E[xp(0)]=∫0t‍pE[xp(a1-b1x-c2y1+α2x+β2y-c3z1+α3x+β3z+12p(p-1)σ12)]ds.
So,
(22)dE[xp(t)]dt=pE[xp(a1-b1x-c2y1+α2x+β2y-c3z1+α3x+β3z+12p(p-1)σ12)]⩽p{(a1+12(p-1)σ12)E[xp(t)]-b1E[xp+1(t)]}⩽p{(a1+12(p-1)σ12)E[xp(t)]-b1E[xp(t)](p+1)/p}=pE[xp(t)]{[(a1+12(p-1)σ12)]-b1E[xp(t)]1/p}.
Let X(t)=E[xp(t)], then we have
(23)dX(t)dt⩽pX(t)[(a1+12(p-1)σ12)-b1X1/p(t)].
From (17), we know
(24)0<b1X1/p(0)=b1x(0)<a1+12(p-1)σ12,
which by the standard comparison argument shows that
(25)(E[xp(t)])1/p=X1/p⩽a1+(1/2)(p-1)σ12b1,
that is,
(26)E[xp(t)]⩽(a1+(1/2)(p-1)σ12b1)p.
Define the functionV2=yp, for (x(t),y(t),z(t))∈R+3 and p>0. By Itô’s formula we get
(27)dV2=pyp-1dy+12p(p-1)yp-2(dy)2=pyp[(a2-b2y-h2yf2+g2z+d2x1+α2x+β2y)dt(a2-b2y-h2yf2+g2z+d2x1+α2x+β2y)+σ2dB2(t)+12(p-1)σ22dt]=pyp[(a2-b2y-h2yf2+g2z+d2x1+α2x+β2y(a2-b2y-h2yf2+g2z+d2x1+α2x+β2y+12(p-1)σ22)dt+σ2dB2(t)].
Integrating from 0 to t and taking expectations yields
(28)E[yp(t)]-E[yp(0)]=∫0t‍pE[yp(a2-b2y-h2yf2+g2z+d2x1+α2x+β2y)(a2-b2y-h2yf2+g2z+d2x1+α2x+β2y)+12(p-1)σ22]ds.
So,
(29)dE[yp(t)]dt=pE[yp(a2-b2y-h2yf2+g2z+d2x1+α2x+β2y)(a2-b2y-h2yf2+g2z+d2x1+α2x+β2y)+12(p-1)σ22]⩽p{(a2+d2α2+12(p-1)σ22)E[yp(t)](a2+d2α2+12(p-1)σ22)-b2E[yp+1(t)]}⩽p{(a2+d2α2+12(p-1)σ22)E[yp(t)](a2+d2α2+12(p-1)σ22)-b2E[yp(t)](p+1)/p}=pE[yp(t)]{(a2+d2α2+12(p-1)σ22)(a2+d2α2+12(p-1)σ22)-b2E[yp(t)]1/p}.
Let Y(t)=E[yp(t)], then we have
(30)dY(t)dt⩽pY(t)[a2+d2α2+12(p-1)σ22-b2Y1/p(t)].
From (17), we know that
(31)0<b2Y1/p(0)=b2y(0)<a2+d2α2+12(p-1)σ22,
which by the standard comparison argument shows that
(32)(E[yp(t)])1/p=Y1/p(t)⩽a2+(d2/α2)+(1/2)(p-1)σ22b2,
that is,
(33)E[yp(t)]⩽(a2+(d2/α2)+(1/2)(p-1)σ22b2)p.
Similarly, we can show that
(34)E[zp(t)]⩽(a3+(d3/α3)+(1/2)(p-1)σ32b3)p.
This completes the proof.Theorem 5.
Assume that Assumption A holds, the solutions of system (4) with initial value (x0,y0,z0)∈R+3 are stochastically ultimately bounded.Proof.
If(x(t),y(t),z(t))∈R3, its norm here is denoted by |x(t),y(t),z(t)|=(x2(t)+y2(t)+z2(t))1/2, then
(35)|(x(t),y(t),z(t))|p⩽3p/2(|x(t)|p+|y(t)|p+|z(t)|p)
by Lemma 4, E[|x(t),y(t),z(t)|p]⩽K(p), t∈(0,+∞). K(p) is dependent on (x0,y0,z0) and is defined by K(p)=3p/2(K1(p)+K2(p)+K3(p)). By virtue of Chebyshev inequality, we can easily obtain that the solution (x(t),y(t),z(t)) of system (4) is stochastically ultimately bounded.


## 4. Stochastic Permanence

Definition 6 (see [10]).
The solution(x(t),y(t),z(t)) of system (4) is said to be stochastically permanent, if for any ϵ∈(0,1), there exist a pair of positive constants δ=δ(ϵ) and χ=χ(ϵ) such that for any initial value (x0,y0,z0)∈R+3, the solution (x(t),y(t),z(t)) of system (4) has the properties that
(36)lim⁡inf⁡t⇀∞P{|(x(t),y(t),z(t))|⩾δ}⩾1-ϵ,lim⁡inf⁡t⇀∞P{|(x(t),y(t),z(t))|⩽χ}⩾1-ϵ.Assumption B.
One has
1
/
2
(
max
⁡
{
σ
1
,
σ
2
,
σ
3
}
)
2
<
min
⁡
{
a
1
,
a
2
-
(
c
2
/
α
2
)
,
a
3
-
(
c
3
/
α
3
)
}.Theorem 7.
Under Assumption B, for any initial value(x0,y0,z0)∈R+3, the solution (x(t),y(t),z(t)) of system (4) satisfies that
(37)limsup⁡t→∞E[1|(x(t),y(t),z(t))|θ]⩽H,
where θ is an arbitrary positive constant satisfying
(38)θ+12(max⁡{σ1,σ2,σ3})2<min⁡{a1,a2-c2α2,a3-c3α3},
and k is an arbitrary positive constant satisfying
(39)θmin⁡{a1,a2-c2α2,a3-c3α3}-θ(θ+1)2(max⁡{σ1,σ2,σ3})2-k>0.Proof.
DefineV(x,y,z)=x+y+z for (x,y,z)∈R+3, then
(40)dV(x,y,z)=x(a1-b1x-c2y1+α2x+β2y-c3z1+α3x+β3z)dt+y(a2-b2y-h2yf2+g2z+d2x1+α2x+β2y)dt+z(a3-b3z-h3zf3+g3y+d3x1+α3x+β3z)dt+σ1xdB1(t)+σ2ydB2(t)+σ3zdB3(t).
Define U(x,y,z)=(1/(V(x,y,z))) for (x,y,z)∈R+3, by Itô’s formula, we get
(41)dU=-U2[x(a1-b1x-c2y1+α2x+β2y-c3z1+α3x+β3z)+y(a2-b2y-h2yf2+g2z+d2x1+α2x+β2y)+z(a3-b3z-h3zf3+g3y+d3x1+α3x+β3z)]dt+U3(σ12x2+σ22y2+σ32z2)dt-U2(σ1xdB1(t)+σ2ydB2(t)+σ3zdB3(t))=LUdt-U2(σ1xdB1(t)+σ2ydB2(t)+σ3zdB3(t)),
where
(42)LU=-U2[x(a1-b1x-c2y1+α2x+β2y-c3z1+α3x+β3z)+y(a2-b2y-h2yf2+g2z+d2x1+α2x+β2y)+z(a3-b3z-h3zf3+g3y+d3x1+α3x+β3z)]+U3(σ12x2+σ22y2+σ32z2).
Under Assumption B, choosing a positive constantθ such that it satisfies (38). By Itô’s formula, we get
(43)d(1+U)θ=[θ(1+U)θ-1LU+θ(θ-1)2U4(1+U)θ-2θ(θ-1)2×(σ12x2+σ22y2+σ32z2)]dt-θU2(1+U)θ-1(σ1xdB1(t)+σ2ydB2(t)+σ3zdB3(t)),
where
(44)L(1+U)θ=θ(1+U)θ-1LU+θ(θ-1)2U4×(1+U)θ-2(σ12x2+σ22y2+σ32z2),
then choosing a positive constant k such that it obeys (39), by Itô’s formula,
(45)dekt(1+U)θ=ektd(1+U)θ+kekt(1+U)θdt,
where
(46)Lekt(1+U)θ=ektL(1+U)θ+kekt(1+U)θ=ekt(1+U)θ-2×{k(1+U)2+θ(θ+1)2U4(σ12x2+σ22y2+σ32z2)+θU3(σ12x2+σ22y2+σ32z2)-θ(1+U)U2×[x(a1-b1x-c2y1+α2x+β2y-c3z1+α3x+β3z)+y(a2-b2y-h2yf2+g2z+d2x1+α2x+β2y)+z(a3-b3z-h3zf3+g3y+d3x1+α3x+β3z)]}⩽ekt(1+U)θ-2×{k(1+U)2+θ(θ+1)2U4(σ12x2+σ22y2+σ32z2)+θU3(σ12x2+σ22y2+σ32z2)-θU2[(a1x+(a2-c2α2)y+(a3-c3α3)z)-(b1x2+(b2+h2f2)y+(b3+h3f3)z)]-θU3[(a1x+(a2-c2α2)y+(a3-c3α3)z)-(b1x2+(b2+h2f2)y+(b3+h3f3)z)]}⩽ekt(1+U)θ-2×[(k+θmax⁡{b1,b2+h2f2,b3+h3f3})+(2k-θmin⁡{a1,a2-c2α2,a3-c3α3}+θmax⁡{b1,b2+h2f2,b3+h3f3}{a1,a2-c2α2,a3-c3α3}+θ(max⁡{σ1,σ2,σ3})2)U+(k-θmin⁡{a1,a2-c2α2,a3-c3α3}+θ(θ+1)2(max⁡{σ1,σ2,σ3})2)U2],
where
(47)U3(σ12x2+σ22y2+σ32z2)⩽(max⁡{σ1,σ2,σ3})2U,θ(θ+1)2U4(σ12x2+σ22y2+σ32z2)⩽θ(θ+1)2(max⁡{σ1,σ2,σ3})2U2.
Hence, it implies that there exists a positive constant K such that
(48)Lekt(1+U)θ⩽Kekt.
Then we have
(49)dekt(1+U)θ⩽Kektdt-ektθU2(1+U)θ-1×(σ1xdB1(t)+σ2ydB2(t)+σ3zdB3(t)).
Integrating both sides of the above inequality from 0 to t and then taking the expectations leads to
(50)E[ekt(1+U(t))θ]⩽(1+U(0))θ+Kkekt=(1+U(0))θ+H1ekt,
where H1=K/k. So
(51)limsup⁡t→∞E[U(t)θ]⩽limsup⁡t→∞E[(1+U(t))θ]⩽H1.
Since that (x+y+z)θ⩽3θ(x2+y2+z2)θ/2=3θ|x,y,z|θ, where (x,y,z)∈R+3, obviously
(52)limsup⁡t→∞E[1|(x(t),y(t),z(t))|θ]⩽3θlimsup⁡t→∞E[U(t)θ]⩽3θH1=:H,
as required.Theorem 8.
Under Assumption B, system (4) is stochastically permanent.Proof.
By Theorem5, we know that
(53)limsup⁡t→∞E|(x(t),y(t),z(t))|p⩽K(p).
Now, for any ϵ>0, let χ=(K(p)/ϵ)1/p. Then by Chebyshev’s inequality, we can obtain the conclusion easily.


## 5. Global Asymptotic Stability

Definition 9.
Let(x1(t),y1(t),z1(t)) be a positive solution of system (4). If we say that (x1(t),y1(t),z1(t)) is globally asymptotically stable in expectation, it means that any other solution (x2(t),y2(t),z2(t)) of system (4) has t⩾0 and that we have initial value (x0,y0,z0)∈R+3. That is
(54)P{lim⁡t⇀∞E[|(x1(t),y1(t),z1(t))-(x2(t),y2(t),z2(t))|]=0lim⁡t⇀∞E[|(x1(t),y1(t),z1(t))}=1.Lemma 10 (see [11]).
Suppose that an n-dimensional stochastic processX(t) on t⩾0 satisfies the condition
(55)E|X(t)-X(s)|α⩾c|t-s|1+β,0⩽s,t<∞,
for some positive constants α, β, and c. Then there exists a continuous modification X~(t) of X(t) which has the property that for every ϑ∈(0,β/α) there is a positive random variable h(ω) such that
(56)P{ω:sup⁡0<|t-s|<h(ω),0⩽s,t<∞|X~(t,ω)-X~(s,ω)||t-s|ϑ⩽21-2-ϑ:sup⁡0<|t-s|<h(ω),0⩽s,t<∞|X~(t,ω)-X~(s,ω)||t-s|ϑ}=1.
In other words, almost every sample path of X~(t) is locally but uniformly Hölder continuous with exponent ϑ.Lemma 11.
Let(x(t),y(t),z(t)) be a solution of system (4) on t⩾0 with initial value (x0,y0,z0)∈R+3, then almost every sample path of (x(t),y(t),z(t)) is uniformly continuous on t⩾0.Proof.
From system (4), we have the following stochastic integral equation
(57)x(t)=x(0)+∫0t‍x(s)×(a1-b1x(s)-c2y(s)1+α2x(s)+β2y(s)-c3z(s)1+α3x(s)+β3z(s))ds+∫0t‍σ1x(s)dB1(s).
Let  f(s)=x(s)(a1-b1x(s)-(c2y(s)/(1+α2x(s)+β2y(s))) − (c3z(s)/(1+α3x(s)+β3z(s)))), g(s) = σ1x(s), notice that
(58)E|f(t)|p=E[|x(a1-b1x-c2y1+α2x+β2y-c3z1+α3x+β3z)|p]⩽12E|x|2p+12E×[|a1-b1x-c2y1+α2x+β2y-c3z1+α3x+β3z|2p]⩽12E|x|2p+12E[|a1+b1x+c2y+c3z|2p]⩽12E|x|2p+42p-12×(a12p+b12pE|x|2p+c2E|y|2p+c3E|z|2p)⩽12K1(2p)+42p-12×(a12p+b12pK1(2p)+c2K2(2p)+c3K3(2p))=:K4(p).
On the other hand, by the moment inequality (see [12]) for stochastic integrals, we have that for 0⩽t1<t2<∞ and p>2,
(59)E|∫t1t2‍g(s)dB1(s)|p⩽[p(p-1)2]p/2(t2-t1)(p-2)/2∫t1t2‍E|g(s)|pds⩽[p(p-1)2]p/2(t2-t1)(p-2)/2∫t1t2‍σ1pk1(p)ds=[p(p-1)2]p/2(t2-t1)p/2G1(p),
where G1(p)=σ1pk1(p). Let 0⩽t1<t2<∞, t2-t1⩽1, and (1/p)+(1/q)=1, we obtain
(60)E|x(t2)-x(t1)|p=E|∫t1t2‍f(s)ds+∫t1t2‍g(s)dB1(s)|p⩽2p-1E|∫t1t2‍f(s)ds|p+2p-1E|∫t1t2‍g(s)dB1(s)|p⩽2p-1(∫t1t2‍1qds)p/qE(∫t1t2‍|f(s)|pds)+2p-1[p(p-1)2]p/2(t2-t1)p/2G1(p)⩽2p-1(t2-t1)pK4(p)+2p-1[p(p-1)2]p/2(t2-t1)p/2G1(p)⩽2p-1(t2-t1)p/2{(t2-t1)p/2+[p(p-1)2]p/2}K5(p)⩽2p-1{1+[p(p-1)2]p/2}K5(p)(t2-t1)p/2,
where K5(p):=max⁡{K4(p),G1(p)}. Then, we have that almost every sample path of x(t) is locally but uniformly Hölder continuous with exponent ϑ for every ϑ∈(0,(p-2)/2p) and therefore almost every sample path of x(t) is uniformly continuous on t⩾0. Similarly, we can show that almost every sample path of y(t) and z(t) is uniformly continuous on t⩾0.Lemma 12 (see [13]).
Letf(t) be a nonnegative function defined on [0,∞) such that f(t) is integrable on [0,∞) and is uniformly continuous on t⩾0. Then lim⁡t→∞f(t)=0.Theorem 13.
If(61)A:=b1-c2α2β2-c3α3β3-d2-d3>0,B:=b2-c2-2h2f2-d2β2α2-h3g3(a2+d2/α2+σ22)b2f32>0,C:=b3-c3-2h3f3-d3β3α3-h2g2(a3+d3/α3+σ32)b3f22>0,
then system (4) is globally asymptotically stable.Proof.
Define(62)V(t)=|ln⁡x1(t)-ln⁡x2(t)|+|ln⁡y1(t)-ln⁡y2(t)|+|ln⁡z1(t)-ln⁡z2(t)|,
then V(t) is a continuous positive function on t⩾0. A direct calculation of the right differential d+V(t) of V(t), and then applying Itô’s formula, we have
(63)d+V(t)=sgn⁡(x1(t)-x2(t))×{[dx1(t)x1(t)-(dx1(t))22x12(t)]-[dx2(t)x2(t)-(dx2(t))22x22(t)]}+sgn⁡(y1(t)-y2(t))×{[dy1(t)y1(t)-(dy1(t))22y12(t)]-[dy2(t)y2(t)-(dy2(t))22y22(t)]}+sgn⁡(z1(t)-z2(t))×{[dz1(t)z1(t)-(dz1(t))22z12(t)]-[dz2(t)z2(t)-(dz2(t))22z22(t)]}=sgn⁡(x1(t)-x2(t))×[-b1(x1(t)-x2(t))-c3z2(t)1+α3x2(t)+β3z2(t))-(c2y1(t)1+α2x1(t)+β2y1(t)-c2y2(t)1+α2x2(t)+β2y2(t))-(c3z1(t)1+α3x1(t)+β3z1(t)-c3z2(t)1+α3x2(t)+β3z2(t))]dt+sgn⁡(y1(t)-y2(t))×[-b2(y1(t)-y2(t))-d2x2(t)1+α2x2(t)+β2y2(t)-(h2y1(t)f2+g2z1(t)-h2y2(t)f2+g2z2(t))+(d2x1(t)1+α2x1(t)+β2y1(t)-d2x2(t)1+α2x2(t)+β2y2(t))]dt+sgn⁡(z1(t)-z2(t))×[-b3(z1(t)-z2(t))-d3x2(t)1+α2x2(t)+β2z2(t)-(h3z1(t)f3+g3y1(t)-h3z2(t)f3+g3y2(t))+(d3x1(t)1+α3x1(t)+β3z1(t)-d3x2(t)1+α2x2(t)+β2z2(t))]dt.
Integrating from 0 to t and taking expectations yields
(64)E[V(t)]-E[V(0)]=E[c2y1(s)1+α2x1(s)+β2y1(s)∫0t‍[-d3x2(s)1+α2x2(s)+β2z2(s)))sgn⁡(x1(s)-x2(s))×(-b1(x1(s)-x2(s))-c3z2(s)1+α3x2(s)+β3z2(s)-(c2y1(s)1+α2x1(s)+β2y1(s)-c2y2(s)1+α2x2(s)+β2y2(s))-(c3z1(s)1+α3x1(s)+β3z1(s)-c3z2(s)1+α3x2(s)+β3z2(s)))+sgn⁡(y1(s)-y2(s))×(-b2(y1(s)-y2(s))d2x2(s)1+α2x2(s)+β2y2(s)-(h2y1(s)f2+g2z1(s)-h2y2(s)f2+g2z2(s))+(d2x1(s)1+α2x1(s)+β2y1(s)-d2x2(s)1+α2x2(s)+β2y2(s)))+sgn⁡(z1(s)-z2(s))×(-b3(z1(s)-z2(s))+sgn⁡(z1(t)-z2(t))+sgn⁡(z1(t)-z2(t))d3x2(s)1+α2x2(s)+β2z2(s)-(h3z1(s)f3+g3y1(s)-h3z2(s)f3+g3y2(s))+(d3x1(s)1+α3x1(s)+β3z1(s)-d3x2(s)1+α2x2(s)+β2z2(s)))]ds].
So
(65)dE[V(t)]dt=E[sgn⁡(x1(t)-x2(t))d3x2(t)1+α2x2(t)+β2z2(t)×(-b1(x1(t)-x2(t))d3x2(t)1+α2x2(t)+β2z2(t)-(c2y1(t)1+α2x1(t)+β2y1(t)-c2y2(t)1+α2x2(t)+β2y2(t))-(c3z1(t)1+α3x1(t)+β3z1(t)-c3z2(t)1+α3x2(t)+β3z2(t)))+sgn⁡(y1(t)-y2(t))×(-b2(y1(t)-y2(t))d3x2(t)1+α2x2(t)+β2z2(t)-(h2y1(t)f2+g2z1(t)-h2y2(t)f2+g2z2(t))+(d2x1(t)1+α2x1(t)+β2y1(t)-d2x2(t)1+α2x2(t)+β2y2(t)))+sgn⁡(z1(t)-z2(t))×(-b3(z1(t)-z2(t))d3x2(t)1+α2x2(t)+β2z2(t)-(h3z1(t)f3+g3y1(t)-h3z2(t)f3+g3y2(t))+(d3x1(t)1+α3x1(t)+β3z1(t)-d3x2(t)1+α2x2(t)+β2z2(t)))]⩽-b1E[|x1(t)-x2(t)|]-c2E[sgn⁡(x1(t)-x2(t))d3x2(t)1+α2x2(t)+β2z2(t)×(y1(t)1+α2x1(t)+β2y1(t)-y2(t)1+α2x2(t)+β2y2(t))]-c3E[sgn⁡(x1(t)-x2(t))d3x2(t)1+α2x2(t)+β2z2(t)×(z1(t)1+α3x1(t)+β3z1(t)-z2(t)1+α3x2(t)+β3z2(t))]-b2E[|y1(t)-y2(t)|]-h2E[sgn⁡(y1(t)-y2(t))d3x2(t)1+α2x2(t)+β2z2(t)×(y1(t)f2+g2z1(t)-y2(t)f2+g2z2(t))]+d2E[sgn⁡(y1(t)-y2(t))d3x2(t)1+α2x2(t)+β2z2(t)×(x1(t)1+α2x1(t)+β2y1(t)-x2(t)1+α2x2(t)+β2y2(t))]-b3E[|z1(t)-z2(t)|]-h3E[sgn⁡(z1(t)-z2(t))d3x2(t)1+α2x2(t)+β2z2(t)×(z1(t)f3+g3y1(t)-z2(t)f3+g3y2(t))]+d3E[sgn⁡(z1(t)-z2(t))d3x2(t)1+α2x2(t)+β2z2(t)×(x1(t)1+α3x1(t)+β3z1(t)-x2(t)1+α2x2(t)+β2z2(t))]⩽-b1E[|x1(t)-x2(t)|]+c2E[|y1(t)1+α2x1(t)+β2y1(t)-y2(t)1+α2x2(t)+β2y2(t)|]+c3E[|z1(t)1+α3x1(t)+β3z1(t)-z2(t)1+α3x2(t)+β3z2(t)|]-b2E[|y1(t)-y2(t)|]+h2E[|y1(t)f2+g2z1(t)-y2(t)f2+g2z2(t)|]+d2E[|x1(t)1+α2x1(t)+β2y1(t)-x2(t)1+α2x2(t)+β2y2(t)|]-b3E[|z1(t)-z2(t)|]+h3E[|z1(t)f3+g3y1(t)-z2(t)f3+g3y2(t)|]+d3E[|x1(t)1+α3x1(t)+β3z1(t)-x2(t)1+α2x2(t)+β2z2(t)|]⩽-b1E[|x1(t)-x2(t)|]+c2E[|y1(t)-y2(t)|]+c2α2β2E[|x1(t)-x2(t)|]+c3E[|z1(t)-z2(t)|]+c3α3β3E[|x1(t)-x2(t)|]-b2E[|y1(t)-y2(t)|]+h2f2E[|y1(t)-y2(t)|]+h2g2f22E[y1(t)]E[|z1(t)-z2(t)|]+h2f2E[|y1(t)-y2(t)|]+d2E[|x1(t)-x2(t)|]+d2β2α2E[|y1(t)-y2(t)|]-b3E[|z1(t)-z2(t)|]+h3f3E[|z1(t)-z2(t)|]+h3g3f32E[z1(t)]E[|y1(t)-y2(t)|]+h3f3E[|z1(t)-z2(t)|]+d3E[|x1(t)-x2(t)|]+d3β3α3E[|z1(t)-z2(t)|]⩽(-b1+c2α2β2+c3α3β3+d2+d3)×E[|x1(t)-x2(t)|]+(-b2+c2+2h2f2+d2β2α2+h3g3f32(E[z13(t)])1/3)×E[|y1(t)-y2(t)|]+(-b3+c3+2h3f3+d3β3α3+h2g2f22(E[y13(t)])1/3)×E[|z1(t)-z2(t)|].
By Lemma 4(66)(E[x13(t)])1/3⩽a1+σ12b1,(E[y13(t)])1/3⩽a2+(d2/α2)+σ22b2,(E[z13(t)])1/3⩽a3+(d3/α3)+σ32b3.
Thus
(67)dE[V(t)]dt⩽(-b1+c2α2β2+c3α3β3+d2+d3)E[|x1(t)-x2(t)|]+(-b2+c2+2h2f2+d2β2α2+h3g3(a2+d2/α2+σ22)b2f32)×E[|y1(t)-y2(t)|]+(-b3+c3+2h3f3+d3β3α3+h2g2(a3+d3/α3+σ32)b3f22)×E[|z1(t)-z2(t)|]=AE[|x1(t)-x2(t)|]+BE[|y1(t)-y2(t)|]+CE[|z1(t)-z2(t)|].
Integrating both sides leads to
(68)E[V(t)]⩽E[V(0)]+∫0t‍AE[|x1(s)-x2(s)|]+BE[|y1(s)-y2(s)|]+CE[|z1(s)-z2(s)|]ds.
Therefore
(69)0⩽E[V(t)]+∫0t‍AE[|x1(s)-x2(s)|]+BE[|y1(s)-y2(s)|]+CE[|z1(s)-z2(s)|]ds⩽E[V(0)]<∞.
Then we get
(70)E[|(x1(t),y1(t),z1(t))-(x2(t),y2(t),z2(t))|]=E[+|z1(t)-z2(t)|2)1/2(|x1(t)-x2(t)|2+|y1(t)-y2(t)|2+|z1(t)-z2(t)|2)1/2]⩽E[|x1(t)-x2(t)|]+E[|y1(t)-y2(t)|]+E[|z1(t)-z2(t)|]∈L1[0,∞).
Therefore from Lemmas 11 and 12, we have
(71)lim⁡t→+∞E[|(x1(t),y1(t),z1(t))-(x2(t),y2(t),z2(t))|]=0.
This completes the proof of Theorem 13.


## 6. Extinction

Theorem 14.
For any initial value(x0,y0,z0)∈R+3, the solution (x(t),y(t),z(t)) of system (4) obeys(72)limsup⁡t→∞ln⁡x(t)t⩽a1-σ122,limsup⁡t→∞ln⁡y(t)t⩽a2+d2α2-σ222,limsup⁡t→∞ln⁡z(t)t⩽a3+d3α3-σ322.Proof.
Define Lyapunov functionsln⁡x, ln⁡y, and ln⁡z, respectively, hence by Itô’s formula, we get
(73)ln⁡x(t)=ln⁡x0+(a1-σ122)t+∫0t‍-b1x(s)-c2y(s)1+α2x(s)+β2y(s)-c3z(s)1+α3x(s)+β3z(s)ds,ln⁡y(t)=ln⁡y0+(a2-σ222)t+∫0t‍-b2y(s)-h2y(s)f2+g2z(s)ds+∫0t‍d2x(s)1+α2x(s)+β2y(s)ds,ln⁡z(t)=ln⁡z0+(a3-σ322)t+∫0t‍-b3z(s)-h3z(s)f3+g3y(s)ds+∫0t‍d3x(s)1+α3x(s)+β3z(s)ds.
Then we have
(74)ln⁡x(t)⩽ln⁡x0+(a1-σ122)t+σ1B1(t),ln⁡y(t)⩽ln⁡y0+(a2+d2α2-σ222)t+σ2B2(t),ln⁡z(t)⩽ln⁡z0+(a3+d3α3-σ322)t+σ3B3(t).
Dividing t on the both sides and letting t→∞, we can derive
(75)limsup⁡t→∞ln⁡x(t)t⩽a1-σ122,(76)limsup⁡t→∞ln⁡y(t)t⩽a2+d2α2-σ222,(77)limsup⁡t→∞ln⁡z(t)t⩽a3+d3α3-σ322,
as required.So we can obtain that ifa1-(σ12/2)<0, a2+(d2/α2)-(σ22/2)<0, and a3+(d3/α3)-(σ32/2)<0 hold, then for any initial value (x0,y0,z0)∈R+3, the solution (x(t),y(t),z(t)) of system (4) will be extinct exponentially with probability one.


## 7. Numerical Simulations

At last, we numerically simulate the solution(x(t),y(t),z(t)) of system (4) to substantiate the analytical findings. By the method mentioned in [14], we consider the discretization equation:
(78)xk+1=xk+xk×(a1-b1xk-c2yk1+α2xk+β2yk-c3zk1+α3xk+β3zk)Δt+σ1xkξkΔt+12σ12xk(ξk2-1)Δt,yk+1=yk+yk×(a2-b2yk-h2ykf2+g2zk+d2xk1+α2xk+β2yk)Δt+σ2ykηkΔt+12σ22yk(ηk2-1)Δt,zk+1=zk+zk×(a3-b3zk-h3zkf3+g3yk+d3xk1+α3xk+β3zk)Δt+σ3zkζkΔt+12σ32zk(ζk2-1)Δt,
where ξk, ηk, and ζk are Gaussian random variables that follow N(0,1).For example, in Figure1, we choose initial value (x0,y0,z0)=(0.5,0.5,0.5) and parameters satisfying conditions of Theorem 7; the system is stochastically permanent. In Figure 2, (a) is the solution to the deterministic system (4), (b) is the solution to the stochastic system (4), and (c) is the phase diagram of the stochastic system (4). Comparing (a) with (b), we find that with small environmental noise, the stochastic system is getting more similar to the deterministic. In Figure 3, (a) is the solution to the deterministic system (4) and (b) is the solution to the stochastic system (4). Comparing (a) with (b), we find that the sufficiently large environmental noise make the stochastic system extinct. The phase diagram of (x,y,z) is displayed in (c).Figure 1
Solutions of system (4) for (x0,y0,z0) = (0.5,0.5,0.5),a1 = 2;a2=1;a3 = 1.5;b1=3;b2=2;b3=2;c2 = 0.2;c3=0.1;d2=1;d3 = 1.5;α2=1;β2=0.5;α3=0.6;β3=0.8;h2=0.5;h3=1;f2=1;f3=0.3;g2 = 0.5;g3=1;σ1=0.05;σ2 = 0.05;σ3=0.05.Solutions of system (4) for (x0,y0,z0) = (0.6,0.5,0.5),a1=2;a2=1;a3 = 1.5;b1=2.5;b2 = 2;b3=1.5;c2=0.2;c3 = 0.1;d2=1.5;d3=1;α2 = 1;β2=0.5;α3=1;β3=0.8;h2=0.5;h3=1;f2=1;f3=0.3;g2=1;g3=0.5;σ1=0.05;σ2 = 0.05;σ3=0.05.
(a)
(b)Solutions of system (4) for (x0,y0,z0)=(1.5,1,1),a1 = 1.1;a2=1.2;a3=1.5;b1=0.8;b2=1.1;b3=1.2;c2=0.02;c3 = 0.01;d2=1.2;d3=1;α2=0.7;β2=0.5;α3 = 0.8;β3=0.5;h2=0.8;h3=0.5;f2=1;f3=0.3;g2=0.8;g3=0.5;σ1 = 1.8;σ2=1.8;σ3=1.8.
(a)
(b)
(c)


## 8. Discussion

This paper has been devoted to dynamics of a stochastic cooperative predator-prey system with Beddington-DeAngelis functional. Firstly we show that, although the coefficients in the model neither satisfy the linear growth condition nor local Lipschitz continuity, the stochastic model has a globally positive solution. Then we know that the positive solution is stochastically bounded. Moreover, under some conditions, we analyze global asymptotic stability of the positive solutions. We can find that the stochastic model will preserve this nice property provided that the noise is sufficiently small. Some meaningful questions deserve further investigation. One way we can consider colored noise in the models owing to sudden environmental changes caused by seasons or other reasons. Moreover, it is worth to rebuild our model with some parameters, not the ones studied above, which are also subject to stochastic excitation. The ecosystem we considered has some limitations, for example, both predators(y and z) with the common prey (x) and we only established a model with a mutual cooperation between both predators. It is worth to discuss mechanics of the population contained predators with a mutual competition in other ecosystems.



---

*Source: 927217-2014-07-13.xml*