The effects of fear and delay on a predator-prey model with Crowley-Martin functional response and stage structure for predator

1.   Introduction

In the study of ecological dynamics, the diversity of populations is a very important index for the continuous development of populations. The dynamics of such systems as the predator-prey system, competition system, co-operation system and more are a vital research area ^[1,2,3].

In the real world, the prey-predator system is the most popular, in which the predation is an important feature in evaluating the interactivity between prey and predator since the predator's survival depends on the successful consumption on prey populations, as well as regulating the overall population size. The functional response is such an index that reflects the consumption rate of predator populations on prey populations. It is characterized as the predation efficiency of the prey populations consumed by predator populations and it is usually decided by the biomass of prey populations, predator populations or both of them. Many authors have applied the Holling-type functional response to depict the predation relationship^[4,5,6]. The Holling-type functional response is only dependent on the prey's biomass, but when it is at a low level, the number of prey predated by predator usually does not follow a linearly increasing curve with the prey's biomass. Some researches have proved that the consumption of prey by predator depends not only on the prey's biomass, but also on the predator's biomass^{[7,8,9,10,11,12]}. As a result, many predator-dependent functional responses have been introduced, such as Crowley-Martin^[9,10], Beddington-DeAngelis^[11] and Hassell-Varley^[12], which are more appropriate choices in mathematical modeling by experimental evidence. Taking into account the interference among predator populations, the predation rate is often declined when the predator populations are higher, even if the prey populations are higher, too. In this situation, the interactions between predator populations should be considered and the Crowley-Martin functional response is better to describe this kind of interplay than the other two, as it is more realistic in reflecting this biological background^[13].

In the ecological world, almost all kinds of prey and predator species have two stages of structure: The immature stage and the mature stage. In different stages, the individuals exhibit very different biological behaviors, such as the ability of predating and reproduction and susceptibility to predator. Hence, in mathematical modeling, the assumption that every predator has the same skill to predate prey is very impractical, whereas a great quantity of research works were performed on the assumption that every individual undergoes only one stage in its whole living history^{[14,15,16,17]}. Compared with those single-stage systems, the authors in ^[14,15] showed that it was more realistic for people to predict the system dynamics if the stage structure was incorporated in the mathematical models. Taking the stage structure of species into account, recently some nice results of prey-predator systems with stage structure have been reported^[16,17,18]. For example, the authors in ^[18] presented a predator-prey model where the predator population was divided into two stages and the immature predator was not capable of predation. Many researchers have shown that the parameters representing the stage structure of individuals can change the behaviors of species from stable status to fluctuating status, and vice versa^[19]. Therefore, it is necessary to consider the effect of the stage structure of species on the system dynamics.

Considering the Crowley-Martin type functional response of predator to prey and the stage structure of predators, we have the following system:

where $P_1, P_2$ and $P_3$ (the simplicity of $P_1(t)$ $P_2(t)$ and $P_3(t)$) are the biomass of prey, immature predator and mature predator at time $t$, respectively. $r > 0$ is the growth rate and $d_1$ is the natural death rate of prey. $c_1$ is the intra-specific competition rate between individuals. The immature predator has no predation ability with natural death rate $d_2$, and the transversion rate of predator from immature to mature is $c_2$. The predator functional response to prey is supposed to be the Crowley-Martin type. Parameters $b_1, b_2$ and $b_3$ describe the effect of handling time, capture rate and the magnitude of interference among predators, respectively. $\eta$ is the conversion rate of predators. The intra-specific competition rate between adult predator populations is $c_3$, and the natural death rate is $d_3$.

The functional response of predator to prey is relatively easily observed^[20]. However, besides the direct predation, more and more experimental evidences show that predators can create fear in prey populations by their voice or smell. This non-consumptive action of predators often induces some significant impacts on prey and forces them to change their behaviors ^[21]. In fact, in some cases, the fear effect of predators can work more efficiently than direct killing, resulting in the reduction of the number of prey or even the extinction of prey, even if the direct predation is absent. The fear effect from predators is a manifestation of sustained psychological stress because the prey species are always worried about the possible attack of a predator. To protect themselves, the prey will take some measures to prevent the predator's attack and present some anti-predator behaviors. Due to the effect of fear, some important aspects of prey behaviors are varied like hunting, foraging, reproduction ability and so on^{[22,23,24,25]}. Empirical observations show that the prey shifts to a relatively safer place with less food and is accompanied by starvation to avoid the predator's attack^[26], forages in a less amount due to the predator's fear ^[27], spends less time hunting to prevent the predation risk of mountain lions^[28] or even changes their reproductive physiology due to the predation risk of wolves^[29].

Therefore, it is interesting to take into account the fear effect as well as the counter-predation sensitiveness. We then get the modified version of (1.1) below.

where $f$ denotes the fear effect induced by the predator and $k$ is the sensitiveness of prey to predation risk. That is, if the value of parameter $k$ is higher, i.e., the counter-predation level of prey is higher, then the predation rate is to be lower instead.

On the other hand, except for the effect of direct and indirect predation of a predator, time delay is another important factor that impacts the dynamics of prey-predator interactions^[30]. Time delay emerges in almost all ecological processes. For example, there is a time lag for the immature individual to become mature, which is called maturation delay^[31]. For predators, they absorb energy by consuming prey to keep their survival, but the conversion of energy is not finished instantly, and there exists a digestion delay^[30]. The gestation delay of the newborn is more popular ^[30,32]. In addition, due to the fear effect induced by the predator, the reproduction efficiency of prey is reduced by many manners, which is all invisible immediately. It takes a relatively long period of time to manifest the impact of a predator's fear on the birth rate of prey. The fear delay refers to the time when prey population starts to change their physiological characteristics due to fear of predators. That is, there is always a time lag for the fear effect on the system dynamics^{[33,34,35,36]}. The delayed biological system has been extensively studied by use of the theory of delay differential equations. Many results indicate that time delays can destroy the system stability. For instance, the authors in ^[37] considered a predator-prey system with mature delay and established the conditions of global stability and bifurcation. Mortoja et al.^[38] considered the dynamics of a predator-prey model with mature delay and disease in prey population. Currently, the studies of how multi-delays affect the system dynamics have been executed and lots of helpful results have been obtained. Particularly the bifurcation thresholds of delays have been derived ^[35,39,40]. By study of a delayed predator-prey system with no intra-specific competitions, the authors in ^[39] observed that the system was stable when the time delay was low, while as for the increasing of time delays, the stability was lost and bifurcation appeared, and when the delay was large enough, then chaotic behavior occurred. Beretta et al.^[40] involved a kind of geometric method into a delayed differential equation to determine whether the stability was changed or not. All these reveal that time delays play a key role in the analysis of system stability.

