Modeling the fear effect in the predator-prey dynamics with an age structure in the predators

1.   Introduction

Interactions between the predators and the preys are diverse and complex in ecology. The predators increase the mortality rate of preys by direct predation. In existing literatures, many predator-prey models only involve direct predation for predator-prey interactions ^[1,2,3,4]. However, besides the population loss caused by direct predation, the prey will modify their behavior, psychology and physiology in response to the predation risk. This is defined as the fear effect by Cannon ^[5]. Zanette et al. ^[6] investigated the variation of song sparrows offspring reproduction when the sounds and calls of predators were broadcasted to simulate predation risk. They discovered that the fear effect alone can cause a significant reduction of song sparrows offspring reproduction. Inspired by this experimental result, Wang et al. ^[7] incorporated the fear effect into the predator-prey model. In their theoretical analysis, linear and Holling type Ⅱ functional response are chosen respectively. According to their results, the structure of the equilibria will not be affected by the fear effects, but the stability of equilibria and Hopf bifurcation are slightly different from models with no fear effects.

Soon afterward, Wang and Zou ^[8] considered a model with the stage structure of prey (juvenile prey and adult prey) and a maturation time delay. Additionally, Wang and Zou ^[9] pointed out anti-predation behaviors will not only decrease offspring reproduction of prey but also increase the difficulty of the prey being caught. Based on this assumption, they derived an anti-predation strategic predator-prey model.

Recently, the aforementioned ordinary or delayed differential equation models were extended to reaction diffusion equation ^[10,11,12]. Wang et al. ^[11] used several functional responses to study the effect of the degree of prey sensitivity to predation risk on pattern formation. Following this work, Wang et al. ^[10] introduced spatial memory delay and pregnancy delay into the model. Their numerical simulation presented the effect of some biological important variables, including the level of fear effect, memory-related diffusion, time delay induced by spatial memory and pregnancy on pattern formation. Moreover, Dai and Sun ^[13] incorporated chemotaxis and fear effect into predator-prey model, and investigated the Turing-Hopf bifurcation by selecting delay and chemotaxis coefficient as two analysis parameters.

Denote by $u_1(x, t)$ and $u_2(x, t)$, the population of the prey and adult predator at location $x$ and time $t$, respectively. We suppose juvenile predators are unable to prey on. By choosing simplest linear functional response in the model of Wang and Zou ^[7], the equation for prey population is given by

where $d_1$ is a diffusion coefficient for prey, $\mu_1$ is the mortality rate for prey, $m$ is the intraspecies competition coefficient, $p$ is predation rate, $r$ is reproduction rate for prey and $g(K, u_2)$ represents the cost of anti-predator defense induced by fear with $K$ reflecting response level. We assume $g(K, u_2)$ satisfies the following conditions.

($\mathbf{H}_1$) $g(0, u_2) = g(K, 0) = 1$, $\lim\limits_{u_2\to\infty}g(K, u_2) = \lim\limits_{K\to\infty}g(K, u_2) = 0$, $\frac{\partial g(K, u_2)}{\partial K} < 0$ and $\frac{\partial g(K, u_2)}{\partial u_2} < 0$.

Considering the maturation period of the predator, we set $b(x, t, a_1)$ be the density of the predator at age $a_1$, location $x$ and time $t$. Establish the following population model with spatial diffusion and age-structure

where $x\in\Omega$, a bounded spatial habitat with the smooth boundary $\partial\Omega$, $t, a_1 > 0$, $c$ is the conversion rate of the prey to predators, $\tau > 0$ be the maturation period of predator and age-specific functions

represent the diffusion rate and mortality rate at age $a_1$, respectively. We introduce the total population of the matured predator as $u_2 = \int_\tau^\infty b(x, t, a_1)da_1$. Thus, (1.1) together with $b(x, t, \infty) = 0$ yields

Let $s_1 = t-a_1$ and $w(x, t, s_1) = b(x, t, t-s_1)$. Along the characteristic line, solving (1.1) yields

Assume linear operator $d_0\Delta-\gamma$ with Neumann boundary conditions yields the $C_0$ semigroups $T_1(t)$. Therefore,

In particular, $\mathcal{G}(x, y, t)$ denotes the kernel function corresponding to $T_1(t)$. Thus

The above equation together with (1.2) yields a nonlocal diffusive predator-prey model with fear effect and maturation period of predators

Since the spatial movement of mature predators is much bigger than that of juvenile predators, we assume that diffusion rate of juvenile predators $d_0$ approaches zero. Hence, the kernel function becomes $\mathcal{G} = e^{-\gamma\tau}f(x-y)$ with a Dirac-delta function $f$. Thus, the equation of $u_2(x, t)$ in (1.3) becomes

It follows from the properties of Dirac-delta function that

where $B_\epsilon(x)$ is the open ball of radius $\epsilon$ centered at $x$. Therefore, model (1.3) equipped with nonnegative initial conditions and Neumann boundary conditions is

If $r\le \mu_1$, then as $t\rightarrow \infty$, we have $(u_1, u_2)\to(0, 0)$ for $x\in\Omega$, namely, two species will extinct. Throughout this paper, suppose $r > \mu_1$, which ensures that the prey and predator will persist.

This paper is organized as follows. We present results on well-posedness and uniform persistence of solutions and prove the global asymptotic stability of predator free equilibrium in Section 2. The nonexistence of nonhomogeneous steady state and steady state bifurcation are proven in Section 3. Hopf bifurcation analysis is carried out in Section 4. In Section 5, we conduct numerical exploration to illustrate some theoretical conclusions and further explore the dynamics of the nonlocal model numerically. We sum up our paper in Section 6.

2.   Preliminary

Denote by $\mathbf{C}: = C([-\tau, 0], X^2)$ the Banach space of continuous maps from $[-\tau, 0]$ to $X^2$ equipped with supremum norm, where $X = L^2(\Omega)$ is the Hilbert space of integrable function with the usual inner product. $\mathbf{C}_+$ is the nonnegative cone of $\mathbf{C}$. Let $u_1(x, t)$ and $u_2(x, t)$ be a pair of continuous function on $\Omega\times[-\tau, \infty)$ and $(u_{1t}, u_{2t})\in\mathcal{C}$ as $(u_{1t}(\vartheta), u_{2t}(\vartheta)) = (u_1(\cdot, t+\vartheta), u_2(\cdot, t+\vartheta))$ for $\vartheta\in[-\tau, 0]$. By using [14,Corollary 4], we can prove model (1.4) exists a unique solution. Note (1.4) is mixed quasi-monotone ^[15], together with comparison principle implies that the solution of (1.4) is nonnegative.

