# Qualitative Analysis of a Spatiotemporal Prey-Predator Model with Additive Allee Effect and Fear Effect

Authors: Changcheng Ke; Ming Yi; Yanfeng Guo
Journal: Complexity
 (2022)

---

## Abstract

A diffusive predator-prey system with both the additive Allee effect and the fear effect in the prey subject to Neumann boundary conditions is considered in this paper. Firstly, non-negative and non-trivial solution a priori estimations are shown. Furthermore, for specific parameter ranges, the absence of non-constant positive solutions is demonstrated. Secondly, we use the linearized theory to investigate the stability of non-negative constant solutions. The spatially homogeneous and non-homogeneous periodic solutions, as well as non-constant steady state solutions, are next investigated by using Allee effect parameters as the bifurcation parameter. Finally, numerical simulation is used to illustrate some theoretical results.


---

## Body

## 1. Introduction

The biodynamics of ecosystems are current hot issues in biology and ecology. The intense effort to understand the pattern formation and mechanisms of spatial diffusion during the late 20th century, especially in the context of biological and ecological contexts, has gradually raised more and more concerns. Especially, in biochemical reactions characterized by interactions of different species, the study on predator-prey types has been studied widely in [1–4].Recently, Allee effect, which was initially introduced by Allee in 1931 [5], has been studied extensively [6–10]. With the development of the theory for reaction-diffusion equations, many scholars have done many mathematical research to better describe the relationship between different species. Especially, introducing the Allee effect into the model makes the dynamic behavior of the model closer to reality. The spatiotemporal complexity of a delayed predator-prey model with double Allee effect was given by [11]. In [12], P. J. Pal and S. Tapan consider a system with a double Allee effect in prey population growth, which are very sensitive to parameter perturbations and position of initial conditions. H. Molla and S. Sarwardi developed a predator-prey model that combines these phenomena, considering variable prey refuge with additive Allee effect on the prey species, and also investigated the appearance of Hopf bifurcations in a neighborhood of the unique interior equilibrium point of the dynamical system [13]. The rich behaviour of the dynamics suggests that both prey refuge and a strong Allee affect are important factors in ecological complexity. For a reaction-diffusion system with double Allee effect induced by fear factors subject to Neumann boundary conditions, for details, please refer to [2]. The dynamical behavior of a reaction-diffusion-advection model with weak Allee effect type growth has been studied in [9]. Han and Dai investigated the spatiotemporal pattern formation and selection driven by nonlinear cross-diffusion of a toxic-phytoplankton-zooplankton model with Allee effect. By taking cross-diffusion rate as bifurcation parameter, amplitude equations under nonlinear cross-diffusion are derived that describe the spatiotemporal dynamics [14].Some researchers have indicated that predators can not only capture prey directly but also affect the behavior of prey, even that it could affect the prey more influential than predation [15, 16]. In fact, all animals show various kinds of antipredator responses, such as feeling of fear, habitat changes, vigilance, foraging, and different physiological changes ([17–21]).The cost of fear is objective, and it should be taken into consideration when establishing predation and predation models. For example, Jana et al. [22] have explored the influence of habitat complexity on a predator-prey system under fear effect by incorporating self-diffusion. Tiwari et al. analyzed a predator-prey interaction model with Beddington-DeAngelis functional response (BDFR) and incorporating the cost of fear into prey reproduction. For the spatial system, the Hopf bifurcation around the interior equilibrium, stability of homogeneous steady state, direction, and stability of spatially homogeneous periodic orbits have been established [23]. For a plankton-fish model with both the zooplankton refuge and the fear effect, the local and global dynamics of such a model have been investigated in [24]. Moreover, the investigation in [25] has revealed the threshold behavior of a stochastic predator-prey system with fear effect, prey refuge, and non-constant mortality rate. Sasmal and Takeuchi studied the dynamics of a prey-predator interaction model using Monod–Haldane type functional response and provided detailed mathematical results, including basic dynamical properties, existence of positive equilibria, asymptotic stability of all equilibria, Hopf bifurcation, direction, and stability of bifurcated periodic solutions [26]. Furthermore, they also investigated the role of predation fear and its carry-over effects in the prey-predator model. Basic dynamical properties, as well as the global stability of each equilibrium, have been discussed [27].Allee effect comes in different forms, including multiplicative Allee effect and additive Allee effect. Furthermore, Dennis [6] first proposed the equation incorporating additive Allee effect:(1)dudt=ru1−uk−mu+a,where m and a are constants, which reflect the degree of Allee effect; m/u+a denotes the additive Allee effect; r is the intrinsic growth rate of prey; k presents capacity. We note that if 0<m<a, then (1) has the weak Allee effect and if m>a, then it has the strong Allee effect.Motivated by the previous works above, we further consider the following reaction-diffusion system with fear effect and additive Allee effect:(2)∂u∂t=d1Δu+ru1−u−mu+a11+fv−buv,x∈Ω,t>0,∂v∂t=d2Δv+cbuv−dv,x∈Ω,t>0,∂νu=∂νv=0,x∈∂Ω,t>0,ux,0=u0x≥0,x∈Ω,vx,0=v0x≥0,x∈Ω,where Δ is the Laplace operator on domains. d1>0,d2>0 meanss the diffusion coefficients. The homogeneous Neumann boundary condition is imposed so that there is no population flow across the boundary, ν denotes the outward normal to the boundary ∂Ω. u, v stand for the density of the prey and predator, respectively; m and a are constants, which reflect the degree of Allee effect; f is a constant, which reflects the degree of fear effect; 1/1+fv and m/u+a denote the fear effect and additive Allee effect, respectively; b represents the modified capture rate; c is the conversion coefficient; r is the intrinsic growth rate of prey; d is the death rate of predator. Then, the steady-state system corresponding to (2) is(3)d1Δu+ru1−u−mu+a11+fv−buv=0,x∈Ω,d2Δv+cbuv−dv=0,x∈Ω,∂νu=∂νv=0,x∈∂Ω.The remainder of the paper is structured as follows. In Section2, we carry out a priori estimates for (3) and the requirements for the nonexistence of non-constant positive solutions. In Section 3, we consider the stability of non-negative constant steady state solutions for system (3). In Section 4, we demonstrate the existence of Hopf bifurcation and steady state bifurcation. In Section 5, we show how the parameters affect the dynamical behavior of the system. Furthermore, we verify the analysis results with the numerical simulation results. In section 6, the paper ends with some conclusions.


