Bifurcation analysis in a modified Leslie-Gower predator-prey model with fear effect and multiple delays

1.   Introduction

Predation is a fundamental interaction between two species, wherein one species benefits by acquiring resources at the expense of the other. Over the years, numerous mathematical models have been developed, drawing from various biological contexts, to examine the dynamics of predator-prey interactions ^[1,2,3]. Investigating the dynamical behavior of these models enhances our comprehension of the regulatory mechanisms underlying predation, thereby enabling more precise predictions of predator and prey population dynamics.

In general, the dynamic behavior of predator-prey models can be influenced by various factors, including the functional response function, birth rate function, mortality function, and others ^[4,5]. Recent studies have highlighted that certain prey species may exhibit fearful behavior in response to predators, resulting in decreased birth rates and reduced outdoor foraging frequency ^[6,7,8]. Zanette et al. ^[9] observed a forty percent reduction in the reproduction rate of song sparrows due to fear induced by predators. Similarly, Elliott et al. ^[10] demonstrated that the presence of mantis odor in the environment increased the likelihood of a group of drosophilas going extinct sevenfold. Furthermore, some scholars have shown that predator-induced fear can lead to a decrease in prey populations, with this effect often surpassing that of direct predation ^[6,11,12,13]. Taking this fear effect into consideration, Wang et al. ^[14] introduced a fear factor function $\frac{1}{1+ky}$ into their model. Their results suggest that high levels of fear (or strong anti-predator behavioral responses) can stabilize the system by preventing the presence of periodic solutions. They also demonstrated that low levels of fear effects may lead to a limit cycle through Hopf bifurcation. Subsequently, the fear factor function has been widely incorporated into other predator-prey models ^[15,16]. However, Wang and Zou considered only the disadvantages of prey's anti-predation responses, overlooking the potential benefits. In reality, prey's anti-predation responses obviously decrease the likelihood of being captured by predators, as reflected in the model from ^[17].

Leslie and Gower introduced the classic Leslie-Gower model in ^[18], where they assumed that the carrying capacity of the predator is directly proportional to the number of its preferred prey. This classical model sets an upper limit on the relative growth rate of predators, which holds more biological significance. Consequently, numerous researchers have conducted extensive studies based on this model, leading to several important theoretical findings ^[19,20,21]. While this assumption holds true under certain conditions, the predation dynamics depicted are not universally applicable. Predators often need to diversify their prey sources when their preferred prey is scarce ^[22,23]. Mathematically, this situation is typically represented by augmenting the denominator of the Leslie-Gower term with a constant $a$, indicating the predator's reliance on alternative prey in the absence of its favored prey. We introduce the fear effect of prey on predators into a modified Leslie-Gower predator-prey model, resulting in the following two-dimensional model:

In model (1.1), $x$ and $y$ represent the densities of prey and predator species, respectively. Parameters $r$ and $s$ are the intrinsic growth rates of prey and predators, while $d_1$ denotes the natural (density-independent) death rate, and $d_2$ represents the density-dependent death rate due to intraspecific competition. The prey's fear-induced anti-predation response is characterized by $\frac{1}{1+ky(t)}$, where parameter $k$ denotes the perceived level of fear-based response by the prey ^[14]. Conversely, the prey's active anti-predation response, aimed at reducing the chances of being captured, is described by $\frac{1}{1+ck}$, where $c$ is the extent of predation reduction relative to the response level $k$ ^[17]. Additionally, $m$ denotes the capture rate, and $q$ represents the handling time for each prey captured. The term $\frac{y}{nx}$ is referred to as the Leslie-Gower term, while $\frac{y}{nx+a}$ represents a modified version of the growth law of the predator population in the Leslie-Gower model. Here, $nx+a$ denotes the carrying capacity of predators, where $n$ measures the quality of prey as food for the predator, and $a$ represents the availability of alternative food sources for the predator. Compared to the classical Leslie-Gower model, the modified Leslie-Gower model holds greater practical significance.

Furthermore, to enhance the realism of the model, constructing time-delayed differential equation models is a well-established technique used to describe delayed activity at specific stages. The inclusion of time delays often alters the dynamic behaviour of the models. In recent years, an increasing body of literature has incorporated relevant delays into predator-prey models and interacting population models, such as those involving incubation, gestation, maturation, and digestive biological processes ^{[24,25,26,27]}. Typically, researchers introduce discrete time delays into existing ordinary differential equation models to construct corresponding delayed models and explore their dynamics ^[28,29,30]. The introduction of time delays may destabilise the stability of the positive equilibrium and induce periodic oscillations. These oscillations, induced by delay, can be periodic, quasi-periodic, or even chaotic, representing a characteristic feature of time-delayed predator-prey models. Research on predator-prey models with time delays has mostly focused on the instability induced by time delays, the derivation of thresholds, and the stability of oscillatory solutions. In reality, after perceiving certain cues, prey require time to assess the risk. Hence, fear of predation risk does not respond immediately to changes in prey populations but involves a certain time lag. Consequently, we introduce a delay $\tau$, termed the fear response delay ^[31,32]. On the other hand, there is a time lag between predators consuming their prey and reproducing offspring. Therefore, it is reasonable to introduce a delay $\sigma$, known as the gestation delay, into model (1.1). This delay represents the time lag between consuming prey and reproducing offspring. It is assumed that the rate of change of predator species is determined by the number of prey and predators existing at time $(t-\sigma)$, and that the delay $\sigma$ is only considered in the numerical response ^[32,33,34]. Considering these two time delays, we obtain the following model:

