# Dynamics of a Predator-Prey Model with Fear Effect and Time Delay

Authors: Junli Liu; Pan Lv; Bairu Liu; Tailei Zhang
Journal: Complexity
 (2021)

---

## Abstract

In this paper, we propose a time-delayed predator-prey model with Holling-type II functional response, which incorporates the gestation period and the cost of fear into prey reproduction. The dynamical behavior of this system is both analytically and numerically investigated from the viewpoint of stability, permanence, and bifurcation. We found that there are stability switches, and Hopf bifurcations occur when the delayτ passes through a sequence of critical values. The explicit formulae which determine the direction, stability, and other properties of the bifurcating periodic solutions are given by using the normal form theory and center manifold theorem. We perform extensive numerical simulations to explore the impact of some important parameters on the dynamics of the system. Numerical simulations show that high levels of fear have a stabilizing effect while relatively low levels of fear have a destabilizing effect on the predator-prey interactions which lead to limit-cycle oscillations. We also found that the model with or without a delay-dependent factor can have a significantly different dynamics. Thus, ignoring the delay or not including the delay-dependent factor might result in inaccurate modelling predictions.


---

## Body

## 1. Introduction

Predator-prey interaction is a central topic in ecology and evolutionary biology which has been extensively studied from different aspects by many researchers over the last few decades. Mathematical models have played a vital role in understanding these complex scenarios. Most of the existing predator-prey models are based upon classical Lotka–Volterra formalism, and it is usually assumed that the predators can impact the prey population only through direct killing as direct predation is relatively easy to observe in nature. However, the presence of predator may significantly change the behaviour and physiology of prey to such an extent that it could affect the prey population more effective than direct predation [1–3].In the presence of predators, many prey species change their behaviour due to predation risk and show a variety of antipredator responses, which incorporates foraging activity [4], habitat alterations, vigilance [5] and some physiological changes, etc. [3, 6]. Many researchers observed that due to fear of predation risk, the reproduction of the scared prey decreases. For example, in the Greater Yellowstone Ecosystem, wolves influence the reproductive physiology of Elks [7]. Snowshoe hares shift to relatively safe but less profitable microhabitats [8]. Birds react to the sound of predator species with antipredator defenses, and they flee from their nests at the first sight of danger [1]. Other pieces of evidence suggest that fear can also affect populations like mule deer [9] and dugongs [10].In 2011, Zanette et al. conducted a manipulation on song sparrows during a whole breeding season and showed that there is a40% reduction in the offspring production of the song sparrows due to the fear from predator [11]. Very few researchers have investigated the impact of fear effect in a predator-prey model with the help of mathematical modelling. Wang et al. first proposed a predator-prey model by incorporating the cost of fear into prey production. Their analysis showed that fear effect can stabilize the oscillation of the system [12]. Wang and Zou formulated a predator-prey model with the cost of fear and adaptive avoidance of predators. The prey population is divided into juvenile and mature stages; this model is represented by a system of delay differential equations with maturation delay [13]. They theoretically analyzed the model dynamics for constant defense level and adaptive defense level. Numerical simulations showed that both strong adaptation of adult prey and the large cost of fear have destabilizing effect, but large population of predators has a stabilizing effect on the system. Sasmal and Takeuchi studied a predator-prey system with fear effect, group defense, and Holling-type IV functional response; they observed that the model can exhibit multistability [14]. Panday et al. studied the impact of fear in a tritrophic food chain model with Holling-type II functional response and observed that fear can stabilize the system from chaos to stable focus through the period-halving phenomenon, and chaotic dynamics can be controlled by the fear factors [15]. Sarkar and Khajanchi [16] proposed and analyzed a prey-predator system introducing the cost of fear into prey reproduction with Holling-type II functional response. Their mathematical analysis suggests that strong antipredator responses can stabilize the prey-predator interactions by ignoring the existence of periodic behaviors.Inspired by the aforementioned works, the aim of the current work is to develop a predator-prey model, which incorporates the fear effects of predators on prey as well as the time delay. Further, we investigate the impacts of important parameters on the dynamics of system. The rest of the paper is organized as follows. In the next section, we will give the underlying assumptions and formulate the model. In the following section, we present some results on the stability of equilibria and existence of Hopf bifurcation. In Section4, we investigate the properties of the bifurcating periodic solutions by using the normal form theory and center manifold theorem. We show the permanence of the system in Section 5. In Section 6, we provide some numerical simulation results which reveal some potential roles that the fear effect and the time delay may play in predator-prey interactions. A brief discussion completes the paper.


## 2. Model Formulation

