# Dynamical Analysis of a Stochastic Predator-Prey Model with an Allee Effect

Authors: Feng Rao
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

We present and analyze a modified Holling type-II
predator-prey model that includes some important factors such as
Allee effect, density-dependence, and environmental noise. By
constructing suitable Lyapunov functions and applying Itô
formula, some qualitative properties are given, such as the
existence of global positive solutions, stochastic boundedness, and
the global asymptotic stability. A series of numerical simulations
to illustrate these mathematical findings are presented.


---

## Body

## 1. Introduction

The study of the interrelationship between species and their environment, ecology, in such areas as predator-prey and competition interactions, renewable resource management, evolution of pesticide resistant strains, ecological and genetically engineered control of pests, multispecies societies, and so on, is now an enormous field [1]. Predator-prey interaction is one of the basic interspecies relations for ecological and social models, and it is also a basic block of more complicated food chain, food web and biochemical network structure [2]. In recent years, one of the important predator-prey models is the Holling type-II model, originally due to Holling [3, 4], which has been extensively studied in many articles [5–11].In ecological populations, species growth is often damped by an Allee effect. The so-called Allee effect is described by the positive relationship between population density and the per capita growth rate [12, 13]. The logistic model assumes that per capita growth rate declines monotonically with density. But for population subject to an Allee effect, per capita growth rate shows a humped curve increasing (from negative to positive) at low density, up to a maximum at intermediate density and then declining [14]. The Allee effect may arise from a number of sources such as difficulties in finding mates when population density is low, social dysfunction at small population sizes, food exploitation, inbreeding depression, and predator avoidance or defense [15]. This effect usually saturates or disappears as population get larger [12]. The Allee effect has been modeled in different ways [16–18], such as multiplicative Allee effect and additive Allee effect. In an ecological point of view, Allee effect has been modeled into strong and weak cases. The strong Allee effect introduces a population threshold, and the population must surpass this threshold to grow. In contrast, population with a weak Allee effect do not have a threshold [19]. the Allee effect’s strong potential impact on population dynamics has been attracting much attention recently. Detailed investigations relating to the Allee effect may be found in [12–15, 17–28].In reality, there are many benefits to be gained in using stochastic models because real life is full of randomness and stochasticity. A large number of stochastic differential equation models are considered as they provide some additional degree of realism compared to their deterministic counterpart [2]. The environment factors are time dependent randomly varying and should be taken as stochastic [29, 30]. The basic mechanism and factors of population growth like the resources and vital rates-birth, death, immigration and emigration change nondeterministically due to random environment [31]. Most natural phenomena do not follow strictly deterministic laws but rather oscillate randomly about some average values, so that the population density never attains a fixed value with the advancement of time [32]. In predator-prey model, the random fluctuations are also undeniably arising from either environmental variability or internal species. For example, the birth and death processes of individuals are intrinsically stochastic fluctuations which become especially pronounced when the number of individuals is small [33]. Recent advances in stochastic differential equations enable a lot of authors to introduce stochastic factor into the model of physical phenomena, whether it is a random noise in the system of differential equations or environmental fluctuations in parameters [34–45].Based upon these factors, to the best of our knowledge, a small amount of work has been done with stochastic perturbation on a modified Holling type-II predator-prey model which has an Allee effect on prey population and density-dependence of predator population. The purpose of this paper is to study that the stochastic factor has a significant effect on the dynamics of Holling type-II predator-prey model with an Allee effect and density-dependence. The organization of this paper is as follows. In the next section, we present the formulation of mathematical model with Allee effect, density-dependence, and environmental noise. Moreover, we present a qualitative analysis of the stochastic model. In Section3, we use numerical simulations to reveal the influence of noise on the dynamical behavior of the model. A brief discussion is given in Section 4.


## 2. The Model and Analysis

In this paper, we present a Holling type-II predator-prey model that includes some important factors such as an Allee effect on prey population, density-dependence of predator population, and environmental noise. Assuming that the random factors in the environment would display themselves mainly as fluctuations on the intrinsic growth rates of prey and predator, then the stochastic modified Holling type-II predation model takes the following form:(1)dN=N(1-N-mb+N-wPa+N)(rdt+αdη1(t)),dP=P(Na+N-d-cP)(sdt+βdη2(t)).N(0)=N0>0,P(0)=P0>0.Here, we replace the intrinsic growth ratesr (>0) and s (>0) of prey and predator by
(2)r⟶r+αη˙1(t),s⟶s+βη˙2(t).In the above model,N=N(t) and P=P(t) represent prey and predator densities at time t, respectively. mN/(b+N) is the term of additive Allee effect and m∈(0,1) and b∈(0,1) are Allee effect constants. If m<b, the function of prey growth rate with the Allee effect of the prey is called weak Allee effect; if m>b, the modified function of prey growth rate is strong Allee effect. cP stands for the density-dependence of predator population and c>0. w (>0) is the predation rate, a (>0) the half-saturation density of Holling type-II functional response, and d (>0) the predator death rate. The positive constants α and β, respectively, are the coefficients of the effect of the environmental stochastic perturbation on the prey and on the predator population. η˙i(t)(i=1,2) is a standard white noise with increment Δηi(t)=ηi(t+Δt)-ηi(t) that is independent Gaussian random variables N(0,Δt). That is to say, ηi(t)  (i=1,2) is a Brownian motion defined on a complete probability space (Ω,ℱ,P) with a filtration {ℱt}t∈ℝ+ satisfying the usual conditions; namely, it is right continuous and increasing while ℱ0 contains all P-null sets [46].Consider model (1) as an Itô stochastic differential system of the type
(3)dU(t)=F(t,U(t))dt+G(t,U(t))dη(t),U(t0)=U0,
where
(4)U(t)=(N(t),P(t))T,η(t)=(η1(t),η2(t))T,F=(rN(1-N-mb+N-wPa+N)sP(Na+N-d-cP)),G=(αN(1-N-mb+N-wPa+N)00βP(Na+N-d-cP)),
and the solution U(t)  (t>0) is an Itô process, F is a slowly varying continuous component or drift coefficient, and G is the rapidly varying continuous random component or diffusion coefficient [46].


