Dynamics in a predator-prey model with memory effect in predator and fear effect in prey

1.   Introduction

The relationship between predator and prey is an important research content in ecosystem, and many scholars have studied this interaction by differential equation models ^{[1,2,3,4,5,6]}. The direct relationship between predator and prey is the consumption of prey. One of the indirect relationship between predator and prey is the fear of predator. When the prey gets the predator's signal (chemical/vocal), they will increase their vigilance time and reduce their foraging ^[7], such as mule deer v.s. mountain lions ^[8], elk v.s. wolves ^[9]. Consider the fear effect, Panday et al. ^[10] proposed the following model

All parameters are positive. The biological interpretation of parameters is given in Table 1. They incorporated fear effect by modifying the prey intrinsic growth rate $R_0$ as $\frac{R_0}{1+Kv}$. By the scaling

model (1.1) is changed to (drop the bars)

Panday et al. ^[10] introduced time delay of perceiving predator signals to model (1.3), and mainly studied the boundedness, persistence, local and global behavior of the delayed model.

In the real world, in addition to the fear effect of the prey, the clever predator also has spatial-memory and cognition ^[11], which is often ignored in modeling the predator-prey interaction. For example, blue whales rely on memory for migration, which is presented by B. Abrahms et al. ^[12] and W. F. Fagan ^[13]. As another example, animals in polar regions usually determine their spatial movement by judging footprints, which record the history of species distribution and movement, involving time delay ^[14]. Obviously, highly developed animals can even remember the historical distribution or cluster of species in space. Much progress has been made in implicitly integrating spatial cognition or memory ^{[14,15,16,17,18]}. To incorporate the memory effect, Shi et al. proposed a single specie model with spatial memory by introducing a additional delayed diffusion term ^[14]. They supposed that in addition to the negative gradient of the density distribution function at the present time, there is a directed movement toward the negative gradient of the density distribution function at past time ^[14]. After this pioneering work, some recent works ^{[19,20,21,22,23]} about the population model with memory effect have emerged. In ^[23], Song et al. obtained a computing method for the normal forms of the Hopf bifurcations in the diffusive predator-prey model with memory effect, which is friendly to use.

Inspired by the above work, we suppose the predator has spatial-memory diffusion and the prey has fear effect, then modified the model (1.3) as follow

All parameters are positive. The biological description of parameters is given in Table 1. The term $-d\nabla (v \nabla u(t-\tau))$ represents the memory-based diffusion effect of the predator. The Neumann boundary conditions is used. The aim of this paper is to study the effect of predator's memory-based diffusion and prey's fear on the model (1.4).

The rest of this paper is organized as follows. In Section 2, the stability of coexisting equilibrium and the existence of Hopf bifurcation are considered. In Section 3, the property of Hopf bifurcation is studied. In Section 3, some numerical simulations are given to analyze the effect of spatial-memory and fear effect. In Section 4, a short conclusion is given.

2.   Stability analysis

For simplicity, we choose $\Omega = (0, l\pi)$. Denote $\mathbb{N}$ as positive integer set, and $\mathbb{N}_0$ as nonnegative integer set. It is easy to obtain $(0, 0)$ and $(1, 0)$ are two boundary equilibria of model (1.4). Next, we will give the existence of coexisting equilibrium.

Lemma 2.1. If $a > (1+b) \mu$, model (1.4) has one unique coexisting equilibrium $(u_*, v_*)$ where $u_* = \frac{\mu }{a-b \mu }$, $v_* = \frac{1}{2 c}\left(-1+\sqrt{1+\frac{4 c [a-(1+b) \mu ]}{(a-b \mu)^2}}\right)$.

Proof. The coexisting equilibrium of (1.4) is a positive root of the following equations

From the second equation, we have $u = \frac{\mu }{a-b \mu }$. Substitute it into the first equation, we have

Then $v = \frac{1}{2 c}\left(-1\pm \sqrt{1+4 c \frac{(a-(1+b) \mu)}{(a-b \mu)^2}}\right)$. Obviously, the conclusion holds.

In this paper, we mainly study the stability of coexisting equilibrium $E_*(u_*, v_*)$. Linearize model (1.4) at $E_*(u_*, v_*)$, we have

where