where $b = \frac{m}{1+ck}$, $\tau$ represents the fear response delay, while $\sigma$ is the time delay caused by the predator's gestation. All of the parameters are positive. Our aim of this study is to analyze the impact of the two time delays on the dynamic behavior of the model. Moreover, model (1.2) satisfies the non-negative conditions $x(\theta) = \phi_1(\theta)\geq0, y(\theta) = \phi_2(\theta)\geq0, \theta\in[-\mu, 0], \mu = \max\{\tau, \sigma\}, \phi_i(\theta)\in C([-\mu, 0]\rightarrow R_+)$, $i = 1, 2.$

The remainder of this paper is structured as follows. In the following section, we first verify the positivity and boundedness of the solutions to model (1.2). In Section 3, we investigate the existence of Hopf bifurcation at the positive equilibrium, which is caused by the two time delays, as well as the local stability of every possible equilibrium. In Section 4, we calculate the normal form to derive the properties of Hopf bifurcation. Numerical simulations are performed to substantiate our analytical results in Section 5. A summary as well as an outlook for future research of the work are given in conclusion.

2.   Preliminaries

The positivity and boundedness of the model solutions are prerequisites for the study. Positivity signifies species survival and boundedness is considered to be a result of limited resources.

Lemma 2.1. The solution $(x(t), y(t))$ of model (1.2) with initial conditions $(x(0), y(0))\in R_+^2$ is always positive.

Proof. Model (1.2) can be rewritten as

Integrating the above differential equation in the interval $[0, t]$, for $x(0) > 0$, we have

Then $x(t)$ is positive for all $t$ and $x(0) > 0$. Same process as above, we get

Similarly, $y(t)$ is positive for all $t$ and $y(0) > 0$.

Therefore, this lemma is provided a thorough proof.

Lemma 2.2. Given any non-negative initial value, the solution of model (1.2) is bounded. Moreover, the set $\Omega = \{(x(t), y(t))\in R_+^2\ |\ 0\leq x(t)\leq K, \ 0\leq y(t)\leq L\}$ is a positive invariant, where $K = \frac{r}{d_2}$, $L = (nK+a)e^{s\sigma}$.

Proof. From the first equation of model (1.2), we get

Hence, we further have

From the second equation of model (1.2) we obtain

which implies that

Furthermore, for any $\epsilon > 1$, there exists a positive $T_\epsilon$ such that $x(t) < \epsilon K$ for $t > T_\epsilon$. Thus, for $t > T_\epsilon+\sigma$, one can obtain that

which means $\limsup\limits_{t\rightarrow +\infty} y(t)\leq L_\epsilon,$ where $L_\epsilon = (\epsilon nK+a)e^{s\sigma}$. This lemma is supported by letting $\epsilon\rightarrow 1$.

3.   Stability analysis and Hopf bifurcation

In this section, we will analyze the existence and local stability of each possible equilibrium of model (1.2), and discuss the existence of the Hopf bifurcation induced by the delays $\tau$ and $\sigma$.

3.1.   Existence of equilibria

Obviously, model (1.2) always has the trivial equilibrium $E_0 = (0, 0)$ and the predator-only equilibrium $E_1 = (0, a)$, and if $r > d_1$, the model also has a prey-only equilibrium $E_2 = \left(\frac{r-d_1}{d_2}, 0\right)$. Any positive equilibrium $E^* = (x^*, y^*)$, if exists, must satisfy

from which we know that $y^* = nx^*+a$ and $x^*$ is a positive root of the following cubic algebraic equation:

where $a^* = d_2knq > 0, b^* = ad_2kq+bkn^2+d_1knq+d_2kn+d_2q > 0, c^* = 2abkn+ad_1kq+ad_2k+d_1kn+bn+d_1q-qr+d_2, d^* = a^2bk+ad_1k+ab+d_1-r.$

Obviously, if $F(0) = d^* < 0$, Eq (3.1) has one unique positive root, and therefore model (1.2) has one unique positive equilibrium. Denote

If $F(0) = d^* > 0$, according to Shengjin formula, regarding the quantity and presence of positive equilibrium, we can arrive to the following conclusions:

(i) If $c^* < 0$,

a) There exist two positive equilibria when $\Delta < 0$ (see Figure 1(a));

b) There exists a double positive equilibrium when $\Delta = 0$ (see Figure 1(b));

c) There exists no positive equilibrium when $\Delta > 0$ (see Figure 1(c)).

(ii) If $c^*\geqslant 0$, there exists no positive equilibrium (see Figure 1(d)).

3.2.   Stability and Hopf bifurcation

Theorem 3.1. When they exist, the equilibria $E_0 = (0, 0)$ and $E_2 = \left(\frac{r-d_1}{d_2}, 0\right)$ of model (1.2) are always unstable.

Proof. The characteristic equations of model (1.2) at $E_0$ and $E_2$ are

and

respectively. It is clear that both of them have a positive root $\lambda = s$, indicating that $E_0$ and $E_2$ are always unstable.

