Hopf bifurcation induced by fear: A Leslie-Gower reaction-diffusion predator-prey model

1.   Introduction

Recently, one of the hot issues in the ecosystem being studied is the impact of fear effect on the dynamics of prey and predators in predator-prey models. The fear effect is an inherent psychological reaction of the organisms to increase alertness and respond to danger^[1]. It can trigger various anti-predation responses, such as changing the reproduction capacity and strategies ^[2], changing foraging behaviors and selecting new habitats ^[3], and reducing the birth and survival rate of offspring ^[4]. In 2020, Sarkara and Khajanchi ^[5] proposed the following predator-prey model that introduces the cost of fear into the prey

where $u$ and $v$, respectively, represent the densities of prey and predator. The meaning of the parameters can refer to ^[5]. In addition, in ^[5], $g(\eta, \alpha, v) = \eta+\frac{\alpha(1-\eta)}{\alpha+v}$ stands for a fear function which describes that the prey is affected by the fear of predator. So, we know some characteristics of $g(\eta, \alpha, v)$ showing $\lim_{\alpha\rightarrow \infty}g(\eta, \alpha, v) = 1$, $\lim_{v\rightarrow \infty}g(\eta, \alpha, v) = \eta$, $g(0, \alpha, v) = \frac{\alpha }{\alpha+v}$, $g(\eta, 0, v) = \eta$, $g(\eta, \alpha, 0) = 1$, $g(1, \alpha, v) = 1$, $\frac{\partial g}{\partial \eta} > 0$, $\frac{\partial g}{\partial \alpha} > 0$, $\frac{\partial g}{\partial v} < 0$, which imply that the inhibitory effect of fear on the birth rate of prey will increase with the increase of predator population and will decrease with the enhancement of the ability to identify the capture of predator.

In model (1.1), $\frac{a u}{1+bu}$ stands for the Holling Ⅱ functional response ^[6], which expresses the prey consumed by each predator per unit time, which only depends on the density of prey without being disturbed by the predator. It is of great significance to use different functional response functions to describe the relationship between prey and predator, which is caused by the difficulty of capturing the prey by the predator and the different absorption conversion rates of the predator. As we all know, the Beddington-DeAngelis (B-D) functional response function ^[7,8] is described as $\frac{u}{1+au+bv}$, which depends on both the density of prey and predator populations. Assume that $b = 0$, the B-D functional response becomes the Holling Ⅱ functional response, which implies that it is more comprehensive and accurate in describing the interference and handling of populations. In 1960, Leslie and Gower ^[9] proposed a novel functional response $\frac{v}{cu}$ to describe the conversion rate of prey to predator, called the Leslie-Gower (L-G) functional response. Compared with the Holling Ⅱ functional response, the L-G functional response is affected by prey and interfered with by predators ^[10]. Therefore, the Holling Ⅱ functional response function in model (1.1) can be replaced by the B-D functional response and L-G functional response, respectively, to study further the impact of different functional response functions on populations. Therefore, we consider the following predator-prey model, in which prey is subject to B-D functional response and predator is subject to modified L-G functional response

where the meaning of the parameters are the same as model (1.1), $b$ is the growth rate of predator, and $\delta$ stands for the strength of prey refuge ^[11] with $\delta\in[0, 1)$. $\frac{u}{a_1+u+ev}$ represents the B-D functional response. $\frac{cv}{a_2+ u}$ represents a modified L-G functional response, where $c$ is the maximum value of the diminishment rate of predator due to prey, and $a_2$ measures the extent to which environment provides predator.

To explore how the fear effect acts on the populations, many scholars have constructed and studied a large number of models with different functional response functions. For example, Wang et al. ^[12] as well as Sarkara and Khajanchi ^[5], both proposed a predator-prey model with fear and Holling Ⅱ. They all found that prey and predator populations were affected by fear effect. Pal et al. ^[13] constructed a B-D predator-prey model with fear, which found that fear has a destructive effect on stability. Wang et al. ^[14] investigated an improved L-G predator-prey model with fear. They found that as the degree of fear increases, it will lead to a decrease in population density and the extinction of prey.

For simplicity, let $x = \frac{c}{a_1}u, y = \frac{c}{a_2b}v, \tau = bt$, and model (1.1) can transform into the following simplified model (for simplicity, $u$, $v$, $t$ represent $x$, $y$, $\tau$ again, respectively):

where $r = \frac{r_0\eta-d_0}{b}$, $\gamma = \frac{a_1\beta}{bc}$, $\theta = \frac{r_0(1-\eta)}{b}$, $K = \frac{a_2b}{c\alpha}$, $p = \frac{a_1c}{a_2(1-\delta)}$, $q = \frac{a_1(1-\delta)}{a_2c}$, $h = \frac{a_1}{a_2}$, $m = be$. Here, $r$ can be positive or negative. Now, we give an analysis of the nature of $r$ as follows:

(1) When $r_0\eta > d_0$, that is, $1 > \eta > \frac{d_0}{r_0}$, we have $r > 0$, which means that when the cost of minimum fear is high, the birth rate of prey affected by fear is higher than the mortality rate of prey.

(2) When $r_0\eta < d_0$, that is, $\frac{d_0}{r_0} > \eta > 0$, one has $r < 0$, which means that when the cost of minimum fear is small, the birth rate of prey affected by fear effect is lower than the mortality rate of prey.