and $a_1 = u_*\left(\frac{a b v_*}{(1+b u_*)^2}-\frac{1}{1+c v_*}\right)$, $a_2 = -\frac{u_*(1-u_*)(1+2 c v_*)}{v_*(1+c v_*)^2} < 0$, $b_1 = \frac{a v_*}{(1+b u_*)^2} > 0$. The characteristic equations are

where

2.1.   $\tau=0$

When $\tau = 0$, the characteristic Eq (2.3) become

where $B_n+C_n = -a_2 b_1-(a_2 d v_*+a_1 d_2)\frac{ n^2}{l^2}+d_1 d_2 \frac{n^4}{l^4}$. Make the following hypothesis

Theorem 2.1. For model (1.4) with $\tau = 0$, $E_*(u_*, v_*)$ is locally asymptotically stable under the hypothesis ${\bf (H_1)}$.

Proof. If ${\bf (H_1)}$ holds, we can easily obtain that $A_n > 0$ and $B_n+C_n > 0$. Then the characteristic roots of (2.4) all have negative real parts. Then $E_*(u_*, v_*)$ is locally asymptotically stable.

2.2.   $\tau > 0$

In the following, we assume ${\bf (H_1)}$ holds. Let $\mbox{i}\omega$ ($\omega > 0$) be a solution of Eq (2.3), then we have

We can obtain ${\rm{cos}}\omega\tau = \frac{\omega ^2-B_n}{C_n}$, ${\rm{sin}}\omega\tau = \frac{A_n\omega }{C_n} > 0$ under hypothesis ${\bf (H_1)}$. It leads to

Let $p = \omega^2$, then (2.5) becomes

and the roots of (2.6) are $p_n^{\pm} = \frac{1}{2}[-\left(A_n^2-2 B_n\right) \pm \sqrt{\left(A_n^2-2 B_n\right)^2-4(B_n^2-C_n^2)}].$ By direct computation, we have

and $B_n+C_n > 0$ under hypothesis ${\bf (H_1)}$. Define $z_{\pm} = \frac{-(a_2 d v_*-a_1 d_2)\pm\sqrt{(a_2 d v_*-a_1 d_2)^2-4 d_1 d_2 (-a_2 b_1)}}{2 d_1 d_2}$, $d_* = \frac{a_1 d_2}{a_2v_*}+\frac{2}{v_*} \sqrt{-\frac{b_1 d_1 d_2}{a_2}}$, and $\mathbb{M} = \{n|\frac{n^2}{l_2}\in (z_-, z_+), \; n\in\mathbb{N}_0 \}$. Then we can obtain that

The existence of purely imaginary roots of Eq (2.3) can be divided into the following two cases.

${\bf Case\; 1:}\; a_1^2+2 a_2 b_1 > 0$. We can obtain $A_n^2-2B_n > a_1^2+2 a_2 b_1 > 0$. For $d > d_*$ and $n\in \mathbb{M}$, Eq (2.3) has a pair of purely imaginary roots $\pm i\omega^+_n$ at $\tau^{j, +}_n$ for $j \in \mathbb{N}_0$ and $n\in\mathbb{M}$. Otherwise, Eq (2.3) does not have characteristic roots with zero real parts.

${\bf Case\; 2:}\; a_1^2+2 a_2 b_1 < 0$. This case can be divided into the following two subcases.

$\bullet$ For $d \leq d_*$ and $n\in \mathbb{M}_1: = \{n|A_n^2-2 B_n < 0, \; \left(A_n^2-2 B_n\right)^2-4(B_n^2-C_n^2) > 0, \; n\in\mathbb{N}_0 \}$, Eq (2.3) has two pairs of purely imaginary roots $\pm i\omega^\pm_n$ at $\tau^{j, \pm}_n$ for $j \in \mathbb{N}_0$ and $n\in\mathbb{M}_1$. Otherwise, Eq (2.3) does not have characteristic roots with zero real parts.