Remark 1. Theorem 3.1 demonstrates that the boundary equilibria $E_0 = (0, 0)$ and $E_2 = \left(\frac{r-d_1}{d_2}, 0\right)$ are always unstable, independent of the two time delays.

Theorem 3.2. When $r < (ka+1)(ab+d_1)$ and $0\leq\sigma < \sigma_0$, the predator-only equilibrium $E_1 = (0, a)$ of model (1.2) is locally asymptotically stable; Otherwise it is unstable.

Proof. The characteristic equation of model (1.2) at $E_1$ is

Obviously, $\lambda_1 = \frac{r}{ka+1}-d_1-ab$ is one root of Eq (3.2) and other characteristic roots satisfy

If $\sigma = 0$, then the root of (3.3) is $\lambda_2 = -s < 0$. When $\sigma > 0$, assuming $i\omega_0 \; (\omega_0 > 0)$ is a root of Eq (3.3), and plugging it into the equation, we obtain

Further calculations yield $\omega_0 = s$ and

Moreover, we can easily verify the transversality condition:

Therefore, $E_1$ is locally asymptotically stable if $r < (ka+1)(ab+d_1)$ and $0\leq\sigma < \sigma_0$, and unstable if $r > (ka+1)(ab+d_1)$ or $\sigma\geq\sigma_0$. Thus this theorem is proved.

Remark 2. Theorem 3.2 indicates that the predator gestation delay $\sigma$ affects the stability of equilibrium $E_1$. That is to say, when $r < (ka+1)(ab+d_1)$, the prey will eventually die out, and the predator can rely on additional prey to maintain its growth provided $\sigma < \sigma_0$.

Next, we discuss the local stability of the positive equilibrium. We make the following assumption:

$(\textbf{H1})$ (i) $d^* < 0$ or (ii) $c^* < 0$ and $\Delta < 0$.

Based on the analysis in Subsection 3.1, we know that model (1.2) has at least one positive equilibrium $E^* = (x^*, y^*)$ if (H1) holds. By Linearizing the model at $E^*$, we can obtain that:

where $I_2$ is a two-dimensional unit matrix and

where

Therefore, the characteristic equation at $E^*$ is

where $p_1 = -l_{11}-l_{22}, p_0 = l_{11}l_{22}, n_1 = -m_{22}, n_0 = l_{11}m_{22}-l_{12}m_{21}, q_0 = -m_{12}m_{21}.$ To analyze the separate and combined effects of two time delays on the stability of the positive equilibrium $E^*$, we discuss them in the following five subsections.

3.2.1.   $\tau=0, \sigma=0$

In this assumption, Eq (3.6) reduces to

Therefore, if

${(\textbf{H2})} \; p_1+n_1 > 0, n_0+q_0 > 0$

holds, all the roots of Eq (3.7) have negative real parts. Thus, we have the following conclusion.

Theorem 3.3. In the absence of time delays, if both (H1) and (H2) hold, then the positive equilibrium $E^*$ is locally asymptotically stable.

In the following, to explore the effect of time delays on the stability of the positive equilibrium, we always assume that the positive equilibrium is present and locally stable in the absence of delays, i.e., both (H1) and (H2) are satisfied.

3.2.2.   $\tau > 0, \sigma=0$

In this subsection, we ignore the gestation delay of the predator ($\sigma = 0$) and only consider the effect of the fear delay ($\tau > 0$) on the stability of $E^*$. The corresponding characteristic Eq (3.6) becomes

Applying $\lambda = i\omega\ (\omega > 0)$ into Eq (3.8) and dividing the imaginary and real components:

which leads to

where $\nu = \omega^2$. Notice that $n_0^2-q_0^2$ has the same sign as $n_0-q_0$ since $n_0+q_0 > 0$. Next, we discuss the distribution of the roots of the Eq (3.10) in three situations.

(I) If condition

${(\textbf{H3})} \; n_0 > q_0, (p_1+n_1)^2 > 2\left(n_0-\sqrt{n_0^2-q_0^2}\right)$

holds, then Eq (3.10) has no positive root, implying Eq (3.8) has no pure virtual root. Consequently, the positive equilibrium $E^*$ is locally asymptotically stable for any $\tau > 0$ under conditions (H2) and (H3).

(II) If condition

${(\textbf{H4})} \; n_0 > q_0, (p_1+n_1)^2 < 2\left(n_0-\sqrt{n_0^2-q_0^2}\right)$

holds, then Eq (3.10) has two different positive roots

Accordingly, the critical values of delay $\tau$ are

Since $\omega_1^+ > \omega_1^-$, then $\tau_1^{(j+1)+}-\tau_1^{j+} = \frac{2\pi}{\omega_1^+} < \frac{2\pi}{\omega_1^-} = \tau_1^{(j+1)-}-\tau_1^{j-}$, for $j = 0, 1, 2, \cdots$ ^[35]. From Eq (3.8), we obtain

Now, we establish the stability switch due to delay as follows:

