Dynamics analysis of a predator-prey model incorporating fear effect in prey species

1.   Introduction

Predator-prey models constitute a fundamental class of ecological systems, extensively employed to characterize the complex interaction dynamics between prey and predator populations ^[1,2,3,4]. Typically, in a predator-prey system, the prey species will be consumed by the predators; in other words, the growth of the predator will dependent on the density of the prey species. Nowadays, scholars have developed diverse mathematical frameworks capturing the intrinsic behavioral traits of both populations owing to a basic logistic growth equation framework. Of particular interest in our investigation, we consider the following non-dimensional Leslie-Gower type predator-prey system with a Holling Type Ⅲ functional response:

where $w$ and $v$ represent the densities of the prey and predators, the term $\frac{a w^2}{w^2+sv^2}$ is the well-known ratio-dependent Holling Type Ⅲ functional response ^[5,6], and $\frac{v}{u}$ corresponds to the Leslie-Gower-type trophic interaction ^[7,8,9], where the parameters $a, \tau$, and $s$ are positive constants. For prior dynamic analyses of system (1.1), we direct the readers to foundational studies ^[10,11].

Fear effects, as ubiquitous non consumptive mechanisms in predator-prey interactions, significantly modulate prey species' intrinsic growth rates through predator-induced physiological stress responses. In fact, the integration of fear effects into predator-prey systems has been systematically investigated through various modeling frameworks. According to a predator-prey model with the fear factor, Tripathi et al. ^[12] obtained the stability of the system's equilibria, the existence of the Hopf bifurcation, and the Turing instability, as well as pattern formation. By incorporating fear effects and prey refuges for a Filippov prey-predator model, Hamdallah and Arafa ^[13] analyzed the stability and they obtained bifurcation sets of pseudo equilibrium and local/global sliding bifurcations. By utilizing an ecological model with fear, group defense, and the Allee effect, Kumar et al. ^[14] investigated the complex dynamical properties of the predator-prey system via saddle-node, Hopf, homoclinic, Bautin, and Bogdanov-Takens bifurcation. Cao et al. ^[15] reported the Hopf bifurcation, steady-state bifurcation, and the stability of the bifurcating periodic solution of a reaction-diffusion predator-prey system with fear effects. For comprehensive analyses of fear effect dynamics in ecological modeling, see the existing studies in ^{[16,17,18,19,20]}. Building upon these theoretical foundations, we propose the following modified system through fear effect coupling:

where $\frac{1}{1+\theta v}$ models the fear effect with the positive constant $\theta$, which describes the level of fear. It is clear that $|\frac{w}{1+\theta v}|\leq|w|$. This implies that the fear effect directly impacts the prey population's growth rate, i.e., a larger value of the fear effect control parameter $\theta$ corresponds to a lower actual growth rate of the prey population. Meanwhile, if the fear factor disappears in the system (1.2), namely, $\theta = 0$, then the system (1.2) will reduce to system (1.1).

Both systems (1.1) and (1.2) are governed by ordinary differential equations (ODEs), which exclusively describe the temporal dynamics of predator-prey interactions. These formulations capture solely temporal population dynamics without accounting for the spatial distribution mechanisms. However, natural population movements necessitate the incorporation of diffusion processes for both species within the spatial domains, thereby extending the framework to reaction-diffusion systems. Such systems exhibit enhanced dynamic complexity characterized by spatiotemporal pattern formation. For foundational studies on ecological reaction-diffusion systems, see ^[21,22,23]. Consequently, if we consider the diffusion effect in the system (1.2), then we get

where $d_1\Delta w$ and $d_2\Delta v$ describe the random diffusion progress of the prey $w$ and predator $v$ with the diffusion rates $d_1$ and $d_2$, respectively; the notation $\Delta$ is the Laplacian operator; $\Omega = (0, \pi)$ is a bounded domain with the smooth boundary $\partial\Omega$; $\mathbf{n}$ is the outward unit normal vector; and $\partial \mathbf{n}$ is the operator of the directional derivative along the direction $\mathbf{n}$. In addition, $w_0(x)$ and $v_0(x)$ represent the initial densities of the prey $w$ and predator $v$, respectively. All parameters involved in the system (1.3) are set to be positive constants.