Letxt and yt be the densities of prey and predator populations at time t, respectively. We assume that the prey population follows a logistic growth without any predator population, which can be split into three different parts, namely, the birth rate, the natural death rate, and the density dependent death due to intraspecific competition among preys. Let r be the birth rate of prey, d1 be the natural death rate of prey, and d2 be the death rate due to intraprey competition. Since the fear effect will reduce the production, we will multiply the production term of prey rx by a factor fk,y=1/1+ky which accounts for the cost of antipredator defense due to predation fear. Also, it is reasonable to assume that the production of predators after predation is not instantaneous; for the sake of simplicity, we suppose that the time lag required for gestation of predator is τ. Therefore, our predator-prey model with fear effect is given by(1)dxdt=rx1+ky−d1x−d2x2−mxy1+ax,dydt=mce−dτxt−τyt−τ1+axt−τ−dy.Here, the parameter k refers to the level of fear, which is due to the antipredator response of prey. In this study, we adopt the widely used Holling-type II functional response of predators [17], where m/a is the predation rate and 1/a is the half-saturation constant. c is the conversion efficiency from prey biomass to predator biomass and d is the death rate of predator. Since predator goes through a gestation stage, and it endures the death rate, then a factor e−dτ appears in the first term of the second equation of system (1). The parameters r,d1,d2,m,a,c, and d are all positive constants, and k and τ are nonnegative constants.The initial conditions take the form of(2)xθ,yθ=ϕ1θ,ϕ2θ∈C−τ,0,R+2,where C−τ,0,R+2 denotes the space of continuous functions mapping the interval −τ,0 into R+2 and denotes the norm of an element ϕ=ϕ1,ϕ2 in C−τ,0,R+2 by ϕ=sup−τ≤θ≤0ϕ1θ,ϕ2θ.It is straightforward to show (see [18]) the following.Lemma 1.
The solution of system (1) with an initial condition (2) is unique and nonnegative for all t≥0.
Note that from the first equation in (1), we obtain dx/dt≤r−d1x; by a comparison argument, we see that if r<d1, then limt⟶∞xt=0, and applying the theory of asymptotically autonomous systems to the second equation in (1), we also have limt⟶∞yt=0. That is, when r<d1, both prey and predator species go to extinction. So, throughout this paper, we assume that r>d1.Lemma 2.
There exists anM>0 such that for any solution of system (1) with an initial condition (2), there must be a T>0 such that xt≤M and yt≤M for t≥T.Proof.
We first consider the system(3)dzdt=r−d1z−d2z2,whose solution satisfies(4)lim supt⟶∞zt≤r−d1d2.
Sincedx/dt≤r−d1x−d2x2, then by the comparison principle, we have(5)lim supt⟶∞xt≤r−d1d2.
Hence, for anyε>0 small enough, there exists a t∗>0 sufficiently large such that(6)xt≤r−d1d2+εfort>t∗.
ForLt=ce−dτxt−τ+yt and t>t∗+τ, we have(7)dLtdt=ce−dτrxt−τ1+kyt−τ−d1xt−τ−d2x2t−τ−dyt,≤ce−dτr−d1xt−τ−dyt=ce−dτr−d1+dxt−τ−dLt≤ce−dτr−d1+dr−d1d2+ε−dLt.
It follows that(8)lim supt⟶∞Lt≤ce−dτr−d1+dr−d1/d2+εd.
This completes the proof. □Lemma 3 (see [19], Lemma 2.1).
Consider the following equation:(9)x˙t=axt−τ−bxt,wherea,b,τ>0;xt>0for−τ≤t≤0. We have(i)
ifa<b, then limt⟶+∞xt=0(ii)
ifa>b, then limt⟶+∞xt=+∞
It is easy to check that (1) has a trivial extinction equilibrium E0=0,0. This system also has a semitrivial equilibrium E1=r−d1/d2,0 under the condition r>d1. To find the positive (coexistence) equilibrium E∗=x∗,y∗, letting the right-hand side be zero, we then get x∗=d/mce−dτ−da, which exists if mce−dτ>da, and y∗ is a positive root of the equation(10)a0y2+a1y+a2=0,where a0=mk,a1=m+kd1+d2x∗1+ax∗,a2=d1+d2x∗−r1+ax∗. Clearly, (10) has one positive root y∗ if and only if d1+d2x∗−r<0; this inequality is equivalent to(11)mce−dτ>da+dd2r−d1.Then, we conclude that (1) has a unique positive equilibrium E∗=x∗,y∗ if (11) holds, and the expression for y∗ is as follows:(12)y∗=−a1+a12−4a0a22a0.


## 3. Stability Analysis and Hopf Bifurcation