$\bullet$ For $d > d_*$ and $n\in \mathbb{M}_2: = \{n|A_n^2-2 B_n < 0, \; \left(A_n^2-2 B_n\right)^2-4(B_n^2-C_n^2) > 0, \; n\in\mathbb{N}_0, \; n\notin\mathbb{M} \}$, Eq (2.3) has two pairs of purely imaginary roots $\pm i\omega^\pm_n$ at $\tau^{j, \pm}_n$ for $j \in \mathbb{N}_0$ and $n\in\mathbb{M}_1$. For $d > d_*$ and $n\in \mathbb{M}$, Eq (2.3) has a pair of purely imaginary roots $\pm i\omega^+_n$ at $\tau^{j, +}_n$ for $j \in \mathbb{N}_0$ and $n\in\mathbb{M}$. Otherwise, Eq (2.3) does not have characteristic roots with zero real parts.

The $\omega^\pm_n$ and $\tau^{j, \pm}_n$ are defined as follow

Define

We have the following lemma.

Lemma 2.2. Assume ${\bf (H_1)}$ holds. Then$\mathit{{\rm{Re}}}(\frac{d \lambda} {d \tau})|_{\tau = \tau^{j, +}_n} > 0$, $\mathit{{\rm{Re}}} (\frac{d \lambda} {d \tau})|_{\tau = \tau^{j, -}_n} < 0$ for $\tau^{j, \pm}_n \in \mathbb{S}$ and $j \in \mathbb{N}_0$.

Proof. By Eq (2.3), we have

Then

Therefore ${\rm{Re}}(\frac{d \lambda} {d \tau})|_{\tau = \tau^{j, +}_n} > 0$, ${\rm{Re}} (\frac{d \lambda} {d \tau})|_{\tau = \tau^{j, -}_n} < 0$.

Denote $\tau_* = min\{\tau^{0, \pm}_n|\; \tau^{0, \pm}_n\in\mathbb{S}\}$. We have the following theorem.

Theorem 2.2. Assume ${\bf (H_1)}$ holds, then the following statements are true for model (1.4).

$\bullet$ $E_*(u_*, v_*)$ is locally asymptotically stable for $\tau > 0$ when $\mathbb{S} = \varnothing$.

$\bullet$ $E_*(u_*, v_*)$ is locally asymptotically stable for $\tau \in [0, \tau_*)$ when $\mathbb{S}\neq \varnothing$.

$\bullet$ $E_*(u_*, v_*)$ is unstable for $\tau \in (\tau_*, \tau_*+\varepsilon)$ for some $\varepsilon > 0$ when $\mathbb{S}\neq \varnothing$.

$\bullet$ Hopf bifurcation occurs at$(u_*, v_*)$ when $\tau = \tau^{j, +}_n$ $(\tau = \tau^{j, -}_n)$, $j\in \mathbb{N}_0$, $\tau^{j, \pm}_n \in \mathbb{S}$.

Remark 2.1. In the Theorem 2, if $E_*(u_*, v_*)$ is locally asymptotically stable, the densities of prey and predator will tend to the equilibrium state in the whole region when the initial densities of prey and predator is near $E_*(u_*, v_*)$. When Hopf bifurcation occurs at$(u_*, v_*)$, then the densities of prey and predator will produce periodic oscillation. Especially, spatially homogeneous periodic oscillations may occur when $\tau$ near the critical value $\tau = \tau^{j, +}_0$ or $\tau = \tau^{j, -}_0$, and spatially inhomogeneous periodic oscillations may occur when $\tau$ near the critical value $\tau = \tau^{j, +}_n$ or $\tau = \tau^{j, -}_n$ ($n > 0$).

3.   Property of Hopf bifurcation

In this section, we use the algorithm in ^[23] to compute the normal form of Hopf bifurcation. We denote the critical value of Hopf bifurcation as $\tilde{\tau}$ and the purely imaginary roots as $\pm i \omega_n$ of Eq (2.3). Let $\bar{u}(x, t) = u(x, \tau t)-u_*$ and $\bar{v}(x, t) = v(x, \tau t)-v_*$. Drop the bar, the model (1.4) can be written as

Define the real-valued Sobolev space $\mathbb{X} = \left\{U = (u, v)^T\in W^{2, 2}(0, l\pi)^2, \; (\frac{\partial{u}}{\partial{x}}, \frac{\partial{v}}{\partial{x}})|_{x = 0, l\pi} = 0 \right\}$, the inner product

and $\mathbb{C} = C([-1, 0];\mathbb{X})$. Set $\tau = \tilde{\tau}+\varepsilon$, where $\varepsilon$ is small perturbation. Then system (3.1) is rewritten as

