A functioning ecosystem depends on the framework of food webs that it supports [1]. Food webs involve many types of predator–prey interactions, which are fundamental interactions for sustaining species [2]. Interaction between the prey and the predator can be affected by factors such as refuge, disease, stage structure, competition, and fear [3–7]. Intraguild predation (IGP), on the other hand, is described as predator–prey interactions among consumers who may be fighting for limited resources. In natural communities, there is a growing body of literature emphasizing the relevance of intraguild predation [8, 9]. Three species are involved in the simplest intraguild predation model: a superpredator (IG predator), a specialist predator (IG prey), and a basal prey. Bai et al. recently suggested a three-species IGP food web model with the IG predator, IG prey, and basal prey, in which the basal prey grows logistically with a large Allee effect [10]. They looked into the model's local and global dynamics, focusing on the impact of the Allee effect and discovered that the intraguild predation food web model has rich and complicated dynamic behavior and that a large Allee effect in the basal prey raises the danger of extinction for not only the basal prey but also the IG prey or/and IG predator.

Many predators in food webs are superpredators, who may not restrict their diets to a specific prey species but feed also on other predators [11]. Therefore, superpredators are expected to compete not only with other predators for food and space but in many cases also through intraguild predation [5, 12, 13]. Predators induce indirect effects such as fear in prey that can change the prey's behavior [14]. Fear takes the form of sustained psychological stress on the prey, as prey species are always wary of possible attack [15]. Suraci et al. experimentally showed that fear of large carnivores reduces mesocarnivore foraging, which benefits the mesocarnivore's prey [16].

A handful of mathematical models have studied the effect of fear on food webs. Panday et al. investigated the impact of fear in a tri-trophic food chain model, with prey fear in response to both predators [15], from the middle predator to the generalist predator [17] and with delays [18]. Cong et al. introduced a fear-adjusted birth rate to a three-species food web [19]. Hossain et al. limited the growth of the prey due to fear in an intraguild predation model [20]. Mukerjee incorporated interspecific competition and fear affecting the death rate of the prey [21]. Ibrahim et al. limited prey growth due to fear of the generalist predator and showed that fear could have a stabilizing effect on the system [22]. Mondal et al. showed that the prey's fear of predators was responsible for the increase in intraspecific competition among the prey species [23]. Roy et al. showed that fear could play a destablizing role if it caused a reduction in the birth rate of susceptible prey, whereas the levels of fear responsible for the increase in the intraspecies competition of susceptible prey and eradication of the disease prevalence could stabilize an otherwise unstable system [24]. Maity et al. considered time-varying fear effects, showing that periodic solutions could arise [25]. Hossain et al. showed that perceived fear of predators could reduce the prey birth rate [26]. Tiwari et al. showed that seasonal variations in the level of prey fear generated higher-order periodic solutions [27].

Here, we examine the impact of fear on the dynamical behavior of an IGP food web system in which the prey responds with fear to the specialist predator but not the superpredator; because the superpredator has alternative food sources, it will be less dangerous to the basal prey compared to the IG prey. To the best of our knowledge, this is the first model to examine this effect.

Our food web consists of a prey at the first level, a specialist predator at the second level, and a superpredator at the third level. Let x(t), y(t), and z(t) be the population densities at time t for the prey, specialist predator, and superpredator, respectively. The prey grows logistically in the absence of predators, while it has a fear property of predation in the presence of the specialist predator. Hence, the intrinsic growth rate of prey becomes r1+ky, which is a monotonic decreasing function of both k and y, where k represents the fear rate [28]. The food transport attack rates are given by the parameters a₁, a₂, and a₃, with conversion rates e₁, e₂, and e₃. Finally, the predators face natural death rates d₁ and d₂ for the specialist predator and superpredator, respectively. The dynamics of the food web with fear can be represented mathematically by the following set of differential equations: (1)dxdt=rx(1-x)(11+ky)-a1xy-a2xz=xf1(x,y,z),dydt=e1a1xy-a3yz-d1y=yf2(x,y,z),                     dzdt=e2a2xz+e3a3yz-d2z=zf3(x,y,z).                   Initial conditions satisfy x(t) ≥ 0, y(t) ≥ 0, and z(t) ≥ 0, and all parameters are assumed to be positive.