If $\tau_1^{0+} < \tau_1^{1+} < \tau_1^{0-}$. All the roots of Eq (3.8) have negative real parts when $0\leq\tau < \tau_1^{0+}$, and there is a pair of pure imaginary roots at $\tau = \tau_1^{0+}$. As $\tau$ increases and approaches $\tau_1^{0-}$, the above pure imaginary roots cross the imaginary axis and are accompanied by positive real parts so that $E^*$ becomes unstable. Similarly, when $\tau > \tau_1^{1+}$, there exist two pairs of eigenvalues with positive real parts and $E^*$ is unstable. Since there is only one pair of eigenvalues crossing the imaginary axis at each $\tau = \tau_1^{j+}$ or $\tau = \tau_1^{j-}$, and two consecutive $\tau_1^{j-}$ are not possible, instability always exists when $\tau > \tau_1^{0+}$.

If $0 < \tau_1^{0+} < \tau_1^{0-} < \tau_1^{1+} < \tau_1^{1-} < \cdots < \tau_1^{n+} < \tau_1^{(n+1)+} < \tau_1^{n-} < \cdots$ for some positive integer $n$. Then when $\tau\in(0, \tau_1^{0+})\cup(\tau_1^{0-}, \tau_1^{1+})\cup\cdots\cup(\tau_1^{(n-1)-}, \tau_1^{n+})$, all the roots of Eq (3.8) have negative real parts; when $\tau\in(\tau_1^{0+}, \tau_1^{0-})\cup(\tau_1^{1+}, \tau_1^{1-})\cup\cdots\cup(\tau_1^{(n-1)+}, \tau_1^{(n-1)-})$, Eq (3.8) has at least a pair of conjugate complex roots with positive real parts. Once the delay satisfy $\tau_1^{n+} < \tau_1^{(n+1)+} < \tau_1^{n-}$, similar to the analysis above, $E^*$ becomes unstable for all $\tau > \tau_1^{n+}$ ^[35,36,37].

(III) If condition

${(\textbf{H5})} \; n_0 < q_0$

holds, then Eq (3.10) has a unique positive root

further calculations we have

Let $\tau_{10} = \min{\tau_1^{(0)}}$, then

and

The following theorem can be derived from the analysis above.

Theorem 3.4. Assume that (H1) and (H2) hold. Then for $\sigma = 0$ we have

(I) If (H3) holds, $E^*$ is always locally asymptotically stable for any $\tau\geq0$.

(II) If (H4) holds, for $\tau_1^{0+} < \tau_1^{1+} < \tau_1^{0-}$, $E^*$ is locally asymptotically stable for $0\leq\tau < \tau_1^{0+}$ and unstable for $\tau > \tau_1^{0+}$, model (1.2) undergoes a Hopf bifurcation at $E^*$ when $\tau = \tau_1^{0+}$. For $0 < \tau_1^{0+} < \tau_1^{0-} < \tau_1^{1+} < \tau_1^{1-} < \cdots < \tau_1^{n+} < \tau_1^{(n+1)+} < \tau_1^{n-} < \cdots$, $E^*$ is locally asymptotically stable for $\tau\in(0, \tau_1^{0+})\cup(\tau_1^{0-}, \tau_1^{1+})\cup\cdots\cup(\tau_1^{(n-1)-}, \tau_1^{n+})$ and unstable for $\tau\in(\tau_1^{0+}, \tau_1^{0-})\cup(\tau_1^{1+}, \tau_1^{1-})\cup\cdots\cup(\tau_1^{(n-1)+}, \tau_1^{(n-1)-})\cup(\tau_1^{n+}, +\infty)$. In addition, model (1.2) undergoes a Hopf bifurcation at $E^*$ when $\tau = \tau_1^{j\pm}$.

(III) If (H5) holds, $E^*$ is locally asymptotically stable for $0\leq\tau < \tau_{10}$ and unstable for $\tau > \tau_{10}$. Moreover, model (1.2) undergoes a Hopf bifurcation at $E^*$ when

3.2.3.   $\tau=0, \sigma > 0$

In this subsection, we only consider the effect of the gestation delay of the predator $\sigma$ on the stability of $E^*$, characteristic Eq (3.6) becomes

We substitute $\lambda = i \omega$ into (3.12) and then obtain

Separating the real and imaginary parts leads to

Squaring and adding both above equations lead to

Making $\omega^2 = v$, the above equation can be re-written by

Obviously, $(n_0+q_0)^2 > 0$ when (H2) holds. Therefore (3.13) has a unique positive root

By calculation, we obtain

Define the first critical value as $\sigma_{10} = \min\limits{\sigma_{1}^{(j)}}$. Moreover, we can count

Then

based on the fact $e^{-\lambda \sigma} = -\left(\frac{\lambda^2+p_1\lambda}{n_1\lambda+n_0+q_0}\right)$. Thus

We arrived at the following conclusions according to the aforementioned analysis.

Theorem 3.5. Assume that (H1) and (H2) hold. Then for $\tau = 0$, $E^*$ is locally asymptotically stable for $0\leq\sigma < \sigma_{10}$ and is unstable for $\sigma > \sigma_{10}$. Moreover, the model undergoes a Hopf bifurcation at $E^*$ when $\sigma = \sigma_{10}$, where

Remark 3. Theorem 3.5 shows that when the gestation delay $\sigma$ crosses a threshold value, the stability of the positive equilibrium of model (1.2) changes, accompanied by the occurrence of a Hopf bifurcation.

3.2.4.   $\sigma > 0, \tau\in (0, \tau_{10})$

This subsection explores the impact of delay $\sigma$ on the stability of equilibrium $E^*$, fixing $\tau$ in the interval $(0, \tau_{10})$. The corresponding characteristic equation at $E^*$ is