(3) When $\eta = 1$, i.e., without the effect of the fear effect, we have $r_0 > d_0$, which means that the birth rate is higher than the mortality rate, which is consistent with the actual situation of biological species in the ecosystem.

(4) $\theta+r$ represents the natural growth rate of prey, which is greater than 0, meaning that prey will survive for a long time.

Then, some natural questions arise from model (1.3):

● What are the conditions for the existence and stability of the equilibria?

● What are the conditions for the occurrence of Hopf bifurcation, and, if so, how to determine its stability and direction?

On the other hand, species not only evolve on the timescale but also move randomly on the spatial scale. So, it is inevitable to consider the issues of species in time and space. In 1952, Turing ^[15] first described the movement of species on time and space scales via using the reaction-diffusion equations. He found that the steady state equilibrium in the spatial model is stable in the absence of diffusion but transforms unstable in the presence of diffusion, which means that diffusion can induce instability of populations ^[16,17]. In addition to the instability driven by diffusion, the reaction-diffusion system also can be triggered by other mechanisms, such as steady states solutions ^[18], Hopf bifurcation ^[19], pattern formation and etc ^[20]. Han et al. ^[21] considered a modified L-G predator-prey model with cross-diffusion and indirect predation effect, which shows that cross-diffusion can drive Turing instability and pattern formation. Tiwari et al. ^[22] proposed and analyzed a B-D predator-prey interaction model with fear. They investigated some properties of bifurcation, such as Hopf bifurcation and pitchfork bifurcation.

Motivated by these works, we construct a reaction-diffusion predator-prey model as follows

where $d_{1}\geq 0$ and $d_{2}\geq 0$ are, respectively, the diffusion rate of prey and predator. Laplacian operator $\Delta = \frac{\partial^2}{\partial x^2}$ is in the one-dimensional diffusion, $\Delta = \frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}$ is in the two-dimensional diffusion, $\Delta = \frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2}$ is in the three-dimensional diffusion, $\Omega\subset\mathbb{R}^n$ is a bounded domain with smooth boundary $\partial \Omega$, and $\mathbf{n}$ is the outward unit normal vector of the boundary $\partial \Omega$ as in ^[23]. There is no population flux across the boundaries owing to homogeneous Neumann boundary conditions. Then, we ask:

● What are the critical conditions that determine Turing instability?

● What are the conditions for the occurrence of the spatially homogeneous and spatial inhomogeneous periodic solutions?

● If spatial homogeneous and spatial inhomogeneous periodic solutions occur, what are the conditions for determining the stability and direction?

The rest of the paper is organized as follows. In Section 2, we focus on discussing the existence and stability of the equilibria and give the Hopf bifurcation analysis of the nonspatial model (1.3). In Section 3, we investigate the Hopf bifurcation analysis of the spatial model (1.4). In Section 4, we present a series of numerical simulations to reveal the theoretical analysis.

2.   Analysis of the nonspatial model (1.3)

2.1.   Positivity and boundedness

From the nonspatial model (1.3), one has

We know $u(t), v(t)\geq0$ based on the above two expressions. Then, $\mathbb{R}_{+}^2 = \{u > 0, v > 0\}$ is positively invariant for the nonspatial model (1.3).

Lemma 2.1. All solutions $(u(t), v(t))$ of model (1.3) are contained in the region $U = \bigg\{(u, v) \in \mathbb{R}^2_+:0\leq u(t)\leq \frac{r+\theta}{\gamma}, 0\leq v(t)\leq 1+\frac{q(r+\theta)}{\gamma} \bigg\}$ as $t\rightarrow +\infty$.

2.2.   Existence of equilibria

Define thresholds

From model (1.3), we denote

Clearly, model (1.3) has three boundary equilibrium points: $E_0 = (0, 0)$, $E_1 = \left(\frac{\theta +r}{\gamma}, 0\right)$, $E_2 = \left(0, 1\right)$, where $E_0$ and $E_2$ always exist, and $E_1$ existes when $R_0 > 0$.

By solving (2.1), we obtain that $v = 1+qu$ and

where

It is obvious that Eq (2.2) is a third-order algebraic equation, which has one, two, or three positive roots. Hence, discussing the number of positive equilibria of model (1.3) is equivalent to discussing the number of positive roots of Eq (2.2).

First, let

From Eq (2.2), we know $A_0 > 0$, and the signs of $A_1$, $A_2$ and $A_3$ are uncertain. First, we discuss the sign of $A_3$.

(1) When $R_0 > 0$ and $r > 0$ hold, that is, $\theta +r > 0$ and $1 > \eta > \frac{d_0}{r_0}$, if $\frac{1}{p+m} > r$, we know that $1-r(p+m) > 0$, which can be obtained that $A_3 > 0$ when $\theta < \theta^*$, $A_3 < 0$ when $\theta > \theta^*$, or $A_3 = 0$ when $\theta = \theta^*$; if $\frac{1}{p+m}\leq r$, we know that $1-r(p+m)\leq0$, which we can get that $A_3 < 0$.