Lemma 2.1. For any initial condition $(u_{10}(x, \vartheta), u_{20}(x, \vartheta))\in \mathbf{C}_+$, model (1.4) possesses a unique solution $(u_1(x, t), u_2(x, t))$ on the maximal interval of existence $[0, t_{max})$. If $t_{max} < \infty$, then $\limsup\limits_{t\rightarrow t_{max}^{-}}(\|u_1(\cdot, t)\|+\|u_2(\cdot, t)\|) = \infty$. Moreover, $u_1$ and $u_2$ are nonnegative for all $(x, t)\in \overline\Omega\times[-\tau, t_{max})$.

We next prove $u_1$ and $u_2$ are bounded which implies that $t_{max} = \infty$.

Theorem 2.2. For any initial condition $\varphi = (u_{10}(x, \vartheta), u_{20}(x, \vartheta))\in \mathbf{C}_+$, model (1.4) possesses a global solution $(u_1(x, t), u_2(x, t))$ which is unique and nonnegative for $(x, t)\in \overline\Omega\times[0, \infty)$. If $u_{10}(x, \vartheta)\ge 0(\not\equiv0), u_{20}(x, \vartheta)\ge 0(\not\equiv0)$, then this solution remains positive for all $(x, t)\in \overline\Omega\times(0, \infty)$. Moreover, there exists a positive constant $\xi$ independent of $\varphi$ such that $\limsup\limits_{t\rightarrow \infty}u_1\le \xi, \ \limsup\limits_{t\rightarrow \infty}u_2\le \xi$ for all $x\in\Omega.$

Proof. Consider

Clearly, $w_1(x, t)$ of (2.1) is a upper solution to (1.4) due to ${\partial u_1/\partial t}\le d_1 \Delta u_1+ru_1-\mu_1 u_1-m u_1^2$. By using Lemma 2.2 in ^[16], $(r-\mu_1)/m$ of (2.1) is globally asymptotically stable in $C(\overline\Omega, \mathbb{R}^+)$. This together with comparison theorem indicates

Thus, there exists $\widetilde{\xi} > 0$ which is not dependent on initial condition, such that $\|u_1\|\le \widetilde{\xi}$ for all $t > 0$. $T_2(t)$ denotes the $C_0$ semigroups yielded by $d_2\Delta-\mu_2$ with the Neumann boundary condition. Then from (1.4),

Let $-\delta < 0$ be the principle eigenvalue of $d_2\Delta-\mu_2$ with the Neumann boundary condition. Then, $\|T_2(t)\|\le e^{-\delta t}$. The above formula yields

by choosing constants $\widetilde{B}_1\ge\|u_2(\cdot, 0)\|+ cp\tau \widetilde{\xi}\sup\limits_{a\in[-\tau, 0]}\|u_2(\cdot, a)\|$ and $\widetilde{B}_2\ge cp\widetilde{\xi}$. Using Gronwall's inequality yields $\|u_2(\cdot, t)\|\le \widetilde{B}_1 e^{\widetilde{B}_2t}$ for all $0\le t < t_{max}$. Lemma 2.1 implies $t_{max} = \infty$. So $(u_1, u_2)$ is a global solution. Moreover, if $u_{10}(x, \vartheta)\ge 0(\not\equiv0), u_{20}(x, \vartheta)\ge 0(\not\equiv0)$, then by [17,Theorem 4], this solution is positive for all $t > 0$ and $x\in\overline\Omega$.

We next prove $(u_1, u_2)$ is ultimately bounded by a constant which is not dependent on the initial condition. Due to (2.2), there exist $t_0 > 0$ and $\xi_0 > 0$ such that $u_1(x, t)\le \xi_0$ for any $t > t_0$ and $x\in\overline\Omega$. Let $z(x, t) = c u_1(x, t-\tau)+u_2(x, t)$, $\mu = \min\{\mu_1, \mu_2\}$ and $I_1 = \int_\Omega z(x, t)dx$. We integrate both sides of (1.4) and add up to obtain

Comparison principle implies

Especially, there exist $t_1 > t_0$ and $\xi_1 > 0$ such that $\|u_2\|_{1}\le \xi_1$ for all $t\ge t_1$.

Now, we define $V_{l}(t) = \int_{\Omega}(u_2(x, t))^ldx$ with $l\ge1$, and estimate the upper bound of $V_{2}(t)$. For $t > t_1$, the second equation of (1.4) and Young's inequality yield

The Gagliardo-Nirenberg inequality states:

We obtain

where $C_1 = \hat c\epsilon^{-n/2-1}2d_2\xi_1^2 > 0$, $C_2 = cp\xi_0 > 0$ and $C_3 = 2(d_2/\epsilon-C_2) > 0$ with small $\epsilon\in(0, d_2/C_2)$. Using comparison principle again yields $\limsup\limits_{t\to\infty}V_2(t)\le C_1/C_3$, which implies there exist $t_2 > t_1$ and $\xi_2 > 0$ such that $V_2\le \xi_2$ for $t\ge t_2$.

Let $L_l = \limsup\limits_{t\to\infty} V_l(t)$ with $l\ge1$, we want to estimate $L_{2l}$ with the similar method of estimation for $L_2$. Multiply the second equation in (1.4) by $2lv^{2l-1}$ and integrate on $\Omega$. Young's inequality implies

Then

via Gagliardo-Nirenberg inequality. Since $L_l = \limsup\limits_{t\to\infty}V_{l}(t)$, there exists $t_l > 0$ such that $V_{l}\le 1+L_l$ when $t > t_l$. Hence,

with $C_4 = 2cp\xi_0$. We choose $\epsilon^{-1} = (2C_4+1)l/(2d_2)$ and $C_5 = 2d_2\hat c[(2 C_4+1)/(2d_2)]^{n/2+1}$. Then for $t > t_l$, we obtain

By comparison principle, the above inequality yields $L_{2l}\le C_5 l^{n/2}(1+L_l)^2$, with a constant $C_5$ which is not dependent on $l$ and initial conditions. Finally, prove $L_{2^s} < \infty$ for all $s\in\mathbb{N}_0$. Let $B = 1+C_5$ and $\{b_s\}_{s = 0}^\infty$ be an infinite sequence denoted by $b_{s+1} = B^{(1/2)^{(s+1)}}2^{sn((1/2)^{(s+2)})}b_s$ with the first term $b_0 = L_1+1$. Clearly, $L_{2^s}\le(b_s)^{2^s}$ and

Therefore,

The above inequality leads to $\limsup\limits_{t\rightarrow \infty}u_1\le \xi$ and $\limsup\limits_{t\rightarrow \infty}u_2\le \xi$ for all $x\in\Omega$.

