Bifurcations, stability switches and chaos in a diffusive predator-prey model with fear response delay

1.   Introduction

The interaction between predator and prey is one of the most important topics in mathematical biology and theoretical ecology ^[1]. The direct interaction, which is reflected by predation, has been extensively studied ^[2,3,4,5]. Recently, a field experimental study ^[6] provided evidence that show that indirect effect (e.g., fear effect) cannot be ignored. Even though there is no direct killing between predators and prey, the presence of predators cause a reduction in prey population ^[6,7]. To explore the impact that fear can have on population dynamics, Wang et al. ^[8] formulated the following model incorporating the cost of fear

where $u$ and $v$ are the population densities of the prey and predator, respectively. Function $f(k, v(t))$ reflects the cost of anti-predation response due to fear, where $k$ measures the level of fear. $r_{0}$ is the reproduction rate of prey in the absence of predator, $d$ and $m$ represent the natural death rate, $a$ reflects the death rate due to intra-specific competition, the positive constant $c$ is the efficiency in biomass transfer, $g(u(t))$ is the functional response. After their work, other biological phenomena, such as the Allee effect, harvesting, cooperation hunting, prey refuge, group defense and so on, were introduced into predator-prey models with the fear effect ^{[1,9,10,11,12,13,14,15]}.

In fact, besides the cost in the reproduction of prey due to fear, there are also some benefits for an anti-predation response. Wang and Zou ^[16] described such benefits with $g(u(t), k)$ instead of $g(u(t))$, and considered both linear functional response and Holling type II functional response. In the model with the Holling type II functional response, they choose $g(u(t), k) = \frac{1}{1+c_{1}k}\cdot\frac{pu(t)}{1+qu(t)}$, $f(k, v) = \frac{1}{1+c_{2}kv(t)}$ and get the following model

where $c_{1}, c_{2}$ represent the decreasing rate of reproduction and predation, respectively.

To obtain more resources, species tend to migrate from a high population density area to a low population density area. Therefore, spatial diffusion should be considered when we model the interaction between predator and prey. Many diffusive models have been proposed to investigate the influence of the cost of fear on the spatial distribution of species ^{[17,18,19,20]}. Wang and Zou ^[21] proposed a reaction-diffusion-advection predator prey model, in which conditions of spatial pattern formation are obtained.

In reality, there are time delays in almost every process of predator and prey interaction. Many kinds of delays have been incorporated into predator-prey models with the fear effect ^{[16,17,18,22,23,24,25]}. Since the biomass transfer is not instantaneous after the predation of prey, the biomass transfer delay is considered in ^[16]. A generation time delay in prey has been considered in ^[23]. In fact, the reproduction of prey will not respond immediately to fear, but will rather reduce after a time lag. Such a fear response delay has been considered in ^[18,24,25].

Motivated by ^[16] and the previous work, we propose the following model

where $d_{1}, d_{2} > 0$ represent the diffusion coefficients of prey and predator, respectively, and $\tau$ is the fear response delay in prey reproduction.

Although there have been some literature on a predator-prey model with fear effect and delay, no work has been done to explore the joint impact of fear level and fear response delay on the population dynamics. In this paper, we aim to reveal how these two parameters $k$ and $\tau$ jointly affect the dynamics of system (1.3) from the view of bifurcation analysis. The characteristic equation may have no, one or two pairs of purely imaginary roots under different conditions. Correspondingly, the fear response delay may have no effect, both stabilizing and destabilizing effect, or destabilizing effect on the stability of the positive equilibrium. By the discussion of the monotonicity of the critical values of Hopf bifurcation, we get the order of all the Hopf bifurcation values. We prove the existence of stability switches for system (1.3) under suitable condition. Through numerical simulations, we find that the level of fear $k$ has a crucial role in the stability of positive equilibrium and the occurrence of Hopf bifurcation induced by fear response delay. To reveal the complex phenomena induced by $k$ and $\tau$, double Hopf bifurcation analysis is carried out, which can induce complex spatio-temporal dynamics ^[26,27]. It provides a qualitative classification of dynamical behaviors on the $(k, \tau)$ plane, which can help us to explicitly observe the different dynamics corresponding to different values of $k$ and $\tau$, including quasi-periodic oscillations on two or three dimensional torus, and even chaos.

The rest of this paper is organized as follows. In Section 2, we analyze the local stability of equilibria and obtain the condition of Hopf bifurcation in three different cases. In Section 3, we give the condition for double Hopf bifurcation and derive the normal form of the double Hopf bifurcation. In Section 4, numerical simulations are presented to verify our theoretical results. Finally, conclusions and discussions are given in Section 5.

2.   Local stability and Hopf bifurcation

In this section, we discuss the local stability of equilibria. System (1.3) has three possible equilibria. The trivial equilibrium $E_{0}(0, 0)$ always exists. When $r_{0}-d > 0$, system (1.3) has a predator-free equilibrium $E_{u}(\frac{r_{0}-d}{a}, 0)$. Moreover, if

holds, there is a unique positive equilibrium $E^{*}(u^{*}, v^{*})$, where

with

It is well known that the laplacian operator $\Delta$ has eigenvalues $-n^{2}/l^{2}\; (n = 0, 1, 2, \cdots)$.

Now, we study the local stability for each of equilibria. The linearization of system (1.3) at an equilibrium $(\bar{u}, \bar{v})$ can be written as

where

with

The characteristic equation is given by

From (2.4), we get the corresponding characteristic equation at $E_{0}(0, 0)$ as

Thus, we have

Obviously, $\lambda_{2, n} < 0$. If $r_{0} < d$, $\lambda_{1, n} < 0$ for all $n\in\mathbb{N}_{0}$, and if $r_{0} > d$, $\lambda_{1, 0} > 0$. Thus, if $r_{0} < d$, $E_{0}$ is locally asymptotically stable. If $r_{0} > d$, $E_{0}$ is unstable.

