# Population Dynamic Study of Prey-Predator Interactions with Weak Allee Effect, Fear Effect, and Delay

Authors: Ye Xuan Li; Hua Liu; Yu Mei Wei; Ming Ma; Gang Ma; Jing Yan Ma
Journal: Journal of Mathematics
 (2022)

---

## Abstract

In this study, a predator-prey model with the Allee effect and fear effect is established. We use the comparison principle to prove boundedness. The zero equilibrium point and nonzero equilibrium point of the model are calculated, and the local stability conditions are obtained. Next, according to the Sotomayor theorem, the cross-sectional conditions of transcritical bifurcation and Hopf bifurcation are obtained. The conditions for the Hopf bifurcation to be supercritical or subcritical can be calculated by the normal form theory. Then, to make the model more realistic, we introduce the gestation delay in the proposed mathematical model. Stability and Hopf bifurcation are also analyzed. Finally, several numerical simulations are presented to verify the conclusions. Our results demonstrate that the Allee effect, fear effect, and delay play significant roles in population dynamics. The Allee effect and delay destabilize the originally stable model, after which Hopf bifurcation occurs. However, the fear effect can enhance stable coexistence.


---

## Body

## 1. Introduction

In mathematical ecology, the interaction between the predator and prey is important and complex and will continue to be predominant in ecology because it is widespread globally. Since the establishment of the well-known Lotka–Volterra model solved the key problems of ecological processes [1, 2], applied mathematicians usually use mathematical models as tools to study the complex dynamics between predators and prey. It is an effective method to study the dynamic behavior of the population by establishing realistic mathematical models. Different types of models are proposed to investigate the interactions in mathematical biology. The classic predator-prey model can be expressed by the following first-order ordinary differential equations:(1)dxdt=xlx−gx,yy,dydt=−δy+hx,yy,where x and y are the population of prey and predator, respectively, lx is the birth rate of prey, and δ is the predator death rate. Here, the role of the functional response term is to connect the predator and the prey, that is, the prey biomass consumed by each predator in unit time, which is described by gx,y, while hx,y is called the numerical response of the predator [3], which means the growth rate of the predator realized by consuming the prey. Holling proposed several different functional response types. Because this method can reflect the interspecies relationship between the predator and prey in the model, it has wide applicability. Many scholars incorporate this method into mathematical ecology modeling. Holling type-I, type-II, and type-III functional responses have been extensively studied [3–5]. Usually, gx,y is called a prey-dependent functional response, which is a function that considers the density of the prey. We choose a simple linear function gx=ax, that is, the Holling type-I functional response in this paper, where a is the predator predation rate. In the traditional system, g⋅=h⋅; however, in a dissipative system, hx,y=cgx,y, where c is the prey-to-predator conversion rate and 0<c<1. As is well known, restricted by the competition of the same species, environmental factors, and other mechanisms, populations in nature cannot grow indefinitely. Allowing the prey population to grow logically is a common modeling method and has the form xlx=xr1−x/K, where r is the intrinsic growth rate of prey and K is the carrying capacity. In summary, the interaction between the predator and prey can be expressed by the following equations:(2)dxdt=rx1−xK−axy,dydt=caxy−δy.Taking the long-term view, predators can impact the population density of prey only through direct killing. However, the mere presence of predators can change the physiology of prey and even reduce its birth rate. This phenomenon can exert an effect on the population density of prey that is even more obvious than direct predation and is called the fear effect [6–8]. Some animals show a series of antipredator responses, such as habitat migration and foraging behavior change. For example, birds have changed habitats because of a fear of predators in the same habitat. Although new habitats are safer, the relocated animals are unable to adapt to new habitats, and a decline in population ensues [9, 10]. Wang et al. [11] established a model with a fear effect and Holling functional response. They first formulated the fear-effect term c and combined the fear factor with a predator-prey model widely [12]. Here, f represents the fear level that depends on the predator y. The fear effect is introduced to further change the prey birth rate by the following equation:(3)dxdt=nf,yxlx,and nf,y must meet the following conditions before it can be used as a fear function [11]:(4)n0,y=1,nf,0=1,limf⟶∞nf,y=0,limy⟶∞nf,y=0,∂nf,y∂f<0,∂nf,y∂y<0.A simple functionnf,y=1/1+fy satisfies all the above conditions.In population dynamics, the cluster lifestyle of biological populations is beneficial to the growth of populations, but the excessive density of clusters often causes the growth of biological populations to be suppressed due to competition for resources. The population density is too sparse or too dense to be unfavorable for the survival and development of the population. For each population, there must be an independent optimal density for growth and reproduction; this mechanism is called the Allee effect [13–15]. The Allee effect is generally divided into two types: strong and weak [16, 17]. When the population has a weak Allee effect, although the population density grows slowly, the increase in population density is positive. However, when a species with a sparse population density has a strong effect, the population will eventually become extinct. Since the Allee effect is widely present in natural populations, taking it into account in modeling will make the model more complete. Considering the weak Allee effect ϕx=x/x+AA>0 [18–20] and the fear effect, a further modification is required.(5)dxdt=rx1+fy1−xKxx+A−axy,dydt=caxy−δy.To intuitively understand the impact of fear on the model, the relation betweendx/dt and f is shown in Figure 1 [21]. With the increase of the fear parameter, the value of dx/dt decreases.Figure 1
Relation betweendx/dt and f of model (5) (y=1, r=2.65, K=10, A=5, a=0.02, c=0.22, and δ=1.04).There are many types of time delay in the biological models of predators and prey. It is considered herein that the reproduction of the prey population will not be completed instantly and that there is a gestation delayτ required for gestation of prey [22, 23]. To make model (5) closer to reality, gestation delay is introduced as follows:(6)dxtdt=rxt1+fyt1−xt−τKxtxt+A−axtyt,dytdt=caxtyt−δyt.Human beings have a long history of using mathematics to explore the world of biological populations and have made tremendous progress. Nowadays, mathematical models have been widely used in many fields. Among them, there are many works on the application of mathematical models to study the dynamics of populations and infectious diseases [24–26]. Establishing mathematical models is a widely used and effective method of studying the dynamic behavior of populations. The advantage of this method is that it is convenient to use the computer to process the main variables and parameters of a model, save time and cost, and make predictions based on past and present information [27]. The fear effect and Allee effect are ecological effects that actually exist in populations, so incorporating rich ecological effects into the mathematical model and constructing a suitable model to study the impact of ecological effects on the dynamic behavior of the model are our research motivation. Although our model is not complicated, an overly complicated model will not be able to facilitate a detailed analysis due to the limitations of research methods. The objective of this study is to explore how these ecological effects will affect the dynamic behavior of the population when there are Allee effects, fear effects, and time delay in the model. In addition, it will be interesting to explore the bifurcation caused by the Allee effect and time delay.


