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

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.

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 Section 3, 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−𝑁−−𝑏+𝑁𝑤𝑃𝑎+𝑁𝑟𝑑𝑡+𝛼𝑑𝜂₁(𝑡),𝑁𝑑𝑃=𝑃𝑎+𝑁−𝑑−𝑐𝑃𝑠𝑑𝑡+𝛽𝑑𝜂₂(𝑡).𝑁(0)=𝑁₀>0,𝑃(0)=𝑃₀>0. (1)

Here, we replace the intrinsic growth rates 𝑟 (>0) and 𝑠 (>0) of prey and predator by 𝑟⟶𝑟+𝛼𝜂̇₁(𝑡),𝑠⟶𝑠+𝛽𝜂̇₂(𝑡). (2)

In the above model, 𝑁=𝑁(𝑡) and 𝑃=𝑃(𝑡) represent prey and predator densities at time 𝑡, respectively. 𝑚𝑁/(𝑏+𝑁) is the term of additive Allee effect and 𝑚∈(0,1) and 𝑏∈(0,1) are Allee effect constants. If 𝑚<𝑏, the function of prey growth rate with the Allee effect of the prey is called weak Allee effect; if 𝑚>𝑏, the modified function of prey growth rate is strong Allee effect. 𝑐𝑃 stands for the density-dependence of predator population and 𝑐>0. 𝑤 (>0) is the predation rate, 𝑎 (>0) the half-saturation density of Holling type-II functional response, and 𝑑 (>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. 𝜂̇_𝑖(𝑡)(𝑖=1,2) is a standard white noise with increment Δ𝜂_𝑖(𝑡)=𝜂_𝑖(𝑡+Δ𝑡)−𝜂_𝑖(𝑡) that is independent Gaussian random variables 𝐍(0,Δ𝑡). That is to say, 𝜂_𝑖(𝑡)  (𝑖=1,2) is a Brownian motion defined on a complete probability space (Ω,ℱ,𝐏) with a filtration {ℱ_𝑡}_𝑡∈ℝ+ satisfying the usual conditions; namely, it is right continuous and increasing while ℱ₀ contains all 𝐏-null sets [46].

Consider model (1) as an Itô stochastic differential system of the type 𝑑𝑈(𝑡)=𝐹(𝑡,𝑈(𝑡))𝑑𝑡+𝐺(𝑡,𝑈(𝑡))𝑑𝜂(𝑡),𝑈𝑡₀=𝑈₀, (3) where 𝑈(𝑡)=(𝑁(𝑡),𝑃(𝑡))^𝑇,𝜂(𝑡)=𝜂₁(𝑡),𝜂₂(𝑡)^𝑇,𝑚𝐹=𝑟𝑁1−𝑁−−𝑏+𝑁𝑤𝑃𝑁𝑎+𝑁𝑠𝑃,𝑎+𝑁−𝑑−𝑐𝑃𝐺=^{𝑚𝛼𝑁1−𝑁−−𝑏+𝑁𝑤𝑃0}_𝑁^𝑎+𝑁_{0𝛽𝑃𝑎+𝑁−𝑑−𝑐𝑃}, (4) and the solution 𝑈(𝑡)  (𝑡>0) is an Itô process, 𝐹 is a slowly varying continuous component or drift coefficient, and 𝐺 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 (𝑁(𝑡),𝑃(𝑡)) for 𝑡∈[0,𝜏_𝑒) to model (1) a.s. for any given initial value (𝑁₀,𝑃₀)∈ℝ²₊={𝑁,𝑃∈ℝ²∣𝑁>0,𝑃>0}, where 𝜏_𝑒 is the explosion time.

Proof. Set 𝑝(𝑡)=ln𝑁(𝑡), 𝑞(𝑡)=ln𝑃(𝑡) and we consider the equations 𝑑𝑝=𝑟1−𝑒^𝑝−𝑚𝑏+𝑒^𝑝−𝑤𝑒^𝑞𝑎+𝑒^𝑝−𝛼²21−𝑒^𝑝−𝑚𝑏+𝑒^𝑝−𝑤𝑒^𝑞𝑎+𝑒^𝑝2𝑑𝑡+𝛼1−𝑒^𝑝−𝑚𝑏+𝑒^𝑝−𝑤𝑒^𝑞𝑎+𝑒^𝑝𝑑𝜂₁(𝑡),𝑠𝑒𝑑𝑞=^𝑝𝑎+𝑒^𝑝−𝑑−𝑐𝑒^𝑞−𝛽²2𝑒^𝑝𝑎+𝑒^𝑝−𝑑−𝑐𝑒^𝑞2𝑒𝑑𝑡+𝛽^𝑝𝑎+𝑒^𝑝−𝑑−𝑐𝑒^𝑞𝑑𝜂₂(𝑡),𝑝(0)=ln𝑁₀,𝑞(0)=ln𝑃₀ (2.1) on 𝑡≥0. Then we know that the coefficients of (2.1) satisfy the local Lipschitz condition, and there is a unique local solution (𝑝(𝑡),𝑞(𝑡)) on [0,𝜏_𝑒). By Itô’s formula, we can see that 𝑁(𝑡)=𝑒^𝑝(𝑡), 𝑃(𝑡)=𝑒^𝑞(𝑡) is the unique positive local solution to model (1) with the initial value 𝑁₀>0, 𝑃₀>0.

Next, we show that the solution of model (1) is global; that is, 𝜏_𝑒=∞ a.s.

Theorem 2. For model (1) and any given initial value (𝑁₀,𝑃₀)∈ℝ²₊, there is a unique solution (𝑁(𝑡),𝑃(𝑡)) on 𝑡≥0 and the solution will remain in ℝ²₊ with probability one.

Proof. Let 𝑛₀>0 be sufficiently large for 𝑁₀ and 𝑃₀ lying within the interval [1/𝑛₀,𝑛₀]. For each integer 𝑛≥𝑛₀, let us define a sequence of stopping times by 𝜏_𝑛=inf𝑡∈0,𝜏_𝑒(𝑡)∉1∶𝑁𝑛(𝑡)∉1,𝑛or𝑃𝑛,𝑛, (6) where we set inf∅=∞ (∅ represents the empty set) in this paper. Since 𝜏_𝑛 is nondecreasing as 𝑛→∞, there exists the limit 𝜏_∞=lim_𝑛→∞𝜏_𝑛. (7) Then 𝜏_∞≤𝜏_𝑒a.s. (8) Now, we need to show 𝜏_∞=∞ a.s. If this statement is violated, then there exist 𝑇>0 and 𝜀∈(0,1) such that 𝐏𝜏_∞≤𝑇>𝜀. (9) Thus, there is an integer 𝑛₁≥𝑛₀ such that 𝐏𝜏_𝑛≤𝑇≥𝜀∀𝑛≥𝑛₁. (10)
Define a ℂ²-function 𝑉∶ℝ²₊→ℝ₊ by 𝑉(𝑁,𝑃)=𝑁+𝑃−ln𝑁𝑃−2, (11) which is a nonnegative function. If (𝑁(𝑡),𝑃(𝑡))∈ℝ²₊, by using Itô’s formula, we compute ()=()𝑚𝑑𝑉𝑁,𝑃2𝑟𝑁−11−𝑁−−𝑏+𝑁𝑤𝑃()𝑁𝑎+𝑁+2𝑠𝑃−1+𝛼𝑎+𝑁−𝑑−𝑐𝑃₂2𝑚1−𝑁−−𝑏+𝑁𝑤𝑃𝑎+𝑁²+𝛽²2𝑁𝑎+𝑁−𝑑−𝑐𝑃²(𝑚𝑑𝑡+𝛼𝑁−1)1−𝑁−−𝑏+𝑁𝑤𝑃𝑎+𝑁𝑑𝜂₁(𝑁+𝛽𝑃−1)𝑎+𝑁−𝑑−𝑐𝑃𝑑𝜂₂≤()𝑚2𝑟1−𝑁+𝑏+𝑁𝑤𝑃+𝛼𝑎+𝑁₂2𝑚1−𝑁−−𝑏+𝑁𝑤𝑃𝑎+𝑁²+𝑑𝑡2𝑠𝑁𝑃+𝛽𝑎+𝑁+𝑑+𝑐𝑃₂2𝑁𝑎+𝑁−𝑑−𝑐𝑃²()𝑚𝑑𝑡+𝛼𝑁−11−𝑁−−𝑏+𝑁𝑤𝑃𝑎+𝑁𝑑𝜂₁()𝑁+𝛽𝑃−1𝑎+𝑁−𝑑−𝑐𝑃𝑑𝜂₂≤𝑀₁+𝑀₂()×𝑚+𝛼𝑁−11−𝑁−−𝑏+𝑁𝑤𝑃𝑎+𝑁𝑑𝜂₁(𝑁+𝛽𝑃−1)𝑎+𝑁−𝑑−𝑐𝑃𝑑𝜂₂, (12) where 𝑀₁, 𝑀₂ are positive numbers. Integrating both sides of the above inequality from 0 to 𝜏_𝑛∧𝑇 and taking the expectations give rise to 𝐄𝑉𝑁𝜏_𝑛∧𝑇,𝑃𝜏_𝑛∧𝑇≤𝑉𝑁₀,𝑃₀+𝑀₁+𝑀₂𝑇. (13) Set Ω_𝑛={𝜏_𝑛≤𝑇} and by inequality (10), we get 𝐏(Ω_𝑛)≥𝜀. For every 𝜇∈Ω_𝑛, 𝑁(𝜏_𝑛,𝜇) and 𝑃(𝜏_𝑛,𝜇) equal either 𝑛 or 1/𝑛, hence 𝑉𝑁𝜏_𝑛,𝜇,𝑃𝜏_𝑛1,𝜇≥min𝑛−ln𝑛−1,𝑛+ln𝑛−1. (14) It then follows from (13) that 𝑉𝑁₀,𝑃₀+𝑀₁+𝑀₂𝑇≥𝐄1_Ω𝑛𝜇𝑉𝑁𝜏_𝑛,𝑃𝜏_𝑛1≥𝜀min𝑛−ln𝑛−1,𝑛+ln𝑛−1, (15) where 1_Ω𝑛 is the indicator function of Ω_𝑛. Letting 𝑛→∞ leads to the contradiction ∞>𝑉(𝑁₀,𝑃₀)+(𝑀₁+𝑀₂)𝑇=∞. This completes the proof.

2.2. Stochastic Boundedness

Theorem 3. For any 𝜃∈(0,1), there is a positive constant 𝐶=𝐶(𝜃), which is independent of initial value (𝑁₀,𝑃₀)∈ℝ²₊, such that the solution of model (1) has the property that limsup_𝑡→∞𝐄|(𝑁(𝑡),𝑃(𝑡))|^𝜃≤𝐶. (16)

Proof. Define a function by 𝑉(𝑁,𝑃)=𝑁^𝜃+𝑃^𝜃,(𝑁,𝑃)∈ℝ²₊. (17) For any 𝑡>0, using Itô formula gives that 𝐋𝑉(𝑁,𝑃)=𝑟𝜃𝑁^𝜃𝑚1−𝑁−−𝑏+𝑁𝑤𝑃𝑎+𝑁+𝑠𝜃𝑃^𝜃𝑁+𝛼𝑎+𝑁−𝑑−𝑐𝑃₂2𝜃(𝜃−1)𝑁^𝜃𝑚1−𝑁−−𝑏+𝑁𝑤𝑃𝑎+𝑁²+𝛽²2𝜃(𝜃−1)𝑃^𝜃𝑁𝑎+𝑁−𝑑−𝑐𝑃², (18) as 0<𝜃<1; then, 𝐋𝑉(𝑁,𝑃)≤𝑟𝜃𝑁^𝜃()−𝛼1−𝑁²2𝜃(1−𝜃)𝑁^𝜃×𝑚1−𝑁−−𝑏+𝑁𝑤𝑃𝑎+𝑁²+𝑠𝜃𝑃^𝜃𝑁−𝛽𝑎+𝑁₂2𝜃(1−𝜃)𝑃^𝜃𝑁𝑎+𝑁−𝑑−𝑐𝑃²≤𝑟𝜃𝑁^𝜃−𝛼²2𝜃(1−𝜃)𝑁^𝜃×𝑚²(𝑏+𝑁)²+(1−𝑁)1−𝑁−2𝑚𝑏+𝑁−1+𝑠𝜃𝑃^𝜃−𝛽²2𝜃(1−𝜃)𝑃^𝜃𝑐²𝑃²−2𝑐(1−𝑑)𝑃=𝑟𝜃𝑁^𝜃−𝛼²2𝜃(1−𝜃)𝑁^𝜃×𝑚²(𝑏+𝑁)²+()×1−𝑁1−𝑁−2𝑚𝑏+𝑁−1+𝑁^𝜃+𝑠𝜃𝑃^𝜃−𝛽²2𝜃(1−𝜃)𝑃^𝜃×𝑐²𝑃²−2𝑐(1−𝑑)𝑃+𝑃^𝜃−𝑉(𝑁,𝑃)≤𝐶₀−𝑉(𝑁,𝑃), (19) where 𝐶₀ is an integer. Hence, 𝑑𝑉(𝑁,𝑃)=𝐋𝑉(𝑁,𝑃)𝑑𝑡+𝛼𝜃𝑁^𝜃𝑚1−𝑁−−𝑏+𝑁𝑤𝑃𝑎+𝑁𝑑𝜂₁+𝛽𝜃𝑃^𝜃𝑁𝑎+𝑁−𝑑−𝑐𝑃𝑑𝜂₂≤𝐶₀−𝑉(𝑁,𝑃)𝑑𝑡+𝛼𝜃𝑁^𝜃𝑚1−𝑁−−𝑏+𝑁𝑤𝑃𝑎+𝑁𝑑𝜂₁+𝛽𝜃𝑃^𝜃𝑁𝑎+𝑁−𝑑−𝑐𝑃𝑑𝜂₂. (20)
Applying Itô formula yields 𝑑𝑒^𝑡𝑉(𝑁,𝑃)=𝑒^𝑡(𝑉(𝑁,𝑃)𝑑𝑡+𝑑𝑉(𝑁,𝑃))≤𝑒^𝑡𝐶₀𝑑𝑡+𝑒^𝑡𝛼𝜃𝑁^𝜃×𝑚1−𝑁−−𝑏+𝑁𝑤𝑃𝑎+𝑁𝑑𝜂₁+𝑒^𝑡𝛽𝜃𝑃^𝜃𝑁𝑎+𝑁−𝑑−𝑐𝑃𝑑𝜂₂. (21) Using the sequence of stopping times {𝜏_𝑛}^∞_𝑛=1 defined in the proof of Theorem 2 and from (21), we have 𝐄𝑒^{𝑡∧𝜏𝑛}𝑉𝑁𝑡∧𝜏_𝑛,𝑃𝑡∧𝜏_𝑛≤𝑉𝑁₀,𝑃₀+𝐶₀𝐄𝑒^{𝑡∧𝜏𝑛}−1∀𝑡≥0. (22) Letting 𝑛→∞ in the latter inequality with a fact that 𝑉(𝑁(𝑡∧𝜏_𝑛),𝑃(𝑡∧𝜏_𝑛))>0 and 0<𝑒^{𝑡∧𝜏𝑛}≤𝑒^𝑡 a.s. and using Fatou’s lemma, we obtain 𝑒^𝑡𝐄𝑉(𝑁,𝑃)≤𝑉𝑁₀,𝑃₀+𝐶₀𝑒^𝑡−1. (23) Therefore, limsup_𝑡→∞𝐄𝑉(𝑁(𝑡),𝑃(𝑡))≤𝐶. (24)
On the other hand, since |𝑁,𝑃|^𝜃≤2^𝜃max𝑁^𝜃,𝑃^𝜃≤2^𝜃𝑉(𝑁,𝑃), (25) we obtain limsup_𝑡→∞𝐄|(𝑁(𝑡),𝑃(𝑡))|^𝜃≤2^𝜃limsup_𝑡→∞(𝑁(𝑡)(𝑡≤𝐄𝑉,𝑃))2^𝜃𝐶≜𝐶(𝜃). (26) The proof is complete.

Theorem 4. The solutions of model (1) are stochastically ultimately bounded for any initial value (𝑁₀,𝑃₀)∈ℝ²₊.

Proof. The solutions of model (1) will remain in ℝ²₊ for all 𝑡≥0 with probability one from Theorem 2. And by Theorem 3, let 𝜃=1/2; there is a 𝐶₁>0 such that limsup_𝑡→∞𝐄|(𝑁(𝑡),𝑃(𝑡))|^1/2≤𝐶₁. (27) For any 𝜀>0, set 𝑋=𝐶²₁/𝜀². Then by Chebyshev inequality, 𝐏𝑁(𝑡)(𝑡}≤𝐄𝑁(𝑡)(𝑡{|(,𝑃))|>𝑋|(,𝑃))|^1/2√𝑋. (28) That is, limsup_𝑡→∞𝐏{|(𝑁(𝑡),𝑃(𝑡))|>𝑋}≤𝜀. (29) 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 𝐻=𝐻(𝜃)>0; let (𝑁(𝑡),𝑃(𝑡)) be the solution of model (1) with any initial value (𝑁₀,𝑃₀); then, limsup_𝑡→∞1𝑡^𝑡₀𝐄𝑁^𝜃+2(𝜏)+𝑃^𝜃+2(𝜏)𝑑𝜏≤𝐻. (30)

Proof. Set 𝑉(𝑁,𝑃)∶ℝ²₊→ℝ₊, and from Theorem 3, we have 𝐋𝑉(𝑁,𝑃)≤𝑟𝜃𝑁^𝜃−𝛼²2𝜃(1−𝜃)𝑁^𝜃×𝑁²+𝑚²(𝑏+𝑁)²−2𝑚(1−𝑁)𝑏+𝑁−2𝑁+𝑠𝜃𝑃^𝜃−𝛽²2𝜃(1−𝜃)𝑃^𝜃𝑐²𝑃²−2𝑐(1−𝑑)𝑃, (31) so, ()+𝛼𝐋𝑉𝑁,𝑃²𝜃(1−𝜃)4𝑁^𝜃+2+𝑐²𝛽²𝜃(1−𝜃)4𝑃^𝜃+2≤𝑟𝜃𝑁^𝜃−𝛼²𝜃(1−𝜃)4𝑁^𝜃+2+𝛼²𝜃(1−𝜃)𝑁^𝜃𝑚(1−𝑁)𝑚𝑏+𝑁+𝑁−²2(𝑏+𝑁)²+𝑠𝜃𝑃^𝜃−𝑐²𝛽²𝜃(1−𝜃)4𝑃^𝜃+2+𝑐(1−𝑑)𝛽²𝜃(1−𝜃)𝑃^𝜃+1≤𝐾, (32) where 𝐾 is a positive number. Then ()+𝑘𝐋𝑉𝑁,𝑃²𝛾²𝜃(1−𝜃)4𝑁^𝜃+2+𝑃^𝜃+2()+𝛼≤𝐋𝑉𝑁,𝑃²𝜃(1−𝜃)4𝑁^𝜃+2+𝑐²𝛽²𝜃(1−𝜃)4𝑃^𝜃+2≤𝐾, (33) where 𝑘²𝛾²=min{𝛼²,𝑐²𝛽²}; we can obtain (𝑁(𝑡)(𝑡+𝑘𝐄𝑉,𝑃))²𝛾²𝜃(1−𝜃)4×∫^𝑡₀𝐄𝑁^𝜃+2(𝜏)+𝑃^𝜃+2(𝜏)𝑁𝑑𝜏≤𝑉₀,𝑃₀+𝐾𝑡. (34) That is, ^𝑡₀𝐄𝑁^𝜃+2(𝜏)+𝑃^𝜃+2(𝜏)4𝑉𝑁𝑑𝜏≤⁰,𝑃₀+𝐾𝑡𝑘²𝛾²𝜃().1−𝜃 (35) Set 𝐻=4𝐾/(𝑘²𝛾²𝜃(1−𝜃)); therefore, limsup_𝑡→∞1𝑡^𝑡₀𝐄𝑁^𝜃+2(𝜏)+𝑃^𝜃+2(𝜏)𝑑𝜏≤𝐻. (36)

2.3. Stochastic Asymptotic Stability

From Theorem 2, 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: 𝑁𝑚1−𝑁−−𝑏+𝑁𝑤𝑃𝑃𝑁𝑎+𝑁=0,𝑎+𝑁−𝑑−𝑐𝑃=0. (37) From (37), the positive equilibrium (𝑁,𝑃) satisfies the equation 𝑐𝑁⁴+(2𝑎𝑐−𝑐+𝑐𝑏)𝑁³+𝑎²()()(𝑁𝑐−𝑐𝑏−𝑚−2𝑎𝑐1−𝑏+𝑤1−𝑑))²+𝑎²𝑐()()())𝑁𝑏−1−2𝑎𝑐𝑏−𝑚+𝑤𝑏1−𝑑−𝑤𝑎𝑑−𝑎²𝑐(𝑏−𝑚)−𝑤𝑎𝑏𝑑=0,𝑁>𝑎𝑑1−𝑑>0. (38) As an example, we only give the proof of a positive equilibrium (𝑁^∗,𝑃^∗) of model (1).

