Dynamics of a predator-prey model with fear effects and gestation delays

1.   Introduction

Predator-prey interactions play an important role in shaping community structure and maintaining ecological diversity. Over the past few decades, food chain models have been extensively studied and a large number of papers have been published. Researchers have used mathematical models to understand the complexity of ecosystems and predict the consequences of ecosystems ^{[1,2,3,4,5,6,7,8]}.

Generally, predators affect prey populations by killing them directly, but many recent studies ^[9,10,11] have shown that fear effects play a crucial role in ecology and evolutionary biology. In ^[9], fear was found to cause prey species to exhibit many anti-predator behaviors and defense mechanisms. In ^[10], prey birth and death rates were affected by fear of predation risk. In ^[11], higher levels of fear were found to result in lower equilibrium levels of coexistence and global asymptotic stability of species near equilibrium. Their results all point to the importance of fear in predator-prey systems. The presence of predators forces prey to be more alert and spend less time foraging or choosing new habitats due to the fear of being killed, while fear also reduces the mating rate of prey populations, resulting in a decrease in prey population size. In 2011, Zanette et al. ^[12] used song sparrows as subjects for their experiments. They left the population of song sparrows in a predator-free area and then manipulated artificial predation risk by playing the sounds of predators. Through experiments, they observed that the song sparrows' perception of predation risk led to a reduction in their breeding population by about $40\%$ compared to the previous population. So this experiment confirms that fear reduces the reproductive rate of prey. In another experiment involving elk, Creel et al. ^[13] observed that the risk of predation by wolves caused elk to remain alert and move to areas with low predation risk. This suggests that the physiological behavior of elk is influenced by anti-predator behavior. Suraci et al. ^[14] used sound playback of large carnivores to induce fear in medium carnivores and found that fear of large carnivores increased vigilance and greatly reduced foraging behavior in medium carnivores. These experiments all confirm the important effects of fear on the dynamic behavior of species. In recent years, many researchers have proposed that the predator-prey system should incorporate a fear factor. In 2016, Wang et al. ^[15] proposed a mathematical model with fear and confirmed the importance of fear on the growth of prey.

In ecology, how predators feed on prey is a very important factor, and the way they interact with each other affects the dynamics of the entire ecosystem. This interaction is called the "functional response", which gives the per capita rate of the predator feeding on prey and is a key component of the mathematical model, allowing us to study the dynamics of the model in a better way. The most widely used functional responses include Holling type Ⅱ, Holling type Ⅲ, Holling type Ⅳ ^[16,17,18], Monod-Haldane functional response, Beddington-DeAngelis functional response, Crowley-Martin functional response and the ratio-dependent functional responses, etc. ^{[19,20,21,22,23]}. The Holling type Ⅱ functional response is the more commonly used functional response, but this functional response depends only on prey density and is not applicable for ecosystems with large population sizes. To achieve the elimination of this situation, we introduce a ratio-dependent functional response, which is determined simultaneously by the relationship between prey density and predator density. In particular, when predators have to go in search of food, they will face a situation of sharing and competition with other predators. Therefore, ratio-dependent functional responses will be more closely related to practicality. Also, considering that predators are disturbed by other predators while handling or searching for prey, the Beddington-DeAngelis functional response deserves further study. Therefore, we will consider introducing ratio-dependent functional response and Beddington-DeAngelis functional response in the model.

Delay differential equations (DDEs) have a long history in modeling predator-prey systems, where delays usually represent the required reaction time, gestation period, maturation, feeding time, etc. ^{[24,25,26,27,28,29]}. In terms of mathematical modeling, delayed systems ^[30,31] make more sense than non-delayed systems. The introduction of time delays in a predator-prey model causes the model to exhibit more complex dynamics. Among them, the gestation delay represents the time delay that occurs when biomass is transferred from prey to predator. For example, Tripathyi et al. ^[5] investigated the degree of mutual interference between density-dependent predators in a model incorporating discrete gestation delays and Beddington-DeAngelis-type functional responses. They observed that when the delay parameter exceeded a certain threshold, the system experienced Hopf bifurcation and the system lost its stability. Based on this finding, they investigated two different time delay effects in a proportionally dependent predator-prey system and studied the conditions under which Hopf-bifurcation occurred when both delays were present in the model. From the above example, we can determine the importance of delay in stability analysis. Therefore, we add the time lag parameters to the model.