## 2. Preliminaries

In this section, we first present some properties of equilibrium solutions of (3) including a priori estimates. Then, we discuss the nonexistence of non-constant positive solutions for certain parameter range. It is an essential part for analysis of the existence of non-constant positive steady states and the global bifurcation. We first recall the maximum principle in [28].Lemma 1 (see [28]).
We suppose thatFx,w∈CΩ¯×R. If w∈C2Ω∩C1Ω¯ satisfies(4)Δwx+Fx,wx≥0,x∈Ω∂w∂ν≤0,x∈∂Ω,and wx0=maxΩ¯w, then Fx0,wx0≥0. Similarly, if the two inequalities are reversed and wx0=minΩ¯w, then Fx0,wx0≤0.
We note thatΩ is a bounded domain in RN with smooth boundary. Let λi,i=0,1,2,… be the eigenvalues of −Δ under Neumann boundary condition.
By Lemma1, we have a prior estimates as follows:Theorem 1.
Letux,νx be non-negative and non-trivial solution of (3); we assume that cr/1+ad1+a4/4−m+cd1/d2>0. Then, ux,νxsatisfies(5)0<ux≤1,0<ν≤cr1+ad1+a24−m+cd1d2,where d,r,c,d1,d2,a,m>0.Proof.
From the strong maximum principle, we haveu>0 and v>0. Then, by Lemma 1, it follows 0<ux≤1. The first equation of (3) is multiplied by c and adding the two equations of (3), we obtain(6)−cd1Δu+d2Δv=cru1−u−mu+a11+fv−dv=cruu+a1−uu+a−m11+fv−dv,≤cruu+a1+a24−m11+fv−dv,≤cr1+a1+a24−m−dv,≤cr1+a1+a24−m+dd1cd2−ed2cd1u+d2v,
which leads to(7)Δcd1u+d2v+cr1+a1+a24−m+dd1cd2−dd2cd1u+d2v≥0,under the condition of m≤a. Then, by Lemma 1, we obtain (8)cd1u+d2v≤crd21+ad1+a24−m+cd1,
which implies(9)v≤cr1+ad1+a24−m+cd1d2.Theorem 2.
For any fixedd,r,a,b,c,f>0, there exists d∗r,b,c,d,m,a,f,Ω such that if mind1,d2>d∗, then (3) has no non-constant positive solution.Proof.
Letu,v be a non-negative solution of (3). We denote (10) as(10)u¯=Ω−1∫Ωudx,v¯=Ω−1∫Ωvdx,Fu=ru1−u−mu+a.
Then,(11)∫Ωu−u¯dx=∫Ωv−v¯dx=0.
Multiplying the first equation of (3) by u−u¯ and integrating on Ω , applying Theorem 1 that(12)d1∫Ω∇u−u¯2dx=∫Ωu−u¯Fu11+fvdx−∫Ωbuvu−u¯dx=∫Ωu−u¯Fu11+fv−Fu¯11+fv¯dx−∫Ωbvu−u¯2dx−∫Ωbvu¯u−u¯dx≤∫Ωu−u¯Fu11+fv−Fu¯11+fv+Fu¯11+fv−Fu¯11+fv¯−∫Ωbvu¯u−u¯dx≤r1+ma2∫Ωu−u¯2dx+∫ΩfFu¯u−u¯v−v¯1+fv1+fv¯dx−∫Ωbvu¯u−u¯dx.
Similarly, multiplying the second equation of (3) by v−v¯, we obtain(13)d2∫Ω∇v−v¯2dx=cb∫Ωuvv−v¯dx−∫Ωdvv−v¯dx=cb∫Ωuvv−v¯dx−∫Ωdv−v¯v−v¯dx≤cb∫Ωuvv−v¯−uv¯v−v¯+uv¯v−v¯−u¯v¯v−v¯dx=cb∫Ωuv−v¯2+u−u¯v−v¯v¯dx≤cb∫Ωv−v¯2dx+cb∫Ωuv¯v−v¯dx.
Multiplying the first equation of (3) by c, added to the second equations of (3), and integrating on Ω, we obtain(14)−∫Ωcd1Δu+d2Δvdx=∫Ωcru1−u−mu+a11+fv−dvdx.
Subject it to the boundary conditions, we have(15)d∫Ωvdx=∫Ωcru1−u−mu+a11+fvdx≤Ωcr4.
Hence,(16)v¯=1Ω∫Ωvdx≤cr4d.
By and (16), it follows from Theorem 1 and Young inequality that(17)∫Ωuv¯v−v¯dx=∫Ωv¯u−u¯v−v¯≤cr4d∫Ωu−u¯v−v¯dx≤cr8d∫Ωu−u¯2dx+cr8d∫Ωv−v¯2dx.
Similarly, we have(18)∫Ω−vu¯u−u¯dx=∫Ωv¯−vu¯u−u¯dx≤12∫Ωu−u¯2dx+12∫Ωv−v¯2dx(19)∫ΩfFu¯u−u¯v¯−v1+fv1+fv¯dx=∫Ωfru¯1−u¯−mu¯+au−u¯v¯−v1+fv1+fv¯dx≤rf8+mrf2a∫Ωu−u¯2dx+rf8+mrf2a∫Ωv−v¯2dx.
From (12), (13), (16)–(19) and the Poincaré inequality, we obtain that(20)d1∫Ω∇u−u¯2dx+d2∫Ω∇v−v¯2dx≤1λ1A∫Ω∇u−u¯2dx+B∫Ω∇v−v¯2dx,where(21)A=r1+ma2+bc2r8d+rf8+mrf2a+12,B=bc2r8d+rf8+mrf2a+12+cb.
This shows that if(22)mind1,d2>1λ1maxA,B:=d∗,
then(23)∇u−u¯=∇v−v¯=0.and u,v must be a constant solution.


