0Liu and Wei Advances in Difference Equations (2016) 2016:21 DOI 10.1186/s13662-015-0698-x
0 Advances in Difference Equations
a SpringerOpen Journal
RESEARCH
Open Access
Dynamics of a stochastic cooperative predator-prey system with Beddington-DeAngelis functional response
CrossMark
Xiaoyan Liu and Yiqiang Wei*
Correspondence: wei_yiqiang@163.com College of Mathematics, Taiyuan University of Technology, 79 Yingze W St., Taiyuan, Shanxi 030024, China
Abstract
In this paper, a stochastic cooperative predator-prey system with Beddington-DeAngelis functional response is proposed. By constructing Lyapunov functions, we demonstrate the existence of global positive solution, which is stochastically bounded and permanent. In addition, under some specific conditions, the solution of this stochastic system is globally asymptotically stable, which means that the properties of the solution will not be changed when the stochastic perturbation is small. Finally, the influence of the stochastic perturbation is demonstrated by simulation.
Keywords: stochastic differential equation; Beddington-DeAngelis functional response; stochastic perturbation; global asymptotic stability
1 Introduction
As an important branch of ecology science, population ecology has become a systematic discipline where mathematics is thoroughly applied. The relationship between predator and prey is one of the basic relationships among species. Many significant functional responses are constructed to model various situations. Most of the functional response are prey-dependent, which fail to model the interference among predators. In fact, when predators have to search, share, and compete for food, the functional response should be predator dependent; such a scenario occurs more frequently in nature and laboratory (see [1] and references therein).
In predator-prey systems, there are three classical predator-dependent functional responses: the Hassell-Varley, Beddington-DeAngelis, and Crowley-Martin responses (see [1] and references therein). By comparing statistical evidence from 19 predator-prey systems, Skalski and Gilliam claimed that the above three predator-dependent functional responses provided a better description of a predator feeding over a range of predator-prey abundances [2]. Amongst the Beddington-DeAngelis type functional responses some cases fitted better [2]. Beddington and DeAngelis et al. in 1975 first introduced the Beddington-DeAngelis type predator-prey model (see [3] and references therein), which is
dt =t(ri ailX l+fitly y )'
ÉL = y(_r„ + a2i* )
dt y( '2 + l+fix+yy ''
ft Spri
ringer
© 2016 Liu and Wei. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use,distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
where x = x(t) and y = y(t) represent prey and predator densities, respectively; r1 stands for the average growth rate of prey, rl is the death rate of predator, and an is the predator density-dependence rate. ri, aj, ¡3, and y are positive constants for i,j = 1, I. Because of its practical significance, many scholars have studied the system in recent years [4-11].
In the study of predator-prey systems, the cooperation element is always neglected. Nevertheless, the cooperative system is a sometimes rudimentary and sometimes important ecological system in mathematical biology [13,14]. To describe the mutual cooperation between two species, May [15] proposed the following equations:
, dx = r1x(1- aA-y- c1x), (l)
, d = riy(1- ai+biy- ciy),
where ri represents the growth rate, and ci = j (Ki is the carrying capacity). ri, ai, bi, and ci are positive constants for i, j = 1, I.
For each predator, various species of prey coexist in nature, so the multiple species predator-prey system simulates the actual situation better. However, this has been neglected in many existing studies [1,16]. In this paper, we investigate the following cooperative predator-prey system of one predator feeding on two prey species with Beddington-DeAngelis functional response:
dx = x(a1 - b1x - f+gy - ) dt,
dy = y(ai- b-y -fhy- -I_ll+3-) ^ (3)
dz = z(a3 - b3z + T+Odxkz + 1+all+3iz) dt,
where the species x, y are the preys of z; x and y are cooperative species. ai stands for the average growth rate, and bi is the density-dependence rate. All the parameters in system (3) are positive constants.
In reality, a population system is inevitably affected by environmental perturbation. May [15] have claimed that owing to the environmental perturbation, the birth rates in the system should be stochastic. At the same time, the natural growth of many species vary with t, e.g. owing to the seasonality. However, environmental fluctuations have been neglected in many existing studies. References [1, 3,17] consider the effect of environmental noise. We introduce stochastic perturbation into the intrinsic growth rate. The intrinsic growth rate can be written as an average growth rate with some small random perturbed terms. In general, by the central limit theorem, the small terms follow some normal distributions, thus we can approximate the error term by a white noise aiBi(t), where of is the intensity of the noise and Bi(t) is a standard white noise. Bi(t) is a Brownian motion defined on a complete probability space (fi, F, P). The growth rates ai, al, and a3 are disturbed to ai + a1BB 1(t), al + alBl(t), a3 + a3B3(t), respectively. Finally we obtain the following stochastic system:
dx = x(a1 - b1x - f+gxy - 1+ax+Az) dt + a1xdB1(t),
dy = y(ai - b-y - f+glx - ) dt + aiydBi(D, (4)
dz = z(a3 - b3z + T+Oxkrz + 1+udy+3lz) dt + a3zdB3(t),
where b3z denotes the density dependence of the predator population. All the parameters in system (4) are positive constants.
In this paper, we will study the stochastic predator-prey system of one predator feeding on two prey species with Beddington-DeAngelis functional response. The rest of this paper is organized as follows. Section 2, we show the existence of global positive solution. Section 3, the stochastic boundedness of solution is studied. In Section 4, we prove that the system is stochastically permanent. In Section 5, we investigated the global attractivity of system (4). Finally, in Section 6, we present numerical simulation to verify our analytical results.
2 Global positive solutions
In model (4), as x(t), y(t), z(t) represent predator and prey densities, the solutions of model (4) should be non-negative. In order for a stochastic differential equation to have a unique global solution for any given initial value, the coefficients of the equation are generally required to satisfy the linear growth condition and the local Lipschitz condition (see [18]), whereas the coefficients of (4) neither comply with the linear growth condition nor with the local Lipschitz condition. In the following sections, we will construct Lyapunov functions to demonstrate the existence and uniqueness of positive solutions.
Lemma 1 For any initial value xo > 0, yo > 0, zo > 0, system (4) has a unique positive local solution (x(t), y(t), z(t)) for t e [0, re) almost surely (a.s.), where xe is the explosion time.
Proof Consider the following equations:
du = a - al- be -fhguv - irafW^w) dt + a! ^(t),
dv = (a2-a4-b2ev-f+2u-17oCw) dt + dB2(t) (5)
dw =(a3 - # - b3ew + + ) dt + a3 dB3(t)
on t > 0 with initial value u(0) = lnx0, v(0) = lny0, w(0) = lnz0. The coefficients of (5) satisfy the local Lipschitz condition, so there is a unique local solution (u(t), v(t), z(t)) on t e [0, Te). By Ito's formula, we can see that x(t) = eu(t), y(t) = ev(t), z(t) = ew(t) is the unique positive local solution to (5) with initial value x0 > 0, y0 > 0, z0 > 0. □
Theorem 1 Consider system (4), for any given initial value (x0,y0, z0) e R+, there is a unique solution (x(t),y(t), z(t)) for t > 0 and the solution will remain in R+ with probability 1.
Proof Based on Lemma 1, we only need to show the xe = to. Let m0 > 0 be sufficiently large for x0, y0, z0 lying within the interval [m■, m0]. For each integer m > m0, define the stopping times
Tm = inf jt e [0, Te]: min{x(t),y(t),z(t)} < — or max{x(t),y(t),z(t)} > m^. (6)
Clearly, t— is increasing as m ^to. Set tto = limm^TO t—, whence tto< t—. If we can prove = to, then Te = to and (x0, y0, z0) e R+ a.s. for all t > 0. If this statement is false, there is
a pair of constants T > 0 and e e (0,1) such that P[xm < T} > e. So there exists an integer m1 > m0 such that
P[Tm < T}>e, m > m1. (7)
Define V(x,y, z) = (x - 1 - lnx) + (y -1 - lny) + (z - 1 - lnz). Since u - 1 - ln u > 0 for all u > 0, we see that V(x, y, z) is non-negative. Applying Ito's formula, we have
nr , r h\x c1 z \ , al
dV = (x-1) a1 - b1x —---dt + — at
\ J1+ g1y 1 + «1x + ftz/ I
hiy ciz aii
+ (y -1) a- - biy —----— )dt + — dt
\ Ji+ gix 1 + a-y + ¡izj I
d1x diy a3i
+ (z -1) a3 - b3z + --— + --— )dt + — dt
\ 1 + a1x + ¡¡1z 1 + a2y + p2z/ I
+ a1 (x -1) dB1 (t) + ai(y - 1) dBi(t)+ a3(z -1) dB3(t), (8)
h1x c1z
LV = (x -1)1 a1 - b1x -
+ (y -1) ^ a- - b-y -+ (z -1) ( a3 - b3z +
J1+ g1y 1 + a1x + ¡1z hiy ciz
fi+ gix 1 + a-y + ¡izy
d1x diy \ a1i + a2l + a|
< x(a1 - b1x) + y(a2 - b2y) + z( a3 - b3z +
1 + a1x + ¡1z 1 + a2y + ¡2z/ 2
d1x d2y
1 + a1x + ¡1z 1 + a2y + ¡¡2z
222 + a12 + a22 + a32
, , , , , , / , d1 dA a12 + a2 + a32 < x(a1 - b1x) + y(a2 - b2y) + ^ a3 - b3z + — + — I +-
«1 «2,
< K, (9)
where K is a positive number. Consequently,
dV < Kdt+(x - 1)a1 dB1(t) + (y - 1)a2 dB2(t) + (z - 1)a3 dB3(t). (10)
Integrating both sides of the inequality from 0 to xm A T and then taking the expectations, we get
EV(x(tm A T),y(Tm A T), z(tm A T)) < V(x0,y0, z0) + KT. (11)
Set = [rm < T} for m > m1, then we get P(fim) > e by (7). For every w e , there is at least one of x ( Tm, w), y(Tm, w), z(Tm, w), which equals either m or m, therefore V(x(xm A T),y(Tm A T),z(Tm A T)) is no less than m -1 - lnm or m -1-ln -m;. Consequently, we have
V(x(tm A T),y(Tm A T),z(tm A T)) > (m -1-lnm) A\ — - 1 - ln — ). (12)
It follows from (11) that
V(x0,yo, zo) + KT > E[lam V(x(xm A T),y(xm A T), z(xm A T))]
> e(m — 1 — lnm) a( —-1-ln —), (13)
\m m /
where is the indicator function of Qm. Letting m ^ro, it results in the contradiction that ro > V(x0,y0, z0) + KT = ro, thus we draw a conclusion that Tro = ro a.s. □
Theorem 1 only tells us the solution of model (4) will remain in R3 with probability 1. Next, we will discuss how the solution varies in R+ more detail.
3 Stochastic boundedness
Definition 1 (see [3] and references therein) The solution (x(t),y(t), z(t)) of system (4) is said to be stochastically ultimately bounded, if for any e e (0,1), there is a positive constant S = 5(e).Suchthatforanyinitialvalue(x0,y0,z0) e R+,thesolution(x(t),y(t),z(t)) ofsystem (4) has the property that
limsupP{|x(t),y(t),z(t) | = ^x2(t) + y2(t) + z2(t)> 5} < e. (14)
Assumption 1 For any initial value (x0,y0, z0) e R+, there existsp > 1 such that
al + p-1a12 , N
x0 <-, (15)
a2 + , ,
y0 <-T~ , (6)
d1 d2 P-1 O2
a3 + —i---+ -^-On
„ . 3 a1 a2 2 3 /1-7-1
z0 <-:-. (7)
Lemma 2 If Assumption 1 holds, let (x(t), y(t), z(t)) be a solution to system (4) with initial value (x0,y0,z0) e R+. Then, for allp >1,
E[xp(t)] < K1(p), (18)
E[/(t)] < K2(p), (19)
E[zp(t)] < K3(p), (20)
( a1 + p-1o12 \p
*i(p):=( 1 bi2 1 ] , (1) ( a2 + p-1o22\p
K2(p) := 2 b2 2j , (2)
( a3 + ^ + ^ + p-1o32 \P K3(p):=( 3 a1 ba2 2 M . (3)
Proof Define V1 = xp for x e R+ andp > 0. By Ito's formula, we have 1
dV1 = pxp-1 dx + -p(p - 1)xp-2(dx)2
= pxP-1
x( a1 - b1x-
J1 + g1y 1 + «1x + ¡1 z
dt + a1xdB1(t)
+ 2P(P - 1)xpa12 dt
a1 - b1x -
c1z 1 2
+ 7t(p -1)a12
dt + a1 dB1(t)[. (24)
J1 + g1y 1 + «1x + ¡1z 2 Integrating both sides of the equality from 0 to t and then taking the expectations, we get
E[xp(t)] - E[xp(0)]
a1 - b1x -
h1x c1z 1 2
+ -(p -^a2
J1 + g1y 1 + «1x + ¡1z 2
dE[xP(t)] dt
= pE xp
a1 - b1x -
J1 + g1y 1 + a1x + ¡1z ' 2
c1z 1 2
+ 7t(p -1)a2
P - b1xP+1 + -(p - 1)a12xp
a1 + t(p -1)a12
E[xp(t)] - b1E[xP+1(t)]J
:pE[xP(t)]j
a1 + 1(p - 1)a2 - b1E[xP(t)]p [.
Let X(t) = E[xp(t)], then
^ = PX(t)ia1 + 2(p - 1)a2 - b1Xp(t) [.
By Assumption 1, we know 0 < b1Xp (0) = b1x(0) < a1 + P-1 a12. Applying the standard comparison argument,
E[xp(t)]p < ■
E[xp(t)] < K1(p) =
a1 + P-1 a2 \p
Similarly, we can prove that
E[y(t)] <
a2 + P-1 a22 xp
E[zp(t)] <
d.2 U2
Theorem 2 Assume that Assumption 1 holds, the solutions of system (4) with initial value (x0, yo, z0) e R+ are stochastically ultimately bounded.
Proof If (x(t),y(t),z(t)) e R+, its norm here is denoted by
|x(t), y(t), z(t) | = Vx2(t)+y2(t)+z2(t). Then |x(t), y(t), z(t)|p = 3p(|xp(t)| + |yp(t)| + |zp (t)|). By Lemma 2, E[|x(t), y(t), z(t)|p] < L(p),
t > 0.L(p) is dependent on (x0,y0,z0) e R+ and defined by L(p) = 32 K1(p)+K2(p)+K3(p)]. By Chebyshev's inequality, the above result is straightforward. □
4 Stochastic permanence
Definition 2 (see [3] and references therein) The solution (x(t), y(t), z(t)) of system (4) is said to be stochastically permanent, if for any e e (0,1), there is a positive constant S = 5(e) and x = X (e), such that for any initial value (x0,y0, z0) e R+, system (4) has the properties that
liminfP{|x(t),y(t),z(t) I > 5} > 1 - e, (2)
t—1 '
liminf P{ |x(t), y(t), z(t)| < x} > 1 - e. (3)
Assumption 2
1 i 2 2 2 1 i ci c2 1 . . - max a-f, ai, ai\ < mim a1, a2, a3 + — + — >. 04)
2 ( a a2 J
Theorem 3 If Assumption 2 holds, for any initial value (x0,y0, z0) e R+, the solution (x(t), y(t), z(t)) sc stant satisfying
(x(t), y(t), z(t)) satisfies limsupt—TO E[ |x(t) y(1) z(t) ] < H, where 0 is an arbitrary positive con-
0 + 3 ( 2 2 21 i c1 c2 I
-maxaf, ai, a32! < mim a1, a2, a3 + — + — >, (35)
2 ( a a2 J
where k is an arbitrary positive constant satisfying
0(0 + 3) . 2 2 91 i c1 c2 1
k +-max a1, a2, a3 ! < 0 mim a1, a2, a3 + — + — k (36)
2 ( a a2 J
Proof Define U(x, y, z)= x + y + z. Applying Ito's formula, we have
(h^x c1Z \
a1 - b1x —---I dt + a1xdB1(t)
f1+ g1y 1 + a1x + ftz/
( , h2y c2z \ , , . .
+ yU2 - b2y ------— dt + a2ydB2(t)
\ f2+ g2x 1 + a2y + P2Z/
+ z(a3- b3Z + --d1x -d2y )dt + a3zdB3(t). (37)
\ 1 + a1x + p1z 1 + a2y + p2z)
Define V(x,y,z) = Utyz). Using Ito's formula, we obtain
hix ciz
dV (x, y, z) = -V2
x( a1 - b1x -
+ y^2 - b2y -+ z( a3 - b3z +
f1+ g1y 1 + a1x + ftz h2y C2z
f2+ g1x 1 + a2y + ftz
d1x d2y
1 + a1x + fa1z 1 + a2y + fa2z + V3(o2x2 + o22y2 + o32z2) dt - V2(o1xdB1(t) + O2ydB2(t) + O3zdB3(t)),
LV = -V
xya1 - b1 x -+ y(^2 - b2y -+ z( a3 - b3z +
f1+ g1y 1 + a1x + ftz h2y C2z
f2+ g1x 1 + a2y + fa2z
d1x d2y
1 + a1x + fa1z 1 + a2y + fa2z
, 1/3/ 2„2 , „2.2 .
+ V (o1 x + o2y + o3 z j.
Under Assumption 2, choosing a positive constant 0 such that it satisfies (35), then by Ito's formula, we have
d(1 + V)0 = 0(1 + V)0-1 dV + 10(0 -1)(1 + V)0-2(dV)2
= |0(1 + V)0-1LV + ^0(0 - 1)(1 + V)0-2 V4(o2x2 + o22y2 + o32z2) j dt - 0 V2(1 + V)0-1( O1xdB1(t) + O2ydB2(t) + O3zdB3(t)).
L(1 + V)0 = 0(1 + V)0-1LV + ^0(0 -1)(1 + V)0-2 V4(o2x2 + o22y2 + o32z2).
Then choosing a positive constant k such that it satisfies (36). By Ito's formula, we get
dekt (1 + V)0 = ektd(1 + V)0 + kek (1 + V)0 dt
= ektL(1 + V)0 dt + kekt (1 + V)0 dt - ekt 0 V2(1 + V)0 -1[oixdB1(t) + O2ydB2(t) + 03 zdB3(t)],
Lekt (1 + V)0
= ekt(1 + V)0 Ak(1 + V)2 + 10(0 - 1) V4(o12x2 + o22y2 + o32z2)
+ 0 (1 + V )V 3(o12x2 + o22y2 + o32z2)
- 0 (1 + V )V2 + y^a2 - b2y -+ z( a3 - b3z +
x( a1 - b1x -h2y
f1+ g1y 1 + a1x + ftz
f2+ g2x 1 + a2 y + foz
d1x d2y
1 + a1x + fa1z 1 + a2y + fa2z
< ekt(1 + V)0"2Jk(1 + V)2 + ^0(0 + 1)V4max{o12,o22,o32}(x2 + y2 + z2)
+ 0 (1 + V )V3 max{ o12,o22,o32}(x2 + y2 + z2)
h1x C1z
- 0 (1 + V )V2
x( a1 - b1x -
f1 a1x
+ y ( a2 - b2y - h2y - — ) + z(a3 - b3z)
f2 a2y.
< ekt(1 + V)02jk(1 + V)2 + ±0(0 + 1)V4ma^jo12,o22,o32}(x2 + y2 + z2) + 0 (1 + V )V3 max{ o12,o22,o32}(x2 + y2 + z2)
a1x + a2y + a3-- --
a1 «2.
- 0 (1 + V )V2
-(b1 + b. + |) y2- b3z2
< ekt(1 + V)0 Ak(1 + V)2 + -0(0 + 1)V4max{o12,o22,o32^-1r
2 1 2 3 V2
+ 0 (1 + V) V3 maxj „2, o22, o32} — - 0 (1 + V) V2
- rnaxj b1 + b2 + j2, b3\ VT2j
min a1, a2,
- c^ll
a U2] V
< ekt (1 + V)0
k + 0 max-j b1 + h-1, b2 + h-2, b3 f1 f2
2k + 0 max{o2, o|, „32} - 0 min ja1, a2, a3 -
+ 0 maxj b1 + h1, b2 + h2, b3 f1 f2
k + 10(0 + 3) max{o12, o2:, „3} - 0 mini a1, a2, a3 - — - — 2 [ a a2
V4(o2x2 + o22y2 + o32z2) < V2ma^jo2,o2,o32},
( C1 C2 \ . i C1
a1x + a2y +1 a3----Iz > mim a1,a2,a3--
\ a a2 / I a1
Hence, we get Lekt(1 + V)0 < K4ekt. Then
CAU a a.2) V'
(3) (44)
dek(1 + V)0 < K4ektdt-0ektV2(1 + V)0^[o^dB^t) + o2ydB2(t) + 03zdB3(t)]. (45)
Integrating both sides of the inequality from 0 to t and then taking the expectations, we can see that
Eekt(1 + V)0 < (1 + V(0))0 + -jekt = (1 + V(0))0 +, (46)
where K5 = K4, so
limsupE[V(t)]° < limsupE[1 + V(t)]0 < K5. (7)
Since (x + y + z)0 < 30 (x2 + y2 + z2) 2 = 30 |x + y + z|0 and V(x,y, z) = x+y+z,
1 < 30--= 30 V0 (t). (48)
|x + y + z|0 (x + y + z)0 Clearly,
limsup E
t—TO
< 30 limsupE[V0(t)] < 30K5 = K6, (49)
t—TO
Jx + y + z|0_
which is the desired assertion. □
Theorem 4 Under Assumption 2, system (4) is stochastically permanent.
Proof By Theorem 2, we know that limsupt—TOE|x + y + z|p < L(p). Let x = (Lep))p, for all e > 0. By Chebyshev's inequality, we can obtain the required assertion. □
5 Global asymptotic stability
Definition 3 (see [3] and references therein) 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) e R+. That is,
P{t——TOE[|(x1(t),y1(t),z1(t)) - (*2(t),y2(t),z2(t))|] =0 = 1. (50)
Lemma 3 [19] Suppose that an n-dimensional stochastic process X(t) ont > 0 satisfies the condition
E|X(t)-X(s)|a < c|t - s|1+^, 0 < s, t < to, (1)
for some positive constants a, f$, and c. There exists a continuous modification X (t) ofX(t) which has the property that for every & e (0, a) there is a positive random variable h(a>) such that
X(t, a) - X(s, to) 2
0<|t-s|<h(ffl),0<s,t<TO |t - s|& 1 - 2-
sup = 1. (52)
In other words, almost every sample path ofX(t) is locally but uniformly Holder continuous with exponent &.
Lemma 4 Let (x(t), y(t), z(t)) be a positive solution of system (4) ont > 0 with initial value (x0,y0, z0) e R+. Then almost every sample path of (x(t),y(t), z(t)) is uniformly continuous ont > 0.
Proof The first equation of system (4) is equivalent to the following stochastic integral equation:
hix(s)
ciz(s)
x(t) = x(0) + i x(s)(ai - bix(s) - r\ar\
J0 \ fi+ giy(s) 1 + a-x(s) +P-z(s)
+ / aix(s) dBi(s).
We estimate
x ai - bix -
fi+ giy i+aix + P-z
i 9„ i
- -E|x|2p + -E
i 2p i
- -E|x|2p + -E
- 2 | | 2
ai - bix -
fi+ giy i+aix + P-z
hix ci ai + bix + — + — Ji Pi
32p-i 2
E|x|2p + la2/ + ( bi + h^) E|x|2p +
hi ^ fi, hi^
- 2Ki(2p) + ^ jaf + bi + ^ Ki(2p) + ^p =: K5 (p).
In addition, applying the moment inequality for stochastic integrals [18], we have, for 0 < t1 < t2 < to and p >2,
a1x(s) dB1(s)
- (--Y
p(p -1) ' 2
p(p -1)
(t2- ti)"^ E|x(s)|pds
'(t2- ti)2Ki(p).
Let t2 -1, — 1 and - + - = 1, we obtain
E|x(t2)-x(ti)|
hix(s)
ciz(s)
- 2p-1E + 2p-1E
x(s)( ai- bix(s)-,. , , . ,\a,\ ti V fi+ giy(s) 1 + a-x(s) +Piz(s)
ds + / a1x(s) dB,(s)
hix(s)
ciz(s)
x(s)( ai- bix(s)-,. , , . ,\a,\ ti V fi+ giy(s) 1 + a-x(s) +Piz(s)
a, x(s) dB,(s)
- 2p-1(t2- ti)q E ti
x(s) I a, - b,x(s) -
hix(s)
ciz(s)
fi + giy(s) 1 + a-x(s) + P-z(s),
+ 2p - >2) P
p(p -1)
(t2 - tl)2Ki(p)
< 2p-1(t2 - h)PK5(p) + 2p-1(t2 - ti)p (a2)p p(p - 1)
p(p - 1)
< 2p- 1(t2- ti)^ |l + < 2p- 1(t2- ti)^1 +
p(p - 1)
max{Ks(p), (a2)PKi(p)} K6(p),
where K6(p) = max{K5(p), (a12)PKL(p)j. Then according to Lemma 3, we know that almost every sample path of x(t) is locally but uniformly Holder continuous with exponent & for every & e (0, P-2). Therefore almost every sample path of x(t) is uniformly continuous on t > 0. In the same way, we can demonstrate that almost every sample path ofy(t), z(t) is uniformly continuous on t > 0. □
Lemma 5 [1] Letf be a non-negative function defined on R+ such thatf is integrable on R+ and is uniformly continuous ont > 0. Then limt^TOf (t) = 0.
Theorem 5 If
A - h h h1 c№ 9 > h2g2 [a2 +
A = b1 - h1 - f. - — - 2d--bf-> 0,
J) - h h h2 C2a2 r¡i h1g1[a1 + n B = b2 - h2 - h2 - "/j2 2d2--b1f2-> 0,
C = b3 - 2c1 - 2c2 -
ftdj ^2d2
then system (4) is globally attractive.
Proof Define W(t) = | lnxi(t) - lnx2(t)| + | lnyi(t) - lny2(t)| + | lnzi(t) - lnz2(i)|. We can prove W(t) is continuous positive function on t > 0. By a direct calculation of the right differential d+W(t) of W(t), and then applying Ito's formula, we get
dx2 (dx2)2n x2 2x22
d+ W (t) = sgn(x1 - X2) + sgn(y1- y2) + sgn(z1 - Z2)
dx1 (dx1)2
x1 2x2 -
'dy1 (dyO2" ' dy 2 (dy2)2]
. y1 2y2 . . y 2 2y2 J
dz1 (dz1)21 ' dz2 (dz2)21
. z1 2z2 . . z2 2z2 J
= sgn(x1 - x2H -b1(x1 - x2) -
h1x1 h1x2
1+ gy f1+ g1y2
1 + a1x1 + jj1z1 1 + a1x2 + j1z2
+ sgn(y1 -y2H -h(y1 -y2) -
2 + g2x1 f2+ g2x2
1 + a2y1 + j2z1 1 + a2y2 + j2z2
+ sgn(z1 - z2 H -b3 (z1 - z2) +
1 + a1x1 + a1z1 1 + a1x2 + a1z2,
, 1 + a2y1 + ^2z1 1 + a2y2 + faz2, Integrating both sides of the equality from 0 to t and then taking expectations
E[W(t) - W(0)]
_sgn(x1(s)-x2(s^{-b^x1(s)-x2(s^ - "fT^)
c1z1(s) c1z2(s)
1 + a1x1 (s) + (s) 1 + a1x2 (s) + ^1z2 (s),
h2y1(s)
■ sgn(y1(s) -y2(s)^-b1(y1(s) -y2(s)) - f
h2y2(s) 2 + g2x1(s) f2+ g2x2(s)
c2z1(s)
c2z2(s)
, 1 + a2y1 (s) + (s) 1 + a2y2 (s) + ^2 (s)
+ sgn (z1 (s) - z2 (s)) j -b3 (z1 (s) - z2 (s))
d1x1(s) d1x2(s)
1 + a1x1(s) + p1z1{s) 1 + a1x2 (s) + 0^2 (s) d2y1(s) d2y2(s)
1 + a2y1 (s) + (s) 1 + a2y2 (s) + fcz2 (s)
dE[W (t)] dt
_sgn(x1(t) -b1(x1(t) -x2(t)) - f1^ -f^y^k)
c1 z1(t)
c1z2(t)
1 + a1x1 (t) + fa (t) 1 + a1x2 (t) + Az2 (t),
h2y1(t)
sgn(y1(t) -y2(t)){-b2(y1(t) -y2(t)) - f
h2y2(t) 2 + g2x1(t) f2+ g2x2(t)
c2z1(t)
c2z2(t)
, 1 + a2y1 (t) + ^2z1 (t) 1 + a2y2 (t) + ^2 (t)
+ sgn(z1(t)-z2(t))| -b3(z1(t)-z2(t))
d1 x1(t) d1x2(t)
1 + a1x1 (t) + Az1(t) 1 + a1x2 (t) + ftz2 (t) d2y1(t) d2y2(t)
1 + a2y1 (t) + (t) 1 + a2y2 (t) + 0^2z2 (t) h1x1(t) h1x2(t)
< -b1E|x1(t) - x2(t) | + E c1z1(t)
f1+ g1y1(t) f1+ g1y2(t) c1z2(t)
1 + a1x1 (t) + 01z1 (t) 1 + a1x2 (t) + 01z2 (t)
h2y1(t) h2y2(t)
- b2E|yl(t)-y2(t)| -E
f2+ g2x1(t) f2+ g2x2(t)
C2 Zi(t)
C2Z2(t)
1 + a2yi (t) + ^2Zl (t) 1 + a2y2 (t) + P2Z2 (t) d1x1(t)
-b3E\zi(t) - Z2(t)| + E d2yi(t)
di X2 (t)
i + aixi (t) + fiiZi (t) i + aiX2 (t) + ^iZ2 (t) d2y2(t)
i + a2yi (t) + ^2Zi (t) i + a2y2 (t) + P2Z2 (t)
< - biE\xi(t) - X2(t)\ + hi(i + j )E\xi(t) - X2(t)\
E\xi(t)\E\yi(t) - y2(t) \ + 2ciE\Zi(t) - Z2(t)\
higi + fi2
+ C^E\xi(t) -x2(t)\ -b2E\yi(t) -y2(t)\ Pi
+ J^j E\yi(t)-y2(t)\ + h2f2E\yi(t)\E\xi(t)-x2(t)\
+ 2c2E\Zi(t)-Z2(t)\ + C2a2E\yi(t) -y2(t)\ P2
- b3E\Zi(t)-Z2(t)\ + 2diE\xi(t) -x2(t)\
+ PidiE\Zi(t)-Z2(t)\ + 2d2E\yi(t) -y2(t)\ + — E\Zi(t)-Z2(t)\, (60)
dE[W (t)] dt
bi - hi- ^ - Cai-2di- h-22 E\yi(t)\l E\xi(t)-x2(t)\
fi Pi f22
b2 - -2 - h2 - CP02 -2d2 - ^E\xi(t)\lE\yi(t) -y2(t)\
f2 P2 fi2
b3 - 2ci - 2C2 -
Pidi P2d2 ai a2
E\Zi(t)-Z2(t)\.
By Lemma 2,
E\xi(t)\ = E\x3(t)\3 = {E[x3(t)]}3 < fli+bygi ,
E\yi(t)\ = E\y3(t)\3 = {E[y3(t)]}3 < ^^p-,
di d2 p-i 2
E\Zi(t)\ = E\Z3(t)\3 = {E[Z3(t)]}3 < 3 ai ¡2 2°\
dE[W (t)] dt
hi ciai ^ , h2g2(a-2 + V^2)
bi - hi -------2di -
, , h2 c2a2 b2 - h2 - —----2d2 -
higi(ai + p^of)1
b3 - 2ci - 2c2 -
Pidi P2 d2
bifi2 E\Zi(t) - Z2(t)\
E\xi(t) - x2(t)\ E\yi(t) - y2(t)\
< -AE\xi(t) - x2(t)\-BE\yi(t) - y2(t)\ - CE\Zi(t) - Z2(t)\,
A b b h1 ca , h2g2 (a2 + Pp-1 a22) . . A = b1- h1--------2d1--—2-, (66)
f1 bf2 h2 c2a2 h1g1(a1 + p-1 a12)
J-~ftT-2d2 f
B = b2- h2 - -r —- 2d2 ———, 2 , (67)
run o a1d1 °2d2
C = b3 - 2c1 - 2c2----. (8)
Integrating both sides, we get E[W(t)] < E[W(0)]
- f [AE^s) -x2(s)| + BE|y1(s) -y2(s)| + CE|z1(s)-z2(s)|]ds. (69)
Consequently,
W(t)+ i [AE|x1(s)-x2(s)|+ BE|y1(s)-y2(s)| + CE|z1(s)-z2(s)|] ds < W(0) < to. (70)
Recall that W(t) > 0 and A > 0, B > 0, C >0, yielding
k(t)-x2(t)| eL1[0, +to), |y1(t)-y2(t)| eL1[0,+to), 1 2 1 2 (71)
|z1(t)-z2(t)| e L [0, +to). Then from Lemmas 4 and 5, we get the desired assertion. □
6 Numerical simulations
In this section, to substantiate the analytical results; the dynamics of system (4) with and without environmental noises are illustrated by the Milstein method [20]. Consider the discrete equations:
/ , h1xk c1zk \ . ^
xk+1 = xk + xk a1 - b1xk - T-----— At
\ f1+ g1yk 1 + a1xk + 01 zk )
+ a1xkHkVAt + 2a12xk(Hk - 1) At,
( U h2yk c2zk \
yk+1 = yk + yA a2 - b2yk - 7-----— At
\ f2 + g2xk 1 + a2xk + 02zk /
+ a2yknkVA + 2a%yk(nl - 1) At,
d1xk d2yk
zk+1 = zk + zk a3 - b3zk + ---— + ---— At
\ 1 + a1xk + fa 1 + a2xk + 02 zk /
+ a3zkZkVAt + 2a^zk (Zk -1) At,
where Hk, nk, and Zk are Gaussian random variables that follow N(0,1).
In Figure 1, we choose the initial value (x0,y0, z0) = (0.5,0.5,0.5) and all the parameters satisfying the conditions of Theorem 3, the system (4) is stochastically permanent.
Figure 1 Solutions of system (4)for(x0,y0,z0) = (0.5, 1.3 0.5,0.5), a1 =2, a2 = 1, a3 = 1.5, b1 =2, b2 = 2, b3 = 3, c1 = 0.2, c2 = 0.1, d1 = 1, d2 = 1, a1 = 1, fa = 0.5, -| 2
a2 = 0.6, fa = 0.8, h1 = 0.5, h2 = 0.8, f1=1, f2 = 0.3, '
g1 =0.5, g2 = 1, a1 =0.05, =0.05, = 0.05.
0 100 200 t
sqrt(x2+y2+z2)
(a) (b)
Figure 2 Solutions of system (4)for(x0,y0,z0) = (0.5,0.5,0.4), a1 =2, a2 = 2, a3 = 1, b1 =2.5, b2 = 2.5, b3 = 2.6, c1 =0.1, c2 = 0.15, d1 = 0.8, d2 = 0.8, a1 = 1, fa = 0.8, a2 = 1, fa = 0.8, h1 = 0.25, h2 = 0.3, f1 = 1, f2 = 0.8, g1 = 0.5, g2 = 0.6, a1 = 0.05, a2 = 0.05, a2 = 0.05.
0.5 Population x
Population y
— Population z
Population x
— Population y
— Population z
Figure 3 Solutions of system (4) for (x0,y0,z0) = (1.5,1.2,1.0), a1 =1.1, a2 = 1.1, a3 = 1,b1 = 0.8, b2 = 0.9, b3 = 1.1, c1 = 0.02, c2 = 0.01, d1 =1.2, d2 = 1, a1 = 0.8, fa = 0.5, a2 = 0.7, fa = 0.5, h1 = 0.5, h2 = 0.3, f1 =1, f2 = 0.5, g1 = 0.5, g2 = 0.6, a1 = 1.8, a2 = 1.8, a2 = 2.
In Figure 2, we choose the initial value (x0,y0,z0) = (0.5,0.5,0.5) and all the parameters satisfying conditions of Theorem 5. Using Matlab, we see that, under a small perturbation, the solution of system (4) is fluctuating in a small neighborhood, at this time, the stochastic system is getting more similar to the deterministic.
In Figure 3, suffering sufficiently large white noise, system (4) gets extinct, while none of the species in the deterministic system die out.
7 Conclusion
In this paper, we studied a stochastic cooperative predator-prey system with Beddington-DeAngelis functional response. We first demonstrate the existence and uniqueness of global positive solution. We then investigate that the solution is stochastically bounded and permanent. Under some specific conditions, the solution of this stochastic system is globally asymptotically stable, which is useful to estimate the risk of extinction of species in the system. Some interesting topics deserve further study. For example, we can try to study three or more trophic levels, which simulates the actual situation better. Moreover, we can consider the extinction of system (4).
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
XL proposed the model and drafted the main part of the manuscript. YW contributed to the revision of the manuscript
and polished the language. All authors read and approved the final manuscript.
Acknowledgements
This work was supported by Shanxi Scholarship Council of China (No. 800104-02090206).
Received: 11 May 2015 Accepted: 16 November 2015 Published online: 26 January 2016
References
1. Qiu, H, Liu, M, Wang, K, Wang, Y: Dynamics of a stochastic predator-prey system with Beddington-DeAngelis functional response. Appl. Math. Comput. 219,2303-2312 (2012)
2. Skalski, GT, Gilliam, JF: Functional responses with predator interference: viable alternatives to the Holling type II model. Ecology 82,3083-3092 (2001)
3. Huang, D, Li, Y, Guo, Y: Dynamics of a stochastic cooperative predator-prey system. Abstr. Appl. Anal. 2014, Article ID 927217(2014)
4. Cosner, C, DeAngelis, DL, Ault, JS, Olson, DB: Effects of spatial grouping on the functional response of predators. Theor. Popul. Biol. 56,65-75 (1999)
5. Cantrell, RS, Cosner, C: On the dynamics of predator-prey models with the Beddington-DeAngelis functional response. J. Math. Anal. Appl. 257,206-222 (2001)
6. Hwang, TW: Uniqueness of limit cycles of the predator-prey system with Beddington-DeAngelis functional response. J. Math. Anal. Appl. 290,113-122 (2004)
7. Dimitrov, D, Kojouharov, H: Complete mathematical analysis of predator-prey models with linear prey growth and Beddington-DeAngelis functional response. Appl. Math. Comput. 162, 523-538 (2005)
8. Cui, J, Takeuchi, Y: Permanence, extinction and periodic solution of predator-prey system with Beddington-DeAngelis functional response. J. Math. Anal. Appl. 317,464-474 (2006)
9. Lin, G, Hong, Y: Delay induced oscillation in predator-prey system with Beddington-DeAngelis functional response. Appl. Math. Comput. 190,1296-1311 (2007)
10. Ko, W, Ryu, K: Analysis of diffusive two-competing-prey and one-predator systems with Beddington-DeAngelis functional response. Nonlinear Anal. 71,4185-4202 (2009)
11. Li, H, Takeuchi, Y: Dynamics of the density dependent predator-prey system with Beddington-DeAngelis functional response. J. Math. Anal. Appl. 374,644-654 (2011)
12. Baek, H: Qualitative analysis of Beddington-DeAngelis type impulsive predator-prey models. Nonlinear Anal., Real World Appl. 11,1312-1322 (2010)
13. Goh, BS: Stability in models of mutualism. Am. Nat. 113, 261-275 (1979)
14. Yang, PH, Xu, R: Global asymptotic stability of periodic solution in n-species cooperative system with time delays. J. Biomath. 13,841-862 (1998)
15. May, RM: Stability and Complexity in Model Ecosystems. Princeton University Press, Princeton (2001)
16. Jost, C, Arditi, R: From pattern to process: identifying predator-prey interactions. Popul. Ecol. 43,229-243 (2001)
17. Li, X, Mao, X: Population dynamical behavior of nonautonomous Lotka-Volterra competitive system with random perturbation. Discrete Contin. Dyn. Syst., Ser. A 24, 523-545 (2009)
18. Mao, X: Stochastic Differential Equations and Applications. Horwood, Chichester (2007)
19. Karatzas, I, Shreve, SE: Brownian Motion and Stochastic Calculus. Springer, Berlin (1991)
20. Higham, DJ: An algorithmic introduction to numericalsimulation of stochastic differentialequations. SIAM Rev. 43, 525-546 (2001)
Submit your manuscript to a SpringerOpen journal and benefit from:
► Convenient online submission
► Rigorous peer review
► Immediate publication on acceptance
► Open access: articles freely available online
► High visibility within the field
► Retaining the copyright to your article
Submit your next manuscript at ► springeropen.com