Incorporating the response delay of prey to the fear of predator and mature delay of predator, (1.2) becomes

The initial conditions are

For biological justifications, we assume all parameters in the models are positive and $\varphi_i(\theta)\geq0.$ In biological research, the system stability is one of the most important issues. We aim to study the local asymptotical stability, parameter bifurcations and the cost and benefit of fear and time delays on the stability of model (1.3) both theoretically and numerically.

The innovations and contributions of this paper are as follows:

● The different predation ability of juvenile and adult predators, the mutual interference between predators, the fear induced by predator on prey as well as the counter-predator sensitivity of prey, and the fear and mature delays are incorporated simultaneously in the model.

● The effects of the conversion rate of prey to predator, stage structure, fear and time delays on the system stability are investigated. Particularly, the critical values of some bifurcation parameters are obtained by numerical simulations.

● The complicated ecological phenomena like bubble and bistability are presented intuitively.

The rest of the article is constructed as follows. The existence, boundedness and extinction of species are presented in section two. The local stability and bifurcation of equilibrium points of (1.2) and (1.3) are executed in sections three and four, respectively. Simulations and numerical analysis are given in section five. Finally, a brief conclusion and discussion are given to end this paper in section six.

2.   Preparations

By the theory of the delayed functional differential equation^[30], if the initial values of system (1.3) are positive, then there is a unique positive solution. That is, system (1.3) is positive invariant. Next, we discuss the boundedness of solutions. From the first equation of (1.3), we have

then

By the boundedness of $P_1$, together with the second equation of (1.3), we have

Then $P_2\leq \frac{\eta b_2(r-d_1)}{kb_3c_1 (d_2+c_2)}$ and the boundedness of $P_2$ is derived. By the third equation of (1.3) as well as the boundedness of $P_2$, we obtain

Then, $P_3\leq \frac{c_2\eta b_2(r-d_1)}{kb_3c_1 d_3(d_2+c_2)}$ and the boundedness of $P_3$ is derived. Therefore, the boundedness of all species is derived.

We proceed with the discussion of the extinction of species. It is clear that if $r < d_1$, then $P_1$ is extinct, which results in the extinction of the juvenile and adult predator. Now, we study the case of $r > d_1$. Define $P(t) = P_2(t)+P_3(t)$, differentiating $P(t)$ w.r.t $t$ along the solution of (1.3), and we have

that is,

Consequently, we derive that

Hence,

Therefore, we can obtain that the adult predator $P_3$ is extinct if $\frac{\eta b_2}{b_1} < d_3$, which means the juvenile predator is also extinct. Here, the prey is stable with $\lim\limits_{t\rightarrow \infty} P_1(t) = \frac{r-d_1}{c_1}.$

Finally, in this part we give the equilibrium points of (1.2). The equilibrium points can be classified into the following three classes:

● The trivial equilibrium point $E_0(0, 0, 0),$ which always exists.

● The predator-free equilibrium point $\overline{E}(\overline {P}_1, 0, 0)$ with $\overline {P}_1 = \frac{r-d_1}{c_1}$, which exists under $r > d_1$.

● The positive interior equilibrium $\widetilde{E}(\widetilde{P}_1, \widetilde{P}_2, \widetilde{P}_3),$ which conditionally exists such that the right sides of (1.2) are equivalent to zero, which cannot be expressed by explicit function and is to be discussed numerically in section five.

3.   Dynamics of (1.2)

For the discussion of local asymptotic stability (LAS) of system (1.2), we start with the linearized system of (1.2). For simplicity, let $g(P_1, P_3) = 1+b_1P_1+kb_3P_3+kb_1b_3P_1P_3$, and rewrite (1.2) as

then we can linearize system (1.2) as $P'(t) = JP(t)$ with Jacobian matrix

where

For convenience of using the Sotomayor's bifurcation theorem^[41] latter, we compute the partial derivatives of $F$ on viable $P_i(i = 1, 2, 3)$ as follows:

With these preparations, we analyze the dynamics of (1.2) near the equilibrium points.

(1) Dynamics near $E_0$

At $E_0(0, 0, 0)$, the Jacobian matrix $J_0$ of system (1.2) is,

It is not difficult to compute that the roots of characteristic equation are $\lambda_1 = r-d_1$, $\lambda_2 = -d_2-c_2 < 0$, $\lambda_3 = -d_3 < 0.$ Thus, if $r < d_1$, then it is LAS, otherwise it is unstable^[42]. We discuss the existence of parameter bifurcation $d_1$ at $r-d_1 = 0$. Denote $d_1^T$ the value of $d_1-r = 0$. The eigenvectors of $\lambda_1 = 0$ corresponding to $J_0$ and $J_0^T$ are $v_1 = (1, 0, 0)^T, v_2 = (1, 0, 0)^T,$ respectively. The partial derivatives of $F(P_1, P_2, P_3)$ on viable $d_1$ and $F_{d_1}$ on viable $P_i(i = 1, 2, 3)$ are

Now, we verify the conditions of bifurcation theorem point by point.

Consequently, by use of Sotomayor's bifurcation theorem, there exists a transcritical bifurcation near $E_0$ at the threshold $d_1^T = r.$

(2) Dynamics near $\overline{E}$

We linearize (1.2) and get the Jacobian matrix $J_1$ at $\overline{E}$ as below,

By computation, one characteristic root is $\lambda_1 = -(r-d_1),$ and the other two meet

namely,