## 2. Analysis of the Predator-Prey Model

Model (5) has three equilibrium points E0=0,0, EK=K,0, and E∗=x∗,y∗. Here, x∗=δ/ac and y∗=1+4frx∗K−x∗/Kax∗+A−1/2f. According to the actual meaning of the parameters, δ/ac<K is obtained.


### 2.1. Boundedness

Theorem 1.
The solutionxt,yt of the model satisfies(7)limt⟶∞supxt+1cyt≤K4rδδ+r2.Proof.
DefineΛ=xt+1/cyt.(8)dΛdt=dxdt+1cdydt=rx1+fy1−xKxx+A−axy+1ccaxy−1cδy=rx1+fy1−xKxx+A−1cδy.
Then, for allt>0,(9)dΛdt+δΛ=rx1+fy1−xKxx+A−1cδy+δx+1cδy=rx1+fy1−xKxx+A+δx≤rx1−xK+δx=xr+δ−rxK≤K4rδ+r2.
Thus,(10)limt⟶∞supxt+1cyt≤K4rδδ+r2.


### 2.2. Stability Analysis

Theorem 2.
(1)
Trivial equilibrium pointE0 is stable.(2)
Ifa>δ/cK, the axial equilibrium point EK is a saddle point and is stable, otherwise.(3)
Interior equilibriumE∗ is locally asymptotically stable for(11)A<rx∗2K−3x∗+r2x∗22K−3x∗2+4ay∗K1+fy∗rx∗2x∗−K2ay∗K1+fy∗−δac,and E∗ is unstable, otherwise.Proof.
Let(12)Fx,y=rx1+fy1−xKxx+A−axy,Gx,y=caxy−δy.
The Jacobian matrix of model (3) is given by(13)J=∂F∂x∂F∂y∂G∂x∂G∂y,where(14)∂F∂x=Krx2K−3x1+fyx+A+rx2Kx−K1+fyK1+fyx+A2−ay,∂F∂y=rx2fx−KKx+A1+fy2−ax,∂G∂x=cay,∂G∂y=cax−δ.
The community matrix atE0=0,0, EK=K,0, and E∗=x∗,y∗ is given by(15)JE0=000−δ,JEK=−rK1+fyK+A−aK0caK−δ,JE∗=Krx∗2K−3x∗x∗+A+rx∗2Kx∗−KK21+fy∗x∗+A2−ay∗rx∗2fx∗−KKx∗+A1+fy∗2−δccay∗0.
It is found thatE0 is stable. The first eigenvalue of JEK is negative. If a<δ/cK, EK is stable, and if a>δ/cK, EK is a saddle point.
Next, the community matrix at the interior equilibriumE∗ is considered. The characteristic polynomial of JE∗ is(16)Hλ=λ2−TrJE∗λ+DetJE∗,where(17)TrJE∗=Krx∗2K−3x∗x∗+A+rx∗2Kx∗−KK21+fy∗x∗+A2−ay∗,DetJE∗=δay∗−rx∗2fx∗−Kcay∗Kx∗+A1+fy∗2.
One must obtain the conditions ofTrJE∗<0DetJE∗>0; that is, if(18)A<rx∗2K−3x∗+r2x∗22K−3x∗2+4ay∗K1+fy∗rx∗2x∗−K2ay∗K1+fy∗−δac,E∗ is locally asymptotically stable. If(19)A>rx∗2K−3x∗+r2x∗22K−3x∗2+4ay∗K1+fy∗rx∗2x∗−K2ay∗K1+fy∗−δac,E∗ is unstable.


### 2.3. Bifurcation Analysis

#### 2.3.1. Transcritical Bifurcation

Theorem 3.
Ata=a0, where a0=δ/cK, model (5) enters into transcritical bifurcation around EK.Proof.
The eigenvalues ofJEK are −rK/1+fyK+A and λ2=caK−δ. λ1 is negative when a=δ/cK=a0 and λ2=0. Thus, JEK has a zero eigenvalue. V and W are the eigenvectors corresponding to the zero eigenvalue of the matrices JEK and JEKT, where V=−a1+fyK+A/r,1T and W=0,1T.(20)JEKT=−rK1+fyK+A0−aK0,WTfax¯,y¯,a¯=0,WTDfax¯,y¯,a¯V=cK≠0,WTD2fax¯,y¯,a¯V,V=−2ca21+fyK+Ar≠0.
By the Sotomayor theorem [28], the model attains transcritical bifurcation around EK at a=a0.


#### 2.3.2. Hopf Bifurcation