## 3. Non-Negative Constant Steady-State Solutions

In this section, the stability of non-negative constant steady state solutions of (3) will be investigated by the standard linearization theory. By [17], under particular situations, (3) has the non-negative constant steady state solutions as follows.(1)
the trivial solutionE0=0,0 always exists.(2)
ifa∈0,1, there is no boundary constant solution when a<a+12/4<m.(3)
ifa∈0,1, then E11−a/2,0 is unique boundary equilibria when a<m=a+12/4.(4)
ifa∈0,1, there is two boundary constant solution E21−a−a+12−4m/2,0 and E31−a+a+12−4m/2,0 when a<m<a+12/4.(5)
ifa∈0,1, there is unique boundary constant solution E3 under the condition of 0<m≤a<a+12/4.(6)
ifa=1, there is unique boundary constant solution E41−m,0 only when 0<m<1.(7)
ifa>1, there is unique boundary constant solution E51−a+a+12−4m/2,0 when 0<m<a<a+12/4.(8)
there is unique positive constant solutionE∗d/cb,−b+Δ/2bf=u∗,v∗ with Δ=b2+4bfr1−u∗−m/u∗+a when 1−u∗−m/u∗+a>0.Under the no-flux boundary condition,−Δ has eigenvalues 0=λ0<λ1≤λ2≤… and limi⟶∞λi=∞. Let Xλi be the eigenspace generated by the eigenfunctions corresponding to λi. Let mi be the algebraic multiplicity of λi. Let ϕiji≥0,1≤j≤mi be the normalized eigenfunctions corresponding to λi. Then, the set ϕij1≤j≤mi forms a complete orthonormal basis in L2Ω.Next, we consider the stability of constant steady state solutions.Theorem 3.
For all constantsa,b,c,d,r,f,d1,d2>0, we have that(1)
For trivial solutionE0, ifm>a, thenE0is locally asymptotically stable; ifm<a, thenE0is unstable(2)
If0<a<m=a+12/4<1, thenE1is unstable(3)
If0<a<m<a+12/4<1, thenE2is unstable(4)
Ifcbuj−d<0,Ejj=3,4,5is stable and ifcbuj−d>0,Ejis unstable(5)
E∗exists if and only if−m/a+u∗−u∗+1>0. Ifa>m−u∗, thenE∗is stable. Ifa<m−u∗, thenE∗is unstableProof.
We rewrite (3) as(24)d1Δu+F1u,v=0,x∈Ω,d2Δv+F2u,v=0,kx∈Ω,∂νu=∂νv=0,x∈∂Ω.
The linearization matrix of (3) at a constant solution E = (u0,V0) can be expressed by(25)J=∂uF1u0,v0+d1Δ∂vF1u0,v0∂uF2u0,v0∂vF2u0,v0+d2Δ.where(26)∂uF1u0,v0=−bv+ru−1+m/a+u21+fv+r1−u−m/a+u1+fv,∂vF1u0,v0=−bu−fru1−u−m/a+u1+fv2,∂uF2u0,v0=cbv,∂vF2u0,v0=cbu−d.
We define thatXij=a⋅ϕij:a∈R2, Xi=⊕j=1miXij, and X=⊕i=1∞Xi. Let Φx,Ψx be a pair of eigenfunction of J corresponding to an eigenvalue λ. Then, we have(27)JΦΨ=fu+d1Δfvgugv+d2ΔΦΨ=λΦΨ.
We set(28)Φ=∑0≤i≤∞,1≤j≤miaijϕij,Ψ=∑0≤i≤∞,1≤j≤mibijψij.
Then, we obtain(29)∑0≤i≤∞,1≤j≤mifu+d1λifvgugv+d2λiaijbijϕij=∑0≤i≤∞,1≤j≤miPiaijbijϕij=λaijbijϕij.
From the chapter 5 of [29, 30], we know that if all the eigenvalues of J have negative real parts, then the constant solution E is locally asymptotically stable; J is unstable if there is an eigenvalue of J with positive real part; if all the eigenvalues have non-positive real parts while some eigenvalues have zero real parts, then the stability of E cannot be determined by the linearization. Furthermore, λ is an eigenvalue of J if and only if λ is an eigenvalue of the matrix λI−Pi for some i≥0. We have.(30)λI−Pi=λ2−Tiλ+Di,where(31)Ti=−d1+d2λi+fu+gv,Di=d1d2λi2−d2fu+d1gvλi+fugv−fvgu.(1)
For trival solutionE0=0,0,(32)Ti=−d1+d2λi+r1−ma−d,Di=d1d2λi2+dd1+rma−1d2λi+drma−1.Ifm>a, then for all eigenvalues λ, we have Ti<0 and Di>0, which leads to Reλ<0. Hence, E0 is locally asymptotically stable. If m<a, then for i=0, there exists a positive eigenvalue rm/a−1, which implies that E0 is unstable. In addition, if m=a=1, E0 is stable, else if m=a≠1, E0 is unstable.(2)
ForE11−a/2,0, with 0<a<m=a+12/4<1,(33)Ti=−d1+d2λi+cba−12−d,Di=d1d2λi2+d1d+bc1−a2λi.For corresponding ordinary system,E1 is unstable, so for any d1,d2>0, E1 is unstable.(3)
ForE21−a−a+12−4m/2,0=u2,v2, with 0<a<m<a+12/4<1,(34)Ti=−d1+d2λi+ru2−1+mu2+a2+cbu2−d,Di=d1d2λi2−d2ru2−1+mu2+a2+d1cbu2−dλi+d1d2cbu2−d−ru2+mru2u2+a2.Fori=0, there exists a positive eigenvalue ru2−1+m/u2+a2. So, E2 is always unstable.(4)
Forj=3,4,5, Ej is stable when cbuj−d<0, and in this case, Ti<0,Di>0 for any i≥0. Additionally, in other cases, Ej is unstable, so for any d1,d2>0, Ej is unstable.(35)Ti=−d1+d2λi+rui−1+mui+a2+cbui−d,Di=d1d2λi2−d2rui−1+mui+a2+d1cbui−dλi+rui−1+mui+a2cbui−d.(5)
For positive constant solution,E∗u∗,v∗=d/cb,1/2bf−b+b2+4bfr1−u∗−m/u∗+a. The Jacobi matrix of (3) at E∗ is(36)JE∗=ma+u∗2−1ru∗fv∗+1−d1λi−−ma+u∗−u∗+1fru∗fv∗+12−bu∗cbv∗−d2λi.
It is noted that(37)Ti=−d1+d2λi+ru∗ma+u∗2−11fv∗+1,Di=d1d2λi2−d2ru∗ma+u∗2−11fv∗+1λi+cbv∗−ma+u∗−u∗+1fru∗fv∗+12+bu∗.
ForE∗ exists if and only if −m/a+u∗−u∗+1>0, so it is easy to conclude that Ti<0,Di>0 if m/a+u∗2−1<0a>m−u∗, which implies that E∗ is stable. If a<m−u∗, for i=0, we obtain that Ti>0 and Di>0, so it follows that there exist two of the eigenvalues with positive real parts, which implies that E∗ is unstable.


