# The Long Time Behavior of Equilibrium Status of a Predator-Prey System with Delayed Fear in Deterministic and Stochastic Scenarios

Authors: Weili Kong; Yuanfu Shao
Journal: Journal of Mathematics
 (2022)

---

## Abstract

In view of a time lag between the time that the prey perceive predator signals and make some changes or behavior responses, we establish a predator–prey system with direct and indirect predation in this paper. First, we investigate the existence, boundedness, and global asymptotical stability of the positive equilibrium status. Next, by perturbing the mortality rates of prey species and predator species, we stretch the deterministic system to the stochastic scenario and investigate the existence of stochastic process and the global asymptotical stability of the equilibrium status. The analytical findings show that fear affects the value of the equilibrium status, and stochastic disturbance affects its stability, but time delay has no effect if some conditions are satisfied, which are verified by some examples and numerical simulations.


---

## Body

## 1. Introduction

For predator and prey interaction model, the direct predation (the depletion of prey species by predator species) is very popular and is easy to be detected in this field. In mathematical modeling, many classes of functional responses are introduced to reveal the impact of straight predation by predator on the dynamics of ecological models, see [1–5] and references cited therein. Different from the direct predation, the indirect predation (the fear from predator on prey) is comparatively difficult to be observed and, hence, is studied very little, whereas many experiments imply that the fear from predator on prey sometimes can change the prey’s demography and cannot be ignored. For example, when there is no straight killing, Zanette et al. [6] did an experiment on songbirds so as to survey the impact of fear on songbirds’s reproduction, which was performed during a total breeding period. Their results show that due to the fear from predator, the female songbirds laid and hatched fewer eggs resulting in decreasing the survival rate of the descendants significantly. Particularly, the number of descendants of the songbirds decline 40%. In real world, when prey populations perceive some predation risks, they will always show a few anti-predator actions, such as vigilance, foraging actions, and some psychological changes, which cut down the growth speed of the prey species and the survival speed of adults is also impacted accordingly. Recently, many predator-prey population systems containing fear of prey species are studied and some valuable results are reported [7–9].Time delay is very popular in biological systems. For populations, they will always take some time to complete such processes as digestion of food, maturation, prey hunting, disease transmission [10,11]. Similarly, when fear from predators appears, there is a time delay for prey to make response to the predation danger after they perceive predator signals through chemical or vocal cues [12].The deterministic systems are very popular in real world, but they are exposed to the survival environments without exception, and hence, the environmental white noise will bring large influence to their dynamics. The coefficients in the deterministic model cannot catch the impact of stochastic environmental disturbances, and it is necessary to spread the deterministic situation into a stochastic scenario. Many scholars incorporated white noise to investigate the impacts of population disturbance in the system dynamics of the interconnected populations [13–17].Stimulated by above analysis, we establish a deterministic predator–prey system together with Holling-II type functional response and fear effect of prey from predator species, where there is a time delay for prey’s response between their perceiving predation risk and the behavioral changes. Then by considering multiple environmental white noise, we extend it to stochastic scenario. To the author’s knowledge, besides direct predation, there are few findings reported on the dynamics of such stochastic models together with fear factor and delayed response.The paper is structured as follows: The boundedness of solutions, existence, and stability of the equilibrium status of the deterministic system (2) are discussed in Section 2. The existence of a unique stochastic process and the long time behavior of the equilibrium status of stochastic model (21) are investigated in Section 3. Then some examples and simulation analysis are carried out to validate the theoretical findings in Section 4. Finally, a brief conclusion is presented to conclude this work in Section 5.


## 2. Deterministic Model

The following predator–prey system model with direct functional response of predator (Holling-II type) and indirect predation (fear from predator) is proposed in Reference [9]:(1)dxtdt=xtr1+fyt−d1−bxt−pyt1+axt,dytdt=yt−d2−hyt+cpxt1+axt.The xt denotes the biomass of prey and yt denotes the biomass of predator population at time t. Coefficient f denotes the fear level of the prey induced by predator population, r denotes the birth rate of prey population, d1 denotes the mortality rate of prey population, b denotes the mortality rate of prey population due to intraspecies competition, p denotes the predation rate, p/a is the maximum number of prey eaten per predator per unit of time, d2 denotes the common mortality rate of predator population, h denotes the mortality rate of predator population due to intraspecies competition, c is the conversion coefficient. More biological backgrounds of this model are referred to [9].In view of a time lag between the time of perceiving predator signals and behavioral responses in the prey population, the above model becomes(2)dxtdt=xtr1+fyt−τ−d1−bxt−pyt1+axt,dytdt=yt−d2−hyt+cpxt1+axt.The initial conditions of (2) are taken as(3)xθ=ϕ1θ>0,yθ=ϕ2θ>0,−τ≤θ≤0,where the biological meanings of all parameters are the same as given before, τ is a time delay of the prey’s response to the fear induced by predator, ϕθ=ϕ1θ,ϕ2θT is a continuous mapping from −τ,0 to R+2. For any ϕ, we let ϕ=supθ∈−τ,0ϕθ, where ϕθ is the usual distance in R2. All coefficients in system (2) are assumed to be positive for biological backgrounds.


### 2.1. Positivity and Boundedness of Solutions

For assurance that model (2) is biologically well-behaved, we want to show that solutions of (2) are positive and bounded. For the first equation of (2), we integrate and get(4)xt=x0exp∫0tr1+fys−τ−d1−bxs−pys1+axsds>0.Applying the same integrating to the second equation of (2) yields(5)yt=y0exp∫0t−d2−hys+c2xs1+axsds>0.Therefore, all solutions of system (2) are positive. Now, we are to discuss the boundedness by use of comparison method. We can obtain from the first equation of (2) that(6)dxtdt≤xtr−d1−bxt.Integrating the above inequality and by use of comparison theory, then we have(7)xt≤r−d1er−d1t1+ber−d1t,which implies(8)lim supt⟶∞xt≤r−d1b.DefineVt=cxt+yt, then(9)V′+μV≤cxr−d1−bx−yd2+hy≤cxr+μ−d1−bx−yd2−μ≤cr+μ−d124b≐ϱ,where μ is a constant such that μ<d2. By the differential inequality theory, we obtain that 0≤V≤ϱ1−eμt/μ+Vx0,y0e−μt. Letting t⟶∞, then 0≤V≤ϱ/μ holds. Therefore, the boundedness of solutions of (2) are obtained.


