Prey fear of a specialist predator in a tri-trophic food web can eliminate the superpredator

We propose an intraguild predation ecological system consisting of a tri-trophic food web with a fear response for the basal prey and a Lotka–Volterra functional response for predation by both a specialist predator (intraguild prey) and a generalist predator (intraguild predator), which we call the superpredator. We prove the positivity, existence, uniqueness, and boundedness of solutions, determine all equilibrium points, prove global stability, determine local bifurcations, and illustrate our results with numerical simulations. An unexpected outcome of the prey's fear of its specialist predator is the potential eradication of the superpredator.

Introduction

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.

Food-web model formulation

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 $\frac{r}{1+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:

$\begin{array}{ll}\begin{array}{l}\frac{dx}{dt}=rx(1-x)\left(\frac{1}{1+ky}\right)-a_{1}xy-a_{2}xz=x{f}_1(x,y,z),\\ \frac{dy}{dt}=e_1{a}_{1}xy-a_{3}yz-d_{1}y=y{f}_2(x,y,z),\\ \frac{dz}{dt}=e_2{a}_{2}xz+e_3{a}_{3}yz-d_{2}z=z{f}_3(x,y,z).\end{array}& (1)\end{array}$

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:

$\begin{array}{l}x(t)=x(0)e^{{\int }_0^t\left[r(1-x(s))(\frac{1}{1+ky(s)})-a_{1}y(s)-a_{2}z(s)\right]ds}\ge 0\end{array}$

$\begin{array}{l}y(t)=y(0)e^{{\int }_0^t\left[e_1{a}_{1}x(s)-a_{3}z(s)-d_1\right]ds}\ge 0\end{array}$

$\begin{array}{l}z(t)=z(0)e^{{\int }_0^t\left[e_2{a}_{2}x(s)+e_3{a}_{3}y(s)-d_2\right]ds}\ge 0.\end{array}$

From the first equation of system (1):

$\begin{array}{l}\frac{dx}{dt}\le rx(1-x).\end{array}$

Then, it is easy to verify that $x(t)\le \frac{r}{4}$ for all t.

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

$\frac{dW}{dt}\le \frac{r}{4}-M\left[x+y+z\right]$,

where M = min{1, d₁, d₂}. Solving the following linear differential inequality

$\begin{array}{l}\frac{dW}{dt}+MW\le \frac{r}{4},\end{array}$

we obtain that, as t → ∞,

$\begin{array}{l}W(t)\le \frac{r}{4M}.\end{array}$

Existence of equilibrium points

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, $E_2=(\bar{x},0,\bar{z})$, where

$\begin{array}{ll}\bar{x}=\frac{d_2}{e_2{a}_2}\text{ and }\bar{z}=\frac{r(1-\bar{x})}{a_2},& (\mathrm{2a})\end{array}$

exists provided that

$\begin{array}{ll}{d}_2<e_2{a}_2.& (\mathrm{2b})\end{array}$

The superpredator-free equilibrium point, $E_3=(\widehat{x},\widehat{y},0)$, where

$\begin{array}{ll}\widehat{x}=\frac{d_1}{e_1{a}_1}\text{ and }\widehat{y}=-\frac{1}{2k}+\frac{1}{2a_{1}k}\sqrt{a_1^2+4kr\frac{(e_1{a}_1-d_1)}{e_1}},& (\mathrm{3a})\end{array}$

exists, provided that

$\begin{array}{ll}{d}_1<e_1{a}_1.& (\mathrm{3b})\end{array}$

There is also an interior equilibrium point, $E_4=(x^{*},y^{*},z^{*})$, where

$\begin{array}{ll}{y}^{{}^{*}}=\frac{d_2-e_2{a}_2{x}^{{}^{*}}}{e_3{a}_3}\text{ and }{z}^{{}^{*}}=\frac{e_1{a}_1{x}^{{}^{*}}-d_1}{a_3},& (\mathrm{4a})\end{array}$

with x^* a positive root of the following second-order polynomial equation:

$\begin{array}{ll}{\beta }_1{x}^2+{\beta }_{2}x+{\beta }_3=0,& (\mathrm{4b})\end{array}$

where

$\begin{array}{l}{\beta }_1=e_2{a}_1{a_2}^{2}k\left[e_1{e}_3-e_2\right],\end{array}$

$\begin{array}{l}\begin{array}{c}{\beta }_2=e_2{e}_3{a}_1{a}_2{a}_3-r{e_3}^2{a_3}^2+2e_2{a}_1{a}_{2}k{d}_2-e_1{e}_3{a}_1{a}_2{a_3}^2\\ -e_1{e}_3{a}_1{a}_{2}k{d}_2-e_2{e}_3{a_2}^{2}k{d}_1,\end{array}\end{array}$

$\begin{array}{l}{\beta }_3=r{e_3}^2{a_3}^2-e_3{a}_1{a}_3{d}_2-a_{1}k{d_2}^2+e_3{a}_2{a_3}^2{d}_1+e_3{a}_{2}k{d}_1{d}_2.\end{array}$

Therefore, there is a unique interior equilibrium point in the interior of ${ℝ}_{+}^3$ provided that the following conditions hold:

$\begin{array}{ll}\frac{d_1}{e_1{a}_1}<x^{{}^{*}}<\frac{d_2}{e_2{a}_2}& (\mathrm{4c})\end{array}$

$\begin{array}{ll}\begin{array}{c}{\beta }_1>0;{\beta }_3<0\\ OR\\ {\beta }_1<0;{\beta }_3>0\end{array}\}.& (\mathrm{4d})\end{array}$

Persistence

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

$\begin{array}{ll}\begin{array}{l}\frac{dx}{dt}=x\left[\frac{r(1-x)}{1+ky}-a_{1}y\right]=x{f}_1(x,y)\\ \frac{dy}{dt}=y\left[e_1{a}_{1}x-d_1\right]=y{g}_1(x,y).\end{array}& (5)\end{array}$

The second subsystem can be written as follows:

$\begin{array}{ll}\begin{array}{l}\frac{dx}{dt}=x\left[r(1-x)-a_{2}z\right]=x{f}_2(x,z)\\ \frac{dz}{dt}=z\left[e_2{a}_{2}x-d_2\right]=z{g}_2(x,z).\end{array}& (6)\end{array}$

Subsystem (5) has a unique positive equilibrium point in the interior of xy-plane given by $(\widehat{x},\widehat{y})\equiv E_3$, while subsystem (6) has a unique positive equilibrium point in the interior of xz-plane given by $(\bar{x},\bar{z})\equiv E_2$. 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)=\frac{1}{xy}$. Then, we have:

$\begin{array}{l}\Delta =\frac{\partial (B{f}_1)}{\partial x}+\frac{\partial (B{g}_1)}{\partial y}=-\frac{r}{y(1+ky)}.\end{array}$

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)=\frac{1}{xz}$.

Theorem (3): System (1) is uniformly persistent in the interior of ${ℝ}_{+}^3,$ provided that

$\begin{array}{ll}{e}_1{a}_1>d_1,& (\mathrm{7a})\end{array}$

$\begin{array}{ll}{e}_2{a}_2>d_2,& (\mathrm{7b})\end{array}$

$\begin{array}{ll}{e}_1{a}_1\bar{x}>a_3\bar{z}+d_1,& (\mathrm{7c})\end{array}$

$\begin{array}{ll}{e}_2{a}_2\widehat{x}+e_3{a}_3\widehat{y}>d_2.& (\mathrm{7d})\end{array}$

Proof. Consider the function $(x,y,z)=x^{b_1}{y}^{b_2}{z}^{b_3}$, 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

$\frac{dU}{dt}=\frac{\partial U}{\partial x}\frac{dx}{dt}+\frac{\partial U}{\partial y}\frac{dy}{dt}+\frac{\partial U}{\partial z}\frac{dz}{dt}$,

with

$\begin{array}{l}\frac{\partial U}{\partial x}=b_1{x}^{(b_1-1)}{y}^{b_2}{z}^{b_3},\end{array}$

$\begin{array}{l}\frac{\partial U}{\partial y}=b_2{x}^{b_1}{y}^{(b_2-1)}{z}^{b_3},\end{array}$