Model (5) undergoes bifurcation when(21)A=rx∗2K−3x∗+r2x∗22K−3x∗2+4ay∗K1+fy∗rx∗2x∗−K2ay∗K1+fy∗−δac.The model (5) undergoes a Hopf bifurcation will be proved in this section.The parameterA0 is chosen to analyze the Hopf bifurcation at E∗=x∗,y∗, where(22)A0=rx∗2K−3x∗+r2x∗22K−3x∗2+4ay∗K1+fy∗rx∗2x∗−K2ay∗K1+fy∗−δac.WhenA=A0, one has TJE∗=0, and λ=±i2Kfx∗+Acay∗/K1+fy∗x∗+A2+δay∗ are imaginary eigenvalues of the Jacobian matrix JE∗. Let λ=αA±iβA be the roots of λ2−TJE∗+DJE∗=0,(23)α2−β2−aT+D=0,2αβ−Tβ=0.Then,(24)α=T2,β=4D−T22,dαdAA=A0=12ddATrJE∗A=A0≠0.The cross-sectional condition is satisfied, so model (3) undergoes a Hopf bifurcation at E∗=x∗,y∗ when A=A0. Next, the detailed nature of the Hopf bifurcation is analyzed. Transforming the equilibrium point E∗ into the origin by setting x=X+x∗ and y=Y+y∗, one has the following model:(25)dXdt=rX+x∗1+fY+y∗1−X+x∗kX+x∗X+x∗+A−aX+x∗Y+y∗,dYdt=caX+x∗Y+y∗−δY+y∗.Expanding Taylor’s series of the above model atX,Y=0,0,(26)X˙=c10X+c01Y+c20X2+c11XY+c02Y2+c30X3+c21X2Y+c12XY2+c03Y3+o,Y˙=d10X+d01Y+d20X2+d11XY+d02Y2+d30X3+d21X2Y+d12XY2+d03Y3+o,where(27)c10=Krx∗2K−3x∗1+fy∗x∗+A+rx∗2Kx∗−K1+fy∗K1+fy∗x∗+A2−ay∗,c01=rx∗2fx∗−KKx∗+A1+fy∗2−δc,c20=r2K−6x∗1+fy∗2x∗+A2−rx∗2K−3x∗1+fy∗2x∗+A−2rx∗2x∗−K1+fy∗22K1+fy∗x∗+A3,c11=−frx∗2K−3x∗x∗+A+frx∗2x∗−12K1+fy∗x∗+A2−a,c02=−2rx∗2f2x∗−KKx∗+A1+fy∗3,c30=r15x∗+6A−2Kx∗+A2+2rx∗x∗−2Ax∗−K−rx∗2K−3x∗12Kx∗+A41+fy∗+r2x∗+A2K−3x∗−2rx∗212Kx∗+A31+fy∗,c21=−fr2K−6x∗x∗+A21+fy∗+2frx∗2x∗−K1+fy∗2K1+fy∗x∗+A3++fr3x∗−22K−3x∗2K1+fy∗x∗+A2,c12=frx∗2K−3x∗x∗+A2+frx∗2x∗−Kx∗+A6K1+fy∗x∗+A3,c03=6rx∗2f3x∗−KKx∗+A1+fy∗4,d10=cay∗,d01=0,d20=0,d11=ca2,d02=0,d30=0,d21=0,d12=0,d03=0.Let(28)P1X,Y=c20X2+c11XY+c02Y2+c30X3+c21X2Y+c12XY2+c03Y3,P2X,Y=d20X2+d11XY+d02Y2+d30X3+d21X2Y+d12XY2+d03Y3.Now,(29)X˙=c10X+c01Y+P1X,Y,Y˙=d10X+d01Y+P2X,Y.Taking the transformationu=X and v=−1/βc10X+c01Y,(30)dudt=−βv+G1u,v,dvdt=βu+G2u,v,where(31)G1u,v=P1u,−1c01c10u+βv,G2u,v=−1βc10P1u,−1c01c10u+βv+c01P2u,−1c01c10u+βv.Therefore,(32)G1u,v=c20u2+c11u−1c01c10u+βv+c02−1c01c10u+βv2+c30u3+c21u2−1c01c10u+βv+c12u−1c01c10u+βv2+c03−1c01c10u+βv3,G2u,v=−1βc10c20u2+c11u−1c01c10u+βv+c02−1c01c10u+βv2+c30u3+c21u2−1c01c10u+βv+c12u−1c01c10u+βv2+c03−1c01c10u+βv3+c01d11u−1c01c10u+βv.The following is set:(33)ψ=116∂3G1∂u3+∂3G1∂u∂v2+∂3G2∂u2∂v+∂3G2∂v3+116β∂2G1∂u∂v∂2G1∂u2+∂2G1∂v2−∂2G2∂u∂v∂2G2∂u2+∂2G2∂v2−∂2G1∂u2∂2G2∂u2+∂2G1∂v2∂2G2∂v2,where(34)∂3G1∂u3=6c30−c10c21c01+c102c12c012−c103c03c013,∂3G1∂u∂v2=2β2c12c012−2c10c03c013−c10c03c013,∂2G1∂u∂v=β2c10c02c012−c11c01−2βuc21c01−2c10c12c012+c102c03c013+2c102c03c013+2β2vc12c012−2c10c03c013−c10c03c013,∂2G1∂u2=2c20−c10c11c01+c102c02c012+6uc30−c10c21c01+c102c12c012−c103c03c013+2βv2c10c12c012−c21c01−c102c03c013−2c102c03c013,∂2G1∂v2=2β2uc12c012−2c10c03c013−c10c03c013+2c02c012β2−6c03c013β3v,∂3G2∂v3=6c10c03c013β2,∂3G2∂u2∂v=2c10c21c01−2c102c12c012+c103c03c013+2c103c03c013,∂2G2∂u∂v=2uc10c21c01−2c102c12c012+c103c03c013+2c103c03c013+2βvc102c03c013−c10c12c012+2c102c03c013+c10c11c01−2c102c02c012−d11β,∂2G2∂u2=6uc102c21c01β−c10c30β−c103c12c012β+c104c03c013β+2vc10c21c01−2c102c12c012+c103c03c013+2c103c03c013+2c102c11c01β−c10c20β−c103c02c012β−c10d11,∂2G2∂v2=2βuc102c03c013−c10c12c012+2c102c03c013−2c10c02c012β+6c10c03c013β2v.Ifψ<0, the Hopf bifurcation is supercritical, and the Hopf bifurcation is subcritical if ψ>0.