### 2.1. Existence of Global Positive Solutions

For a stochastic differential equation, in order to have a unique global solution (i.e., no explosion in a finite time) for any given initial value, the coefficients of the equation are generally required to satisfy the linear growth condition and local Lipschitz condition [46]. To show that model (1) has a global positive solution, let us firstly prove that the model has a positive local solution by making the change of variables.Theorem 1.
There is a unique local solution(N(t),P(t)) for t∈[0,τe) to model (1) a.s. for any given initial value (N0,P0)∈ℝ+2={N,P∈ℝ2∣N>0,P>0}, where τe is the explosion time.Proof.
Setp(t)=lnN(t), q(t)=lnP(t) and we consider the equations
(2.1)dp=[r(1-ep-mb+ep-weqa+ep)-α22(1-ep-mb+ep-weqa+ep)2]dt+α(1-ep-mb+ep-weqa+ep)dη1(t),dq=[s(epa+ep-d-ceq)-β22(epa+ep-d-ceq)2]dt+β(epa+ep-d-ceq)dη2(t),p(0)=lnN0,q(0)=lnP0
on t≥0. Then we know that the coefficients of (2.1) satisfy the local Lipschitz condition, and there is a unique local solution (p(t),q(t)) on [0,τe). By Itô’s formula, we can see that N(t)=ep(t), P(t)=eq(t) is the unique positive local solution to model (1) with the initial value N0>0, P0>0.Next, we show that the solution of model (1) is global; that is, τe=∞ a.s.Theorem 2.
For model (1) and any given initial value (N0,P0)∈ℝ+2, there is a unique solution (N(t),P(t)) on t≥0 and the solution will remain in ℝ+2 with probability one.Proof.
Letn0>0 be sufficiently large for N0 and P0 lying within the interval [1/n0,n0]. For each integer n≥n0, let us define a sequence of stopping times by
(6)τn=inf{t∈[0,τe]:N(t)∉(1n,n)orP(t)∉(1n,n)},
where we set inf∅=∞ (∅ represents the empty set) in this paper. Since τn is nondecreasing as n→∞, there exists the limit
(7)τ∞=limn→∞τn.
Then
(8)τ∞≤τea.s.
Now, we need to show τ∞=∞ a.s. If this statement is violated, then there exist T>0 and ɛ∈(0,1) such that
(9)P{τ∞≤T}>ɛ.
Thus, there is an integer n1≥n0 such that
(10)P{τn≤T}≥ɛ∀n≥n1.
Define aℂ2-function V:ℝ+2→ℝ+ by
(11)V(N,P)=N+P-lnNP-2,
which is a nonnegative function. If (N(t),P(t))∈ℝ+2, by using Itô’s formula, we compute
(12)dV(N,P)=[β22(Na+N-d-cP)22r(N-1)(1-N-mb+N-wPa+N)+2s(P-1)(Na+N-d-cP)+α22(1-N-mb+N-wPa+N)2+β22(Na+N-d-cP)2]dt+α(N-1)(1-N-mb+N-wPa+N)dη1+β(P-1)(Na+N-d-cP)dη2≤[α22(1-N-mb+N-wPa+N)22r(1-N)(mb+N+wPa+N)+α22(1-N-mb+N-wPa+N)2]dt+[β22(Na+N-d-cP)22s(NPa+N+d+cP)+β22(Na+N-d-cP)2]dt+α(N-1)(1-N-mb+N-wPa+N)dη1+β(P-1)(Na+N-d-cP)dη2≤M1+M2+α(N-1)×(1-N-mb+N-wPa+N)dη1+β(P-1)(Na+N-d-cP)dη2,
where M1, M2 are positive numbers. Integrating both sides of the above inequality from 0 to τn∧T and taking the expectations give rise to
(13)EV(N(τn∧T),P(τn∧T))≤V(N0,P0)+(M1+M2)T.
Set Ωn={τn≤T} and by inequality (10), we get P(Ωn)≥ɛ. For every μ∈Ωn, N(τn,μ) and P(τn,μ) equal either n or 1/n, hence
(14)V(N(τn,μ),P(τn,μ))≥min{n-lnn-1,1n+lnn-1}.
It then follows from (13) that
(15)V(N0,P0)+(M1+M2)T≥E[1Ωn(μ)V(N(τn),P(τn))]≥ɛmin{n-lnn-1,1n+lnn-1},
where 1Ωn is the indicator function of Ωn. Letting n→∞ leads to the contradiction ∞>V(N0,P0)+(M1+M2)T=∞. This completes the proof.