The interaction functions f_i (i = 1, 2, 3) are continuous and have continuous partial derivatives and are thus Lipschitzian, so system (1) has a unique solution. Furthermore, for any initial condition in ℝ+3, the solution of system (1) is positive and uniformly bounded as shown in the following theorem. Hence, system (1) will be a dissipative system.

Theorem (1): The domain of system (1), ℝ+3, is positively invariant, and all solutions of system (1) starting in ℝ+3 are uniformly bounded.

Proof. Let x(t), y(t), and z(t) be any solution of system (1). Since the solution (x(t), y(t), z(t)) of the system (1) with initial condition in ℝ+3 exists and is unique on [0 , δ), where 0 < δ ≤ +∞, we have: x(t)=x(0)e∫0t[r(1-x(s))(11+ky(s))-a1y(s)-a2z(s)]ds≥0 y(t)=y(0)e∫0t[e1a1x(s)-a3z(s)-d1]ds≥0 z(t)=z(0)e∫0t[e2a2x(s)+e3a3y(s)-d2]ds≥ 0. From the first equation of system (1): dxdt≤rx(1-x ). Then, it is easy to verify that x(t)≤r4 for all t.

Let = x + y + z. Then, since e_i ∈ (0, 1); i = 1, 2, 3, we have

dWdt≤r4-M[x+y+z ],

where M = min{1, d₁, d₂}. Solving the following linear differential inequality dWdt+MW≤ r4, we obtain that, as t → ∞, W(t)≤r4M.

System (1) has at most five nonnegative biologically feasible equilibrium points. The trivial equilibrium E₀ = (0, 0, 0) always exists. The axial equilibrium E₁ = (1, 0, 0) always exists on the boundary of the first octant. The specialist-free equilibrium point, E2=(x̄,0,z̄), where (2a)x̄=d2e2a2 and z̄=r(1-x̄)a2, exists provided that (2b)d2<e2a2. The superpredator-free equilibrium point, E3=(x^,y^,0), where (3a)x^=d1e1a1 and y^=-12k+12a1ka12+4kr(e1a1-d1)e1, exists, provided that (3b)d1<e1a1. There is also an interior equilibrium point, E4=(x*,y*,z*), where (4a)y*=d2-e2a2x*e3a3 and z*=e1a1x*-d1a3, with x^* a positive root of the following second-order polynomial equation: (4b)β1x2+β2x+β3=0, where β1=e2a1a22k[e1e3-e2],  β2=e2e3a1a2a3-re32a32+2e2a1a2kd2-e1e3a1a2a32-e1e3a1a2kd2-e2e3a22kd1,  β3=re32a32-e3a1a3d2-a1kd22+e3a2a32d1+e3a2kd1d2. Therefore, there is a unique interior equilibrium point in the interior of ℝ+3 provided that the following conditions hold: (4c)d1e1a1<x*<d2e2a2 (4d) β1>0;β3<0ORβ1<0;β3>0}.

Next, we determine the requirements that ensure persistence in the system (1). Because the types of attractors available in the boundary planes affect the creation of our IGP food web model's persistence conditions, an analysis of the dynamics in the boundary planes of the IGP food web system is conducted. System (1) has two subsystems: the first occurs in the absence of the superpredator (IG predator) and the second occurs in the absence of the specialist predator (IG prey). The first subsystem is (5)dxdt=x[r(1-x)1+ky-a1y]=xf1(x,y)dydt=y[e1a1x-d1]=yg1(x,y). The second subsystem can be written as follows: (6)dxdt=x[r(1-x)-a2z]=xf2(x,z)dzdt=z[e2a2x-d2]=zg2(x,z).       Subsystem (5) has a unique positive equilibrium point in the interior of xy-plane given by (x^,y^)≡E3, while subsystem (6) has a unique positive equilibrium point in the interior of xz-plane given by (x̄,z̄)≡E2. The interior equilibrium point of subsystem (5) corresponds to system (1)'s superpredator-free equilibrium point; similarly, the interior equilibrium point of subsystem (6) corresponds to system (1)'s specialist-predator-free equilibrium point. The dynamics around these equilibrium points can be described in the following theorem.

Theorem (2): There are no periodic dynamics in the interior of xy-plane or the xz-plane.