## 2.1. Boundedness

Theorem 1.
The solutionxt,yt of the model satisfies(7)limt⟶∞supxt+1cyt≤K4rδδ+r2.Proof.
DefineΛ=xt+1/cyt.(8)dΛdt=dxdt+1cdydt=rx1+fy1−xKxx+A−axy+1ccaxy−1cδy=rx1+fy1−xKxx+A−1cδy.
Then, for allt>0,(9)dΛdt+δΛ=rx1+fy1−xKxx+A−1cδy+δx+1cδy=rx1+fy1−xKxx+A+δx≤rx1−xK+δx=xr+δ−rxK≤K4rδ+r2.
Thus,(10)limt⟶∞supxt+1cyt≤K4rδδ+r2.


## 2.2. Stability Analysis

Theorem 2.
(1)
Trivial equilibrium pointE0 is stable.(2)
Ifa>δ/cK, the axial equilibrium point EK is a saddle point and is stable, otherwise.(3)
Interior equilibriumE∗ is locally asymptotically stable for(11)A<rx∗2K−3x∗+r2x∗22K−3x∗2+4ay∗K1+fy∗rx∗2x∗−K2ay∗K1+fy∗−δac,and E∗ is unstable, otherwise.Proof.
Let(12)Fx,y=rx1+fy1−xKxx+A−axy,Gx,y=caxy−δy.
The Jacobian matrix of model (3) is given by(13)J=∂F∂x∂F∂y∂G∂x∂G∂y,where(14)∂F∂x=Krx2K−3x1+fyx+A+rx2Kx−K1+fyK1+fyx+A2−ay,∂F∂y=rx2fx−KKx+A1+fy2−ax,∂G∂x=cay,∂G∂y=cax−δ.
The community matrix atE0=0,0, EK=K,0, and E∗=x∗,y∗ is given by(15)JE0=000−δ,JEK=−rK1+fyK+A−aK0caK−δ,JE∗=Krx∗2K−3x∗x∗+A+rx∗2Kx∗−KK21+fy∗x∗+A2−ay∗rx∗2fx∗−KKx∗+A1+fy∗2−δccay∗0.
It is found thatE0 is stable. The first eigenvalue of JEK is negative. If a<δ/cK, EK is stable, and if a>δ/cK, EK is a saddle point.
Next, the community matrix at the interior equilibriumE∗ is considered. The characteristic polynomial of JE∗ is(16)Hλ=λ2−TrJE∗λ+DetJE∗,where(17)TrJE∗=Krx∗2K−3x∗x∗+A+rx∗2Kx∗−KK21+fy∗x∗+A2−ay∗,DetJE∗=δay∗−rx∗2fx∗−Kcay∗Kx∗+A1+fy∗2.
One must obtain the conditions ofTrJE∗<0DetJE∗>0; that is, if(18)A<rx∗2K−3x∗+r2x∗22K−3x∗2+4ay∗K1+fy∗rx∗2x∗−K2ay∗K1+fy∗−δac,E∗ is locally asymptotically stable. If(19)A>rx∗2K−3x∗+r2x∗22K−3x∗2+4ay∗K1+fy∗rx∗2x∗−K2ay∗K1+fy∗−δac,E∗ is unstable.


## 2.3. Bifurcation Analysis

### 2.3.1. Transcritical Bifurcation

Theorem 3.
Ata=a0, where a0=δ/cK, model (5) enters into transcritical bifurcation around EK.Proof.
The eigenvalues ofJEK are −rK/1+fyK+A and λ2=caK−δ. λ1 is negative when a=δ/cK=a0 and λ2=0. Thus, JEK has a zero eigenvalue. V and W are the eigenvectors corresponding to the zero eigenvalue of the matrices JEK and JEKT, where V=−a1+fyK+A/r,1T and W=0,1T.(20)JEKT=−rK1+fyK+A0−aK0,WTfax¯,y¯,a¯=0,WTDfax¯,y¯,a¯V=cK≠0,WTD2fax¯,y¯,a¯V,V=−2ca21+fyK+Ar≠0.
By the Sotomayor theorem [28], the model attains transcritical bifurcation around EK at a=a0.


### 2.3.2. Hopf Bifurcation