In this paper, we aim to study the effect of time lag parameters on the dynamics of the system. The rest of the paper is organized as follows. In the next section, we construct a predator-prey model with fear effects and multiple delays. In Section 3, we analyze the positivity and boundedness of the solutions of the system. In Section 4, we present the criterion of the existence of equilibrium points. In Section 5, we study the stability of the equilibrium point and the conditions of the occurrence of Hopf bifurcation with time lags as the bifurcation parameter. In Section 6, we use numerical simulations to verify our main analytical results. Finally, we conclude the whole paper.

2.   Model formulation

Let $x_1(t)$, $x_2(t)$, and $x_3(t)$ be the population densities of prey, intermediate predator, and top predator at time $t$, respectively. We construct the mathematical model through the following assumptions.

$(1)$ In the absence of predators, the prey population increases logically. However, in the presence of intermediate predators, the prey perceives fear, and fear affects the reproductive rate of the prey.

$(2)$ Intermediate predators prey on preys based on a ratio-dependent functional response. However, the interaction between intermediate and top predators follows a Beddington-DeAngelis type functional response. Among them, the expression form of the ratio-dependent functional response is $f(x_1, x_2) = \frac{\alpha x_1}{x_1+mx_2}$, and the expression form of the Beddington-DeAngelis type functional response is $f(x_2, x_3) = \frac{\theta x_2}{v+px_2+qx_3}$.

$(3)$ Most biological processes involve delays, and here we consider two gestation delays, a time interval between the predator preying on its prey and the prey breeding offspring, denoted by $\tau_1$, $\tau_2$.

Based on these assumptions, we obtain the following model:

where $r, k, b, a$ and $\alpha$ are the intrinsic growth rate of prey, intermediate predator-induced fear, mortality, intraspecific competition rate and prey consumption rate, respectively. Meanwhile, $\beta, d, \theta, \eta$ and $e$ are prey conversion rate, intermediate predator mortality and consumption rate, intermediate predator conversion rate and top predator mortality rate, respectively. Here $m$ is the half-saturation constant of the ratio-dependent functional response, $v$ is the saturation constant of the Beddington-DeAngelis-type functional response. $p$ and $q$ are the abundance of intermediate predators and disturbance of top predators, respectively. $\tau_1$ and $\tau_2$ represent gestation time delays.

The initial conditions are

where $\tau = \max \left\lbrace \tau_1, \tau_2 \right\rbrace$, $\mathit{\Phi} = (\phi _1, \phi_2, \phi_3)$, and $\mathit{\Phi} = (\phi _1, \phi_2, \phi_3)$ belongs to the Banach space of continuous function $\mathit{\Phi}:[-\tau, 0]\to \mathit{R^3}$ with

From the biological perspective, we assume $\phi _i(s)\ge 0$ and $\phi _i(0) > 0, \ i = 1, 2, 3$.

3.   Positivity and boundedness of solutions

Theorem 3.1. All solutions of system (2.1) with initial conditions (2.2) are positive for all $t > 0$.

Proof. System (2.1) can be written in the following matrix form:

where $W = \begin{pmatrix} x_1\\ x_2\\ x_3 \end{pmatrix}$, and $F(W) = \begin{pmatrix} F_1(W)\\ F_2(W)\\ F_3(W) \end{pmatrix} = \begin{pmatrix} \frac{rx_1}{1+kx_2}-bx_1-ax_1^2-\frac{\alpha x_1x_2}{x_1+mx_2}\\ \frac{\beta x_1(t-\tau_1)x_2(t-\tau_1)}{x_1(t-\tau_1)+mx_2(t-\tau_1)}-dx_2-\frac{\theta x_2x_3}{v+px_2+qx_3}\\ \frac{\eta x_2(t-\tau_2)x_3(t-\tau_2)}{v+px_2(t-\tau_2)+qx_3(t-\tau_2)}-ex_3 \end{pmatrix}$ with the initial condition

We have easily checked in system (2.1) that whenever choosing $W(s)\in \mathit{R^3_+}$, then

Using the Lemma $4$ given in ^[32], we obtain that every solution $W(t)$ of the system (2.1) with initial conditions (2.2) exists in the region $\mathit{R^3_+}$ and is positive and invariant for all $t > 0$.