### 2.2. Globally Asymptotic Stability of the Equilibrium Points

To begin with, we denote the equilibrium status of model (2) by E∗x∗,y∗. In this subsection, we aim to discuss the existence and global stability of E∗x∗,y∗. By the definition of equilibrium status, E∗x∗,y∗ satisfies the following equalities,(10)r1+fy∗=d1+bx∗+py∗1+ax∗,cpx∗1+ax∗=d2+hy∗.An easy computation yields that the positive equilibrium pointx∗ is(11)x∗=d2+hy∗cp−ad2+hy∗,ifcp−ad2+hy∗>0,and the equilibrium point y∗ satisfies the following equation: (12)pfy∗2+p+fd1+bx∗1+ax∗y∗+1+ax∗bx∗+d1−r=0.By verification, ifr>d1+bx∗, then the solution y∗>0 will always exists.By the monotonicity of the variablef, we derive that the increasing of fear f results in the decreasing of the value of equilibrium point x∗, which leads to a lower value of equilibrium point y∗, while if the conditions cp−ad2+hy∗>0 and r>d1+bx∗ hold, then the equilibrium points E∗x∗,y∗ always exist. In later discuss, we always assume that system (2) has a unique positive equilibrium points E∗x∗,y∗.For readers’s convenience, we denote(13)D1=b−fr21+fy∗−pay∗1+ax∗,D2=1+ax∗ch−fr21+fy∗.Our finding of the long time behavior ofE∗x∗,y∗ of system (2) is as follows:Theorem 1.
Under the conditionDi>0,i=1,2, then for any solutionxt,ytof system (2) with positive initial data, we havelimt⟶∞xt=x∗,limt⟶∞yt=y∗,that is, the equilibrium status of (2) is globally asymptotically stable.Proof.
By the functional differential equation theory, we only need to find a Lyapunov functionalVxt,yt such that D+Vt≤0, where Vxt,yt is a nonnegative function on R+2, D+Vt presents the Dini upper-right derivative. For this purpose, we construct the below two functions,(14)V1x=∫0x−x∗ξdξξ+x∗,V2y=∫0y−x∗ξdξξ+y∗,where x∗,y∗ is the equilibrium status of (2). It is not difficult to verify that both V1x and V2y are nonnegative on R+=x∈R,x>0. If x,y∈R+2, we compute the Dini upper right derivatives of V1x and V2y along system (2), respectively, then(15)D+V1x=x−x∗r1+fyt−τ−d1−bx−py1+ax=x−x∗r1+fyt−τ−d1−bx−py1+ax−r1+fy∗+d1+bx∗+py∗1+ax∗=x−x∗−fryt−τ−y∗1+fyt−τ1+fy∗−bx−x∗−py1+ax∗−py∗1+ax1+ax1+ax∗≤−bx−x∗2+fr1+fy∗x−x∗yt−τ−y∗−p1+ax∗1+ax1+ax∗x−x∗y−y∗+pay∗1+ax∗x−x∗2≤−bx−x∗2+fr21+fy∗x−x∗2+yt−τ−y∗2−p1+ax∗1+ax1+ax∗x−x∗y−y∗+pay∗1+ax∗x−x∗2=−b+fr21+fy∗+pay∗1+ax∗x−x∗2+fr21+fy∗yt−τ−y∗2−p1+ax∗1+ax1+ax∗x−x∗y−y∗,(16)D+V2y=y−y∗−d2−hyt+cpx1+ax=y−y∗−hyt+cpx1+ax+hy∗−cpx∗1+ax∗=−hy−y∗2+cp1+ax1+ax∗x−x∗y−y∗.
In order to determinate the coefficient term ofyt−τ−y∗2, we define a nonnegative function defined on R+ as follows: (17)V3t=∫t−τtyξ−y∗2dξ.
Let(18)Vt=V1x+1+ax∗cV2y+f21+fy∗V3t,then Vt is nonnegative on R+. An easy computation yields(19)D+Vt=−b+fr21+fy∗+pay∗1+ax∗x−x∗2+−1+ax∗ch+fr21+fy∗y−y∗2=−D1x−x∗2−D2y−y∗2.
By the conditions thatD1>0,D2>0, one can easily obtain that, apart from the trail of the equilibrium status, D+Vt<0 always holds, and hence, the differential stability theory (see Reference [18]) implies that E∗x∗,y∗ of system (2) is globally asymptotically stable. The proof is complete.Remark 1.
In Theorem1, by the monotonicity of function Ff=fr/21+fy∗ on the variable f, we find that if the fear f is lower, then the conditions are easier to be satisfied, and if the fear from predators is absent, that is, f=0, then the conditions change to the form of b−pay∗/1+ax∗>0, which means that without fear from predators, the condition is the easiest to be satisfied. In addition, Theorem 1 indicates that the time lag has no influence on the global stability of the equilibrium points of system (2). All these will be elucidated by numerical simulations in Section 4.


## 2.1. Positivity and Boundedness of Solutions