### 2.2.  Stochastic Boundedness

Theorem 3.
For anyθ∈(0,1), there is a positive constant C=C(θ), which is independent of initial value (N0,P0)∈ℝ+2, such that the solution of model (1) has the property that
(16)limsupt→∞E|(N(t),P(t))|θ≤C.Proof.
Define a function by(17)V(N,P)=Nθ+Pθ,(N,P)∈ℝ+2.
For any t>0, using Itô formula gives that
(18)LV(N,P)=rθNθ(1-N-mb+N-wPa+N)+sθPθ(Na+N-d-cP)+α22θ(θ-1)Nθ(1-N-mb+N-wPa+N)2+β22θ(θ-1)Pθ(Na+N-d-cP)2,
as 0<θ<1; then,
(19)LV(N,P)≤rθNθ(1-N)-α22θ(1-θ)Nθ×(1-N-mb+N-wPa+N)2+sθPθNa+N-β22θ(1-θ)Pθ(Na+N-d-cP)2≤rθNθ-α22θ(1-θ)Nθ×(m2(b+N)2+(1-N)(1-N-2mb+N)-1)+sθPθ-β22θ(1-θ)Pθ(c2P2-2c(1-d)P)=rθNθ-α22θ(1-θ)Nθ×(m2(b+N)2+(1-N)×(1-N-2mb+N)-1m2(b+N)2)+Nθ+sθPθ-β22θ(1-θ)Pθ×(c2P2-2c(1-d)P)+Pθ-V(N,P)≤C0-V(N,P),
where C0 is an integer. Hence,
(20)dV(N,P)=LV(N,P)dt+αθNθ(1-N-mb+N-wPa+N)dη1+βθPθ(Na+N-d-cP)dη2≤(C0-V(N,P))dt+αθNθ(1-N-mb+N-wPa+N)dη1+βθPθ(Na+N-d-cP)dη2.
Applying Itô formula yields(21)d(etV(N,P))=et(V(N,P)dt+dV(N,P))≤etC0dt+etαθNθ×(1-N-mb+N-wPa+N)dη1+etβθPθ(Na+N-d-cP)dη2.
Using the sequence of stopping times {τn}n=1∞ defined in the proof of Theorem 2 and from (21), we have
(22)E[et∧τnV(N(t∧τn),P(t∧τn))]≤V(N0,P0)+C0(Eet∧τn-1)∀t≥0.
Letting n→∞ in the latter inequality with a fact that V(N(t∧τn),P(t∧τn))>0 and 0<et∧τn≤et a.s. and using Fatou’s lemma, we obtain
(23)etEV(N,P)≤V(N0,P0)+C0(et-1).
Therefore,
(24)limsupt→∞EV(N(t),P(t))≤C.
On the other hand, since(25)|N,P|θ≤2θmax{Nθ,Pθ}≤2θV(N,P),
we obtain
(26)limsupt→∞E|(N(t),P(t))|θ≤2θlimsupt→∞EV(N(t),P(t))≤2θC≜C(θ).
The proof is complete.Theorem 4.
The solutions of model (1) are stochastically ultimately bounded for any initial value (N0,P0)∈ℝ+2.Proof.
The solutions of model (1) will remain in ℝ+2 for all t≥0 with probability one from Theorem 2. And by Theorem 3, let θ=1/2; there is a C1>0 such that
(27)limsupt→∞E|(N(t),P(t))|1/2≤C1.
For any ɛ>0, set X=C12/ɛ2. Then by Chebyshev inequality,
(28)P{|(N(t),P(t))|>X}≤E|(N(t),P(t))|1/2X.
That is,
(29)limsupt→∞P{|(N(t),P(t))|>X}≤ɛ.
This completes the proof.Next, we study the asymptotic properties of the moment solutions of model (1).Theorem 5.
For any givenθ∈(0,1), there is a H=H(θ)>0; let (N(t),P(t)) be the solution of model (1) with any initial value (N0,P0); then,
(30)limsupt→∞1t∫0t‍E(Nθ+2(τ)+Pθ+2(τ))dτ≤H.Proof.
SetV(N,P):ℝ+2→ℝ+, and from Theorem 3, we have
(31)LV(N,P)≤rθNθ-α22θ(1-θ)Nθ×(N2+m2(b+N)2-2m(1-N)b+N-2N)+sθPθ-β22θ(1-θ)Pθ(c2P2-2c(1-d)P),
so,
(32)LV(N,P)+α2θ(1-θ)4Nθ+2+c2β2θ(1-θ)4Pθ+2≤rθNθ-α2θ(1-θ)4Nθ+2+α2θ(1-θ)Nθ(m(1-N)b+N+N-m22(b+N)2)+sθPθ-c2β2θ(1-θ)4Pθ+2+c(1-d)β2θ(1-θ)Pθ+1≤K,
where K is a positive number. Then
(33)LV(N,P)+k2γ2θ(1-θ)4(Nθ+2+Pθ+2)≤LV(N,P)+α2θ(1-θ)4Nθ+2+c2β2θ(1-θ)4Pθ+2≤K,
where k2γ2=min{α2,c2β2}; we can obtain
(34)EV(N(t),P(t))+k2γ2θ(1-θ)4×∫0t‍E(Nθ+2(τ)+Pθ+2(τ))dτ≤V(N0,P0)+Kt.
That is,
(35)∫0t‍E(Nθ+2(τ)+Pθ+2(τ))dτ≤4(V(N0,P0)+Kt)k2γ2θ(1-θ).
Set H=4K/(k2γ2θ(1-θ)); therefore,
(36)limsupt→∞1t∫0t‍E(Nθ+2(τ)+Pθ+2(τ))dτ≤H.