Theorem 3.2. The solutions of system (2.1) which starts in $\mathit{R^3_+}$ are uniformly bounded.

Proof. Let $U(t) = x_1(t-\tau_1-\tau_2)+\frac{\alpha }{\beta }x_2(t-\tau_2)+\frac{\alpha \theta }{\beta \eta } x_3(t)$.

Calculating the derivative of $U(t)$ along the solution of system (2.1), we have

Let $z = \min \left \{ b, d, e \right \}$ and $Q = \frac{r^2}{4a}$. Using the theory of differential inequality ^[33], we obtain

When $t\longrightarrow \infty$, we get $0 < U\le \frac{Q}{z}$. Therefore, all the solutions of system (2.1) starting from $R^3_+$ are restricted in the region:

This completes the proof.

4.   Equilibrium points and their criteria of existence

The possible equilibrium points of the model are as follows:

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

$(2)$ The predator free axial equilibrium point $E_1\left(\frac{r-b}{a}, 0, 0\right)$ exists if $r > b$.

$(3)$ The top predator free equilibrium point $E_2(\bar{x}_1, \bar{x}_2, 0)$ where $\bar{x}_2 = \frac{\beta -d}{md}\bar{x}_1$ and $\bar{x}_1$ is a positive root of $-ac(1+cm)x_1^2-\left [(1+cm)(cb+a)+c^2\alpha \right]x_1+(r-b)(1+cm)-c\alpha = 0$, where $c = \frac{\beta -d}{md}$. According to Descarte's rule of sign, the equation must have a unique positive root if $(r-b)(1+cm) > c\alpha$. Thus, the equilibrium $E_2$ exists if $(r-b)(1+cm) > c\alpha$.

$(4)$ The interior equilibrium point $E^*(x_1^*, x_2^*, x_3^*)$ exists if the following two isoclines have a positive intersection:

where

As $x_2\longrightarrow 0$, then the above two isoclines become:

where the number of positive roots of the $f(x_1, 0) = 0$ depends on the sign of $r-b$, then the above two isoclines have a unique intersection on the first quadrant (the positive quadrant) if the following sufficient condition holds:

5.   Local stability and Hopf bifurcation

In this section, our goal is to determine the local behavior of the system (2.1) around the equilibrium point. First, system (2.1) can be expressed as

where $W(t) = [x_1(t), x_2(t), x_3(t)]^T$.

We take $E(x_1, x_2, x_3)$ as an arbitrary equilibrium point of the system (2.1). Assume $\tilde{x}_1(t) = x_1(t)-x_1$, $\tilde{x}_2(t) = x_2(t)-x_2$ and $\tilde{x}_3(t) = x_3(t)-x_3$. Then, the linearized system (2.1) around the equilibrium point $E$ is

where

and $G(t) = [\tilde{x}_1(t), \tilde{x}_2(t), \tilde{x}_3(t)]^T$. The Jacobian matrix of the system (2.1) at $E$ is given by

Then

where

According to the Jacobian matrix, we analyze the stability of the equilibrium point $E_0$, $E_1$ and $E_2$, and obtain the following theorems:

Theorem 5.1. The trivial equilibrium point $E_0(0, 0, 0)$ is locally asympototically stable if $r < b$ and becomes unstable if $r > b$.

Theorem 5.2. When $\tau_1\ge 0$, the predator free axial equilibrium point $E_1\left(\frac{r-b}{a}, 0, 0\right)$ is locally asympototically stable if $\beta < d$; otherwise it is unstable.

Theorem 5.3. For system (2.1), according to the value of $\tau_1$, there are the following two cases:

(1) When $\tau_1 = 0$, the top predator free equilibrium point $E_2 = (\bar{x}_1, \bar{x}_2, 0)$ is locally asympototically stable for any $\tau_2\ge 0$ if $\bar{x}_2 < \frac{ve}{\eta -pe}$, $a\bar{x}_1+d+\frac{\beta \bar{x}_1^2}{(\bar{x}_1+m\bar{x}_2)^2} > \frac{\alpha \bar{x}_1\bar{x}_2}{(\bar{x}_1+m\bar{x}_2)^2}$ and $\left (-a\bar{x}_1+\frac{\alpha \bar{x}_1\bar{x}_2}{(\bar{x}_1+m\bar{x}_2)^2} \right)\left(\frac{\beta \bar{x}_1^2 }{(\bar{x}_1+m\bar{x}_2)^2} -d \right)$$+\frac{\beta m\bar{x}_2^2}{(\bar{x}_1+m\bar{x}_2)^2}\left(\frac{kr\bar{x}_1}{(1+k\bar{x}_2)^2}+\frac{\alpha \bar{x}_1^2}{(\bar{x}_1+m\bar{x}_2)^2} \right) > 0$; otherwise unstable.