For assurance that model (2) is biologically well-behaved, we want to show that solutions of (2) are positive and bounded. For the first equation of (2), we integrate and get(4)xt=x0exp∫0tr1+fys−τ−d1−bxs−pys1+axsds>0.Applying the same integrating to the second equation of (2) yields(5)yt=y0exp∫0t−d2−hys+c2xs1+axsds>0.Therefore, all solutions of system (2) are positive. Now, we are to discuss the boundedness by use of comparison method. We can obtain from the first equation of (2) that(6)dxtdt≤xtr−d1−bxt.Integrating the above inequality and by use of comparison theory, then we have(7)xt≤r−d1er−d1t1+ber−d1t,which implies(8)lim supt⟶∞xt≤r−d1b.DefineVt=cxt+yt, then(9)V′+μV≤cxr−d1−bx−yd2+hy≤cxr+μ−d1−bx−yd2−μ≤cr+μ−d124b≐ϱ,where μ is a constant such that μ<d2. By the differential inequality theory, we obtain that 0≤V≤ϱ1−eμt/μ+Vx0,y0e−μt. Letting t⟶∞, then 0≤V≤ϱ/μ holds. Therefore, the boundedness of solutions of (2) are obtained.


## 2.2. Globally Asymptotic Stability of the Equilibrium Points

To begin with, we denote the equilibrium status of model (2) by E∗x∗,y∗. In this subsection, we aim to discuss the existence and global stability of E∗x∗,y∗. By the definition of equilibrium status, E∗x∗,y∗ satisfies the following equalities,(10)r1+fy∗=d1+bx∗+py∗1+ax∗,cpx∗1+ax∗=d2+hy∗.An easy computation yields that the positive equilibrium pointx∗ is(11)x∗=d2+hy∗cp−ad2+hy∗,ifcp−ad2+hy∗>0,and the equilibrium point y∗ satisfies the following equation: (12)pfy∗2+p+fd1+bx∗1+ax∗y∗+1+ax∗bx∗+d1−r=0.By verification, ifr>d1+bx∗, then the solution y∗>0 will always exists.By the monotonicity of the variablef, we derive that the increasing of fear f results in the decreasing of the value of equilibrium point x∗, which leads to a lower value of equilibrium point y∗, while if the conditions cp−ad2+hy∗>0 and r>d1+bx∗ hold, then the equilibrium points E∗x∗,y∗ always exist. In later discuss, we always assume that system (2) has a unique positive equilibrium points E∗x∗,y∗.For readers’s convenience, we denote(13)D1=b−fr21+fy∗−pay∗1+ax∗,D2=1+ax∗ch−fr21+fy∗.Our finding of the long time behavior ofE∗x∗,y∗ of system (2) is as follows:Theorem 1.
Under the conditionDi>0,i=1,2, then for any solutionxt,ytof system (2) with positive initial data, we havelimt⟶∞xt=x∗,limt⟶∞yt=y∗,that is, the equilibrium status of (2) is globally asymptotically stable.Proof.
By the functional differential equation theory, we only need to find a Lyapunov functionalVxt,yt such that D+Vt≤0, where Vxt,yt is a nonnegative function on R+2, D+Vt presents the Dini upper-right derivative. For this purpose, we construct the below two functions,(14)V1x=∫0x−x∗ξdξξ+x∗,V2y=∫0y−x∗ξdξξ+y∗,where x∗,y∗ is the equilibrium status of (2). It is not difficult to verify that both V1x and V2y are nonnegative on R+=x∈R,x>0. If x,y∈R+2, we compute the Dini upper right derivatives of V1x and V2y along system (2), respectively, then(15)D+V1x=x−x∗r1+fyt−τ−d1−bx−py1+ax=x−x∗r1+fyt−τ−d1−bx−py1+ax−r1+fy∗+d1+bx∗+py∗1+ax∗=x−x∗−fryt−τ−y∗1+fyt−τ1+fy∗−bx−x∗−py1+ax∗−py∗1+ax1+ax1+ax∗≤−bx−x∗2+fr1+fy∗x−x∗yt−τ−y∗−p1+ax∗1+ax1+ax∗x−x∗y−y∗+pay∗1+ax∗x−x∗2≤−bx−x∗2+fr21+fy∗x−x∗2+yt−τ−y∗2−p1+ax∗1+ax1+ax∗x−x∗y−y∗+pay∗1+ax∗x−x∗2=−b+fr21+fy∗+pay∗1+ax∗x−x∗2+fr21+fy∗yt−τ−y∗2−p1+ax∗1+ax1+ax∗x−x∗y−y∗,(16)D+V2y=y−y∗−d2−hyt+cpx1+ax=y−y∗−hyt+cpx1+ax+hy∗−cpx∗1+ax∗=−hy−y∗2+cp1+ax1+ax∗x−x∗y−y∗.
In order to determinate the coefficient term ofyt−τ−y∗2, we define a nonnegative function defined on R+ as follows: (17)V3t=∫t−τtyξ−y∗2dξ.
Let(18)Vt=V1x+1+ax∗cV2y+f21+fy∗V3t,then Vt is nonnegative on R+. An easy computation yields(19)D+Vt=−b+fr21+fy∗+pay∗1+ax∗x−x∗2+−1+ax∗ch+fr21+fy∗y−y∗2=−D1x−x∗2−D2y−y∗2.
By the conditions thatD1>0,D2>0, one can easily obtain that, apart from the trail of the equilibrium status, D+Vt<0 always holds, and hence, the differential stability theory (see Reference [18]) implies that E∗x∗,y∗ of system (2) is globally asymptotically stable. The proof is complete.Remark 1.
In Theorem1, by the monotonicity of function Ff=fr/21+fy∗ on the variable f, we find that if the fear f is lower, then the conditions are easier to be satisfied, and if the fear from predators is absent, that is, f=0, then the conditions change to the form of b−pay∗/1+ax∗>0, which means that without fear from predators, the condition is the easiest to be satisfied. In addition, Theorem 1 indicates that the time lag has no influence on the global stability of the equilibrium points of system (2). All these will be elucidated by numerical simulations in Section 4.


## 3. Stochastic Model