$\begin{array}{l}\frac{\partial U}{\partial z}=b_3{x}^{b_1}{y}^{b_2}{z}^{(b_3-1)}.\end{array}$

Therefore,

$\begin{array}{l}\begin{array}{l}\text{Ψ}(x,y,z)=\frac{(\frac{dU}{dt})}{U(x,y,z)}=\frac{b_1}{x}\frac{dx}{dt}+\frac{b_2}{y}\frac{dy}{dt}+\frac{b_3}{z}\frac{dz}{dt},\\ =b_1\left[r(1-x)\left(\frac{1}{1+ky}\right)-a_{1}y-a_{2}z\right]\\ +b_2\left[e_1{a}_{1}x-a_{3}z-d_1\right]+b_3\left[e_2{a}_{2}x+e_3{a}_{3}y-d_2\right].\end{array}\end{array}$

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

$\begin{array}{l}\text{Ψ}(E_1)=b_2\left[e_1{a}_1-d_1\right]+b_3\left[e_2{a}_2-d_2\right],\end{array}$

$\begin{array}{l}\text{Ψ}(E_2)=b_2\left[e_1{a}_1\bar{x}-a_3\bar{z}-d_1\right],\end{array}$

$\begin{array}{l}\text{Ψ}(E_3)=b_3\left[e_2{a}_2\widehat{x}+e_3{a}_3\widehat{y}-d_2\right].\end{array}$

Consequently, conditions (7a)–(7d) satisfy Ψ(E_i) > 0 for each i = 1, 2, 3.

Global stability analysis

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:

$\begin{array}{ll}{e}_2<e_1{e}_3<\frac{d_2}{a_2}.& (8)\end{array}$

Proof. Consider the following positive-definite, real-valued function around E₁:

$\begin{array}{l}{V}_1=c_1\left[x-1-lnx\right]+c_{2}y+c_{3}z.\end{array}$

Then, we have:

$\begin{array}{l}\begin{array}{l}\frac{d{V}_1}{dt}=-\frac{c_{1}r}{1+ky}{(x-1)}^2-(c_1-c_2{e}_1)a_{1}xy-(c_1-c_3{e}_2)a_{2}xz\\ -(c_2-c_3{e}_3)a_{3}yz-(c_2{d}_1-c_1{a}_1)y-(c_3{d}_2-c_1{a}_2)z.\end{array}\end{array}$

Then, by choosing c₁ = e₁, c₂ = 1 and $c_3=\frac{1}{e_3}$, we obtain

$\begin{array}{l}\frac{d{V}_1}{dt}=-\frac{e_{1}r}{1+ky}{(x-1)}^2-\left(e_1-\frac{e_2}{e_3}\right)a_{2}xz-(d_1-e_1{a}_1)y\\ -\left(\frac{d_2}{e_3}-e_1{a}_2\right)z.\end{array}$

Clearly, under condition (8), $\frac{d{V}_1}{dt}$ 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 $E_2=(\bar{x},0,\bar{z})$ is globally asymptotically stable provided that the following sufficient conditions hold:

$\begin{array}{ll}{a}_1\bar{x}<\frac{d_1}{e_1}+\frac{e_3{a}_3}{e_2}\bar{z},& (\mathrm{9a})\end{array}$

$\begin{array}{ll}{e}_1{e}_3<e_2.& (\mathrm{9b})\end{array}$

Proof. Consider the following positive-definite, real-valued function around E₂:

$\begin{array}{l}{V}_2=c_1\left[x-\bar{x}-\bar{x}ln\frac{x}{\bar{x}}\right]+c_{2}y+c_3\left[z-\bar{z}-\bar{z}ln\frac{z}{\bar{z}}\right].\end{array}$

Then, we have:

$\begin{array}{l}\begin{array}{l}\frac{d{V}_2}{dt}=-c_{1}r{(x-\bar{x})}^2-a_2\left[c_1-c_3{e}_2\right](x-\bar{x})(z-\bar{z})\\ -(c_1-c_2{e}_1)a_{1}xy\\ -(c_2{d}_1-c_1{a}_1\bar{x}+c_3{e}_3{a}_3\bar{z})y-(c_2-c_3{e}_3)a_{3}yz.\end{array}\end{array}$

Then, by choosing c₁ = 1, $c_2=\frac{1}{e_1}$ and $c_3=\frac{1}{e_2}$, we have:

$\begin{array}{l}\frac{d{V}_2}{dt}=-r{(x-\bar{x})}^2-\left(\frac{d_1}{e_1}+\frac{e_3{a}_3}{e_2}\bar{z}-a_1\bar{x}\right)y\\ -\left(\frac{1}{e_1}-\frac{e_3}{e_2}\right)a_{3}yz.\end{array}$

Obviously, under conditions (9a)–(9b), $\frac{d{V}_2}{dt}$ is negative definite. Moreover, since V₂ is radially unbounded, the specialist-free equilibrium point $E_2=(\bar{x},0,\bar{z})$ is globally asymptotically stable.

Theorem (6): The superpredator-free equilibrium point $E_3=(\widehat{x},\widehat{y},0)$is globally asymptotically stable, provided that, in addition to condition (9b), the following sufficient conditions hold:

$\begin{array}{ll}{a}_2\widehat{x}+\frac{a_3}{e_1}\widehat{y}<\frac{d_2}{e_2},& (\mathrm{10a})\end{array}$

$\begin{array}{ll}(x-\widehat{x})(y-\widehat{y})>0.& (\mathrm{10b})\end{array}$

Proof. Consider the following positive-definite, real-valued function around E₃:

$\begin{array}{l}{V}_3=c_1\left[x-\widehat{x}-\widehat{x}ln\frac{x}{\widehat{x}}\right]+c_2\left[y-\widehat{y}-\widehat{y}ln\frac{y}{\widehat{y}}\right]+c_{3}z.\end{array}$

Then, we have

$\begin{array}{l}\begin{array}{l}\frac{d{V}_3}{dt}=-\frac{c_{1}r{(x-\widehat{x})}^2}{R}-\frac{c_{1}rk(1-\widehat{x})(x-\widehat{x})(y-\widehat{y})}{R\widehat{R}}\\ -(c_1-c_2{e}_1)a_1(x-\widehat{x})(y-\widehat{y})-(c_1-c_3{e}_2)a_{2}xz\\ -(c_2-c_3{e}_3)a_{3}yz-(c_3{d}_2-c_1{a}_2\widehat{x}-c_2{a}_3\widehat{y})z,\\ \end{array}\end{array}$

where R = (1 + ky) and $\widehat{R}=(1+k\widehat{y})$. Then, by choosing c₁ = 1, $c_2=\frac{1}{e_1}$ and $c_3=\frac{1}{e_2}$, we obtain:

$\begin{array}{l}\frac{d{V}_3}{dt}\le -\frac{r{(x-\widehat{x})}^2}{R}-\frac{rk(1-\widehat{x})(x-\widehat{x})(y-\widehat{y})}{R\widehat{R}}\\ -\left(\frac{1}{e_1}-\frac{e_3}{e_2}\right)a_{3}yz-\left(\frac{d_2}{e_2}-a_2\widehat{x}-\frac{a_3}{e_1}\widehat{y}\right)z.\end{array}$

Under conditions (10a)–(10b), $\frac{d{V}_3}{dt}$ is negative definite. Since V₃ is radially unbounded, the superpredator-free equilibrium point $E_3=(\widehat{x},\widehat{y},0)$ is globally asymptotically stable.

Theorem (7): The coexistence equilibrium point $E_4=(x^{*},y^{*},z^{*})$ is globally asymptotically stable provided that, in addition to condition (9b), the following sufficient condition holds:

$\begin{array}{ll}\left(x-x^{{}^{*}}\right)\left(y-y^{{}^{*}}\right)>0.& (11)\end{array}$

Proof. Consider the following positive-definite, real-valued function around E₄:

$\begin{array}{l}{V}_4=c_1\left[x-x^{{}^{*}}-x^{{}^{*}}ln\frac{x}{x^{{}^{*}}}\right]+c_2\left[y-y^{{}^{*}}-y^{{}^{*}}ln\frac{y}{y^{{}^{*}}}\right]\\ +c_3\left[z-z^{{}^{*}}-z^{{}^{*}}ln\frac{z}{z^{{}^{*}}}\right].\end{array}$