(2) When $\tau_1 > 0$ and $\bar{x}_2 < \frac{ve}{\eta -pe}$, if $H_2^2-N_2^2 < 0$ and $2\varphi _0^2+H_1-N_1^2-2H_2 > 0$ hold, the super predator free equilibrium point $E_2 = (\bar{x}_1, \bar{x}_2, 0)$ is locally asympototically stable for $\tau_1\in \left [0, \hat{\tau_1} \right)$ and unstable for $\tau_1 > \hat{\tau_1}$. Moreover, the system (2.1) undergoes a Hopf bifurcation when $\tau_1 = \hat{\tau_1}$, where

and $\varphi _0$ is defined in the proof.

Next, we study the local stability of $E^* = (x_1^*, x_2^*, x_3^*)$. The Jacobian matrix at $E^*$ is

where

The corresponding characteristic equation of the above Jacobian matrix is

where

We divide the following analysis into five cases according to the different ranges of values of the two time delays. Based on the characteristic equation (5.1), we analyze the stability of the interior equilibrium point $E^*$ and the conditions for the occurrence of Hopf-bifurcation, and obtain the following theorems.

Case 1: $\tau_1 = \tau_2 = 0$.

Theorem 5.4. In absence of both delays, the interior equilibrium point $E^*$ is locally asymptotically stable if $P_1+Q_1+R_1 > 0$, $P_3+Q_3+R_3+S_2 > 0$ and $(P_1+Q_1+R_1)(P_2+Q_2+R_2+S_1) > P_3+Q_3+R_3+S_2$.

Case 2: $\tau_1 > 0, \tau_2 = 0$.

Theorem 5.5. For system (2.1), when $\tau_1 > 0, \tau_2 = 0$, if $(P_3+R_3)^2 < (Q_3+S_2)^2$ holds, the interior equilibrium point $E^*$ is locally asymptotically stable for all $\tau_1\in \left [0, \tau_{10} \right)$ and unstable for $\tau_1 > \tau_{10}$. Furthermore, when $3\varphi _{0}^4+(2E_{11}^2-4E_{12}-2F_{11}^2)\varphi _{0}^2+(E_{12}^2-2E_{11}E_{13}+2F_{11}F_{13}-F_{12}^2) > 0$, then Hopf-bifurcation occurs at $\tau_1 = \tau_{10}$, where

and $\varphi _{k}$ are defined in the proof.

Case 3: $\tau_1 = 0, \tau_2 > 0$.

Theorem 5.6. For system (2.1), when $\tau_1 = 0, \tau_2 > 0$, if $(P_3+Q_3)^2 < (R_3+S_2)^2$ holds, the interior equilibrium point $E^*$ is locally asymptotically stable for all $\tau_2\in \left [0, \tau_{20} \right)$ and unstable for $\tau_2 > \tau_{20}$. Furthermore, when $3\varphi_{0}^4+(2E_{21}^2-4E_{22}-2F_{21}^2)\varphi_{0}^2+(E_{22}^2-2E_{21}E_{23}+2F_{21}F_{23}-F_{22}^2) > 0$, then Hopf-bifurcation occurs at $\tau_2 = \tau_{20}$, where

and $\varphi _{k}$ are defined in the proof.

Case 4: $\tau_1 > 0, \tau_2\in (0, \tau_{20})$.

Theorem 5.7. For system (2.1), when $\tau_1 > 0, \tau_2\in (0, \tau_{20})$, if $(P_3+R_3)^2 < (Q_3+S_2)^2$ holds, the interior equilibrium point $E^*$ is locally asymptotically stable for all $\tau_1\in \left [0, \hat{\tau}_{10} \right)$ and unstable for $\tau_1 > \hat{\tau}_{10}$. Furthermore, when $M_1M_3+M_2M_4 > 0$, then Hopf-bifurcation occurs at $\tau_1 = \hat{\tau}_{10}$, where