To get the impact of environmental fluctuation on the dynamics of model (2), we generalize the deterministic model (2) to a stochastic scenario. Since the system randomness mainly brings influence to the natural mortality rates of the interaction species, we structure a stochastic model by perturbing the mortality rate of prey species and predator species respectively [19, 20]. By the known central limit theorem, the error term follows a normal distribution, which sometimes depends on the difference of population scale from the equilibrium status, that is,(20)−d1⟶−d1+σ1xt−x∗tB˙1t,−d2⟶−d2+σ2yt−y∗tB˙2t,where σi2 denotes the magnitude of the environmental noise, Bitt>0 is defined on Ω,ℱ,ℱt≥0,P, which denotes a complete probability space, ℱt≥0 is a filtration satisfying the usual conditions, i=1,2. We assume B1t and B2t are standard and mutually independent Brownian motion. Then we obtain the following stochastic version of (2):(21)dxt=xtr1+fyt−τ−d1−bxt−pyt1+axtdt+σ1xtxt−x∗tdB1t,dyt=yt−d2−hyt−cpxt1+axtdt+σ2ytyt−y∗tdB2t.Supposext be a stochastic process of the below differential equation with n-dimensions,(22)dxt=ft,xtdt+gt,xtdBt.For any Lyapunov functionalVxt, we define the Lyapunov operator L of Vxt as follows:(23)LVxt=Vxxtft,xt+12tracegTt,xtVxxxtgt,xt.Next, we aim to study the existence and uniqueness of positive stochastic process and asymptotic stability of the equilibrium status of stochastic system (21).


### 3.1. Existence and Uniqueness of Solutions

The existence and uniqueness of stochastic solution of (21) is necessary for us to analyze the dynamical behaviors of system (21). In this regard, we have the following finding.Theorem 2.
For system (21) with initial dataϕiθ>0,θ∈−τ,0i=1,2, there exists a positive solutionxt,yt, which is unique and global almost surely (a.s. for abbreviation).Proof.
One can easily find that the coefficients of (21) are all local Lipschitz, then by use of the existence theory of solutions for functional differential equation [21], system (21) has a local solution on −τ,τe which is positive and unique, where τe denotes the explosion time of solutions. As to the proof of global solutions, we only want to show that τe=∞. For this purpose, we take a sufficiently large constant k0 subject to(24)1k0<min−τ≤θ≤0ϕiθ≤max−τ≤θ≤0ϕiθ<k0,i=1,2.
For anyk>k0, we define the stopping time as follows:(25)τk=inft∈−τ,τe:xt∉1k,koryt∉1k,k.
Denoteτ∞=limt⟶∞τk. The definition of τk implies that τ∞≤τe a.s., and hence, it is sufficient for us to prove τ∞=∞ a.s., which means Pτk≤T⟶0 as k⟶∞ for all T>0.
Define(26)Vx,y=x−12lnx−1+y−12lny−1,where xt=x and yt=y for simplicity. Obviously, the function Vx,y is positive when x,y∈R+2. We compute the derivative of Vx,y by applying Itô’s formula along the solution of system (21), then(27)dVx,y=x−12r1+fyt−τ−d1−bxt−py1+axdt+σ122−x81+fyt−τ2x−x∗2dt+x−121+fyt−τx−x∗dB1t+y−12−d2−hyt−cpx1+axdt+σ122−y8y−y∗2dt+y−12y−y∗dB2t≤12r−d1x1/2−b2x3/2−12px1/21+ax+12bx+12py+σ1241+fyt−τ2x2+σ1241+fyt−τ2x∗2−σ1221+fyt−τ2x∗x−σ1281+fyt−τ2x5/2−σ1281+fyt−τ2x∗2x1/2+σ1241+fyt−τ2x∗x3/2+12d2+h2y−d22y1/2−h2y3/2+12y1/2qx+σ224y2++σ224y∗2−σ224y∗y−σ228y3/2−σ228y∗2y1/2+σ224y∗y3/2+x−121+fyt−τx−x∗dB1t+y−12y−y∗dB2t≤−σ1281+fM2x5/2−b2−σ12x∗4x3/2+σ124x2+12bx+12r−d1x1/2+σ12x∗24−σ228y5/2−h2−σ224y∗y3/2+σ224y2+12p+h2−σ224y∗y−d22+σ228y∗2−12cpy1/2+σ224y∗2+12d2+x−121+fyt−τx−x∗dB1t+y−12y−y∗dB2t≤K1+K2dt+x−121+fyt−τx−x∗dB1t+y−12y−y∗dB2t,where(28)K1=supx⟶∞−σ1281+fM2x5/2−b2−σ12x∗4x3/2+σ124x2+12bx+12r−d1x1/2+σ12x∗24,K2=supy⟶∞−σ228y5/2−h2−σ224y∗y3/2+σ224y2+12p+h2−σ224y∗y−d2+σ228y∗2−12cpy1/2+σ224y∗2+12d2.
Integrating the last inequality from 0 toτk∧T yields(29)Vxτk∧T,yτk∧T≤K1+K2τk∧T+∫0τk∧Tx−121+fyt−τx−x∗dB1s+∫0τk∧Ty−12y−y∗dB2s.
Taking expectation in the above inequality leads to(30)EVxτk∧T,yτk∧T≤K1+K2τk∧T.
For everyω∈τk≤T,xτk,ω∉1/k,k and yτk,ω∉1/k,k. Therefore,(31)K1+K2T≥EVxτk∧T,yτk∧T=Pτk≤TVxτk,yτk+Pτk≥TVxτk,yτk≥Pτk≤Tk−12lnk−1∧1k+12lnk−1.
It is not difficult to find thatlimk⟶∞k−1/2lnk−1∧1/k+1/2lnk−1=∞. Then by the upper argumentation, we get limk⟶∞Pτk≤T=0. Due to the arbitrariness of T>0, then Pτ∞<∞=0. That is to say, Pτ∞=∞=1. The proof is complete.