If $r_{0} > d$, the predator-free equilibrium $E_{u}(\frac{r_{0}-d}{a}, 0)$ exists. The corresponding characteristic equation is given by

Since $r_{0}-d > 0$, $\lambda_{1, n} = -d_{1}\frac{n^{2}}{l^{2}}-(r_{0}-d) < 0$. Now we consider the sign of $\lambda_{2, n} = -[d_{2}\frac{n^{2}}{l^{2}}+m-\frac{cp(r_{0}-d)}{(1+c_{1}k)[a+q(r_{0}-d)]}]$. If $(\mathrm{H1})$ holds, $\lambda_{2, 0} = -m+\frac{cp(r_{0}-d)}{(1+c_{1}k)[a+q(r_{0}-d)]} > 0$, which indicates that $E_{u}$ is unstable. If $\frac{m(1+c_{1}k)}{pc-mq(1+c_{1}k)} > \frac{r_{0}-d}{a} > 0$ holds, $\lambda_{2, n} < 0$ for all $n\in\mathbb{N}_{0}$, and $E_{u}$ is locally asymptotically stable.

Theorem 1. (i) When $r_{0} < d$, there is only the trivial equilibrium $E_{0}$ for system (1.3), which is locally asymptotically stable; when $r_{0} > d$, $E_{0}$ is unstable and there is a predator-free equilibrium $E_{u}$.

(ii) When $\frac{m(1+c_{1}k)}{pc-mq(1+c_{1}k)} > \frac{r_{0}-d}{a} > 0$, $E_{u}$ is locally asymptotically stable; When $(\mathrm{H1})$ holds, $E_{u}$ is unstable and there is a positive equilibrium $E^{*}$.

For the positive equilibrium $E^{*}(u^{*}, v^{*})$ of system (1.3), we have

The characteristic equation becomes

where

When $\tau = 0$, the characteristic equation (2.6) becomes

Obviously, $D_{0}+M = -J_{21}J_{12}-J_{21}K_{12} > 0$. Hence, if

holds, we have $\ T_{n} > 0$, $D_{n}+M > 0$ for all $n\in\mathbb{N}_{0}$, thus all roots of Eq (2.7) have negative real parts. If $J_{11} > 0$, we have $T_{0} < 0$ and $D_{0}+M > 0$, the roots of (2.7) have positive real parts when $n = 0$, which means that $E^{*}$ is unstable in the absence of diffusion. To sum up, Turing instability induced by diffusion will not occur.

In the following, we assume that $(\mathrm{H2})$ always holds, which ensures $\ T_{n} > 0$, $D_{n} > 0$ and $D_{n}+M > 0$. Thus, a change of stability at $E^{*}(u^{*}, v^{*})$ can only happen when there is at least one root of Eq (2.6) across the imaginary axis on the complex plane. After excluding Turing instability induced by diffusion, now we seek the critical values of $\tau$ where the roots of Eq (2.6) will cross the imaginary axis from the left half plane to the right half plane. Plugging $\lambda = \mathrm{i}\omega_{n} (\omega_{n} > 0)$ into Eq (2.6), we obtain

Separating the real and imaginary parts, we obtain

Squaring and adding both equations of (2.8), we have

The number of positive roots of Eq (2.9) is relevant to the signs of $T^{2}_{n}-2D_{n}$ and $D^{2}_{n}-M^{2}$. For the convenience of discussion, we make the following assumptions.

Before the discussion of these cases, we need to investigate some properties of $T^{2}_{n}-2D_{n}$ and $D^{2}_{n}-M^{2}$.

Lemma 1. Suppose that $(\mathrm{H1})$ and $(\mathrm{H2})$ hold. $T^{2}_{n}-2D_{n}$ and $D^{2}_{n}-M^{2}$ are both monotonically increasing with respect to $n$.

Proof. When $(\mathrm{H2})$ holds, we have

which means that $T^{2}_{n}-2D_{n}$ is monotonically increasing with respect to $n$. Similarly, when $J_{11} < 0$, we have $D_{n} > 0$, and

□

2.1.   Case I when $(\mathrm{H3})$ holds

We first consider the case when $(\mathrm{H3})$ holds.

Lemma 2. If $(\mathrm{H1})$, $(\mathrm{H2})$ and $(\mathrm{H3})$ hold, then all roots of Eq (2.6) have negative real parts for all $\tau\geq0$.

Proof. Noticing that $D_{n}+M > 0$ for all $n\in\mathbb{N}_{0}$, thus $D_{0}-M\geq0$ leads to $D^{2}_{0}-M^{2}\geq0$. From Lemma 1 and $(\mathrm{H3})$, we can deduce that $T^{2}_{n}-2D_{n}\geq0$ and $D^{2}_{n}-M^{2}\geq0$ for all $n\in\mathbb{N}_{0}$. Therefore, Eq (2.9) has no positive roots, and Eq (2.6) has no imaginary roots. □

Theorem 2. If $(\mathrm{H1})$, $(\mathrm{H2})$ and $(\mathrm{H3})$ hold, the positive equilibrium $E^{*}$ is locally asymptotically stable for all $\tau\geq0$.

2.2.   Case II when $(\mathrm{H4})$ holds

Lemma 3. If $(\mathrm{H1})$, $(\mathrm{H2})$ and $(\mathrm{H4})$ hold, there exists $n_{1}\in\mathbb{N}_{0}$ such that the following conclusions hold.

(i) If either $n > n_{1}$ or $0\leq n\leq n_{1}$ and $\Delta = (T^{2}_{n}-2D_{n})^{2}-4(D^{2}_{n}-M^{2}) < 0$, then all roots of Eq (2.6) have negative real parts for all $\tau\geq0$.