Theorem 1.
The trivial extinction equilibriumE0=0,0is always an unstable saddle point. The semitrivial equilibriumE1=r−d1/d2,0is globally asymptotically stable inC˜=ϕ1θ,ϕ2θ∈C−τ,0,R+2|ϕ10>0if (11) is reversed and is unstable if (11) holds; meanwhile,E∗exists.Proof.
For equilibriumE0=0,0, the characteristic equation is(13)λ+d1−rλ+d=0.
Obviously, (13) has one positive root r−d1 and one negative root −m. Hence, E0 is a saddle point.
The characteristic equation atE1=r−d1/d2,0 is as follows:(14)λ+r−d1λ+d−mce−dτr−d1d2+ar−d1e−λτ=0.
Letfλ=λ+d−mce−dτr−d1/d2+ar−d1e−λτ; then it is clear that (14) has one root d1−r<0 and the other roots are determined by the equation fλ=0.
If (11) is reversed, then d>mce−dτr−d1/d2+ar−d1. Let λ=u+iv with u,v∈ℝ be a root of fλ=0. Then, we have(15)u+d2+v2=mce−dτr−d1e−uτd2+ar−d12.
Ifu≥0, then(16)u+d2+v2>mce−dτr−d1d2+ar−d12≥mce−dτr−d1e−uτd2+ar−d12.
A Contradiction. This shows that all roots of fλ=0 must have negative real parts. Therefore, E1 is locally asymptotically stable.
Next, we will show thatE1 is globally attractive if (11) is reversed. Since mce−dτ<da+dd2/r−d1, we can choose small enough ε>0 such that(17)mce−dτr−d1+d2εd2+ar−d1+d2ε<d.
From (6), we know that(18)dydt≤mce−dτr−d1/d2+ε1+ar−d1/d2+εyt−τ−dyt=mce−dτr−d1+d2εd2+ar−d1+d2εyt−τ−dyt.for t>t∗. But for the equation(19)dVdt=mce−dτr−d1+d2εd2+ar−d1+d2εVt−τ−dVt,by (17) and Lemma 3, it follows that limt⟶∞Vt=0. Therefore, an application of the standard comparison argument yields limt⟶∞yt=0. Then, the limiting equation for xt is(20)dxdt=rx−d1x−d2x2,x0>0.
Therefore,limt⟶∞xt=r−d1/d2, and E1 is globally asymptotically stable in C˜.
Note also that if (11) holds, then f0=d−mce−dτr−d1/d2+ar−d1<0, and f+∞=∞. Hence, fλ=0 has at least one positive root and E1 is unstable.Theorem1 implies that E1 loses its stability to the occurrence of the positive equilibrium E∗ when (11) holds. The local stability of E∗ is addressed in the following theorem.When (11) holds, linearizing system (1) at the positive equilibrium E∗=x∗,y∗, we get the characteristic equation as follows:(21)λ2+m1λ+m2+e−λτm3λ+m4=0,where(22)m1=d+d2x∗−max∗y∗1+ax∗2,m2=dd2x∗−max∗y∗1+ax∗2,m3=−d,m4=d−d2x∗+max∗y∗1+ax∗2+mce−dτy∗1+ax∗2rkx∗1+ky∗2+mx∗1+ax∗.Whenτ=0, (21) yields(23)λ2+m1+m3λ+m2+m4=0.It is easy to see thatm2+m4>0 and m1+m3=m2/d. To ensure the local stability of E∗, we need m2>0.Letλ=iω, ω>0 is a root of (21); then, (21) implies(24)ω2−m2=m3ωsinωτ+m4cosωτ,m1ω=m4sinωτ−m3ωcosωτ.Squaring and adding both equations of (24), it follows that(25)z2+m12−m32−2m2z+m22−m42=0,where z=ω2. A straightforward computation shows that(26)m12−m32−2m2=m22d2.Now, we can summarize the above result on the local stability of the positive equilibriumE∗.Theorem 2.
Assume that (11) holds; thenE∗=x∗,y∗exists. Moreover, (i) ifm2>0andm2>m4, then (21) has no purely imaginary root, andE∗is locally asymptotically stable for allτ≥0. (ii) Ifm4>m2>0, then (21) has a unique pair of purely imaginary roots±iω0, andE∗may lose stability asτincreases from zero, whereω0=z0andz0=−m12−m32−2m2/2+m12−m32−2m22−4m22−m42/2is the unique positive root of (25).To further investigate the stability ofE∗ under condition (ii) of Theorem 2, we need to check the transversality condition at λ=±iω0. Since x∗ and y∗ are both time delay-dependent, it is hard to discuss the transversality condition, if we suppose that the coefficient in conversing prey into predator ce−dτ=c0 is a constant, independent of delay τ, then we can extend the result in Theorem 2 as follows.Theorem 3.
Ifce−dτ=c0is a constant,mc0>da+dd2/r−d1andm4>m2>0, then there exists a sequence of values ofτ:(27)τ0<τ1<⋯<τn<⋯such thatE∗=x∗,y∗is locally asymptotically stable whenτ∈0,τ0and unstable whenτ>τ0. Furthermore, system (1) undergoes a Hopf bifurcation atE∗whenτ=τn,n=0,1,2,…, whereτnis defined by (28).Proof.
From (24), we know that τn corresponding to ω0 is(28)τn=1ω0arccosm4ω02−m2−m1m3ω02m32ω02+m42+2nπω0,n=0,1,2,….
To complete the proof of Theorem3, we only need to show that(29)dReλdτ|τ=τn>0.
Differentiating (21) with respect to τ, we have(30)2λ+m1dλdτ−τe−λτm3λ+m4dλdτ−λe−λτm3λ+m4+m3e−λτdλdτ=0.
It follows that(31)dλdτ−1=2λ+m1−τe−λτm3λ+m4+m3e−λτλe−λτm3λ+m4=2λ+m1λe−λτm3λ+m4+m3λm3λ+m4−τλ=λ2−m2−λ2λ2+m1λ+m2−m4λ2m3λ+m4−τλ,and hence,(32)signdReλdτ|τ=τn=signRedλdτ−1|τ=τn=signReλ2−m2−λ2λ2+m1λ+m2|τ=τn−Rem4λ2m3λ+m4|τ=τn=signReλ2−m2−λ2λ2+m1λ+m2|τ=τn−Rem4λ2m3λ+m4|τ=τn=signω04−m22ω02m2−ω022+m12ω02+m42ω02m42+m32ω02=signω04−m22+m42ω02m42+m32ω02.
Then, we havedReλ/dτ|τ=τn>0. Therefore, the transversality condition holds, and the Hopf bifurcation occurs at τ=τn.