### 3.2. Long-Time Behavior of the Equilibrium Status

Now, we begin to study the long-time behavior of the equilibrium status of (21) in this part.Theorem 3.
We assume that(32)D1>σ12x∗2,D2>σ22y∗2,then for system (21) with initial data ϕiθ>0i=1,2,θ∈−τ,0, all solutions xt,yt of (21) converge to the equilibrium status E∗x∗,y∗ almost surely, that is, limt⟶∞xt=x∗,limt⟶∞yt=y∗,a.s.Proof.
Define two identical functions as given in Theorem 2.2, that is,(33)V1x=∫0x−x∗ξdξξ+x∗,V2y=∫0y−x∗ξdξξ+y∗.
They are nonnegative onR+. By using the Itô’s formula to V1x and V2y, respectively, then(34)LV1x=x−x∗r1+fyt−τ−d1−bx−py1+ax+σ12x∗2x−x∗2≤x−x∗r1+fyt−τ−d1−bx−py1+ax−r1+fy∗+d1+bx∗+py∗1+ax∗+σ12x∗2x−x∗2=x−x∗−fryt−τ−y∗1+fyt−τ1+fy∗−bxt−x∗−py1+ax∗−py∗1+ax1+ax1+ax∗σ12x∗2x−x∗2≤σ12x∗2−bx−x∗2+fr1+fy∗x−x∗yt−τ−y∗−p1+ax∗1+ax1+ax∗x−x∗y−y∗+pay∗1+ax∗x−x∗2≤σ12x∗2−bx−x∗2+fr21+fy∗x−x∗2+yt−τ−y∗2−p1+ax∗1+ax1+ax∗x−x∗y−y∗+pay∗1+ax∗x−x∗2=σ12x∗2−b+fr21+fy∗+pay∗1+ax∗x−x∗2+fr21+fy∗yt−τ−y∗2−p1+ax∗1+ax1+ax∗x−x∗y−y∗,(35)LV2y=y−y∗−d2−hyt+cpx1+ax+σ22y∗2y−y∗2=y−y∗−d2−hyt+cpx1+ax+d2+hy∗−cpx∗1+ax∗+σ22y∗2y−y∗2=σ22y∗2−hy−y∗2+cp1+ax1+ax∗x−x∗y−y∗.
DefineV3t=∫t−τtyξ−y∗2dξ so as to eliminate the term of yt−τ−y∗2. Let(36)Vt=V1x+1+ax∗cV2y+fr21+fy∗V3t.
Easily, we getVt≥0 on t>0. By computation,(37)LVt=σ12x∗2−b+fr21+fy∗+pay∗1+ax∗x−x∗2+σ22y∗2−1+ax∗ch+fr21+fy∗y−y∗2=−D1−σ12x∗2x−x∗2−D2−σ12x∗2y−y∗2.
By the condition thatDi>σi2x∗/2,i=1,2, then we obtain LV<0. It follows from the stochastic differential equation theory [21] that all solutions converge to E∗x∗,y∗, and it is globally asymptotically stable. Hence, the result holds.Remark 2.
The conditions of Theorem3 indicate that the stochastic environmental noises play some crucial influences on the long-time behavior of the positive equilibrium status E∗x∗,y∗ of (21), while under these conditions, time delay τ has no effect on it.


## 3.1. Existence and Uniqueness of Solutions

The existence and uniqueness of stochastic solution of (21) is necessary for us to analyze the dynamical behaviors of system (21). In this regard, we have the following finding.Theorem 2.
For system (21) with initial dataϕiθ>0,θ∈−τ,0i=1,2, there exists a positive solutionxt,yt, which is unique and global almost surely (a.s. for abbreviation).Proof.
One can easily find that the coefficients of (21) are all local Lipschitz, then by use of the existence theory of solutions for functional differential equation [21], system (21) has a local solution on −τ,τe which is positive and unique, where τe denotes the explosion time of solutions. As to the proof of global solutions, we only want to show that τe=∞. For this purpose, we take a sufficiently large constant k0 subject to(24)1k0<min−τ≤θ≤0ϕiθ≤max−τ≤θ≤0ϕiθ<k0,i=1,2.
For anyk>k0, we define the stopping time as follows:(25)τk=inft∈−τ,τe:xt∉1k,koryt∉1k,k.
Denoteτ∞=limt⟶∞τk. The definition of τk implies that τ∞≤τe a.s., and hence, it is sufficient for us to prove τ∞=∞ a.s., which means Pτk≤T⟶0 as k⟶∞ for all T>0.
Define(26)Vx,y=x−12lnx−1+y−12lny−1,where xt=x and yt=y for simplicity. Obviously, the function Vx,y is positive when x,y∈R+2. We compute the derivative of Vx,y by applying Itô’s formula along the solution of system (21), then(27)dVx,y=x−12r1+fyt−τ−d1−bxt−py1+axdt+σ122−x81+fyt−τ2x−x∗2dt+x−121+fyt−τx−x∗dB1t+y−12−d2−hyt−cpx1+axdt+σ122−y8y−y∗2dt+y−12y−y∗dB2t≤12r−d1x1/2−b2x3/2−12px1/21+ax+12bx+12py+σ1241+fyt−τ2x2+σ1241+fyt−τ2x∗2−σ1221+fyt−τ2x∗x−σ1281+fyt−τ2x5/2−σ1281+fyt−τ2x∗2x1/2+σ1241+fyt−τ2x∗x3/2+12d2+h2y−d22y1/2−h2y3/2+12y1/2qx+σ224y2++σ224y∗2−σ224y∗y−σ228y3/2−σ228y∗2y1/2+σ224y∗y3/2+x−121+fyt−τx−x∗dB1t+y−12y−y∗dB2t≤−σ1281+fM2x5/2−b2−σ12x∗4x3/2+σ124x2+12bx+12r−d1x1/2+σ12x∗24−σ228y5/2−h2−σ224y∗y3/2+σ224y2+12p+h2−σ224y∗y−d22+σ228y∗2−12cpy1/2+σ224y∗2+12d2+x−121+fyt−τx−x∗dB1t+y−12y−y∗dB2t≤K1+K2dt+x−121+fyt−τx−x∗dB1t+y−12y−y∗dB2t,where(28)K1=supx⟶∞−σ1281+fM2x5/2−b2−σ12x∗4x3/2+σ124x2+12bx+12r−d1x1/2+σ12x∗24,K2=supy⟶∞−σ228y5/2−h2−σ224y∗y3/2+σ224y2+12p+h2−σ224y∗y−d2+σ228y∗2−12cpy1/2+σ224y∗2+12d2.
Integrating the last inequality from 0 toτk∧T yields(29)Vxτk∧T,yτk∧T≤K1+K2τk∧T+∫0τk∧Tx−121+fyt−τx−x∗dB1s+∫0τk∧Ty−12y−y∗dB2s.
Taking expectation in the above inequality leads to(30)EVxτk∧T,yτk∧T≤K1+K2τk∧T.
For everyω∈τk≤T,xτk,ω∉1/k,k and yτk,ω∉1/k,k. Therefore,(31)K1+K2T≥EVxτk∧T,yτk∧T=Pτk≤TVxτk,yτk+Pτk≥TVxτk,yτk≥Pτk≤Tk−12lnk−1∧1k+12lnk−1.
It is not difficult to find thatlimk⟶∞k−1/2lnk−1∧1/k+1/2lnk−1=∞. Then by the upper argumentation, we get limk⟶∞Pτk≤T=0. Due to the arbitrariness of T>0, then Pτ∞<∞=0. That is to say, Pτ∞=∞=1. The proof is complete.