In this current paper, we want to explore the spatiotemporal dynamics of the modified systems. To be specific, we will investigate the stability of the positive equilibrium, the Hopf bifurcation, and its natural property for the ODE system (1.2) by virtue of the first Lyapunov coefficient technique; for the reaction-diffusion system (1.3), we will examine the precise conditions of the stability of the equilibrium, the existence of the Turing instability, and the Hopf bifurcation. Especially, we will determine the stability of the bifurcating periodic solution of the reaction-diffusion system (1.3) by using normal form theory and center manifold reduction ^{[24,25,26,27]}. Collectively, the complex spatiotemporal dynamics of the reaction-diffusion system (1.3) can be presented by incorporating the fear factor in the system. The outline of this paper is as follows. In Section 2, we study the stability of the equilibrium, and the occurrence conditions of the Hopf bifurcation and its direction for the ODE system. In Section 3, we conduct a stability analysis and explore the existence of the Turing instability and the Hopf bifurcation for the reaction-diffusion system. Section 4 displays the computational validation through numerical simulations. In Section 5, we summarize this paper with some conclusions.

2.   Hopf bifurcation for the local system

2.1.   Stability analysis

To yield the equilibria of system (1.2), we set

Now, we look for the positive equilibrium, say $E_* = (w_*, v_*)$, of the system (1.2) by setting $f(w, v) = g(w, v) = 0$ as $w > 0$ and $v > 0$. Utilizing (2.1), one has

where $A: = a/(1+s) > 0$. By using (2.2), we can see that $L(w_*)\geq\theta w_*^2+(1+\theta A)w_* > 0$ as $a\geq1+s$. This implies that we must ensure that $0 < a < 1+s$ for the possibility of $L(w_*) = 0$. Now, if $0 < a < 1+s$ is valid, we have the existence criterion of the roots of (2.2) is $\Delta(w_*): = (\theta A-1)^2+4\theta > 0$. That is to say, $L(w_*) = 0$ must have two real roots. Note that $0 < a < 1+s$, so one has $0 < A < 1$. Therefore, $L(w_*) = 0$ has a unique positive root

where $A = a/(1+s)$. Hence, there is the unique positive equilibrium $E_* = (w_*, v_*) = (w_*, w_*)$ of the system (1.2) when $0 < a < 1+s$.

We have the following result concerning the system (1.2).

Theorem 2.1. Suppose that $0 < a < 1+s$ is valid.

(1) The positive equilibrium $E_*$ is locally asymptotically stable if one of the following conditions is satisfied:

(C1) $s\geq1$;

(C2) $0 < s < 1, \; \frac{a(1-s)}{(1+s)^2}\leq w_*$;

(C3) $0 < s < 1, \; \frac{a(1-s)}{(1+s)^2} > w_*, \; \tau > \frac{a(1-s)}{(1+s)^2}-w_*$.

However, it becomes unstable when

(C4) $0 < s < 1, \; \frac{a(1-s)}{(1+s)^2} > w_*, \; 0 < \tau < \frac{a(1-s)}{(1+s)^2}-w_*$.

(2) If $0 < s < 1, \; \frac{a(1-s)}{(1+s)^2} > w_*$, the system (1.2) experiences the Hopf bifurcation as $\tau = \tau_0^H$, where

Proof. At the positive equilibrium $E_* = (w_*, v_*)$, the Jacobian matrix takes the form:

We obtain the characteristic equation as follows:

where

(1) If $s\geq1$ is satisfied, then one has $T_{E_*}(\tau) < 0$ and $D_{E_*}(\tau) > 0$, and it is concluded that the positive equilibrium $E_*$ is locally asymptotically stable. (2) We can infer that if $0 < s < 1$ and $\frac{a(1-s)}{(1+s)^2}\leq w_*$, then one has $T_{E_*}(\tau) < 0$ and $D_{E_*}(\tau) > 0$. This means that the positive equilibrium $E_*$ is locally asymptotically stable. (3) If $0 < s < 1, \; \frac{a(1-s)}{(1+s)^2} > w_*$, and $\tau > \frac{a(1-s)}{(1+s)^2}-w_*$, we can see that $T_{E_*}(\tau) < 0$ and $D_{E_*}(\tau) > 0$. So, the positive equilibria $E_*$ is locally asymptotically stable. (4) If $0 < s < 1, \; \frac{a(1-s)}{(1+s)^2} > w_*$, and $0 < \tau < \frac{a(1-s)}{(1+s)^2}-w_*$, this immediately yields $T_{E_*}(\tau) > 0$, and $D_{E_*}(\tau) > 0$. So, the positive equilibrium $E_*$ is unstable. Finally, solving $\lambda$ for (2.3) yields