(2) When $R_0 > 0$ holds, if $r < 0$, that is, $\frac{d_0}{r_0} > \eta > 0$, which implies that $1-r(p+m) > 0$, then we get that $A_3 > 0$ when $\theta < \theta^*$, $A_3 < 0$ when $\theta > \theta^*$, or $A_3 = 0$ when $\theta = \theta^*$.

(3) When $R_0 > 0$ and $r = 0$ hold, we know that $A_3 > 0$ when $\frac{1+K}{p+m} > \theta > 0$, that is, $\theta < \theta^*$, $A_3 < 0$ when $\theta > \theta^*$, or $A_3 = 0$ when $\theta = \theta^*$.

Therefore, we get the following lemma:

Lemma 2.2. Regarding the sign of $A_3$, we have

(1) $A_3\geq0$ if $(\mathbf{S_1}): R_0 > 0, \; \theta\leq\theta^*, \; \frac{1}{p+m} > r;$

(2) $A_3 < 0$ if $(\mathbf{S_2}): R_0 > 0, \; \theta > \theta^*, \; \frac{1}{p+m} > r;$

(3) $A_3 < 0$ if $(\mathbf{S_3}): R_0 > 0, \; r > 0, \; \frac{1}{p+m}\leq r.$

From Lemma 2.2(1), we know that $A_3\geq0$ when $(\mathbf{S_1})$ holds, which means that Eq (2.2) has at most two positive roots (see Figure 1). Next, by discussing the sign of $\Delta_2$, we obtain the number of positive roots of Eq (2.2), seeing the following lemma:

Lemma 2.3. Suppose that $(\mathbf{S_1})$ holds. Then, we have

(1) if $\Delta_2 < 0$, then Eq (2.2) has no positive root;

(2) if $\Delta_2 = 0$, then Eq (2.2) has a unique positive root;

(3) if $\Delta_2 > 0$, then Eq (2.2) has two different positive roots.

From Lemma 2.2(2), (3), we know that $A_3 < 0$ when $(\mathbf{S_2})$ or $(\mathbf{S_3})$ holds, which means that Eq (2.2) has at least one positive root and at most three positive roots (see Figure 2). Then, by discussing the sign of $\Delta_1$ and $\Delta_2$, one has the following lemma:

Lemma 2.4. Suppose that $(\mathbf{S_2})$ or $(\mathbf{S_3})$ holds. We obtain

(1) if $\Delta_2 < 0$, then Eq (2.2) has a unique positive root;

(2) if $\Delta_2 = 0$, and

(i) $\Delta_1 = 0$, then Eq (2.2) has a unique positive root;

(ii) $\Delta_1 > 0$, then Eq (2.2) has two positive roots;

(3) if $\Delta_2 > 0$, then Eq (2.2) has three different positive roots.

Therefore, summarizing the above lemmas, we obtain the following theorem (for convenience, we summarize Theorem 2.1 in Table 1).

Theorem 2.1. (Ⅰ) Model (1.3) has a unique trivial equilibrium $E_0 = (0, 0)$ and one semi-trivial equilibrium $E_2 = \left(0, 1\right)$;

(Ⅱ) if $R_0 > 0$, model (1.3) also has a semi-trivial equilibrium $E_1 = (u_1, 0) = \left(\frac{\theta +r}{\gamma}, 0\right)$;

(Ⅲ) when $(\mathbf{S_1})$ holds, model (1.3) has at most two positive equilibria. Moreover,

(1) if $\Delta_2 > 0$, then model (1.3) has two different positive equilibria $E^*_i(i = 2, 3)$;

(2) if $\Delta_2 = 0$, then model (1.3) has a unique positive equilibrium $E_*$;

(3) if $\Delta_2 < 0$, then model (1.3) has no positive equilibrium;

(Ⅳ) when $(\mathbf{S_2})$ or $(\mathbf{S_3})$ holds, model (1.3) has at least one positive equilibrium and at most three positive equilibria. Moreover,

(1) if $\Delta_2 > 0$, then model (1.3) has three different positive equilibria $E^*_i(i = 1, 2, 3)$;

(2) if $\Delta_2 = 0$, and

(i) $\Delta_1 > 0$, then model (1.3) has two positive equilibria $E_*$ and $E^*_1(or E^*_3)$;

(ii) $\Delta_1 = 0$, then model (1.3) has a unique positive equilibrium $E^*_* = \left(-\frac{A_1}{3A_0}, 1-\frac{q A_1}{3A_0}\right)$;

(3) if $\Delta_2 < 0$, then model (1.3) has a unique positive equilibrium $E^*_3$.

2.3.   Stability of equilibria

For model (1.3), there is at least one positive equilibrium and at most three. Consequently, there are numerous rich and complex dynamical characteristics, which increases the difficulty in studying the dynamical properties of this model. Therefore, in order to alleviate the complexity of this study, we will only consider the case with a unique positive equilibrium, denoted as $E^*$. So in the subsequent research, we consistently assume that model (1.3) has a unique positive equilibrium, namely, $E^*(u^*, v^*)$.

Theorem 2.2. (Ⅰ) $E_0$ and $E_1$ are always hyperbolic saddle;

(Ⅱ) (1) if $R_0 > 0$, $\theta\leq\theta^*$, $\frac{1}{p+m} > r$, then $E_2$ is a stable node;

(2) if one of the conditions hold