## 3.2. Long-Time Behavior of the Equilibrium Status

Now, we begin to study the long-time behavior of the equilibrium status of (21) in this part.Theorem 3.
We assume that(32)D1>σ12x∗2,D2>σ22y∗2,then for system (21) with initial data ϕiθ>0i=1,2,θ∈−τ,0, all solutions xt,yt of (21) converge to the equilibrium status E∗x∗,y∗ almost surely, that is, limt⟶∞xt=x∗,limt⟶∞yt=y∗,a.s.Proof.
Define two identical functions as given in Theorem 2.2, that is,(33)V1x=∫0x−x∗ξdξξ+x∗,V2y=∫0y−x∗ξdξξ+y∗.
They are nonnegative onR+. By using the Itô’s formula to V1x and V2y, respectively, then(34)LV1x=x−x∗r1+fyt−τ−d1−bx−py1+ax+σ12x∗2x−x∗2≤x−x∗r1+fyt−τ−d1−bx−py1+ax−r1+fy∗+d1+bx∗+py∗1+ax∗+σ12x∗2x−x∗2=x−x∗−fryt−τ−y∗1+fyt−τ1+fy∗−bxt−x∗−py1+ax∗−py∗1+ax1+ax1+ax∗σ12x∗2x−x∗2≤σ12x∗2−bx−x∗2+fr1+fy∗x−x∗yt−τ−y∗−p1+ax∗1+ax1+ax∗x−x∗y−y∗+pay∗1+ax∗x−x∗2≤σ12x∗2−bx−x∗2+fr21+fy∗x−x∗2+yt−τ−y∗2−p1+ax∗1+ax1+ax∗x−x∗y−y∗+pay∗1+ax∗x−x∗2=σ12x∗2−b+fr21+fy∗+pay∗1+ax∗x−x∗2+fr21+fy∗yt−τ−y∗2−p1+ax∗1+ax1+ax∗x−x∗y−y∗,(35)LV2y=y−y∗−d2−hyt+cpx1+ax+σ22y∗2y−y∗2=y−y∗−d2−hyt+cpx1+ax+d2+hy∗−cpx∗1+ax∗+σ22y∗2y−y∗2=σ22y∗2−hy−y∗2+cp1+ax1+ax∗x−x∗y−y∗.
DefineV3t=∫t−τtyξ−y∗2dξ so as to eliminate the term of yt−τ−y∗2. Let(36)Vt=V1x+1+ax∗cV2y+fr21+fy∗V3t.
Easily, we getVt≥0 on t>0. By computation,(37)LVt=σ12x∗2−b+fr21+fy∗+pay∗1+ax∗x−x∗2+σ22y∗2−1+ax∗ch+fr21+fy∗y−y∗2=−D1−σ12x∗2x−x∗2−D2−σ12x∗2y−y∗2.
By the condition thatDi>σi2x∗/2,i=1,2, then we obtain LV<0. It follows from the stochastic differential equation theory [21] that all solutions converge to E∗x∗,y∗, and it is globally asymptotically stable. Hence, the result holds.Remark 2.
The conditions of Theorem3 indicate that the stochastic environmental noises play some crucial influences on the long-time behavior of the positive equilibrium status E∗x∗,y∗ of (21), while under these conditions, time delay τ has no effect on it.


## 4. Numerical Analysis