## 4. Existence of Non-Constant Positive Solutions

In this section, we consider the existence of non-constant positive solutions to (3) in Ω=0,lπ. First, the existence of spatially homogeneous and non-homogeneous periodic solutions is studied by taking m as the bifurcation parameter. Then, the structure and the stability of the bifurcation solutions that bifurcate from u∗,v∗ are shown. From Theorem 3, the stability of u∗,v∗ is determined by the trace and determinant of J. Furthermore, we will restrict −m/a+u∗−u∗+1>0. To put out our discussion into the context of the Hopf bifurcation, we convert (3) into the following system by u˜=u−u∗ and v˜=v−v∗ and drop “∼” for simplicity. We have(38)∂u∂t=d1∂2u∂x2+ru+u∗1−u+u∗−mu+u∗+a11+fv+v∗−bu+u∗v+v∗,x∈0,lπ,t>0,∂v∂t=d2∂2v∂x2+cbu+u∗v+v∗−dv+v∗,x∈0,lπ,t>0,∂νu0,t=∂νulπ,t=∂νv0,t=∂νvlπ,t=0,t>0.Firstly, we define the real-valued Sobolev space(39)X=u,v∈H20,lπ×H20,lπ:∂u0,t∂x=∂ulπ,t∂x=∂v0,t∂x=∂vlπ,t∂x=0,and the corresponding complexification space is given by XC=X⊕iX=a+ib:a,b∈X.The linearized operator of the steady state system of (39) evaluated at m,0,0 is(40)Lm=ru∗m/a+u∗2−11/fv∗+1+d1∂2/∂x2−fru∗−m/a+u∗−u∗+11/fv∗+12−bu∗cbv∗d2∂2/∂x2,where XC is the domain of Lm.The adjoint operator ofLm is defined by(41)L∗m=ru∗m/a+u∗2−11/fv∗+1+d1∂2/∂x2cbv∗−fru∗−m/a+u∗−u∗+11/fv∗+12−bu∗d2∂2/∂x2,where the domain of L∗m is XC.The following condition in [31] is crucial to ensure that the Hopf bifurcation occurs.(H1) There exists a neighborhoodO of m0 such that for m∈O, Lm has a pair of complex, simple, conjugate eigenvalues αm±iωm, continuously differentiable in m, with αm0=0,ωm0=ω0>0, and α′m0≠0, all other eigenvalues of Lm have non-zero real parts for m∈O.Motived by [31], we apply the Hopf bifurcation theory to analyze our system. For the eigenvalue problem(42)−φ′′=λφ,x∈0,lπ,φ′0=φ′lπ=0,we know that the corresponding (42) eigenvalues are λn=n2/l2n=0,1,…, with corresponding eigenfunctions φnx=cosnx/l. Let(43)ϕψ=∑n=0∞anbncosnxl,be a pair of eigenfunctions ofLm corresponding to an eigenvalue ρm, that is, Lmϕ,ψT=ρmϕ,ψT. By a straightforward analysis, we have(44)Lnmanbn=ρmanbn,n=0,1,…,where(45)Lnm=ru∗m/a+u∗2−11/fv∗+1−d1n2/l2−fru∗−m/a+u∗−u∗+11/fv∗+12−bu∗cbv∗−d2n2/l2.Hence, the eigenvalues ofLm are given by the eigenvalues of Lnm,n=0,1,…. The characteristic equation of Lnm is(46)ρ2−Tnmρ+Dnm=0,n=0,1,…,where(47)Tnm=−d1+d2n2l2+ru∗ma+u∗2−11fv∗+1,Dnm=d1d2n4l4−d2ru∗ma+u∗2−11fv∗+1n2l2+cbv∗fru∗−ma+u∗−u∗+11fv∗+12+bu∗.(48)T0=ru∗ma+u∗2−11fv∗+1,D0=cbv∗fru∗−ma+u∗−u∗+11fv∗+12+bu∗>0.Therefore, the eigenvalues are determined by(49)ρm=Tnm±Tn2m−4Dnm2,n=0,1,…If the condition (H1) holds,La has a pair of simple purely imaginary ±iω0 at a=a0, if and only if there exists a unique n∈N such that ±iω0 are the purely imaginary eigenvalues of Lnm. The related eigenvector is denoted by q=qn=an,bnTcosnx/l, with an,bn∈C, such that Lm0q=iω0q.We identify the Hopf bifurcation pointm0 which satisfies the condition (H1): there exists n∈N such that(50)Tnm0=0,Dnm0>0,Tjm0≠0,Djm0≠0forj≠n,and for the unique pair of complex eigenvalues near the imaginary axisαm±iωm(51)α′m0≠0.It is easy to obtainTnm<0 and Dnm>0 if 0<m<a+u∗2, which implies that the steady state u∗,v∗ is locally asymptotically stable. Hence, any potential bifurcation points must be in the interval a+u∗2,a+u∗1−u∗. This means that u∗<1−a/2 is essential for bifurcation condition. For any Hopf bifurcation point m0 in a+u∗2,a+u∗1−u∗, αm±iωm are the eigenvalues of Lnm, where(52)αm=ru∗m/a+u∗2−11/fv∗+12−d1+d2n22l2,ωm=Dnm−α2m,(53)α′m>0,form in a+u∗2,a+u∗1−u∗. Hence, the transversality condition is always satisfied.From the discussion above, the determination of Hopf bifurcation points reduces to describing the set(54)Γ≔m∈a+u∗2,a+u∗1−u∗:for somen∈N,which satisfies the conditionH1.when a set of parameters d1,d2,l,a,b,c,d,f,r are given.In the following, ford1,d2,a,b,c,d,f,r>0 and 0<m<a+u∗1−u∗ fixed, we choose l appropriately. mH=a+u∗2 is always an element of Γ for any l>0 because of T0m0H=0,Tjm0H<0 for any j>1, and Dkm0H>0 for any k∈N. This corresponds to the Hopf bifurcation of spatially homogeneous periodic solution. Apparently, m0H is also the unique value m for the Hopf bifurcation of spatially homogeneous periodic solution for any l>0.In the following, we search for spatially non-homogeneous Hopf bifurcation points forn≥1. As T0m0H=0 and T0′m>0 for m∈m0H,a+u∗1−u∗, we obtain that 0<T0m<T0a+u∗1−u∗=ru∗1−a−2u∗/a+u∗fv∗+1:=M∗ for m∈m0H,a+u∗1−u∗. We define(55)ln=nd1+d2M∗,n∈N∗.Then forln<l≤ln+1, and 1≤j≤n, we derive the root of T0m=d1+d2j2/l2 as mjH such that m0H<mjH<a+u∗1−u∗. Moreover, by T0′m>0 in m0H,a+u∗1−u∗, we derive (56) and (57)(56)0<m0H<m1H<m2H<…<mnH<a+u∗1−u∗(57)TjmjH=0,TimjH≠0fori≠j.SinceDjmjH>0, now we discuss a condition to verify DnmjH≠0 for j≠n. For m∈m0H,a+u∗1−u∗, we have(58)Dim=d1d2i4l4−d2T0mi2l2+D0m:=τi2l2.The quadratic functionτi2/l2 is positive for all l∈R if the discriminant of τi2/l2=0 is negative, which means that (60)(59)d22T02m−4d1d2D0=d22ru∗ma+u∗2−11fv∗+12,−4d1d2cbv∗fru∗−ma+u∗−u∗+11fv∗+12+bu∗.We note that(60)f¯=r2d22u∗4m/a+u∗2−1416d12c21−u∗−m/u∗+a2−b214br1−u∗−m/u∗+a.Fora+u∗2≤m<a+u∗1−u∗, we can choose f>f¯ such that the discriminant of τi2/l2=0 is negative. Then, τi2/l2>0 for i∈N such that DimjH>0.We summarize our analysis above and apply Theorem2 in [31]. The existence of both spatially homogeneous and non-homogeneous periodic solutions bifurcation from u∗,v∗ can be obtained as follows:Theorem 4.
For anylinln,ln+1andf>f¯, system (2) undergoes Hopf bifurcation at eachm=mjH1≤j≤n. Moreover, the bifurcation periodic solutions nearm,u,v=mjH,u∗,v∗can be parameterized asms,us,vsso thatms∈L∞in the form ofms=mjH+osfors∈0,δfor some smallδ>0, (61) and (62)(61)ust,x=u∗+sanei2πt/Ts+a¯ne−i2πt/Tscosnxl+os2,vst,x=v∗+sbnei2πt/Ts+b¯ne−i2πt/Tscosnxl+os2,where an,bn is the corresponding eigenvector, and(62)Ts=2πω01+τ2s2+os4,τ2=−1ω0Imc1mjH−Rec1mjHα′mjHω0′mjH,(63)T′′0=4πω0τ2=−4πω02Imc1mjH−Rec1mjHα′mjHω0′mjH.
Furthermore, we notice that(1)
The bifurcating periodic orbits fromm=m0H are spatially homogeneous, which coincide with the periodic orbits of the corresponding ODE system(2)
The bifurcating periodic orbits fromm=mjH are spatially non-homogeneous.Then, we consider the direction and stability of spatially homogeneous Hopf bifurcation.Theorem 5.
For system (2), if all other eigenvalues of Lnm0H have negative real parts and Rec1m0H<0resp.>0, the spatially homogeneous periodic solutions bifurcating from m=m0H are locally asymptotically stable (resp. unstable). Moreover, the Hopf bifurcation at m0H is supercritical (resp. subcritical) if 1/α′m0HRec1m0H<0resp.>0.Proof.
Here, the notations and calculations in [31] are used in the same way. For the sake of simplicity, we denote 1−u∗−m/u∗+a by M. Then, we introduce(64)q≔a0b0=1−−cb+cb2+4bfrM/2fω0iq∗≔a0∗b0∗=1/2πl−ω0f/−cb+cb2+4bfrMπli,
such thatq∗,q=1,q∗,q¯=0, Lm0Hq=iω0q and L∗m0Hq∗=−iω0q∗, where(65)ω0=2bu∗cMb2+4bfrMb+b2+4bfrM.
Andu,v=∫0lπu¯Tvdx denotes the inner product in L20,lπ×L20,lπ. Then, we get the derivatives at u∗,v∗,m0H as follows:(66)fuu=−2ra3+3a2u∗−am+3au∗2+u∗3a+u∗3fv∗+1,guu=0,fuv=fra22u∗−1+am+2u∗2u∗−1+u∗22u∗−1a+u∗2fv∗+12−b,guv=cb,fvv=2f2ru∗M1+fv∗3,gvv=0,fuuu=−6amra+u∗41+fv∗,guuu=0,fuuv=2fra3+3a2u∗−am+3au∗2+u∗3a+u∗3fv∗+12,guuv=0,fuvv=−2f2ra22u∗−1+am+2u∗2u∗−1+u∗22u∗−1a+u∗2fv∗+13,guvv=0,fvvv=−6f3ru∗M1+fv∗4,gvvv=0.
In addition, we note(67)Qqq=cndncos2nxl,Qqq¯=enfncos2nxl,Qqqq¯=gnhncos3nxl,where cn,dn,en,fn,gn,hn are defined as the same with [31].(68)cn=fuuan2+2fuvanbn+fvvbn2,dn=guuan2+2guvanbn+gvvbn2,en=fuuan2+fuvanb¯n+a¯nbn+fvvbn2,fn=guuan2+guvanb¯n+a¯nbn+gvvbn2,gn=fuuuan2an+fuuv2an2bn+an2b¯n+fuvv2bn2an+bn2a¯n+fvvvbn2bn,hn=guuuan2an+guuv2an2bn+an2b¯n+guvv2bn2an+bn2a¯n+gvvvbn2bn.
Forn=0, by calculation, we derive (69)(69)c0=−2cf2r3M2fv∗+134bfrM+b2+2rm/a+u∗2−1fv∗+1−2mru∗a+u∗3fv∗+1−i2crM−frM/fv∗+12−fru∗m/a+u∗2−1/fv∗+12−bcu∗M4bfrM+b2,d0=−i2c2brMcu∗Mb4frM+b,e0=2cf2r3M2fv∗+134bfrM+b2−2rmu∗a+u∗3fv∗+1+2rm−a+u∗2a+u∗21+fv∗,g0=r32fv∗+14192ic3f3r3u∗M4cu∗Mb2+4bfrM3/2+64cf2r2M2ma+2u∗−1a+u∗2fv∗+1u∗a+u∗a+u∗u∗−1+mb2+4bfrM+64icfrMa+u∗3−amfv∗+12a+u∗3cu∗Mb2+4bfrM+192mu∗fv∗+13a+u∗4−192mfv∗+13a+u∗3,f0=h0=0.
Then, we can obtain (70)(70)q∗,Qqq=bc−cf2r3M2fv∗+13b2+4bfrM−mru∗a+u∗3fv∗+1+rm−a+u∗2a+u∗21+fv∗+icrMb−fra+u∗22u∗−1+am/a+u∗21+fv∗2cu∗Mb2+4bfrM,q∗,Qqq¯=cf2r3M2fv∗+13b2+4bfrM−mru∗a+u∗3fv∗+1+rm−a+u∗2a+u∗2fv∗+1,q¯∗,Qqq=−bc−cf2r3M2fv∗+13b2+4bfrM−mru∗a+u∗3fv∗+1+rm−a+u∗2a+u∗21+fv∗+icrMb−fra+u∗22u∗−1+am/a+u∗21+fv∗2cu∗Mb2+4bfrM,q¯∗,Qqq¯=cf2r3M2fv∗+13b2+4bfrM−mru∗a+u∗3fv∗+1+rm−a+u∗2a+u∗2fv∗+1.q¯∗,Cqqq¯=r64fv∗+14192ic3f3r3u∗M4cu∗Mb2+4bfrM3/2+64cf2r2M2ma+2u∗−1a+u∗2fv∗+1u∗a+u∗a+u∗u∗−1+mb2+4bfrM+64icfrMa+u∗3−amfv∗+12a+u∗3cu∗Mb2+4bfrM+192mu∗fv∗+13a+u∗4−192mfv∗+13a+u∗3.
And, we note(71)ω20=2iω0I−Lm0H−1H20,ω11=−Lm0H−1H11,H20=c0d0−q∗,Qqqa0b0−q¯∗,Qqqa¯0b¯0=0,H11=e0f0−0q∗,Qqq¯a0b0−q¯∗,Qqq¯a¯0b¯0=0.
Hence,ω20=ω11=0, q∗,Qω20q¯=q∗,Qω11q=0. By further calculation, we obtain that(72)Rec1m0H=Rei2ω0q∗,Qqq⋅q∗,Qqq¯+12q∗,Cqqq¯=−12ω0crMb−fra+u∗22u∗−1+am/a+u∗21+fv∗2cu∗Mb2+4bfrM⋅cf2r3M2fv∗+13b2+4bfrM−mru∗a+u∗3fv∗+1+rm−a+u∗2a+u∗2fv∗+1+11281+fv∗464cf2r2M2ma+2u∗−1a+u∗2fv∗+1u∗a+u∗a+u∗u∗−1+mb2+4bfrM+192mu∗fv∗+13a+u∗4−192mfv∗+13a+u∗3.
From above analysis, we know thatα′m0H>0. Hence, by Theorem 2 in [31], the bifurcating solutions bifurcated from m0H,u∗,v∗ are locally asymptotically stable (resp. unstable) if Rec1m0H<0resp.>0 and Tjm0H<0, Djm0H>0 for j≥1, and the Hopf bifurcation at m0H is supercritical (resp. subcritical) if 1/α′m0HRec1m0H<0resp.>0. The proof is complete.
Inspired by [31, 32], we take m as the bifurcation parameter and also restrict a+u∗2≤m<a+u∗1−u∗. We suppose that Ω=0,lπ. The non-negative steady state solutions of (72) satisfy the elliptic problem corresponding to(73)d1∂2u∂x2+ru+u∗1−u+u∗−mu+u∗+a11+fv+v∗−bu+u∗v+v∗=0,x∈0,lπ,t>0,d2∂2v∂x2+cbu+u∗v+v∗−dv+v∗=0,x∈0,lπ,t>0,∂νu0,t=∂νulπ,t=∂νv0,t=∂νvlπ,t=0.
From Theorem3, we know that u∗,v∗ is locally asymptotically stable for 0<m<a+u∗2 and unstable for a+u∗2<m<a+u∗1−u∗.
The steady state bifurcation pointm0 satisfies the steady state bifurcation condition (H2) in [31]:(74)Dnm0=0,Tnm0≠0,andTjm0≠0,Djm0≠0forj≠n∈N0.(75)ddmDnm0≠0.
It means that the potential steady state bifurcation pointsm shall satisfy conditions (74) and (75). Recall that for m∈0,m0H, Tnm≤0 and Dnm>0. Then, any potential bifurcation point m0 must be in the interval m0H,a+u∗1−u∗. Hence, the steady state bifurcation points reduces to the set(76)Λ:=m∈m0H,a+u∗1−u∗:for somen∈N,74and75are satisfied,
when a set of parametersd1,d2,a,b,c,d,f,r,l are fixed.
Recall thatDnm=d1d2ρ2−d2T0mρ+D0, where ρ=n2/l2. By solving Dnm=0 , we have(77)ρ=ρ±m:=d2T0m±d22T02m−4d1d2D0m2d1d2.
We define that(78)zm=d2r2u∗2ma+u∗2−12−4d1cr1−u∗−mu∗+ab2+4bfr1−u∗−mu∗+a,Bm=d22T02m−4d1d2D0m=4b2d2u∗b+b2+4bfr1−u∗−m/u∗+a2zm.
Forzm, we have z′m>0 and za+u∗2<0,za+u∗1−u∗>0 for a+u∗2≤m<a+u∗1−u∗.
Hence, there exists a unique root ofzm=0 denoted by mB, which implies that BmB=0 and ρ±m>0 exists only for mB≤m<a+u∗1−u∗. Therefore, the potential steady state bifurcation points reduces to the set (79).(79)Θ:=m∈mB,a+u∗1−u∗:forsomen∈N,74and75are satisfied.
Then, the properties ofρ±m can be summarized as follows:Lemma 2.
We assume thatd1,d2,f>0,u∗<1−a/2. Then, for any m∈mB,a+u∗1−u∗, ρ±m exists. Moreover, ρ+m is increasing and ρ−m is decreasing.(80)limm⟶mBρ+m=limm⟶mBρ−m=T0mB2d1,limm⟶mBρ+′m=+∞,limm⟶mBρ−′m=−∞,ρ+a+u∗1−u∗=1d1ru∗1−u∗a+u∗−1,ρ−a+u∗1−u∗=0.Proof.
The first limit equation is trivial, so we omit here. We mainly analyze the monotonicity result onρ±a+u∗1−u∗ with respect to m for m∈mB,a+u∗1−u∗.
DifferentiatingDnm with respect to m, it follows that (81)2d1d2ρ±mρ±′m−d2T0mρ±m−d2T0ρ±′m+D0′=0.
Hence,ρ±′m=d2T0′mρ±m−D0′m/2d1d2ρ±−d2T0m. It is easy to get 2d1d2ρ+m−d2T0m>0 and 2d1d2ρ−m−d2T0m<0 from (77). In addition, by calculation, we obtain that for m∈mB,1−u∗a+u∗, d2T0′mρ±m−D0′m>0. The proof is completed.
It follows from Lemma1 that the curve m,ρ± forms a smooth connected curve which connects m,ρ=mB,T0mB/2d1, 1−u∗a+u∗,1/d1ru∗1−u∗/a+u∗−1, and m,ρ=m,0.
By the properties ofρ±, if (82)0<n2l2<1d1ru∗1−u∗a+u∗−1,
then there existsmnB∈mB,1−u∗a+u∗ such that ρ+mnB=n2/l2 or ρ−mnB=n2/l2, and thus DnmnB=0. We define l¯=n/1/d1ru∗1−u∗/a+u∗−1. Then, for any l>l¯, there exists a mnB such that DnmnB=0.
Next, we verifydDnmnB/dm≠0. We recall that Dn′m=−d2n2/l2T0′m+D0′m. Moreover, we know that T0′m>0 and D0′m<0 . It follows that dDnmnB/dm<0.