### 2.3. Stochastic Asymptotic Stability

From Theorem2, model (1) has a global positive solution. By constructing some Lyapunov functions, we analyze stochastic asymptotic stability of the positive equilibria of the stochastic model. The equilibria of the stochastic model (1) we consider are solutions of the following model:
(37)N(1-N-mb+N-wPa+N)=0,P(Na+N-d-cP)=0.
From (37), the positive equilibrium (N,P) satisfies the equation
(38)cN4+(2ac-c+cb)N3+(a2c-c(b-m)-2ac(1-b)+w(1-d))N2+(a2c(b-1)-2ac(b-m)+wb(1-d)-wad)N-a2c(b-m)-wabd=0,N>ad1-d>0.
As an example, we only give the proof of a positive equilibrium (N*,P*) of model (1).Theorem 6.
Let(39)θ=w(a+N*)(r-α2N*)as+β2acP*,
and
(40)A=r(mb(b+N*)+wP*a(a+N*)-1)+α2N*2(mb(b+N*)+wP*a(a+N*)-1)2+β2a2wP*(r-α2N*)2a3(a+N*)(s+β2cP*),B=α2w2N*2a2+cw(a+N*)(r-α2N*)(β2cP*-2s)2a(s+β2cP*),C=α2mwN*ab(b+N*)+α2w2N*P*a2(a+N*);
if A<0, C2-4AB<0, then the equilibrium point (N*,P*) of model (1) is stochastically asymptotically stable in ℝ+2; that is, for any initial value (N0,P0), the solution of model (1) has the property that
(41)limt→∞N(t)=N*,limt→∞P(t)=P*,
almost surelyProof.
Based on the theory of stability of stochastic differential equations [46], we only need to construct a Lyapunov function V(U) satisfying LV≤0 and the identity holds if and only if U=U*, where U=U(t) is the solution of the stochastic differential equation (3) and U* is the equilibrium point of model (3).
Define Lyapunov functions(42)V1(N)=N-N*-N*lnNN*,V2(P)=P-P*-P*lnPP*,
The nonnegativity of this function can be observed from x-1-lnx≥0 on x>0. If (N(t),P(t))∈ℝ+2, applying Itô’s formula, we have
(43)LV1=r(N-N*)(1-N-mb+N-wPa+N)+α2N*2(1-N-mb+N-wPa+N)2=[+wP*(a+N*)(a+N)-1)2r(m(b+N*)(b+N)+wP*(a+N*)(a+N)-1)+α2N*2(m(b+N*)(b+N)+wP*(a+N*)(a+N)-1)2]×(N-N*)2+α2w2N*2(a+N)2(P-P*)2-[rwa+N+α2wN*a+N×(m(b+N*)(b+N)+wP*(a+N*)(a+N)-1)rwa+N+α2wN*a+N]×(N-N*)(P-P*)≤[+α2N*2(mb(b+N*)+wP*a(a+N*)-1)2r(mb(b+N*)+wP*a(a+N*)-1)+α2N*2(mb(b+N*)+wP*a(a+N*)-1)2]×(N-N*)2+α2w2N*2a2(P-P*)2+α2wN*-rwa+N(N-N*)(P-P*)+(α2mwN*ab(b+N*)+α2w2N*P*a2(a+N*))|N-N*||P-P*|.
Similarly, we obtain
(44)LV2=s(P-P*)(Na+N-d-cP)+β2P*2(Na+N-d-cP)2=s(P-P*)(a(a+N)(a+N*)(N-N*)+c(P*-P)a(a+N)(a+N*))+β2P*2(a2(a+N)2(a+N*)2(N-N*)2+c2(P-P*)2+2ac(a+N)(a+N*)×(N-N*)(P-P*)a2(a+N)2(a+N*)2)=β2a2P*2(a+N)2(a+N*)2(N-N*)2+(β2c2P*2-sc)(P-P*)2+as+acβ2P*(a+N)(a+N*)(N-N*)(P-P*)≤β2a2P*2a2(a+N*)2(N-N*)2+(β2c2P*2-sc)(P-P*)2+as+β2acP*(a+N)(a+N*)(N-N*)(P-P*).
Define(45)V(N,P)=V1(N)+θV2(P),
and choose
(46)θ=w(a+N*)(r-α2N*)as+β2acP*.
We compute(47)LV(N,P)=LV1(N)+LθV2(P)≤[β2a2wP*(r-α2N*)2a2(a+N*)(as+β2acP*)r(mb(b+N*)+wP*a(a+N*)-1)+α2N*2(mb(b+N*)+wP*a(a+N*)-1)2+β2a2wP*(r-α2N*)2a2(a+N*)(as+β2acP*)](N-N*)2+[α2w2N*2a2+cw(a+N*)(r-α2N*)(β2cP*-2s)2(as+β2acP*)]×(P-P*)2+(α2mwN*ab(b+N*)+α2w2N*P*a2(a+N*))|N-N*|×|P-P*|=A(N-N*)2+B(P-P*)2+C|N-N*||P-P*|.
Let |U-U*|=(|N-N*|,|P-P*|)T; then, we have
(48)LV(N,P)⩽12|U-U*|T(2ACC2B)|U-U*|.
Clearly, if A<0 and C2-4AB<0 hold, then the above inequality implies LV(N,P)<0 along all trajectories in the first quadrant except (N*,P*). Then the desired assertion (41) follows immediately.