(ii) If $\Delta > 0$ for $0\leq n\leq n_{1}$, then Eq (2.6) has a pair of imaginary roots $\pm\mathrm{i}\omega_{n}^{+}(\pm\mathrm{i}\omega_{n}^{-}, \ \mathrm{respectively})$ for $0\leq n\leq n_{1}$.

Proof. Similar as the proof in Lemma 2, we have $D^{2}_{n}-M^{2} > 0$ for all $n\in\mathbb{N}_{0}$. When $T^{2}_{0}-2D_{0} < 0$, from Lemma 1, there exists $n_{1}\in\mathbb{N}_{0}$ such that

It means that Eq (2.9) has no positive roots for $n > n_{1}$. For $0\leq n\leq n_{1}$, if $\Delta = (T^{2}_{n}-2D_{n})^{2}-4(D^{2}_{n}-M^{2}) < 0$, Eq (2.9) has no positive roots, and if $\Delta > 0$, Eq (2.9) has two positive roots $\omega^{\pm}_{n}$, where

Correspondingly, Eq (2.6) has a pair of purely imaginary roots $\pm\mathrm{i}\omega^{+}_{n}(\pm\mathrm{i}\omega^{-}_{n}, \; \mathrm{respectively})$ for $0\leq n\leq n_{1}$. □

Since $(\mathrm{H2})$ holds, we have $T_{n} > 0$. Noticing that $M = -J_{21}K_{12} > 0$, from (2.8), we have $S_{n}(\omega^{\pm}_{n}) > 0$, and we get

Lemma 4. Assume that $(\mathrm{H1})$, $(\mathrm{H2})$ and $(\mathrm{H4})$ hold. If $\Delta > 0$, then $\mathrm{Re}\left[\frac{\mathrm{d}\lambda}{\mathrm{d}\tau}\right]_{\tau = \tau_{n}^{j+ }} > 0$, $\mathrm{Re}\left[\frac{\mathrm{d}\lambda}{\mathrm{d}\tau}\right]_{\tau = \tau_{n}^{j- }} < 0$, for $0\leq n\leq n_{1}$, $j = 0, 1, 2, \cdots$, where $n_{1}$ is defined in (2.10).

Proof. Differentiating two sides of (2.6) with respect to $\tau$, we obtain

The sign of $\mathrm{Re}\left[\frac{\mathrm{d}\lambda}{\mathrm{d}\tau}\right]_{\tau = \tau_{n}^{j\pm }}$ is the same as that of $\mathrm{Re}\left[\frac{\mathrm{d}\lambda}{\mathrm{d}\tau}\right]^{-1}_{\tau = \tau_{n}^{j\pm }}$, thus we have $\mathrm{Re}\left[\frac{\mathrm{d}\lambda}{\mathrm{d}\tau}\right]_{\tau = \tau_{n}^{j+ }} > 0$, $\mathrm{Re}\left[\frac{\mathrm{d}\lambda}{\mathrm{d}\tau}\right]_{\tau = \tau_{n}^{j- }} < 0$. □

Now we consider the order of the critical values Hopf bifurcation.

Lemma 5. Assume that $(\mathrm{H1})$, $(\mathrm{H2})$ and $(\mathrm{H4})$ hold, and $\Delta > 0$ for $0\leq n\leq n_{1}$. $\tau_{n}^{j+}$ is monotonically increasing and $\tau_{n}^{j-}$ is monotonically decreasing with respect to $n$ ($0 < n\leq n_1$), where ${n_1}$ and $\tau _{n}^ {j\pm}$ are defined in (2.10) and (2.12), respectively.

Proof. See Appendix A. □

Theorem 3. Assume that $(\mathrm{H1})$, $(\mathrm{H2})$ and $(\mathrm{H4})$ hold, and $n_{1}$ is defined in (2.10).

(i) If $\Delta < 0$ for $0\leq n\leq n_{1}$, $E^{*}$ is locally asymptotically stable for all $\tau > 0$.

(ii) If $\Delta > 0$ for $0\leq n\leq n_{1}$, system (1.3) undergoes a Hopf bifurcation at $E^{*}$ when $\tau = \tau^{j\pm}_{n} \; (0\leq n\leq n_{1}, \; j = 0, 1, 2, \cdots)$. Moreover, $E^{*}$ is locally asymptotically stable when $\tau\in[0, \tau_{0}^{0+})\cup(\tau_{0}^{0-}, \tau_{0}^{1+})\cup\cdots\cup(\tau_{0}^{(s-1)-}, \tau_{0}^{s+})$, and it is unstable when $\tau\in(\tau_{0}^{0+}, \tau_{0}^{0-})\cup(\tau_{0}^{1+}, \tau_{0}^{1-})\cup\cdots\cup(\tau_{0}^{s+}, +\infty)$.

Proof. (ⅰ) From Lemma 3, when $\Delta < 0$, Eq (2.6) has no purely imaginary roots, thus the stability of $E^*$ can not be changed by delay.

(ⅱ) When $\Delta > 0$, from Lemma 5, we have $\tau_{0}^{j+} < \tau_{1}^{j+} < \cdots < \tau_{n_{1}}^{j+}$, and $\tau_{n_{1}}^{j-} < \cdots < \tau_{0}^{j-}, \; j = 0, 1, 2, \cdots.$ From $\omega_{n_{1}}^{+} > \omega_{n_{1}}^{-}$, we have $\arccos \left[\frac{(\omega_{n_{1}}^{+})^{2}-D_{n_{1}}}{M}\right] < \arccos\left[\frac{(\omega_{n_{1}}^{-})^{2}-D_{n_{1}}}{M}\right]$, and $\tau_{n_{1}}^{j+} < \tau_{n_{1}}^{j-}$. Therefore, $\tau_{0}^{j+} < \cdots < \tau_{n_{1}}^{j+} < \tau_{n_{1}}^{j-} < \cdots < \tau_{0}^{j-}$.