If $r-d_1 > 0$ and $(d_2+c_2)d_3- \frac{\eta c_2b_2(r-d_1)}{c_1+b_1(r-d_1)} > 0$, then all characteristic roots are negative, and system (1.2) is LAS. Otherwise, it is unstable.

If $\lambda_1 = r-d_1 = 0,$ by the same manner as before, we can derive that (1.2) undergoes a transcritical bifurcation near $\overline{E}$ at $d_1^T = r.$

If $\lambda_1 = -(r-d_1) < 0,$ $(d_2+c_2)d_3- \frac{\eta c_2b_2(r-d_1)}{c_1+b_1(r-d_1)} = 0,$ then $\lambda_2 = -d_2-d_3-c_2 < 0, \lambda_3 = 0.$ Choose $\eta$ as a bifurcation parameter and analyze whether the bifurcation happens or not. Let $\eta = \eta^T$ such that $(d_2+c_2)d_3- \frac{\eta c_2b_2(r-d_1)}{c_1+b_1(r-d_1)} = 0.$ The eigenvector of zero characteristic value corresponding to $J_1$ is $v_1 = (\overline{\varrho}_1, \overline{\varrho}_2, 1)^T,$ and the eigenvector corresponding to $J_1^T$ is $v_2 = (0, \overline{\varrho}_3, 1)^T,$ where $\overline{\varrho}_1 = \frac{1}{r-d_1}, \overline{\varrho}_2 = \frac{d_3}{c_2}, \overline{\varrho}_3 = \frac{d_3(c_1+b_1(r-d_1))}{\eta b_2(r-d_1)}.$ Compute the partial derivatives about parameter $\eta$, then

The conditions of the bifurcation theorem are validated as follows:

$\begin{array}{lll} (i)\quad v_2^T[F_\eta(\overline{E}, \eta = \eta^T] = 0. \end{array}$

$\begin{array}{lll} (ii)\quad v_1^T[DF_\eta(\overline{E}, \eta = \eta^T)v_2] & = & \left(\overline{\varrho}_1, \overline{\varrho}_2, 1\right) \left(\begin{array}{ccc} 0 & 0 & 0 \\ 0 & 0 & \frac{ b_2(r-d_1)}{c_1+b_1(r-d_1)} \\ 0 & 0 & 0 \\ \end{array} \right) \left(\begin{array}{c} 0 \\ \overline{\varrho}_3 \\ 1 \\ \end{array} \right) \\ & = & \frac{ b_2d_3(r-d_1)}{c_2(c_1+b_1(r-d_1))}\neq0. \end{array}$

$\begin{array}{lll} (iii) & D^2F(\overline{E}, \eta = \eta^T)(v_1, v_1) = F_{P_1P_1}(E_1, \eta = \eta^T)\overline{\varrho}_1^2 +2F_{P_1P_2}(E_1, \eta = \eta^T)\overline{\varrho}_1\overline{\varrho}_2 \\ & \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; +2F_{P_1P_3}(E_1, \eta = \eta^T)\overline{\varrho}_1 +F_{P_2P_2}(E_1, \eta = \eta^T)\overline{\varrho}_2^2 \\ & \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; + 2F_{P_2P_3}(E_1, \eta = \eta^T)\overline{\varrho}_2+F_{P_3P_3}(E_1, \eta = \eta^T) \\\\ & = \left(\begin{array}{c} -2c_1\\ 0 \\ 0 \end{array} \right)\overline{\varrho}_1^2 + 2\left(\begin{array}{c} -rkf- \frac{b_2}{(1+b_1\overline{P}_1)^2} \\ \frac{\eta b_2}{(1+b_1\overline{P}_1)^2}\\ 0 \end{array} \right) \overline{\varrho}_1 +\left(\begin{array}{c} 2rk^2f^2P_1^1+ \frac{2kb_2b_3\overline{P}_1}{1+b_1\overline{P}_1} \\ - \frac{2\eta kb_2b_3\overline{P}_1}{1+b_1\overline{P}_1}\\ -2c_3 \end{array} \right). \end{array}$

Consequently,

Therefore, Sotomayor's bifurcation theorem implies that system (1.2) has a transcritical bifurcation near $E_1$ at $\eta = \eta^T$. Similarly we can obtain the bifurcation thresholds of other parameters (for example, $c_2$).

(3) Dynamics near $\widetilde{E}$

Next, we research the stability and bifurcation around the coexistence equilibrium state $\widetilde{E}$. The variational matrix $J_2$ at $\widetilde{E}$ is

where $\quad \widetilde{a}_{11} = \frac{r}{1+kf\widetilde{P}_3 }-d_1-2c_1\widetilde{P}_1-\frac{b_2 \widetilde{P}_3(1+kb_3\widetilde{P}_3)}{g^2(\widetilde{P}_1, \widetilde{P}_3)}, \quad\widetilde{a}_{13} = \frac{-rkf\widetilde{P}_1}{(1+kf\widetilde{P}_3)^2} -\frac{b_2\widetilde{P}_1(1+b_1\widetilde{P}_1)}{g^2(\widetilde{P}_1, \widetilde{P}_3)}, \\ \widetilde{a}_{21} = \frac{\eta b_2\widetilde{P}_3(1+kb_3\widetilde{P}_3)}{g^2(\widetilde{P}_1, \widetilde{P}_3)}, \quad a_{22} = -d_2-c_2, \quad\widetilde{a}_{23} = \frac{\eta b_2\widetilde{P}_1(1+b_1\widetilde{P}_1)}{g^2(\widetilde{P}_1, \widetilde{P}_3)}, \quad a_{32} = c_2, \quad a_{33} = -d_3-2c_3\widetilde{P}_3.$

The characteristic equation is

where $\delta_1 = -(\widetilde{a}_{11}+\widetilde{a}_{22}+\widetilde{a}_{33}), \delta_2 = \widetilde{a}_{11}\widetilde{a}_{22}+\widetilde{a}_{11}\widetilde{a}_{33} +\widetilde{a}_{22}\widetilde{a}_{33}-\widetilde{a}_{23}\widetilde{a}_{32}, \delta_3 = \widetilde{a}_{11}\widetilde{a}_{23}\widetilde{a}_{32} -\widetilde{a}_{11}\widetilde{a}_{22}\widetilde{a}_{33}-\widetilde{a}_{13}\widetilde{a}_{21}\widetilde{a}_{32}.$ By the Rowth-Hurwitz theorem, if $\delta_1 > 0, \delta_3 > 0, \delta_1\delta_2 > \delta_3$, then it is LAS, otherwise it is unstable.

If $\delta_3 = 0$, then there is a zero eigenvalue of (3.1). Choose $c_2$ as a bifurcation parameter and denote $c_2^S$ such that $\delta_3 = 0$, then we study the bifurcation about parameter $c_2$. The eigenvectors of the zero eigenvalue corresponding to matrix $J_2$ and $J_2^T$ are $\widetilde{v}_1 = (\widetilde{\varrho}_1, \widetilde{\varrho}_2, 1)^T$ and $\widetilde{v}_2 = (\widetilde{\varrho}_3, 1, \widetilde{\varrho}_4)^T$, respectively, where $\widetilde{\varrho}_1 = -\frac{\widetilde{a}_{13}}{\widetilde{a}_{11}}, \widetilde{\varrho}_2 = -\frac{\widetilde{a}_{33}}{\widetilde{a}_{32}}, \widetilde{\varrho}_3 = -\frac{\widetilde{a}_{21}}{\widetilde{a}_{11}}, \widetilde{\varrho}_4 = -\frac{\widetilde{a}_{22}}{\widetilde{a}_{32}}.$ An easy computation results in

Then, with the help of Matlab2014a, we compute that

$(i)\; \quad \widetilde{v}_2^T[F_{c_2}(\widetilde{E}, c_2 = c_2^S]\neq0$.

$(ii)\quad\widetilde{v}_2^T[D^2F(\widetilde{E}, c_2 = c_2^S)(\widetilde{v}_1, \widetilde{v}_1)]\neq0.$

Therefore, by the bifurcation theorem (Theorem 1 in Page 338 of ^[41]), (1.2) has a saddle-node bifurcation near $\widetilde{E}$ at $c_2 = c_2^S$.

Next, we investigate the Hopf bifurcation of (1.2) near $\widetilde{E}$. Suppose (1.2) is LAS; that is, $\delta_1 > 0$, $\delta_3 > 0, \delta_1\delta_2 > \delta_3$. If $\delta_3 = \delta_1\delta_2,$ then the characteristic equation becomes

It is clear that there is a negative root $\lambda_1 = -\delta_1.$ Another two characteristic roots of (3.2) are $\lambda_{2, 3} = \pm\sqrt{\delta_2}i.$ Choose $k$ as the bifurcation parameter, then the roots of (3.2) can be written as

Denote $k = k^H$ satisfying $\delta_3-\delta_1\delta_2 = 0$. Now we verify the tranversality condition $\frac{d}{dk}Re(\lambda_i(k))|_{k = k^H}\neq0$, $i = 2, 3.$ Putting $\lambda_2 = \mu(k)+ i\nu(k)$ into (3.1) and calculating the derivative w.r.t $k$, noting $\mu(k) = 0$ and $\nu(k) = \sqrt{\delta_2(k)},$ we have

where $\wedge_1(k) = 2\delta_2(k), \wedge_2(k) = 2\delta_1(k)\sqrt{\delta_2(k)}, \wedge_3(k) = \delta_1'(k)\delta_2(k) -\delta_3(k), \wedge_4(k) = -\delta_2(k)\sqrt{\delta_2(k)}.$ If

then we compute that

Note that $\lambda_1(k) = -\delta_1\neq0,$ so by Theorem 2 in Page 353 of ^[41], a Hopf bifurcation occurs near $\widetilde{E}$ at $k = k^H.$ Similarly, we can obtain the bifurcation thresholds of other parameters. We summarize our findings as follows.

Theorem 3.1. For system (1.2), there are the following properties.

(i) For equilibrium point $E_0(0, 0, 0)$, if $r < d_1$, then system (1.2) is LAS, and there exists a transcritical bifurcation near $E_0$ at the threshold $d_1^T = r.$

(ii) For the predator-free equilibrium state $\overline{E}(\overline{P}_1, 0, 0)$, if $r-d_1 > 0$ and $(d_2+c_2)d_3- \frac{\eta c_2b_2(r-d_1)}{c_1+b_1(r-d_1)} > 0$, then system (1.2) is LAS. If $r-d_1 = 0$ or $(d_2+c_2)d_3- \frac{\eta c_2b_2(r-d_1)}{c_1+b_1(r-d_1)} = 0$, then (1.2) undergoes transcritical bifurcations near $\overline{E}$ with threshold $d_1^T = r$ or $\eta = \eta^T$, respectively.

(iii) For the interior equilibrium point $\widetilde{E}(\widetilde{P}_1, \widetilde{P}_2, \widetilde{P}_3)$, if $\delta_1 > 0, \delta_3 > 0, \delta_1\delta_2 > \delta_3$, then system (1.2) is LAS, otherwise it is unstable. If $\delta_3 = 0$, then a saddle-node bifurcation occurs near $\widetilde{E}$ at $c_2 = c_2^S$. If $\delta_1\delta_2 = \delta_3, \delta_2^2(k)+\delta_1(k)(\delta'_1(k)\delta_2(k)-\delta_3(k))\neq0,$ then a Hopf bifurcation occurs near $\widetilde{E}$ at $k = k^H$.

Remark 3.1. The instability of $E_0$ ensures the existence of $\overline{E}$ and vice versa, but the relationship between the instability of $\overline{E}$ and the existence of $\widetilde{E}$ is not clear. In addition, by use of similar deduction, the bifurcation thresholds of other parameters, such as, the functional response parameters like $b_i$ $(i = 1, 2, 3)$, can be obtained similarly.

4.   Dynamics of (1.3)

As stated before, the time delay usually plays a crucial role in the system dynamics, so we investigate the delayed case now. By applying Taylor's formula, we linearize the delayed system (1.3) and get

with Jacobian matrix

where $\breve{a}_{13} = -\frac{b_2P_1(1+b_1P_1)}{g^2(P_1, P_3)}, \hat{a}_{13} = \frac{-rkfP_1}{(1+kfP_3)^2},$ and other parameters are the same as before. The variational matrix of delayed system (1.3) is

Thus the characteristic equation of (1.3) is

i.e.,

where $\Upsilon_1 = -(a_{11}+a_{22}+a_{33}), \Upsilon_2 = a_{11}a_{22}+a_{11}a_{33}+a_{22}a_{33}-a_{23}a_{32}e^{-\lambda\tau_2}$, $\Upsilon_3 = -a_{11}a_{23}a_{32}$, $\Upsilon_4 = a_{11}a_{23}a_{32}-\breve{a}_{13}a_{21}a_{32}$, $\Upsilon_5 = -\hat{a}_{13}a_{21}a_{32}.$ Clearly, if $\tau_1 = \tau_2 = 0$, it is accordant with (3.1).

For $E_0$ and $\overline{E}$, the dynamics are the same as before, so we only discuss the dynamics near the coexistence equilibrium state $\widetilde{E}(\widetilde{P}_1, \widetilde{P}_2, \widetilde{P}_3)$.

We suppose (1.3) is LAS in the absence of delays, i.e., conditions $\delta_1 > 0, \delta_2 > 0, \delta_1\delta_2 > \delta_3^2$ hold. According to the scenarios of one delay or two delay appearance, we classify our discussion into the following two cases.

$(1)$ The scenario of one delay.

● $\tau_1 > 0, \tau_2 = 0$.

For system (1.3), if there is only one delay $\tau_1$ and the mature delay is free ($\tau_2 = 0$), then we have the following conclusion.

Theorem 4.1. Suppose the delayed system (1.3) is LAS in the absence of delays. If $\tau_1 > 0, \tau_2 = 0$, then there exists a Hopf bifurcation for system (1.3) near $\widetilde{E}$ at $\tau_1 = \widetilde{\tau}_1$ under the condition that

where $\widetilde{\tau}_1, \Xi_i (i = 1, 2, 3, 4)$ are defined later in the proof.

Proof. At $\tau_2 = 0$, the characteristic Eq (4.1) becomes

Take $\tau_1$ as a bifurcation parameter and assume the characteristic root of (4.2) is in the general form $\lambda(\tau_1) = \alpha(\tau_1)+i\beta(\tau_1)$. Due to the LAS of (1.2) with no delays, that is, $\alpha(\tau_1) < 0$, by the continuity of functions there exists a $\widetilde{\tau}_1$ such that $\alpha(\widetilde{\tau}_1) = 0$. When $\tau_1 < \widetilde{\tau}_1$, it remains the LAS, while if $\tau_1 > \widetilde{\tau}_1$, it turns unstable from stable and may undergo bifurcations. Based on the analysis, we substitute the characteristic root $\lambda$ into (4.2) and then separate the real part and imaginary part respectively, then

If $\widetilde{E}$ alters its stability from stable to unstable, then the necessary condition is that the roots of (4.2) are in the following forms: $\lambda(\tau_1) = i\beta(\tau_1)$, i.e., $\alpha(\tau_1) = 0,$ then

Squaring both sides of (4.3) and adding them leads to

Denote $\gamma_1 = \Upsilon_1^2-2\Upsilon_2, \gamma_2 = \Upsilon_2-2\Upsilon_1(\Upsilon_3+\Upsilon_4), \gamma_3 = -\Upsilon_5^2$. For simplicity we rewrite (4.4) in the below form,

Let $h(x) = x^6+\gamma_1x^4+\gamma_2x^2+\gamma_3$, then $h(0) = \gamma_3 = -\Upsilon_5^2 < 0, h(\infty) = \infty.$ By the continuity of $h(x)$, we obtain that (4.5) owns one positive solution at least and three positive solutions at most.

Without loss of generality, we suppose it has three solutions denoted by $\beta_i, i = 1, 2, 3,$ then for some invariant $\beta_i$, by (4.3), we can seek out a sequence of $\tau_1^{(i, k)}$ such that

Let $\widetilde{\tau}_1 = \min\limits_{i, k}\tau_1^{(i, k)}$, then by Butler's lemma, $\widetilde{E}$ remains stable when $\tau_1 < \widetilde{\tau}_1$. Next we verify

i.e., the tranversality condition of Hopf bifurcation holds. Differentiating (4.2) w.r.t. $\tau_1$ and substituting $\tau_1 = \widetilde{\tau}_1, \alpha(\widetilde{\tau}_1) = 0, \beta_i(\widetilde{\tau}_1) = \widetilde{\beta}_i$, we have

where

Under the condition that $\Xi_1\Xi_4-\Xi_2\Xi_3 > 0,$ solving (4.6) results in

That is, the tranversality condition holds and, hence, a Hopf bifurcation occurs near $\widetilde{E}$ at $\tau_1 = \widetilde{\tau}_1$. □

● $\tau_1 = 0, \tau_2 > 0$.

In this case, we take $\tau_2$ as a bifurcation parameter, then the characteristic equation is

Repeating the previous process, we have the following result.

Theorem 4.2. Suppose the delayed system (1.3) is LAS in the absence of delays. If $\tau_1 = 0, \tau_2 > 0$, then there exists a Hopf bifurcation for system (1.3) near $\widetilde{E}$ at $\tau_2 = \widetilde{\tau}_2$ under the condition that $\xi_1\xi_4-\xi_2\xi_3 >$$0$, where

$\widetilde{\beta}_i = \beta_i(\widetilde{\tau}_2), \beta_i (i = 1, 2, 3)$ is the root of the following equation:

$(2)$ The scenario of two delays.

If there are two delays for (1.3), i.e., $\tau_i > 0(i = 1, 2),$ then we fix one delay in the stable interval and choose the other as the Hopf bifurcation parameter. There are the following two scenarios.

● $\tau_2\in(0, \widetilde{\tau}_2)$ and $\tau_1 > 0$.

We fix $\tau_2$ at some arbitrary point $\tau_2 = \breve{\tau}_2\in(0, \widetilde{\tau}_2)$, and take $\tau_1$ as a bifurcation parameter, then we obtain the criterion of the existence of Hopf bifurcation of system (1.3).

Theorem 4.3. Suppose the delayed system (1.3) is LAS in the absence of delays. If $\tau_1 > 0, \tau_2\in(0, \widetilde{\tau}_2)$, then there exists a Hopf bifurcation for system (1.3) near $\widetilde{E}$ at $\tau_1 = \widehat{\tau}_1$ under the conditions that

where $\Upsilon_i(i = 3, 4, 5)$ remains unchanged as those in (4.2), and $\Omega_i(i = 1, 2, 3, 4)$ are defined in the proof.

Proof. We fix $\tau_2 = \breve{\tau}_2\in(0, \widetilde{\tau}_2)$, then by (4.1), the characteristic equation is

Based on the same method, let $\lambda(\tau_1) = \alpha(\tau_1)+i\beta(\tau_1)$ be a characteristic root of (4.8), then putting it into (4.8) and separating the real and imaginary parts leads to

The stability near $\widetilde{E}$ changes from stable to unstable under the condition that there is a pair of purely imaginary roots for (4.8), i.e., $\alpha(\tau_1) = 0,$ and whence (4.9) becomes

Squaring the left and right sides of (4.10) and then adding them, we have

where

Let $h(x) = x^6+\Lambda_1x^4+\Lambda_2x^2+\Lambda_3\sin\breve{\tau}_2x+\Lambda_4\cos \breve{\tau}_2x+\Lambda_5$. It is clear that $h(\infty) = \infty,$ then under the condition that $h(0) = (\Upsilon_3+\Upsilon_4)^2-\Upsilon_5^2 < 0$, by the continuity of $h(x)$, (4.13) has at least one positive root and at most a finite number of positive solutions denoted by $\beta_j, j = 1, \cdots, N.$ Thus, for some fixed solution $\beta_j(j\in{1, 2, \cdots, N})$, by (4.10) we can find a sequence of ${\tau_1^{(j, k)}}$ such that

Let $\widehat{\tau}_1 = \min\limits_{j, k}\tau_1^{(j, k)}$, then by Butler's lemma, $\widetilde{E}$ remains stable when $\tau_1 < \widehat{\tau}_1$. Again, we verify the following tranversality condition

Differentiating (4.9) w.r.t. $\tau_1$ and putting $\tau_1 = \widehat{\tau}_1, \alpha(\widehat{\tau}_1) = 0, \widehat{\beta}_j = {\beta}_j(\widehat{\tau}_1)$ for some $j\in\{1, 2, \cdots, N\}$, we have

where

$\begin{array}{lll} \Omega_1 = -3\widehat{\beta}_j^2+\Upsilon_2-\Upsilon_4\breve{\tau}_2\cos\widehat{\beta}_j\breve{\tau}_2 -\Upsilon_5(\breve{\tau}_2+\widehat{\tau}_1)\cos\widehat{\beta}_j(\widehat{\tau}_1+\breve{\tau}_2), \\ \Omega_2 = 2\Upsilon_1\widehat{\beta}_j+\Upsilon_4\breve{\tau}_2\sin\widehat{\beta}_j\breve{\tau}_2+\gamma_5(\tau_1+\breve{\tau}_2) \sin\widehat{\beta}_j(\widehat{\tau}_1+\breve{\tau}_2), \\ \Omega_3 = \gamma_5\widehat{\beta}_j\sin\widehat{\beta}_j(\widehat{\tau}_1+\breve{\tau}_2), \qquad \Omega_4 = \gamma_5\widehat{\beta}_j\cos\widehat{\beta}_j(\widehat{\tau}_1+\breve{\tau}_2). \end{array}$

Solving them results in

provided that $\Omega_1\Omega_3-\Omega_2\Omega_4 > 0.$ That is, the tranversality condition holds and, hence, a Hopf bifurcation occurs near $\widetilde{E}$ at $\tau_1 = \widehat{\tau}_1$. □

● $\tau_1\in(0, \widetilde{\tau}_1)$ and $\tau_2 > 0$.

Fix $\tau_1$ at $\tau_1 = \breve{\tau}_1\in(0, \widetilde{\tau}_1)$ and take $\tau_2$ as a free parameter, then we have the similar theorem.

Theorem 4.4. Suppose the delayed system (1.3) is LAS in the absence of delays. If $\breve{\tau}_1\in(0, \widetilde{\tau}_1), \tau_2 > 0$, there exists a Hopf bifurcation for system (1.3) near $\widetilde{E}$ at $\tau_2 = \widehat{\tau}_2$ under the conditions that

where

and $\widehat{\beta}_j = \beta_j(\widehat{\tau}_2)$, $\beta_j(j = 1, 2, \cdots, N)$ is the root of the following equation:

5.   Numerical analysis

In this part, we give some examples to simulate and validate the previous theoretical analysis and take the discussion further to discuss how the fear and time delays affect the stability of the equilibrium points. For this purpose, we fix a set of parameters and change the values of fear level $f$ induced by the predator, the counter-predation parameter $k$ of prey, the transversion level $\eta$ of predator from prey and the maturing efficiency $c_2$ from immature to mature predator. Without the others stated, the rest of the parameters are chosen as follows:

We make the numerical simulations from the following four scenarios.

5.1.   Dynamics of the equilibrium states

For system (1.2), let $r = 0.08, f = k = \eta = c_2 = 0$ and other parameter values are the same as (5.1), then it is not difficult to find that all species will become extinct and the trivial equilibrium point is $E_0(0, 0, 0)$. Keep the parameters unchanged and $r = 3$, then the prey can survive while the predator remains extinct. The equilibrium state is $\overline{E}(29, 0, 0)$. By Theorem 3.1, both $E_0$ and $\overline{E}$ are LAS. See Figures 1 and 2, respectively.

Due to the complexity of the equations of equilibrium points, next we numerically validate that there exists a coexistence equilibrium point. Choose $f = 0.5, k = 1, \eta = 0.5, c_2 = 0.8$ and the rest of the parameters are taken from (5.1), then by computation we get the unique interior equilibrium point as $\widetilde{E}(3.6540, 0.9939, 4.2827)$. It is LAS by Theorem 3.1 (see Figure 3). That is to say, all the species can coexist together.

5.2.   Effects of functional response and stage structure

In order to find whether the functional response and the stage structure affect the system dynamics, we let $\eta = 0.9$ or $c_2 = 0.3$, and the other parameter values are identical with Figure 3. By simulation we get the time series graph of prey and predator species (see Figure 4). Comparing Figure 3 with Figure 4, we find that when the parameter values of $\eta$ or $c_2$ are changed, the interior equilibrium point of (1.2) becomes unstable from stable, i.e., the conversion rate of prey to predator as well as the stage structure of predators affects the system stability.

5.3.   Effect of fear and anti-predation sensitivity

Due to the predator's fear on prey, the prey presents some counter-predation features, then it is necessary to explore how the anti-predation sensitivity $k$ affects the stability of the system. In this subsection, we mainly pay our attention on the effect of parameters $k$ and $f$ on the system dynamics. By resetting the parameter values of $f = 2$ or $k = 0.2$ and keeping the rest the same as Figure 3, we get the time series diagrams of model (1.2) (Figure 5). Simulation figures reveal that system (1.2) becomes oscillatory accompanied with a period fluctuation, which means that the system bifurcation about parameters $k$ and $f$ may occur, so we further explore the system bifurcation phenomenon.

(1) Bifurcation analysis of parameter $k$.

In view of the importance of parameter $k$, which indirectly reveals how the fear level affects the whole system dynamics, we start with the bifurcation analysis of $k$.

(i) Bifurcation of $k$ for different fear level $f$.

For system (1.2) with the parameter values given in Figure 5(a), by use of Matlab and Moncont ^[43], we depict the continuous trajectories of the coexistence equilibrium point of prey with $k$ in different regions (see Figure 6). We find from Figure 6 that when $k = 0.002041$ or $k = 1.8719$, the Hopf bifurcation appears, and when $k = 0.002169$, the saddle-node bifurcation appears. That is, when $k$ is smaller, the system dynamics are more complex.

To better visualize the effect of $k$ on the stability for different fear $f$, we change the value of $f$ and draw the bifurcation picture of prey about parameter $k$ with $f = 0.5$ (see Figure 7). Compared with Figure 6(a), we find that the saddle-node bifurcation region is enlarged when the fear is reduced, which indicates that the fear indirectly affects the prey stability by the anti-predation sensitivity $k$.

Figures 6 and 7 show that if the fear level of prey is high but the anti-predation sensitivity is very low, then the system changes its stability from stable to unstable and even makes the prey extinct. When the anti-predation sensitivity is high enough, like $k > 1.8719$, then the system becomes stable again and the prey and predator can coexist. If the fear level is lower, then the system remains a stable and larger region. The survival region of prey is also enlarged.

(ii) Two parameter bifurcation.

$\blacklozenge$ The interplay of $k$ and $f$ on the dynamics of (1.2).

As previously stated, both the parameters $k$ and $f$ represent the effects of fear on the system dynamics, and we now study their joint effects. By Matcont, we get the bi-parameter bifurcation picture (see Figure 8). Figure 8(a) is with $k\in[0, 2]$ and Figure 8(b) is the enlarged part of Figure 8 (a) with $k\in[0, 0.8],$ where $R_1$ is the region of stable $\overline{E}$, $R_2$ is the region of bistability between stable $\overline{E}$ and oscillatory $\widetilde{E}$, $R_3$ is the region of bistability between stable $\overline{E}$ and $\widetilde{E}$, $R_4$ is the region of oscillatory $\widetilde{E}$, $R_5$ is the region of stable $\widetilde{E}$, $R_6$ is the region of bistability between two limit cycles and $R_7$ is the region of bistability between stable $\widetilde{E}$ and oscillatory $\widetilde{E}$. Figure 7 shows that as the changes of $k$ and $f$ occur, the cush point, homoclinic bifurcation curve (HC)), saddle-node bifurcation of limit cycle (LPC) and generalized Hopf point appear. The stability of $\widetilde{E}$ is changed and many bistability regions appear.

To better visualize the dynamics of system switching, keep the parameter values given in (5.1) unchanged. Take $f = 0.5, \eta = 0.5, c_2 = 0.8$ and change the value of $k$, or keep $k = 1, \eta = 0.5, c_2 = 0.8$ and change the value of $f$, then we draw the phase graph of the system with three groups of different initial data (see Figure 9) where the initial data for the blue line, green line and red line are (5, 4, 2), (1, 2, 0.5) and (1, 1, 4), respectively.

Figures 8 and 9 show that the fear level and anti-predation sensitivity can jointly affect the equilibrium stability of (1.2). Specifically, when the anti-predation sensitivity is low, the coexistence state of prey and predator will lose. The system will jump from one equilibrium point to another and oscillate between stability and instability frequently, even leading to bistability. The system dynamics are very complex.

$\blacklozenge$ The interplay of $k$ and $c_2$ on the dynamics.

As a stage structure system, the transversion rate is a very important indictor that decides on the whole system dynamics, so we investigate the interrelationship of anti-predation parameter $k$ and the transversion rate $c_2$. By simulation we get the graph of interplay effects between $k$ and $c_2$ on the system dynamics (see Figure 10) where $R_1$ is the region of stable $\widetilde{E}$, $R_2$ is the region of bistable between stable $\widetilde{E}$ and oscillatory $\widetilde{E}$ and $R_3$ is the region of oscillatory $\widetilde{E}$. Clearly, when $k$ is small, the dynamics of (1.2) are very rich. The cush point produces at the point $(k, c_2) = (0.008273, 0.355558)$ and the bistability regions between oscillatory $\widetilde{E}$ and stable $\widetilde{E}$ are presented. It shows that the anti-predation behaviors of prey affect the predation rate of the predator, resulting in a change of the transversion rate from immature predator to mature predator, and the system presents some complicated dynamics.

(2) Bubble phenomenon.

In the biological system dynamics, we often see a closed loop-like structure where oscillation appears from one Hopf bifurcation point and disappears from another Hopf bifurcation point. From the diagram it seems like a bubble in the real world, so we call it a bubble phenomenon ^[32,35]. By drawing the Hopf bifurcation picture w.r.t. fear level $f$, we see the interesting bubble phenomenon (see Figure 11). It shows that as the change of parameter $f$ occurs, the system (1.2) changes its stability from stable to unstable, then from unstable to stable again. That is, if the fear induced by the predator is low, then prey and predator can coexist and the system remains stable, but when it becomes large ($f > 0.649$), then the system becomes unstable. If the fear level is sufficiently large, for example, $f > 21.104$, then the prey and predator coexist again and the system remains stable.

(3) Bistability.

For a system with fixed parameter values, if the initial data is changed, then the solution curves of the system is convergent with different attractors and we say the system has a bistability phenomenon. From the bifurcation graphs of two parameters (Figures 8 and 10), we find that (1.2) has some different bistability behaviors. Given some different initial data, we depict the diagrams of time series and phase of system (1.2) to better visualize the bistability phenomenon (see Figure 12) where the parameter values are same as given in (5.1), and $f = 2, k = 1.87, \eta = 0.5$ and $c_2 = 0.8.$ Figure 12(a) shows that the trajectories with different initial data converge to different cycles respectively. The trajectories of Figure 12(c) corresponding to different initial data converges to a cycle and a stable equilibrium respectively. However, the trajectories of Figure 12(e) converge to the same cycle although the initial data is different. Biologically, it means that under different conditions, the prey and predator can experience multiple stable coexistence states, which is valuable for people to keep the species alive and serve the economic development.

In a word, the fear of predator brings some important influences to system (1.2), which leads to richer dynamical properties of the system. It makes the prey produce anti-predation behaviors and affects the predation rate of the predator, as well as the transversion rate of mature predator from immature predator. Particularly, the change of anti-predation sensitivity brings large influences to the system dynamics, such as changing the system stability from stable to unstable or from unstable to stable, accompanied by bifurcations, limit cycles and the bubble or bistability phenomena.

5.4.   Effect of time delays

Now we investigate how the time delay affects the system dynamics of system (1.3). For system (1.3), we keep the parameter values the same as given in Figure 3; that is, it is LAS when the delays are absent. Then, by choosing two small time delays as $\tau_1 = \tau_2 = 0.08$, we find it is still stable (see Figure 13(a) and (b)). However, if we increase the delay values, for example, let $\tau_1 = \tau_2 = 5$, then we observe that it is unstable and presents periodic fluctuation behaviors (see Figure 13(c) and (d)).

By the simulation Figure 13, we conjecture that there may be a Hopf bifurcation. First, we consider the case of one delay. Let $\tau_2 = 0$ (or $\tau_1 = 0$) and depict the bifurcation picture w.r.t. $\tau_1$ (or $\tau_2$), then we get the bifurcation pictures (see Figure 14(a) or Figure 14(b)). Second, we consider the case of two delays. Keep one delay fixed in the range of assuring the system stability, and take the other delay as bifurcation parameter, then we get the delay bifurcation pictures (see Figure 14(c) and (d)). Similarly, as the change of one delay, Hopf bifurcations occur and system stability are changed from stable to unstable. If they are large enough, then the system will become unstable and produce chaos (see Figure 14(e) and (f)).

All simulation figures indicate that if the prey can respond timely to the fear of predator and the mature period of the predator is very short, then the system can remain stable. Meanwhile, if the fear delay or mature delay is large enough, then the system stability is lost and accompanied with periodic fluctuation. If one of them becomes too large, then chaos occurs and the system becomes unstable. In a word, the fear and mature delays bring some significant influences to the system stability.

6.   Discussion and conclusions

Combining the real ecological environment, we formulated a delayed predator-prey system with a Crowley-Martin functional response and stage structure for predators. The existence of equilibrium points, the locally asymptotical stability and bifurcation behaviors of the system were investigated. All theoretical results were summarized in Theorems 3.1 and 4.1–4.4. In order to validate the theoretical results and intuitively reveal the influence on system dynamics of fear, stage structure, anti-predation sensitivity and time delays, some numerical analyses were given. Figures 1–3 validated the stability of equilibrium points. Figure 4 showed the functional response and stage structure of how predators affected the system stability. Figures 5–12 revealed how the fear and anti-predation sensitivity affected the system stability numerically. Figures 13 and 14 explained the effects of fear delay and mature delay on the system stability. We summarize our findings as follows:

$(1)$ The conversion rate from prey to predator and the stage structure of predator bring some influence on the stability of the equilibrium point of the system.

$(2)$ The fear induced by the predator seriously affects the stability of the interior equilibrium point. It makes the prey produce anti-predator behaviors, which affects the predation rate and transversion rate of mature predator from immature predator. The fear and anti-predation sensitivity jointly affects the system stability and produces complex dynamical phenomena such as bubble, or bistability.

$(3)$ The fear delay and mature delay will bring some important influence on the system dynamics. They can change the system stability from stable to unstable and even lead to chaos if they are large enough.

By these findings, from the angle of protecting the species from extinction to maintaining the system permanence and the biological balance, we should make some measures such as

● Constructing some refuge zones for prey to prevent the fear of predators and increase the anti-predation sensitivity of prey.

● Supply some additional food for predators and decrease the predation rate to decline the negative effect of predators fear on prey.

● Shorten the mature delay and accelerate the mature speed of predators from immature predator by using some biological control strategies, like changing their genes.

There are still some topics to be investigated. For example, the environments are always fluctuated by random noise, so incorporating the effects of environmental stochastic noise on the system dynamics is very interesting^[44]. On the other hand, how the refuge and additional food supplement affects the stability is another interesting topic^[45,46]. In addition, whether the mature delay can be efficiently controlled by some reasonable ways is also necessary to be considered. All of these are left for our future research.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research is supported by the National Science Foundation of China (No. 11861027, No. 12161062) and the Science and Technology Program of Inner Mongolia Autonomous Region (No. 2023YFHH0024, No. 2022YFHH0017, No. 2022YFHH0063).

The authors would like to thank the editor and the anonymous reviewers for their valuable comments, which helped to improve the presentation of the paper.

Conflict of interest

There is no competing interests in this paper.