Then, we have:

$\begin{array}{l}\begin{array}{l}\frac{d{V}_4}{dt}=-\frac{c_{1}r{\left(x-x^{{}^{*}}\right)}^2}{R}-\frac{c_{1}rk\left(1-x^{{}^{*}}\right)\left(x-x^{{}^{*}}\right)\left(y-y^{{}^{*}}\right)}{R{R}^{{}^{*}}}\\ \begin{array}{c}-\left(c_1-c_2{e}_1\right)a_1\left(x-x^{{}^{*}}\right)\left(y-y^{{}^{*}}\right)\\ -\left(c_2-c_3{e}_3\right)a_3\left(y-y^{{}^{*}}\right)\left(z-z^{{}^{*}}\right)\\ -\left(c_1-c_3{e}_2\right)a_2\left(x-x^{{}^{*}}\right)\left(z-z^{{}^{*}}\right),\\ \end{array}\end{array}\end{array}$

where R^* = 1 + ky^*. Then, by choosing c₁ = e₂, c₂ = e₃, and c₃ = 1, we have

$\begin{array}{l}\frac{d{V}_4}{dt}=-\frac{e_{2}r{\left(x-x^{{}^{*}}\right)}^2}{R}-\frac{e_{2}rk\left(1-x^{{}^{*}}\right)\left(x-x^{{}^{*}}\right)\left(y-y^{{}^{*}}\right)}{R{R}^{{}^{*}}}\\ -\left(e_2-e_1{e}_3\right)a_1\left(x-x^{{}^{*}}\right)\left(y-y^{{}^{*}}\right).\end{array}$

Then, under condition (11) with (9b), $\frac{d{V}_4}{dt}$ is negative definite. Since V₄ is radially unbounded, the coexistence equilibrium point $E_4(x^{*},y^{*},z^{*})$ is globally asymptotically stable.

Bifurcation analysis

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:

$\begin{array}{l}\frac{dX}{dt}=F(X),X={\left(x,y,z\right)}^T\text{and }F={\left(x{f}_1,y{f}_2,z{f}_3\right)}^T.\end{array}$

The second derivative of F with respect to X can be written as:

$\begin{array}{l}{D}^{2}F(U,U)=\\ \left(\begin{array}{c}-\frac{2r{u_1}^2}{1+ky}-\frac{2kr(1-2x)u_1{u}_2}{{(1+ky)}^2}-2a_1{u}_1{u}_2-2a_2{u}_1{u}_3+\frac{2k^{2}rx(1-x){u_2}^2}{{(1+ky)}^3}\\ 2e_1{a}_1{u}_1{u}_2-2a_3{u}_2{u}_3\\ 2e_2{a}_2{u}_1{u}_3+2e_3{a}_3{u}_2{u}_3\end{array}\right),& (12)\end{array}$

where $U={(u_1,u_2,u_3)}^T$ is a general vector.

Theorem (8): Assume that $d_2=e_2{a}_2(\equiv {d_2}^{*})$. 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

$\begin{array}{ll}{e}_1{a}_1<d_1.& (13)\end{array}$

Proof. The Jacobian matrix of system (1) at the axial equilibrium point E₁ with $d_2=e_2{a}_2(\equiv {d_2}^{*})$ can be written as:

$\begin{array}{l}{J}_1=J(E_1,{d_2}^{*})=\left[\begin{array}{c}-r-a_1-a_2\\ 0e_1{a}_1-d_{1}0\\ 000\\ \end{array}\right].\end{array}$

J₁ has eigenvalues ${{\text{λ}}_{11}}^{*}=-r<0$, ${{\text{λ}}_{12}}^{*}=e_1{a}_1-d_1<0$ under condition (13), and ${{\text{λ}}_{13}}^{*}=0$. Hence, the necessary but not sufficient condition for a bifurcation is satisfied, and E₁ is a non-hyperbolic point.

Let ${\Phi }_1={(v_{11},v_{12},v_{13})}^T$ be the eigenvector of J₁ corresponding to the eigenvalue ${{\text{λ}}_{13}}^{*}=0$. Straightforward computation gives ${\Phi }_1={({\beta }_1{v}_{13},0,v_{13})}^T$, where v₁₃ represents any nonzero real number and ${\beta }_1=-\frac{a_2}{r}<0$.

Let ${\text{Ψ}}_1={({\psi }_{11},{\psi }_{12},{\psi }_{13})}^T$ be the eigenvector of ${J_1}^T$ corresponding to the eigenvalue ${{\text{λ}}_{13}}^{*}=0$. Direct calculation shows that ${\text{Ψ}}_1={(0,0,{\psi }_{13})}^T$, where ψ₁₃ is any nonzero real number. Because $\frac{\partial F}{\partial d_2}=F_{d_2}={(0,0,-z)}^T$, we obtain that $F_{d_2}(E_1,{d_2}^{*})={(0,0,0)}^T$, which yields:

$\begin{array}{l}{{\text{Ψ}}_1}^T\left[F_{d_2}\left(E_1,{d_2}^{{}^{*}}\right)\right]=0.\end{array}$

By Sotomayor's theorem, system (1) at E₁ with $d_2={d_2}^{*}$ does not experience a saddle-node bifurcation. Moreover, we have

$\begin{array}{l}{{\text{Ψ}}_1}^T\left[D{F}_{d_2}\left(E_1,{d_2}^{{}^{*}}\right){\Phi }_1\right]=-v_{13}{\text{Ψ}}_{13}\ne 0,\end{array}$

where DF_d2 represents the derivative of F_d2 with respect to X. Applying equation (12) at $(E_1,{d_2}^{*})$ with the eigenvector Φ₁, we obtain that

$\begin{array}{l}{{\text{Ψ}}_1}^T\left[D^{2}F\left(E_1,{d_2}^{{}^{*}}\right)({\Phi }_1,{\Phi }_1)\right]=2e_2{a}_2{\beta }_1{v_{13}}^2{\psi }_{13}\ne 0.\end{array}$

By Sotomayor's theorem, system (1) near the equilibrium point E₁ with $d_2={d_2}^{*}$ 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 $d_1=e_1{a}_1\bar{x}+a_3\bar{z}(\equiv {d_1}^{*})$, but neither a saddle-node nor a pitchfork bifurcation can occur, provided

$\begin{array}{ll}{e_1{e}_3{a}_{1}a}_2+r{e}_3{a}_3\ne e_2{a}_2\left[kr(1-\bar{x})+a_1\right].& (14)\end{array}$

Proof. The Jacobian matrix of system (1) at the specialist-free equilibrium point E₂ with $d_1={d_1}^{*}$ takes the form:

$\begin{array}{l}{J}_2=J(E_2,{d_1}^{{}^{*}})=\left[\begin{array}{ccc}-r\bar{x}& -kr\bar{x}(1-\bar{x})-a_1\bar{x}& -a_2\bar{x}\\ 0& 0& 0\\ e_2{a}_2\bar{z}& e_3{a}_3\bar{z}& 0\\ \end{array}\right].\end{array}$

J₂ has eigenvalues

$\begin{array}{l}{{\text{λ}}_{21}}^{{}^{*}}=-\frac{r\bar{x}}{2}+\frac{1}{2}\sqrt{{(r\bar{x})}^2-4e_2{a_2}^2\bar{x}\bar{z}},\\ {{\text{λ}}_{23}}^{{}^{*}}=-\frac{r\bar{x}}{2}-\frac{1}{2}\sqrt{{(r\bar{x})}^2-4e_2{a_2}^2\bar{x}\bar{z}},\end{array}$

while ${{\text{λ}}_{22}}^{*}=0$; hence, the necessary but not sufficient condition for a bifurcation is satisfied, and E₂ is a non-hyperbolic point.

Let ${\Phi }_2={(v_{21},v_{22},v_{23})}^T$ be the eigenvector of J₂ corresponding to the zero eigenvalue. Straightforward computation gives ${\Phi }_2={({\alpha }_1{v}_{22},v_{22},{\alpha }_2{v}_{22})}^T$, where v₂₂ represents any nonzero real number, with