Therefore, if $0 < s < 1, \; \frac{a(1-s)}{(1+s)^2} > w_*$, and $\tau: = \tau_0^H = \frac{a(1-s)}{(1+s)^2}-w_*$, we get

This implies that the characteristic Eq (2.3) has a pair of purely imaginary roots. On the other hand, we can compute

Consequently, the system (1.2) presents the Hopf bifurcation as $\tau = \tau_0^H$. The proof readily follows.

2.2.   Direction of Hopf bifurcation for system (1.2)

Utilizing (2) of Theorem 2.1, we know that the system (1.2) experiences the Hopf bifurcation as $\tau = \tau_0^H$, where

Now, our goal is to investigate the direction of this Hopf bifurcation. To this end, consider the translation $\tilde{w} = w-w_*$, $\tilde{v} = v-v_*$ and still denote $\tilde{w}$ and $\tilde{v}$ as $w$ and $v$, respectively. As such, the system (1.2) becomes

Rewrite the system (2.5) as follows:

where

with

and

Define a matrix

where

and $\zeta(\tau)$ and $\eta(\tau)$ can be found in (2.4). Let

Substituting it into (2.6), we have

where

The Taylor series expansion of $\widehat{F}(\tau, \widehat{w}, \widehat{v})$ and $\widehat{G}(\tau, \widehat{w}, \widehat{v})$ shows that

where

and $k_{ij} = \frac{\zeta(\tau)-J_{11}(\tau)}{\eta(\tau)}j_{ij}-\frac{J_{12}(\tau)}{\eta(\tau)}\widetilde{k}_{ij}$ for $i, j = 0, 1, 2, \cdots$ with

When $\tau = \tau_0^H$, we define the first Lyapunov coefficient as follows:

where

and where

and $k_{ij}^* = -\frac{J_{11}(\tau_0^H)}{\eta(\tau_0^H)}j_{ij}^*-\frac{J_{12}(\tau_0^H)}{\eta(\tau_0^H)}\widetilde{k}_{ij}^*$ for $i, j = 0, 1, 2, \cdots$ with

By some direct calculations, one has

Now we summarize the analysis above and yield the following result.

Theorem 2.2. Suppose that $0 < a < 1+s, \; 0 < s < 1, \; \frac{a(1-s)}{(1+s)^2} > w_*$ are true. Then, the Hopf bifurcation is supercritical (respectively, subcritical) when $\mathit{\mbox{Re}}(\ell_1(\tau_0^H)) < 0$ (respectively, $\mathit{\mbox{Re}}(\ell_1(\tau_0^H)) > 0$) and the bifurcating periodic solution is stable (respectively, unstable).

3.   Hopf bifurcation for the diffusive system

3.1.   An estimate

In this subsection, we want to give an estimate with respect to the classic solution $(w, v)$ of the reaction-diffusion system (1.3).

Theorem 3.1. Suppose that $0 < a < \frac{2\sqrt{s}}{1+\theta}$ is fulfilled. Then the solution $(w(x, t), v(x, t))$ of the system (1.3) admits

In addition, one obtains

Proof. By employing the first equation of the system (1.3), we get

As such, the comparison principle for the parabolic equations show that $t_1\gg1$ and $0 < \varepsilon_1\ll1$ exist such that $w(x, t)\leq 1+\varepsilon_1$ for $\forall x\in\overline{\Omega}, \; t\geq t_1$. Keeping this result in mind and utilizing the second equation of the system (1.3), one derives

According to the comparison principle of parabolic equations, we get $v(x, t)\leq (1+\varepsilon_1)+\varepsilon_2$ for $\forall x\in\overline{\Omega}, \; t\geq t_2$. In the sequel, let us explore the lower-bounds of the classic solution $(w(x, t), v(x, t))$ of the system (1.3). By the way, the lower bounds of the classic solution $(w(x, t), v(x, t))$ imply that the system (1.3) is persistence. Using the first equation of the system (1.3) again, we obtain

