# The Asymptotic Behavior of a Stochastic Predator-Prey System with Holling II Functional Response

Authors: Zhenwen Liu; Ningzhong Shi; Daqing Jiang; Chunyan Ji
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

We discuss a stochastic predator-prey system with Holling II functional response. First, we show that this system has a unique positive solution as this is essential in any population dynamics model. Then, we deduce the conditions that there is a stationary distribution of the system, which implies that the system is permanent. At last, we give the conditions for the system that is going to be extinct.


---

## Body

## 1. Introduction

One of the most popular predator-prey model is the one with Michaelis-Menten type (or Holling Type II) functional response [1, 2]:
(1.1)x˙(t)=x(t)(a-bx(t)-αy(t)1+βx(t)),y˙(t)=y(t)(-e+kαx(t)1+βx(t)),
where x(t) and y(t) are the population densities of prey and predator at time t, respectively. The constants a, b/a, α, β, e, and k are positive constants that stand for prey intrinsic growth rate, carrying capacity, the maximum ingestion rate, half-saturation constant, predator death rate, and the conversion factor, respectively. This model exhibits the well-known but highly controversial “paradox of enrichment” observed by Hairston et al. [3] and by Rosenzweig [4] which is rarely reported in nature. It is very important to study the existence and asymptotical stability of equilibria and limit cycle for autonomous predator-prey systems with Holling II functional response. If kaαβ>aeβ2+kbα+beβ, then system (1.1) has a unique limit cycle which is stable. If akα>aeβ+be, then system (1.1) has a unique positive equilibrium:
(1.2)x*=ekα-eβ,y*=kα(akα-aeβ-be)(kα-eβ)2,
which is a stable node or focus (see [5]).However, countless organisms live in seasonally or diurnally forced environments. Hence, authors considered models with periodic ecological parameters or perturbations. For example, Liu and Chen [6] introduced periodic constant impulsive immigration of predator into system (1.1) and gave conditions for the system to be extinct and permanence, respectively. Zhang and Chen [7] studied a Holling II functional response food chain model with impulsive perturbations. Zhang et al. [8] further considered system (1.1) with periodic constant impulsive immigration of predator and periodic variation in the intrinsic growth rate of the prey.On the other hand, the white noise is always present, and we cannot omit the influence of the white noise to the system. May [9] pointed out that due to continuous fluctuation in the environment, the birth rates, death rates, carrying capacity, competition coefficients, and all other parameters involved with the model exhibit random fluctuation to a great lesser extent, and as a result the equilibrium population distribution never attains a steady value, but fluctuates randomly around some average value. Many authors studied the effect of the stochastic perturbation to the predator-prey system with different functional responses, such as [10–14]. Therefore, in this paper, we also introduce stochastic perturbation system (1.1) and obtain the following stochastic system:
(1.3)dx(t)=x(t)(a-bx(t)-αy(t)1+βx(t))dt+σ1x(t)dB1(t),dy(t)=y(t)(-e+kαx(t)1+βx(t))dt+σ2y(t)dB2(t),
where B1(t) and B2(t) are mutually independent Brownian motion with B1(0)=B2(0)=0, and σ12, σ22 are intensities of the white noise.The aim of this paper is to discuss the long time behavior of system (1.3). As the deterministic population models, we are also interested in the permanence and extinction of the system. The global stability of the positive equilibrium means that the system is permanence. But, for the stochastic system, there is no positive equilibrium. Hence, it is impossible that the solution of system (1.3) will tend to a fixed point. In this paper, we show that there is a stationary distribution of system (1.3) mainly according to the theory of Has’meminskii [15], if the white noise is small. While if the white noise is large, based on the techniques developed in [16, 17], we prove that the predator population will die out a.s. and the prey population will either extinct or its distribution converges to a probability measure. It does not happen that both the prey population and the predator population in system (1.3) will die out, which is brought by large white noise, such as weather, epidemic disease. From this point, we say the stochastic model is more realistic than the deterministic model.The rest of this paper is organized as follows. In Section2, we show that there is a unique nonnegative solution of system (1.3). In Section 3, we show that there is a stationary distribution under small white noise. While in Section 4, we consider the situation when the white noise is large. We prove that the system will be extinct. Finally, we give an appendix containing the stationary distribution theory used in Section 3.


## 2. Existence and Uniqueness of the Nonnegative Solution