We claim that there exists a positive integer $s$ such that $\tau_{0}^{0+} < \tau_{0}^{0-} < \tau_{0}^{1+} < \tau_{0}^{1-} < \cdots < \tau_{0}^{s+} < \tau_{0}^{(s+1)+} < \tau_{0}^{s-}$. Since $\omega_{0}^{-} < \omega_{0}^{+}$, $\tau_{0}^{j+}-\tau_{0}^{(j-1)+} = \frac{2\pi}{\omega_{0}^{+}} < \tau_{0}^{j-}-\tau_{0}^{(j-1)-} = \frac{2\pi}{\omega_{0}^{-}}$, which means that the alternation for $\tau_{0}^{j+}$ and $\tau_{0}^{j-}$ cannot persist for the entire sequence, and there exists a positive integer $s$ such that $\tau_{0}^{(s-1)+} < \tau_{0}^{(s-1)-} < \tau_{0}^{s+} < \tau_{0}^{(s+1)+} < \tau_{0}^{s-}$. From Lemma 4, we can get the stable and unstable interval for $\tau$. □

2.3.   Case III when $(\mathrm{H5})$, $(\mathrm{H6})$ or $(\mathrm{H7})$ holds

First, we consider the case when $(\mathrm{H5})$ holds. If $(\mathrm{H5})$ holds, Eq (2.9) has a positive root $\omega_{0}^{+} = 2D_{0}-T_{0}^{2}$ when $n = 0$. Similar as Lemma 3, there exists $n_{1}\in\mathbb{N}_{0}$ such that (2.10) holds. If $\Delta > 0$ for $0 < n\leq n_{1}$, Eq (2.6) has a pair of imaginary roots $\pm\mathrm{i}\omega_{n}^{+}(\pm\mathrm{i}\omega_{n}^{-}, \mathrm{respectively})$ for $0 < n\leq n_{1}$.

Lemma 6. If $(\mathrm{H1})$, $(\mathrm{H2})$ and $(\mathrm{H5})$ hold, there exists $n_{1}\in\mathbb{N}_{0}$ defined in (2.10).

$(i)$ If either $n > n_{1}$ or $0 < n\leq n_{1}$ and $\Delta < 0$, then all roots of Eq (2.6) have negative real parts for all $\tau\geq0$.

$(ii)$ If $n = 0$, then Eq (2.6) has a pair of imaginary roots $\pm\mathrm{i}\omega_{0}^{+}$.

$(iii)$ If $\Delta > 0$ for $0 < n\leq n_{1}$, then Eq (2.6) has a pair of imaginary roots $\pm\mathrm{i}\omega_{n}^{+}(\pm\mathrm{i}\omega_{n}^{-}, \ \mathrm{respectively})$ for $0 < n\leq n_{1}$.

Remark 1. Similar as the proof of Lemma 5 and Theorem 3, we can get $\tau_{0}^{j+} < \tau_{1}^{j+} < \cdots < \tau_{n_{1}}^{j+} < \tau_{n_{1}}^{j-} < \cdots < \tau_{1}^{j-}$, where $\tau_n^{j\pm}$ is defined in (2.12). Similar as Lemma 4, we can verify that $\mathrm{Re}\left[\frac{\mathrm{d}\lambda}{\mathrm{d}\tau}\right]_{\tau = \tau_{n}^{j+ }} > 0$ for $0\leq n\leq n_{1}$ and $\mathrm{Re}\left[\frac{\mathrm{d}\lambda}{\mathrm{d}\tau}\right]_{\tau = \tau_{n}^{j- }} < 0$ for $0 < n\leq n_{1}$. Being different from case II, stability switches can not happen.

Now we consider the case when $(\mathrm{H6})$ holds. From Lemma 1, and $T^{2}_{0}-2D_{0}\geq0$, we have $T^{2}_{n}-2D_{n}\geq0$ for all $n\in\mathbb{N}_{0}$. From Lemma 1, combining $D_{0}-M < 0$ with $D_{n}+M > 0$, we get

Equation (2.9) has only one positive root $\omega^{+}_{n}$ for $0\leq n\leq n_{2}$, where

Lemma 7. Assume that $(\mathrm{H1})$, $(\mathrm{H2})$ and $(\mathrm{H6})$ hold, and $n_{2}\in\mathbb{N}_{0}$ is defined in (2.13).

(i) If $n > n_{2}$, then all roots of Eq (2.6) have negative real parts for all $\tau\geq0$.

(ii) If $0\leq n\leq n_{2}$, then Eq (2.6) has a pair of imaginary roots $\pm\mathrm{i}\omega_{n}^{+}$ for $0\leq n\leq n_{2}$.

It is easy to get $S_{n}(\omega^{+}_{n}) > 0$, from (2.8), we have

When $(\mathrm{H7})$ holds, from Lemma 1, we can deduce that there exists $n_{1}$ and $n_{2}$ such that

If $n_{1}\leq n_{2}$, Eq (2.9) has only one positive root $\omega^{+}_{n}$ for $0\leq n\leq n_{2}$, where

and Eq (2.9) has no positive roots for $n > n_{2}$.

If $n_{2}\leq n_{1}$, Eq (2.9) has one positive root $\omega^{+}_{n}$ for $0\leq n\leq n_{2}$, and has two positive roots $\omega^{\pm}_{n}$ for $n_{2}+1 < n\leq n_{1}$, and has no positive roots for $n > n_{1}$. When $n = n_{2}+1$, if $D_{n_{2}+1}^{2}-M^{2} > 0$, Eq (2.9) has two positive roots $\omega^{\pm}_{n_{2}+1}$, and if $D_{n_{2}+1}^{2}-M^{2} = 0$, Eq (2.9) has one positive root $\omega^{+}_{n_{2}+1}$. Similar as the proof in Lemma 5, we have