## 4. Properties of Hopf Bifurcations

In this section, we will study the properties of the Hopf bifurcations by using the normal theory and the center manifold theorem due to Hassard et al. [20].Letx1=x−x∗,x2=y−y∗,τ=τn+v,x¯it=xiτt and dropping the bars for simplification of notations, system (1) is transformed into a functional differential equation in C=C−1,0,ℝ2 as(33)x˙t=Lvxt+Fv,xt,where xt=x1t,x2tT∈ℝ2, and Lv:C⟶ℝ2,F:ℝ×C⟶ℝ2 are given, respectively, by(34)Lvϕ=τn+va11a120a22ϕ10ϕ20+τn+v00b21b22ϕ1−1ϕ2−1,and(35)Fv,ϕ=τn+vF1,F2T,where(36)a11=r1+ky∗−d1−2d2x∗−my∗1+ax∗2,a12=−krx∗1+ky∗2−mx∗1+ax∗,a22=−d,b21=c0my∗1+ax∗2,b22=c0mx∗1+ax∗,F1=K11ϕ120+K12ϕ10ϕ20+K13ϕ220+K14ϕ130+K15ϕ120ϕ20+K16ϕ10ϕ220+K17ϕ230+⋯,F2=K21ϕ12−1+K22ϕ1−1ϕ2−1+K23ϕ13−1+K24ϕ12−1ϕ2−1+⋯,K11=may∗1+ax∗3−d2,K12=−kr1+ky∗2−m1+ax∗2,K13=k2rx∗1+ky∗3,K14=−ma2y∗1+ax∗4,K15=ma1+ax∗3,K16=k2r1+ky∗3,K17=−k3rx∗1+ky∗4,K21=−c0may∗1+ax∗3,K22=c0m1+ax∗2,K23=c0ma2y∗1+ax∗4,K24=−c0ma1+ax∗3.By the Riesz representation theorem, there exists a functionηθ,v of bounded variation for θ∈−1,0, such that(37)Lvϕ=∫−10dηθ,vϕθ,forϕ∈C.In fact, we can choose(38)ηθ,v=τn+va11a120a22δθ−τn+v00b21b22δθ+1,where δ denotes the Dirac delta function. For ϕ∈C1−1,0,ℝ2, definex(39)Avϕ=dϕθdθ,θ∈−1,0,∫−10dηs,vϕs,θ=0,Rvϕ=0,θ∈−1,0,Fv,ϕ,θ=0.Then, system (1) is equivalent to(40)x˙t=Avxt+Rvxt,where xtθ=xt+θ for θ∈−1,0. For ψ∈C10,1,ℝ2∗, define(41)A∗ψs=−dψsds,s∈0,1,∫−10dηTt,0ψ−t,s=0,and a bilinear inner product(42)ψs,ϕθ=ψ¯0ϕ0−∫−10∫ξ=0θψ¯ξ−θdηθϕξdξ,where ηθ=ηθ,0. Then, A0 and A∗ are adjoint operators. By the discussion in Section 3, we know that ±iω0τn are eigenvalues of A0. Hence, they are also eigenvalues of A∗. We first need to compute the eigenvectors of A0 and A∗ corresponding to iω0τn and −iω0τn, respectively.Supposeqθ=1,αTeiω0τnθ is the eigenvector of A0 corresponding to iω0τn; then A0qθ=iω0τnqθ. Then, from the definition of A0, (34), (37), and (38), we have(43)α=iω0−a11a12.Similarly, we can obtain the eigenvectorq∗s=D1,α∗eiω0τns of A∗ corresponding to −iω0τn, where(44)α∗=−iω0+a11b21eiω0τn.ChoosingD as D=1+α¯α∗+α∗τneiω0τnb21+α¯b22−1, then by (42), we see q∗s,qθ=1.Now, we will use the algorithms given in [20] and using a computation process similar to that in [21–23], we get the coefficients used in determining the qualities of bifurcating periodic solutions:(45)g20=2τnD¯K11+αK12+α2K13+α¯∗K21+αK22e−2iω0τn,g11=2τnD¯K11+K12Reα+αα¯K13+α¯∗K21+K22Reα,g02=2τnD¯K11+α¯K12+α¯2K13+α¯∗K21+α¯K22e2iω0τn,g21=2τnD¯W2010+2W1110K11+W1120+12W2020+α¯2W2010K12+2αW1120+α¯W2020K13+3K14+2α+α¯K15+α2α¯+αK16+3α2α¯K17+α¯∗2W111−1e−iω0τn+W201−1eiω0τnK21+W112−1+αW111−1e−iω0τn+12W202−1+α¯W201−1eiω0τnK22+3K23+α¯+2αK24e−iω0τn,where(46)W20θ=ig20ω0τnq0eiω0τnθ+ig¯023ω0τnq¯0e−iω0τnθ+E1e2iω0τnθ,W11θ=−ig11ω0τnq0eiω0τnθ+ig¯11ω0τnq¯0e−iω0τnθ+E2,E11=M11M1,E12=M12M1,E21=M21M2,E22=M22M2,M1=2iω0−a11−a12−b21e−2iω0τn2iω0−a22−b22e−2iω0τn,M11=2K11+αK12+α2K13−a12K21+αK222iω0−a22−b22e−2iω0τn,M12=22iω0−a11K11+αK12+α2K13−b21e−2iω0τnK21+αK22,M2=−a11−a12−b21−a22−b22,M21=2K11+K12Reα+αα¯K13−a12K21+K22Reα−a22−b22,M22=2−a11K11+K12Reα+αα¯K13−b21K21+K22Reα.Thus, we can determineW20θ and W11θ. Furthermore, we can compute g21 by (45). Then, we can compute the following values:(47)c10=i2ω0τng20g11−2g112−g0223+g212,μ2=−Rec10Reλ′τn,T2=−Imc10+μ2Imλ′τnω0τn,β2=2Rec10.The properties of bifurcating periodic solution in the center manifold at the critical valueτn are determined by the values in (47). From [20], we know that μ2 determines the directions of the Hopf bifurcation: if μ2>0μ2<0, then the Hopf bifurcation is supercritical (subcritical); β2 determines the stability of the bifurcating periodic solutions: the bifurcating periodic solutions are stable (unstable) if β2<0β2>0; and T2 determines the period of the bifurcating periodic solutions: the period increase (decrease) if T2>0T2<0.