$\begin{array}{l}{\alpha }_1=-\frac{e_3{a}_3}{e_2{a}_2}\text{ and }{\alpha }_2=\frac{r{e}_3{a}_3-\left[kr(1-\bar{x})+a_1\right]e_2{a}_2}{e_2{a_2}^2}.\end{array}$

Let ${\text{Ψ}}_2={({\psi }_{21},{\psi }_{22},{\psi }_{23})}^T$ be the eigenvector of ${J_2}^T$ corresponding to the zero eigenvalue. Direct calculation shows that ${\text{Ψ}}_2={(0,{\psi }_{22},0)}^T$, where ψ₂₂ is any nonzero real number. Because $\frac{\partial F}{\partial d_1}=F_{d_1}={(0,-y,0)}^T$, we obtain that $F_{d_1}(E_2,{d_1}^{*})={(0,0,0)}^T$, which yields

$\begin{array}{l}{{\text{Ψ}}_2}^T\left[F_{d_1}\left(E_2,{d_1}^{{}^{*}}\right)\right]=0.\end{array}$

Thus, by Sotomayor's theorem, system (1) at E₂ with $d_1={d_1}^{*}$ does not experience a saddle-node bifurcation. Moreover, we have:

$\begin{array}{l}{{\text{Ψ}}_2}^T\left[D{F}_{d_1}\left(E_2,{d_1}^{{}^{*}}\right){\Phi }_2\right]=-v_{22}{\psi }_{22}\ne 0,\end{array}$

where DF_d1 represents the derivative of F_d1 with respect to X. Using equation (12) at $(E_2,{d_1}^{*})$ with the eigenvector Φ₂, we obtain that

$\begin{array}{l}{{\text{Ψ}}_2}^T\left[D^{2}F\left(E_2,{d_1}^{{}^{*}}\right)({\Phi }_2,{\Phi }_2)\right]=(2e_1{a}_1{\alpha }_1-2a_3{\alpha }_2){v_{22}}^2{\psi }_{22}.\end{array}$

It is easy to verify that ${{\text{Ψ}}_2}^T\left[D^{2}F(E_2,{d_1}^{*})({\Phi }_2,{\Phi }_2)\right]\ne 0$ due to condition (14). Hence, by Sotomayor's theorem, system (1) near the equilibrium point E₂ with $d_1={d_1}^{*}$ 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 $d_2=e_2{a}_2\widehat{x}+e_3{a}_3\widehat{y}(\equiv {d_2}^{*})$, but neither a saddle-node nor a pitchfork bifurcation can occur, provided

$\begin{array}{l}{e}_2{a}_2\left[kr(1-\widehat{x})+a_1{(1+k\widehat{y})}^2\right]\ne \\ e_3(1+k\widehat{y})\left[a_{3}r+e_1{a}_1{a}_2(1+k\widehat{y})\right].& (15)\end{array}$

Proof. The Jacobian matrix of system (1) at the superpredator-free equilibrium point E₃ with $d_2={d_2}^{*}$ is

$\begin{array}{l}{J}_3=J\left(E_3,{d_2}^{{}^{*}}\right)=\left[\begin{array}{ccc}-\frac{r\widehat{x}}{1+k\widehat{y}}& -\frac{kr\widehat{x}(1-\widehat{x})}{{(1+k\widehat{y})}^2}-a_1\widehat{x}& -a_2\widehat{x}\\ e_1{a}_1\widehat{y}& 0& -a_3\widehat{y}\\ 0& 0& 0\\ \end{array}\right].\end{array}$

J₃ has eigenvalues

$\begin{array}{l}{{\text{λ}}_{31}}^{{}^{*}}=-\frac{r\widehat{x}}{2(1+k\widehat{y})}\\ +\frac{1}{2}\sqrt{{\left(\frac{r\widehat{x}}{1+k\widehat{y}}\right)}^2-4e_1{a}_1\widehat{y}\left(\frac{kr\widehat{x}(1-\widehat{x})}{{(1+k\widehat{y})}^2}+a_1\widehat{x}\right)},\end{array}$

$\begin{array}{l}{{\text{λ}}_{32}}^{{}^{*}}=-\frac{r\widehat{x}}{2(1+k\widehat{y})}\\ -\frac{1}{2}\sqrt{{\left(\frac{r\widehat{x}}{1+k\widehat{y}}\right)}^2-4e_1{a}_1\widehat{y}\left(\frac{kr\widehat{x}(1-\widehat{x})}{{(1+k\widehat{y})}^2}+a_1\widehat{x}\right)},\end{array}$

while ${{\text{λ}}_{33}}^{*}=0$; hence, the necessary but not sufficient condition for bifurcation is satisfied, and E₃ is a non-hyperbolic point.

Let ${\Phi }_3={(v_{31},v_{32},v_{33})}^T$ be the eigenvector of J₃ corresponding to the zero eigenvalue. Straightforward computation gives ${\Phi }_3={(s_1{v}_{33},s_2{v}_{33},v_{33})}^T$, where v₃₃ represents any nonzero real number, while

$\begin{array}{l}{s}_1=\frac{a_3}{e_1{a}_1}\text{ and }{s}_2=-\left[\frac{r{a}_3(1+k\widehat{y})+e_1{a}_1{a}_2{(1+k\widehat{y})}^2}{e_1{a}_1\left[kr(1-\widehat{x})+a_1{(1+k\widehat{y})}^2\right]}\right].\end{array}$

Let ${\text{Ψ}}_3={({\psi }_{31},{\psi }_{32},{\psi }_{33})}^T$ be the eigenvector of ${J_3}^T$ corresponding to the zero eigenvalue. Direct calculation shows that ${\text{Ψ}}_3={(0,0,{\psi }_{33})}^T$, where ψ₃₃ is any nonzero real number. Because $\frac{\partial F}{\partial d_2}=F_{d_2}={(0,0,-z)}^T$, we obtain $F_{d_2}(E_3,{d_2}^{*})={(0,0,0)}^T$, which yields

$\begin{array}{l}{{\text{Ψ}}_3}^T\left[F_{d_2}\left(E_3,{d_2}^{{}^{*}}\right)\right]=0.\end{array}$

Thus, by Sotomayor's theorem, system (1) at E₃ with $d_2={d_2}^{*}$ does not experience a saddle-node bifurcation. Moreover, we have:

$\begin{array}{l}{{\text{Ψ}}_3}^T\left[D{F}_{d_2}\left(E_3,d_2\right){\Phi }_3\right]=-v_{33}{\psi }_{33}\ne 0,\end{array}$

where DF_d2 represents the derivative of F_d2 with respect to X. By using equation (12) at $(E_3,{d_2}^{*})$ with the eigenvector Φ₃, we obtain that:

$\begin{array}{l}{{\text{Ψ}}_3}^T\left[D^{2}F\left(E_3,{d_2}^{{}^{*}}\right)({\Phi }_3,{\Phi }_3)\right]=(2e_2{a}_2{s}_1+2e_3{a}_3{s}_2){v_{33}}^2{\psi }_{33}.\end{array}$

Straightforward computation shows that ${{\text{Ψ}}_3}^T\left[D^{2}F(E_3,{d_2}^{*})({\Phi }_3,{\Phi }_3)\right]\ne 0$ due to condition (15). Thus, by Sotomayor's theorem, system (1) near the equilibrium point E₃ with $d_2={d_2}^{*}$ undergoes a transcritical bifurcation, but a pitchfork bifurcation cannot occur.

Theorem (11): System (1) undergoes a saddle-node bifurcation at the coexistence equilibrium point E₄, but neither a transcritical nor a pitchfork bifurcation can occur when e₃ passes through the value ${e_3}^{*}=\frac{b_{12}{b}_{23}{b}_{31}}{(b_{11}{b}_{23}-b_{13}{b}_{21})a_3{z}^{*}}$, provided the following conditions hold:

$\begin{array}{ll}\left(\frac{kr{x}^{{}^{*}}\left(1-x^{{}^{*}}\right)}{{\left(1+k{y}^{{}^{*}}\right)}^2}+a_1{x}^{{}^{*}}\right)e_2{a}_2>\left(\frac{r{x}^{{}^{*}}}{1+k{y}^{{}^{*}}}\right){e_3}^{{}^{*}}{a}_3,& (\mathrm{16a})\end{array}$