where

For either $n_{1}\leq n_{2}$ or $n_{2}\leq n_{1}$, we have

From the previous discussion, we can get the following conclusions.

Theorem 4. Assume that $(\mathrm{H1})$ and $(\mathrm{H2})$ hold, and $n_{1}$, $n_{2}$ are defined in (2.15). If $(\mathrm{H5})$, $(\mathrm{H6})$ or $(\mathrm{H7})$ holds, $E^{*}$ is locally asymptotically stable when $\tau\in[0, \tau_{0}^{0+})$, and $E^{*}$ is unstable when $\tau > \tau_{0}^{0+}$.

Proof. If $(\mathrm{H5})$, $(\mathrm{H6})$ or $(\mathrm{H7})$ holds, we easily know that $\tau_{0}^{0+}$ is the smallest critical value of Hopf bifurcation values. Similar as Lemma 4, we can get the transversality condition, which can lead to the results. □

3.   Double Hopf bifurcation

In order to figure out the joint effect of fear level $k$ and fear response delay $\tau$ on the dynamical behavior of system (1.3), we carry out double Hopf bifurcation analysis. We first verify the condition for the occurrence of double Hopf bifurcation.

Remark 2. Assume that $(\mathrm{H1})$, $(\mathrm{H2})$ and $(\mathrm{H4})$ hold. If there exists $k = k^{*}$ such that $\tau^{j_{1}+}_{n_{1}} = \tau^{j_{2}-}_{n_{2}} = \tau^{*}$. Then system (1.3) may undergo a double Hopf bifurcation at the positive equilibrium $E^{*}(u^{*}, v^{*})$ when $(k, \tau) = (k^{*}, \tau^{*})$.

Define the real-valued Hilbert space $X: = \left\{(u, v)\in H^{2}(0, l\pi)\times H^{2}(0, l\pi):(u_{x}, v_{x})\big|_{x = 0, l\pi} = 0\right\}$, and the corresponding complex space $X_{\mathbb{C}}: = \left\{x_{1}+\mathrm{i}x_{2}, x_{1}, x_{2}\in X\right\}$. Let $\widetilde{u} = u(x, \tau t)-u^{*}$, $\widetilde{v} = v(x, \tau t)-v^{*}$, and drop the hats for the convenience of notation. Denote $U(t) = (u(x, t), v(x, t))^{T}$, and $U_{t}(\theta) = U(t+\theta), -1\leq\theta\leq0$. System (1.3) can be written as

where the coefficients of nonlinear terms $F_{uu}, F_{uv}$, etc., are in Appendix B, and

Here $J_{11}(k), J_{12}(k), J_{21}(k)$, and $K_{12}(k)$ are defined in (2.5), in which $u^{*}, v^{*}$ are functions of $k$ defined in (2.1).

Let $k = k^{*}+\alpha_{1}$, $\tau = \tau^{*}+\alpha_{2}$, then $(\alpha_{1}, \alpha_{2}) = (0, 0)$ is a double Hopf bifurcation point. System (3.1) can be written as

where

with

and

We leave the detailed calculation and some expressions in Appendix C.

Define an enlarged phase space $\mathcal{BC}$ by

Equation (3.2) can be written as an abstract ordinary differential equation in $\mathcal{BC}$

where $A\varphi = \dot{\varphi}+X_{0}[D_{0}\Delta\varphi(0)+L_{0}(\varphi)-\dot{\varphi}(0)]$, with

It is obvious that $\{\pm\mathrm{i}\omega^{+}_{n_{1}}\tau^{*}, \pm\mathrm{i}\omega^{-}_{n_{2}}\tau^{*}\}$ are the pure imaginary eigenvalues of the operator $A$ and the corresponding eigenfunctions in $\mathcal{BC}$ are $\{\phi_{1}(\theta)\gamma_{n_{1}}, \phi_{2}(\theta)\gamma_{n_{1}}, \phi_{3}(\theta)\gamma_{n_{2}}, \phi_{4}(\theta)\gamma_{n_{2}}\}$, respectively, where

and $\phi_{1}(\theta) = (1, p_{12})^{T}e^{\mathrm{i}\omega^{+}_{n_{1}}\tau^{*}\theta}$, $\phi_{3}(\theta) = (1, p_{32})^{T}e^{\mathrm{i}\omega^{-}_{n_{2}}\tau^{*}\theta}$, $\phi_{2}(\theta) = \overline{\phi_{1}(\theta)}$, $\phi_{4}(\theta) = \overline{\phi_{3}(\theta)}$, with

A direct calculation derived that $\psi_{1}(s) = D_{1}(1, q_{12})e^{-\mathrm{i}\omega^{+}_{n_{1}}\tau^{*}s}$, $\psi_{3}(s) = D_{3}(1, q_{32})e^{-\mathrm{i}\omega^{-}_{n_{2}}\tau^{*}s}$ are the formal adjoint eigenvectors of $\phi_{1}(\theta)$, $\phi_{3}(\theta)$ and satisfy $(\phi_{1}, \psi_{1})_{1} = 1$, $(\phi_{3}, \psi_{3})_{2} = 1$. Here

and $(\cdot, \cdot)_{k}$ are the adjoint bilinear form

with a function of bounded variation $\eta_{k}(\theta, \tau^{*}), \; \theta \in[-1, 0]$ are given by

Denote

where $\psi_{2}(s) = \overline{\psi_{1}(s)}$, $\psi_{4}(s) = \overline{\psi_{3}(s)}$.