In Theorem 2.2, we proved that the solution of model (1.4) is uniformly bounded for any nonnegative initial condition, this implies the boundedness of the population of two species. Clearly model (1.4) exists two constant steady states $E_0 = (0, 0)$ and $E_1 = ((r-\mu_1)/m, 0)$, where $E_0$ is a saddle. Define the basic reproduction ratio ^[18] by

Thus model (1.4) possesses exactly one positive constant steady state $E_2 = (u_1^*, u_2^*)$ if and only if ${R}_0 > 1$, which is equivalent to $cp(r-\mu_1) > m\mu_2$ and $0\le\tau < \tau_{max}: = {1\over \gamma}\ln {cp(r-\mu_1)\over m\mu_2}$. Here,

The linearization of (1.4) at the positive constant steady state $(\tilde u_1, \tilde u_2)$ gives

where domain $Y: = \{(u_1, u_2)^T:u_1, u_2\in C^2(\Omega)\cap C^1(\bar\Omega), (u_1)_\nu = (u_2)_\nu = 0\ \text{on}\ \partial\Omega\}$, $W = (u_1(x, t), u_2(x, t))^T$, $\mathcal{D} = \text{diag}(d_1, d_2)$ and a bounded linear operator $\mathcal{L}:\mathbf{C}\to X^2$ is

with

The characteristic equation of (2.4) gives

or equivalently

Here, $\sigma_n$ is the eigenvalue of $-\Delta$ in $\Omega$ with Neumann boundary condition with respect to eigenfunction $\psi_n$ for all $n\in\mathbb{N}_0$, and

Theorem 2.3. (i) The trivial constant steady state $E_0 = (0, 0)$ is always unstable.

(ii) If $R_0 > 1$, then $E_1 = ((r-\mu_1)/m, 0)$ is unstable, and model (1.4) possesses a unique positive constant steady state $E_2 = (u_1^*, u_2^*)$.

(iii) If $R_0\le1$, then $E_1$ is globally asymptotically stable in $\mathbf{C}_+$.

Proof. (ⅰ) Note, (2.5) at $E_0$ takes form as $(\rho+\sigma_n d_1-r+\mu_1)(\rho+\sigma_n d_2+\mu_2) = 0$ for all $n\in\mathbb{N}_0$. Then $r-\mu_1 > 0$ is a positive real eigenvalue, namely, $E_0$ is always unstable.

(ⅱ) The characteristic equation at $E_1$ gives

Note that one eigenvalue $\rho_1 = -\sigma_n d_1-r+\mu_1$ remains negative. Hence, we only need to consider the root distribution of the following equation

According to [19,Lemma 2.1], we obtain that (2.8) exists an eigenvalue $\rho > 0$ when $R_0 > 1$, namely, $E_1$ is unstable when $R_0 > 1$.

(ⅲ) By using [19,Lemma 2.1] again, any eigenvalue $\rho$ of (2.8) satisfies $Re(\rho) < 0$ when $R_0 < 1$, namely, $E_1$ is locally asymptotically stable when $R_0 < 1$. Now, consider $R_0 = 1$. $0$ is an eigenvalue of (2.7) for $n = 0$ and all other eigenvalues satisfy $Re(\rho) < 0$. To prove the stability of $E_1$, we shall calculate the normal forms of (1.4) by the algorithm introduced in ^[20]. Set

Obviously, $\Upsilon = \{0\}$ when $R_0 = 1$. System (1.4) satisfies the non-resonance condition relative to $\Upsilon$. Denote $\overline{u}_1 = (r-\mu_1)/m$, and let $\mathbf{w} = ({w}_1, {w}_2)^T = (\overline{u}_1-u_1, u_2)^T$ and (1.4) can be written as

Here linear operator $\mathcal{A}_0$ is given by $(\mathcal{A}_0 \varphi)(\vartheta) = (\varphi(\vartheta))'$ when $\vartheta\in[-\tau, 0)$ and

and the nonlinear operator $\mathcal{F}_0$ satisfies $[\mathcal{F}_0(\varphi)](\vartheta) = 0$ for $-\tau\le\vartheta < 0$. By Taylor expansion, $[\mathcal{F}_0(\varphi)](0)$ can be written as

Define a bilinear form