In this section, some examples and numerical analysis are presented to verify the main theoretical findings by using Matlab 2014a. For stochastic system (21), we approximate the solution by Milstein method mentioned in Reference [22]. First, we choose the following biologically feasible values of the corresponding parameters,(38)r=0.8;f=0.1;b=0.2;p=0.8;c=0.8;d1=0.2;d2=0.4;h=0.15;a=0.4;τ=0.5;σ1=0.2;σ2=0.2.The initial values are taken asϕ1θ=1,ϕ2θ=0.5 unless other stated. By verification, there exists a coexisting equilibrium status E∗x∗,y∗, where(39)x∗=1.1043,y∗=0.6014.Now, we focus on the long-time behaviors of the equilibrium statuses of (2) and (21) respectively. By an easy computation, we find the criteria of Theorem 1 and Theorem 3 are all satisfied, then the equilibrium status E∗x∗,y∗ of the deterministic system (2) and the stochastic system (21) is globally asymptotically stable, respectively, see Figure 1, which is a numerical testification of the theory results. For low intensity of noise, Figure 1 implies that E∗x∗,y∗ is globally asymptotically stable.Figure 1
The global asymptotical stability ofE∗x∗,y∗, where xt and xi∼t represent the intensities of prey, yit and yi∼t,i=1,2,3 represent the intensities of predator of (2) and (21), respectively. (a) The time series graph of prey and predator of (2). (b) The time series graph of prey and predator of (21).
(a)(b)Now, we are in the position to explore the effect on the long-time behavior ofE∗x∗,y∗. We concentrate on the factors of fear from predator, delayed response to the predation danger, and environmental perturbations numerically.(i)
The effect of fear induced by predator on the equilibrium status. To clarify the effect of fear from predator on the dynamical behaviors of (2) and (21), we choose different values of fii=1,2,3 as f1=0,f2=1 and f3=2, while keeping all remaining parameter values as given above. For different fi, we compute the corresponding equilibrium status. By use of Matlab 2014a, we obtain three pairs of equilibrium status Ei∗xi∗,yi∗,i=1,2,3 as follows:(40)x1∗=1.1414,y1∗=0.6768;x2∗=0.9836,y2∗=0.3450;x3∗=0.9408,y3∗=0.2499.Obviously, different fearfi leads to different equilibrium status Ei∗. By verification, for every value of fear fi, the inequalities in Theorem 1 and Theorem 3 are all gratified respectively, hence, the equilibrium state Ei∗ is globally asymptotically stable, where i=1,2,3. The numerical simulations are shown in Figure 2, where xit and xi∼ti=1,2,3 represent the intensities of prey species, and yit and yi∼ti=1,2,3 represent the intensities of predator species of the deterministic system (2) and the stochastic system (21), respectively. Figure 2 shows that the increasing of the value of fear parameter leads to the decreasing of value of corresponding equilibrium status, but the prey and predator species still maintain global asymptotical stability around the equilibrium status. Numerically, we observe that the fear f induced by predator plays an important role in the equilibrium status.Figure 2
The effect of fear from predator onEi∗xi∗,yi∗ of the deterministic system (2) and the stochastic system (21) with f1=0,f2=1 and f3=2, respectively, where xit and xi∼t,i=1,2,3 represent the intensities of prey species, and yit and yi∼ti=1,2,3 represent the intensities of predator species of (2) and (21), respectively. (a) The time series graph of prey of (2) and (21), respectively. (b) The time series graph of predator of (2) and (21) respectively. (c) The phase graph of prey and predator of (2) with f3=2. (d) The phase graph of prey and predator of (21) with f3=2.
(a)(b)(c)(d)


### 4.1. The Effect of Delay on the Equilibrium Status

In order to illustrate the influence of time delay on the equilibrium points of (2) and (21) through numerical simulations, we keep changing the parameter value of τ and keeping the remaining parameter values as given above. Take τ1=0,τ2=10 and τ3=50, then the equilibrium points are still E∗x∗,y∗. We solve the corresponding equations and plot the time series graph, respectively, see Figure 3, which shows that the time delay has no influences on the equilibrium points E∗x∗,y∗. Whereas if we take a lower initial data ϕ1θ=ϕ2θ=0.1 and a larger time delay as τ=100, then by computer generated simulation, we get the time series graph of the solutions of (2) and (21), see Figures 3(c) and 3(d), respectively. Figures (c) and (d) indicate that the time delay brings no influence on the global stability of the equilibrium points E∗x∗,y∗, but when the time t<100, the solutions are not stable around the equilibrium points. We only plot the part of t<100 of Figures 3(c) and 3(d) and get Figures 3(e) and 3(f) so as to clearly present the influence of time delay. The last two figures of Figure 2 show that the time delay brings important impact on the local stability of the equilibrium points. Comparison of Figures 3(c) and 3(d) implies that the white noises are useful for the solutions staying around the equilibrium points. That is, the white noises are helpful to locally stabilize the equilibrium status.Figure 3
The influence ofτ on the stability of E∗x∗,y∗ of (2) and (21) with τ1=0,τ2=10 and τ3=50, where xit and xi∼t,i=1,2,3 represent the intensities of prey species, and yit and yi∼t,i=1,2,3 represent the intensities of predator species of (2) and (21), respectively. (a) The time series graph of (2) with different τii=1,2,3. (b) The time series graph of (21) with different τii=1,2,3. (c) The time series graph of (2) with ϕ1θ=ϕ2θ=0.1 and τ=100. (d) The time series graph of (21) with ϕ1θ=ϕ2θ=0.1 and τ=100. (e) The local part of (c) as t∈0,100. (f) The local part of (d) as t∈0,100.
(a)(b)(c)(d)(e)(f)


### 4.2. The Effect of Stochastic Parameters on the Equilibrium Status

Analogously, we keep changing the parameter values ofσii=1,2 and keeping other parameter values as fixed. For different white noise with σi=0.2 or σi=1 or σi=5,i=1,2, we solve (21) and get the numerical simulations, see Figure 4(a) and Figure 4(b). It shows that when the level of stochastic parameter value is higher, the equilibrium points E∗x∗,y∗ of (21) still keep globally asymptotically stable, but if the level of stochastic noise is large enough, for example, σ1=20,σ2=26, then the global stability of the equilibrium points E∗x∗,y∗ of (21) is failed, see Figure 4(c) and Figure 4(d). Numerically, the stochastic parameter values have very important influences to the stability of the equilibrium state of (21).Figure 4
The impact ofσi,i=1,2 on the equilibrium state of (2) and (21) with different σi=0.2,1 and σi=5 for i=1,2, where xi∼t and yi∼t,i=1,2,3 represent the intensities of prey and predator species of (21), respectively. (a) The time series graph of prey and predator of (21) with different stochastic parameters. (b) The phase graph of prey and predator of (21) with σ1=σ2=5. (c) The time series graph of prey of (21) with σ1=20,σ2=26. (d) The time series graph of predator of (21) with σ1=20,σ2=26.
(a)(b)(c)(d)