## 2.1. Existence of Global Positive Solutions

For a stochastic differential equation, in order to have a unique global solution (i.e., no explosion in a finite time) for any given initial value, the coefficients of the equation are generally required to satisfy the linear growth condition and local Lipschitz condition [46]. To show that model (1) has a global positive solution, let us firstly prove that the model has a positive local solution by making the change of variables.Theorem 1.
There is a unique local solution(N(t),P(t)) for t∈[0,τe) to model (1) a.s. for any given initial value (N0,P0)∈ℝ+2={N,P∈ℝ2∣N>0,P>0}, where τe is the explosion time.Proof.
Setp(t)=lnN(t), q(t)=lnP(t) and we consider the equations
(2.1)dp=[r(1-ep-mb+ep-weqa+ep)-α22(1-ep-mb+ep-weqa+ep)2]dt+α(1-ep-mb+ep-weqa+ep)dη1(t),dq=[s(epa+ep-d-ceq)-β22(epa+ep-d-ceq)2]dt+β(epa+ep-d-ceq)dη2(t),p(0)=lnN0,q(0)=lnP0
on t≥0. Then we know that the coefficients of (2.1) satisfy the local Lipschitz condition, and there is a unique local solution (p(t),q(t)) on [0,τe). By Itô’s formula, we can see that N(t)=ep(t), P(t)=eq(t) is the unique positive local solution to model (1) with the initial value N0>0, P0>0.Next, we show that the solution of model (1) is global; that is, τe=∞ a.s.Theorem 2.
For model (1) and any given initial value (N0,P0)∈ℝ+2, there is a unique solution (N(t),P(t)) on t≥0 and the solution will remain in ℝ+2 with probability one.Proof.
Letn0>0 be sufficiently large for N0 and P0 lying within the interval [1/n0,n0]. For each integer n≥n0, let us define a sequence of stopping times by
(6)τn=inf{t∈[0,τe]:N(t)∉(1n,n)orP(t)∉(1n,n)},
where we set inf∅=∞ (∅ represents the empty set) in this paper. Since τn is nondecreasing as n→∞, there exists the limit
(7)τ∞=limn→∞τn.
Then
(8)τ∞≤τea.s.
Now, we need to show τ∞=∞ a.s. If this statement is violated, then there exist T>0 and ɛ∈(0,1) such that
(9)P{τ∞≤T}>ɛ.
Thus, there is an integer n1≥n0 such that
(10)P{τn≤T}≥ɛ∀n≥n1.
Define aℂ2-function V:ℝ+2→ℝ+ by
(11)V(N,P)=N+P-lnNP-2,
which is a nonnegative function. If (N(t),P(t))∈ℝ+2, by using Itô’s formula, we compute
(12)dV(N,P)=[β22(Na+N-d-cP)22r(N-1)(1-N-mb+N-wPa+N)+2s(P-1)(Na+N-d-cP)+α22(1-N-mb+N-wPa+N)2+β22(Na+N-d-cP)2]dt+α(N-1)(1-N-mb+N-wPa+N)dη1+β(P-1)(Na+N-d-cP)dη2≤[α22(1-N-mb+N-wPa+N)22r(1-N)(mb+N+wPa+N)+α22(1-N-mb+N-wPa+N)2]dt+[β22(Na+N-d-cP)22s(NPa+N+d+cP)+β22(Na+N-d-cP)2]dt+α(N-1)(1-N-mb+N-wPa+N)dη1+β(P-1)(Na+N-d-cP)dη2≤M1+M2+α(N-1)×(1-N-mb+N-wPa+N)dη1+β(P-1)(Na+N-d-cP)dη2,
where M1, M2 are positive numbers. Integrating both sides of the above inequality from 0 to τn∧T and taking the expectations give rise to
(13)EV(N(τn∧T),P(τn∧T))≤V(N0,P0)+(M1+M2)T.
Set Ωn={τn≤T} and by inequality (10), we get P(Ωn)≥ɛ. For every μ∈Ωn, N(τn,μ) and P(τn,μ) equal either n or 1/n, hence
(14)V(N(τn,μ),P(τn,μ))≥min{n-lnn-1,1n+lnn-1}.
It then follows from (13) that
(15)V(N0,P0)+(M1+M2)T≥E[1Ωn(μ)V(N(τn),P(τn))]≥ɛmin{n-lnn-1,1n+lnn-1},
where 1Ωn is the indicator function of Ωn. Letting n→∞ leads to the contradiction ∞>V(N0,P0)+(M1+M2)T=∞. This completes the proof.