## 5. Numerical Simulations

In this section, in order to reveal the influence of fear effect, Allee effect, and other factors on the predator-prey model, a numerical method is used to analyze the effect of parameters on the asymptotic behavior of system (2) so as to verify and supplement the theoretical results mentioned before.In Figure1, we choose d1=0.1,d2=0.1,a=0.5,b=1,c=1,d=0.2,r=1,f=15. Varying the parameter m and choosing the initial data near u∗,v∗, we indicate the following numerical results on the effects of parameter m:(1)
Take0<m=0.3<a. Since m<d/bc+a2, u∗,v∗≈0.2000,0.1275 is locally asymptotically stable by Theorem 3. The simulation results indicate that system (2) converges to the equilibrium (see Figures 1(a) and 1(b)).(2)
Takingm=0, there is no Allee effect on prey. Since m<d/bc+a2, u∗,v∗≈0.2000,0.2000 is locally asymptotically stable by Theorem 3. The simulation results indicate that system (2) converges to the equilibrium (see Figures 1(c) and 1(d)).(3)
Taking0<d/bc+a2<m=0.495<a<d/bc+a1−d/bc, it satisfies the condition of the weak Allee effect and Hopf bifurcation condition by Theorem 4. The simulation results indicate that system (2) undergoes Hopf bifurcation (see Figures 1(e) and 1(f)).(4)
Taking0<d/bc+a2<a<m=0.51<d/bc+a1−d/bc, it satisfies the condition of the strong Allee effect and Hopf bifurcation condition by Theorem 4. The simulation results indicate that system (2) undergoes Hopf bifurcation (see Figure 1(g) and 1(h)).Figure 1
The effects of parameterm for T=800,l=2. The values of parameter m are as follows: (a, b) m=0.3; (c, d) m=0; (e, f) m=0.495; (g, h) m=0.51. The right column is the L1 norm of u and v.
(a)(b)(c)(d)(e)(f)(g)(h)In Figure2, we choose d1=0.1,d2=0.1,a=0.5,b=1,c=1,d=0.2,r=1,m=0.51. Varying the parameter f and choosing the initial data near u∗,v∗, it shows that when f is small, system (2) has obviously periodic oscillation (see Figure 2(a)). When f increases, the maximum L1 norms of u are almost the same. However, the maximum L1 norms of v decrease with f increasing (see Figures 2(a)–2(d)). This means that the fear has a negative impact on predators. Moreover, the period of periodic solutions becomes larger as f increases. (see Figures 2(a)–2(d)).Figure 2
The effects of parameterf for T=800 , l=2. The values of parameter f are as follows: (a) f=5; (b) f=50; (c) f=500; (d) f=5000.
(a)(b)(c)(d)In Figure3, we choose d1=0.1,d2=0.1,a=0.5,b=1,c=1,r=1,f=15,m=0.51. Varying the parameter d and choosing the initial data near u∗,v∗, it indicates that when d increases, the period of periodic solution is decreasing. Furthermore, with d increasing, the amplitude of periodic solutions is also decreasing. (see Figures 3(a)–3(c)). As d continues to increase, system (2) converges to an equilibrium (see Figure 3(d)).Figure 3
The effects of parameterd for T=800,l=2. The values of parameter d are as follows: (a) d=0.18; (b) d=0.2; (c) d=0.22; (d) d=0.3.
(a)(b)(c)(d)