Theorem 6. Let 𝑤𝜃=𝑎+𝑁^∗𝑟−𝛼²𝑁^∗𝑎𝑠+𝛽²𝑎𝑐𝑃^∗, (39) and 𝑚𝐴=𝑟𝑏(𝑏+𝑁^∗)+𝑤𝑃^∗𝑎(𝑎+𝑁^∗)+𝛼−1₂𝑁^∗2𝑚𝑏(𝑏+𝑁^∗)+𝑤𝑃^∗𝑎(𝑎+𝑁^∗)−1²+𝛽²𝑎²𝑤𝑃^∗𝑟−𝛼²𝑁^∗2𝑎³(𝑎+𝑁^∗)𝑠+𝛽²𝑐𝑃^∗,𝛼𝐵=²𝑤²𝑁^∗2𝑎²+𝑐𝑤𝑎+𝑁^∗𝑟−𝛼²𝑁^∗𝛽²𝑐𝑃^∗−2𝑠2𝑎𝑠+𝛽²𝑐𝑃^∗,𝛼𝐶=²𝑚𝑤𝑁^∗𝑎𝑏(𝑏+𝑁^∗)+𝛼²𝑤²𝑁^∗𝑃^∗𝑎²(𝑎+𝑁^∗); (40) if 𝐴<0, 𝐶²−4𝐴𝐵<0, then the equilibrium point (𝑁^∗,𝑃^∗) of model (1) is stochastically asymptotically stable in ℝ²₊; that is, for any initial value (𝑁₀,𝑃₀), the solution of model (1) has the property that lim_𝑡→∞𝑁(𝑡)=𝑁^∗,lim_𝑡→∞𝑃(𝑡)=𝑃^∗, (41) almost surely