## 5. Permanence

In order to prove the uniform persistence of system (1), we need the uniform persistence theorem for infinite dimensional systems from [24]. Let X be a complete metric space. Suppose that X0 is open and dense in X, with X0∪X0=X, X0∩X0=∅. Assume that Tt is a C0−semigroup on X satisfying(48)Tt:X0⟶X0,Tt:X0⟶X0.LetTbt=Tt|X0 and let Ab be the global attractor for Tbt. Assume further that(i)
There is at0≥0 such that Tt is compact for t>t0;(ii)
Tt is point dissipative in X;(iii)
A˜b=∪x∈Abωx is isolated and has an acyclic covering M, where M=∪i=1kMk;(iv)
WsMi∩X0=∅ for i=1,2,…,k.Then,X0 is a uniform repeller with respect to X0, i.e., there is an ϵ>0 such that for any x∈X0,liminft⟶+∞dTtx,X0≥ε, where d is the distance of Ttx from X0.Now, we consider the persistence of system (1).Theorem 4.
Let (11) hold; then there exists anε>0such that every solutionxt,ytof system (1) with initial conditionsϕ1θ≥0,ϕ2θ≥0andϕ1θ0ϕ2θ0≠0for someθ0∈−τ,0satisfies(49)liminft⟶∞xt≥ε,liminft⟶∞yt≥ε.Proof.
Let(50)X=C−τ,0,R+2,X0=ϕ1,ϕ2∈X|ϕ1θ∗ϕ2θ∗>0for someθ∗∈−τ,0,X0=X∖X0=ϕ1,ϕ2∈X|ϕ1θϕ2θ≡0for allθ∈−τ,0.
It suffices to show that there exists anε>0 such that for any solution ut of system (1) initiating from X0,lim inft⟶∞dut,X0≥ε. To this end, we verify below that the conditions of the above uniform persistence theorem are satisfied. Firstly, we show X0 is positively invariant. By (1), we have(51)xt=xθ∗exp∫θ∗tr1+kyρ−d1−d2xρ−myρ1+axρdρ,and(52)yt=y0+∫0tmxρ−τyρ−τe−dτ−ρ1+axρ−τdρe−dt.
Then,xt>0 and yt>0 for all t>0ifΦ=ϕ1,ϕ2∈X0. This implies that X0 is positively invariant.
Similarly, we can show thatX0 is positively invariant, and condition (48) is satisfied. We have verified the point dissipativeness of the semiflow of system (1) in 2.
Denote theω−limit set of the solution of system (1) starting in ϕ1,ϕ2∈X by ωϕ1,ϕ2. Let(53)Ω=∪ωϕ1,ϕ2|ϕ1,ϕ2∈X0.
For any givenϕ1,ϕ2∈X0, we have ϕ1θϕ2θ≡0 for all θ∈−τ,0. Ifϕ1θ˜=0 for some θ˜∈−τ,0, then xt≡0 for all t≥0, and yt⟶0 as t⟶∞. Ifϕ2θ^=0 for some θ^∈−τ,0, we have the following two cases.Case 1.
y0=0. Then from (52), we have yt≡0 for all t≥0, and x˙=rx−d1x−d2x2, it follows that xt⟶0orxt⟶r−d1/d2 as t⟶∞.Case 2.
y0>0. Then we must have x0=0; from the first equation of (1), we have xt≡0 for all t≥0; again from the second equation of (1), we have yt⟶0 as t⟶∞.
Finally, we see that for any givenϕ1,ϕ2∈X0,xt,yt⟶E0 as t⟶∞, or xt,yt⟶E1 as t⟶∞. Thus, Ω=E0,E1. If invariant sets E0 and E1 are isolated, then Ω is isolated and is an acyclic covering. The isolated invariance of E0 and E1 will follow from the following proof.
We now show that(54)WsE0∩X0=∅,WsE1∩X0=∅,where WsE0(WsE1) denotes the stable manifold of E0(E1). Suppose WsE0∩X0≠∅; then there exists a solution xt,yt∈X0,t≥0, of (1) with initial conditions in X0, such that(55)limt⟶∞xt=0,limt⟶∞yt=0.
Sincer>d1, we can choose sufficiently small η1>0 so that(56)r1+kη1−d1−mη1>0.
Lett0>0 sufficiently large such that(57)0<yt<η1,t>t0.
Then, fort>t0, we have(58)dxdt≥r1+kη1x−d1x−d2x2−mη1x,xt0>0.
Consider the following equation:(59)dWdt=r1+kη1W−d1W−d2W2−mη1W,Wt0=xt0>0.
By the comparison principle, we havext≥Wt for all t≥t0. It is easy to see that (59) has a positive equilibrium W∗=r/1+kη1−d1−mη1/d2, which is globally attractive. Note that xt≥Wt for all t≥t0 and limt⟶∞Wt=W∗; this is a contradiction to (55).
SupposeWsE1∩X0≠∅; then there exists a solution xt,yt∈X0,t≥0, of (1) with initial conditions in X0, such that(60)limt⟶∞xt=r−d1d2,limt⟶∞yt=0.
Then, for anyr−d1/d2>η2>0, by (60), there exists a t1>0 such that(61)r−d1d2−η2<xt<r−d1d2+η2,for allt>t1−τ.
Hence, by (1),(62)dytdt≥mce−dτd2−η21+ar−d1/d2−η2yt−τ−dyt,=mce−dτr−d1−d2η2d2+ar−d1−d2η2yt−τ−dyt,t>t1.
Since (11) holds, we can choose η2>0 sufficiently small so that mce−dτr−d1−d2η2/d2+ar−d1−d2η2>d. Then, by Lemma 3, the solution of the linear delay differential equation u˙t=mce−dτr−d1−d2η2/d2+ar−d1−d2η2ut−τ−dut with positive initial condition must converge to infinity as t⟶∞. Therefore, by the comparison principle, we have yt⟶∞ as t⟶∞, which is a contradiction to (60). Thus, (54) holds. This completes the proof.In Lemma2, we have shown that the solutions of system (1) are ultimately bounded; together with Theorem 4, we now give the main result in this section.Theorem 5.
System1is permanent provided11holds and the initial conditions are shown inTheorem 4.