(i) $R_0 > 0, \; \theta > \theta^*, \; \frac{1}{p+m} > r$;

(ii) $R_0 > 0, \; \frac{1}{p+m}\leq r$,

then $E_2$ is a hyperbolic saddle;

(3) if $R_0 > 0$, $\theta = \theta^*$, $\frac{1}{p+m} > r$, then $E_2$ is a degenerate equilibrium.

Proof. The corresponding Jacobian matrix for equilibria $E = (u, v)$ of model (1.3) is easily obtained as follows:

where

We can obtain the following characteristic equation:

where

(Ⅰ) For equilibrium $E_0$ and $E_1$, it is obvious that one of the eigenvalues is $\lambda = 1 > 0$, then $E_0$ and $E_1$ are hyperbolic saddle.

(Ⅱ) For equilibrium $E_2$, we obtain

It is obvious that the sign of $\mathrm{det}(J_{E_2})$ is equivalent to the sign of $A_3$. Then, from Lemma 2.2, we know the sign of $\mathrm{det}(J_{E_2})$. Hence, when $R_0 > 0$, $\theta\leq\theta^*$, $\frac{1}{p+m} > r$, we have $\mathrm{tr}(J_{E_2}) < 0$ and $\mathrm{det}(J_{E_2}) > 0$, which imply that $E_2$ is a stable node. When $(\mathbf{S_2})$ or $(\mathbf{S_3})$ holds, we have $\mathrm{det}(J_{E_2}) < 0$, which implies that $E_2$ is a hyperbolic saddle. When $R_0 > 0$, $\theta = \theta^*$, $\frac{1}{p+m} > r$, we have $\mathrm{det}(J_{E_2}) = 0$, which implies that $E_2$ is a degenerate equilibrium. □

Theorem 2.3. Assume that $E^*$ exists, and

(1) if $\chi > \max\{\chi^*, \chi^*+1-qb_{12}\}$, then $E^*$ is stable;

(2) if $\chi < \min\{\chi^*, \chi^*+1-qb_{12}\}$, then $E^*$ is unstable;

(3) if $\chi = \chi^*$ and $qb_{12} > 1$, then Hopf bifurcation occurs.

Proof. For $E^*(u^*, v^*)$, we obtain

where

Then we obtain the characteristic equation of ${E}^*$

where

By using the Routh-Hurwitz criterion, the conditions for the stability of $E^*$ are established, that is to say, $E^*$ is stable when $\chi > \max\{\chi^*, \chi^*+1-qb_{12}\}$, while $E^*$ is unstable when $\chi < \min\{\chi^*, \chi^*+1-qb_{12}\}$.

When $\chi = \chi^*$ and $qb_{12} > 1$, we have $\mathrm{tr}(J_{E^*}) = 0$ and $\mathrm{det}(J_{E^*}) > 0$. Then, we have roots $\lambda_{1, 2} = \pm i\sqrt{\mathrm{det}(J_{E^*})}$ at $\chi = \chi^*$. Choosing $\chi$ as a bifurcation parameter (in fact, $\theta$ is the control parameter), we have

Hence, model (1.3) undergoes a Hopf bifurcation at $\chi = \chi^*$. □

Remark 2.1. For positive equilibrium $E^*$, assume that $qb_{12} > 1$, and

(1) if $\chi > \chi^*$, then $E^*$ is stable;

(2) if $\chi < \chi^*$, then $E^*$ is unstable;

(3) if $\chi = \chi^*$, then $E^*$ emerges Hopf bifurcation.

That means that $\chi^*$ can determine the stability of $E^*$.

Remark 2.2. (1) When $m = 0$ and $p = 0$, the B-D functional response of model (1.3) changes to linear approximations. In this case, the model has a unique positive equilibrium and exhibits limit cycles under fear interference.

(2) When $m = 0$, the B-D functional response of model (1.3) becomes the Holling Ⅱ functional response. In this case, the model still has a unique positive equilibrium and exhibits limit cycles under fear interference, but shows complexity compared to case (1). Readers can refer to the research results in ^[14].

Theorem 2.4. Suppose that $R_0 > 0$, $\gamma > \frac{h}{pm}$, and $\sqrt{\frac{4}{1+q}\left(\gamma-\frac{h}{pm}\right)} > \left(\frac{1}{p}+q+\frac{r q^2}{\gamma}\right)$ hold. If

then $E^*$ is globally asymptotically stable.

Proof. Denote

Then, we obtain

where $B_1 = \gamma-\frac{h}{pm}$, $B_2 = \frac{1}{1+q}$, $B_3 = {K\theta}+\frac{1}{p}+q+\frac{q^2(\theta+r)}{\gamma}$.

Since conditions $R_0 > 0$, $\gamma > \frac{h}{pm}$, and $\sqrt{\frac{4}{1+q}\left(\gamma-\frac{h}{pm}\right)} > \left(\frac{1}{p}+q+\frac{r q^2}{\gamma}\right)$ hold, we obtain $B_2 > \frac{B_3^2}{4B_1}$ when

which imply that $\frac{ \rm{d}V}{ \rm{d}t} < 0$. Then, the proof is finished. □

2.4.   Hopf bifurcation