## 6. Conclusion

In this paper, a diffusive predator-prey model with additive Allee effects induced by fear factors is considered, in which prey can represent the antipredator behavior due to fear factors. Analytical results indicate that the upper bound ofv depends on the diffusion rates d1, d2, the death rate of the predators d, the Allee effects parameters m,a, the conversion rate c, and the level of fear f. There exists d∗, which depends on a,b,c,d,r,m,f,Ω, such that if mind1,d2>d∗, the system have only constant positive solution. Furthermore, the dynamic behavior near E∗=u∗,v∗ is of more concerned. Then, we indicate the existence of non-constant positive solutions. Taking m as a bifurcation parameter, system undergoes Hopf bifurcation at each m=mjH0≤j≤n. Furthermore, for d1,d2,a,b,c,d,r,f>0, 0<m<1−d/bcd/bc+a are fixed, and there is a smooth curve Γn of non-constant positive solutions bifurcating from u∗,v∗.We observe that the Allee effect is essential to the dynamical behavior of system (2) by numerical simulations. The amplitude in the strong Allee effect is increasing with m increasing. On the other hand, numerical simulations reveal that the fear effect have an impact on the dynamical behavior of system (2). With the fear effect increasing, the period of periodic solutions is increasing, but the maximum L1 norm of u is almost the same. On the contrary, the maximum L1 norm of v decreases with f increasing. From a biological standpoint, the prey survives by adopting antipredator behavior as a result of the fear effect, and the predator is impacted by the prey’s antipredator behavior. At last, we show how the death rate d affects system (2) with the strong Allee effect. As d increases, the amplitude of periodic solution is decreasing, and the period of periodic solution is also decreasing.It is extremely important to construct animal interaction models in the incorporation of these different types of factors. Considering the different ways of introducing the Allee effect and the interaction with the fear effect, further analysis of the bifurcating solutions of (2) remains a challenging problem. From our discussion before, we conjecture that Turing-Hopf bifurcation, Hopf-Hopf bifurcation is likely to exist in the system, which reveals more complex dynamic behavior and potential biological significance.



---

*Source: 5715922-2022-06-02.xml*