Model (5) undergoes bifurcation when(21)A=rx∗2K−3x∗+r2x∗22K−3x∗2+4ay∗K1+fy∗rx∗2x∗−K2ay∗K1+fy∗−δac.The model (5) undergoes a Hopf bifurcation will be proved in this section.The parameterA0 is chosen to analyze the Hopf bifurcation at E∗=x∗,y∗, where(22)A0=rx∗2K−3x∗+r2x∗22K−3x∗2+4ay∗K1+fy∗rx∗2x∗−K2ay∗K1+fy∗−δac.WhenA=A0, one has TJE∗=0, and λ=±i2Kfx∗+Acay∗/K1+fy∗x∗+A2+δay∗ are imaginary eigenvalues of the Jacobian matrix JE∗. Let λ=αA±iβA be the roots of λ2−TJE∗+DJE∗=0,(23)α2−β2−aT+D=0,2αβ−Tβ=0.Then,(24)α=T2,β=4D−T22,dαdAA=A0=12ddATrJE∗A=A0≠0.The cross-sectional condition is satisfied, so model (3) undergoes a Hopf bifurcation at E∗=x∗,y∗ when A=A0. Next, the detailed nature of the Hopf bifurcation is analyzed. Transforming the equilibrium point E∗ into the origin by setting x=X+x∗ and y=Y+y∗, one has the following model:(25)dXdt=rX+x∗1+fY+y∗1−X+x∗kX+x∗X+x∗+A−aX+x∗Y+y∗,dYdt=caX+x∗Y+y∗−δY+y∗.Expanding Taylor’s series of the above model atX,Y=0,0,(26)X˙=c10X+c01Y+c20X2+c11XY+c02Y2+c30X3+c21X2Y+c12XY2+c03Y3+o,Y˙=d10X+d01Y+d20X2+d11XY+d02Y2+d30X3+d21X2Y+d12XY2+d03Y3+o,where(27)c10=Krx∗2K−3x∗1+fy∗x∗+A+rx∗2Kx∗−K1+fy∗K1+fy∗x∗+A2−ay∗,c01=rx∗2fx∗−KKx∗+A1+fy∗2−δc,c20=r2K−6x∗1+fy∗2x∗+A2−rx∗2K−3x∗1+fy∗2x∗+A−2rx∗2x∗−K1+fy∗22K1+fy∗x∗+A3,c11=−frx∗2K−3x∗x∗+A+frx∗2x∗−12K1+fy∗x∗+A2−a,c02=−2rx∗2f2x∗−KKx∗+A1+fy∗3,c30=r15x∗+6A−2Kx∗+A2+2rx∗x∗−2Ax∗−K−rx∗2K−3x∗12Kx∗+A41+fy∗+r2x∗+A2K−3x∗−2rx∗212Kx∗+A31+fy∗,c21=−fr2K−6x∗x∗+A21+fy∗+2frx∗2x∗−K1+fy∗2K1+fy∗x∗+A3++fr3x∗−22K−3x∗2K1+fy∗x∗+A2,c12=frx∗2K−3x∗x∗+A2+frx∗2x∗−Kx∗+A6K1+fy∗x∗+A3,c03=6rx∗2f3x∗−KKx∗+A1+fy∗4,d10=cay∗,d01=0,d20=0,d11=ca2,d02=0,d30=0,d21=0,d12=0,d03=0.Let(28)P1X,Y=c20X2+c11XY+c02Y2+c30X3+c21X2Y+c12XY2+c03Y3,P2X,Y=d20X2+d11XY+d02Y2+d30X3+d21X2Y+d12XY2+d03Y3.Now,(29)X˙=c10X+c01Y+P1X,Y,Y˙=d10X+d01Y+P2X,Y.Taking the transformationu=X and v=−1/βc10X+c01Y,(30)dudt=−βv+G1u,v,dvdt=βu+G2u,v,where(31)G1u,v=P1u,−1c01c10u+βv,G2u,v=−1βc10P1u,−1c01c10u+βv+c01P2u,−1c01c10u+βv.Therefore,(32)G1u,v=c20u2+c11u−1c01c10u+βv+c02−1c01c10u+βv2+c30u3+c21u2−1c01c10u+βv+c12u−1c01c10u+βv2+c03−1c01c10u+βv3,G2u,v=−1βc10c20u2+c11u−1c01c10u+βv+c02−1c01c10u+βv2+c30u3+c21u2−1c01c10u+βv+c12u−1c01c10u+βv2+c03−1c01c10u+βv3+c01d11u−1c01c10u+βv.The following is set:(33)ψ=116∂3G1∂u3+∂3G1∂u∂v2+∂3G2∂u2∂v+∂3G2∂v3+116β∂2G1∂u∂v∂2G1∂u2+∂2G1∂v2−∂2G2∂u∂v∂2G2∂u2+∂2G2∂v2−∂2G1∂u2∂2G2∂u2+∂2G1∂v2∂2G2∂v2,where(34)∂3G1∂u3=6c30−c10c21c01+c102c12c012−c103c03c013,∂3G1∂u∂v2=2β2c12c012−2c10c03c013−c10c03c013,∂2G1∂u∂v=β2c10c02c012−c11c01−2βuc21c01−2c10c12c012+c102c03c013+2c102c03c013+2β2vc12c012−2c10c03c013−c10c03c013,∂2G1∂u2=2c20−c10c11c01+c102c02c012+6uc30−c10c21c01+c102c12c012−c103c03c013+2βv2c10c12c012−c21c01−c102c03c013−2c102c03c013,∂2G1∂v2=2β2uc12c012−2c10c03c013−c10c03c013+2c02c012β2−6c03c013β3v,∂3G2∂v3=6c10c03c013β2,∂3G2∂u2∂v=2c10c21c01−2c102c12c012+c103c03c013+2c103c03c013,∂2G2∂u∂v=2uc10c21c01−2c102c12c012+c103c03c013+2c103c03c013+2βvc102c03c013−c10c12c012+2c102c03c013+c10c11c01−2c102c02c012−d11β,∂2G2∂u2=6uc102c21c01β−c10c30β−c103c12c012β+c104c03c013β+2vc10c21c01−2c102c12c012+c103c03c013+2c103c03c013+2c102c11c01β−c10c20β−c103c02c012β−c10d11,∂2G2∂v2=2βuc102c03c013−c10c12c012+2c102c03c013−2c10c02c012β+6c10c03c013β2v.Ifψ<0, the Hopf bifurcation is supercritical, and the Hopf bifurcation is subcritical if ψ>0.


## 2.3.1. Transcritical Bifurcation

Theorem 3.
Ata=a0, where a0=δ/cK, model (5) enters into transcritical bifurcation around EK.Proof.
The eigenvalues ofJEK are −rK/1+fyK+A and λ2=caK−δ. λ1 is negative when a=δ/cK=a0 and λ2=0. Thus, JEK has a zero eigenvalue. V and W are the eigenvectors corresponding to the zero eigenvalue of the matrices JEK and JEKT, where V=−a1+fyK+A/r,1T and W=0,1T.(20)JEKT=−rK1+fyK+A0−aK0,WTfax¯,y¯,a¯=0,WTDfax¯,y¯,a¯V=cK≠0,WTD2fax¯,y¯,a¯V,V=−2ca21+fyK+Ar≠0.
By the Sotomayor theorem [28], the model attains transcritical bifurcation around EK at a=a0.


## 2.3.2. Hopf Bifurcation