## 2.2.  Stochastic Boundedness

Theorem 3.
For anyθ∈(0,1), there is a positive constant C=C(θ), which is independent of initial value (N0,P0)∈ℝ+2, such that the solution of model (1) has the property that
(16)limsupt→∞E|(N(t),P(t))|θ≤C.Proof.
Define a function by(17)V(N,P)=Nθ+Pθ,(N,P)∈ℝ+2.
For any t>0, using Itô formula gives that
(18)LV(N,P)=rθNθ(1-N-mb+N-wPa+N)+sθPθ(Na+N-d-cP)+α22θ(θ-1)Nθ(1-N-mb+N-wPa+N)2+β22θ(θ-1)Pθ(Na+N-d-cP)2,
as 0<θ<1; then,
(19)LV(N,P)≤rθNθ(1-N)-α22θ(1-θ)Nθ×(1-N-mb+N-wPa+N)2+sθPθNa+N-β22θ(1-θ)Pθ(Na+N-d-cP)2≤rθNθ-α22θ(1-θ)Nθ×(m2(b+N)2+(1-N)(1-N-2mb+N)-1)+sθPθ-β22θ(1-θ)Pθ(c2P2-2c(1-d)P)=rθNθ-α22θ(1-θ)Nθ×(m2(b+N)2+(1-N)×(1-N-2mb+N)-1m2(b+N)2)+Nθ+sθPθ-β22θ(1-θ)Pθ×(c2P2-2c(1-d)P)+Pθ-V(N,P)≤C0-V(N,P),
where C0 is an integer. Hence,
(20)dV(N,P)=LV(N,P)dt+αθNθ(1-N-mb+N-wPa+N)dη1+βθPθ(Na+N-d-cP)dη2≤(C0-V(N,P))dt+αθNθ(1-N-mb+N-wPa+N)dη1+βθPθ(Na+N-d-cP)dη2.
Applying Itô formula yields(21)d(etV(N,P))=et(V(N,P)dt+dV(N,P))≤etC0dt+etαθNθ×(1-N-mb+N-wPa+N)dη1+etβθPθ(Na+N-d-cP)dη2.
Using the sequence of stopping times {τn}n=1∞ defined in the proof of Theorem 2 and from (21), we have
(22)E[et∧τnV(N(t∧τn),P(t∧τn))]≤V(N0,P0)+C0(Eet∧τn-1)∀t≥0.
Letting n→∞ in the latter inequality with a fact that V(N(t∧τn),P(t∧τn))>0 and 0<et∧τn≤et a.s. and using Fatou’s lemma, we obtain
(23)etEV(N,P)≤V(N0,P0)+C0(et-1).
Therefore,
(24)limsupt→∞EV(N(t),P(t))≤C.
On the other hand, since(25)|N,P|θ≤2θmax{Nθ,Pθ}≤2θV(N,P),
we obtain
(26)limsupt→∞E|(N(t),P(t))|θ≤2θlimsupt→∞EV(N(t),P(t))≤2θC≜C(θ).
The proof is complete.Theorem 4.
The solutions of model (1) are stochastically ultimately bounded for any initial value (N0,P0)∈ℝ+2.Proof.
The solutions of model (1) will remain in ℝ+2 for all t≥0 with probability one from Theorem 2. And by Theorem 3, let θ=1/2; there is a C1>0 such that
(27)limsupt→∞E|(N(t),P(t))|1/2≤C1.
For any ɛ>0, set X=C12/ɛ2. Then by Chebyshev inequality,
(28)P{|(N(t),P(t))|>X}≤E|(N(t),P(t))|1/2X.
That is,
(29)limsupt→∞P{|(N(t),P(t))|>X}≤ɛ.
This completes the proof.Next, we study the asymptotic properties of the moment solutions of model (1).Theorem 5.
For any givenθ∈(0,1), there is a H=H(θ)>0; let (N(t),P(t)) be the solution of model (1) with any initial value (N0,P0); then,
(30)limsupt→∞1t∫0t‍E(Nθ+2(τ)+Pθ+2(τ))dτ≤H.Proof.
SetV(N,P):ℝ+2→ℝ+, and from Theorem 3, we have
(31)LV(N,P)≤rθNθ-α22θ(1-θ)Nθ×(N2+m2(b+N)2-2m(1-N)b+N-2N)+sθPθ-β22θ(1-θ)Pθ(c2P2-2c(1-d)P),
so,
(32)LV(N,P)+α2θ(1-θ)4Nθ+2+c2β2θ(1-θ)4Pθ+2≤rθNθ-α2θ(1-θ)4Nθ+2+α2θ(1-θ)Nθ(m(1-N)b+N+N-m22(b+N)2)+sθPθ-c2β2θ(1-θ)4Pθ+2+c(1-d)β2θ(1-θ)Pθ+1≤K,
where K is a positive number. Then
(33)LV(N,P)+k2γ2θ(1-θ)4(Nθ+2+Pθ+2)≤LV(N,P)+α2θ(1-θ)4Nθ+2+c2β2θ(1-θ)4Pθ+2≤K,
where k2γ2=min{α2,c2β2}; we can obtain
(34)EV(N(t),P(t))+k2γ2θ(1-θ)4×∫0t‍E(Nθ+2(τ)+Pθ+2(τ))dτ≤V(N0,P0)+Kt.
That is,
(35)∫0t‍E(Nθ+2(τ)+Pθ+2(τ))dτ≤4(V(N0,P0)+Kt)k2γ2θ(1-θ).
Set H=4K/(k2γ2θ(1-θ)); therefore,
(36)limsupt→∞1t∫0t‍E(Nθ+2(τ)+Pθ+2(τ))dτ≤H.