## 6. Numerical Simulations

In order to investigate the role of time delay and fear in our model (1), we perform some numerical simulations to illustrate the analytical results observed in the previous sections. All the simulations are carried out in MATLAB.First, we choose the parameter values as(63)r=0.1,k=0.3,d1=0.01,d2=0.01,m=0.5,c=0.4,d=0.22,a=0.6,τ=0.1,where most of the parameters are taken from [12].For the above set of parameter values, in Figure1, we observe that the system shows a limit cycle around the interior equilibrium E∗. If we further increase the value of τ, for τ=0.4, we see that the solutions approach the coexistence equilibrium; see Figure 2.Figure 1
The figure shows that the dynamics of system (1) of prey and predator populations. (a) represents the stable and periodic oscillations around the interior equilibrium E∗ and the phase plane diagram. (b) represents the stable limit cycle around E∗ when τ=0.1. The parameter values are given in (63).
(a)(b)Figure 2
The figure shows that the dynamics of the model system (1) of prey and predator populations. (a) represents the coexistence equilibrium E∗ that is locally asymptotically stable and (b) represents the phase-portrait of E∗ that is locally asymptotically stable when τ=0.4. The other parameter values are the same as in Figure 1.
(a)(b)To observe the long-term behavior for a range of values ofτ, in Figure 3, we draw the bifurcation diagram with respect to τ. From Figure 3, we can see that as we increase the value of τ, the prey and predator populations show oscillatory coexistence for the small values of τ; then, both populations show stable coexistence. This indicates that increasing delay τ may induce a transition from the state where the populations of the prey and predator oscillate periodically to a stability situation. And also, we note that the prey population increases, while the predator population decreases as we gradually increase τ, and finally the predator population goes to extinction; the system becomes stable around the semitrivial equilibrium E1.Figure 3
The figure shows the bifurcation diagrams of system (1) of prey and predator populations with respect to the parameter τ, which indicates a transition from the state where the populations of the prey and predator oscillate periodically to a stability situation. The other parameters are given in (63).
(a)(b)To investigate the impact of fear when the system is unstable around the positive equilibrium, we fixτ=0.2. We draw the bifurcation diagram considering k as a bifurcation parameter. Figure 4 shows that as we increase the level of the fear effect in prey, the coexistence of prey-predator changes from stable periodic oscillations to stable equilibrium. Figure 4 indicates that the prey population tends to a steady state. The fear effect will not affect the prey population over the long term. Biologically, after a certain level of fear with prey population, the fear has no effect due to physiological impact when the prey populations are habituated with fear from predator species, while the predator population eventually will decrease with increasing k.Figure 4
The figure shows the bifurcation diagrams of system (1) of prey and predator populations with respect to the parameter k, which indicates that high levels of fear can stabilize the predator-prey system by excluding the existence of periodic solutions, where τ=0.2 and other parameters are given in (63).
(a)(b)If we choose parameter values as those in Figure1 and let m vary in the interval 0.3,0.65, Figure 5 gives the bifurcation diagram. From this figure, we can see that the predator population can survive only for the large values of m. If the value of the predation rate is low, then the predator population cannot survive. As we increase the value of m, the coexistence of prey-predator changes at stable equilibrium to stable oscillatory coexistence.Figure 5
(a) and (b) are the bifurcation diagrams of prey and predators with respect to parameterm, respectively, which indicate the emergence of Hopf bifurcations as k increases. The parameters are given in (63). In (c), τ=0.1, m=0.45, and the positive equilibrium is asymptotically stable. In (d), τ=0.1, m=0.6, there occurs a limit cycle.
(a)(b)(c)(d)To investigate the properties of the bifurcating periodic solution, we neglect the delay-dependent factore−dτ; that is, we consider the following system:(64)dxdt=rx1+ky−d1x−d2x2−mxy1+ax,dydt=mc0xt−τyt−τ1+axt−τ−dy.We use the following parameter values:(65)r=0.4,k=0.3,d1=0.01,d2=0.04,m=0.5,c0=0.4,d=0.22,a=0.01.Then, (64) has a positive equilibrium E∗=1.1122,0.579 and satisfies the conditions indicated in Theorem 3; then, E∗=1.1122,0.579 is asymptotically stable when τ=0. Take n=0; we then have ω0=0.2504, τ0=0.7003, μ2=0.1746>0, β2=−0.0081<0, and T2=0.2383>0. Thus, E∗ is stable for τ<τ0 as illustrated in Figure 6. When τ passes through the critical value τ0, E∗ will lose its stability and a Hopf bifurcation will occurs; that is, a family of periodic solutions bifurcate from E∗. Since μ2>0, β2<0, then the Hopf bifurcation is supercritical and its direction is τ>τ0. These bifurcating periodic solutions from E∗ at τ0 are stable; see Figure 7.Figure 6
The figure shows that the dynamics of the system (64) of prey and predators. (a) represents the positive equilibrium E∗=1.1122,0.579 that is locally asymptotically stable and (b) represents the phase-portrait of E∗ that is locally asymptotically stable when τ=0.6<τ0=0.7003. The parameters are given in (65).
(a)(b)Figure 7
(a) represents the stable and periodic oscillations around the positive equilibriumE∗=1.1122,0.579 and the phase plane diagram (b) represents the stable limit cycle around E∗ when τ=0.75>τ0=0.7003. The parameters are given in (65).
(a)(b)In Figure8, we explore the effects of delay τ on the dynamics of system (64). Figure 8 shows that as τ increases, at first the positive equilibrium is stable; then when τ is greater than a certain critical value, periodic solution occurs; this is coincident with Theorem 3, which implies that delay destabilizes the equilibrium and induces oscillatory dynamics. By comparing Figure 3 and Figure 8, we see that the difference is big in the dynamics when the delay-dependent factor is incorporated or not.Figure 8
The figure shows the bifurcation diagrams of system (64) of prey and predator populations, which indicates the emergence of Hopf bifurcations as parameter τ increases, and other parameters are given in (65).
(a)(b)In Figure9, taking k as a bifurcation parameter, we found that the fear effects can stabilize the system. In Figure 10, we draw the bifurcation diagram with respect to m; the shapes of the bifurcation diagrams for system (64) in Figure 10 and the original system (1) in Figure 5 are very similar.Figure 9
The figure shows the bifurcation diagrams of system (64) of prey and predator populations with respect to the parameter k, which shows that high levels of fear can stabilize the predator-prey system. However, relatively low levels of fear can induce oscillations via Hopf bifurcations, where τ=0.8 and other parameters are given in (65).
(a)(b)Figure 10
This figure shows the bifurcation diagrams of system (64) of prey and predators with respect to parameter m, which indicates the emergence of Hopf bifurcations as k increases, where τ=0.8 and other parameters are given in (65).
(a)(b)