Model (5) undergoes bifurcation when(21)A=rx∗2K−3x∗+r2x∗22K−3x∗2+4ay∗K1+fy∗rx∗2x∗−K2ay∗K1+fy∗−δac.The model (5) undergoes a Hopf bifurcation will be proved in this section.The parameterA0 is chosen to analyze the Hopf bifurcation at E∗=x∗,y∗, where(22)A0=rx∗2K−3x∗+r2x∗22K−3x∗2+4ay∗K1+fy∗rx∗2x∗−K2ay∗K1+fy∗−δac.WhenA=A0, one has TJE∗=0, and λ=±i2Kfx∗+Acay∗/K1+fy∗x∗+A2+δay∗ are imaginary eigenvalues of the Jacobian matrix JE∗. Let λ=αA±iβA be the roots of λ2−TJE∗+DJE∗=0,(23)α2−β2−aT+D=0,2αβ−Tβ=0.Then,(24)α=T2,β=4D−T22,dαdAA=A0=12ddATrJE∗A=A0≠0.The cross-sectional condition is satisfied, so model (3) undergoes a Hopf bifurcation at E∗=x∗,y∗ when A=A0. Next, the detailed nature of the Hopf bifurcation is analyzed. Transforming the equilibrium point E∗ into the origin by setting x=X+x∗ and y=Y+y∗, one has the following model:(25)dXdt=rX+x∗1+fY+y∗1−X+x∗kX+x∗X+x∗+A−aX+x∗Y+y∗,dYdt=caX+x∗Y+y∗−δY+y∗.Expanding Taylor’s series of the above model atX,Y=0,0,(26)X˙=c10X+c01Y+c20X2+c11XY+c02Y2+c30X3+c21X2Y+c12XY2+c03Y3+o,Y˙=d10X+d01Y+d20X2+d11XY+d02Y2+d30X3+d21X2Y+d12XY2+d03Y3+o,where(27)c10=Krx∗2K−3x∗1+fy∗x∗+A+rx∗2Kx∗−K1+fy∗K1+fy∗x∗+A2−ay∗,c01=rx∗2fx∗−KKx∗+A1+fy∗2−δc,c20=r2K−6x∗1+fy∗2x∗+A2−rx∗2K−3x∗1+fy∗2x∗+A−2rx∗2x∗−K1+fy∗22K1+fy∗x∗+A3,c11=−frx∗2K−3x∗x∗+A+frx∗2x∗−12K1+fy∗x∗+A2−a,c02=−2rx∗2f2x∗−KKx∗+A1+fy∗3,c30=r15x∗+6A−2Kx∗+A2+2rx∗x∗−2Ax∗−K−rx∗2K−3x∗12Kx∗+A41+fy∗+r2x∗+A2K−3x∗−2rx∗212Kx∗+A31+fy∗,c21=−fr2K−6x∗x∗+A21+fy∗+2frx∗2x∗−K1+fy∗2K1+fy∗x∗+A3++fr3x∗−22K−3x∗2K1+fy∗x∗+A2,c12=frx∗2K−3x∗x∗+A2+frx∗2x∗−Kx∗+A6K1+fy∗x∗+A3,c03=6rx∗2f3x∗−KKx∗+A1+fy∗4,d10=cay∗,d01=0,d20=0,d11=ca2,d02=0,d30=0,d21=0,d12=0,d03=0.Let(28)P1X,Y=c20X2+c11XY+c02Y2+c30X3+c21X2Y+c12XY2+c03Y3,P2X,Y=d20X2+d11XY+d02Y2+d30X3+d21X2Y+d12XY2+d03Y3.Now,(29)X˙=c10X+c01Y+P1X,Y,Y˙=d10X+d01Y+P2X,Y.Taking the transformationu=X and v=−1/βc10X+c01Y,(30)dudt=−βv+G1u,v,dvdt=βu+G2u,v,where(31)G1u,v=P1u,−1c01c10u+βv,G2u,v=−1βc10P1u,−1c01c10u+βv+c01P2u,−1c01c10u+βv.Therefore,(32)G1u,v=c20u2+c11u−1c01c10u+βv+c02−1c01c10u+βv2+c30u3+c21u2−1c01c10u+βv+c12u−1c01c10u+βv2+c03−1c01c10u+βv3,G2u,v=−1βc10c20u2+c11u−1c01c10u+βv+c02−1c01c10u+βv2+c30u3+c21u2−1c01c10u+βv+c12u−1c01c10u+βv2+c03−1c01c10u+βv3+c01d11u−1c01c10u+βv.The following is set:(33)ψ=116∂3G1∂u3+∂3G1∂u∂v2+∂3G2∂u2∂v+∂3G2∂v3+116β∂2G1∂u∂v∂2G1∂u2+∂2G1∂v2−∂2G2∂u∂v∂2G2∂u2+∂2G2∂v2−∂2G1∂u2∂2G2∂u2+∂2G1∂v2∂2G2∂v2,where(34)∂3G1∂u3=6c30−c10c21c01+c102c12c012−c103c03c013,∂3G1∂u∂v2=2β2c12c012−2c10c03c013−c10c03c013,∂2G1∂u∂v=β2c10c02c012−c11c01−2βuc21c01−2c10c12c012+c102c03c013+2c102c03c013+2β2vc12c012−2c10c03c013−c10c03c013,∂2G1∂u2=2c20−c10c11c01+c102c02c012+6uc30−c10c21c01+c102c12c012−c103c03c013+2βv2c10c12c012−c21c01−c102c03c013−2c102c03c013,∂2G1∂v2=2β2uc12c012−2c10c03c013−c10c03c013+2c02c012β2−6c03c013β3v,∂3G2∂v3=6c10c03c013β2,∂3G2∂u2∂v=2c10c21c01−2c102c12c012+c103c03c013+2c103c03c013,∂2G2∂u∂v=2uc10c21c01−2c102c12c012+c103c03c013+2c103c03c013+2βvc102c03c013−c10c12c012+2c102c03c013+c10c11c01−2c102c02c012−d11β,∂2G2∂u2=6uc102c21c01β−c10c30β−c103c12c012β+c104c03c013β+2vc10c21c01−2c102c12c012+c103c03c013+2c103c03c013+2c102c11c01β−c10c20β−c103c02c012β−c10d11,∂2G2∂v2=2βuc102c03c013−c10c12c012+2c102c03c013−2c10c02c012β+6c10c03c013β2v.Ifψ<0, the Hopf bifurcation is supercritical, and the Hopf bifurcation is subcritical if ψ>0.