Proof. Based on the theory of stability of stochastic differential equations [46], we only need to construct a Lyapunov function 𝑉(𝑈) satisfying 𝐋𝑉≤0 and the identity holds if and only if 𝑈=𝑈^∗, where 𝑈=𝑈(𝑡) is the solution of the stochastic differential equation (3) and 𝑈^∗ is the equilibrium point of model (3).
Define Lyapunov functions 𝑉₁(𝑁)=𝑁−𝑁^∗−𝑁^∗𝑁ln𝑁^∗,𝑉₂(𝑃)=𝑃−𝑃^∗−𝑃^∗𝑃ln𝑃^∗, (42) The nonnegativity of this function can be observed from 𝑥−1−ln𝑥≥0 on 𝑥>0. If (𝑁(𝑡),𝑃(𝑡))∈ℝ²₊, applying Itô’s formula, we have 𝐋𝑉₁=𝑟𝑁−𝑁^∗𝑚1−𝑁−−𝑏+𝑁𝑤𝑃+𝛼𝑎+𝑁₂𝑁^∗2𝑚1−𝑁−−𝑏+𝑁𝑤𝑃𝑎+𝑁²=𝑟𝑚(𝑏+𝑁^∗)+)(𝑏+𝑁𝑤𝑃^∗(𝑎+𝑁^∗)+𝛼)(𝑎+𝑁−1₂𝑁^∗2𝑚(𝑏+𝑁^∗)+)(𝑏+𝑁𝑤𝑃^∗(𝑎+𝑁^∗)(𝑎+𝑁)−1²×𝑁−𝑁^∗2+𝛼²𝑤²𝑁^∗2(𝑎+𝑁)²𝑃−𝑃^∗2−𝑟𝑤+𝛼𝑎+𝑁²𝑤𝑁^∗×𝑚𝑎+𝑁(𝑏+𝑁^∗)+)(𝑏+𝑁𝑤𝑃^∗(𝑎+𝑁^∗)×)(𝑎+𝑁−1𝑁−𝑁^∗𝑃−𝑃^∗≤𝑟𝑚𝑏(𝑏+𝑁^∗)+𝑤𝑃^∗𝑎(𝑎+𝑁^∗)+𝛼−1₂𝑁^∗2𝑚𝑏(𝑏+𝑁^∗)+𝑤𝑃^∗𝑎(𝑎+𝑁^∗)−1²×𝑁−𝑁^∗2+𝛼²𝑤²𝑁^∗2𝑎²𝑃−𝑃^∗2+𝛼²𝑤𝑁^∗−𝑟𝑤𝑎+𝑁𝑁−𝑁^∗𝑃−𝑃^∗+𝛼²𝑚𝑤𝑁^∗𝑎𝑏(𝑏+𝑁^∗)+𝛼²𝑤²𝑁^∗𝑃^∗𝑎²(𝑎+𝑁^∗)𝑁−𝑁^∗𝑃−𝑃^∗. (43) Similarly, we obtain 𝐋𝑉₂=𝑠𝑃−𝑃^∗𝑁+𝛽𝑎+𝑁−𝑑−𝑐𝑃₂𝑃^∗2𝑁𝑎+𝑁−𝑑−𝑐𝑃²=𝑠𝑃−𝑃^∗𝑎(𝑎+𝑁)(𝑎+𝑁^∗)𝑁−𝑁^∗+𝑐𝑃^∗+𝛽−𝑃₂𝑃^∗2𝑎²(𝑎+𝑁)²(𝑎+𝑁^∗)²𝑁−𝑁^∗2+𝑐²𝑃−𝑃^∗2+2𝑎𝑐(𝑎+𝑁)(𝑎+𝑁^∗)×𝑁−𝑁^∗𝑃−𝑃^∗=𝛽²𝑎²𝑃^∗2(𝑎+𝑁)²(𝑎+𝑁^∗)²𝑁−𝑁^∗2+𝛽²𝑐²𝑃^∗2−𝑠𝑐𝑃−𝑃^∗2+𝑎𝑠+𝑎𝑐𝛽²𝑃^∗(𝑎+𝑁)(𝑎+𝑁^∗)𝑁−𝑁^∗𝑃−𝑃^∗≤𝛽²𝑎²𝑃^∗2𝑎²(𝑎+𝑁^∗)²𝑁−𝑁^∗2+𝛽²𝑐²𝑃^∗2−𝑠𝑐𝑃−𝑃^∗2+𝑎𝑠+𝛽²𝑎𝑐𝑃^∗(𝑎+𝑁)(𝑎+𝑁^∗)𝑁−𝑁^∗𝑃−𝑃^∗. (44)
Define 𝑉(𝑁,𝑃)=𝑉₁(𝑁)+𝜃𝑉₂(𝑃), (45) and choose 𝑤𝜃=𝑎+𝑁^∗𝑟−𝛼²𝑁^∗𝑎𝑠+𝛽²𝑎𝑐𝑃^∗. (46)
We compute 𝐋𝑉(𝑁,𝑃)=𝐋𝑉₁(𝑁)+𝐋𝜃𝑉₂(𝑃)≤𝑟𝑚𝑏(𝑏+𝑁^∗)+𝑤𝑃^∗𝑎(𝑎+𝑁^∗)+𝛼−1₂𝑁^∗2𝑚𝑏(𝑏+𝑁^∗)+𝑤𝑃^∗𝑎(𝑎+𝑁^∗)−1²+𝛽²𝑎²𝑤𝑃^∗𝑟−𝛼²𝑁^∗2𝑎²(𝑎+𝑁^∗)𝑎𝑠+𝛽²𝑎𝑐𝑃^∗𝑁−𝑁^∗2+𝛼²𝑤²𝑁^∗2𝑎²+𝑐𝑤𝑎+𝑁^∗𝑟−𝛼²𝑁^∗𝛽²𝑐𝑃^∗−2𝑠2𝑎𝑠+𝛽²𝑎𝑐𝑃^∗×𝑃−𝑃^∗2+𝛼²𝑚𝑤𝑁^∗𝑎𝑏(𝑏+𝑁^∗)+𝛼²𝑤²𝑁^∗𝑃^∗𝑎²(𝑎+𝑁^∗)𝑁−𝑁^∗×𝑃−𝑃^∗=𝐴𝑁−𝑁^∗2+𝐵𝑃−𝑃^∗2+𝐶𝑁−𝑁^∗𝑃−𝑃^∗. (47) Let |𝑈−𝑈^∗|=(|𝑁−𝑁^∗|,|𝑃−𝑃^∗|)^𝑇; then, we have ()⩽1𝐋𝑉𝑁,𝑃2𝑈−𝑈^∗𝑇2𝐴𝐶𝐶2𝐵𝑈−𝑈^∗. (48) Clearly, if 𝐴<0 and 𝐶²−4𝐴𝐵<0 hold, then the above inequality implies 𝐋𝑉(𝑁,𝑃)<0 along all trajectories in the first quadrant except (𝑁^∗,𝑃^∗). 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 𝑁_𝑖+1=𝑁_𝑖+𝑟𝑁_𝑖1−𝑁_𝑖−𝑚𝑏+𝑁_𝑖−𝑤𝑃_𝑖𝑎+𝑁_𝑖Δ𝑡+𝛼𝑁_𝑖1−𝑁_𝑖−𝑚𝑏+𝑁_𝑖−𝑤𝑃_𝑖𝑎+𝑁_𝑖√Δ𝑡𝜉_1𝑖+𝛼²2𝑁²_𝑖1−𝑁_𝑖−𝑚𝑏+𝑁_𝑖−𝑤𝑃_𝑖𝑎+𝑁_𝑖²𝜉²_1𝑖𝑃−1Δ𝑡,_𝑖+1=𝑃_𝑖+𝑠𝑃_𝑖𝑁_𝑖𝑎+𝑁_𝑖−𝑑−𝑐𝑃_𝑖Δ𝑡+𝛽𝑃_𝑖𝑁_𝑖𝑎+𝑁_𝑖−𝑑−𝑐𝑃_𝑖√Δ𝑡𝜉_2𝑖+𝛽²2𝑃²_𝑖𝑁_𝑖𝑎+𝑁_𝑖−𝑑−𝑐𝑃_𝑖²𝜉²_2𝑖−1Δ𝑡, (49) where 𝜉_1𝑖 and 𝜉_2𝑖(𝑖=1,2,…,𝑛) are the Gaussian random variables 𝐍(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 𝑁 and predator population 𝑃 of model (1) with initial value (𝑁₀,𝑃₀)=(0.2,1.4) and other parameters values 𝑟=1.5, 𝑠=0.6, 𝑎=0.5, 𝑏=0.6, 𝑐=0.1, 𝑑=0.12, and 𝑤=0.25. In Figure 1(a), with the Allee effect 𝑚=0.15, one can see that a positive equilibrium point (𝑁^∗,𝑃^∗)=(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 𝑚 to 𝑚=0.45, a stable limit cycle occurs which indicates that prey population coexist with predator population with oscillator balance behavior.

(a)

(b)

Figure 1 

(left) Time-series plots; (right) phase portraits for model (1) without noise (i.e., 𝛼=𝛽=0) and with 𝑟=1.5, 𝑠=0.6, 𝑎=0.5, 𝑏=0.6, 𝑐=0.1, 𝑑=0.12, and 𝑤=0.25 and initial condition (𝑁0,𝑃0)=(0.2,1.4). (a) 𝑚=0.15 and (b) 𝑚=0.45.

Figure 2 shows time-series plots for model (1) with 𝑟=1.5, 𝑠=0.6, 𝑎=0.5, 𝑏=0.6, 𝑐=0.1, 𝑑=0.12, 𝑚=0.15, and 𝑤=0.25 and initial condition (𝑁₀,𝑃₀)=(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 (𝑁^∗,𝑃^∗) is stochastically asymptotically stable in ℝ²₊. In Figure 2(a), we choose 𝛼=𝛽=0.5 and the conditions of Theorem 6 are satisfied. That is to say, the equilibrium point (𝑁^∗,𝑃^∗) is stochastically asymptotically stable. In Figure 2(b), choosing 𝛼=2.4 and 𝛽=0.8, one can see that two species 𝑁 and 𝑃 will die out. In other words, when the environmental noise intensities are large and violate the conditions of Theorem 6, the positive equilibrium (𝑁^∗,𝑃^∗) is no longer globally stable and not permanent.

(a)

(b)

Figure 2 

Time-series plots for model (1) with 𝑟=1.5, 𝑠=0.6, 𝑎=0.5, 𝑏=0.6, 𝑐=0.1, 𝑑=0.12, 𝑚=0.15, and 𝑤=0.25 and initial condition (𝑁0,𝑃0)=(0.2,1.4). (a) 𝛼=𝛽=0.5 and (b) 𝛼=2.4 and 𝛽=0.8.

By comparing Figure 1(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.

Acknowledgments

The author would like to thank the editor and the anonymous referees for their helpful comments and suggestions.