and $\varphi _{k}$ are defined in the proof.

Case 5: $\tau_1\in (0, \tau _{10}), \tau_2 > 0$.

Theorem 5.8. For system (2.1), when $\tau_1\in (0, \tau _{10}), \tau_2 > 0$, if $(P_3+Q_3)^2 < (R_3+S_2)^2$ holds, the interior equilibrium point $E^*$ is locally asymptotically stable for all $\tau_2\in \left [0, \hat{\tau} _{20} \right)$ and unstable for $\tau_2 > \hat{\tau} _{20}$. Furthermore, when $M_5M_7+M_6M_8 > 0$, then Hopf-bifurcation occurs at $\tau_2 = \hat{\tau} _{20}$, where

and $\varphi _{k}$ are defined in the proof, where

Remark 1. The proofs of above theorems are given in the Appendix.

6.   Numerical simulations

In this section, we use MATLAB R2020a to perform some numerical simulations to illustrate the results of the analysis in the previous sections and plot the corresponding graphs of system (2.1). Since the problem is not a species-specific case study and no real data are available, some hypothetical data are taken here for the simulation.

We make the parameters values $r = 4.5$, $k = 6$, $b = 0.2$, $a = 0.06$, $\alpha = 1$, $m = 1$, $\beta = 1.3$, $d = 0.03$, $\theta = 4$, $v = 3.1$, $p = 1.14$, $q = 1.13$, $\eta = 2.2$, $e = 0.2$, and the initial value is $(\phi_1(s), \phi_2(s), \phi_3(s)) = (2, 2, 2), \ s\in[-\tau, 0]$. According to the above parameters, we can obtain all equilibrium points, respectively $E_0(0, 0, 0)$, $E_1(71.67, 0, 0)$, $E_2(0.0111, 0.4699, 0)$ and $E^*(14.5, 1.73, 0.51)$. In the following, we mainly show several cases of interior equilibrium points.

For $\tau_1 = 0, \tau_2 = 0$, the interior equilibrium point $E^*$ is locally asymptotically stable (see Figure 1).

For $\tau_1 > 0, \tau_2 = 0$, we plotted the Hopf-bifurcation diagram Figure 2 concerning $\tau_1$. If we change the value of $\tau_1$ from $0$ to $14$, then the system undergoes Hopf-bifurcation. Hence, the interior equilibrium $E^*$ is locally asymptotically stable for $\tau_1 = 1 < \tau_{10}\approx 10.88$ (see Figure 3) and unstable for $\tau_1 = 13 > \tau_{10}\approx 10.88$ (see Figure 4).