To investigate the dynamical behavior, the first concern is the global existence of the solutions. Hence in this section we show that the solution of system (1.3) is global and nonnegative. It is not difficult to check the uniqueness and global existence of solutions if the coefficients of the equation satisfy the linear growth condition and local Lipschitz condition (cf. [18]). However, the coefficients of system (1.3) do not satisfy the linear growth condition, but locally Lipschitz continuous, so the solution of system (1.3) may explode at a finite time. In this section, by changing variables, we first show that system (1.3) has a local solution, then show that this solution is global.Theorem 2.1.
For any initial value(x(0),y(0))∈R+2, there is a unique solution (x(t),y(t)) of system (1.3) on t≥0, and the solution will remain in R+2 with probability 1.Proof.
First, consider the following system, by changing variables,x(t)=eu(t), y(t)=ev(t),
(2.1)du(t)=(a-σ122-beu(t)-αev(t)1+βeu(t))dt+σ1dB1(t),dv(t)=(-e-σ222+kαeu(t)1+βeu(t))dt+σ2dB2(t).
It is clear that the coefficients of system (2.1) are locally Lipschitz continuous for the given initial value (logx(0),logy(0))∈R2 there is a unique local solution (u(t),v(t)) on t∈[0,τe), where τe is the explosion time (see [18]). Hence, by Itô formula, we know (eu(t),ev(t)),t∈[0,τe) is a unique positive local solution of system (1.3). To show that this solution is global, we need to show that τe=∞ a.s. Let m0≥1 be sufficiently large so that x(0),y(0) all lie within the interval [1/m0,m0]. For each integer m≥m0, define the stopping time:
(2.2)τm=inf{t∈[0,τe):min{x(t),y(t)}≤1mormax{x(t),y(t)}≥m},
Where, throughout this paper, we set inf∅=∞ (as usual ∅ denotes the empty set). Clearly, τm is increasing as m→∞. Set τ∞=limm→∞τm, whence τ∞≤τe a.s. If we can show that τ∞=∞ a.s., then τe=∞ and (x(t),y(t))∈R+2 a.s. for all t≥0. In other words, to complete the proof all we need to show is that τ∞=∞ a.s. If this statement is false, then there is a pair of constants T>0 and ϵ∈(0,1) such that
(2.3)P{τ∞≤T}>ϵ.
Hence there is an integer m1≥m0 such that
(2.4)P{τm≤T}≥ϵ∀m≥m1.
Define a C2-function V:R+2→R-+ by
(2.5)V(x,y)=(x-c-clogxc)+1k(y-1-logy),
where c is a positive constant to be determined later. The nonnegativity of this function can be seen from u-1-logu≥0,forallu>0. Using Itô's formula, we get
(2.6)dV:=LVdt+σ1(x-c)dB1(t)+σ2k(y-1)dB2(t),
where
(2.7)LV=(x-c)(a-bx-αy1+βx)+cσ122+1k(y-1)(-e+kαy1+βx)+σ222k=-ac+cσ122+ek+σ222k+(a+bc)x-eky-bx2+αcy1+βx≤-ac+cσ122+ek+σ222k+(a+bc)x-bx2-(ek-αc)y.
Choose c=e/αk such that e/k-αc=0, then
(2.8)LV≤-ac+cσ122+ek+σ222k+(a+bc)x-bx2≤K,
where K is a positive constant. Therefore
(2.9)∫0τm∧TdV(x(t),y(t))≤∫0τm∧T‍Kdt+∫0τm∧T‍σ1(x(s)-c)dB1(s)+σ2k(y(s)-1)dB2(s),
which implies that,
(2.10)E[V(x(τm∧T),y(τm∧T))]≤V(x(0),y(0))+E∫0τm∧TKdt≤V(x(0),y(0))+KT.
Set Ωm={τm≤T} for m≥m1, then by (2.4), we know that P(Ωm)≥ϵ. Note that for every ω∈Ωm, there is at least one of x(τm,ω) and y(τm,ω) equals either m or 1/m, then
(2.11)V(x(τm),y(τm))≥(m-c-clogmc)∧(1m-c+clog(cm))∧1k(m-1-logm)∧1k(1m-1+logm).
It then follows from (2.4) and (2.10) that
(2.12)V(x(0),y(0))+KT≥E[1Ωm(ω)V(x(τm),y(τm))]≥ϵ(m-c-clogmc)∧(1m-c+clog(cm))∧1k(m-1-logm)∧1k(1m-1+logm),
where 1Ωm(ω) is the indicator function of Ωm. Letting m→∞ leads to the contradiction that ∞>V(x(0),y(0))+KT=∞. So we must therefore have τ∞=∞ a.s.


## 3. Permanence