## 3. Stability and Hopf Bifurcation Analysis of the Delay Model

LettingXt=xt−x∗ and Yt=yt−y∗, after linearization, model (6) could be expressed as the following matrix:(35)ddtXtYt=A1XtYt+A2xt−τyt−τ,where(36)A1=2rx∗1+fy∗x∗+A−rx∗2K1+fy∗1+fy∗x∗+A2−2rx∗2K1+fy∗x∗+A−rx∗31+fy∗K1+fy∗x∗+A2−ay∗−rfx∗2x∗+A1+fy∗x∗+A2+rKfx∗2x∗+Ak1+fy∗x∗+A2−δccay∗0,A2=−rx∗2k1+fy∗x∗+A000.Let(37)G=2rx∗1+fy∗x∗+A−rx∗2K1+fy∗1+fy∗x∗+A2−2rx∗2K1+fy∗x∗+A−rx∗31+fy∗K1+fy∗x∗+A2−ay∗,H=rx∗2k1+fy∗x∗+A,I=rKfx∗2x∗+Ak1+fy∗x∗+A2−rfx∗2x∗+A1+fy∗x∗+A2.The characteristic polynomial is(38)Hλ=λ2−G−He−λtλ+aδy∗−Icay∗.Theorem 4.
IfH2−G2+2F−2IL>0 and −F2+2IFL−I2L2>0, the interior equilibrium E∗ is locally asymptotically stable, and Hopf bifurcation occurs when τ passes the value(39)τ¯=1ω0arctanω02+IL−FGω0.Proof.
The characteristic equation is(40)λ2−G−He−λtλ+aδy∗−Icay∗=0.
Supposing thatλτ0=iω0 is a solution of Hλ, one has(41)−ω02−G−He−iω0τiω0+aδy∗−Icay∗=0.
It is known that(42)e−ω0τ=cosω0τ−isinω0τ,so one obtains(43)−ω02−G−Hcosω0τ−isinω0τiω0+aδy∗−Icay∗=0.
LetF=aδy∗ and L=Icay∗; then,(44)−ω02−iω0G+iω0Hcosω0τ+Hω0sinω0τ+F−IL=0.
Next, separating real and imaginary parts,(45)Hω0sinω0τ=ω02+IL−F,ω0Hcosω0τ=ω0G,one obtains(46)ω04+ω02H2−G2+2F−2IL−F2+2IFL−I2L2=0.
Assuming that(47)W=ω02,one has(48)W2+WH2−G2+2F−2IL−F2+2IFL−I2L2=0.
IfH2−G2+2F−2IL>0 and −F2+2IFL−I2L2>0, all roots of the equation have negative real parts, which means that the interior equilibrium E∗ is locally asymptotically stable:(49)W=−N±N2−4−F2+2IFL−I2L22,where(50)N=H2−G2+2F−2IL.
If−N>0,N2=4−F2+2IFL−I2L2, and −F2+2IFL−I2L2>0, the equation has a unique positive solution. The interior equilibrium E∗ is stable. If −N>0,N2>4−F2+2IFL−I2L2, and −F2+2IFL−I2L2>0, the equation has two positive solutions.
One has(51)Hω0sinω0τ=ω02+IL−F,ω0Hcosω0τ=ω0G,so(52)sinω0τ=ω02+IL−FHω0,cosω0τ=GH,and one obtains(53)tanω0τ=ω02+IL−FGω0.
This shows that if−N<0 and −F2+2IFL−I2L2<0, the characteristic equation(54)λ2−G−He−λtλ+aδy∗−Icay∗=0has a pair of imaginary eigenvalues λ=±iω0, ητ¯=0 when τ¯j±=1/ω0±arctanω02+IL−F/Gω0+2πj/ω0±,j=0,1,….
Next, the cross-sectional condition holds and is proved.
The characteristic equation is(55)λ2−G−He−λtλ+F−IL=0at j=0 and τ=τ¯=τ¯j.
One obtains(56)2λdλdτ−dλdτG−He−λτ+λGe−λτ−τdλdτ−λ=0,dλdτ−1=τλ−Hλ2G+1λ2e−λτ−2λGe−λτ,dReλdτ−1≠0.


## 4. Numerical Simulations