## 4.1. The Effect of Delay on the Equilibrium Status

In order to illustrate the influence of time delay on the equilibrium points of (2) and (21) through numerical simulations, we keep changing the parameter value of τ and keeping the remaining parameter values as given above. Take τ1=0,τ2=10 and τ3=50, then the equilibrium points are still E∗x∗,y∗. We solve the corresponding equations and plot the time series graph, respectively, see Figure 3, which shows that the time delay has no influences on the equilibrium points E∗x∗,y∗. Whereas if we take a lower initial data ϕ1θ=ϕ2θ=0.1 and a larger time delay as τ=100, then by computer generated simulation, we get the time series graph of the solutions of (2) and (21), see Figures 3(c) and 3(d), respectively. Figures (c) and (d) indicate that the time delay brings no influence on the global stability of the equilibrium points E∗x∗,y∗, but when the time t<100, the solutions are not stable around the equilibrium points. We only plot the part of t<100 of Figures 3(c) and 3(d) and get Figures 3(e) and 3(f) so as to clearly present the influence of time delay. The last two figures of Figure 2 show that the time delay brings important impact on the local stability of the equilibrium points. Comparison of Figures 3(c) and 3(d) implies that the white noises are useful for the solutions staying around the equilibrium points. That is, the white noises are helpful to locally stabilize the equilibrium status.Figure 3
The influence ofτ on the stability of E∗x∗,y∗ of (2) and (21) with τ1=0,τ2=10 and τ3=50, where xit and xi∼t,i=1,2,3 represent the intensities of prey species, and yit and yi∼t,i=1,2,3 represent the intensities of predator species of (2) and (21), respectively. (a) The time series graph of (2) with different τii=1,2,3. (b) The time series graph of (21) with different τii=1,2,3. (c) The time series graph of (2) with ϕ1θ=ϕ2θ=0.1 and τ=100. (d) The time series graph of (21) with ϕ1θ=ϕ2θ=0.1 and τ=100. (e) The local part of (c) as t∈0,100. (f) The local part of (d) as t∈0,100.
(a)(b)(c)(d)(e)(f)


## 4.2. The Effect of Stochastic Parameters on the Equilibrium Status

Analogously, we keep changing the parameter values ofσii=1,2 and keeping other parameter values as fixed. For different white noise with σi=0.2 or σi=1 or σi=5,i=1,2, we solve (21) and get the numerical simulations, see Figure 4(a) and Figure 4(b). It shows that when the level of stochastic parameter value is higher, the equilibrium points E∗x∗,y∗ of (21) still keep globally asymptotically stable, but if the level of stochastic noise is large enough, for example, σ1=20,σ2=26, then the global stability of the equilibrium points E∗x∗,y∗ of (21) is failed, see Figure 4(c) and Figure 4(d). Numerically, the stochastic parameter values have very important influences to the stability of the equilibrium state of (21).Figure 4
The impact ofσi,i=1,2 on the equilibrium state of (2) and (21) with different σi=0.2,1 and σi=5 for i=1,2, where xi∼t and yi∼t,i=1,2,3 represent the intensities of prey and predator species of (21), respectively. (a) The time series graph of prey and predator of (21) with different stochastic parameters. (b) The phase graph of prey and predator of (21) with σ1=σ2=5. (c) The time series graph of prey of (21) with σ1=20,σ2=26. (d) The time series graph of predator of (21) with σ1=20,σ2=26.
(a)(b)(c)(d)


## 5. Concluding Remarks

Considering the fear from predator species and a delayed response of prey to the danger outside, we construct a delayed predator–prey mathematical model. Then we spread the model to stochastic environment and get a stochastic predator–prey model. Using Lyapunov functional theory, we investigate the existence and long-time behavior of the coexisting equilibrium status. We find that under some conditions, the equilibria will keep stable in the long term.Then we verify our main results by some numerical examples. Numerically, we show that such factors as fear from predator, the delayed response to the danger outside and stochastic environments bring some significant influences to the long-time behavior of solutions of both system (2) and system (21). Figure 2 implies that the level of fear will change the value of coexisting equilibrium, but all solutions still converge to the equilibrium state from long time aspect. Biologically high level of fear induced from predator leads to low birth speed of prey species, and hence, the equilibrium status of prey is lessened. As a result of low birth speed of prey species, the predator is also affected and the equilibrium state of predator is lessened consequently, which is accordance with the real biological explanation. The delayed fear does not affect the long-time behavior of deterministic system (2) and the stochastic system (21) under some conditions, which is reasonable. For stochastic system, when the stochastic parameter value is not too large, the equilibrium status of (21) keeps globally asymptotically stable. As the increasing of the stochastic noises, the long-time behavior of solutions will be changed from stable to unstable. If the noise is sufficiently large, then equilibrium status will be unstable. In a word, numerical simulations represent that fear from predator and stochastic environment can bring much influence on the values and long-time behaviors of the equilibrium status; however, under some certain conditions, it will not be affected by the delayed fear of prey.Whereas in practice, different functional responses [3,4] often appear and exhibit more different dynamical behaviors. For environmental perturbation, Lévy jump [23] often appears and attracts many authors’ attention. For these scenarios, we believe they are interesting for us to study in our coming work.



---

*Source: 3214358-2022-06-29.xml*