where for $\varphi = (\varphi_, \varphi_2)^T \in \mathbb{C}$, $d(\varepsilon)\Delta$, $L(\varepsilon):\mathbb{C}\rightarrow \mathbb{X}$, $F:\mathbb{C}\times \mathbb{R}^2\rightarrow \mathbb{X}$. They are defined as

and

Denote $L_0(\varphi) = \tilde{\tau}A\varphi(0)$, and rewrite (3.2) as

where $\tilde{F}(\varphi, \varepsilon) = \varepsilon A\varphi(0)+F(\varphi, \varepsilon)+F^d(\varphi, \varepsilon)$. The characteristic equation for the linearized equation $\frac{dU(t)}{dt} = d_0\Delta (U_t)+L_0 (U_t)$ is $\tilde{\Gamma}_n(\lambda) = det\left(\tilde{M}_n((\lambda))\right)$, where

The eigenvalue problem

has eigenvalues $\frac{n^2}{l^2}$ and normalized eigenfunctions

Set $\beta_n^{(j)} = z_n(x)e_j$, $j = 1, 2$, where $e_1 = (1, 0)^T$ and $e_2 = (0, 1)^T$. Define $\eta_n(\theta)\in BV([-1, 0], \mathbb{R}^2)$, such that

$C = C([-1, 0], \mathbb{R}^2)$, $C^* = C([0, 1], \mathbb{R}^{2*})$, and

Let $\bigwedge = \{i\tilde{\omega}, -i\tilde{\omega} \}$, the eigenspace $P$, and corresponding adjoint space $P^*$. Decompose $C = P\oplus Q$, where $Q = \{\varphi\in C: < \psi, \varphi > = 0, \forall \psi\in P^* \}$. Choose $\Phi(\theta) = (\phi(\theta), \bar{\phi}(\theta))$, $\Psi(\theta) = col(\psi^T(s), \bar{\psi}^T(s))$, where

and

Then $\phi(\theta)$ and $\psi(s)$ are the bases of $P$ and $P^*$, respectively, and such that $< \phi, \psi > = I_2$.

By direct computation, we have

where $f^{(1)}_{20} = \frac{2 a b v_*}{(1+b u_*)^3}-\frac{2}{1+c v_*}$, $f^{(1)}_{11} = -\frac{a}{(1+b u_*)^2}+\frac{c (-1+2 u_*)}{(1+c v_*)^2}$, $f^{(1)}_{02} = \frac{2 c^2 (1-u_*) u_*}{(1+c v_*)^3}$, $f^{(1)}_{30} = -\frac{6 a b^2 v_*}{(1+b u_*)^4}$, $f^{(1)}_{21} = \frac{2 a b}{(1+b u_*)^3}+\frac{2 c}{(1+c v_*)^2}$, $f^{(1)}_{12} = \frac{2 c^2 (1-2 u_*)}{(1+c v_*)^3}$, $f^{(1)}_{03} = -\frac{6 c^3 (1-u_*) u_*}{(1+c v_*)^4}$, $f^{(2)}_{20} = -\frac{2 a b v_*}{(1+b u_*)^3}$, $f^{(2)}_{11} = \frac{a}{(1+b u_*)^2}$, $f^{(2)}_{02} = 0$, $f^{(2)}_{30} = \frac{6 a b^2 v_*}{(1+b u_*)^4}$, $f^{(2)}_{21} = -\frac{2 a b}{(1+b u_*)^3}$, $f^{(2)}_{12} = 0$, $f^{(2)}_{03} = 0$. We can compute the following parameters

and $\widetilde{A}_{j_1 j_2} = A_{j_1 j_2}-2\frac{n^2}{l^2}A^d_{j_1 j_2}$ for $j_1, j_2 = 0, 1, 2, \; j_1+j_2 = 2$. In addition, $h_{0, 20}(\theta) = \frac{1}{l\pi}(\widetilde{M}_0(2i\tilde{\omega}))^{-1}A_{20}e^{2i\tilde{\omega} \theta}$, $h_{0, 11}(\theta) = \frac{1}{l\pi}(\widetilde{M}_0(0))^{-1}A_{11}$, $h_{2n, 20}(\theta) = \frac{1}{2l\pi}(\widetilde{M}_{2n}(2i\tilde{\omega}))^{-1}\widetilde{A}_{20}e^{2i\tilde{\omega} \theta}$, $h_{2n, 11}(\theta) = \frac{1}{l\pi}(\widetilde{M}_{2n}(0))^{-1}\widetilde{A}_{11}$.