In this section, the realization of numerical simulation uses Matlab 2016a (MathWorks, USA). Through extensive numerical simulations, several results that verify the above analytical analysis conclusion of models (5) and (6) are obtained. The parameters are fixed as r = 2.65, f = 0.1, K = 900, c = 0.215, and δ  = 1.06, and a, A, and τ are considered controlling parameters. Letting a=0.002<δ/cK, EK is stable. Considering the parameters a = 0.02 and A = 5, which satisfy the conditions of Theorem 2, the interior equilibrium E∗246.5,26.4 is locally asymptotically stable. Let A = 149, the other parameters are unchanged, and the interior equilibrium becomes unstable. A series of calculations show that τ¯=7.752, and E∗ loses its stability when τ crosses through the critical τ¯.First, the stability conditions ofEK and E∗ are obtained. By choosing appropriate parameters for the numerical simulations, the stable phase portraits of EK and E∗ are obtained, as shown in Figures 2 and 3. The stable phase portrait of E∗ without the Allee effect is also obtained (Figure 4, A = 0). By comparing Figures 4–6, it is found that the stability of the equilibrium point E∗ will change as parameter A changes. A bifurcation diagram is provided, so the results look more intuitive. As Figure 7 shows, a bifurcation occurs when the bifurcation parameter value is approximately 65, which is consistent with this study. The population dynamic changes from stable to unstable with an increasing value of parameter A. Finally, the dynamical behaviors of model (6) are studied, and mainly, the effect of gestation delay on the stability of the interior equilibrium point is simulated. Figures 8 and 9 show that the equilibrium point is stable at τ=0.9 and unstable at τ=8. In Figure 10, when f=0, the other parameters remain unchanged, and the model is simulated without the fear effect. The phase portrait is unstable and shows periodicity with increasing time τ=0.9. Upon comparison with Figure 8, the existence of the fear effect promotes the stability of the phase portrait.Figure 2
Phase portrait of model (5); EK is stable (r = 2.65, f = 0.1, K = 900, A = 5, a = 0.002, c = 0.215, and δ = 1.06).Figure 3
Phase portrait of model (5); E∗ is stable (r = 2.65, f = 0.1, K = 900, A = 5, a = 0.02, c = 0.215, and δ = 1.06).Figure 4
Phase portrait of model (5); E∗ is stable (r = 2.65, f = 0.1, K = 900, A = 0, a = 0.02, c = 0.215, and δ = 1.06).Figure 5
Phase portrait of model (5); E∗ is stable (r = 2.65, f = 0.1, K = 900, A = 5, a = 0.02, c = 0.215, and δ = 1.06).Figure 6
Phase portrait of model (5); E∗ is unstable (r = 2.65, f = 0.1, K = 900, A = 149, a = 0.02, c = 0.215, and δ = 1.06).Figure 7
Bifurcation diagram with parameterA (r = 2.65, f = 0.1, K = 900, a = 0.02, c = 0.215, and δ = 1.06).Figure 8
Phase portrait of model (6); E∗ is stable (r = 2.65, f = 0.1, K = 900, A = 0.01, a = 0.02, c = 0.215, δ = 1.06, and τ = 0.9).Figure 9
Phase portrait of model (6); E∗ is unstable (r = 2.65, f = 0.1, K = 900, A = 0.01, a = 0.02, c = 0.215, δ = 1.06, and τ = 8).Figure 10
Phase portrait of model (6); E∗ is unstable (r = 2.65, f = 0, K = 900, A = 0.01, a = 0.02, c = 0.215, δ = 1.06, and τ = 0.9).


## 5. Conclusions

Mathematical models play an important role in many fields, such as ecology [29–31], dynamics [32], physics [33], algorithms [34, 35], and epidemiology [36–38]. In this study, a predator-prey model is established with a fear effect and weak Allee effect, and the effect of the weak Allee effect on the model is studied. All equilibrium points of the model were solved, and stability analysis was performed. Furthermore, the conditions required for the equilibrium points to be stable were determined. The boundedness and stability of the model were analyzed. Model (5), which undergoes a transcritical bifurcation at the axial equilibrium point EK and undergoes a Hopf bifurcation at the interior equilibrium point E∗, was also analyzed. Moreover, the direction of the Hopf bifurcation was determined. After completing these proofs, gestation delay was introduced into model (5) to obtain model (6). Time delay can change the stability of model (6); that is, the stability of the interior equilibrium point of model (6) changes with the change in delay. Finally, parameters were selected that satisfy the conditions for numerical simulation. Stable phase portraits of the equilibrium points EK and E∗ were obtained. By increasing the value of parameter A, it was found that the stability of the equilibrium point E∗ changed. Upon comparison with the situation without the Allee effect, it was found that the Allee effect has a significant impact on the stability of the equilibrium point. Then, phase portraits and a bifurcation diagram with A as the parameter were generated to verify the study. The effect of gestation delay on the stability of the interior equilibrium point was mainly studied using model (6). By giving unstable and stable phase portraits, the gestation delay will change the stability of the original model. When the model has no fear effect, the stable phase portrait becomes unstable. In other words, the fear effect promotes stability.Dynamic complexities are the common characteristics in a variety of systems. The complexity of an ecosystem is reflected in both time and space [39, 40]. From a biological point of view, this research shows that a weak Allee effect and gestation delay can destroy the stability of a species and lead to a decrease in population density. Moreover, the introduction of some ecological effects makes the model closer to representing the complex relationships between populations. Appropriate model construction can capture accurate characteristics of interacting population dynamics. This is an important modeling idea to help us better understand population dynamics. Generally speaking, considering the weak Allee effect, fear effect, and gestation delay, the model proposed in this study is more biologically reasonable than some previous predator-prey models. Our finding suggests that the Allee effect and delay can cause unexpected dynamic predictions of interacting populations, such as stability switching. This result is more or less similar to some previous viewpoints. Ye et al. [4] studied the dynamic behavior of a predator-prey model and found that the Allee effect, searching delay, and digestion delay can all cause Hopf bifurcation. Huang et al. [41] showed that the Allee effect, fear effect, and prey refuge will complicate the dynamic behavior of a predator-prey model. Another important result is that the stability of the model with weak Allee effect and gestation delay is very sensitive to changes in the fear level. Therefore, the model with a Holling type-I functional response and weak Allee effect makes it difficult to resist changes in the fear effect. Hence, a low level of fear will have a significant impact on the stability of the model. Overall, it is found that a weak Allee effect, fear effect, and delay all play significant roles in population dynamics. The previous studies show that the Allee effect and delay can induce Hopf bifurcation and that the fear effect will promote the system stability [4, 22, 42, 43]. Wang et al. [11] considered the fear effect term in the interacting prey-predator model and completed a comprehensive work on fear factor. Qiao et al. [21] showed that the fear effect reduces the population density of prey and causes Hopf bifurcation. Overall, these findings suggest that a weak Allee effect, fear effect, and delay all play significant roles in population dynamics. According to the existing references, the present work covers a variety of ecological effects and obtained appropriate results.



---

*Source: 8095080-2022-02-09.xml*