$\begin{array}{l}\frac{a_1{a}_2}{a_3}\left(\frac{e_2}{{e_3}^{{}^{*}}}-e_1\right)-\frac{r}{\left(1+k{y}^{{}^{*}}\right)}\times \\ \left[1-\frac{k\left(1-2x^{{}^{*}}\right)e_2{a}_2}{\left(1+k{y}^{{}^{*}}\right){e_3}^{{}^{*}}{a}_3}-\frac{k^2{x}^{{}^{*}}\left(1-x^{{}^{*}}\right){\left(e_1{a}_1\right)}^2}{{\left(1+k{y}^{{}^{*}}\right)}^2{a_3}^2}\right]\ne 0.& (\mathrm{16b})\end{array}$

Proof. Direct computation shows that system (1) at the coexistence equilibrium point and ${e_3=e_3}^{*}$ has Jacobian matrix

$\begin{array}{r}{J}_4=J(E_4,{e_3}^{*})=\left[\begin{array}{ccc}-\frac{r{x}^{{}^{*}}}{1+k{y}^{*}}& -\frac{kr{x}^{*}(1-x^{*})}{{(1+k{y}^{*})}^2}-a_1{x}^{*}& -a_2{x}^{*}\\ e_1{a}_1{y}^{*}& 0& -a_3{y}^{*}\\ e_2{a}_2{z}^{*}& {e_3}^{*}{a}_3{z}^{*}& 0\end{array}\right]\\ ={\left[b_{ij}\right]}_{3\times 3}.\end{array}$

Straightforward computation shows that |J₄| = 0, under condition (16a). Hence, J₄ has two eigenvalues with negative real parts, and the third one is λ^* = 0. It follows that E₄ becomes a non-hyperbolic equilibrium point.

Let ${\Phi }_4={(v_{41},v_{42},v_{43})}^T$ be the eigenvector of J₄ corresponding to the eigenvalue λ^* = 0. Straightforward computation gives ${\Phi }_4={(v_{41},{\delta }_1{v}_{41},{\delta }_2{v}_{41})}^T$, where v₄₁ represents any nonzero real number, ${\delta }_1=-\frac{b_{31}}{{b_{32}}^{*}}<0$ and ${\delta }_2=-\frac{b_{21}}{b_{23}}>0$.

Let ${\text{Ψ}}_4={({\psi }_{41},{\psi }_{42},{\psi }_{43})}^T$ be the eigenvector of ${J_4}^T$ corresponding to the eigenvalue λ^* = 0. Direct calculation shows that ${\text{Ψ}}_4={({\psi }_{41},{\mu }_1{\psi }_{41},{\mu }_2{\psi }_{41})}^T$, where ψ₃₃ is any nonzero real number, with ${\mu }_1=-\frac{b_{13}}{b_{23}}<0$ and ${\mu }_2=-\frac{b_{12}}{b_{32}}>0$.

Since $\frac{\partial F}{\partial e_3}=F_{e_3}={(0,0,a_{3}yz)}^T$, we have $F_{e_3}(E_4,{e_3}^{*})={(0,0,a_3{y}^{*}{z}^{*})}^T$, which yields

$\begin{array}{l}{{\text{Ψ}}_4}^T\left[F_{e_3}\left(E_4,{e_3}^{{}^{*}}\right)\right]=a_3{\mu }_2{y}^{{}^{*}}{z}^{{}^{*}}{\psi }_{41}\ne 0.\end{array}$

By Sotomayor's theorem, transcritical and pitchfork bifurcations cannot occur while the first condition of the saddle-node bifurcation is satisfied. Moreover, from equation (12) with $E_4,{e_3}^{*}$ and Φ₄, we obtain that:

$\begin{array}{l}\begin{array}{c}{D}^{2}F\left(E_4,{e_3}^{{}^{*}}\right)\left({\Phi }_4,{\Phi }_4\right)={2v_{41}}^2\times \\ \left(\begin{array}{c}-\frac{r}{1+k{y}^{{}^{*}}}-\frac{kr\left(1-2x^{{}^{*}}\right){\delta }_1}{{\left(1+k{y}^{{}^{*}}\right)}^2}-a_1{\delta }_1-a_2{\delta }_2+\frac{k^{2}r{x}^{{}^{*}}\left(1-x^{{}^{*}}\right){{\delta }_2}^2}{{\left(1+k{y}^{{}^{*}}\right)}^3}\\ e_1{a}_1{\delta }_1-a_3{\delta }_1{\delta }_2\\ e_2{a}_2{\delta }_2+{e_3}^{{}^{*}}{a}_3{\delta }_1{\delta }_2\end{array}\right).\\ \end{array}\end{array}$

Hence,

$\begin{array}{l}\begin{array}{c}{\text{Ψ}}_4^T\left[D^{2}F\left(E_4,{e_3}^{{}^{*}}\right)\left({\Phi }_4,{\Phi }_4\right)\right]={2v_{41}}^2{\psi }_{41}\times \\ [-\frac{r}{1+k{y}^{{}^{*}}}-\frac{kr\left(1-2x^{{}^{*}}\right){\delta }_1}{{\left(1+k{y}^{{}^{*}}\right)}^2}-a_1{\delta }_1-a_2{\delta }_2\\ +\frac{k^{2}r{x}^{{}^{*}}\left(1-x^{{}^{*}}\right){{\delta }_2}^2}{{\left(1+k{y}^{{}^{*}}\right)}^3}+e_1{a}_1{\delta }_1{\mu }_1-a_3{\delta }_1{\delta }_2{\mu }_1+e_2{a}_2{\delta }_2{\mu }_2\\ +{e_3}^{{}^{*}}{a}_3{\delta }_1{\delta }_2{\mu }_2].\end{array}\end{array}$

Further computation shows that

$\begin{array}{l}\begin{array}{l}{{\text{Ψ}}_4}^T\left[D^{2}F\left(E_4,{e_3}^{{}^{*}}\right)\left({\Phi }_4,{\Phi }_4\right)\right]={2v_{41}}^2{\psi }_{41}[\frac{a_1{a}_2}{a_3}\left(\frac{e_2}{{e_3}^{{}^{*}}}-e_1\right)\\ -\frac{r}{\left(1+k{y}^{{}^{*}}\right)}\left[1-\frac{k\left(1-2x^{{}^{*}}\right)e_2{a}_2}{\left(1+k{y}^{{}^{*}}\right){e_3}^{{}^{*}}{a}_3}-\frac{k^2{x}^{{}^{*}}\left(1-x^{{}^{*}}\right){\left(e_1{a}_1\right)}^2}{{\left(1+k{y}^{{}^{*}}\right)}^2{a_3}^2}\right]].\end{array}\end{array}$

Using condition (16b), we have ${{\text{Ψ}}_4}^T\left[D^{2}F(E_4,{e_3}^{*})({\Phi }_4,{\Phi }_4)\right]\ne 0$. Hence, system (1) has a saddle-node bifurcation at E₄ when $e_3={e_3}^{*}$.

Theorem (12): System (2) undergoes a Hopf bifurcation around E₄ when $e_3={e_3}^{**}$, if and only if the following conditions hold:

$\begin{array}{ll}\left(\frac{kr{x}^{*}\left(1-x^{*}\right)}{{\left(1+k{y}^{*}\right)}^2}+a_1{x}^{*}\right)e_2{a}_2<\left(\frac{r{x}^{*}}{1+k{y}^{*}}\right)e_3{a}_3,& (\mathrm{17a})\end{array}$

$\begin{array}{ll}{\left[A({e_3}^{{}^{**}})B({e_3}^{{}^{**}})\right]}^{{}^{\prime }}<C^{{}^{\prime }}({e_3}^{{}^{**}}),& (\mathrm{17b})\end{array}$

where

$\begin{array}{l}A=-b_{11}>0,\end{array}$

$\begin{array}{l}B=-b_{12}{b}_{21}-b_{13}{b}_{31}-b_{23}{b}_{32}^{*}>0,\end{array}$