Let $\lambda = i\omega\ (\omega > 0)$ be a root of (3.16), then we have

which yields

After squaring both sides of (3.18) and adding them together, we have

Obviously, there are $H(0) = -(n_0+q_0)^2 < 0$ and $H(+\infty) = +\infty$. Eq (3.19) is a transcendental equation, and the prediction of the roots is a difficult task. Now, it is considered that the equation has at least a finite number of positive roots $\omega_{21}, \omega_{22}, \cdots, \omega_{2k}$. After a complex calculation process, we get

for $k = 1, 2, 3, 4; j = 0, 1, 2, \cdots,$ and $\pm\omega_{2k}$ is a pair of pure imaginary roots of (3.16).

Next we aim to verify the transversality condition. By plugging $\lambda(\sigma) = \alpha(\sigma)+i\omega(\sigma)$ into (3.16) and separating the real and imaginary part, we get

Let $\sigma_{20} = \min_{k = 1, 2, 3, 4} \{{\sigma_{2k}^{(0)}}\}, \omega_{20} = \omega_{2k}^{(0)}$. Further, we differentiate (3.16) with respect to $\sigma$ and substitute $\sigma = \sigma_{20}$, then

where

For Eq (3.21), we have

If $AC-BD\neq0$, then the transversal condition is satisfied.

Therefore, we have the following conclusion.

Theorem 3.6. Suppose that the conditions in Theorem 3.3 hold. For a given $\tau \in (0, \tau_{10})$, $E^*$ is locally asymptotically stable when $0\leq\sigma < \sigma_{20}$ and unstable when $\sigma > \sigma_{20}$. Moreover, model (1.2) undergoes a Hopf bifurcation at $E^*$ when $\sigma = \sigma_{20}$, where

Remark 4. Theorem 3.6 shows that by fixing the fear response delay to a stable interval and taking the predator's gestation delay as a bifurcation parameter, it is discovered that the stability of the positive equilibrium is altered by the time delay $\sigma$ at the critical value. The above conclusion implies that both time delays affect the kinetic properties of the population.

3.2.5.   $\tau > 0, \sigma\in (0, \sigma_{10})$

This situation is similar to $\sigma > 0, \tau\in (0, \tau_{10})$, and we can easily get the following results using the same method as in the last subsection.

Theorem 3.7. Suppose that the conditions in Theorem 3.3 hold. For a given $\sigma \in (0, \sigma_{10})$, $E^*$ is locally asymptotically stable when $0\leq\tau < \tau_{20}$ and unstable when $\tau > \tau_{20}$. Moreover, model (1.2) undergoes a Hopf bifurcation at $E^*$ when $\tau = \tau_{20}$, where

The proof of Theorem 3.6 is comparable to this one. Here, we leave it out.

4.   Direction of the Hopf bifurcation

In previous section, we proved that the model (1.2) generates a Hopf bifurcation at $E^*$ when $\sigma = \sigma_{20}, \tau\in[0, \tau_{10})$. Here, we apply the normal form theorem and the center manifold theorem by Hassard et al. ^[38] to discuss the direction of Hopf bifurcation and the stability of the periodic solutions.

We need to only calculate the coefficients $\mu_2, \beta_2$ and $T_2$. Since the derivation is standard and lengthy, we place them in the Appendix. Consequently, the coefficients $g_{ij}$ can be solved, and we can obtain the values:

which decides the properties of Hopf bifurcations at the critical value $\sigma = \sigma_{20}$.

Theorem 4.1. For model (1.2), the Hopf bifurcation is supercritical (resp., subcritical) if $\mu_2 > 0$ (resp., $\mu_2 < 0$). The bifurcating periodic solutions are stable (resp., unstable) if $\beta_2 < 0$ (resp., $\beta_2 > 0$). The period increases (resp., decreases) if $T_2 > 0$ (resp., $T_2 < 0$).

5.   Numerical simulations

In this section, we conduct numerical simulations to demonstrate the theoretical findings presented in the preceding sections. Initially, we simulate the impact of time delays on the stability of $E^*$ in model (1.2) under various scenarios. Additionally, we supplement our analysis by simulating the effects of fear intensity and the presence of additional prey on $E^*$.

5.1.   The effect of two time delays on Hopf bifurcation

We fix the other parameters as follows:

Then $\frac{r}{ka+1}-d_1-ab = 1.344615385 > 0$, and model (1.2) has three unstable boundary equilibria $E_0 = (0, 0), E_1 = (0, 0.1)$ and $E_2 = (89.5, 0)$ by Theorems 3.1 and 3.2. Moreover, conditions (H1) and (H2) hold, model (1.2) exists a unique positive equilibrium $E^* = (4.0856, 2.1428)$, which is locally asymptotically stable for $\tau = \sigma = 0$ (see Figure 2).