## 2.3. Stochastic Asymptotic Stability

From Theorem2, model (1) has a global positive solution. By constructing some Lyapunov functions, we analyze stochastic asymptotic stability of the positive equilibria of the stochastic model. The equilibria of the stochastic model (1) we consider are solutions of the following model:
(37)N(1-N-mb+N-wPa+N)=0,P(Na+N-d-cP)=0.
From (37), the positive equilibrium (N,P) satisfies the equation
(38)cN4+(2ac-c+cb)N3+(a2c-c(b-m)-2ac(1-b)+w(1-d))N2+(a2c(b-1)-2ac(b-m)+wb(1-d)-wad)N-a2c(b-m)-wabd=0,N>ad1-d>0.
As an example, we only give the proof of a positive equilibrium (N*,P*) of model (1).Theorem 6.
Let(39)θ=w(a+N*)(r-α2N*)as+β2acP*,
and
(40)A=r(mb(b+N*)+wP*a(a+N*)-1)+α2N*2(mb(b+N*)+wP*a(a+N*)-1)2+β2a2wP*(r-α2N*)2a3(a+N*)(s+β2cP*),B=α2w2N*2a2+cw(a+N*)(r-α2N*)(β2cP*-2s)2a(s+β2cP*),C=α2mwN*ab(b+N*)+α2w2N*P*a2(a+N*);
if A<0, C2-4AB<0, then the equilibrium point (N*,P*) of model (1) is stochastically asymptotically stable in ℝ+2; that is, for any initial value (N0,P0), the solution of model (1) has the property that
(41)limt→∞N(t)=N*,limt→∞P(t)=P*,
almost surelyProof.
Based on the theory of stability of stochastic differential equations [46], we only need to construct a Lyapunov function V(U) satisfying LV≤0 and the identity holds if and only if U=U*, where U=U(t) is the solution of the stochastic differential equation (3) and U* is the equilibrium point of model (3).
Define Lyapunov functions(42)V1(N)=N-N*-N*lnNN*,V2(P)=P-P*-P*lnPP*,
The nonnegativity of this function can be observed from x-1-lnx≥0 on x>0. If (N(t),P(t))∈ℝ+2, applying Itô’s formula, we have
(43)LV1=r(N-N*)(1-N-mb+N-wPa+N)+α2N*2(1-N-mb+N-wPa+N)2=[+wP*(a+N*)(a+N)-1)2r(m(b+N*)(b+N)+wP*(a+N*)(a+N)-1)+α2N*2(m(b+N*)(b+N)+wP*(a+N*)(a+N)-1)2]×(N-N*)2+α2w2N*2(a+N)2(P-P*)2-[rwa+N+α2wN*a+N×(m(b+N*)(b+N)+wP*(a+N*)(a+N)-1)rwa+N+α2wN*a+N]×(N-N*)(P-P*)≤[+α2N*2(mb(b+N*)+wP*a(a+N*)-1)2r(mb(b+N*)+wP*a(a+N*)-1)+α2N*2(mb(b+N*)+wP*a(a+N*)-1)2]×(N-N*)2+α2w2N*2a2(P-P*)2+α2wN*-rwa+N(N-N*)(P-P*)+(α2mwN*ab(b+N*)+α2w2N*P*a2(a+N*))|N-N*||P-P*|.
Similarly, we obtain
(44)LV2=s(P-P*)(Na+N-d-cP)+β2P*2(Na+N-d-cP)2=s(P-P*)(a(a+N)(a+N*)(N-N*)+c(P*-P)a(a+N)(a+N*))+β2P*2(a2(a+N)2(a+N*)2(N-N*)2+c2(P-P*)2+2ac(a+N)(a+N*)×(N-N*)(P-P*)a2(a+N)2(a+N*)2)=β2a2P*2(a+N)2(a+N*)2(N-N*)2+(β2c2P*2-sc)(P-P*)2+as+acβ2P*(a+N)(a+N*)(N-N*)(P-P*)≤β2a2P*2a2(a+N*)2(N-N*)2+(β2c2P*2-sc)(P-P*)2+as+β2acP*(a+N)(a+N*)(N-N*)(P-P*).
Define(45)V(N,P)=V1(N)+θV2(P),
and choose
(46)θ=w(a+N*)(r-α2N*)as+β2acP*.
We compute(47)LV(N,P)=LV1(N)+LθV2(P)≤[β2a2wP*(r-α2N*)2a2(a+N*)(as+β2acP*)r(mb(b+N*)+wP*a(a+N*)-1)+α2N*2(mb(b+N*)+wP*a(a+N*)-1)2+β2a2wP*(r-α2N*)2a2(a+N*)(as+β2acP*)](N-N*)2+[α2w2N*2a2+cw(a+N*)(r-α2N*)(β2cP*-2s)2(as+β2acP*)]×(P-P*)2+(α2mwN*ab(b+N*)+α2w2N*P*a2(a+N*))|N-N*|×|P-P*|=A(N-N*)2+B(P-P*)2+C|N-N*||P-P*|.
Let |U-U*|=(|N-N*|,|P-P*|)T; then, we have
(48)LV(N,P)⩽12|U-U*|T(2ACC2B)|U-U*|.
Clearly, if A<0 and C2-4AB<0 hold, then the above inequality implies LV(N,P)<0 along all trajectories in the first quadrant except (N*,P*). Then the desired assertion (41) follows immediately.