$\begin{array}{l}C=b_{23}(b_{11}{b}_{32}^{*}-b_{12}{b}_{31})-b_{13}{b}_{21}{b}_{32}^{*},\end{array}$

with $b_{32}^{*}={e_3}^{**}{a}_3{z}^{*}$ and,

$\begin{array}{r}{e_3}^{**}=\frac{1}{e_1{a}_1{a}_2{a}_3{y}^{{}^{*}}{z}^{*}}[\left[\frac{kr\left(1-x^{*}\right)}{{\left(1+k{y}^{*}\right)}^2}-a_1\right]\times \\ \left(e_1{a}_1{y}^{*}\left(\frac{r{x}^{*}}{1+k{y}^{*}}\right)+e_2{a}_2{a}_3{y}^{*}{z}^{*}\right)+e_2{a_2}^2{z}^{*}\left(\frac{r{x}^{*}}{1+k{y}^{*}}\right)].\end{array}$

Proof. Consider the Jacobian matrix J₄ given in Theorem (11), where $e_3={e_3}^{**}$. It is simple to determine that the characteristic equation is

$\begin{array}{ll}{\text{λ}}^3+A{\text{λ}}^2+B\text{λ}+C=0.& (18)\end{array}$

Obviously, condition (17a) guarantees that C > 0. A straightforward computation shows that AB = C when $e_3={e_3}^{**}$. Hence, the characteristic equation (18) becomes

$\begin{array}{ll}\left({\text{λ}}^2+B\right)(\text{λ}+A)=0.& (19)\end{array}$

Consequently, we obtain that λ₁ = −A and ${\text{λ}}_{2,3}=\pm i\sqrt{B}=\pm i{\delta }_2({e_3}^{**})$. Therefore, the Jacobian matrix has one negative real eigenvalue and two pure imaginary complex conjugates when $e_3={e_3}^{**}$. As a result, the first criterion for having a Hopf bifurcation is met.

Moreover, where e₃ belongs to the neighborhood of ${e_3}^{**}$, the eigenvalues become λ_2,3 = δ₁(e₃) ± iδ₂(e₃).

Next, we substitute δ₁(e₃) + iδ₂(e₃) into equation (18) and take the derivative of the resulting equation with respect to e₃. Equating their real and imaginary parts, we obtain that

$\begin{array}{ll}\begin{array}{l}{H}_1(e_3){\delta }_1^{{}^{\prime }}(e_3)-H_2(e_3){\delta }_2^{{}^{\prime }}(e_3)=-H_3(e_3),\\ H_2(e_3){\delta }_1^{{}^{\prime }}(e_3)+H_1(e_3){\delta }_2^{{}^{\prime }}(e_3)=-H_4(e_3),\end{array}& (20)\end{array}$

where

$\begin{array}{l}{H}_1(e_3)=3{\left[{\delta }_1(e_3)\right]}^2+2A(e_3){\delta }_1(e_3)-3{\left[{\delta }_2(e_3)\right]}^2+B(e_3).\end{array}$

$\begin{array}{l}{H}_2(e_3)=6{\delta }_1(e_3){\delta }_2(e_3)+2A(e_3){\delta }_2(e_3).\end{array}$

$\begin{array}{r}{H}_3(e_3)=A^{{}^{\prime }}{(e_3)\left[{\delta }_1(e_3)\right]}^2-A^{{}^{\prime }}{(e_3)\left[{\delta }_2(e_3)\right]}^2+B^{{}^{\prime }}(e_3){\delta }_1(e_3)\\ +C^{{}^{\prime }}(e_3).\end{array}$

$\begin{array}{l}{H}_4(e_3)=2A^{{}^{\prime }}(e_3){{\delta }_1(e_3)\delta }_2(e_3)+B^{{}^{\prime }}(e_3){\delta }_2(e_3).\end{array}$

Solving system (20) for ${\delta }_1^{{}^{\prime }}(e_3)$ gives

$\begin{array}{l}{\delta }_1^{{}^{\prime }}(e_3)=-\frac{H_1(e_3)H_3(e_3)+H_2(e_3)H_4(e_3)}{{\left[H_2(e_3)\right]}^2+{(e_3)\left[H_2(e_3)\right]}^2}.\end{array}$

The result follows if and only if ${\delta }_1^{{}^{\prime }}(e_3)\ne 0$ or, equivalently, H₁(e₃)H₃(e₃) + H₂(e₃)H₄(e₃) ≠ 0 when $e_3={e_3}^{**}$. Note that we have the following:

$\begin{array}{l}{\delta }_1\left({e_3}^{{}^{**}}\right)=0\text{ and }{\delta }_2\left({e_3}^{{}^{**}}\right)=\sqrt{B\left({e_3}^{{}^{**}}\right)}.\end{array}$

$\begin{array}{l}{H}_1\left({e_3}^{{}^{**}}\right)=-2B\left({e_3}^{{}^{**}}\right).\end{array}$

$\begin{array}{l}{H}_2\left({e_3}^{{}^{**}}\right)=2A\left({e_3}^{{}^{**}}\right)\sqrt{B\left({e_3}^{{}^{**}}\right)}.\end{array}$

$\begin{array}{l}{H}_3\left({e_3}^{{}^{**}}\right)=-A^{{}^{\prime }}\left({e_3}^{{}^{**}}\right)\left[B\left({e_3}^{{}^{**}}\right)\right]+C^{{}^{\prime }}\left({e_3}^{{}^{**}}\right).\end{array}$

$\begin{array}{l}{H}_4\left({e_3}^{{}^{**}}\right)=B^{{}^{\prime }}\left({e_3}^{{}^{**}}\right)\sqrt{B\left({e_3}^{{}^{**}}\right)}.\end{array}$

Consequently,