Theorem 2.5. Assume that $R_0 > 1$, $qb_{12} > 1$, and $\chi = \chi^*$ hold, periodic solutions bifurcated from Hopf bifurcation are stable (unstable) and the direction is subcritical (supercritical) when $\Upsilon(\chi^*) < 0 (>0)$.

Proof. Let $U = u-u^*$ and $V = v-v^*$, and we can transform model (1.3) (for brevity, $u$ and $v$ stand for $U$ and $V$, respectively) to

Rewrite model (2.4) as

where $J_{E^*}$ is denoted in (2.3) and

with

Assuming that $J_{E^*}$ has two characteristic roots, it can be written as $\lambda_{1, 2}(\chi) = \varphi(\chi)\pm i\psi(\chi)$, where

Obviously, eigenvalues $\lambda_1(\chi)$ and $\lambda_2(\chi)$ are complex conjugate if $-\chi^*+\chi-1+qb_{12} > \varphi^2$. Then, the eigenvectors of $J_{E^*}$ corresponding to the eigenvalues of $\lambda(\chi) = i\psi(\chi)$ at $\chi = \chi^*$ are given by

where $\xi_1 = \frac{\chi^*-\chi+2}{2b_{12}}, \; \xi_2 = \frac{\psi(\chi)}{b_{12}}.$

Based on the transformation $(u, v)^T = \left(\begin{array}{cc} 1 & 0\\ \xi_1 & \xi_2 \end{array} \right)(x, y)^T$, model (2.5) becomes

where

Next, we calculate the 1st Lyapunov coefficient as follows:

where

owing to