## 3. Numerical Simulations

In this section, we perform extensive numerical simulations of model (1) to confirm the analytical results by referring to the method mentioned in Higham [47]. Next, we consider the discretization equations
(49)Ni+1=Ni+rNi(1-Ni-mb+Ni-wPia+Ni)Δt+αNi(1-Ni-mb+Ni-wPia+Ni)Δtξ1i+α22Ni2(1-Ni-mb+Ni-wPia+Ni)2(ξ1i2-1)Δt,Pi+1=Pi+sPi(Nia+Ni-d-cPi)Δt+βPi(Nia+Ni-d-cPi)Δtξ2i+β22Pi2(Nia+Ni-d-cPi)2(ξ2i2-1)Δt,
where ξ1i and ξ2i(i=1,2,…,n) are the Gaussian random variables N(0,1).When choosing the environmental noise intensitiesα=β=0 (i.e., without noise), under the weak Allee effect, Figure 1 shows time-series plots and phase portraits of prey population N and predator population P of model (1) with initial value (N0,P0)=(0.2,1.4) and other parameters values r=1.5, s=0.6, a=0.5, b=0.6, c=0.1, d=0.12, and w=0.25. In Figure 1(a), with the Allee effect m=0.15, one can see that a positive equilibrium point (N*,P*)=(0.20589,1.71689) of model (1) is globally asymptotically stable, which illustrates the stable population distribution of both species. In Figure 1(b), increasing m to m=0.45, a stable limit cycle occurs which indicates that prey population coexist with predator population with oscillator balance behavior.(left) Time-series plots; (right) phase portraits for model (1) without noise (i.e., α=β=0) and with r=1.5, s=0.6, a=0.5, b=0.6, c=0.1, d=0.12, and w=0.25 and initial condition (N0,P0)=(0.2,1.4). (a) m=0.15 and (b) m=0.45.
(a)
(b)Figure2 shows time-series plots for model (1) with r=1.5, s=0.6, a=0.5, b=0.6, c=0.1, d=0.12, m=0.15, and w=0.25 and initial condition (N0,P0)=(0.2,1.4). The only difference between conditions of Figures 2(a) and 2(b) is that the values of environmental noise intensities α and β are different. As pointed out in Theorem 6, if those conditions hold, then the positive equilibrium point (N*,P*) is stochastically asymptotically stable in ℝ+2. In Figure 2(a), we choose α=β=0.5 and the conditions of Theorem 6 are satisfied. That is to say, the equilibrium point (N*,P*) is stochastically asymptotically stable. In Figure 2(b), choosing α=2.4 and β=0.8, one can see that two species N and P will die out. In other words, when the environmental noise intensities are large and violate the conditions of Theorem 6, the positive equilibrium (N*,P*) is no longer globally stable and not permanent.Time-series plots for model (1) with r=1.5, s=0.6, a=0.5, b=0.6, c=0.1, d=0.12, m=0.15, and w=0.25 and initial condition (N0,P0)=(0.2,1.4). (a) α=β=0.5 and (b) α=2.4 and β=0.8.
(a)
(b)By comparing Figure1(a) and Figure 2(a), one can realize that if the noise is not strong, the stochastic perturbation does not cause sharp changes of the dynamics of model (1). However, when the environmental noise is sufficiently large (see Figure 2(b)), the noise can force the population to become extinct.


## 4. Conclusions and Remarks

In this paper, we propose a stochastic version of a modified Holling type-II predator-prey model with an Allee effect in prey population and take into account density-dependence of predator population. We perturb the model with respect to white noise around its intrinsic growth rate of the prey and predator population, respectively. By constructing suitable Lyapunov functions and applying Itô formula, some qualitative properties are given, such as the existence of global positive solutions, stochastic boundedness, and stochastic asymptotic stability which means that in the biological point of view, the community consisting of both prey and predator species is a stable biotic community where all species will coexist.To our knowledge, the stability analysis of the stochastic Holling type-II predator-prey model with the Allee effect and density-dependence we illustrated here appears to be rare. And our complete analysis of the stochastic predator-prey model will give some suggestions to other models.



---

*Source: 340980-2013-11-20.xml*