Hopf Bifurcation Analysis of a Predator–Prey Model with Prey Refuge and Fear Effect Under Non-diffusion and Diffusion

Abstract

In this paper, we propose a predator–prey model with prey refuge and fear effect under non-diffusion and diffusion. For the non-diffusion ODE model, we first analyze the existence and stability of equilibria. Then, the existence of transcritical bifurcation, Hopf bifurcation and limit cycle is discussed, respectively. We find that when the cost of minimum fear $\eta$ is taken as the bifurcation parameter, it not only influence the occurrence of Hopf bifurcation but also alters its direction. For diffusion predator–prey model under homogeneous Neumann boundary conditions, we observe that the Turing instability does not occur, but the Hopf bifurcation will manifest near the interior equilibrium. By considering $\eta$ as the bifurcation parameter, the direction and stability of spatially homogeneous periodic orbits are established. At last, the validity of the theoretical analysis are verified by a series of numerical simulations. The results indicate that prey refuge and fear effect play an key role in the stability of populations.

Appendices

Appendix A

Proof of Theorem 2.1

From model (2), denote

$$\begin{aligned} \left\{ \begin{aligned} H_1(u,v)=&ru\left( \eta +\frac{\alpha (1-\eta ) }{\alpha +v}\right) -\mu _1 u-\gamma u^2-\frac{\beta (1-\delta )uv}{1+a(1-\delta ) u},\\ H_2(u,v)=&\frac{\theta \beta (1-\delta ) uv}{1+a(1-\delta ) u}-\mu _2v. \end{aligned}\right. \end{aligned}$$

(A.1)

It is obvious that we have two boundary equilibrium points

$$\begin{aligned} E_0=(0,0), \quad E_1=(u_1,0)=\left( \frac{r-\mu _1}{\gamma },0\right) , \end{aligned}$$

where $E_0$ is always existed in model (2), $E_1$ is existed when $R_0>1$.

For the positive equilibrium point $E^*=(u^*,v^*)$: if conditions $\mathbf {C.0}$ hold, $u^*=\frac{\mu _2}{(\theta \beta -a\mu _2)(1-\delta )}$ exists and $v^*$ is the positive root of following equation

$$\begin{aligned} v^2+B_1 v+B_2=0, \end{aligned}$$

where $B_1= {\mu _2\theta \gamma }/{(\theta \beta -a\mu _2)^2(1-\delta )^2}+{\theta (\mu _1-r\eta )}/{(\theta \beta -a\mu _2)(1-\delta )}+\alpha$, $B_2={\alpha \theta [\mu _2\gamma -(r-\mu _1)(\theta \beta -a\mu _2)(1-\delta )]}/{(\theta \beta -a\mu _2)^2(1-\delta )^2}$.

For $B_2<0$, that is $\delta ^{*}> \delta \ge 0$, the above equation has a unique positive root given by

$$\begin{aligned} v^*=\frac{-B_1+\sqrt{B_1^2-4 B_2}}{2}. \end{aligned}$$

$\square$

Proof of Theorem 2.2

(I) For equilibrium $E_0$: The corresponding Jacobian matrix of the model (2) is

$$\begin{aligned} J_{E_0}:=\left( \begin{array}{cc} r-\mu _1&{}0\\ 0&{} -\mu _2 \end{array} \right) , \end{aligned}$$

which has two eigenvalues $\lambda _1=r-\mu _1$, $\lambda _2=-\mu _2<0$. If $R_0<1$, then the equilibrium $E_0$ of model (2) is locally asymptotically stable node point. But if $R_0>1$, $E_0$ is a saddle.

(II) For equilibrium $E_1$: The corresponding Jacobian matrix of the model (2) is

$$\begin{aligned} J_{E_1}:=\left( \begin{array}{cc} -r+\mu _1&{} -\frac{r(1-\eta )u_1}{\alpha }-\frac{\beta (1-\delta ) u_1}{1+a(1-\delta )u_1} \\ 0&{}\frac{\theta \beta (1-\delta ) u_1}{1+a(1-\delta )u_1}-\mu _2 \end{array} \right) . \end{aligned}$$

Obviously, we have two eigenvalues $\lambda _1=-r+\mu _1$, $\lambda _2={\theta \beta (1-\delta ) u_1}/{1+a(1-\delta )u_1}-\mu _2$. Noting that $\lambda _1<0$ is equivalent to $R_0>1$ and $\lambda _2<0$ is equivalent to

$$\begin{aligned} \begin{aligned} \mu _2>\frac{1}{\gamma }(\beta \theta -a\mu _2)(r-\mu _1)(1-\delta ). \end{aligned} \end{aligned}$$

(A.2)

When $R_1<1$, that is $\beta \theta <a\mu _2$, which means that the right side of the inequality (A.2) is less than zero. Then we obtain $\lambda _2<0$. When $R_1>1$, that is $\beta \theta >a\mu _2$, we obtain $1>\delta >\delta ^{*}$ by calculating the inequality (A.2). Thus, we obtain conditions for locally asymptotically stability of $E_1$, i.e., $R_1<1$ or $R_1>1$ and $1>\delta >\delta ^{*}$.

(III) For positive equilibrium $E^*=(u^*,v^*)$: The corresponding Jacobian matrix of model (2) is

$$\begin{aligned} J_{E^*}=\left( \begin{array}{cc} a_{11}&{} a_{12}\\ a_{21}&{} a_{22} \end{array} \right) =\left( \begin{array}{cc} -\gamma u^*+\frac{a\beta (1-\delta )^2u^*v^*}{(1+a(1-\delta )u^*)^2}&{}-\left( \frac{r\alpha (1-\eta )u^*}{(\alpha +v^*)^2}+\frac{\beta (1-\delta ) u^*}{1+a(1-\delta )u^*}\right) \\ \frac{\theta \beta (1-\delta ) v^*}{(1+a(1-\delta )u^*)^2}&{} 0 \end{array} \right) .\nonumber \\ \end{aligned}$$

(A.3)

We can obtain the following characteristic equation

$$\begin{aligned} \begin{aligned} \lambda ^2-\textrm{tr}(J_{E^*})\lambda +\textrm{det}(J_{E^*})=0, \end{aligned} \end{aligned}$$

(A.4)

where

$$\begin{aligned} \left\{ \begin{aligned} \textrm{tr}(J_{E^*})=&-\gamma u^*+\frac{a\beta (1-\delta )^2u^*v^*}{(1+a(1-\delta )u^*)^2},\\ \textrm{det}(J_{E^*})=&\frac{\theta \beta (1-\delta ) v^*}{(1+a(1-\delta )u^*)^2} \left( \frac{r\alpha (1-\eta )u^*}{(\alpha +v^*)^2}+\frac{\beta (1-\delta ) u^*}{1+a(1-\delta )u^*}\right) >0. \end{aligned}\right. \nonumber \\ \end{aligned}$$

(A.5)

From (A.5), we have $\textrm{det}(J_{E^*})>0$ always hold. And we know that $\textrm{tr}(J_{E^*})<0$ is equivalent to

$$\begin{aligned} \begin{aligned}\alpha \left( \frac{r(1-\eta )}{\frac{\gamma (\theta \beta +a\mu _2)}{a(\theta \beta -a\mu _2)(1-\delta )}-r\eta +\mu _1}-1\right) <v^*=\frac{-B_1+\sqrt{B_1^2-4 B_2}}{2}. \end{aligned} \end{aligned}$$

(A.6)

When conditions $\mathbf {C.0}$ and $\delta ^{*}> \delta \ge \delta ^{**}$ are satisfied, the left side of the inequality (A.6) is less than zero, then we obtain $\textrm{tr}(J_{E^*})<0$. If the condition $\delta ^{**}> \delta \ge 0$ are satisfied, by calculating the inequality (A.6), we obtain $\eta ^*>\eta \ge 0$ which can make $\textrm{tr}(J_{E^*})<0$. Then we know that the positive equilibrium $E^*$ is a locally asymptotically stable node for $\textrm{tr}^2(J_{E^*})>4\textrm{det}(J_{E^*})$ and stable focus for $\textrm{tr}^2(J_{E^*})<4\textrm{det}(J_{E^*})$.

Now, we consider the global asymptotic stability of the positive equilibrium point $E^*$ under conditions $\mathbf {C.0}$ and $\delta ^{*}> \delta \ge \delta ^{**}$ of Theorem 2.2.

Denote a Dulac function $B(u,v)=(\alpha +v)(1+a(1-\delta )u)u^{-1}v^b$, where b is to be specified later, then

$$\begin{aligned} \begin{aligned} D=&\frac{\partial (H_1(u,v)B(u,v))}{\partial u}+\frac{\partial (H_2(u,v)B(u,v))}{\partial v}\\ =&u^{-1}v^b\left( \frac{1}{\alpha }h_1(u)v+h_2(u)\right) , \end{aligned} \end{aligned}$$

where

$$\begin{aligned} \begin{aligned} h_1(u,b)=&-2a\alpha \gamma (1-\delta )u^2+[(1-\delta )[(2+b)(\theta \beta -a\mu _2)\\&+a(r\eta -\mu _1)]-\gamma ]\alpha u-(2+b)\mu _2\alpha ,\\ h_2(u,b)=&-2a\alpha \gamma (1-\delta )u^2+[(1-\delta )[(1+b)(\theta \beta -a\mu _2)\\&+a(r-\mu _1)]-\gamma ]\alpha u-(1+b)\mu _2\alpha . \end{aligned} \end{aligned}$$

When $1-\frac{\theta \beta -a\mu _2}{ar}\ge \eta \ge 0$ holds, we have

$$\begin{aligned} \begin{aligned} h_2(u,b)-h_1(u,b)=[ar(1-\eta )-(\theta \beta -a\mu _2)]\alpha (1-\delta )u+\mu _2\alpha \ge 0 \end{aligned} \end{aligned}$$

for $u\in {\mathbb {R}}_+$.

Thus, we need to prove $h_2(u,b)\le 0$ for all $u\in {\mathbb {R}}_+$ such that $D\le 0$ for all $(u,v)\in {\mathbb {R}}_+^2$.

Since $a\alpha \gamma (1-\delta )>0$ always holds, then we can obtain that

$$\begin{aligned} \begin{aligned} g(b+1)=g_1(b+1)^2+g_2(b+1)+g_3\le 0 \end{aligned} \end{aligned}$$

such that $D\le 0$ for all $(u,v)\in {\mathbb {R}}_+^2$, where

$$\begin{aligned} \begin{aligned} g_1=&\alpha ^2(\beta \theta -a\mu _2)^2(1-\delta )^2,\\ g_2=&2\alpha ^2(1-\delta )[(\beta \theta -a\mu _2)(a(r-\mu _1 )(1-\delta )-\gamma )-4a\gamma \mu _2],\\ g_3=&\alpha ^2(a(r-\mu _1 )(1-\delta )-\gamma )^2. \end{aligned} \end{aligned}$$

That means that there is a $b+1$ such that

$$\begin{aligned} \begin{aligned} g_2^2-4g_1g_3\ge 0, \end{aligned} \end{aligned}$$

which is equivalent to

$$\begin{aligned} \begin{aligned} 2a\gamma \mu _2\ge (a(r-\mu _1 )(1-\delta )-\gamma )(\beta \theta -a\mu _2). \end{aligned} \end{aligned}$$

This simplifies to $\delta >\delta ^{**}$.

Thus, there exists b such that $D\le 0$ for all $(u,v)\in {\mathbb {R}}_+^2$. By the Bendixson-Dulac theorem [56], we obtain that $E^*$ is globally asymptotically stable under both conditions $\mathbf {C.0}$, $\delta ^{*}> \delta \ge \delta ^{**}$ and $1-{\theta \beta -a\mu _2}/{ar}>\eta \ge 0$. $\square$

Appendix B

Proof of Theorem 2.7

Let $\widehat{u}=u-u^*$ and $\widehat{v}=v-v^*$, then we can rewrite the model (2) by the Taylar expansions as follows (for simplicity, u and v represent $\widehat{u}$ and $\widehat{v}$ again, respectively)

$$\begin{aligned} \left\{ \begin{aligned} \frac{\text {d}{u}}{\text {d} t}=&r(u+u^*)\left( \eta +\frac{\alpha (1-\eta ) }{\alpha +(v+v^*)}\right) \\&-\mu _1 (u+u^*)-\gamma (u+u^*)^2-\frac{\beta (1-\delta )(u+u^*)(v+v^*)}{1+a(1-\delta )(u+u^*)},\\ \frac{\text {d}{v}}{\text {d} t}=&\frac{\theta \beta (1-\delta )(u+u^*)(v+v^*)}{1+a(1-\delta )(u+u^*)}-\mu _2(v+v^*). \end{aligned}\right. \end{aligned}$$

(B.1)

The model (B.1) can be rewritten as

$$\begin{aligned} \begin{aligned} \left( \begin{array}{cc} u_t\\ v_t \end{array} \right) =J_{E^*}\left( \begin{array}{cc} u\\ v \end{array} \right) +\left( \begin{array}{cc} g(u,v,\eta )\\ h(u,v,\eta ) \end{array} \right) , \end{aligned} \end{aligned}$$

(B.2)

where

$$\begin{aligned} \left\{ \begin{aligned} g(u,v,\eta )=&g_{20}u^2+g_{11}uv+g_{02}v^2+g_{30}u^3+g_{21}u^2v+g_{12}uv^2+g_{03}v^3+\cdots ,\\ h(u,v,\eta )=&h_{20}u^2+h_{11}uv+h_{02}v^2+h_{30}u^3+h_{21}u^2v+h_{12}uv^2+h_{03}v^3+\cdots ,\\ \end{aligned}\right. \end{aligned}$$

(B.3)

with

$$\begin{aligned} \begin{aligned} g_{20}=&- \gamma + \frac{a\beta (1-\delta )^2v^*}{(1+a(1-\delta )u^*)^3},\; g_{11}=-\frac{r\alpha (1-\eta )}{(\alpha + v^*)^2} - \frac{\beta (1-\delta )}{(1+a(1-\delta )u^*)^2},\\ g_{02}=&\frac{r\alpha (1-\eta )u^*}{(\alpha + v^*)^3},\\ g_{30}=&-\frac{a^2\beta (1-\delta )^3v^*}{(1+a(1-\delta )u^*)^4},\;g_{21}=\frac{a\beta (1-\delta )^2}{(1+a(1-\delta )u^*)^3},\;g_{12}=\frac{r\alpha (1-\eta )}{(\alpha + v^*)^3},\\ g_{03}=&-\frac{r\alpha (1-\eta )u^*}{(\alpha + v^*)^4}, \\ h_{20}=&-\frac{a\theta \beta (1-\delta )^2v^*}{(1+a(1-\delta )u^*)^3},\;h_{11}=\frac{\theta \beta (1-\delta )}{(1+a(1-\delta )u^*)^2},\;h_{02}=0, \\ h_{30}=&\frac{a^2\theta \beta (1-\delta )^3v^*}{(1+a(1-\delta )u^*)^4},\;h_{21}=-\frac{a\theta \beta (1-\delta )^2}{(1+a(1-\delta )u^*)^3},\;h_{12}=0,\;h_{03}=0. \end{aligned} \end{aligned}$$

We know that the characteristic roots of $J_{E^*}$ are

$$\begin{aligned} \lambda _{1,2}=S\pm iT. \end{aligned}$$

When ${-a_{12}a_{21}>S^2}$, the eigenvalues $\lambda _1$ and $\lambda _2$ will be complex conjugate, especially, when $\eta =\eta ^*$, $\lambda _1$ and $\lambda _2$ will be purely imaginary, that is, $S(\eta ^*)=0$ and $\lambda _{1,2}=\pm i T(\eta ^*)$.

The eigenvectors of $J_{E^*}$ corresponding to the eigenvalues of $\lambda _{1,2}=S + i T$ are given by

$$\begin{aligned} X=\left[ \begin{array}{cc} 1\\ Y-iZ \end{array} \right] , \end{aligned}$$

where Y and Z are defined in (5).

Let the matrix

$$\begin{aligned} G=\left[ \begin{array}{cc} 1&{}0\\ Y&{}Z \end{array} \right] . \end{aligned}$$

By using the transformation

$$\begin{aligned} \left[ \begin{array}{cc} u\\ v \end{array} \right] =G\left[ \begin{array}{cc} x\\ y \end{array} \right] , \end{aligned}$$

then, the model (B.2) becomes

$$\begin{aligned} \begin{aligned} \left[ \begin{array}{cc} \frac{dx}{dt}\\ \frac{dy}{dt} \end{array} \right] =J^*_{E^*}\left[ \begin{array}{cc} x\\ y \end{array} \right] +\left[ \begin{array}{cc} \Phi (x,y)\\ \Psi (x,y) \end{array} \right] , \end{aligned} \end{aligned}$$

(B.4)

where

$$\begin{aligned}{} & {} J^*_{E^*}=\left[ \begin{array}{cc} S&{}-T\\ T&{}S \end{array} \right] ,\\{} & {} \Phi (x,y) = \bigg [- \gamma + \frac{a\beta (1-\delta )^2v^*}{(1+a(1-\delta )u^*)^3}- \frac{\beta (1-\delta )Y}{(1+a(1-\delta )u^*)^2}\\{} & {} \quad +\frac{r\alpha (1-\eta )(u^*Y-\alpha -v^*)Y}{(\alpha + v^*)^3}\bigg ]x^2\\{} & {} \quad -\bigg [ \frac{\beta (1-\delta )Z}{(1+a(1-\delta )u^*)^2}+\frac{r\alpha (1-\eta )(2u^*Y+\alpha +v^*)}{(\alpha + v^*)^3}Z\bigg ]\\{} & {} xy+\frac{r\alpha (1-\eta )u^*}{(\alpha + v^*)^3}Z^2y^2\\{} & {} +\bigg [-\frac{a^2\beta (1-\delta )^3v^*}{(1+a(1-\delta )u^*)^4}+\frac{a\beta (1-\delta )^2 Y}{(1+a(1-\delta )u^*)^3}-\frac{r\alpha (1-\eta )(u^*Y-\alpha -v^*)}{(\alpha + v^*)^4}Y^2\bigg ]x^3\\{} & {} -\frac{r\alpha (1-\eta )(3u^*Y-\alpha -v^*)}{(\alpha + v^*)^4}Z^2 xy^2+\cdots , \Psi (x,y)\\{} & {} =\frac{1}{Z}\bigg [-\frac{a\beta (1-\delta )^2v^*(\theta +Y)}{(1+a(1-\delta )u^*)^3}+\gamma Y+ \frac{\beta (1-\delta )(\theta +Y)Y}{(1+a(1-\delta )u^*)^2}\\{} & {} -\frac{r\alpha (1-\eta )(u^*Y-\alpha -v^*)}{(\alpha + v^*)^3}Y^2\bigg ]x^2\\{} & {} +\bigg [\frac{\beta (1-\delta )(\theta +Y)}{(1+a(1-\delta )u^*)^2}-\frac{r\alpha (1-\eta )(2u^*Y+\alpha +v^*)}{(\alpha + v^*)^3}Y\bigg ]xy\\{} & {} -\frac{r\alpha (1-\eta )u^*}{(\alpha + v^*)^3}YZy^2\\{} & {} +\frac{r\alpha (1-\eta )u^*}{(\alpha + v^*)^4}YZ^2y^3+\bigg [-\frac{a\beta (1-\delta )^2(\theta +Y)}{(1+a(1-\delta )u^*)^3}\\{} & {} +\frac{r\alpha (1-\eta )(3u^*Y-2\alpha -2v^*)}{(\alpha + v^*)^4}Y^2\bigg ]x^2y+\cdots . \end{aligned}$$

The model (B.4) can be written in the polar form as

$$\begin{aligned} \begin{aligned}&\frac{dr}{dt}=S(\eta )r+m(\eta )r^3+\cdots ,\\&\frac{d\theta }{dt}=T(\eta )+h(\eta )r^2+\cdots . \end{aligned} \end{aligned}$$

(B.5)

The Taylors series expansion of (B.5) at $\eta =\eta ^*$ gives

$$\begin{aligned} \begin{aligned}&\frac{dr}{dt}=S'(\eta ^*)(\eta -\eta ^*)r+m(\eta ^*)r^3+\cdots ,\\&\frac{d\theta }{dt}=T(\eta )+T'(\eta ^*)(\eta -\eta ^*)+h(\eta ^*)r^2+\cdots . \end{aligned} \end{aligned}$$

where $m(\eta ^*)$ is defined in (4).

Owing to

$$\begin{aligned} \begin{aligned} \textrm{sgn}{\left[ \frac{\partial (\textrm{Re}S)}{\partial \eta }\right] }_{\eta =\eta ^*}=\frac{ra\mu _2}{4\theta \beta }\frac{\sqrt{B_1^2-4 B_2}-B_1}{\sqrt{B_1^2-4 B_2}}|_{\eta =\eta ^*}>0. \end{aligned} \end{aligned}$$

Therefore, we have $\Lambda =-\frac{m(\eta ^*)}{S'(\eta ^*)}$, which determines the stability of Hopf bifurcating periodic solution. $\square$