Now we decompose $\mathcal{BC}$ into the direct sum of the central subspace $\mathcal{P}$ and its complementary space $\mathrm{Ker}\pi$

where $\pi:\mathcal{BC}\rightarrow \mathcal{P}$ is the projection defined by

According to (3.4), $U_{t}$ can be decomposed as

where $(z_{1}, z_{2}, z_{3}, z_{4})^{T} = (\tilde{z}_{1}, \tilde{z}_{2}) = (\Psi_{k}, \langle U_{t}, \beta_{n_{k}}\rangle)_{k}$, and $y_{t}(\theta)\in \mathrm{Ker}\pi$.

Then system (3.3) in $\mathcal{BC}$ is equivalent to the system

where $A_{1}$ is the restriction of $A$ on $\mathcal{Q}^{1}\subset\mathrm{Ker}\pi\to \mathrm{Ker}\pi$, $A_{1}\varphi = A\varphi$ for $\varphi\in\mathcal{Q}^{1}$, and $\widetilde{F}_{2}(z, y, \alpha)$ and $\widetilde{F}_{3}(z, 0, 0)$ are in Appendix D.

According to ^[28], we can obtain the normal form of the double Hopf bifurcation up to the third order as follows

Here

where

with

4.   Numerical simulations

In this section, we carry out some numerical simulations to explain the theoretical results obtained in the previous sections. Symbolic mathematical software Matlab is used to plot numerical graphs. The delayed reaction-diffusion system (1.3) is numerically solved by transforming the continuous system to discrete system using discretization of time and space. In the discrete system, the Laplacian describing diffusion is calculated using finite differences, and the time evolution is solved using the Euler method. Fix the parameters as

4.1.   Three different cases induced by fear response delay $\tau$

From Theorems 2, 3, 4, we have found that the fear response delay may have no effect, both stabilizing and destabilizing effect, or destabilizing effect on the stability of $E^{*}$, respectively. In this section, we explore that the impact level of fear $k$ has on the stability of the positive equilibrium $E^{*}$ and the occurrence of Hopf bifurcation induced by fear response delay.

Choose $k = 0.04$, such that the level of fear is low. We can verify that $(\mathrm{H1})$ holds, thus $E^{*}(0.0659, 0.4483)$ exists. Moreover, $(\mathrm{H2})$ and $(\mathrm{H3})$ hold, which indicates that $E^{*}$ is locally asymptotically stable for all $\tau > 0$ from Theorem 2. This indicates that the fear response delay does not change the stability of the positive equilibrium, which means that the fear response delay cannot support periodic oscillations in prey and predator populations with a low level of fear.

When $k$ is increased to $k = 1$, the positive equilibrium $E^{*}(0.1299, 0.49)$ still exists, and $(\mathrm{H1})$, $(\mathrm{H2})$ and $(\mathrm{H4})$ hold. By calculating the critical value of Hopf bifurcation in Table 1, we can get $\tau^{0+}_{0} < \tau^{0+}_{1} < \tau^{0-}_{1} < \tau^{0-}_{0} < \tau^{1+}_{0} < \tau^{1+}_{1} < \tau^{2+}_{0} < \tau^{2+}_{1}$. From Theorem 3, $E^{*}$ is locally asymptotically stable when $\tau\in[0, \tau^{0+}_{0})\cup(\tau^{0-}_{0}, \tau^{1+}_{0})$ (see Figures 1 and 3), and it is unstable when $\tau\in(\tau^{0+}_{0}, \tau^{0-}_{0})\cup(\tau^{1+}_{0}, +\infty)$ (see Figures 2 and 4). It means that the fear response delay may lead to stability switches with an intermediate level of fear, which is not observed in the model with biomass transfer delay in ^[16].

When $k$ is further increased to $k = 2.5$, we can verify that $(\mathrm{H1})$, $(\mathrm{H2})$ and $(\mathrm{H7})$ hold, and we get $\tau^{0+}_{0} = 7.3868$. From Theorem 4, $E^{*}$ is locally asymptotically stable when $\tau\in[0, \tau^{0+}_{0})$, and it is unstable when $\tau > \tau^{0+}_{0}$.

4.2.   The joint effect of $k$ and $\tau$

In this section, we explore the joint effect of the fear level and fear response delay on the dynamics of system (1.3). We choose $k$ and $\tau$ as bifurcation parameters, and other parameters are the same as in (4.1). From (2.12), we can draw Hopf bifurcation curves with varying $k$ (see Figure 5). Hopf bifurcation curve $\tau_{0}^{1+}$ intersects with $\tau_{0}^{0-}$ at HH $(k^{*}, \tau^{*}) = (1.3677, 58.6677)$. On HH, Eq (2.6) has two pairs of purely imaginary roots $\pm\mathrm{i}\omega_{0}^{-} = \pm 0.0503\ \mathrm{i}, \ \pm\mathrm{i}\omega_{0}^{+} = \pm 0.1146\ \mathrm{i}$. Therefore, HH is a possible double Hopf bifurcation point. By the calculation of normal form of double Hopf bifurcation, we can get the coefficients in (3.5) with

Make the polar coordinate transformation

and denote

Removing the bars, the four-dimensional system (3.5) is transformed into the following two-dimensional amplitude system

By a simple calculation, we find that the system (4.2) admits the following equilibria

According to ^[29], the unfolding of system (4.2) is of type VIa. There are rays $L_{1}-L_{8}$ near the bifurcation point, which divide the phase space into eight parts (see Figure 6(a)), where

The detailed dynamics in $D_{1}-D_{8}$ near $(k^{*}, \tau^{*})$ have shown in Figure 6(b).