Proof. Define a continuously differential function B(x,y)=1xy. Then, we have: Δ=∂(Bf1)∂x+∂(Bg1)∂y=-ry(1+ky).  Clearly, Δ has the same sign and does not equal zero almost everywhere in a simply connected region of the xy-plane. By the Dulac–Bendixson criterion, system (1) has no periodic solutions lying entirely in the interior of xy−plane. The second part follows by using the Dulac function C(x,z)= 1xz.

Theorem (3): System (1) is uniformly persistent in the interior of ℝ+3, provided that (7a)e1a1>d1, (7b)e2a2>d2, (7c)e1a1x̄>a3z̄+d1, (7d)e2a2x^+e3a3y^>d2.

Proof. Consider the function (x,y,z)=xb1yb2zb3, where b_i, i = 1, 2, 3, are positive constants and U(x, y, z) is a C¹ nonnegative function in the interior of ℝ+3. Hence, we have

dUdt=∂U∂xdxdt+∂U∂ydydt+∂U∂zdzdt,

with ∂U∂x=b1x(b1-1 )yb2zb3, ∂U∂y=b2xb1y(b2-1 )zb3, ∂U∂z=b3xb1yb2z(b3-1 ). Therefore, Ψ(x,y,z)=(dUdt)U(x,y,z)=b1xdxdt+b2ydydt+b3zdzdt,                           =b1[r(1-x)(11+ky)-a1y-a2z]                    +b2[e1a1x-a3z-d1]+b3[e2a2x+e3a3y-d2]. Since there are no periodic solutions in the boundary planes of the system (1), it follows that system (1) is uniformly persistent provided that Ψ(E_i) > 0 for each i = 1, 2, 3. We have Ψ(E1)=b2[e1a1-d1]+b3[e2a2-d2 ], Ψ(E2)=b2[e1a1x̄-a3z̄-d1 ], Ψ(E3)=b3[e2a2x^+e3a3y^-d2 ]. Consequently, conditions (7a)–(7d) satisfy Ψ(E_i) > 0 for each i = 1, 2, 3.

Here, we use Lyapunov functions to investigate the global stability of equilibria.

Theorem (4): The axial equilibrium point E₁ = (1, 0, 0) is globally asymptotically stable provided that the following sufficient condition holds: (8)e2<e1e3<d2a2.

Proof. Consider the following positive-definite, real-valued function around E₁: V1=c1[x-1-ln x]+c2y+c3z . Then, we have: dV1dt=-c1r1+ky(x-1)2-(c1-c2e1)a1xy-(c1-c3e2)a2xz-(c2-c3e3)a3yz-(c2d1-c1a1)y-(c3d2-c1a2)z. Then, by choosing c₁ = e₁, c₂ = 1 and c3=1e3, we obtain dV1dt=-e1r1+ky(x-1)2-(e1-e2e3)a2xz-(d1-e1a1)y           -(d2e3-e1a2 )z. Clearly, under condition (8), dV1dt is negative definite. Moreover, since V₁ is radially unbounded, the axial equilibrium point E₁ = (1, 0, 0) is globally asymptotically stable.

Theorem (5): The specialist-free equilibrium point E2=(x̄,0,z̄) is globally asymptotically stable provided that the following sufficient conditions hold: (9a)a1x̄<d1e1+e3a3e2z̄, (9b)e1e3<e2.

Proof. Consider the following positive-definite, real-valued function around E₂: V2=c1[x-x̄-x̄lnxx̄]+c2y+c3[z-z̄-z̄lnzz̄ ]. Then, we have: dV2dt=-c1r(x-x̄)2-a2[c1-c3e2](x-x̄)(z-z̄)                              -(c1-c2e1)a1xy-(c2d1-c1a1x̄+c3e3a3z̄)y-(c2-c3e3)a3yz.  Then, by choosing c₁ = 1, c2=1e1 and c3=1e2, we have: dV2dt=-r(x-x̄)2-(d1e1+e3a3e2z̄-a1x̄)y-(1e1-e3e2 )a3yz. Obviously, under conditions (9a)–(9b), dV2dt is negative definite. Moreover, since V₂ is radially unbounded, the specialist-free equilibrium point E2=(x̄,0,z̄) is globally asymptotically stable.

Theorem (6): The superpredator-free equilibrium point E3=(x^,y^,0)is globally asymptotically stable, provided that, in addition to condition (9b), the following sufficient conditions hold: (10a)a2x^+a3e1y^<d2e2, (10b)(x-x^)(y-y^)>0.