There is no equilibrium of system (1.3). Hence we cannot show the permanence of the system by proving the stability of the positive equilibrium as the deterministic system. In this section we show that there is a stationary distribution of system (1.3).Remark 3.1.
Theorem2.1 shows that there exists a unique positive solution (x(t),y(t)) of system (1.3) with any initial value (x(0),y(0))∈R+2. From the proof of Theorem 2.1, we obtain that LV≤K. Define V~=V+K, then LV~≤V~, and it is clear that V~R=inf(x,y)∈R+2∖DkV~(x,y)→∞ as k→∞, where Dk=(1/k,k)×(1/k,k). Hence by Remark 2 of Theorem 4.1 of Has’meminskii, 1980, page 86 in [15], we obtain that the solution (x(t),y(t)) is a homogeneous Markov process in R+2.Theorem 3.2.
Ifaeβ+be<akα<aeβ+bkα/β and σ1>0,σ2>0 such that σ22<kαx*/(1+βx*) and
(3.1)(12+l2x*)x*σ12+(1+βx*2+l2y*k)y*σ22k<min{12(b-αβy*1+βx*)(x*)2,l22(αx*k(1+βx*)-σ22k2)(y*)2},
where (x*,y*) is the positive equilibrium of system (1.1) and l2 is defined as in the proof. Then system (1.3) has a stationary ergodic solution.Proof.
Sinceakα>aeβ+be, then there is a positive equilibrium (x*,y*) of system (1.1), and
(3.2)a=bx*+αy*1+βx*,e=kαx*1+βx*.
Let
(3.3)V1(x,y)=(x-x*-x*logxx*)+l1(y-y*-y*logyy*),
where l1 is a positive constant to be determined later. Let L be the generating operator of system (1.3). Then
(3.4)LV1=(x-x*)(a-bx-αy1+βx)+x*σ122+l1(y-y*)(-e+kαx1+βx)+l1y*σ222=(x-x*)[-b(x-x*)-α1+βx(y-y*)+αβy*(1+βx*)(1+βx)(x-x*)]+x*σ122+l1(y-y*)kα(x-x*)(1+βx*)(1+βx)+l1y*σ222≤-(b-αβy*1+βx*)(x-x*)2-α1+βx(1-l1k1+βx*)(x-x*)(y-y*)+x*σ122+l1y*σ222.
Choose l1=(1+βx*)/k such that 1-l1k/(1+βx*)=0 and yields
(3.5)LV1≤-(b-αβy*1+βx*)(x-x*)2+x*σ122+(1+βx*)y*σ222k.
Let
(3.6)V2(x,y)=12[(x-x*)+1k(y-y*)]2.
Note that
(3.7)d[(x-x*)+1k(y-y*)]=(ax-bx2-eky)dt+σ1xdB1(t)+σ2kydB2(t)=[-bx(x-x*)2+αy*(x-x*)-x*(y-y*)1+βx*]dt+σ1xdB1(t)+σ2kydB2(t),
then
(3.8)LV2=[(x-x*)+1k(y-y*)][-bx(x-x*)2+αy*(x-x*)-x*(y-y*)1+βx*]+σ122x2+σ222k2y2=-bx(x-x*)2+(αy*1+βx*+by*k)(x-x*)2-bk(x-x*)2y-αx*k(1+βx*)(y-y*)2+(αy*k(1+βx*)-αx*1+βx*-bx*k)(x-x*)(y-y*)+σ122x2+σ222k2y2≤(αy*1+βx*+by*k+σ12)(x-x*)2-(αx*k(1+βx*)-σ22k2)(y-y*)2+(αy*k(1+βx*)-αx*1+βx*-bx*k)(x-x*)(y-y*)+σ12(x*)2+σ22k2(y*)2,
where L is also the generating operator of system (1.3). Note that
(3.9)(αy*k(1+βx*)-αx*1+βx*-bx*k)(x-x*)(y-y*)≤(αy*/(k(1+βx*))-αx*/(1+βx*)-bx*/k)22(αx*/(k(1+βx*))-σ22/k2)(x-x*)2+12(αx*k(1+βx*)-σ22k2)(y-y*)2:=δ(x-x*)2+12(αx*k(1+βx*)-σ22k2)(y-y*)2,
then
(3.10)LV2≤(αy*1+βx*+by*k+σ12+δ)(x-x*)2-12(αx*k(1+βx*)-σ22k2)(y-y*)2+σ12(x*)2+σ22k2(y*)2.
Now define
(3.11)V(x,y)=V1(x,y)+l2V2(x,y),
where l2 is a positive constant to be determined later. Then
(3.12)LV≤-(b-αβy*1+βx*-l2(αy*1+βx*+by*k+σ12+δ))(x-x*)2-l22(αx*k(1+βx*)-σ22k2)(y-y*)2+(12+l2x*)x*σ12+(1+βx*2+l2y*k)y*σ22k.
Choose l2>0 such that (b-αβy*/(1+βx*)-l2(αy*/(1+βx*)+by*/k+σ12+δ))=(1/2)(b-αβy*/(1+βx*)), then it follows from (3.12) that
(3.13)LV≤-12(b-αβy*1+βx*)(x-x*)2-l22(αx*k(1+βx*)-σ22k2)(y-y*)2+(12+l2x*)x*σ12+(1+βx*2+l2y*k)y*σ22k.
Note that
(3.14)(12+l2x*)x*σ12+(1+βx*2+l2y*k)y*σ22k<min{12(b-αβy*1+βx*)(x*)2,l22(αx*k(1+βx*)-σ22k2)(y*)2},
then the ellipsoid
(3.15)-12(b-αβy*1+βx*)(x-x*)2-l22(αx*k(1+βx*)-σ22k2)(y-y*)2+(12+l2x*)x*σ12+(1+βx*2+l2y*k)y*σ22k=0
lies entirely in R+2. We can take U to be a neighborhood of the ellipsoid with U-⊆El=R+2, so that for (x,y)∈U∖El, LV≤-C (C is a positive constant), which implies condition (B.2) in Lemma A.1 is satisfied. Hence the solution (x(t),y(t)) is recurrent in the domain U, which together with Lemma A.3 and Remark 3.1 implies that (x(t),y(t)) is recurrent in any bounded domain D⊂R+2. Besides, for all D, there is an M=min{σ12x2,σ22y2,(x,y)∈D-}>0 such that
(3.16)∑i,j=12λijξiξj=σ12x2ξ12+σ22y2ξ22≥M|ξ2|allx∈D¯,ξ∈R2,
which implies that condition (B.1) is also satisfied. Therefore, system (1.3) has a stable a stationary distribution μ(·) and it is ergodic.Note that(3.17)dxp=pxp(a-bx-αy1+βx)dt+pσ1xpdB1(t)+12p(p-1)σ12xpdt≤pxp(a+p2σ12-bx)dt+pσ1xpdB1(t),
then
(3.18)dE[xp]dt≤p(a+pσ122)E[xp]-bE[xp+1]≤p(a+pσ122)E[xp]-b(E[xp])(p+1)/p.
Hence by comparison theorem, we get
(3.19)limsupt→∞E[xp(t)]≤(a+pσ12/2b)p,
by which together with the continuity of E[xp(t)], we have that there exists a positive constant K=K(p) such that
(3.20)E[xp(t)]≤K(p).By Doob's martingale inequality, together with the (3.20), for δ>0, we have
(3.21)P{ω:sup(n-1)δ≤t≤nδx(t)t>δ}≤E[xp(nδ)](nδ)pδ≤K(p)npδp+1,p>1.
In view of the well-known Borel-Cantelli lemma, we see that for almost all ω∈Ω,
(3.22)sup(n-1)δ≤t≤nδx(t)t≤δ
holds for all but finitely many n. Hence there exists an n0(ω), for all ω∈Ω excluding a P-null set, for which (3.22) holds whenever n≥n0. Consequently, letting δ→0, we have, for almost all ω,
(3.23)limt→∞x(t)t=0.By the ergodic property, for any given constantm>0, we have
(3.24)limt→∞1t∫0t‍(xp(s)∧m)ds=∫R+2‍(z1p∧m)μ(dz1,dz2),a.s.
On the other hand, by dominated convergence theorem and (3.20), we get
(3.25)E[limt→∞1t∫0t‍(xp(s)∧m)ds]=limt→∞1t∫0t‍E[xp(s)∧m]ds≤K(p),
which together with (3.24) implies
(3.26)∫R+2‍(z1p∧m)μ(dz1,dz2)≤K(p).
Letting m→∞, we get
(3.27)∫R+2‍z1pμ(dz1,dz2)≤K(p).
That is to say, the function f1(z)=z1p is integrable with respect to the measure μ. Therefore, by ergodicity property again, we get
(3.28)limt→∞1t∫0t‍xp(s)ds=∫R+2‍z1pμ(dz1,dz2),a.s.
Besides,
(3.29)x(t)-x(0)t=at∫0t‍x(s)ds-bt∫0t‍x2(s)ds-αt∫0t‍x(s)y(s)1+βx(s)ds+σ1t∫0t‍x(s)dB1(s).
Let M1(t)=∫0t‍x(s)dB1(s) which is a martingale with M1(0)=0 and
(3.30)limsupt→∞〈M1,M1〉tt=limt→∞1t∫0t‍x12(s)ds=∫R+2‍z12μ(dz1,dz2)<∞,
then by the strong law of large numbers, we get
(3.31)limt→∞M1(t)t=limt→∞1t∫0t‍x1(s)dB1(s)=0,
which together with (3.23) and (3.28) implies (3.29) that
(3.32)limt→∞1t∫0t‍x(s)y(s)1+βx(s)ds=aα∫R+2‍z1μ(dz1,dz2)-bα∫R+2‍z12μ(dz1,dz2).Hence from these arguments, we get the following result.Theorem 3.3.
Assume the same conditions as in Theorem3.2. Then one has
(3.33)limt→∞1t∫0t‍xp(s)ds=∫R+2‍z1pμ(dz1,dz2),a.s.,limt→∞1t∫0t‍x(s)y(s)1+βx(s)ds=aα∫R+2‍z1μ(dz1,dz2)-bα∫R+2‍z12μ(dz1,dz2).