Under the parameter setting in (5.1), $n_0^2-q_0^2 = -0.004025 < 0$, so condition (H5) holds. In the case of $\sigma = 0$, we can calculate the critical value $\tau_{10} = 6.823$ from (3.11); the corresponding bifurcation diagram is drawn in Figure 3(a). From Figure 4, we can observe that $E^*$ is locally asymptotically stable for $\tau < \tau_{10}$ (see Figure 4(a)); when the fear response delay $\tau$ crosses the critical value $\tau_{10}$, $E^*$ loses original stability and a stable periodic solution develops (see Figure 4(b)). When $\tau = 0$, the critical value of Hopf bifurcation of gestation delay $\sigma$ is $\sigma_{10} = 1.9715$ from (3.15), the corresponding bifurcation diagram is drawn in Figure 3(b). From Figure 5, we can observe $E^*$ is locally asymptotically stable for $\sigma < \sigma_{10}$ (see Figure 5(a)) and becomes unstable when $\sigma > \sigma_{10}$ (see Figure 5(b)), accompanied by a stable periodic solution. For the stability switching phenomenon presented by Theorem 3.4, due to the different parameter values taken, we stay in the following subsection to give the simulation results.

Now, we simulate the cases in which both time delays exist. Taking $\tau = 2.82\in[0, \tau_{10})$, by Eq (3.22), we easily get the Hopf bifurcation point $\sigma_{20} = 1.615$, and the time series diagrams of the two species can be seen in Figure 6. We can find from the figures that $E^*$ is locally asymptotically stable for $\sigma < \sigma_{20}$ and unstable for $\sigma > \sigma_{20}$. Moreover, the biomass of both species showed periodic fluctuations over time. Further, using the algorithms in Section 4, we obtain $c_1(0) = -0.0086-0.0323i, \ \beta_2 = -0.0172 < 0, \ \mu_2 = 0.1660 > 0, \ T_2 = 0.0569 > 0$. Fixing $\sigma = 1.246\in(0, \sigma_{10})$, we can obtain the Hopf bifurcation point $\tau_{20} = 4.291$ from (3.23). The time series diagrams of the two species can be seen in Figure 7. We can find from the figures that the equilibrium $E^*$ is locally asymptotically stable for $\tau < \tau_{20}$ and unstable for $\tau > \tau_{20}$. Similarly, further we calculate $c_1(0) = -0.0134 - 0.0114i, \ \beta_2 = -0.0268 < 0, \ \mu_2 = 1.5812 > 0, \ T_2 = 0.1514 > 0$.

5.2.   Stability switch of $E^*$ induced by fear response delay $\tau$

In this subsection, regardless of the gestation delay $\sigma = 0$, we specifically simulate the stability switching phenomenon induced by fear response delay $\tau$ at equilibrium $E^*$. For this purpose, we take the parameters: $r = 1.8, k = 2, d_1 = 0.01, d_2 = 0.01, b = 0.9, q = 0.2, s = 0.3, n = 0.8, a = 0.1.$

In this case, all the conditions (H1), (H2), and (H4) hold, and then equilibrium $E^* = (0.961, 0.851)$ is locally asymptotically stable for $\tau = 0$. The corresponding critical values can be calculated as follows:

The bifurcation diagrams of model (1.2) with respect to $\tau$ are shown in Figure 8. We can observe that equilibrium $E^*$ is locally asymptotically stable for $\tau\in(0, \tau_1^{0+})\cup(\tau_1^{0-}, \tau_1^{1+})$ and unstable for $\tau\in(\tau_1^{0+}, \tau_1^{0-})\cup(\tau_1^{1+}, +\infty)$, and model (1.2) undergoes a Hopf bifurcation at $E^*$ when $\tau = \tau_1^{j\pm}$. The corresponding time series diagrams of two species can be found in Figure 9, which confirms the occurrence of the stability switch phenomenon as delay $\tau$ increases.

5.3.   The effect of fear intensity $k$ and alternative prey $a$

The simulations in the previous subsections show that the change of fear response delay $\tau$ and gestation delay $\sigma$ have a a tremendous effect on the stability of $E^*$, and the equilibrium presents a stability switching phenomenon in special cases. In this subsection, we demonstrate the effect of fear intensity and alternative prey on model $(1.2)$ using numerical methods.

First, we simulate the effect of fear intensity $k$ on the critical values $\tau_{10}$ and $\sigma_{10}$. Taking $r = 1.8$, $d_1 = 0.01$, $d_2 = 0.02$, $b = 0.3$, $q = 0.8$, $s = 0.38$, $n = 0.5$, $a = 0.1$, we draw the relationship diagram of $\tau_{10}$ and $\sigma_{10}$ and fear intensity $k$ in Figure 10. We can observe that both the critical values depend on fear intensity $k$. Further observation reveals that both critical values increase when the fear intensity $k$ is larger. Biologically, the increased fear effect makes it more difficult for predators to hunt their primary prey, allowing the two species to co-exist in a new state of affairs. Periodic oscillations of the two species are more difficult to occur. That is, the fear intensity strengthens the stability of the positive equilibrium between the two species.

Moreover, we will explore the effect of fear intensity on population density. The relationship between the fear intensity $k$ and the biomass of the two populations is shown in Figure 11. It can be seen that the fear intensity $k$ reduces the biomass of both species overall. The reason for this is that as fear intensity increases, scared prey may adopt more anti-predator behaviors, which will reduce prey reproduction rates, and thus prey densities will subsequently decrease. Therefore, as the biomass of prey decreases, so does the biomass of predators that feed on them.