Thereby, $\varepsilon_3 > 0$ and $t_3 > 0$ exist such that

for $\forall x\in\overline{\Omega}, \; t\geq t_3$. Finally, one can obtain

This means that

for $\forall x\in\overline{\Omega}, \; t\geq t_4$. The proof is completed.

3.2.   Turing instability

Assume that the domain $\Omega$ takes the form $\Omega = (0, \pi)$. Let $\tilde{w} = w-w_*$, $\tilde{v} = v-v_*$ and still denote $\tilde{w}$ and $\tilde{v}$ as $w$ and $v$, respectively. As such, the system (1.3) has the form

If we let $N = (w-w_*, v-v_*)^T$, at the origin, the system (3.1) can be rewritten as

where

and

where $F(\tau, w, v)$ and $G(\tau, w, v)$ have been defined in (2.6). Obviously, the local linearized system of (3.1) admits the form

Now, consider the following solution of (3.3):

where $m\in\mathbb{N}_0 = \{0, 1, 2, \cdot\cdot\cdot\}$, $\lambda$ is the growth rate of perturbation, and $a_m$ and $b_m$ are two nonzero constants. Substituting it into the linear system (3.3) yields

where $I$ is the identity matrix and

Consequently, the characteristic equation can be obtained by setting $|\lambda I-A_m| = 0$ for $m\in\mathbb{N}_0$. In this manner, we can get

where

To obtain the emergence conditions of Turing instability, we require one of Assumptions (C1)–(C3) in Theorem 2.1 is satisfied. This means that $T_m(\tau) < 0$ must hold for any $m\in\mathbb{N}_0$. Therefore, we only need to discuss the sign of $D_m(\tau)$ to explore the occurrence of Turing instability.

We have the following.

Theorem 3.2. Suppose that $0 < a < 1+s$ and one of Assumptions (C1)–(C3) in Theorem 2.1 is satisfied. The positive equilibrium $E_*$ is locally asymptotically stable if one of the following conditions is satisfied:

(C5) $s\geq1$;

(C6) $0 < s < 1, \; \frac{a(1-s)}{(1+s)^2}\leq w_*$;

(C7) $0 < s < 1, \; \frac{a(1-s)}{(1+s)^2} > w_*, \; \tau > \frac{a(1-s)}{(1+s)^2}-w_*, \; d_2\leq d_1$;

(C8) $0 < s < 1, \; \frac{a(1-s)}{(1+s)^2} > w_*, \; \tau > \frac{a(1-s)}{(1+s)^2}-w_*, \; d_2\leq \frac{\tau}{\frac{a(1-s)}{(1+s)^2}-w_*}d_1$;

and

However, $E_*$ is Turing unstable if

Proof. If $s\geq1$ is true, then we have $D_m(\tau)\geq d_1d_2m^4+w_*\tau+\frac{\theta\tau w_*}{(1+\theta w_*)^2} > 0$. So, all the eigenvalues of (3.4) possess negative real parts. This demonstrates that $E_*$ is locally asymptotically stable. Similarly, if $0 < s < 1, \; \frac{a(1-s)}{(1+s)^2}\leq w_*$, we also obtain $D_m(\tau)\geq d_1d_2m^4+w_*\tau+\frac{\theta\tau w_*}{(1+\theta w_*)^2} > 0$. Hence, the equilibrium $E_*$ is locally asymptotically stable. If (C7) is satisfied, we can deduce

Accordingly, the equilibrium $E_*$ is locally asymptotically stable when (C7) is satisfied. If (C8) is satisfied, it is easy to analyze

Henceforth, the equilibrium $E_*$ is locally asymptotically stable when (C8) is satisfied. By utilizing Condition (C9), we have

As such, we get

where "$>$" holds since

Similarly, if Condition (C10) holds, we have

Therefore, one can derive

where "$=$" holds owing to

Benefitting from (3.6), one can see that $D_m(\tau) < 0$ is satisfied, whereas we have $T_m(\tau) < 0$. Hence, we can claim that there at least one eigenvalue of (3.4) with a positive real part. Consequently, the equilibrium $E_*$ is unstable in the Turing sense. The proof is completed.

3.3.   Hopf bifurcation and the direction

For the Hopf bifurcation of the reaction-diffusion system (1.3), we can establish the following.