## 4. Extinction

In this section, we show the situation when system (1.3) will be extinct.Case 1.
(
a
<
σ
1
2
/
2
).
Obviously,(4.1)dx≤x(a-bx)dt+σ1xdB1(t),
then when a<σ12/2,
(4.2)limt→∞x(t)=0,a.s.
That is to say, for all 0<ϵ1<e+σ22/2, there exist T1=T1(ω) and a set Ωϵ1 such that P(Ωϵ1)>1-ϵ1 and kαx(t)≤ϵ1 for t≥T1 and ω∈Ωϵ1. Then
(4.3)-ey(t)dt+σ2y(t)dB2(t)≤dy(t)≤y(t)(-e+ϵ1)dt+σ2y(t)dB2(t),
and so
(4.4)limt→∞y(t)=0,a.s.Case 2.
(
a
>
σ
1
2
/
2
,
e
+
σ
2
2
/
2
>
k
α
/
β
).
Note that(4.5)dy(t)≤y(t)(-e+kαβ)dt+σ2y(t)dB2(t),
then if e+σ22/2>kα/β, we have
(4.6)limt→∞y(t)=0,a.s..
In this situation, for all 0<ϵ2<a-σ12/2, there exist T2=T2(ω) and a set Ωϵ2 such that P(Ωϵ2)>1-ϵ2 and αy(t)≤ϵ2 for t≥T2 and ω∈Ωϵ2. Then
(4.7)x(t)(a-bx(t)-ϵ2)dt+σ1x(t)dB1(t)≤dx(t)≤x(t)(a-bx(t))dt+σ1x(t)dB1(t),(4.8)(a-σ122-bx(t)-ϵ2)dt+σ1dB1(t)≤dx(t)≤(a-σ122-bx(t))dt+σ1dB1(t).
Consider the following equation:
(4.9)dΦ(t)=(a-σ122-beΦ(t))dt+σ1dB1(t).
If a>σ12/2, (4.9) has the density g*(ζ) such that
(4.10)12σ12g*′(ζ)=(a1-σ122-b11eζ)g*(ζ).
Therefore from (4.8) and the arbitrary of ϵ2, we get that the distribution of logx(t) converges weakly to the probability measure with density g*. Thus, from (4.10), we obtain that the distribution of x(t) converges weakly to the probability measure with density f*(ζ)=C0ζ2(a-σ12/2)/σ12-1e-2bζ/σ12, where C0=(2b/σ12)2(a-σ12/2)/σ12/Γ(2(a-σ12/2)/σ12). Besides, from the ergodic theorem and (4.10), it follows that
(4.11)limt→∞1t∫0t‍x(s)ds=∫-∞∞‍eζg*(ζ)dζ=∫-∞∞‍a-σ12/2bg*(ζ)dζ=a-σ12/2ba.s.
Therefore, by the above arguments, we obtain the following.Theorem 4.1.
Let(x(t),y(t)) be the solution of system (1.3) with the initial value (x(0),y(0))∈R+2. Then,(i)
ifa<σ12/2, then
(4.12)limt→∞x(t)=0,limt→∞y(t)=0,a.s.,(ii)
ifa>σ12/2,e+σ22/2>kα/β, then the distribution of x(t) converges weakly to the probability measure with density f*(ζ)=C0ζ2(a-σ12/2)/σ12-1e-2bζ/σ12, where C0=(2b/σ12)2(a-σ12/2)/σ12/Γ(2(a-σ12/2)/σ12), and
(4.13)limt→∞1t∫0t‍x(s)ds=a-σ12/2b,limt→∞y(t)=0,a.s.



---

*Source: 801812-2012-12-31.xml*