Next, we simulate the effect of alternative prey $a$ on the critical values $\tau_{10}$ and $\sigma_{10}$. The parameters are taken as $r = 1.8, k = 3, d_1 = 0.01, d_2 = 0.02, b = 0.3, q = 0.8, s = 0.38, n = 0.5$. As seen in Figure 12, both time delay thresholds increased with the increase of alternative food $a$. This suggests that the alternative food strengthens the stability of the positive equilibrium between the two species. We also notice that the biomass of both populations decreases as alternative food $a$ increases (see Figure 13). This is because of the rapid increase in predators that results from the rise in alternative food. The major prey $x$ is reduced by substantial predation. Therefore, the carrying capacity $nx+a$ of predators is reduced. Predators $y$ decreases due to the lack of a suitable environment for survival.

6.   Conclusions

We have introduced and examined a modified Leslie-Gower predator-prey model incorporating a fear effect. The model accounts for the fear response delay of prey to predators and the gestation delay of predators. We have analyzed the positivity and boundedness of the model solutions, which are essential prerequisites for ensuring the biological relevance of the model. Subsequently, we conducted a local stability analysis for each feasible equilibrium and discussed the Hopf bifurcations induced by the dual delays near the positive equilibrium. Our findings reveal that the trivial equilibrium $E_0 = (0, 0)$ and the prey-only equilibrium $E_2 = \left(\frac{r-d_1}{d_2}, 0\right)$ (if $r > d_1$) are consistently unstable, indicating the persistence of prey populations for any positive initial values. However, the stability of the predator-only equilibrium $E_1 = (0, a)$ is influenced by the gestation delay $\sigma$ of predators.

For the concurrent presence of two time delays, we examined their impact on the stability of equilibrium $E^*$ in four distinct scenarios. When the model incorporates only one time delay, we identified the critical value ($\tau_{10}$ or $\sigma_{10}$) of the delay that destabilizes equilibrium $E^*$. In cases where both time delays occur, with one time delay fixed ($\sigma\in(0, \sigma_{10}$) or $\tau\in(0, \tau_{10}$)), and the other allowed to vary, we also determined the critical value ($\tau_{20}$ or $\sigma_{20}$) of the delay that destabilises equilibrium $E^*$. Upon crossing certain critical values, $E^*$ loses its stability, and a set of periodic solutions emerges from $E^*$. This underscores the significant impact of the two delays on the dynamics of model (1.2). We observed that $\tau_{20} < \tau_{10}$ and $\sigma_{20} < \sigma_{10}$ under the same other conditions, suggesting that the combined effects of the two time delays are more likely to destabilize equilibrium $E^*$ than a single delay (see Figures 6 and 7). Consequently, both delays exert a considerable influence on the dynamics of model (1.2). Utilising the normal form theory and centre manifold theorem, we further computed the values of $\beta_2$, $\mu_2$, and $T_2$. Our results indicate that the Hopf bifurcation is supercritical, and the resulting periodic solutions are stable, with their period increasing. Towards the end of the article, we also examined the effects of fear intensity $k$ and the availability of alternative food $a$ for predators on the biomass and stabilisation intervals of both populations. We found that both $k$ and $a$ enhance equilibrium stability while reducing the biomass of both species. Therefore, to maintain ecological balance within populations, effective measures can be taken by humans to control prey abundance, adjust fear levels, and manipulate the gestation period of the fear response time delay, thereby ensuring the sustained coexistence of both species within the ecosystem.

In biological terms, these findings suggest that relatively short time delays lead to stable levels in both predator and prey populations. However, as the time delays increase, species abundance fluctuates cyclically, indicating that predators and prey coexist in a cyclically oscillating manner. This study enhances our understanding of predation dynamics in nature and facilitates the investigation of population growth, with potential applications in areas such as biodiversity conservation and economic management.

While we focused on the effect of time delays on the dynamics of model (1.2), it is worth considering whether spatial movement influences system dynamics. From a modeling perspective, incorporating more parameters of practical significance could yield more insightful results. On the analytical front, exploring additional dynamical properties, such as consistent boundedness, persistent behaviour, and other types of bifurcation, could provide further insights. Moving forward, we aim to refine this model and investigate its dynamical properties to obtain more robust theoretical results.

Acknowledgments

We are very grateful to Professor Tonghua Zhang for his valuable comments and suggestions, which have greatly improved the quality and presentation of our paper. This work was supported by the Natural Science Foundation of Shanghai (23ZR1445100), the National Natural Science Foundation of China (11671260, 12071293), and the Science and Technology Key Project of Henan Province of China (242102110341).

Conflict of interest

Sanling Yuan is an editorial board member for Mathematical Biosciences and Engineering and was not involved in the editorial review or the decision to publish this article. All authors declare that there are no competing interests.

Appendix

A1.   Computation of the coefficients $\mu_2$, $\beta_2$, and $T_2$

Throughout this section, we compute the coefficients $\mu_2, \beta_2$, and $T_2$ to determine the properties of Hopf bifurcation.

Let $u_1 = x-x^*, u_2 = y-y^*,$ and $\sigma = \sigma_{20}+\mu$, then model (1.2) can be reduced to the following functional differential equation in $C = C\left([-1, 0], R^2\right)$:

where $u(t) = \left(u_1(t), u_2(t)\right)^T \in R^2,$ and $L_\mu: C \rightarrow R^2, F: C \times R \rightarrow R^2$.

and

where