## 7. Discussion

In this paper, taking into account the time delay corresponding to the gestation period and the effect that the fear of predators has on prey, we have proposed a predator-prey model with Holling-type II functional response. We have shown that the boundary equilibriumE1 is globally asymptotically stable if (11) is reversed and pointed out that the instability of E1 implies the permanence of the system and existence of the positive equilibrium E∗. Furthermore, for the system without delay-dependent factor, by choosing the delay τ as a bifurcation parameter, sufficient conditions for the existence of Hopf bifurcation have been presented. By using the theory of normal form and center manifold theorem, the explicit formulae determining the direction, stability, and other properties of bifurcating periodic solutions have been derived.The novelty of this paper is that we have incorporated both the fear effect and a delay-dependent factor that ensures the predators that do not survive the delay interval do not contribute to the growth of the predator population. Our investigation shows that the dynamics of system (1) are affected significantly by not only the biomass transfer delay, but also the level of the fear effect. In [25], the delay-dependent factor was neglected.Numerical simulations are performed to show the roles that time delay and fear effects can play in predator-prey interactions. For the original system (1), we have shown that both large time delay and high levels of fear (or equivalently strong antipredator responses) can stabilize the predator-prey system by excluding the existence of periodic solutions. And also, we observe that if the time delay as well as the fear level increases, the predator population goes to extinction. While for system (64) which neglects the delay-dependent factor, numerical simulations show that the delay has destabilizing effect while the large cost of fear has a stabilizing effect on the predator-prey interactions. Figure 8 indicates that ignoring the delay-dependent factor results in model in (63) that does not predict the extinction of the predators, even if the delay is very large. Therefore, time delay and fear effect significantly influence the dynamics of the system and play important roles in ecosystems stability.In the current paper, we assumed that the fear effects only reduce the birth rate of the prey, but it may also increase the death rate of adult prey due to long-term physiological impacts [11]. To make the model more realistic, some other important mechanisms can also be incorporated; for example, spatial effect is an important factor due to foraging behaviors of both prey and predator. The functional response plays a vital role in determining the different complex dynamical behaviors; different functional response functions can be chosen to reflect the group defense. It would be interesting and challenging to consider the aforementioned possible extensions; we leave these problems for further investigation.



---

*Source: 9184193-2021-04-28.xml*