Select $\alpha = (1, m/(p-r g'_{u_2}(K, 0))$ and $\beta = (0, 1)^T$ to be the right and left eigenfunction of $\mathcal{A}_0$ relative to eigenvalue $0$, respectively. Decompose $\mathbf{w}_t$ as $\mathbf{w}_t = h\alpha+\delta$ and $\langle\beta, \delta\rangle = 0$. Notice $\mathcal{A}_0\alpha = 0$ and $\langle\beta, \mathcal{A}_0\delta\rangle = 0$. Thus,

Moreover,

It follows from the above two equations that

When the initial value is a small perturbation of $E_1$, then $\delta = O(h^2)$, together with Taylor expansion yields

Therefore, we obtain the norm form of (1.4) as follows

Then for any positive initial value, the stability of zero solution of (2.10) implies $E_1$ is locally asymptotically stable when $R_0 = 1$.

Next, it suffices to show the global attractivity of $E_1$ in $\mathbf{C}_+$ when $R_0\le1$. Establish a Lyapunov functional $V:\mathbf{C}_+\to\mathbb{R}$ as

Along solutions of (1.4), taking derivative of $V(\phi_1, \phi_2)$ with respect to $t$ yields

Note $\{E_1\}$ is the maximal invariant subset of $dV/dt = 0$, together with LaSalle-Lyapunov invariance principle ^[21,22] implies $E_1$ is globally asymptotically stable if $R_0\le1$.

In Theorem 2.3, $E_0$ is always unstable which suggests that at leat one species will persist eventually. Moreover, if $R_0\le1$, then $E_1$ is globally asymptotically stable in $\mathbf{C}_+$, which implies that when the basic reproduction ratio is no more than one, the predator species will extinct and only the prey species can persist eventually. Next, we will prove the solution is uniformly persistent. $\Theta_t$ denotes the solution semiflow of (1.4) mapping $\mathbf{C}_+$ to $\mathbf{C}_+$; namely, $\Theta_t\varphi: = (u_1(\cdot, t+\cdot), u_2(\cdot, t+\cdot))\in \mathbf{C}_+$. Set $\zeta^+(\varphi) = \cup_{t\ge0}\{\Theta_t\varphi\}$ be the positive orbit and $\varpi(\varphi)$ be the omega limit set of $\zeta^+(\varphi)$. Denote

$\Gamma_{\partial}$ as the largest positively invariant set in $\partial Z$, and $W^s((\tilde u_1, \tilde u_2))$ as the stable manifold associated with $(\tilde u_1, \tilde u_2)$. We next present persistence result of model (1.4).

Theorem 2.4. Suppose $R_0 > 1$. Then there exists $\kappa > 0$ such that $\liminf\limits_{t\to\infty} u_1(x, t)\ge\kappa$ and $\liminf\limits_{t\to\infty} u_2(x, t)\ge\kappa$ for any initial condition $\varphi\in Z$ and $x\in\overline\Omega$.

Proof. Note $\Theta_t$ is compact, and Theorem 2.2 implies $\Theta_t$ is point dissipative. Then $\Theta_t$ possesses a nonempty global attractor in $\mathbf{C}_+$ ^[23]. Clearly, $\Gamma_{\partial} = \{(\varphi_1, \varphi_2)\in\mathbf{C}_+: \varphi_2\equiv0\}$, and $\varpi(\varphi) = \{E_0, E_1\}$ for all $\varphi\in \Gamma_{\partial}$. Define a generalized distance function $\psi$ mapping $\mathbf{C}_+$ to $\mathbb{R}^+$ by

Following from strong maximum principle ^[24], $\psi(\Theta_t\varphi) > 0$ for all $\varphi\in Z$. Due to $\psi^{-1}(0, \infty)\subset Z$, assumption $(P)$ in [25,Section 3] holds. Then verify rest conditions in [25,Theorem 3].

First, prove $W^s(E_0)\cap \psi^{-1}(0, \infty) = \emptyset$. Otherwise, there exists an initial condition $\varphi\in \mathbf{C}_+$ with $\psi(\varphi) > 0$, such that $(u_1, u_2)\to E_0$ as $t\to\infty$. Thus, for any sufficiently small $\varepsilon_1 > 0$ satisfying $r g(K, \varepsilon_1)-\mu_1 > p\varepsilon_1$, there exists $t_1 > 0$ such that $0 < u_1, u_2 < \varepsilon_1$ for all $x\in \Omega$ and $t > t_1$. Note that $rg(K, 0)-\mu_1 > 0$ and $\partial g(K, u_2)/\partial u_2 < 0$ ensure the existence of small $\varepsilon_1 > 0$. Then the first equation in (1.4) and $(\mathbf{H}_1)$ lead to

Notice

has a globally asymptotically stable positive steady state $(r g(K, \epsilon_1)-\mu_1-p\epsilon_1)/m$ due to [16,Lemma 2.2], together with comparison principle yields $\lim\limits_{t\to\infty}u_1\ge \lim\limits_{t\to\infty}\hat u_1 > 0$. A contradiction is derived, so $W^s(E_0)\cap \psi^{-1}(0, \infty) = \emptyset$.

Next check $W^s(E_1)\cap\psi^{-1}(0, \infty) = \emptyset$. If not, there exists $\varphi\in\mathbf{C}_+$ with $\psi(\varphi) > 0$ such that $(u_1, u_2)$ converges to $E_1$ as $t\to\infty$. According to (2.2), for any small $\varepsilon_2 > 0$ satisfying $cpe^{-\gamma\tau}((r-\mu_1)/m-\epsilon_2) > \mu_2$, there exists $t_2 > 0$ such that $u_1 > (r-\mu_1)/m-\varepsilon_2$ for all $x\in\overline\Omega$ and $t > t_2-\tau$. Note that $R_0 > 1$ ensures the existence of $\varepsilon_2 > 0$. Thus, the second equation of (1.4) yields

In a similar manner, we derive $\lim\limits_{t\to\infty}u_2(x, t) > 0$ by above inequality, $cpe^{-\gamma\tau}((r-\mu_1)/m-\epsilon_2) > \mu_2$ and comparison principle. A contradiction yields again. Hence, it follows from Theorem 3 in ^[25] that, for any $\varphi\in \mathbf{C}_+$, there exists $\kappa > 0$ such that $\liminf\limits_{t\to\infty}\psi(\Theta_t \varphi)\ge\kappa$ uniformly for any $x\in\overline\Omega$.

Theorem 2.4 implies that when the basic reproduction ratio is bigger than one, both the predator species and prey species will persist eventually. We next investigate the stability of $E_2$. The corresponding characteristic equation at $E_2$ gives

with

Characteristic equation (2.11) with $\tau = 0$ is

We observe that $a_{0, n}+b_{0, n} = \sigma_nd_2(\sigma_nd_1+mu_1^*)-\mu_2 u_2^*(rg'_{u_2}(K, u_2^*)-p) > 0$, and $a_{1, n}+b_{1, n} = \sigma_n(d_1+d_2)+mu_1^* > 0$ for all integer $n\ge0$, which yields that any eigenvalue $\rho$ of (2.12) satisfies $Re(\rho) < 0$. Then, local asymptotic stability of $E_2$ is derived when $\tau = 0$ which implies Turing instability can not happen for the non-delay system of (1.4). In addition, $a_{0, n}+b_{0, n} > 0$ for any $n\in\mathbb{N}_0$ leads to that (2.11) can not have an eigenvalue $0$ for any $\tau\ge0$. This suggests we look for the existence of simple $\rho = \pm i\delta$ ($\delta > 0$) for some $\tau > 0$. Substitute $\rho = i\delta$ into (2.11) and then

with

Thus, $a_{0, n}-b_{0, n}\ge 0$ for all $n\in\mathbb{N}_0$ is equivalent to

If $(\mathbf{A}_0)$ holds, then (2.13) admits no positive roots, together with for $\tau = 0$, any eigenvalues $\rho$ of (2.11) satisfies $Re(\rho) < 0$, yields the next conclusion.

Theorem 2.5. Suppose $R_0 > 1$. Then, $E_2$ is locally asymptotically stable provided that $(\mathbf{A}_0)$ holds.

3.   Nonhomogeneous steady states analyses

Now, we consider positive nonhomogeneous steady states. The steady state $(u_1(x), u_2(x))$ of (1.4) satisfies the elliptic equation

From Theorem 2.3, all the solutions converge to $E_1$ when $R_0\le1$ and the positive nonhomogeneous steady state may exist only if $R_0 > 1$. Throughout this section, we assume that $R_0 > 1$. In what follows, the positive lower and upper bounds independent of steady states for all positive solutions to (3.1) are derived.

Theorem 3.1. Assume that $R_0 > 1$. Then any nonnegative steady state of (3.1) other than $(0, 0)$, and $((r-\mu_1)/m, 0)$ should be positive. Moreover, there exist constants $\overline{\mathcal{B}}, \underline{\mathcal{B}} > 0$ which depend on all parameters of (3.1) and $\Omega$, such that $\underline{\mathcal{B}}\le u_1(x), u_2(x)\le\overline{\mathcal{B}}$ for any positive solution of (3.1) and $x\in\overline\Omega$.

Proof. We first show any nonnegative solution $(u_1, u_2)$ other than $E_0$ and $E_1$, should be $u_1 > 0$ and $u_2 > 0$ for all $x\in\overline\Omega$. To see this, suppose $u_2(x_0) = 0$ for some $x_0\in\overline\Omega$, then $u_2(x)\equiv0$ via strong maximum principle and

Thus the above inequality implies $u_1(x)\equiv0$ or $u_1(x)\equiv (r-\mu_1)/m$. Now, we assume $u_2 > 0$ for all $x\in\overline \Omega$. Strong maximum principle yields $u_1 > 0$ for all $x\in\overline \Omega$. Hence, $u_1 > 0$ and $u_2 > 0$ for all $x\in\overline\Omega$.

We now prove $u_1$ and $u_2$ have a upper bound which is a positive constant. Since $-d_1\Delta u_1(x)\le (r-\mu_1-mu_1(x))u_1(x)$, we then obtain from Lemma 2.3 in ^[26] that $u_1(x)\le (r-\mu_1)/m$ for any $x\in\overline\Omega$.

By two equations in (3.1), we obtain

By using [26,Lemma 2.3] again, we conclude

Next, we only need to prove $\|u_1(x)\|$ and $\|u_2(x)\|$ have a positive lower bound which is not dependent on the solution. Otherwise, there exists a positive steady states sequence $(u_{1, n}(x), u_{2, n}(x))$ such that either $\lim\limits_{n\to\infty}\|u_{1, n}\|_{\infty} = 0$ or $\lim\limits_{n\to\infty}\|u_{2, n}\|_{\infty} = 0$. Integrating second equation of (3.1) gives

If $\|u_{1, n}(x)\|_\infty\to0$ as $n\to\infty$, then $cpe^{-\gamma\tau}u_{1, n}(x)-\mu_2 < -\mu_2/2$ for sufficiently large $n$, which yields $u_{2, n}(cpe^{-\gamma\tau}u_{1, n}-\mu_2) < 0$, a contradiction derived. Thus, $\|u_{2, n}(x)\|_{\infty}\to0$ as $n\to\infty$ holds. We then assume that $(u_{1, n}, u_{2, n})\to(u_{1, \infty}, 0)$, as $n\to\infty$ where $u_{1, \infty}\ge0$. Similarly, we obtain that either $u_{1, \infty}\equiv0$ or $u_{1, \infty}\equiv (r-\mu_1)/m$. Obviously, $u_{1, \infty}\not\equiv0$ based on the above argument, thus $u_{1, \infty}\equiv (r-\mu_1)/m$ and $\lim\limits_{n\to \infty}cpe^{-\gamma\tau}u_{1, n}(x)-\mu_2 = \mu_2(R_0-1) > 0$. This again contradicts (3.2). Hence, we have shown $\|u_1(x)\|_\infty$ and $\|u_2(x)\|_\infty$ have a positive constant lower bound independent on the solution. Therefore, $u_1(x)$ and $u_2(x)$ have a uniform positive constant lower bound independent on the solution of (3.1) via Harnack's inequality [26,Lemma 2.2]. This ends the proof.

Theorem 3.2. Suppose $R_0 > 1$. There exists a constant $\chi > 0$ depending on $r, \mu_1, p, c, \gamma, \tau, \mu_2, g$ and $\sigma_1$, such that if $\min\{d_1, d_2\} > \chi$ then model (1.4) admits no positive spatially nonhomogeneous steady states, where $\sigma_1$ is defined in (2.6).

Proof. Denote the averages of the positive solution $(u_1, u_2)$ of system (3.1) on $\Omega$ by

By (2.2), we have $u_1(x)\le \overline u_1$, where $\overline u_1 = (r-\mu_1)/m$, which implies $\widetilde{u_1}\le \overline u_1$. Multiplying the first equation by $ce^{-\gamma\tau}$ and adding two equations of (3.1) lead to

The integration of both sides for the above equation yields

It is readily seen that $\int_{\Omega}(u_1-\widetilde{u_1})dx = \int_{\Omega}(u_2-\widetilde{u_2})dx = 0$. Note $u_1$ and $u_2$ are bounded by two constants $\overline u_1 > 0$ and $\overline u_2: = rc \overline u_1/\min\{\mu_1d_2/d_1, \mu_2\} > 0$ by Theorem 3.1. Denote $M_f = \max\limits_{u_2\in[0, \overline u_2]}|g'_{u_2}(K, u_2)|$ and we then obtain

Set $A_1 = (r-\mu_1)+\left((r M_f+p)\overline u_1+cpe^{-\gamma\tau} M_v\right)/2,$ and $A_2 = \left((r M_f+p)\overline u_1+cpe^{-\gamma\tau}(2\overline u_1+M_v)\right)/2.$ Then, the above inequalities and Poincar$\acute{e}$ inequality yield that

with a positive constant $\chi = \max\{A_1/\sigma_1, A_2/\sigma_1\}$ depending on $r, f, \mu_1, p, c, \gamma, \tau, \mu_2$ and $\sigma_1$. Hence, if $\chi < \min\{d_1, d_2\}$, then $\nabla(u_1-\widetilde{u_1}) = \nabla(u_2-\widetilde{u_2}) = 0$, which implies $(u_1, u_2)$ is a constant solution.

Select $u_1^*: = \nu$ as the bifurcation parameter and study nonhomogeneous steady state bifurcating from $E^*$. Let $u_{2, \nu}$ satisfy $rg(K, u_2)-pu_2 = \mu_1+m\nu$, $E_2 = (\nu, u_{2, \nu})$, and $(\widehat u_1, \widehat u_2) = (u_1-\nu, u_2-u_{2, \nu})$. Drop $\widehat\cdot$. System (3.1) becomes

for $(\nu, u_1, u_2)\in\mathbb{R}^+\times \mathcal{Y}$ with $\mathcal{Y} = \{(u_1, u_2): u_1, u_2\in H^2(\Omega), (u_1)_{\nu} = (u_2)_{\nu} = 0, \ \text{on}\ \partial\Omega\}$. Calculating Fr$\acute{e}$chet derivative of $\mathcal{H}$ gives

Then the characteristic equation follows

where

Obviously, $Q_i > 0$ and $P_i > 0$ for all $\nu\in\mathbb{R}^+$ and $i\in\mathbb{N}_0$. Therefore, (3.3) does not have a simple zero eigenvalue. According to ^[4], we obtain the nonexistence of steady state bifurcation bifurcating at $E_2$.

Theorem 3.3. Model (1.4) admits no positive nonhomogeneous steady states bifurcating from $E_2$.

4.   Hopf bifurcation analyses

Next, the stability switches at $E_2$ and existence of periodic solutions of (1.4) bifurcating from $E_2$ are studied. Suppose $R_0 > 1$, namely, $cp(r-\mu_1) > m\mu_2$ and $0\le\tau < \tau_{max}: = {1\over \gamma}\ln {cp(r-\mu_1)\over m\mu_2}$ to guarantee the existence of $E_2$.

4.1.   Local Hopf bifurcation analyses

Recall the stability of $E_2$ for $\tau = 0$ is proved and $0$ is not the root of (2.11) for $\tau\ge0$. So, we only consider eigenvalues cross the imaginary axis to the right which corresponds to the stability changes of $E_2$. Now, we shall consider the positive root of $G_n(\delta, \tau)$. Clearly, there exists exactly one positive root of $G_n(\delta, \tau) = 0$ if and only if $a_{0, n} < b_{0, n}$ for $n\in\mathbb{N}_0$. More specifically,

is the sufficient and necessary condition to ensure $G_0(\delta, \tau)$ has exactly one positive zero. For some integer $n\ge1$, the assumption $(\mathbf{A}_1)$ is a necessary condition to guarantee $G_n(\delta, \tau)$ exists positive zeros. Set

Implicit function theorem implies $G_n(\delta, \tau)$ has a unique zero

which is a $C^1$ function for $\tau\in J_n$. Hence, $i\delta_n(\tau)$ is an eigenvalue of (2.11), and $\delta_n(\tau)$ satisfies

for $n\in\mathbb{N}_0$. Let

which is the unique solution of $\sin\vartheta_n = h_{1, n}$ and $\cos\vartheta_n = h_{2, n}$ and satisfies $\vartheta_n(\tau)\in(0, 2\pi]$ for $\tau\in I_n$.

According to ^[3,27], we arrive at the next properties.

Lemma 4.1. Suppose that $R_0 > 1$ and $(\mathbf{A}_1)$ holds.

(i) There exists a nonnegative integer $M_1$ such that $J_{n}\ne\emptyset$ for $0\le n\le M_1$, with $J_{M_1}\subset J_{M_1-1}\subset\cdots\subset J_1\subset J_0$, and $J_n = \emptyset$ for $n\ge 1+M_1$, where $J_n$ is defined in (4.1).

(ii) Define

Then, $\mathcal{S}_0^0(0) < 0$; for $0\le n\le M_1$ and $k\in\mathbb{N}_0$, we have $\lim\limits_{\tau\rightarrow \widehat\tau_n^-}\mathcal{S}_n^k(\tau) = -(2k+1)\pi$, where $\widehat\tau_n = \sup J_n$; $\mathcal{S}_n^{k+1}(\tau) < \mathcal{S}_n^{k}(\tau)$ and $\mathcal{S}_n^k(\tau) > \mathcal{S}_{n+1}^{k}(\tau)$.

(iii) For each integer $n\in[0, N_1]$ and some $k\in\mathbb{N}_0$, $\mathcal{S}_n^k(\tau)$ has one positive zero $\overline\tau_n\in J_n$ if and only if (2.11) has a pair of eigenvalues $\pm i\delta_n(\overline\tau_n)$. Moreover,

When $(\mathcal{S}_n^k)'(\overline\tau_n) < 0$, $\pm i\delta_n(\overline\tau_n)$ cross the imaginary axis from right to left at $\tau = \overline\tau_n$; when $(\mathcal{S}_n^k)'(\overline\tau_n) > 0$, $\pm i\delta_n(\overline\tau_n)$ cross the imaginary axis from left to right.

If $\sup\limits_{\tau\in I_0}\mathcal{S}_0^0\le0$, then $\mathcal{S}_n^k < 0$ in $J_n$ holds for any $0\le n\le M_1$ and $k\in\mathbb{N}_0$; or only $\mathcal{S}_0^0$ has a zero with even multiplicity in $J_0$ and $\mathcal{S}_n^k < 0$ for any positive integers $n$ and $k$. Therefore, $E_2$ is locally asymptotically stable for $\tau\in[0, \tau_{max})$. The following assumption ensures Hopf bifurcation may occur at $E_2$.

$(\mathbf{A}_2)$ $\sup\limits_{\tau\in J_0}\mathcal{S}_0^0(\tau) > 0$ and $\mathcal{S}_n^k(\tau)$ has at most two zeros (counting multiplicity) for integer $0\le n\le M_1$ and $k\in \mathbb{N}_0$.

Note $\sup\limits_{\tau\in J_n} \mathcal{S}_n^0$ is strictly decreasing in $n$ due to Lemma 4.1. It then follows from $(\mathbf{A}_2)$, and $\mathcal{S}_n^k(0) < \mathcal{S}_0^0(0) < 0$, $\lim\limits_{\tau\rightarrow \widehat\tau_n^-}\mathcal{S}_n^k(\tau) < 0$ for any integer $0\le n\le M_1$ and $k\in\mathbb{N}_0$, that we can find two positive integers

and

Then $\mathcal{S}_n^k(\tau)$ admits two simple zeros $\tau_n^k$ and $\tau_n^{2K_n-k-1}$ for $k\in[0, K_n-1]$ and no zeros for $k\ge K_n$. The above analysis, together with Lemma 4.1(ⅲ), yields the next result.

Lemma 4.2. Suppose $R_0 > 1$ and $(\mathbf{A}_1)$ and $(\mathbf{A}_2)$ hold. Let $\mathcal{S}_n^k(\tau)$, $M$ and $K_n$ be defined in (4.3), (4.5) and (4.6).

(i) For integer $n\in[0, M]$, there are $2K_n$ simple zeros $\tau_n^j \ (0\le j\le 2K_n-1)$ of $\mathcal{S}_n^i(\tau) \ (0\le i\le K_n-1)$, $0 < \tau_n^0 < \tau_n^1 < \tau_n^2 < \cdots < \tau_n^{2K_n-1} < \widehat\tau_n$, and $d\mathcal{S}_n^i(\tau_n^i)/d\tau > 0$ and $d\mathcal{S}_n^{i}(\tau_n^{2K_n-i-1})/d\tau < 0$ for each $0\le i\le K_n-1$.

(ii) If there exist exactly two bifurcation values $\tau_{n_1}^j = \tau_{n_2}^i: = \tau_*$ with $n_1\ne n_2$ and $(n_1, j), (n_2, i)\in[0, M]\times[0, 2K_n-1]$, then the double Hopf bifurcation occurs at $E_2$ when $\tau = \tau_*$.

Collect all values $\tau_n^j$ with $0\le n\le M$ and $0\le j\le 2K_n-1$ in a set. To ensure Hopf bifurcation occurs, remove values which appear more than once. The new set becomes

Lemma 4.1(ⅱ) implies $\mathcal{S}_0^0(\tau)$ exists two simple zeros $\tau_0 < \tau_{2L-1}$. When $\tau = \tau_i$ with $0\le i\le 2L-1$, the Hopf bifurcation occurs at $E_2$. Moreover, $E_2$ is locally asymptotically stable for $\tau\in[0, \tau_0)\cup(\tau_{2L-1}, \tau_{max})$ and unstable for $\tau\in(\tau_0, \tau_{2L-1})$. Define

Theorem 4.3. Suppose $R_0 > 1$. Let $J_n$, $\mathcal{S}_n^k(\tau)$, $\Sigma$ and $\Sigma_0$ be defined in (4.1), (4.3), (4.7) and (4.8), respectively.

(i) $E_2$ is locally asymptotically stable for all $\tau\in[0, \tau_{max})$ provided that either $J_0 = \emptyset$ or $\sup\limits_{\tau\in J_0} \mathcal{S}_0^0(\tau)\le0$.

(ii) If $(\mathbf{A}_1)$ and $(\mathbf{A}_2)$ hold, then a Hopf bifurcation occurs at $E_2$ when $\tau\in \Sigma$. $E_2$ is locally asymptotically stable for $\tau\in[0, \tau_0)\cup(\tau_{2L-1}, \tau_{max})$, and unstable for $\tau\in(\tau_0, \tau_{2L-1})$. Further, for $\tau\in\Sigma\backslash \Sigma_0$, the bifurcating periodic solution is spatially nonhomogeneous; for $\tau\in\Sigma_0$, the bifurcating periodic solution is spatially homogeneous.

4.2.   Global Hopf bifurcation analyses

Next investigate the properties of bifurcating periodic solutions by global Hopf bifurcation theorem ^[28]. Set $y(t) = (y_1(t), y_2(t))^T = (u_1(\cdot, \tau t)-u_1^*, u_2(\cdot, \tau t)-u_2^*)^T$ and write (1.4) as

where $y_t(\nu) = y(t+\nu)$ for $\nu\in[-1, 0]$, $A = \text{diag}(\tau d_1\Delta-\tau \mu_1, \tau d_2\Delta-\tau \mu_2)$ and

$\{\Psi(t)\}_{t\ge0}$ denotes the semigroup yielded by $A$ in $\Omega$, with Neumann boundary condition. Clearly, $\lim\limits_{t\to\infty}\Psi(t) = 0$. The solution of (4.9) can be denoted by

If $y(t)$ is a $a-$periodic solution of (4.9), then (4.10) yields

since $\Psi(t+na)y(0)\to0$ as $n\to\infty$. Hence we only need to consider (4.11). The integral operator of (4.11) is differentiable, completely continuous, and G-equivariant by [24,Chapter 6.5]. The condition $\min\{d_1, d_2\} > \chi$ ensures (1.4) admits exactly one positive steady state $E_2$. Using a similar argument as in [29,Section 4.2], $(H1)$–$(H4)$ in [24,Chapter 6.5] hold and we shall study the periodic solution.

Lemma 4.4. Suppose $R_0 > 1$, then all nonnegative periodic solutions of (4.9) satisfies $\kappa\le u_1(x, t), u_2(x, t) \le \xi$ for all $(x, t)\in\overline\Omega\times\mathbb{R}^+$, where $\xi$ and $\kappa$ are defined in Theorems 2.2 and 2.4, respectively.

Lemma 4.4 can be obtained by Theorems 2.2 and 2.4. We further assume

This technical condition is used to exclude the $\tau-$periodic solutions. Note assumption $(\mathbf{A}_3)$ holds when $K = 0$. This, together with that $(g(K, u_2)-g(K, u_2^*))/(u_1-u_1^*)$ is continuous in $K$, implies there exists $\varepsilon > 0$ such that $(\mathbf{A}_3)$ holds for $0\le K < \varepsilon$, that is, $(\mathbf{A}_3)$ holds for the model (1.4) with weak fear effect.

Lemma 4.5. Assume that $R_0 > 1$ and $(\mathbf{A}_3)$ holds, then model (1.4) admits no nontrivial $\tau-$periodic solution.

Proof. Otherwise, let $(u_1, u_2)$ be the nontrivial $\tau-$periodic solution, that is, $(u_1(x, t-\tau), u_2(x, t-\tau)) = (u_1(x, t), u_2(x, t))$. Thus, we have

Claim

To see this, establish the Lyapunov functional $\mathbb{L}_1: C(\bar\Omega, \mathbb{R}^+\times\mathbb{R}^+)\to\mathbb{R}$,

Along the solution of system (4.12), the time derivative of $\mathbb{L}_1(\phi_1, \phi_2)$ is

The assumption ($\mathbf{A}_3$) ensures $d\mathbb{L}_1/dt\le0$ for all $(u_1, u_2)\in C(\bar\Omega, \mathbb{R}^+\times\mathbb{R}^+)$. The maximal invariant subset of $d\mathbb{L}_1/dt = 0$ is $\{E_2\}$. Therefore, $E_2$ attracts all positive solution of (4.12) by LaSalle-Lyapunov invariance principle ^[21,22] which excludes the nonnegative nontrivial $\tau-$periodic solution.

To obtain the nonexistence of $\tau-$periodic solution for model (1.4), we must use the condition $(\mathbf{A}_3)$ which is very restrictive. However, in numerical simulations, Lemma 4.5 remains true even $(\mathbf{A}_3)$ is violated. Thus, we conjecture the nonexistence of $\tau$-periodic solution for (1.4).

In the beginning of this section, when $\tau = \tau_i$ with $0\le i\le 2L-1$, $\pm i \delta_n(\tau_i)$ are a pair of eigenvalues of (2.11). Give the next standard notations:

$(i)$ For $0\le i\le 2L-1$, $(E_2, \tau_i, 2\pi/(\delta_n(\tau_i)\tau_i))$ is an isolated singular point.

$(ii)$ $\Gamma = \mathcal{C}l\{(y, \tau, Q)\in X^2\times\mathbb{R}^+\times\mathbb{R}^+:\ y\ \text{is the nontrivial Q-periodic solution of (4.9)}\}$ is a closed subset of $X^2\times\mathbb{R}^+\times\mathbb{R}^+$.

$(iii)$ For $0\le i\le 2L-1$, $\mathcal{C}_i(E_2, \tau_i, Q_i)$ is the connected component of $(E_2, \tau_i, Q_i)$ in $\Gamma$.

$(iv)$ For integer $0\le k\le \max\limits_{n\in[0, M]} K_n-1$, let $\Sigma_H^k = \{\tau\in\Sigma: \mathcal{S}_n^k(\tau) = 0\ \text{for integer}\ 0\le n\le M\}.$

We are ready to present a conclusion on the global Hopf branches by a similar manner in [30,Theorem 4.12].

Theorem 4.6. Assume $R_0 > 1$, $\min\{d_1, d_2\} > \chi$ and $(\mathbf{A}_1)$–$(\mathbf{A}_3)$ hold. Then we have the following results.

(i) The global Hopf branch $\mathcal{C}_i(E_2, \tau_i, Q_i)$ is bounded for $\tau_i\in\Sigma_H^k$ with $k\ge1$ and $i\in[0, 2L-1]$.

(ii) For any $\tau\in(\min\limits_{k}\Sigma_H^k, \max\limits_{k}\Sigma_H^k)$, model (1.4) possesses at least one periodic solution.

(iii) For $\tau_{i_1}\in\Sigma_H^{k_1}$, $\tau_{i_2}\in\Sigma_H^{k_2}$ with $i_1, i_2\in[0, 2L-1]$ and $k_1, k_2\in[0, \max\limits_{n\in[0, M]} K_n-1]$, we have $\mathcal{C}_{i_1}(E_2, \tau_{i_1}, Q_{i_1})\cap\mathcal{C}_{i_2}(E_2, \tau_{i_2}, Q_{i_2}) = \emptyset$ if $k_1\ne k_2$.

5.   Numerical simulation

To verify obtained theoretical results, the numerical simulation is presented in this part. We choose the fear effect function as $g(K, u_2) = e^{-K u_2}$ and let

Figure 1 shows the existence and stability of $E_0$, $E_1$ and $E_2$, and Hopf bifurcation curve of model (1.4) in $K-\tau$ plane. Above the line $\tau = \tau_{max}$, $E_2$ does not exist, $E_1$ is globally asymptotically stable and $E_0$ is unstable; Below the line $\tau = \tau_{max}$, there exist three constant steady states $E_0, E_1$ and $E_2$. In the region which is bounded by $\tau = \tau_{max}$ and $2mu_1^*+u_2^*(rg'_{u_2}(K, u_2^*)-p) = 0$, no Hopf bifurcation occurs and $E_2$ is stable. In the region which is bounded by $\tau = \tau_0$ and $\tau = \tau_1$, there exist periodic solutions through Hopf bifurcation bifurcating at $E_2$

Fix $K = 1.07$, then by simple calculation, we have $\tau_{max}\approx9.944$, $\tau_0\approx0.85$, $\tau_1\approx3.55$, $J_0 = [0, 5.25]$, $J_1 = [0, 3.05]$, $J_n = \emptyset$ for $n\ge2$, $\sup\limits_{\tau\in J_0} \mathcal{S}_0^0 > 0$ and $\sup\limits_{\tau\in J_1}\mathcal{S}_1^0 < 0$. Thus, all Hopf bifurcation values $\tau_0$ and $\tau_1$ are the all zeros of $S_0^k(\tau)$ for integer $k\ge0$. We summarize the dynamics of model (1.4) as follows.

$(i)$ For $\tau\in[\tau_{max}, \infty)$, we obtain $E_1$ is globally asymptotically stable and $E_0$ is unstable, see Figure 2(a).

$(ii)$ For $\tau\in(0, \tau_0)\cup(\tau_1, \tau_{max})$, we obtain $E_2$ is locally asymptotically stable, and two constant steady state $E_0$ and $E_1$ are unstable, as shown in Figure 2(b).

$(iii)$ For $\tau\in(\tau_0, \tau_1)$, we obtain $E_0, E_1$ and $E_2$ are unstable, a periodic solution bifurcates from $E_2$, as shown in Figure 2(c). Further, a Hopf bifurcation occurs at $E_2$ when $\tau = \tau_0, \text{and}\ \tau_1$.

Next, we explore the global Hopf branches by choosing $\Omega = (0, 4\pi)$ and

As shown in Figure 3, we collect all zeros of $\mathcal{S}_n^k(\tau)$ for nonnegative $n, k$ in set $B_i$ with $i = 0, 1, 2, 3$, namely,

From Theorem 4.3, $E_2$ is locally asymptotically stable when $\tau\in(0, 0.06)\cup(14.04, \tau_{max})$ and unstable when $\tau\in(0.06, 14.04)$, at least one periodic solution emerges for $\tau\in(0.06, 14.04)$. Moreover, a spatially homogeneous periodic solution bifurcates from $\tau\in B_0$, see Figure 4(a); a spatially nonhomogeneous periodic solution bifurcates from $\tau\in B_1\cup B_2\cup B_3$, see Figure 4(b).

In model (1.3), the kernel function takes form as $\mathcal{G} = e^{-\gamma\tau}f(x-y)$ by reasonable assumptions and our theoretical results are derived by choosing $f(x-y)$ as Dirac-delta function. Next, we choose

Here, $f$ is the truncated normal distribution. Clearly, $\int_\Omega f(x-y) dy = 1$. Let $\Omega = (0, 4\pi)$ and the parameter values chosen according to (5.1). As shown in Figure 5, when $\tau = 1.03$, a stable nonhomogeneous periodic solution emerges; when $\tau = 1.5$, a homogeneous periodic solution emerges; when $\tau = 16$, the positive constant steady state is stable. Numerical simulation suggests nonlocal interaction can produce more complex dynamics.

6.   Summary

We formulate an age-structured predator-prey model with fear effect. For $R_0\le1$, the global asymptotic stability for predator-free constant steady state is proved via Lyapunov-LaSalle invariance principle. For $R_0 > 1$, we prove the nonexistence of spatially nonhomogeneous steady states and exclude steady state bifurcation. Finally, we carry out Hopf bifurcation analyses and prove global Hopf branches are bounded.

Our theoretical results are obtained by choosing a special kernel function in model (1.3). However, in numerical results, we explore rich dynamics when the nonlocal interaction is incorporated into the delayed term. The theoretical results concerning the nonlocal model are left as an open problem.

Acknowledgments

H. Shu was partially supported by the National Natural Science Foundation of China (No. 11971285), the Fundamental Research Funds for the Central Universities (No. GK202201002), the Natural Science Basic Research Program of Shaanxi (No. 2023-JC-JQ-03), and the Youth Innovation Team of Shaanxi Universities. W. Xu was partially supported by a scholarship from the China Scholarship Council while visiting the University of New Brunswick. P. Jiang was partially supported by the National Natural Science Foundation of China (No. 72274119).

Conflict of interest

The authors declare no conflicts of interest in this paper.