and

Hence, by Riesz representation theorem, there exists a $2 \times 2$ matrix function $\eta(\theta, \mu)$ of bounded variation when $\theta \in[-1, 0]$, which implies

As a matter of fact, we can choose

For $\phi \in C\left([-1, 0], R^2\right)$, we can define

Then system (A1) can be taken the form $\dot{u_t} = A(\mu) u_t+R (\mu) u_t,$ where $u_t = u(t+\theta), \theta \in[-1, 0]$.

For $\psi \in C\left([-1, 0], \left(R^2\right)^*\right)$, where $\left(R^2\right)^*$ is the 2-dimensional space of row vectors, the adjoint operator $A^*$ of $A(0)$ is further defined as

For $\phi \in C\left([-1, 0], R^2\right)$ and $\psi \in C\left([-1, 0], \left(R^2\right)^*\right)$, we define the bilinear form

where $\eta(\theta) = \eta(\theta, 0), A = A(0)$ and $A^*$ are adjoint operators.

In accordance with the previous section, we observe that $\pm i \omega_{20} \sigma_{20}$ are eigenvalues of $A(0)$. Hence, they are the eigenvalues of $A^*$. Suppose that $q(\theta) = (1, \delta)^T e^{i \omega_{20} \sigma_{20} \theta}$ is the eigenvector of $A(0)$ corresponding to $i \omega_{20} \sigma_{20}$ and $q^*(s) = D\left(1, \delta^*\right) e^{i \omega_{20} \sigma_{20} s}$ is the eigenvector of $A^*$ corresponding to $-i \omega_{20} \sigma_{20}.$ Then $A(0)q(\theta) = i\omega_{20} \sigma_{20}q(\theta)$. It follows from the definition of $A(0)$ and $\eta(\theta, \mu)$ that

Thus, we can easily get

Similarly, let $A^*(0){q^*}^T(s) = -i\omega_{20} \sigma_{20}q(s)$, we have

Then we can easily obtain $\delta^* = \frac{-\left(i \omega_{20}+l_{11}\right)}{m_{21} e^{i \omega_{20} \sigma_{20}}}\ \ \text{and}\ \ q^*(0) = D\left(1, \delta^*\right).$ In order to assure $\langle q^*(s), q(\theta)\rangle = 1$, we need to determine the value of $D$. From (A3), we have

This implies that $\bar{D} = \left(1+\delta \bar{\delta}^*+\tau m_{12}\delta e^{-i \omega_{20} \tau}+\sigma_{20} e^{-i \omega_{20} \sigma_{20}}\left(m_{21} \bar{\delta}^*+m_{22} \delta \bar{\delta}^*\right)\right)^{-1}.$ Then we can obtain that $D = \left(1+\bar{\delta} \delta^*+\tau m_{12}\bar{\delta} e^{i \omega_{20} \tau}+\sigma_{20} e^{i \omega_{20} \sigma_{20}}\left(m_{21} \delta^*+m_{22} \bar{\delta} \delta^*\right)\right)^{-1}.$ Next, we use the same notation as Hassard et al. ^[38] and the methods in ^[39,40,41] to describe the center manifold $C_0$ at $\mu = 0$. Let $u_t$ be the solution of (A1) when $\mu = 0$.

Define $z(t) = \langle q^*, u_t\rangle$ and $W(t, \theta) = u_t(\theta)-2 \operatorname{Re}\{z(t) q(\theta)\}$ on the center manifold $C_0$. We have

where $z$ and $\bar{z}$ are the local coordinates for center manifold $C_0$ in the directions of $q^*$ and $\bar{q}^*$. Noting that $W$ is real. Thus when $\mu = 0$, we have $\dot{z}(t) = < q^*, \dot{u}_t > = i \omega_{20} \sigma_{20} z(t)+\bar{q}^*(0) F(z, \bar{z}),$ where $F(z, \bar{z}) = F\left(0, u_t\right)$. That is $\dot{z}(t) = i \omega_{20} \sigma_{20} z(t)+g(z, \bar{z}),$ where

Notice that $u_t\left(u_{1 t}(\theta), u_{2 t}(\theta)\right) = W(t, \theta)+z q(\theta)+\overline{z q}(\theta),$ then

Thus, from the definition of $F(\mu, \mu_t)$, we have

However

Similarly, we have

where $E_i = \left(E_i^{(1)}, E_i^{(2)}\right) \in R^2(i = 1, 2)$ is a constant vector. We now need to seek appropriate $E_1$ and $E_2$. It follows from the definition of $A$ that for $\theta = 0$,

where $\eta(\theta) = \eta(0, \theta)$. Then, we have

where

Noticing that

By further calculation, we have $\left(2 i \omega_{20} \sigma_{20} I-\int_{-1}^0 e^{2 i \omega_{20} \sigma_{20} \theta} d \eta(\theta)\right) E_1 = 2 \sigma_{20}\left(P_1, P_2\right)^T,$ which leads to

Similarly, we can get $\left(\int_{-1}^0 d \eta(\theta)\right) E_2 = 2 \sigma_{20}\left(Q_1, Q_2\right)^T.$ That is

Thus we can determine $W_{20}(\theta)$ and $W_{11}(\theta)$ from (A6) and (A7). Furthermore, $g_{ij}$ in (A5) can be computed.