Theorem 3.3. Suppose that $0 < a < 1+s$ is satisfied. In this case, reaction-diffusion system (1.3) undergoes Hopf bifurcation at the threshold $\tau = \tau_m^H$ as the following condition is true:

where

Proof. Recall that

By using the definition of $\tau_0^H$, we have

Therefore, when $\tau_m^H = \tau_0^H-(d_1+d_2)m^2$ in (C11), one has $T_m(\tau_m^H) = -(d_1+d_2)m^2+\tau_0^H-\tau_m^H = 0$ for all $m\in\mathbb{N}_0$. Furthermore, if $d_2\leq \frac{\tau_m^H}{\tau_0^H}d_1$ holds, we get

Hence, we can conclude that (3.4) must have a pair of purely imaginary roots. On the other hand, by a direct computation, one yields

Consequently, the reaction-diffusion system (1.3) undergoes Hopf bifurcation when Condition (C11) is satisfied. We end the proof.

Remark 3.1. (1) If $m = 0$, the Hopf bifurcation is spatially homogeneous; however, if $m\in\mathbb{N}_0\backslash\{0\}$, the Hopf bifurcation will be spatially inhomogeneous. (2) For the onset of the two types Hopf bifurcation $\tau_m^H$ and $\tau_0^H$, we have $\tau_m^H < \tau_{m-1}^H < \cdots < \tau_2^H < \tau_1^H < \tau_0^H$.

Now our task is to determine the direction of the Hopf bifurcation for the reaction-diffusion system (1.3). For convenience, we consider a Hopf bifurcation point $\tau = \tau_0^H$. Rewrite the system (1.3) as follows:

where $W = (w, v)^T$ and

and the conjugate operator of $L(\tau)$ is given by

In addition,

When $\tau = \tau_0^H$, we have the following critical system:

where

and

Define a inner product

Let

and

By some direct computations, we can obtain $\langle L^*(\tau_0^H)W_{1}, W_{2} \rangle = \langle W_{1}, L(\tau)W_{2} \rangle, \; L(\tau_0^H)y = i \eta(\tau_0^H) y, \; L^*(\tau_0^H)y^{*} = -i \eta(\tau_0^H) y^{*}, \; \langle y^{*}, y\rangle = 1$, and $\langle y^{*}, \overline{y}\rangle = 0$. In the following, we set $\mathbf{B}_c: = \mathbf{B}\oplus i\mathbf{B} = \{x_1+ix_2|x_1, x_2\in\mathbf{B}\}$. In light of the existing literature ^[24], $\mathbf{B} = \mathbf{B}^c\oplus \mathbf{B}^s$ holds with $\mathbf{B}^c: = \{zy+\overline{zy}|z\in\mathbb{C}\}$ and $\mathbf{B}^s = \{q\in\mathbf{B}|\langle y^*, q\rangle = 0\}$, where $z = \langle y^{*}, W\rangle$ with $W = (w, v)^T$. For any $W\in\mathbf{B}$, $z\in\mathbb{C}$ and $q = (q_1, q_2)^T\in\mathbf{B}^s$ exist such that

In this manner, we obtain

Accordingly, we rewrite (3.7) as follows:

where

By employing ^[24], the nonlinear term $F(\tau_0^H, W)$ in (3.7) can be read as:

where $R(W, W)$ and $C(W, W, W)$ are symmetric multi linear forms and

where

for $U = (u_1, u_2)^T, \; V = (v_1, v_2)^T, \; Y = (y_1, y_2)^T$, and $U, \; V, \; Y$ in $H^{2}([0, \pi])\times H^{2}([0, \pi])$ with

and

By some direct but complex calculations, one obtains

So, we can obtain

where

By direct calculations, we get

This means

From ^[24], the system (3.8) has a center manifold and it could be written as

Owing to

it follows that

In this fashion, we know that the reaction-diffusion system (1.3) can be restricted to a center manifold

where

In what follows, we rewrite system (3.9) in Poincaré normal form

where $\zeta(\tau)$ and $\eta(\tau)$ satisfy

and $\delta_{j}(\tau)$ represents complex-valued coefficients. Then

Thereby, when $\tau = \tau_0^H$, we have

It then follows that

A direct computation gives

where

To summarize, we build the following.