$\begin{array}{l}\begin{array}{l}{H}_1\left({e_3}^{{}^{**}}\right)H_3\left({e_3}^{{}^{**}}\right)+H_2\left({e_3}^{{}^{**}}\right)H_4\left({e_3}^{{}^{**}}\right)\\ =2A^{{}^{\prime }}\left({e_3}^{{}^{**}}\right)B^2\left({e_3}^{{}^{**}}\right)-2B\left({e_3}^{{}^{**}}\right)\\ C^{{}^{\prime }}\left({e_3}^{{}^{**}}\right)+2A\left({e_3}^{{}^{**}}\right)B^{{}^{\prime }}\left({e_3}^{{}^{**}}\right)B\left({e_3}^{{}^{**}}\right)\\ =-2B\left({e_3}^{{}^{**}}\right)[C^{{}^{\prime }}\left({e_3}^{{}^{**}}\right)-(A^{{}^{\prime }}\left({e_3}^{{}^{**}}\right)\\ B\left({e_3}^{{}^{**}}\right)+A\left({e_3}^{{}^{**}}\right)B^{{}^{\prime }}\left({e_3}^{{}^{**}}\right)].\end{array}\end{array}$

Condition (17b) ensures that $H_1({e_3}^{**})H_3({e_3}^{**})+H_2({e_3}^{**})H_4({e_3}^{**})\ne 0$.

Hence, the system has a Hopf bifurcation because ${\delta }_1^{{}^{\prime }}(e_3)>0$ under the conditions (17a)–(17b).

Numerical simulations

To illustrate the global dynamics of the system and confirm our analytical findings, we numerically simulated a hypothetical set of parameter values. Consider the following set of parameters:

$\begin{array}{ll}\begin{array}{l}r=1,k=1,a_1=1,a_2=0.5,a_3=0.75,e_1=0.75\\ e_2=0.25,e_3=0.25,d_1=0.05,d_2=0.1.\end{array}& (21)\end{array}$

Using these data, system (1) asymptotically approaches the coexistence equilibrium point E₄ = (0.5, 0.2, 0.43), as shown in Figure 1.

FIGURE 1

Figure 1. Globally asymptotically stable coexistence equilibrium point using data (17) and multiple initial conditions. (A) 3D-Phase portrait of the system (1). (B) Time series for trajectories of x. (C) Time series for trajectories of y. (D) Time series for trajectories of z.

Next, we varied specific parameters to understand the effect of that parameter on the dynamical behavior of the system. For the values of the parameter r satisfying r ≥ 1.84 and other parameters as in equation (17), system (1) approaches the specialist-free equilibrium point, as shown in Figure 2A for the typical value r = 2. For r ≤ 0.76 and other parameters as in equation (17), system (1) asymptotically approaches the superpredator-free equilibrium point, as shown in Figure 2B.

FIGURE 2

Figure 2. Trajectories of system (1) for different values of r, with other parameters as in equation (17). (A) Time series for trajectories of x, y, and z for r = 2, which approach E₂ = (0.8, 0, 0.8). (B) Time series for trajectories of x, y, and z for r = 0.3, which approach E₃ = (0.06, 0.22, 0).

For different values of the fear rate k with the rest of the data as given in equation (17), system (1) is solved numerically and illustrated in Figure 3.

FIGURE 3

Figure 3. Trajectories of system (1) using data given in equation (17) with different values of k. (A) 3D-Phase portrait of system (1) for different values of k. (B) Time series for the trajectory of x, in which x decreases but remains positive as k increases. (C) Time series for the trajectory of y, in which y increases as k increases. (D) Time series for the trajectory of z, in which z decreases as k increases and approaches zero for k > 2.

As shown in Figure 2, the superpredator decreases as k increases, with extinction for k > 2. We also varied the parameter a₁ with other parameters still fixed as in equation (17). For the range a₁ ≥ 1.3, system (1) asymptotically approaches the superpredator-free equilibrium, as shown in Figure 4A for the typical value a₁ = 1.75. Conversely, system (1) approaches asymptotically to the specialist-free equilibrium point, as shown in Figure 4B for the typical value a₁ = 0.4.

FIGURE 4

Figure 4. Trajectories of system (1) for different values of a₁ and other parameters as in equation (17). (A) Time series for trajectories of x, y, and z for a₁ = 1.75, which approach E₃ = (0.03, 0.39, 0). (B) Time series for trajectories of x, y, and z for a₁ = 0.4, which approach E₂ = (0.8, 0, 0.4).

System (1) still asymptotically approaches the coexistence equilibrium point for the data given in equation (17), when e₁ or e₃ increases or when e₂ decreases. However, for e₁ ≤ 0.44, system (1) asymptotically approaches the specialist-free equilibrium point E₂ = (0.8, 0, 0.4), as shown in Figure 4B. For e₂ ≥ 0.3 or e₃ ≤ 0.2 and the other parameters as in equation (17), the trajectories of system (1) asymptotically approach the specialist-free equilibrium point or the superpredator-free equilibrium point as shown in Figure 5.

FIGURE 5

Figure 5. Trajectories of system (1) for different values of e₂ and e₃, with other parameters as in equation (17). (A) Time series for trajectories of x, y, and z for e₂ = 0.4, which approach E₂ = (0.5, 0, 1). (B) Time series for trajectories of x, y, and z for e₃ = 0.1, which approach E₃ = (0.06, 0.58, 0).

Further investigation for the effect of varying parameters a₂ and a₃ while keeping the rest of the parameters as in equation (17) shows that the parameter a₂ has similar effects as the parameter e₂ with a bifurcation point at a₂ = 0.63. However, parameter a₃ has similar effects as parameter r, with two bifurcation points: a₃ = 1.38 and a₃ = 0.62.

Finally, for d₁ ≥ 0.31 and the rest of the parameters as in equation (17), the trajectories of system (1) asymptotically approach the specialist-free equilibrium point, as shown in Figure 6A for the typical value d₁ = 0.4. However, system (1) approaches the coexistence equilibrium point otherwise. On the other hand, for d₂ ≥ 0.12 or d₂ ≤ 0.08 and the other parameters as in equation (17), trajectories of system (1) asymptotically approach either the superpredator-free equilibrium point or the specialist-free equilibrium point, as shown in Figures 6B,C for the typical values d₂ = 0.15 and d₂ = 0.08, respectively.

FIGURE 6

Figure 6. Trajectories of system (1) for different values of d₁ and d₂, with other parameters as in equation (17). (A) Time series for trajectories of x, y, and z for d₁ = 0.4, which approach E₂ = (0.8, 0, 0.4). (B) Time series for trajectories of x, y, and z for d₂ = 0.15, which approach E₃ = (0.06, 0.58, 0). (C) Time series for trajectories of x, y, and z for d₂ = 0.08, which approach E₂ = (0.64, 0, 0.72).

Varying parameters d₁ and d₂ together so that they satisfy conditions (7a) and (7b) while keeping other parameters as in equation (17) makes the trajectories of system (1) asymptotically approach the axial equilibrium point E₁ = (1, 0, 0) as shown in Figure 7.

FIGURE 7

Figure 7. Globally asymptotically stable axial equilibrium point for d₁ = 0.8, d₂ = 0.15 and other parameters as in equation (17). (A) 3D-Phase portrait. (B) Time series for trajectories of x. (C) Time series for trajectories of y. (D) Time series for trajectories of z.

Discussion

Classical intraguild predation food web models have been extensively studied, but relatively little work has been done on the effect of fear on the prey. To the best of our knowledge, this is the first model considering prey fear of only the specialist predator. We determined the equilibria and global stability properties of the model, found bifurcations, and illustrated the theoretical behavior with numerical simulations.

Many previous studies have found that three-species intraguild predation models, in which a superpredator (IG predator) both attacks and competes with a specialist predator (IG prey), are often unstable, either because one consumer is excluded or because long feedback loops produce sustained oscillations [30]. Despite this, many natural IGP systems continue to thrive. Many empirical intraguild predation systems are entrenched in communities with alternative prey species, and standard models of intraguild predation simplify actual systems in significant ways that could affect persistence. Holt and Huxel [30] presented results of theoretical explorations of how alternative prey can influence the persistence and stability of a focal intraguild predation interaction. They reviewed the key conclusions of standard three-species IGP theory and then presented results of theoretical explorations of how alternative prey can influence the persistence and stability of a focal intraguild predation interaction. Conversely, Bai et al. [10] indicated that the intraguild predation food web model has rich and complicated dynamic behavior, and a large Allee effect in the basal prey raises the extinction risk of not only the basal prey but also the IG prey or/and IG predator.

Hossain et al. investigated fear in an intraguild predation model, in which the growth rate of a specialist predator (IG prey) is reduced due to the cost of fear of a superpredator (IG predator), and the growth rate of a basal prey is suppressed due to the cost of fear of both the IG prey and the IG predator [20]. In a three-species food chain system, they found that omnivory can cause chaos in the absence of fear. Fear, on the other hand, can help to keep the chaos under control. They also discovered that the system exhibits bistability between the IG prey-free and IG predator-free equilibrium, as well as bistability between IG prey-free and interior equilibria. Furthermore, they showed that the system can display numerous stable limit cycles for a given set of parameter values.

However, in our study, instilling the fear of a specialist predator (IG prey) in the presence of a superpredator (IG predator) had the unintended consequence of eradicating the superpredator. This demonstrates the often-surprising complexity of food-chain dynamics and has not been observed in previous articles dedicated to the study of fear [23–27]. In contrast, if the prey's intrinsic growth rate is sufficiently high, the IG prey can be eradicated, while the IG predator can be removed if the prey's intrinsic growth rate is sufficiently low. Coexistence or the eradication of both predators are other possibilities. As a result, system (1) has rich dynamical behavior, is sensitive to parameter changes, and at least one bifurcation point exists for each of the parameters. This result is comparable to decreasing the productivity of a resource (comparable to a higher fear effect in the prey in this manuscript) resulting in a reduction in the density of the highest trophic level [31].

The most important parameters that influenced the outcome were, in order, the growth rate r, the fear effect k, predation upon the prey by the specialist predator a₁, the conversion factor between the prey and the superpredator e₂, the conversion factor between the specialist predator and the superpredator e₃, the death rate of the specialist predator d₁, and the natural death rate of the superpredator d₂.

Our model contains several limitations, which should be acknowledged. We only looked at the prey's fear of the specialized predator, not of the superpredator, and we ignored the specialist predator's fear of the superpredator. Our interpretation of such a case is that the pressure of the superpredator on the basal prey is much lower than that of a specialist predator due to the existence of alternative food sources. Mass-action kinetics and a conversion factor are used to model each predator's attack rates, which is a simplification of the genuine dynamics. Finally, the predators were completely reliant on a single food source, which is not always the case.

As a result, the dynamics of a food web in which the prey is afraid of one predator but not the other can produce unexpected results. In the absence of other fear effects, future research will consider the prey's fear of the superpredator but not the specialized predator; the specialist predator's fear of the superpredator could also be incorporated.

Data availability statement

The original contributions presented in the study are included in the article/supplementary files, further inquiries can be directed to the corresponding author.

Author contributions

NF developed the model, wrote the initial draft, conducted analysis, and performed numerical simulations. RN conducted analysis and performed numerical simulations. MH designed the study. SS? rewrote and edited and manuscript. All authors contributed to the article and approved the submitted version.

Funding

SS? is supported by an NSERC Discovery Grant. For citation purposes, please note that the question mark is part of her name.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Thompson RM, Brose U, Dunne JA, Hall RO, Hladyz S, Kitching RL, et al. Food webs: reconciling the structure and function of biodiversity. Trends Ecol Evolut. (2012) 27:689–697. doi: 10.1016/j.tree.2012.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Allesina S, Tang S. Stability criteria for complex ecosystems. Nature. (2012) 483:205–08. doi: 10.1038/nature10832

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Aziz-Alaoui MA. Study of a leslie gower type tri-trophic population model. Chaos Solitons Fractals. (2002) 14:1275–1293. doi: 10.1016/S0960-0779(02)00079-6

CrossRef Full Text | Google Scholar

4. Gakkhar S, Singh B, Naji RK. Dynamical behavior of two predators competing over a single prey. BioSystem. (2007) 90:808–817. doi: 10.1016/j.biosystems.2007.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Gakkhar S, Naji RK. On a food web consisting of a specialist and generalist predator. J Biol Syst. (2003) 11:365–376. doi: 10.1142/S0218339003000956

CrossRef Full Text | Google Scholar

6. Gakkhar S, Naji RK. Order and chaos in a food web consisting of a predator and two independent preys. Commun Nonlinear Sci Numer Simulat. (2005) 10:105–120. doi: 10.1016/S1007-5704(03)00120-5

CrossRef Full Text | Google Scholar

7. Klebanoff A, Hastings A. Chaos in one predator two prey model: general results from bifurcation theory. Math. Biosciences. (1994) 122:221–223. doi: 10.1016/0025-5564(94)90059-0

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Holt RD, Polis GA. A theoretical framework for intraguild predation. Am Natur. (1997) 149:745–764. doi: 10.1086/286018

CrossRef Full Text | Google Scholar

9. Tanabe K, Namba T. Omnivory creates chaos in simple food web models. Ecology. (2005) 86:3411–3414. doi: 10.1890/05-0720

CrossRef Full Text | Google Scholar

10. Bai YK, Ruan S, Wang L. Dynamics of an intraguild predation food web model with strong Allee effect in the basal prey. Nonlinear Anal Real World Appl. (2021) 58:103206. doi: 10.1016/j.nonrwa.2020.103206

CrossRef Full Text | Google Scholar

11. Sabelis MW. Predatory arthropods. In: Natural enemies, the population biology of predators, parasites and diseases, eds M. J. Crawley (Oxford: Blackwell Scientific) (1992). p. 225–264. doi: 10.1002/9781444314076.ch10

CrossRef Full Text | Google Scholar

12. Crowder DW, Snyder WE. Eating their way to the top? Mechanisms underlying the 456 success of invasive generalist predators. Biol Invas. (2010) 12:2857–2876. doi: 10.1007/s10530-010-9733-8

CrossRef Full Text | Google Scholar

13. Schreiber SJ. Generalist and specialist predators that mediate permanence in ecological communities. J. Math. Biol. (1997) 36:133–148 doi: 10.1007/s002850050094

CrossRef Full Text | Google Scholar

14. Lima S, Dill LM. Behavioral decisions made under the risk of predation: a review and prospectus. Can. J. Zool. (1990) 68:619–640. doi: 10.1139/z90-092

CrossRef Full Text | Google Scholar

15. Panday P, Pal N, Samanta S, Chattopadhyay J. Stability and bifurcation analysis of a three-species food chain model with fear. Int J Bifurcation Chaos. (2018) 28:1850009. doi: 10.1142/S0218127418500098

CrossRef Full Text | Google Scholar

16. Clinchy LMDi, Roberts D, Zanette LY. Fear of large carnivores causes a trophic cascade. Nat Commun. (2016) 7:10698. doi: 10.1038/ncomms10698

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Panday NP, Samanta S, Chattopadhyay J. A three species food chain model with fear induced trophic cascade. Int J Appl Comput Math. (2019) 5:1–26. doi: 10.1007/s40819-019-0688-x

CrossRef Full Text | Google Scholar

18. Panday P, Samanta S, Pal N, Chattopadhyay J. Delay induced multiple stability switch and chaos in a predator–prey model with fear effect. Math Comput Simulat. (2020) 172:134–58. doi: 10.1016/j.matcom.2019.12.015

CrossRef Full Text | Google Scholar

19. Cong MF, Zou X. Dynamics of a three-species food chain model with fear effect. Commun Nonlinear Sci Numer Simulat. (2021) 99:105809. doi: 10.1016/j.cnsns.2021.105809

CrossRef Full Text | Google Scholar

20. Hossain NP, Samanta S, Chattopadhyay J. Fear induced stabilization in an intraguild predation model international. J Bifurcation Chaos. (2020) 30:2050053. doi: 10.1142/S0218127420500534

CrossRef Full Text | Google Scholar

21. Mukherjee D. Role of fear in predator–prey system with intraspecific competition. Math Comput Simulat. (2020) 177: 263–275. doi: 10.1016/j.matcom.2020.04.025

CrossRef Full Text | Google Scholar

22. Ibrahim HA, Naji RK. Chaos in beddington–deangelis food chain model with fear. J Phys Conf Ser. (2020) 1591:012082. doi: 10.1088/1742-6596/1591/1/012082

CrossRef Full Text | Google Scholar

23. Mondal B, Roy S, Ghosh U, Tiwari PK. A systematic study of autonomous and nonautonomous predator–prey models for the combined effects of fear, refuge, cooperation and harvesting. Eur Phys J Plus. (2022) 137:724. doi: 10.1140/epjp/s13360-022-02915-0

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Roy S, Tiwari PK, Nayak H, Martcheva M. Effects of fear, refuge and hunting cooperation in a seasonally forced eco-epidemic model with selective predation. Eur Phys J Plus. (2022) 137:1–31. doi: 10.1140/epjp/s13360-022-02751-2

CrossRef Full Text | Google Scholar

25. Maity PKTi, Pal S. An ecoepidemic seasonally forced model for the combined effects of fear, additional foods and selective predation. J Biol Syst. (2022) 12:1–37. doi: 10.1142/S0218339022500103

CrossRef Full Text | Google Scholar

26. Pal PKTi, Pal N. Bifurcations, chaos, and multistability in a nonautonomous predator–prey model with fear Chaos: an interdisciplinary. J Nonlinear Sci. (2021) 31:123–34. doi: 10.1063/5.0067046

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Tiwari PK, Verma M, Pal S, Kang Y, Misra AK. A delay nonautonomous predator–prey model for the effects of fear, refuge and hunting cooperation. J Biol Syst. (2021) 29:927–69. doi: 10.1142/S0218339021500236

CrossRef Full Text | Google Scholar

28. Wang X, Zanette L, Zou X. Modelling the fear effect in predator–prey interactions. J Math Biol. (2016) 73:1179–1204. doi: 10.1007/s00285-016-0989-1

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Perko L. Differential Equations and Dynamical Systems, Volume 7. New York: Springer (1996). doi: 10.1007/978-1-4684-0249-0

CrossRef Full Text | Google Scholar

30. Holt RD, Huxel GR. Alternative prey and the dynamics of intraguild predation: theoretical perspectives. Ecology. (2007) 88:2706–2712. doi: 10.1890/06-1525.1

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Gurney W, Nisbet RM. Ecological Dynamics. Oxford: Oxford University Press. (1998).

Google Scholar