Proof. Consider the following positive-definite, real-valued function around E₃: V3=c1[x-x^-x^lnxx^]+c2[y-y^-y^lnyy^]+c3z. Then, we have dV3dt=-c1r(x-x^)2R-c1rk(1-x^)(x-x^)(y-y^)RR^                                               -(c1-c2e1)a1(x-x^)(y-y^)-(c1-c3e2)a2xz                                               -(c2-c3e3)a3yz-(c3d2-c1a2x^-c2a3y^)z,  where R = (1 + ky) and R^=(1+ky^). Then, by choosing c₁ = 1, c2=1e1 and c3=1e2, we obtain: dV3dt≤-r(x-x^)2R-rk(1-x^)(x-x^)(y-y^)RR^          -(1e1-e3e2)a3yz-(d2e2-a2x^-a3e1y^ )z. Under conditions (10a)–(10b), dV3dt is negative definite. Since V₃ is radially unbounded, the superpredator-free equilibrium point E3=(x^,y^,0) is globally asymptotically stable.

Theorem (7): The coexistence equilibrium point E4=(x*,y*,z*) is globally asymptotically stable provided that, in addition to condition (9b), the following sufficient condition holds: (11)(x-x*)(y-y*)>0.

Proof. Consider the following positive-definite, real-valued function around E₄: V4=c1[x-x*-x*lnxx*]+c2[y-y*-y*lnyy*]                                               + c3[z-z*-z*lnzz* ]. Then, we have: dV4dt=-c1r(x-x*)2R-c1rk(1-x*)(x-x*)(y-y*)RR*               -(c1-c2e1)a1(x-x*)(y-y*)                -(c2-c3e3)a3(y-y*)(z-z*)                     -(c1-c3e2)a2(x-x*)(z-z*),  where R^* = 1 + ky^*. Then, by choosing c₁ = e₂, c₂ = e₃, and c₃ = 1, we have dV4dt=-e2r(x-x*)2R-e2rk(1-x*)(x-x*)(y-y*)RR*             -(e2-e1e3)a1(x-x*)(y-y* ). Then, under condition (11) with (9b), dV4dt is negative definite. Since V₄ is radially unbounded, the coexistence equilibrium point E4(x*,y*,z*) is globally asymptotically stable.

In this section, we investigate the effect of varying parameter values on the dynamics of the system (1) using local bifurcation analysis and the Sotomayor theorem [29]. First, we rewrite system (1) in the vector form as follows: dXdt=F(X),X=(x,y,z)Tand F=(xf1,yf2,zf3 )T. The second derivative of F with respect to X can be written as: (12)D2F(U,U)=(-2ru121+ky-2kr(1-2x)u1u2(1+ky)2-2a1u1u2-2a2u1u3+2k2rx(1-x)u22(1+ky)32e1a1u1u2-2a3u2u32e2a2u1u3+2e3a3u2u3), where U=(u1,u2,u3)T is a general vector.

Theorem (8): Assume that d2=e2a2 (≡d2*). Then system (1) undergoes a transcritical bifurcation at the axial equilibrium point E₁, but neither a saddle node nor a pitchfork bifurcation can occur if (13)e1a1<d1.

Proof. The Jacobian matrix of system (1) at the axial equilibrium point E₁ with d2=e2a2(≡d2*) can be written as: J1=J(E1,d2*)=[-r      -a1      -a20      e1a1-d1      00            0               0 ].

J₁ has eigenvalues λ11*=-r<0, λ12*=e1a1-d1<0 under condition (13), and λ13*=0. Hence, the necessary but not sufficient condition for a bifurcation is satisfied, and E₁ is a non-hyperbolic point.

Let Φ1=(v11,v12,v13)T be the eigenvector of J₁ corresponding to the eigenvalue λ13*=0. Straightforward computation gives Φ1=(β1v13,0,v13)T, where v₁₃ represents any nonzero real number and β1=-a2r< 0.