Notice that $A_{0}$ of (4.2) corresponds to the positive equilibrium $E^{*}$ of the original system (1.3). $A_{1}$ and $A_{2}$ correspond to the periodic solution of (1.3). $A_{3}$ corresponds to quasi-periodic solution on a 2-torus. Periodic orbit of (4.2) corresponds to quasi-periodic solution on a 3-torus. Due to the fact that double Hopf bifurcation point HH is the intersection of 0-mode Hopf bifurcation curves $\tau_{0}^{1+}$ and $\tau_{0}^{0-}$, all periodic or quasi-periodic solutions near HH are spatially homogeneous.

In region $D_{1}$, system (4.2) has only one equilibrium $A_{0}$, which is a saddle. It means that the positive equilibrium $E^{*}(u^{*}, v^{*})$ of system (1.3) is unstable.

In region $D_{2}$, the trivial equilibrium $A_{0}$ of (4.2) is a sink, and $A_{1}$ is unstable. This means that $E^{*}(u^{*}, v^{*})$ of system (1.3) is asymptotically stable as shown in Figure 7, when $k = 0.9$ and $\tau = 50$ are chosen in $D_{2}$. In addition, the spatially homogeneous periodic solution is born, which is unstable.

In region $D_{3}$, system (4.2) has three equilibria: $A_{0}$, $A_{1}$ and $A_{2}$. $A_{2}$ is stable while other equilibria are unstable. When $k = 0.8$ and $\tau = 55$ are chosen in $D_{3}$, system (1.3) has spatially homogeneous periodic solution, which is stable (see Figure 8).

In region $D_{4}$, there are four equilibria of (4.2): $A_{0}$, $A_{1}$, $A_{2}$ and $A_{3}$. $A_{3}$ is stable while other equilibria are unstable. Since the double Hopf bifurcation occurs at the intersection of two Hopf bifurcation curves with $n = 0$, we choose $x = 0$ to demonstrate the rich dynamical phenomena of system (1.3). When parameters are chosen as $k = 1.345$ and $\tau = 59.5$ in $D_{4}$, system (1.3) has a quasi-periodic solution on 2-torus, which is shown in Figure 9.

In region $D_{5}$, there are four equilibria and a periodic orbit of (4.2). $A_{0}$, $A_{1}$, $A_{2}$ and $A_{3}$ are all unstable, and the periodic orbit is stable. This means that system (1.3) has a stable quasi-periodic solution on a 3-torus. Figure 10 illustrates the existence of a stable quasi-periodic solution on a 3-torus when $k = 1.364$ and $\tau = 59.4$ are chosen in $D_{5}$.

In region $D_{6}$, system (4.2) has four equilibria: $A_{0}$, $A_{1}$, $A_{2}$ and $A_{3}$, which are all unstable. When the parameter enters region $D_{6}$, the 3-torus disappears through heteroclinic orbits, and the system (1.3) generates strange attractor according to the "the Ruelle-Takens-Newhouse" scenario. Figure 11 illustrates system (1.3), which demonstrates chaotic phenomena when $k = 1.3801$ and $\tau = 59.978$ are chosen in $D_{6}$. The right figures of Figure 9, 10 and 11 are the results of Poincaré map on a Poincaré section.

In region $D_{7}$, system (4.2) has three equilibria: $A_{0}$, $A_{1}$ and $A_{2}$, which are all unstable. This means that the positive equilibrium $E^{*}(u^{*}, v^{*})$ and the spatially homogeneous periodic solutions of system (1.3) are all unstable.

In region $D_{8}$, system (4.2) has two equilibria: $A_{0}$ and $A_{2}$, which are both unstable. This means that the positive equilibrium $E^{*}(u^{*}, v^{*})$ and the spatially homogeneous periodic solution of system (1.3) are all unstable.

5.   Conclusion and future work

There have been more and more evidences showing that fear from the predator may affect the behavior and physiology of prey, which could reduce the reproduction of prey ^[1,6,7]. Moreover, anti-predation behavior of scared prey can also reduce the chance of prey being caught by predators. Therefore, we consider a predator-prey model with both costs and benefits due to fear based on ^[16]. Since the cost of fear in prey reproduction is not instantaneous, we incorporate a fear response delay into the model. Our aim is to explore how the fear level and fear response delay jointly affect the dynamics of the predator-prey system from the point of view of both codimension-1 and codimension-2 bifurcation analysis.

First, we find that Turing instability does not occur in our model. Since the pioneering work of Turing ^[30], quite a number of literature has revealed that diffusion might destabilize an equilibrium and generate spatial pattern formation. However, not all diffusive systems exhibit such diffusion-driven instability. In some chemical and biological models, Turing instability induced by diffusion are not observed ^[31,32,33].

Next, we explore the role of fear response delay on the dynamics. Through the analysis of the characteristic equation, it turns out that there are three different cases: the fear response delay may have no effect, both stabilizing and destabilizing effect, or destabilizing effect on the stability of $E^{*}$, which are shown in Theorems 2, 3, 4, respectively. Combining with numerical simulations, it is found that the the occurrence of Hopf bifurcation induced by fear response delay has closely connection with the fear level. When the fear level is low, the delay could not change the stability of the positive equilibrium, which is asymptotically stable for all $\tau$. From the point of view of biology, the fear response delay cannot support periodic oscillations in prey and predator populations with a low level of fear. When the fear level is in an intermediate range, the fear response delay can induce stability switches, which indicates that the delay has effects on both destabilizing and stabilizing the dynamics. It is observed that such stability switches cannot be generated by the biomass transfer delay ^[16]. When the fear level is high, the fear response delay can destabilize the equilibrium and induce Hopf bifurcation.