Theorem 3.4. Suppose that $0 < a < 1+s$ is satisfied. Then the Hopf bifurcation is supercritical (respectively subcritical) if $\frac{1}{\zeta^\prime(\tau_0^H)}\mathit{\mbox{Re}}\{\delta_{1}(\tau_0^H)\} < 0\; (\mathit{\mbox{respectively}} > 0)$. Meanwhile, the bifurcating periodic solution is stable (respectively unstable) if $\mathit{\mbox{Re}}\{\delta_{1}(\tau_0^H)\} > 0\; (\mathit{\mbox{respectively}} < 0)$.

4.   Numerical simulation

In this section, we verify our previous theoretical results by using numerical computational experiments.

First, we verify the theoretical results established in Theorem 2.1. We select the following parameters in the system (1.2): $a = 0.5, \; s = 1.5, \; \theta = 0.25, \; \text{and}\ \tau = 1.0$. Under these specific values, we can see that the assumption (C1) is satisfied, yielding the positive equilibrium $E_* = (0.6586, 0.6586)$, which is locally asymptotically stable, as shown in Figure 1(a). A similar stable result of the positive equilibrium $E_*$ could be found in Figure 1(b), where we set $a = 0.5, \; s = 1.5, \; \theta = 0.25, \; \text{and}\ \tau = 2.0$ such that (C1) is fulfilled and the positive equilibrium $E_* = (0.6586, 0.6586)$ is locally asymptotically stable. In what follows, we choose the parameters $a = 0.5, \; s = 0.5, \; \theta = 0.25, \; \text{and}\ \tau = 1.0$ in the system (1.2). As a consequence, we obtain

In this fashion, Assumption (C2) is satisfied. Our numerical experiment shows that the positive equilibrium $E_* = (0.5465, 0.5465)$ is locally asymptotically stable, as demonstrated in Figure 1(c). Now, taking $a = 0.5, \; s = 0.125, \; \theta = 6.25, \; \text{and}\ \tau = 1.0$, we obtain

These mean that Assumption (C3) holds and the positive equilibrium $E_* = (0.1223, 0.1223)$ is locally asymptotically stable (refer to Figure 1(d)). These numerical results validate the theoretical predictions presented in Theorem 2.1.

To reveal the validity of Theorem 2.2, we set $a = 0.5, \; s = 0.125, \; \text{and}\ \theta = 6.25$, then obtain $\tau_0^H = 0.2234, \; E^* = (0.1223, 0.1223)$, and $\mbox{Re}(\ell_1(\tau_0^H)) = -0.1443 < 0$. Our numerical result illustrates that there is a stable periodic solution around the positive equilibrium $E^* = (0.1223, 0.1223)$; see Figure 2(a). We also set the parameters $a = 0.35, \; s = 0.14, \; \text{and}\ \theta = 7.25$. As a result, we obtain $\tau_0^H = 0.0732, \; E^* = (0.1584, 0.1584)$, and $\mbox{Re}(\ell_1(\tau_0^H)) = -0.3673 < 0$. This produces the stable bifurcating periodic solution shown in Figure 2(b).

We now validate the theoretical results presented in Theorem 3.2 through numerical experiments.

If we take the parameters $a = 1.0, \; s = 1.5, \; \theta = 0.25, \; \tau = 1.0, \; d_1 = 1.0, \; \text{and}\ d_2 = 0.5$, then Assumption (C5) is satisfied. It is found that the positive equilibrium $E_* = (0.4907, 0.4907)$ remains stable, as shown in Figure 3.

When setting the parameters $a = 1.0, \; s = 0.5, \; \theta = 0.25, \; \tau = 1.0, \; d_1 = 1.0, \; \text{and}\ d_2 = 0.5$, one has

Hence, Assumption (C6) is satisfied. It is demonstrated that the positive equilibrium $E_* = \; (0.2701, 0.2701)$ is stable, as shown in Figure 4.

We set the parameters $a = 0.5, \; s = 0.13, \; \theta = 6.5, \; \tau = 1.0, \; d_1 = 1.0, \; \text{and}\ d_2 = 1.0$. We obtain