Therefore, one gets the 1st Lyapunov coefficient $\Lambda = -\frac{\Upsilon(\chi^*)}{\varphi'(\chi^*)}$, which determines the stability and direction of Hopf bifurcating periodic solution. □

3.   Analysis of the spatial model (1.4)

3.1.   Stability of $E^*$

Let $0 = \mu_0 < \mu_1 < \mu_2 < \cdots < \mu_i < \cdots$ be the eigenvalues for the operator $-\Delta$ subject to the homogeneous Neumann boundary condition on $\Omega$, where $\mu_i$ has multiplicity $m_i\geq 1$ and $\mu_1$ is the smallest eigenvalue. Denote the real-valued Sobolev space $\mathbf{S} = \left\{(u, v)\in H^2(0, l\pi)\times H^2(0, l\pi):\frac{\partial u}{\partial \mathbf{n}}|_{\partial\Omega} = \frac{\partial v}{\partial \mathbf{n}}|_{\partial\Omega} = 0\right\},$ and let $\mathbf{S}_{\mathbb{C}} = \mathbf{S}\oplus \mathbf{S}\mathrm{i} = \{s_1+s_2\mathrm{i} :s_1, s_2\in \mathbf{S}\}$. Denote $\mathbb{N}_0 = \mathbb{N}\cup\{0\}$.

The linearization equation of model (1.4) at $E^*$ is

where $D = \mathrm{diag}(d_{1}, d_{2})$ and $J_{E^*}$ is denoted in (2.3).

Then, the characteristic equation of (3.1) is given by

where

When condition $\chi > \max\{\chi^*, \chi^*+1-qb_{12}\}$ holds, then we have $\mathbf{Tr}(0) < 0$ and $\mathbf{Det}(0) > 0$, which imply that $\mathbf{Tr}(k) < 0$ always holds.

Next, we analyze the sign of $\mathbf{Det}(k)$. (1) if $\frac{d_1}{d_{2}}\geq\max\{0, \chi^*-\chi+1\}$, we can obtain $\mathbf{Det}(k) > 0$;

(2) if $\frac{d_1}{d_{2}} < \chi^*-\chi+1$, we need to discuss the sign of $\Delta = (d_1-(\chi^*-\chi+1)d_{2})^2+4d_1d_2(\chi^*-\chi+1-qb_{12})$,

(ⅰ) if $\Delta < 0$, which implies that $\mathbf{Det}(k) > 0$;

(ⅱ) if $\Delta > 0$, which implies that $\mathbf{Det}(k)$ must have negative roots, then by calculating, we have $\mathbf{Det}(k) > 0$ when $-(\chi^*-\chi+1)+2qb_{12}-2\sqrt{-qb_{12}(\chi^*-\chi+1-qb_{12})} < \frac{d_1}{d_{2}} < \chi^*-\chi+1,$ and $\mathbf{Det}(k) < 0$ when $-(\chi^*-\chi+1)+2qb_{12}-2\sqrt{-qb_{12}(\chi^*-\chi+1-qb_{12})} > \frac{d_1}{d_{2}} > 0$.

Hence, the Theorem 3.1 is shown by the properties of $E^*$ in model (1.4).

Theorem 3.1. When $\chi > \max\{\chi^*, \chi^*+1-qb_{12}\}$ hold, and

(1) if

or

then $E^*$ is locally asymptotically stable;

(2) if

then $E^*$ is unstable, and the Turing instability occurs.

Remark 3.1. In Figure 3, the blue line represents the equation

and we can see that the red region is the region where Turing instability occurs, and the green region is the region where the steady state solution $E^*$ is stable. On the other hand, it also shows that when the diffusion of predator is fixed, if the diffusion of prey is smaller, it is more difficult to maintain the stability of populations.

Theorem 3.2. When $R_0 > 0$ and $\gamma > \frac{h}{pm}$ hold, if

then $E^*$ is globally asymptotically stable.

Proof. Denote a Lyapunov function

Then

which is used in the inequality $a^2+b^2\geq 2ab$.

Then, under the Neumann boundary conditions and $\frac{1}{1+q} > \frac{1}{2}\left({K\theta}+\frac{1}{p}+q+\frac{q^2(\theta+r)}{\gamma}\right)$, $\gamma > \frac{h}{pm}+\frac{1}{2}\left({K\theta}+\frac{1}{p}+q+\frac{q^2(\theta+r)}{\gamma}\right)$, we obtain that $\frac{ \rm{d}V}{ \rm{d}t}\leq 0$, which implies that the proof of Theorem 3.2 is completed. □

3.2.   Hopf bifurcation

3.2.1.   Existence

Let $X = u-u^*$ and $Y = v-v^*$, and we can transform model (1.4) (for brevity, $u$ and $v$ represent $X$ and $Y$ again, respectively) to:

Then, model (3.3) can be rewritten as

The linearized operator of model (3.3) at $(0, 0, \chi)$ is given by

where the domain of $P{(\chi)}$ is $\mathbb{S}_{\mathbb{C}}$ and

Let $k = \frac{n^2}{l^2}$, where $\frac{n^2}{l^2}$ is the eigenvalue of $-u_{xx}\rightarrow u$ and its corrsponding eigenfunction is $\phi_{n}(x) = \cos{\frac{n}{l}x}$. Then, (3.2) becomes

where

and we obtain the eigenvalues

Let $\chi_0$ be the possible Hopf bifurcation point satisfying conditions ($\mathbf{H_1}$): there exists $n\in\mathbb{N}_0$ such that

and for the unique pair of complex eigenvalues $\rho(\chi)\pm\omega(\chi)\mathrm{i}$ near the imaginary axis

where

From (3.4), when $\chi > \max\{\chi^*, \chi^*+1+qb_{12}\}$, it is obvious that $\mathbf{T}_n(\chi_0) < 0$ and $\mathbf{D}_n(\chi_0) > 0$, which imply that $(u^*, v^*)$ is stable. Thus, any potential bifurcation points shall be in $\chi < \min\{\chi^*, \chi^*+1-qb_{12}\}$.

By calculating, we get

so the transversality condition (3.6) is always held.

($\mathbf{1}$) When $n = 0$, let $\chi_0^H = \chi^*$, then we find that $\chi_0^H$ is always a Hopf bifurcation point for $l > 0$ due to $\mathbf{T}_0(\chi_0^H) = 0$, $\mathbf{T}_j(\chi_0^H) < 0$ for any $j\geq 1$, and $\mathbf{D}_k(\chi_0^H) > 0$ for any $k\in\mathbb{N}$.

($\mathbf{2}$) When $n\geq 1$, since $b_{11}(\chi^H_0) = 0$, $b_{11}'(\chi) < 0$ for $\chi_0 < \chi^*$, we have $0 < b_{11}(\chi) < b_{11}(0): = \chi^*+1$ for $\chi_0 < \chi^*$.

Denote

Then for $l_n < l < l_{n+1}$, and $1\leq j\leq n$, let $\chi^H_j$ be the root of $\chi^*-\chi = (d_1+d_2)\frac{n^2}{l^2}$, then it satisfies $\chi^H_0 > \chi^H_j > 0$. Moreover, by $b_{11}'(\chi) < 0$ for $\chi_0 < \chi^*$, we obtain

and

Now, it is demonstrated that $\mathbf{D}_n(\chi_j^H) > 0$ for $j\neq n$. For $\chi\in(0, \chi_0^H]$, we get

Obviously, $\mathbf{D}_n(\chi) = 0$ is a quadratic function of variable $\frac{n^2}{l^2}$, so we have $\mathbf{D}_n(\chi^H_j) > 0$ when

or

Theorem 3.3. Assume that $-(\chi^*-\chi+1)+2qb_{12}-2\sqrt{-qb_{12}(\chi^*-\chi+1-qb_{12})} < \frac{d_1}{d_{2}} < \chi^*-\chi+1$ or $\frac{d_1}{d_{2}}\geq\max\{0, \chi^*-\chi+1\}$ hold, and for any $l\in(l_n, l_{n+1}]$, there are $n$ points $\chi^H_j(l)$, $j\in [1, n]$ satisfying

such that Hopf bifurcation arises at each $\chi = \chi_j^H$.

3.2.2.   Direction and stability

We apply the normal form theory and center manifold theorem ^[24,25] to prove Theorems 3.4 and 3.5.

Theorem 3.4. Assume the condition of Theorem 3.3 is satisfied, the bifurcating periodic solutions of spatial homogeneous are stable (unstable), and the direction is subcritical (supercritical) when $\mathrm{Re}(c_1(\chi_0^H)) < 0 (>0)$.

Proof. Let

where $\omega_0 = \sqrt{-1+qb_{12}(\chi^H_0)}$, such that ${P}(\chi_0)\mathfrak{q} = \omega_0\mathfrak{q}\mathrm{i}$, $P^*(\chi^H_0)\mathfrak{q}^* = -\omega \mathfrak{q}^*\mathrm{i}$, $\langle \mathfrak{q}^*, \mathfrak{q}\rangle = 1$, and $\langle \mathfrak{q}^*, \overline{\mathfrak{q}}\rangle = 0$, where $\langle \mathfrak{p}, \mathfrak{q}\rangle = \int_0^{l\pi} \overline{\mathfrak{p}}^T\mathfrak{q } \rm{d}x$.

According to ^[25], we obtain that

where

with

Then, through straightforward calculations, we can get that

According to ^[25], we have

From (3.7), we have $\rho'(\chi) < 0$, then the sign of $\mathrm{Re}(c_1(\chi^H_0))$ can determine the direction and stability of Hopf bifurcation periodic solutions. □

Theorem 3.5. Assume the conditions of Theorem 3.3 are satisfied, the bifurcating periodic solutions of spatial inhomogeneous are stable (unstable), and the direction is subcritical (supercritical) when $\mathrm{Re}(c_1(\chi_j^H)) < 0 (>0)$.

Proof. Let

where $\omega_j = \sqrt{-1+qb_{12}(\chi^H_j)-2d_2\frac{j^2}{l^2}-d_2^2\frac{j^4}{l^4}}$, which satisfies ${P}(\chi_j)\mathfrak{q} = \omega_j\mathfrak{q}\mathrm{i}$, $P^*(\chi^H_j)\mathfrak{q}^* = -\omega \mathfrak{q}^*\mathrm{i}$, $\langle \mathfrak{q}^*, \mathfrak{q}\rangle = 1$, and $\langle \mathfrak{q}^*, \overline{\mathfrak{q}}\rangle = 0$.

From ^[25], we obtain that

where

with

According to ^[25], we have

where

Then, we have

Since

where

Then, we get

Then,

where

Note that

Thus

According to ^[25], we have

From (3.7), we have $\rho'(\chi) < 0$, then the direction and stability of the periodic solution of Hopf bifurcation can be determined by the sign of $\mathrm{Re}(c_1(\chi^H_j))$. □

4.   Numerical simulations

We perform a series of numerical simulation on models (1.3) and (1.4) to illustrate our results. Let initial values $(u(0), v(0)) = (2, 3)$ and parameter values in Table 2. Then, we obtain that $R_0 = 4.6 > 0$, $\theta^* = 2.1833 > \theta = 1.2$, $\chi = 0.6467 < \chi^* = 1.4034$, which imply that model (1.3) has a stable node $E_2(23, 0)$, an unstable positive equilibrium $E^*(u^*, v^*) = (4.1418, 5.9701)$, and a stable limit cycle (see Figure 4). When $\theta = 4$, we have $\theta^* = 2.1833 < \theta = 4$ and $\chi = 0.5132 > \chi^* = 0.4843$, which imply that a stable node $E_2(23, 0)$ becomes an unstable and the unstable equilibrium becomes $E^*(u^*, v^*) = (7.1403, 9.5684)$ which is stable; moreover, the limit cycle is disappears (see Figure 5).

According to Figures 4 and 5, we find that as $\theta$($1.2\rightarrow 4$) increases, the positive equilibrium $E^*(u^*, v^*)$ of model (1.3) changes from unstable to stable, and the limit cycle changes from existence to nonexistence, in addition, it also led to an increase in population size. This fully demonstrates that $\theta$ can affect the stability and the population size of model (1.3).

Note that $\theta = \frac{r_0(1-\eta)}{b}$, that is, $\theta$ and $\eta$ are negatively correlated. Therefore, if the minimum fear cost $\eta$ is used as a key parameter to explore the dynamic behavior and the trend of populations, it can be replaced by $\theta$. We find from Figures 6 and 7 that as $\theta$ increases, that is, as the minimum fear cost $\eta$ decreases, the model (1.3) exhibits a limit cycle, which first changes from a large and stable range to a small and unstable range and gradually disappears, while the positive equilibrium gradually changes from an unstable to a globally stable. This means that a lower cost of fear can maintain population stability, causing a decrease in population size but not leading to extinction, while a higher cost of fear can disrupt stability and cause periodic changes in population size.

We draw the bifurcation diagram of model (1.3) with the various parameters $\theta$, $r$, $K$, $h$, $p$, $q$, $m$ (see Figures 7–9) and observe that the various parameters can stabilize the oscillating model via Hopf bifurcation. It is revealed that: (ⅰ) for $\theta$, $r$, $m$, $p$, the small parameter values will destroy the stability and produce Hopf bifurcation; (ⅱ) for $K$, the large parameter value will destroy the stability and produce Hopf bifurcation; (ⅲ) for $h$, $q$, the bubble phenomenon ^[26] appears, which means that the appropriate parameter values will destroy the stability and produce Hopf bifurcation, but the small and large parameter values cannot destroy the stability of populations. Therefore, different parameters can have an impact on the dynamic behavior of model (1.3), mainly interfering with the stability of the equilibria and the existence, stability, and direction of Hopf bifurcations.

For model (1.4), let initial values be $(u_0(x), v_0(x)) = (4.5297+0.2\sin(x), 7.0543+0.2\sin(x))$, $d_2 = 0.2$, $\Omega = (0, 40)$, and change the different parameter values $\theta$ and $d_1$ to discuss the influence of fear and diffusion rate of prey on the dynamics of model (1.4).

Let $\theta = 4$, $d_1 = 0.1$, and according to Theorem 3.1, we get that $\chi = 0.5132 > \chi^* = 0.4843$ and $-(\chi^*-\chi+1)+2qb_{12}-2\sqrt{-qb_{12}(\chi^*-\chi+1-qb_{12})} = 0.0969 < \frac{d_1}{d_{2}} = 0.5 < \chi^*-\chi+1 = 0.9711,$ which imply that the steady state solution of model (1.4) is locally asymptotically stable. If $d_1$ is selected as 0.01 and 0.001, respectively, according to Theorem 3.1, we have

which imply that the Turing instability occurs. This further reflects that diffusion can lead to the occurrence of Turing patterns, providing ideas for studying the morphology of populations and effectively enriching the explosion of species diversity.

Next, select $\theta$ as 6 and 10, respectively, then we can also get that

By Theorem 3.1, the steady state solution of model (1.4) is locally asymptotically stable when $\theta = 6$ or $\theta = 10$.

By comparing Figures 10 and 11 in detail, we find that (ⅰ) the larger $\theta$ and $d_1$, the steady state solution of model (1.4) is more stable, indicating that the smaller the fear effect, or the larger the diffusion rate of prey (predator), the more conducive to the stability of prey (predator) population; (ⅱ) the smaller $\theta$ and $d_1$, the steady state solution of model (1.4) is more unstable, indicating that the greater the fear effect, or the smaller the diffusion rate of prey, the less conducive to the stability of prey(predator) population. Thus, the conclusion of Theorem 3.1 and Remark 3.1 is verified.

In Figure 12, we give the time evolution of model (1.4) with different $\theta$ and $d_1$ at $x = 10$, which corresponds to Figures 10 and 11. These indicated that when the solution of model (1.4) tends to the equilibria or periodic solution, its changing form is calm; but when its changing form undergoes a big sudden change, it means that Turing instability has occurred.

In Figure 13(a)–(c), the dynamics of the solution of model (1.4) is a smooth oscillatory. In Figure 13(d), the approximations have evolved into the spatially homogeneous steady states $u^*$ and $v^*$. Then, we have a conclusion that the larger $\theta$ can keep the solution of model (1.4) to the stationary distribution, while a small $\theta$ can make the solution tend to the smooth oscillatory, that is, the smaller the fear effect, the more stable the prey (predator) population will be.

5.   Conclusions

In this paper, a predator-prey model with B-D functional response for prey and a modified L-G functional response for predator are established to explore how fear and diffusion affect the spatiotemporal dynamics of the model. For the nonspatial model (1.3):

(1) The model exhibits rich dynamic properties. In terms of the existence of equilibrium points, the model has three boundary equilibria and at most three different positive equilibria (see Theorem 2.1 or Table 2).

(2) Discussing the local stability of the positive equilibria (see Theorems 2.2 and 2.3) and the global stability of $E^*$ (see Theorem 2.4), which shows that under the assumption of $qb_{12} > 1$, when $\chi > \chi^*$, $E^*$ is stable; when $\chi < \chi^*$, $E^*$ is unstable; when $\chi = \chi^*$, the Hopf bifurcation will occur.

(3) Exploring the impact of fear on Hopf bifurcation. Fear determines the direction and stability of periodic solutions bifurcated from Hopf bifurcation that are investigated (see Theorem 2.5). The lower cost of fear facilitates the occurrence of Hopf bifurcation, leading to the emergence of limit cycles (see Figures 4–7).

(4) Our results indicate that fear can enrich and destroy the dynamic properties of the model.

For the spatial model (1.4):

(1) The diffusion of populations led to the occurrence of Turing patterns. When the ratio of the diffusion rate of prey to the diffusion rate of predator exceeds a certain critical value, Turing instability occurs (see Theorem 3.1 and Figure 3), which favors the formation of biodiversity.

(2) Discussing the stability and direction of spatially homogeneous and inhomogeneous periodic solutions (see Theorems 3.4 and 3.5), which are mainly influenced by the decisive effects of fear and diffusion (see Figures 10–13).

(3) Our results reveal that the larger fear and lower diffusion rate of prey both can destroy the stability of populations and further promote the emergence of periodic solutions, but to a certain extent, and it is advantageous to improve the survival of the species and the formation of biodiversity.

(4) Our results enrich and develop the dynamics of predator-prey models with diffusion and different functional responses, such as the Holling Ⅱ term ^[5,12,24], L-G term ^[27], and B-D term ^[13,22].

In the future, we will continue to study the effects of fear and diffusion on the three population model, food web model, partial differential population models with cross diffusion, prey-taxis, and spatial memory. In addition, it can also study the influence of external environmental noise, impulsive, and delay on population dynamics.

Use of AI tools declaration

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

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Nos. 12171004, 12301616), the Doctor Research Project Foundation of Liaoning Province (No. 2023-BS-210), the Liaoning Key Laboratory of Development and Utilization for Natural Products Active Molecules (Nos. LZ202302, LZ202303), the Basic scientific research project of Liaoning Provincial Department of Education (No. JYTMS20231706), and the Doctor Research Project Foundation of Anshan Normal University (No. 23b08).

Conflict of interest

The authors declare there is no conflicts of interest.