Finally, to better reveal the joint effect of fear level and fear response delay on the dynamics of system (1.3), we carry out codimension-2 bifurcation analysis. On the $(k, \tau)$ plane, by finding the interactions of Hopf bifurcation curves, we can find the possible critical point of double Hopf bifurcation. After the calculation of the normal form for double Hopf bifurcation, we can obtain the bifurcation set, and figure out all the dynamical behaviors around the critical point. There are periodic solutions, quasi-periodic solutions and even chaos observed near the double Hopf bifurcation point.

In this paper, we assume that the interaction between individuals of a species is local. The individuals at different locations may compete for common resource or communicate visually or by chemical means. Regarding fear effect, the fear of predator depends on both the local and nearby appearance of the predator. Taking these into account, the influence of the nonlocal effect on the spatiotemporal dynamics of the diffusive predator-prey model with fear would be worth further study. In addition, we consider diffusion induced by random movement, however, cognition and memory also have an important influence on the animal movement. Shi et al. ^[34] proposed a memory-based diffusion equation of a single species. Memory-based diffusion can be modified and applied to the model with fear effect, which will provide much more realistic results and bring mathematical challenges.

Use of AI tools declaration

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

Acknowledgments

The work is supported by National Natural Science Foundation of China (Nos. $11901369$, $11971281$ and $12071268$) and PhD research startup foundation of Shaanxi University of Science and Technology (2023BJ-13).

Conflict of interest

The authors declare there is no conflicts of interest.

Appendix A.   The proof of Lemma 5

By direct calculation, we have

From Lemma 1, we have $\frac{\mathrm{d}(T^{2}_{n}-2D_{n})}{\mathrm{d}n} > 0$, $\frac{\mathrm{d}(D_{n}^{2}-M^{2})}{\mathrm{d}n} > 0$. Moreover, from (2.10), we know that $T^{2}_{n}-2D_{n} < 0$ for $0\leq n \leq n_{1}$. Combining $\frac{\mathrm{d}D_{n}}{\mathrm{d}n} > 0$ and $M > 0$, we have $\frac{\mathrm{d}[(\omega_{n}^{+})^{2}-D_{n}]}{\mathrm{d}n} < 0$, and $\frac{\mathrm{d}[\arccos(\frac{(\omega_{n}^{+})^{2}-D_{n}}{M})+2j\pi]}{\mathrm{d}n} > 0$. Similarly, we can also prove that $\frac{\mathrm{d}\omega_{n}^{+}}{\mathrm{d}n} < 0$. Thus, $\tau_{n}^{j+}$ is monotonically increasing with respect to $n$.

Now we consider the monotonicity for $\tau_{n}^{j-}$.

Since $D_{n}^{2}-M^{2} > 0$ for all $n\in\mathbb{N}_{0}$, and $T^{2}_{n}-2D_{n} < 0$ for $0\leq n \leq n_{1}$, then $\sqrt{\Delta}+T^{2}_{n}-2D_{n} < 0$ for $0\leq n \leq n_{1}$. Combining with Lemma 1, we have $-\frac{\mathrm{d}(T^{2}_{n}-2D_{n})}{\mathrm{d}n}\cdot\left[\frac{\sqrt{\Delta}+(T^{2}_{n}-2D_{n})}{\sqrt{\Delta}}\right] > 0$. From $T^{2}_{n}-2D_{n} < 0$, and $D_{n}^{2} > M^{2}$, we have $4D_{n}^{2}-\Delta = 4(D_{n}^{2}-M^{2})+T^{2}_{n}(4D_{n}-T^{2}_{n}) > 0$, and thus $2D_{n}-\sqrt{\Delta} > 0$. From $\frac{\mathrm{d}D_{n}}{\mathrm{d}n} > 0$, we get $\frac{1}{\sqrt{\Delta}}\cdot\frac{\mathrm{d}(D_{n}^{2}-M^{2})}{\mathrm{d}n}-\frac{\mathrm{d}D_{n}}{\mathrm{d}n} = \frac{\mathrm{d}D_{n}}{\mathrm{d}n}\cdot\left[\frac{2D_{n}-\sqrt{\Delta}}{\sqrt{\Delta}}\right] > 0$. Therefore, $\frac{\mathrm{d}[(\omega_{n}^{-})^{2}-D_{n}]}{\mathrm{d}n} > 0$, and $\frac{\mathrm{d}[\arccos(\frac{(\omega_{n}^{-})^{2}-D_{n}}{M})+2j\pi]}{\mathrm{d}n} < 0$. $\frac{\mathrm{d}[(\omega_{n}^{-})^{2}-D_{n}]}{\mathrm{d}n} > 0$ can also lead to $\frac{\mathrm{d}\omega_{n}^{-}}{\mathrm{d}n} > 0$. Thus, $\tau_{n}^{j-}$ is monotonically decreasing with respect to $n$.

Appendix B.   The coefficients of nonlinear terms

Appendix C.   Some terms in $F_{2}(\alpha_{1}, \alpha_{2}, U_{t})$

From the expression of $F_{2}(\alpha_{1}, \alpha_{2}, U_{t})$,

where

with $x = 1+c_{1}k^{*}$, $y = 1+qu^{*}(k^{*})$, $z = 1+c_{2}k^{*}v^{*}(k^{*})$, $h = c_{1}y+2qu^{*^{\prime}}(k^{*})x$,

Appendix D.   $\widetilde{F}_{2}(z, y, \alpha)$ and $\widetilde{F}_{3}(z, 0, 0)$

Here $S_{yz_{i}}(i = 1, 2, 3, 4)$ are linear operators and

where

and

Denote $f_{q_{1}q_{2}q_{3}q_{4}}^{1(k)} = \sqrt{\frac{1}{l\pi}}\psi_{k}(0)F_{q_{1}q_{2}q_{3}q_{4}}, (q_{1}+q_{2}+q_{3}+q_{4} = 2)$.