Then we have

where $b_{2n}^{(1)} = -\frac{n^2}{l^2}$, $b_{2n}^{(2)} = -2\frac{n^2}{l^2}$, $b_{2n}^{(3)} = -4\frac{n^2}{l^2}$. The normal form of the Hopf bifurcation is

where

By the coordinate transformation $z_1 = \omega_1-i\omega_2$, $z_2 = \omega_1+i\omega_2$, and $\omega_1 = \rho {\rm{cos}}\xi$, $\omega_2 = \rho {\rm{sin}}\xi$, the normal form (3.8) can be rewritten as

where $K_1 = \frac{1}{2}{\rm{Re}}(B_1)$, $K_2 = \frac{1}{3!}{\rm{Re}}(B_2)$.

By the work ^[23], we have the following theorem.

Theorem 3.1. If $K_1K_2 < 0 (>0)$, the Hopf bifurcation is supercritical (subcritical), and the bifurcating periodic orbits is stable(unstable) for $K_2 < 0 (>0)$.

Remark 3.1. In the Theorem 3.1, when Hopf bifurcation is supercritical (subcritical), then the bifurcating periodic solutions exist for $\tau > \tilde{\tau}$ ($\tau < \tilde{\tau}$), where $\tilde{\tau}$ is some critical value $\tau = \tau^{j, \pm}_n$. When the periodic solution is stable, the densities of prey and predator will produce periodic oscillation, and finally continue to oscillate.

4.   Numerical simulations

In this section, we give some numerical simulations to analyze the effect of spatial memory in predator and fear in prey on the model (1.4). Fix the following parameters

4.1.   The effect of $d$

If we choose $c = 1$, then model (1.4) has a unique coexisting equilibrium $(u_*, v_*)\approx(0.6667, \; 0.0.6667)$, and $a_1^2+2 a_2 b_1 \approx 0.0352 > 0$, $d_*\approx 0.6206$. To study the effect of memory-based diffusion coefficient $d$ on the model (1.4), we give the bifurcation diagram of model (1.4) with parameter $d$ as in Figure 1. By the Theorem 2.2, we know that $(u_*, v_*)$ is locally stable for $\tau \geq 0$ when $d < d_*$. But when $d > d_*$, the inhomogeneous Hopf bifurcation curves exist. This means that increasing parameter $d$ is not conducive to the stability of the equilibrium $(u_*, v_*)$, and the densities of prey and predator will produce spatially inhomogeneous periodic oscillation.

4.2.   The effect of $\tau$

Choose $d = 0.7$, we have $\mathbb{M} = \{2, \; 3 \}$ and $\tau_* = \tau_2^{0, +} \approx17.4593 < \tau_3^{0, +} \approx19.1380$. When $\tau\in [0, \tau_*)$, $(u_*, v_*)$ is locally stable (Figure 2). By direct calculation, we can obtain $K_1\approx 0.0166$, and $K_2\approx -0.0699$. Then, the Hopf bifurcation is supercritical and the bifurcating periodic solution is stable (Figure 3). At this time, the bifurcating periodic solution is spatially inhomogeneous and with mode-2.

Choose $d = 0.8$, we have $\mathbb{M} = \{2, \; 3 \}$ and $\tau_* = \tau_3^{0, +} \approx8.4754 < \tau_2^{0, +} \approx9.9645$. When $\tau\in [0, \tau_*)$, $(u_*, v_*)$ is locally stable (Figure 4). By direct calculation, we can obtain $K_1\approx 0.1236$, $K_2\approx -0.3994$. Then, the Hopf bifurcation is supercritical and the bifurcating periodic solution is stable (Figure 5). At this time, the bifurcating periodic solution is spatially inhomogeneous with mode-3. When $\tau_* = \tau_3^{0, +} < \tau_2^{0, +} < \tau = 10$, there is an unstably spatially inhomogeneous periodic solution with mode-2 which transitions to the stably spatially inhomogeneous periodic solution with mode-3 (Figure 6).