For $\tau_1 = 0, \tau_2 > 0$, if we continuously increase the value of $\tau_2$, we can find a critical value of $\tau_{20}$, namely $\tau_{20}\approx 0.11$ for which system (2.1) undergoes Hopf-bifurcation (see Figure 5. Hence, $E^*$ is locally asymptotically stable for $\tau_2 = 0.1 < \tau_{20}\approx 0.11$ (see Figure 6) and unstable for $\tau_2 = 0.2 > \tau_{20}\approx 0.11$ (see Figure 7).

For $\tau_1 > 0$, fixing $\tau_2 = 0.1\in (0, \tau_{20})$, we get $\hat{\tau} _{10}\approx 10.28$, so when $\tau_1 = 2 < \hat{\tau} _{10}$, $E^*$ is locally asymptotically stable (see Figure 8), when $\tau_1 = 12 > \hat{\tau} _{10}$, $E^*$ is unstable (see Figure 9). Then Hopf-bifurcation occurs at $\hat{\tau} _{10}\approx 10.28$ (see Figure 10).

For $\tau_2 > 0$, fixing $\tau_1 = 1\in (0, \tau_{10})$, we get $\hat{\tau} _{20}\approx 0.96$, so when $\tau_2 = 0.1 < \hat{\tau} _{20}$, $E^*$ is locally asymptotically stable (see Figure 11), when $\tau_2 = 1.1 > \hat{\tau} _{20}$, $E^*$ is unstable (see Figure 12). Then Hopf bifurcation occurs at $\hat{\tau} _{20}\approx 0.96$ (see Figure 13).

Based on above analysis, we summarize the dynamics of the interior equilibrium $E^*$ of system (2.1) in Table 1.

Finally we investigate how the fear of predators affect the population dynamics. By increasing the fear values and taking $k = 6, 9, 15$ and 20 respectively, we plot the time series graphs of each species (see Figures 14–16).

(1) The effect of fear on prey. In Figure 14, we observe that fear inhibits prey reproduction and excessive fear leads to prey extinction.

(2) The effect of fear on intermediate predators. In Figure 15, we observe that an increase in the value of fear leads to a shorter duration of intermediate predator turbulence, but little change in the number of species, and too much fear leads to a decrease and stabilization of the number of species.

(3) The effect of fear on the top predator. In Figure 16, an increase in fear leads to a decrease in the turbulence time of the highest predator and a decrease in the number of species, and too much fear leads the highest predator to extinction.

7.   Discussion and conclusions

In this paper, we analyze the dynamics of a three-chain model containing indirect predation (fear of predator) and two delays, where the three species are the prey, intermediate predator and top predator, and the two delays represent the gestation delays of the intermediate and top predator, respectively.

First, we investigate the positivity and boundedness of this system, where boundedness can be seen as a natural limit to expansion due to limited resources and positivity implies that species persist. Next, we study the effect of delay on the stability of the model by varying the delay parameters $\tau_1$ and $\tau_2$. According to the theorem in Section 5, we obtain the conditions when the system is stable. By selecting appropriate parameter values, we draw the Hopf-bifurcation diagrams (Figures 2, 5, 10 and 13), that is, if $\tau$ increases continuously, a threshold $\tau_0$ will be obtained. If $\tau < \tau_0$, the system is stable; otherwise, the system becomes unstable, that is, a Hopf-bifurcation occurs at $\tau = \tau_0$. These observations confirm the important role of delay in the system. Through numerical simulations, we find that fear may reduce the number of species and even lead to extinction of species. Thus, fear has an important influence on the dynamics of predator-prey systems.

In nature, in addition to fear, other factors such as refuge, additional food, and human capture influence the predator-prey system. As a part of future research on the model in this paper, the inclusion of refuge and human capture could be considered to make the model more realistic. Moreover, for the discrete-time model ^[6,7] and the infectious disease model ^[8], we believe there will be some intresting findings. All these will be left to our future work.

Acknowledgments

The authors would like to thank the referees for their suggestions and valuable comments, which have greatly impoved the presentation of this paper.

This work is funded by NSF of China (No. 11861027).

Conflict of interest

The authors declare that they have no competing interests in this paper.

Appendix

Here, we give the proofs of Theorems 5.1 to 5.8.

Proof of Theorem 5.1:

Proof. Based on the eigenvalues less than $0$, we obtain that $E_0$ is locally asymptotically stable when $r < b$.

Proof of Theorem 5.2:

Proof. By calculating the Jacobian matrix at the equilibrium $E_1$, we obtain

Obviously, we get $\lambda _1 = -(r-b) < 0$, $\lambda _2 = -e < 0$, and the other eigenvalue is the root of the following equation:

In the absence of delay, i.e. when $\tau_1 = 0$, $E_1$ is locally asymptotically stable if $\beta < d$.

When $\tau_1 > 0$ and $\beta < d$, we assume that $\lambda = \xi +i\varphi (\xi, \varphi \in \mathit{R})$ is the root of Eq (A.1). Then substituting the value of $\lambda$ into Eq (A.1) and separating the real and imaginary parts, we obtain

To eliminate $\tau_1$, we square and add the Eqs (A.2) and (A.3), and obtain the following equation

Let $\xi \ge 0$, then $\left (\beta e^{-\xi \tau_1}-d-\xi \right)\le \left (\beta -d-\xi \right) < 0$. Thus, if $\xi \ge 0$, no such real $\varphi$ exists, which contradicts the previous assumption $\varphi\in \mathrm {R}$. So $\xi < 0$, i.e., the Eq (A.1) contains a negative real root and an imaginary root with a negative real part (if exists).

Therefore, $E_1$ is locally asymptotically stable if $\beta < d$.

Proof of Theorem 5.3:

Proof. The Jacobian matrix at $E_2$ is:

The characteristic equation of above matrix $J_{E_2}$ is given by $det(J_{E_2}-\lambda I) = 0$.

One root of $det(J_{E_2}-\lambda I) = 0$ is given by $-e+\frac{\eta \bar{x}_2}{v+p\bar{x}_2}e^{-\lambda \tau_2}-\lambda = 0$. The analysis process is similar to Eq (A.1). Therefore, if $-e+\frac{\eta \bar{x}_2}{v+p\bar{x}_2} < 0$, i.e., $\bar{x}_2 < \frac{ve}{\eta -pe}$, then the equation contains negative real roots and imaginary roots with negative real parts (if exists).

The other two eigenvalues are the root of the quadratic equation

where

Case 1: $\tau_1 = 0$, then Eq (A.4) becomes

According to Routh-Hurwitz criteria, both roots of the Eq (A.5) are negative real parts if $(H_2+N_2) > 0$ and $(H_1+N_1) > 0$, i.e. $\left (-a\bar{x}_1+\frac{\alpha \bar{x}_1\bar{x}_2}{(\bar{x}_1+m\bar{x}_2)^2} \right)\left(\frac{\beta \bar{x}_1^2 }{(\bar{x}_1+m\bar{x}_2)^2} -d \right)+\frac{\beta m\bar{x}_2^2}{(\bar{x}_1+m\bar{x}_2)^2} \left(\frac{kr\bar{x}_1}{(1+k\bar{x}_2)^2}+\frac{\alpha \bar{x}_1^2}{(\bar{x}_1+m\bar{x}_2)^2} \right) > 0$ and $a\bar{x}_1+d+\frac{\beta \bar{x}_1^2}{(\bar{x}_1+m\bar{x}_2)^2} > \frac{\alpha \bar{x}_1\bar{x}_2}{(\bar{x}_1+m\bar{x}_2)^2}$.

Case 2: $\tau_1 > 0$ and $\bar{x}_2 < \frac{ve}{\eta -pe}$. By substituting $\lambda = \xi +i\varphi$ into Eq (A.4) and separating the real and imaginary parts, we have

Putting $\xi = 0$, Eqs (A.6) and (A.7) become

By squaring and adding, we obtain

Using Descarte's rule of sign, the equation has at least one positive root $\varphi _0$ if $H_2^2-N_2^2 < 0$. By calculating (A.8) and (A.9), we get

Let $\hat{\tau_1} = \underset{j = 0, 1, 2...}{min}{\tau_{1j}}$.

Now the transversality condition $Re\left (\frac{\mathrm{d} \lambda }{\mathrm{d} \tau_1}\right)^{-1}_{\lambda = i\varphi _{0}} > 0$ will be verified.

By differentiating Eq (A.4) with respect to $\tau_1$, we obtain

which leads to

From (A.10), we have

Then

When $2\varphi _0^2+H_1-N_1^2-2H_2 > 0$, $Re\left (\frac{\mathrm{d} \lambda }{\mathrm{d} \tau_1}\right)_{\lambda = i\varphi _0}^{-1} > 0$. Therefore, the transversality condition is satisfied and a Hopf-bifurcation occurs at $E_2$ for $\tau_1 = \hat{\tau_1}$.

Proof of Theorem 5.4:

Proof. The characteristic (5.1) becomes

Therefore, by Routh-Hurwitz Criteria, we obtain all the roots of (A.12) have negative real part, if $P_1+Q_1+R_1 > 0$, $P_3+Q_3+R_3+S_2 > 0$ and $(P_1+Q_1+R_1)(P_2+Q_2+R_2+S_1) > P_3+Q_3+R_3+S_2$. This means that the interior equilibrium $E^*$ is locally asymptotically stable.

Proof of Theorem 5.5:

Proof. The characteristic (5.1) becomes

where $E_{11} = P_1+R_1$, $E_{12} = P_2+R_2$, $E_{13} = P_3+R_3$, $F_{11} = Q_1$, $F_{12} = Q_2+S_1$, $F_{13} = Q_3+S_2$. Let $\lambda = \xi+i\varphi$ be the root of (A.13). Then we obtian

A necessary condition for changing the stability of the equilibrium point $E^*$ is that the real part of the root of the Eq (A.13) changes the sign. To find the stable switching point, we consider $\xi = 0$, then we have

By squaring and adding (A.14) and (A.15), we obtain

where $K_{11} = E_{11}^2-2E_{12}-F_{11}$, $K_{12} = E_{12}^2-2E_{11}E_{13}+2F_{11}F_{13}-F_{12}^2$, $K_{13} = E_{13}^2-F_{13}^2$.

Let $\varphi ^2 = \sigma$, we get

Then, $h_1(0) = K_{13} = E_{13}^2-F_{13}^2 = (P_3+R_3)^2-(Q_3+S_2)^2$ and $\lim_{\sigma \to \infty}h_1(\sigma) = + \infty$. We assume that $h_1(0) < 0\Rightarrow (P_3+R_3)^2 < (Q_3+S_2)^2$. According to Descarte's rule of sign, the Eq (A.17) has at least one positive root and can have at most three positive roots.

Without loss of generality, we assume that it has three positive roots, denoted by $\varphi _{k} = \sqrt{\sigma _{k}}$ and $k = 1, 2, 3$. By calculating (A.14) and (A.15), we obtain

where $k = 1, 2, 3$ and $j = 0, 1, 2, ...$. Let $\tau_{10} = \underset{k, j}{min}\tau_{1k}^{(j)}, \ k = 1, 2, 3, \ j = 0, 1, 2...$.

Now the transversality condition $Re\left (\frac{\mathrm{d} \lambda }{\mathrm{d} \tau_1}\right)^{-1}_{\lambda = i\varphi _{0}} > 0$ will be verified.

By differentiating Eq (A.13) with respect to $\tau_1$, we obtain

which leads to

From (A.16), we have

Then

When $3\varphi _{0}^4+(2E_{11}^2-4E_{12}-2F_{11}^2)\varphi _{0}^2+(E_{12}^2-2E_{11}E_{13}+2F_{11}F_{13}-F_{12}^2) > 0$, $Re\left (\frac{\mathrm{d} \lambda }{\mathrm{d} \tau_1}\right)^{-1}_{\lambda = i\varphi _{0}} > 0$. Therefore, the transversality condition is satisfied and a Hopf-bifurcation occurs around when $\tau_1$ passes through the critical value $\tau_{10}$.

Proof of Theorem 5.6:

Proof. The characteristic (5.1) becomes

where $E_{21} = P_1+Q_1$, $E_{22} = P_2+Q_2$, $E_{23} = P_3+Q_3$, $F_{21} = Q_1$, $F_{22} = R_2+S_1$, $F_{23} = R_3+S_2$.

The calculation process is similar to Theorem 5.5.

Proof of Theorem 5.7:

Proof. We consider that $\tau_2 = \tau_2^0$ in its stable interval, with $\tau_1$ as the parameter.

Let $\lambda = \xi+i\varphi$ be the root of (5.1), we get

A necessary condition for the stability change of the equilibrium point $E^*$ is that the characteristic Eq (5.1) has a pair of purely imaginary roots, i.e., $\xi = 0$. Let $\xi = 0$, and we obtain

Squaring and adding the two equations, and we have

where

Denote

Obviously, $h_2(0) = K_{23}+K_{25}\mid _{\varphi = 0} = (P_3+R_3)^2-(Q_3+S_2)^2$ and $\lim_{\varphi \to \infty}h_2(\varphi) = + \infty$. We assume that $h_2(0) < 0\Rightarrow (P_3+R_3)^2 < (Q_3+S_2)^2$. Then, the Eq (A.21) has at least one positive root.

Assume that Eq (A.21) has finite positive roots $\varphi _{k}(k = 1, 2, ..., l)$. The critical value

where

Let $\hat{\tau}_ {10} = \underset{k, j}{min}\tau_{1k}^{(j)}, \ k = 1, 2, ..., l; j = 0, 1, 2, ...$.

Next, we differentiate both sides of (5.1) concerning $\tau_1$ to verify the transversality condition.

Taking the derivative of $\lambda$ with respect to $\tau_1$ in (5.1) and substituting $\lambda = i\varphi _{0}$, we obtain

where

If $M_1M_3+M_2M_4 > 0$ holds, then $Re\left (\frac{\mathrm{d} \lambda }{\mathrm{d} \tau_1}\right)_{\lambda = i\varphi _{0}}^{-1} > 0$. Therefore, the transversality condition is satisfied and a Hopf-bifurcation occurs at $E^*$ for $\tau_1 = \hat{\tau}_{10}$.

Proof of Theorem 5.8:

Proof. We consider that $\tau_1 = \tau_1^0$ in its stable interval, with $\tau_2$ as the parameter.

The calculation process is similar to Theorem 5.7.