Let Ψ1=(ψ11,ψ12,ψ13)T be the eigenvector of J1T corresponding to the eigenvalue λ13*=0. Direct calculation shows that Ψ1=(0,0,ψ13)T, where ψ₁₃ is any nonzero real number. Because ∂F∂d2=Fd2=(0,0,-z)T, we obtain that Fd2(E1,d2*)=(0,0,0)T, which yields: Ψ1T[Fd2(E1,d2*)]= 0. By Sotomayor's theorem, system (1) at E₁ with d2=d2* does not experience a saddle-node bifurcation. Moreover, we have Ψ1T[DFd2(E1,d2*)Φ1]=-v13Ψ13≠ 0, where DF_d2 represents the derivative of F_d2 with respect to X. Applying equation (12) at (E1,d2*) with the eigenvector Φ₁, we obtain that Ψ1T[D2F(E1,d2*)(Φ1,Φ1)]=2e2a2β1v132ψ13≠ 0. By Sotomayor's theorem, system (1) near the equilibrium point E₁ with d2=d2* undergoes a transcritical bifurcation, but a pitchfork bifurcation cannot occur.

Theorem (9): System (1) at the specialist-free equilibrium point E₂ undergoes a transcritical bifurcation when d1=e1a1x̄+a3z̄ (≡d1*), but neither a saddle-node nor a pitchfork bifurcation can occur, provided (14)e1e3a1a2+re3a3≠e2a2 [kr(1-x̄)+a1].

Proof. The Jacobian matrix of system (1) at the specialist-free equilibrium point E₂ with d1=d1* takes the form: J2=J(E2,d1*)=[-rx̄-krx̄(1-x̄)-a1x̄-a2x̄000e2a2z̄e3a3z̄0 ].

J₂ has eigenvalues λ21* =-rx̄2+12(rx̄)2-4e2a22x̄z ̄,λ23* =-rx̄2-12(rx̄)2-4e2a22x̄z̄, while λ22*=0; hence, the necessary but not sufficient condition for a bifurcation is satisfied, and E₂ is a non-hyperbolic point.

Let Φ2=(v21,v22,v23)T be the eigenvector of J₂ corresponding to the zero eigenvalue. Straightforward computation gives Φ2=(α1v22,v22,α2v22)T, where v₂₂ represents any nonzero real number, with α1=-e3a3e2a2 and α2=re3a3-[kr(1-x̄)+a1]e2a2e2a22. Let Ψ2=(ψ21,ψ22,ψ23)T be the eigenvector of J2T corresponding to the zero eigenvalue. Direct calculation shows that Ψ2=(0,ψ22,0)T, where ψ₂₂ is any nonzero real number. Because ∂F∂d1=Fd1=(0,-y,0)T, we obtain that Fd1(E2,d1*)=(0,0,0)T, which yields Ψ2T[Fd1(E2,d1*)]= 0. Thus, by Sotomayor's theorem, system (1) at E₂ with d1=d1* does not experience a saddle-node bifurcation. Moreover, we have: Ψ2T[DFd1(E2,d1*)Φ2]=-v22ψ22≠ 0, where DF_d1 represents the derivative of F_d1 with respect to X. Using equation (12) at (E2,d1*) with the eigenvector Φ₂, we obtain that Ψ2T[D2F(E2,d1*)(Φ2,Φ2)]=(2e1a1α1-2a3α2 )v222ψ22. It is easy to verify that Ψ2T[D2F(E2,d1*)(Φ2,Φ2)]≠0 due to condition (14). Hence, by Sotomayor's theorem, system (1) near the equilibrium point E₂ with d1=d1* undergoes a transcritical bifurcation but a pitchfork bifurcation cannot occur.

Theorem (10): At the superpredator-free equilibrium point E₃, system (1) undergoes a transcritical bifurcation when d2=e2a2x^+e3a3y^ (≡d2*), but neither a saddle-node nor a pitchfork bifurcation can occur, provided (15)e2a2[kr(1-x^)+a1(1+ky^)2]≠e3(1+ky^)[a3r+e1a1a2(1+ky^)].

Proof. The Jacobian matrix of system (1) at the superpredator-free equilibrium point E₃ with d2=d2* is J3=J(E3,d2*)=[-rx^1+ky^-krx^(1-x^)(1+ky^)2-a1x^-a2x^e1a1y^0-a3y^000 ].

J₃ has eigenvalues λ31* =-rx^2(1+ky^)+ 12(rx^1+ky^)2-4e1a1y^(krx^(1-x^)(1+ky^)2+a1x^ ), λ32*=-rx^2(1+ky^)- 12(rx^1+ky^)2-4e1a1y^(krx^(1-x^)(1+ky^)2+a1x^ ), while λ33*=0; hence, the necessary but not sufficient condition for bifurcation is satisfied, and E₃ is a non-hyperbolic point.