4.3.   The effect of $c$

Next, we will study the effect of fear effect $c$ on the model (1.4). Fix the parameters as (4.1), then model (1.4) has a unique coexisting equilibrium $(u_*, v_*)$. And $a_1^2+2 a_2 b_1 > 0$ when $0 < c < 2.6194$. We give the figure of $d_*$ with parameter $c$ as in Figure 7. Set parameter $d = 0.7$ and $d = 0.8$, we give the bifurcation diagrams of model (1.4) with parameter $c$ as in Figure 8.

When $d = 0.7$ and $\tau = 20$, increasing parameter $c$ can destroy the stability of the coexisting equilibrium $(u_*, v_*)$, and induce spatially inhomogeneous periodic solution (Figure 9). This means that increasing parameter $c$ is not conducive to the stability of the coexisting equilibrium $(u_*, v_*)$.

When $d = 0.8$ and $\tau = 9$, increasing parameter $c$ can destroy the stability of the coexisting equilibrium $(u_*, v_*)$, and induce spatially inhomogeneous periodic solution initially (Figure 10). But when $c$ is larger enough, increasing parameter $c$ can rule out the spatially inhomogeneous periodic oscillation and stabilize the coexisting equilibrium $(u_*, v_*)$ (Figure 10). This means that increasing parameter $c$ is not conducive to the stability of the equilibrium $(u_*, v_*)$, initially. But when $c$ is large, increasing parameter $c$ is conducive to the stability of the coexisting equilibrium $(u_*, v_*)$.

5.   Conclusions

In this paper, we incorporate the memory effect in predator and fear effect in prey into a predator-prey model. By using time delay in the memory of predator as bifurcating parameter, we analyze the local stability of coexisting equilibrium, the existence of Hopf bifurcation, and the property of Hopf bifurcation by the method in ^[23]. Through the numerical simulations, we analyzed the effect of memory effect in predator and fear in prey on the model.

The spatial memory effect plays an important role in the dynamics of the predator-prey model. Through the numerical simulations, we observed that the memory-based diffusion coefficient $d$ has destabilizing effect on the predator-prey model when it is larger than some critical value. In addition. when $d$ crosses the critical value, time delay $\tau$ in the memory of predator can affect the stability of the equilibrium $(u_*, v_*)$. In the numerical simulations, we observe that the first Hopf bifurcation curve is inhomogeneous bifurcation curve, and homogeneous Hopf bifurcation curve does not exist. This is different from the predator-prey model without the spatial memory effect. When $\tau$ crosses the critical value $\tau_*$, the densities of prey and predator will produce spatially inhomogeneous periodic oscillation. When $\tau$ crosses the second critical value, the spatially inhomogeneous periodic oscillations with different modes exist, but the densities of prey and predator will converge to the spatially inhomogeneous periodic solution corresponding to the first bifurcation curve. This shows that the spatially memory effect in predator can destroy the stability of the coexisting equilibrium, and induce spatially inhomogeneous periodic oscillations.

In addition, the fear effect parameter $c$ in prey can also affect the stability of the coexisting equilibrium $(u_*, v_*)$. A small fear effect parameter $c$ means a large birth rate $\frac{1}{1+cv}$, then the large birth rate can support fluctuations. Increasing parameter $c$ can destroy the stability of the coexisting equilibrium $(u_*, v_*)$, and induce spatially inhomogeneous periodic solution. Hence, we observed the destabilizing effect on the the coexisting equilibrium $(u_*, v_*)$. A large fear effect parameter $c$ means a low birth rate, then the low birth rate can not support fluctuations. Increasing parameter $c$ can rule out the spatially inhomogeneous periodic oscillation and stabilize the coexisting equilibrium $(u_*, v_*)$. Hence, we observed the stabilizing effect on the the coexisting equilibrium $(u_*, v_*)$.

Acknowledgements

This research is supported by the Fundamental Research Funds for the Central Universities (Grant No. 2572022BC01) and Postdoctoral Program of Heilongjiang Province (No. LBH-Q21060).

Conflict of interest

The authors have no relevant financial or non-financial interests to disclose.