implying that Assumption (C7) is satisfied. Our numerical experiments illustrate that the positive equilibrium $E_* = (0.1198, 0.1198)$ is stable, as depicted in Figure 5. When we set the parameters $a = 0.5, \; s = 0.13, \; \theta = 6.5, \; \text{and}\ \tau = 1.0$ in Figure 5, but choose $d_1 = 5.0$ and $d_2 = 1.0$, we have

All conditions in (C8) are satisfied, and numerical simulations demonstrate that the positive equilibrium $E_* = (0.1198, 0.1198)$ is stable; refer to Figure 6.

If we set the parameters $a = 0.5, \; s = 0.2, \; \theta = 4.25, \; \tau = 0.8, \; d_1 = 0.35, \; \text{and}\ d_2 = 5.85$, we obtain

and

That is to say, all conditions in (C9) are satisfied. The numerical result demonstrates that the equilibrium $E_* = (0.0355, 0.0355)$ is stable; refer to Figure 7.

When we set the parameters $a = 1.0, \; s = 0.2, \; \theta = 4.25, \; \tau = 0.8, \; d_1 = 0.35, \; \text{and}\ d_2 = 5.85$, we obtain

and

Accordingly, all conditions in (C10) are fulfilled. In this fashion, the equilibrium $E_* = (0.0355, 0.0355)$ becomes unstable in the Turing sense; see Figure 8. If we keep the same parameters in Figure 8 but only change the level of fear, i.e., we set $\theta = 1.85$, then we get

and

These data demonstrate the validity of (C10), leading to Turing instability; refer to Figure 9.

Figures 3–9 validate the theoretical conclusions of Theorem 3.2.

Finally, when setting the parameters $a = 1.0, \; s = 0.3, \; \theta = 3.25, \; d_1 = 2.35, \; \text{and}\ d_2 = 0.85$, one has $E_* = (0.0623, 0.0623), \; \tau_0^H = 0.3519$, and

Therefore, we have $\mbox{Re}\{\delta_{1}(\tau_0^H)\} = 1.5854e+03 > 0$. From Theorem 3.4, we have the supercritical type Hopf bifurcation, and the bifurcating periodic solution is stable. This prediction is validated by employing numerical simulation, as shown in Figure 10.

Our findings demonstrate that all theoretical results have been rigorously validated through numerical experiments.

5.   Conclusions

In this paper, we consider a Leslie-Gower type predator-prey system with a Holling Type Ⅲ functional response and fear effects subject to no-flux boundary conditions. First, we analyze the stability of the unique positive equilibrium and the existence of Hopf bifurcation by considering the ranges of two parameters $a$ and $\tau$ for an ODE system (1.2), as formalized in Theorem 2.1. Furthermore, we rigorously demonstrate the occurrence of Hopf bifurcation at the critical threshold $\tau = \tau_0^H$ and yield its direction; see Theorem 2.2. Next, we shift our focus to the reaction-diffusion system (1.3). An estimate of the classical solution and the stability results of positive equilibrium are outlined in Theorems 3.1 and 3.2, respectively. Notably, we obtain the precise occurrence of Turing instability so that we can observe the spatial pattern formation of the system. Subsequently, by designating $\tau$ as the control parameter, we conclude that spatially inhomogeneous Hopf bifurcation will emerge in the system (1.3) when $\tau = \tau_m^H$ for $m\in\mathbb{N}_0$; see Theorem 3.3. Finally, the normal formal theory of differential equations helps us characterize the supercritical and subcritical nature of the Hopf bifurcation, allowing explicit classification of the stability of bifurcating periodic solutions; see Theorem 3.4. Numerical results illustrate that the rightness of the previous theoretical analysis. The stable positive constant steady states (see Figures 1–7), the unstable positive constant steady states and pattern formations (see Figures 8 and 9), and the stable bifurcating periodic solutions (see Figure 10) are displayed. These results demonstrate that the predator-prey system with Holling Type Ⅲ functional response and fear effects exhibits rich and complex interaction dynamics.

Author contributions

Xiaoyan Zhao: Formal analysis, writing-original draft, software, numerical computation, review and editing; Liangru Yu: Correspondence, review and editing, project administration; Xue-Zhi Li: Supervision, writing-original draft, methodology, review and editing, and project administration. All authors have read and approved the final version of the manuscript for publication.

Use of Generative-AI tools declaration

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

Acknowledgments

This work was supported by the National Natural Science Foundation of China (No. 12271143).

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.