Analysis of a Delayed Reaction-Diffusion Predator–Prey System with Fear Effect and Anti-Predator Behaviour

Abstract

This paper is devoted to studying the dynamics of a delayed reaction-diffusion predator–prey system incorporating the effects of fear and anti-predator behaviour. First, based on its mathematical model, the global attractor is analyzed and the local stability of its positive equilibria is derived. Moreover, the Hopf bifurcation induced by the time delay variable is also investigated. Furthermore, the existence and non-existence of non-constant positive solutions are analyzed. Finally, numerical simulations are presented to validate the theoretical analysis.

1. Introduction

Since it was proposed in the 1920s [1,2], the classic predator–prey model has attracted increasing attention and become a hot topic in both theoretical ecology and biomathematics [3,4,5]. Based on it, some different forms of the Lotka–Volterra model have been proposed by introducing new factors between the predator and the prey. Moreover the corresponding dynamics such as the non-negative solutions, the global stability, ect., have been analyzed in depth.

In most predator–prey models, the predator is always regarded as a winner. However, in the real world, the prey possesses many adaptations developed through evolution to struggle against the predator, which are known as anti-predator behaviors [6]. There exist many different activities operated by the prey in the anti-predator adaptations. For example, in order to avoid detection, the chameleon is able to change its skin color and the purplefrog always lives underground four meters deep. Furthermore, some prey species can defend themselves against predators by communal defence, ejecting noxious materials or even killing juvenile predators. Ives and Dobson [7] found that the increased efficiency of anti-predator behavior increases prey density, decreases the density ratio of predator and prey species, and dampens the corresponding oscillatory dynamics. Tang and Xiao [8] considered anti-predator behaviour in which juvenile predators can be attacked by the prey species. Their results show that the region of the coexistence of the prey and predator populations can be changed by anti-predator behaviour, and moreover, the predator–prey oscillations can also be diminished.

It is known that fear of predators is also a type of anti-predator behaviour and can influence the dynamics of a system. Wang et al. [9] introduced a new term $\frac{1}{1+ky}$ to a Leslie–Gower model, in considering the effects of fear on the reproduction of the prey. Their results described the influence of the lever of fear factor for the species populations. For example, an increase in fear level leads to a decrease in prey density. Zhang et al. [10] studied a similar predator–prey model in which the fear effect and prey refuge are considered, and showed that the fear effect leads to a reduction in predator density as well as to a stabilization of the system. Qi and Meng [11] introduced the fear factor into a stochastic predator–prey model with prey refuge. They studied its threshold behavior and showed that decreasing the cost of fear leads to an increase in the survival rate of prey. Recently, Qiao et al. [12] studied a Beddington–DeAngelis model in considering the fear effect of predators. Das and Samanta [13] studied a prey–predator model with logistic growth of the prey in the absence of the predator, and also considered the fear effect and investigated the impact of fear of the predator on prey when the predator is provided additional food. Panday et al. [14] proposed a delayed predator–prey model with fear in the prey population by considering that the growth rate of the prey population is suppressed due to the fear of predators. Dubey et al. [15] proposed a prey–predator model with fear in prey due to predator and anti-predator behavior by prey against the predator with fear response delay and gestation delay. Yang and Ma [16] investigated the dynamics of a diffusive predator–prey system with anti-predator behaviour and maturation delay with Neumann boundary condition. Zhang et al. [17] considered a delayed diffusive predator–prey model with fear effect.

On the other hand, more and more researchers realize that the dynamics of predator–prey models are significantly influenced by the factors of delay and diffusion. In fact, the predator and prey species are spatially heterogeneous due to the fact that in regions with higher population density [18,19], resources will be depleted more quickly and consequently individuals will move to regions with lower population density [20] in order to obtain sufficient resources. In addition, time delay also plays a key role in predator–prey systems and should not be ignored. In fact, time delay can lead to a loss of stability, the appearance of various oscillations and periodic solutions, etc. Therefore, it is important and reasonable to incorporate the factors of delay and diffusion in the dynamics of predator–prey systems [21,22,23,24,25]. Ma et al. [21] dealt with a delayed three-species Lotka–Volterra food chain model with diffusion effects and homogeneous Neumann boundary conditions. Guo et al. [22] studied the local stability and Hopf bifurcation of a food chain model with two time delays and diffusion. Zhang and Zhao [25] investigated a delayed reaction-diffusion three-species Lotka–Volterra model with interval biological parameters and harvesting. Cai et al. [24] investigated the dynamical outcomes of a host–parasite model incorporating both horizontal and vertical transmissions in a spatial heterogeneous environment analytically and numerically. Zhang et al. [25] develop a diffusive plant invasion model with delay under the homogeneous Neumann boundary condition. However, we noted that the establishment of these models did not consider the impact of the fear effect.

Notice that in real life, the fear effect of a predator on prey is never instantaneous. Thus it should be assumed that there exists a time delay, denoted $\tau$, which affects the growth rate of a prey population. Let u and v denote the prey and predator population densities, respectively. Based on the above discussion, in considering the fear effect and anti-predator behaviour, a delayed reaction-diffusion predator–prey system will be proposed as follows

$$\left\{\begin{array}{cc}{\displaystyle \frac{\partial u}{\partial t}=d_1\Delta u+\frac{ru}{1+bv(t-\tau )}(1-\frac{u}{K})-\frac{\alpha uv}{a+u},}\hfill & (x,t)\in \Omega \times (0,+\infty ),\hfill \\ {\displaystyle \frac{\partial v}{\partial t}=d_2\Delta v+\frac{\beta uv}{a+u}-mv-\eta uv,}\hfill & (x,t)\in \Omega \times (0,+\infty ),\hfill \\ u(x,t)=\phi (x,t),v(x,t)=\psi (x,t),\hfill & (x,t)\in \overline{\Omega }\times [-\tau ,0]\hfill \\ {\displaystyle \frac{\partial u(x,t)}{\partial \nu }=\frac{\partial v(x,t)}{\partial \nu }=0,}\hfill & (x,t)\in \partial \Omega \times (0,+\infty ),\hfill \end{array}\right.$$

(1)

where $\Delta$ denotes the Laplacian operator and $d_1,d_2>0$ are the diffusion coefficients associated with the prey and predator densities, respectively. The term $1/(1+bv(t-\tau )$ is the fear effect term and b is the level of fear, $ru(1-u/K)$ is a logistic growth term, r is the intrinsic growth rate, and K is the carrying capacity of the prey. The term $\frac{uv}{a+u}$ is a Holling type-II functional response and a is a positive constant measuring the ability of a generic predator to kill and consume generic prey, $\alpha$ is the capture rate of the predator, $\beta$ is the conversion rate of prey into predator, and m is the death rate of the predator. The term $\eta uv$ models anti-predator behaviour, and $\eta$ is the rate of anti-predator behaviour of prey for the predator population. Assume that $\Omega$ is bounded and let $\nu$ denote the outward-pointing normal vector of $\partial \Omega$. The homogeneous Neumann boundary conditions indicate that there is no population flux across the boundaries.

Assume the initial conditions

$$\phi (x,t),\psi (x,t)\in C(X,[-\tau ,0\left]\right),$$

while

$$X=\left\{u,v\in W^{2,2}(\Omega ):\frac{\partial u}{\partial \nu }=\frac{\partial v}{\partial \nu }=0,x\in \Omega \right\}.$$

The rest of this paper is organized as follows. The dissipative property of the time-dependent solution of system (1) is studied in Section 2. The global attractor, the local stability and the Hopf bifurcation are analyzed in Section 3. The existence of non-constant positive equilibria is studied in Section 4. In Section 5, numerical simulations are given to validate the theoretical results. Conclusions are given in Section 6.

2. Global Attractor

In this section, we use the upper–lower method and comparison principle to verify the existence and boundedness of the solution of the model (1). First, we give the following result, which was proposed in [26] (Theorem 3.1 by using $L^1$ bounded).

Lemma 1.

Given the initial/boundary value system

$$\left\{\begin{array}{cc}{\displaystyle \frac{\partial u}{\partial t}=\mu \Delta u+u\Phi (x,t)}\hfill & (x,t)\in \Omega \times {\mathbb{R}}_{+},\\ {\displaystyle \frac{\partial u}{\partial \nu }=0,}\hfill & (x,t)\in \partial \Omega \times {\mathbb{R}}_{+}\\ u(x,0)=u_0\left(x\right),\hfill & x\in \overline{\Omega }.\end{array}\right.$$

(2)

where $\Phi (x,t)$ is locally Lipschitz on $\Omega \times {\mathbb{R}}_{+}$ and uniformly bounded. Suppose that $u(x,t)\ge 0$ and ${sup}_{t\ge 0}\left|\right|u(·,t)\left|\right|<M$ for certain $M>0$, then we have

$$\underset{t\ge 0}{sup}{\left|\right|u(·,t)\left|\right|}_{L^{\infty }(\Omega )}<M^{\ast },$$

where the bound $M^{\ast }$ depends on M and $\left|\right|u_0(·){\left|\right|}_{L^{\infty }(\Omega )}$.

Theorem 1

([27]).

(i) 

If $\phi (x,t)$$\psi (x,t)$ are nonnegative and cannot vanish identically, then system (1) possesses a unique positive solution $(u,v)$, for $x\in \overline{\Omega }$ and $t\in (0,\infty )$.

(ii) 

Any solution $(u,v)$ of system (1) satisfies the following conditions

$$\underset{t\to +\infty }{lim\; sup}u(x,t)\le K,{\int }_{\Omega }v(x,t)\mathrm{d}x\le \frac{eK(m+r)}{m}|\Omega |.$$

Moreover, denote by $K_1=max\{K,{max}_{\overline{\Omega }}\phi \left(x\right)\}$; there exists a constant $C^{\ast }>0$ such that

$${\left|\right|u(·,t)\left|\right|}_{C\left(\overline{\Omega }\right)}\le K_1{,\left|\right|v(·,t)\left|\right|}_{C\left(\overline{\Omega }\right)}\le C^{\ast }.$$

Proof. 

(i) Define

$$\left(\begin{array}{c}f(u,v)\\ g(u,v)\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{ru}{1+bv(t-\tau )}(1-\frac{u}{K})-\frac{\alpha uv}{a+u}}\\ {\displaystyle \frac{\beta uv}{a+u}-mv-\eta uv}\end{array}\right),$$

(3)

and

$${\displaystyle \overline{f}(u,v)=ru(1-\frac{u}{K})-\frac{\alpha uv}{a+u}.}$$

• Obviously, $f\le \overline{f}$ holds for any $(x,t)\in \Omega \times (0,\infty )$. Consider now the auxiliary system

  $$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle \frac{\partial u}{\partial t}=d_1\Delta u+\overline{f},}\hfill \\ {\displaystyle \frac{\partial v}{\partial t}=d_2\Delta v+g.}\hfill \end{array}\right.\end{array}$$

  (4)
• Since $u\ge 0,v\ge 0$, it is easy to see that ${\overline{f}}_v\le 0$ and $g_u\ge 0$. Now consider the following initial value problem for ODEs

  $$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle \frac{\mathrm{d}u}{\mathrm{d}t}=ru(1-\frac{u}{K}),}\hfill \\ {\displaystyle \frac{\mathrm{d}v}{\mathrm{d}t}=\frac{\beta uv}{a+u}-mv-\eta uv,}\hfill \\ u\left(0\right)={\parallel \phi (x,t)\parallel }_{\infty },v\left(0\right)={\parallel \psi (x,t)\parallel }_{\infty }.\hfill \end{array}\right.\end{array}$$

  (5)
• System (5) possesses a unique solution $\left(u\right(t;u\left(0\right),v\left(0\right)),v(t;u\left(0\right),v\left(0\right)$. Denote

  $$\underset{\overline{\Omega },t\in [-\tau ,0]}{min}\phi (x,t)={\phi }_m,\underset{\overline{\Omega },t\in [-\tau ,0]}{max}\phi (x,t)={\phi }_M,\underset{\overline{\Omega },t\in [-\tau ,0]}{min}\psi (x,t)={\psi }_m,\underset{\overline{\Omega },t\in [-\tau ,0]}{max}\psi (x,t)={\psi }_M.$$

  (6)
• It is clear that $(u(t;{\phi }_m,{\psi }_M),v(t;{\phi }_M,{\psi }_m))$ and $(u(t;{\phi }_M,{\psi }_m),v(t;{\phi }_m,{\psi }_M))$ are the lower solution and upper solution of system (4), respectively, which means that they are also the lower solution and upper solution of system (1), respectively. Therefore, system (1) possesses a unique solution satisfying

  $$u(t;{\phi }_m,{\psi }_M)\le u(x,t)\le u(t;{\phi }_M,{\psi }_m),v(t;{\phi }_m,{\psi }_M)\le v(x,t)\le v(t;{\phi }_M,{\psi }_m).$$
• Finally, the conclusion that $u(x,t),v(x,t)$ are positive for $(x,t)\in \overline{\Omega }\times (0,+\infty )$ can be obtained by using the strong maximum principle.
• (ii) Define $U\left(t\right)={\int }_{\Omega }u(x,t)\mathrm{d}x$ and $V\left(t\right)={\int }_{\Omega }v(x,t)\mathrm{d}x$. Calculating the derivatives of $U\left(t\right)$ and $V\left(t\right)$, we have

  $$\begin{array}{c}{\displaystyle \frac{\mathrm{d}U}{\mathrm{d}t}={\int }_{\Omega }{u}_t\mathrm{d}x=d_1{\int }_{\Omega }\Delta u\mathrm{d}x+{\int }_{\Omega }[\frac{ru}{1+bv}(1-\frac{u}{K})-\frac{\alpha uv}{a+u}]\mathrm{d}x={\int }_{\Omega }[\frac{ru}{1+bv}(1-\frac{u}{K})-\frac{\alpha uv}{a+u}]\mathrm{d}x,}\hfill \\ {\displaystyle \frac{\mathrm{d}V}{\mathrm{d}t}={\int }_{\Omega }{v}_t\mathrm{d}x=d_2{\int }_{\Omega }\Delta v\mathrm{d}x+{\int }_{\Omega }(\frac{\beta uv}{a+u}-mv-\eta uv)\mathrm{d}x={\int }_{\Omega }(\frac{\beta uv}{a+u}-mv-\eta uv)\mathrm{d}x,}\hfill \end{array}$$
• which leads to

  $$\begin{array}{cc}\hfill {(\frac{\beta }{\alpha }U+V)}_t& =-mV+\frac{r\beta }{\alpha }{\int }_{\Omega }\left[\frac{u}{1+bv}(1-\frac{u}{K})-\eta uv\right]\mathrm{d}x\hfill \\ \hfill & \le -m(\frac{\beta }{\alpha }U+V)+\frac{\beta }{\alpha }(m+r)U.\hfill \end{array}$$

  (7)
• Notice that the first equation of system (1) leads to the following inequality

  $$\frac{\partial u}{\partial t}-d_1\Delta u\le ru(1-\frac{u}{K})-\frac{\alpha uv}{a+u}\le r(1-\frac{u}{K}).$$
• Then, from the comparison principle for parabolic equations, we obtain

  $$\underset{t\to +\infty \overline{\Omega }}{limsupmax}u(x,t)\le K.$$

  (8)
• On the one hand, according to the maximum principle, (8) implies that ${\left|\right|u(.,t)\left|\right|}_{C\left(\overline{\Omega }\right)}\le K_1$, $\forall t\ge 0$. On the other hand, (8) also implies that $\underset{t\to +\infty \overline{\Omega }}{limsupmax}U\left(t\right)\le |\Omega |$, which means that for any $\epsilon >0$, there exists $t^{\ast }>0$ such that

  $$U\left(t\right)\le (K+\epsilon )|\Omega |,t>t^{\ast }.$$

  (9)
• Substituting (9) into (7) and integrating the obtained inequality, we obtain

  $${\int }_{\Omega }v(x,t)dx<\frac{\beta }{\alpha }U\left(t\right)+V\left(t\right)\le \frac{\frac{\beta }{\alpha }(m+r)}{m}(K+\epsilon )|\Omega |:=K_2,t>t^{\ast },$$

  (10)
• i.e., ${\left|\right|v(·,t)\left|\right|}_{L^1(\Omega )}\le =K_2.$ Thus, it can be concluded, by Lemma 1, that ${\left|\right|v(·,t)\left|\right|}_{L^{\infty }(\Omega )}\le K_3$, where $K_3$ depends on $K_2$ and ${\left|\right|\psi \left(x\right)\left|\right|}_{L^{\infty }(\Omega )}$. Consequently, ${\left|\right|v(·,t)\left|\right|}_{L^{\infty }(\Omega )}\le C^{\ast }$ holds for certain positive $C^{\ast }>0$.

□

3. Analysis of Stability and Bifurcation

In this section, we analyze the stability and Hopf bifurcation of the system at the positive equilibrium by analyzing the characteristic equation.

3.1. Stability Analysis

This section will study the equilibrium points of system (1) and their stabilities when $\tau =0$. Obviously, the constant states $E_0=(0,0)$ and $E_1=(K,0)$ are always solutions of system (1), also named the extinct and the predator-free equilibrium point, respectively. In addition, notice that any positive equilibrium point $(u^{\ast },v^{\ast })$, if it exists, must satisfy the following system of equations

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle \frac{r}{1+bv}(1-\frac{u}{K})-\frac{\alpha v}{a+u}=0,}\hfill \\ {\displaystyle \frac{\beta u}{a+u}-m-\eta u=0,}\hfill \end{array}\right.\end{array}$$

which is, clearly, equivalent to the following form

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle \alpha b{v}^2+\alpha v-r(1-\frac{u}{K})(a+u)=0,}\hfill \\ \eta u^2+(m-\beta +a\eta )u+am=0.\hfill \end{array}\right.\end{array}$$

(11)

Theorem 2.

Denote ${\Delta }_0={(m-\beta +a\eta )}^2-4\eta am$, then we have the following statements for system (1):

(i) 

When ${\Delta }_0<0$, there does not exist any positive solution of system (11). Thus, system (1) has no positive equilibrium point;

(ii) 

When ${\Delta }_0=0$ and $m-\beta +a\eta <0$, there exists a unique positive solution, named $u_1^{\ast }$, of the second equation of (11). Moreover, if $u_1^{\ast }<K$, then system (1) possesses a unique positive equilibrium;

(iii) 

When ${\Delta }_0>0$ and $m-\beta +a\eta <0$, there exist two different solutions, called $u_2^{\ast }$ and $u_3^{\ast }$, of the second equation of (11). Moreover, if $u_2^{\ast }<K$ and $u_3^{\ast }<K$, then system (1) possesses two positive equilibria.

Proof. 

The proof is simple and is omitted here. □

Denote $\Gamma \left(t\right)={(u\left(t\right),v\left(t\right))}^T={(u(·,t),v(·,t))}^T$, $D=\mathrm{diag}(d_1,d_2)$ and

$$\mathrm{dom}(\Delta )=\left\{{(u,v)}^T:u,v\in W^{2,2}(\Omega ),\frac{\partial u}{\partial \nu }=\frac{\partial u}{\partial \nu }=0\right\}.$$

Let L be the map from $C(X,[-\tau ,0\left]\right)$ to X, and for any $\phi ={({\phi }_1,{\phi }_2)}^T\in \mathcal{C}(X,[-\tau ,0])$, we have

$$L\left(\phi \right)=\left(\begin{array}{c}(\frac{r}{1+bv}-\frac{2ru}{K(1+bv)}-\frac{\alpha av}{{(a+u)}^2}){\phi }_1\left(0\right)+\frac{bru(\frac{u}{K}-1)}{{(1+bv)}^2}{\phi }_2(-\tau )-\frac{\alpha u}{a+u}{\phi }_2\left(0\right)\\ (\frac{\beta av}{{(a+u)}^2}-\eta v){\phi }_1\left(0\right)-(\frac{\beta u}{a+u}-\eta u-m){\phi }_2\left(0\right)\end{array}\right).$$

Then the linearization of system (1) at an equilibrium point can be expressed as

$$\dot{\Gamma }=D\Delta \Gamma \left(t\right)+L\left({\Gamma }_t\right),$$

(12)

whose characteristic equation can be, clearly, given by

$$\lambda y-D\Delta y-L\left(e^{\lambda }y\right)=0,$$

(13)

with $y\in \mathrm{dom}(\Delta )$, $y\ne 0$.

Note that the domain of definition of $\Delta$ is bounded and thus it has countable eigenvalues

$$0={\mu }_0<{\mu }_1<{\mu }_2<\dots <{\mu }_n<\dots \to \infty .$$

(14)

Suppose that the corresponding eigenfunctions are $Y_i\left(x\right)$, $i=1,2,\dots$ and substituting

$$y=\sum _{n=1}^{\infty }{Y}_n\left(x\right)\left(\begin{array}{c}{\beta }_n^1\hfill \\ {\beta }_n^2\hfill \end{array}\right)$$

(15)

into Equation (13), we obtain

$$\left(\begin{array}{cc}\frac{r}{1+bv}-\frac{2ru}{K(1+bv)}-\frac{\alpha av}{{(a+u)}^2}-d_1{\mu }_n& \frac{bru(\frac{u}{K}-1)}{{(1+bv)}^2}{e}^{-\lambda \tau }-\frac{u}{a+u}\\ \frac{\beta av}{{(a+u)}^2}-\eta v& \frac{\beta u}{a+u}-\eta u-m-d_2{\mu }_n\end{array}\right)\left(\begin{array}{c}{\beta }_n^1\\ {\beta }_n^2\end{array}\right)=\lambda \left(\begin{array}{c}{\beta }_n^1\\ {\beta }_n^2\end{array}\right).$$

Consequently, the characteristic Equation (13) can be transformed into the following form

$$\begin{array}{c}\hfill {\lambda }^2+A_n\lambda +B_n+C_n{e}^{-\lambda \tau }=0,n=0,1,2,\dots ,\end{array}$$

(16)

where

$$\begin{array}{cc}\hfill A_n=& (d_1+d_2){\mu }_n-(\frac{r}{1+bv}-\frac{2ru}{K(1+bv)}-\frac{\alpha av}{{(a+u)}^2})-(\frac{\beta u}{a+u}-\eta u-m),\hfill \\ \hfill B_n=& d_1{d}_2{\mu }_n^2-((\frac{\beta u}{a+u}-\eta u-m)d_1+(\frac{r}{1+bv}-\frac{2ru}{K(1+bv)}-\frac{\alpha av}{{(a+u)}^2})d_2){\mu }_n\hfill \\ \hfill & +(\frac{r}{1+bv}-\frac{2ru}{K(1+bv)}-\frac{\alpha av}{{(a+u)}^2})(\frac{\beta u}{a+u}-\eta u-m)\hfill \\ \hfill & +\frac{u}{a+u}(\frac{\beta av}{{(a+u)}^2}-\eta v),\hfill \\ \hfill C_n=& (\frac{\beta av}{{(a+u)}^2}-\eta v)\frac{bru(1-\frac{u}{K})}{{(1+bv)}^2}.\hfill \end{array}$$

(17)

For $E_0=(0,0)$, Equation (16) becomes

$$(\lambda +d_1{\mu }_n-r)(\lambda +d_2{\mu }_n+m)=0.$$

Therefore, $E_0$ is always unstable.

For $E_1=(K,0)$, Equation (16) becomes

$$(\lambda +d_1{\mu }_n+r)(\lambda +d_2{\mu }_n-\frac{K\beta }{K+a}+K\eta +m)=0.$$

Therefore, when $\frac{K\beta }{K+a}-K\eta -m<0$, $E_1$ is locally asymptotically stable.

For $E^{\ast }=(u^{\ast },v^{\ast })$, Equation (16) becomes

$${\lambda }^2+(d_1{\mu }_n+d_2{\mu }_n-a_1)\lambda +d_2{\mu }_n(d_1{\mu }_n-a_1)+b_1{a}_3+b_1{a}_2{e}^{-\lambda \tau }=0,$$

(18)

where

$$\begin{array}{cc}\hfill & a_1=\frac{r}{1+b{v}^{\ast }}-\frac{2r{u}^{\ast }}{(1+b{v}^{\ast })K}-\frac{\alpha a{v}^{\ast }}{{(a+u^{\ast })}^2},a_2=\frac{br{u}^{\ast }(1-\frac{u^{\ast }}{K})}{{(1+b{v}^{\ast })}^2},\hfill \\ \hfill & a_3=\frac{\alpha u^{\ast }}{a+u^{\ast }},b_1=\frac{a\beta v^{\ast }}{{(a+u^{\ast })}^2}-\eta v^{\ast }.\hfill \end{array}$$

(19)

Obviously, we have the following result for the positive equilibrium $E^{\ast }=(u^{\ast },v^{\ast })$.

Theorem 3.

Assume $\tau =0$. If the conditions $a_1<0$ and $b_1(a_3+a_2)>0$ hold, then $E^{\ast }=(u^{\ast },v^{\ast })$ is locally asymptotically stable.

3.2. Hopf Bifurcation

An interesting question is how the stability changes as the delay $\tau$ varies. In this section, we answer this question. Note that the characteristic Equation (18) can be written into

$${\lambda }^2+B_{1n}\lambda +B_{2n}+B_3{e}^{-\lambda \tau }=0,$$

(20)

where

$$\begin{array}{ccc}{B}_{1n}\hfill & =\hfill & (d_1{\mu }_n+d_2{\mu }_n-a_1),\hfill \\ B_{2n}\hfill & =\hfill & d_2{\mu }_n(d_1{\mu }_n-a_1)+b_1{a}_3,\hfill \\ B_3\hfill & =\hfill & b_1{a}_2.\hfill \end{array}$$

(21)

Substitute $\lambda =i\omega$ into (20) and a simple calculation gives

$$\begin{array}{ccc}{\omega }^2-B_{2n}\hfill & =\hfill & B_{3}cos\left(\omega \tau \right),\hfill \\ B_{1n}\omega \hfill & =\hfill & B_{3}sin\left(\omega \tau \right),\hfill \end{array}$$

(22)

which leads to

$${\omega }^4+P_n{\omega }^2+Q_n=0,$$

(23)

where

$$P_n=B_{1n}^2-2B_{2n},Q_n=B_{2n}^2-B_3^2.$$

Denote ${\Delta }_0=P_n^2-4Q_n^2$, then solving Equation (23) gives ${\omega }^2\left(n\right)=\frac{-P_n\pm \sqrt{{\Delta }_0}}{2}$. Note that, for any solution $\omega \left(n\right)$ of Equation (23), substituting it into (22) and solving for $\tau$, we will obtain the corresponding delay values. Thus, the following results can be obtained easily by analyzing Equations (22) and (23).

Theorem 4.

Assume that $E^{\ast }=(u^{\ast },v^{\ast })$ of system (1) exists, then we have:

(i) 

If ${\Delta }_0<0$, or ${\Delta }_0>0,Q_n\ge 0,P_n\ge 0$, or ${\Delta }_0=0,P_n\ge 0$ holds, then there does not exist any positive root for Equation (23);

(ii) 

If $Q_n<0$, or $P_n=-2\sqrt{Q_n},Q_n>0$, or $Q_n=0,P_n<0$ holds, then there exists only one positive real root, denoted by ${\omega }_1^{+}\left(n\right)$ for Equation (23) and consequently, Equation (16) has purely imaginary roots $\pm i{\omega }_1^{+}\left(n\right)$ while the delay $\tau ={\tau }_{1,j}^{+}\left(n\right)$ is given by

$${\tau }_{1,j}^{+}\left(n\right)=\frac{1}{{\omega }_1^{+}\left(n\right)}\left\{arccos\left[\frac{C_n({\omega }_1^{{+}^2}\left(n\right)-B_n)-{\omega }_1^{{+}^2}\left(n\right)A_n}{\alpha {\omega }_1^{{+}^2}\left(n\right)+\alpha C_n^2}\right]+2j\pi \right\},j=0,1,2,\dots$$

(24)

which can be obtained by solving (22).

(iii) 

If $P_n<0$, $Q_n>0$, ${\Delta }_0>0$, then there exist two positive real roots, denoted by ${\omega }_2^{\pm }\left(n\right)$ for Equation (23) and consequently Equation (16) has purely imaginary roots $\pm i{\omega }_2^{\pm }\left(n\right)$, when $\tau ={\tau }_{2,j}^{\pm }\left(n\right)$, while the delay $\tau ={\tau }_{2,j}^{\pm }\left(n\right)$ is given by

$${\tau }_{2,j}^{\pm }\left(n\right)=\frac{1}{{\omega }_2^{+}\left(n\right)}\left\{arccos\left[\frac{C_n({\omega }_2^{{\pm }^2}\left(n\right)-B_n)-{\omega }_2^{{\pm }^2}\left(n\right)A_n}{\alpha {\omega }_2^{{\pm }^2}\left(n\right)+\alpha C_n^2}\right]+2j\pi \right\},j=0,1,2,\dots$$

(25)

Lemma 2.

Let $\lambda \left(\tau \right)=\xi \left(\tau \right)\pm i\omega \left(\tau \right)$ be the root of Equation (16) which satisfies the conditions

$$\begin{array}{cc}\hfill & \xi \left({\tau }_{1,j}^{+}\left(n\right)\right)=0(\mathrm{resp}.\xi \left({\tau }_{2,j}^{\pm }\left(n\right)\right)=0),\hfill \\ \hfill & \omega \left({\tau }_{1,j}^{+}\left(n\right)\right)={\omega }_1^{+}\left(n\right)(\mathrm{resp}.\omega \left({\tau }_{2,j}^{\pm }\left(n\right)\right)={\omega }_2^{\pm }\left(n\right)).\hfill \end{array}$$

Then the following transversality conditions hold

$$\begin{array}{cc}\hfill & {\left.\frac{\mathrm{d}\left(\mathrm{Re}\lambda \right(\tau \left)\right)}{\mathrm{d}\tau }\right|}_{\tau ={\tau }_{1,j}^{+}\left(n\right)}>0,{\left.\frac{\mathrm{d}\left(\mathrm{Re}\lambda \right(\tau \left)\right)}{\mathrm{d}\tau }\right|}_{\tau ={\tau }_{2,j}^{+}\left(n\right)}>0,{\left.\frac{\mathrm{d}\left(\mathrm{Re}\lambda \right(\tau \left)\right)}{\mathrm{d}\tau }\right|}_{\tau ={\tau }_{2,j}^{-}\left(n\right)}<0.\hfill \end{array}$$

(26)

Proof. 

Differentiating Equation (18) with respect to $\lambda$, we have

$$2\lambda +A_n+\alpha e^{-\lambda \tau }-\alpha (\lambda +C_n)(\tau +\lambda \frac{\mathrm{d}\tau }{\mathrm{d}\lambda })e^{-\lambda \tau }=0,$$

which gives

$$\begin{array}{cc}\hfill & {\left[\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\right]}^{-1}=\frac{(2\lambda +A_n)e^{\lambda \tau }}{\lambda (\alpha \lambda +\alpha C_n)}+\frac{1}{\lambda (\lambda +C_n)}-\frac{\tau }{\lambda }.\hfill \end{array}$$

(27)

Combining with Equation (27), we have

$$\begin{array}{cc}\hfill & {\mathrm{Re}{\left[\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\right]}^{-1}}_{\lambda =\pm i{\omega }_2^{\pm }\left(n\right)}\hfill \\ \hfill & ={\left[\frac{2{\omega }^2+A_n^2-2B_n-{\alpha }^2}{{\alpha }^2{C}_n^2+{\alpha }^2{\omega }^2}\right]}_{\lambda =\pm i{\omega }_2^{\pm }\left(n\right)}={\left[\frac{2{\omega }^2+P_n}{{\alpha }^2{C}_n^2+{\alpha }^2{\omega }^2}\right]}_{\lambda =\pm i{\omega }_2^{\pm }\left(n\right)}=\pm \frac{\sqrt{{\Delta }_0}}{{\alpha }^2{C}_n^2+{\alpha }^2{\omega }_2^{{\pm }^2}\left(n\right)}.\hfill \end{array}$$

Substituting the value of $\tau$ and $\omega$, we have

$$\begin{array}{cc}\hfill & {\left[\frac{\mathrm{d}\left(\mathrm{Re}\lambda \right(\tau \left)\right)}{\mathrm{d}\tau }\right]}_{\begin{array}{c}\tau ={\tau }_{1,j}^{+}\left(n\right),\hfill \\ \omega ={\omega }_1^{+}\left(n\right)\hfill \end{array}}^{-1}=\frac{\sqrt{{\Delta }_0}}{{\alpha }^2{C}_n^2+{\alpha }^2{\omega }_1^{{+}^2}\left(n\right)}>0,\hfill \\ \hfill & {\left[\frac{\mathrm{d}\left(\mathrm{Re}\lambda \right(\tau \left)\right)}{\mathrm{d}\tau }\right]}_{\begin{array}{c}\tau ={\tau }_{2,j}^{+}\left(n\right),\hfill \\ \omega ={\omega }_2^{+}\left(n\right)\hfill \end{array}}^{-1}=\frac{\sqrt{{\Delta }_0}}{{\alpha }^2{C}_n^2+{\alpha }^2{\omega }_2^{{+}^2}\left(n\right)}>0,\hfill \end{array}$$

$$\begin{array}{cc}\hfill & {\left[\frac{\mathrm{d}\left(\mathrm{Re}\lambda \right(\tau \left)\right)}{\mathrm{d}\tau }\right]}_{\begin{array}{c}\tau ={\tau }_{2,j}^{-}\left(n\right),\hfill \\ \omega ={\omega }_2^{-}\left(n\right)\hfill \end{array}}^{-1}=-\frac{\sqrt{{\Delta }_0}}{{\alpha }^2{C}_n^2+{\alpha }^2{\omega }_2^{{-}^2}\left(n\right)}<0.\hfill \end{array}$$

□

Denote

$$\begin{array}{cc}\hfill & \Gamma =\left\{n\in {\mathbb{N}}_0|Q_n<0\mathrm{or}{P}_n=-2\sqrt{Q_n},Q_n>0\mathrm{or}{Q}_n=0,P_n<0\right\}.\hfill \end{array}$$

Note that for fixed n, we have ${\tau }_{1,j_1}^{+}\left(n\right)\le {\tau }_{1,j_2}^{+}\left(n\right)$, for $j_1\le j_2$, and thus ${\tau }_{1,0}^{+}\left(n\right)=\underset{j\in {\mathbb{N}}_0}{min}{\tau }_{1,j}^{+}\left(n\right)$. Define the smallest critical value of delay by

$${\tau }^{\ast }=\underset{n\in \Gamma }{min}\left\{{\tau }_{1,0}^{+}\left(n\right)\right\}.$$

(28)

Recall that $E^{\ast }$ denotes the positive equilibrium point of system (1). The following result can be concluded based on the above discussions.

Theorem 5.

Assume that $a_1<0,b_1(a_3+a_2)>0$ are satisfied, where $a_1,a_2,a_3,b_1$ are defined by (19), then the following statements hold true for system (1):

(i) 

If ${\Delta }_0<0$, or ${\Delta }_0>0,Q_n\ge 0,P_n\ge 0$, or ${\Delta }_0=0,P_n\ge 0$ hold, then, for any delay $\tau \ge 0$, $E^{\ast }$ is always locally asymptotically stable.

(ii) 

If for certain $n\in {\mathbb{N}}_0$, $Q_n<0$, or $P_n=-2\sqrt{Q_n},Q_n>0$ or $Q_n=0,P_n<0$ hold, then when $0\le \tau <{\tau }^{\ast })$, $E^{\ast }$ is locally asymptotically stable, and when $\tau >{\tau }^{\ast }$, $E^{\ast }$ is unstable. Moreover, a Hopf bifurcation occurs at $E^{\ast }$ when $\tau ={\tau }_{1,j}{\left(n\right)}^{+}$ with $j\in {\mathbb{N}}_0$ and $n\in \Gamma$;

(iii) 

If $P_n<0$, $Q_n>0$, ${\Delta }_0>0$ hold, then the stability of $E^{\ast }$ will switch limit times from stable to unstable, to stable, and so on.

Remark 1.

If Equation (28) holds for $n=0$, then the system undergoes spatially homogeneous Hopf bifurcation. If Equation (28) holds for $n>0$, then the system undergoes spatially inhomogeneous Hopf bifurcation and the ODE system cannot occur.

4. Non-Constant Positive Solutions

In this section, with the maximum principle and theory of the Leray–Schauder degree, the nonexistence and existence of nonconstant positive solutions will be studied, and some conditions will be derived. For this purpose, consider the following elliptic system

$$\begin{array}{c}\hfill \left\{\begin{array}{c}-d_1\Delta u=\frac{ru}{1+bv}(1-\frac{u}{K})-\frac{\alpha uv}{a+u}\hfill \\ -d_2\Delta v=\frac{\beta uv}{a+u}-mv-\eta uv.\hfill \end{array}\right.\end{array}$$

(29)

Note that from now on, all the solutions are assumed to be defined in $C^2(\Omega )\cap C^1\left(\overline{\Omega }\right)$.

4.1. A Priori Estimate of the Positive Solution

Recall two useful lemmas [28].

Lemma 3.

Assume that $\varphi \in C^2\left(\overline{\Omega }\right)$ and ${\partial }_v\varphi =0$ on $\partial \overline{\Omega }$. If $\varphi \left(x_0\right)={\max }_{\overline{\Omega }}\varphi$, then we have $-\Delta \varphi \left(x_0\right)\ge 0$.

Lemma 4.

Consider the equation $\Delta \varphi +c\left(x\right)\varphi =0$ with a homogeneous Neumann boundary condition, where $c\left(x\right)\in C\left(\overline{\Omega }\right)$. Assume that $\varphi \in C^2(\Omega )\cap C^1\left(\overline{\Omega }\right)$ is a positive solution, then there exists a constant $C_{\ast }=C_{\ast }{(\parallel c\parallel }_{\infty })>0$ such that

$$\underset{\overline{\Omega }}{max}\varphi ⩽C_{\ast }\underset{\overline{\Omega }}{min}\varphi .$$

The following theorem describes the upper and lower bound of the positive solutions for system (29).

Theorem 6.

Any positive solution $(u,v)$ of system (29) satisfies the following conditions

$$0<\underset{\overline{\Omega }}{max}u\left(x\right)\le K,0<\underset{\overline{\Omega }}{max}v\left(x\right)\le \frac{K{(r+m{d}_1)}^2}{4\alpha d_2^{2}rm}.$$

(30)

Proof. 

Assume that $\left(u\right(x),v(x\left)\right)$ is a solution of (1) satisfying $u\left(x\right)\ge 0$ and $v\left(x\right)\ge 0$. Notice that if for some $x_0\in \Omega$, we have $v\left(x_0\right)=0$, then it can be concluded, by the strong maximum principle, that $v\left(x\right)\equiv 0$, and moreover $u\left(x\right)$ satisfies the following equations

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}-d_1\Delta u=ru(1-\frac{u}{K}),\hfill & x\in \Omega ,\\ \frac{\partial u}{\partial n}=0,\hfill & x\in \partial \Omega .\end{array}\right.\end{array}$$

In the same way, if $u\left(x_0\right)=0$, then $u\left(x\right)\equiv 0$ and consequently, from (29), we obtain that $v\left(x\right)\equiv 0$. Otherwise, for $x\in \overline{\Omega }$, we have $u\left(x\right)>0$ and $v\left(x\right)>0$.

By Lemma 3, we have $u\left(x\right)\le K$ and thus, from the strong maximum principle, it can be concluded that $u\left(x\right)<K$, $x\in \overline{\Omega }$. A direct calculation gives

$$\begin{array}{cc}\hfill -(\frac{\beta }{\alpha }{d}_1\Delta u+d_2\Delta v)& =\frac{\beta }{\alpha }\frac{ru}{1+bv}(1-\frac{u}{K})-mv-\eta uv\hfill \\ \hfill & \le (\frac{\beta r}{\alpha }+\frac{m\beta d_1}{\alpha d_2})u-\frac{\beta r}{\alpha K}{u}^2-\frac{m}{d_2}(\frac{\beta d_1}{\alpha }u+d_{2}v)\hfill \\ \hfill & \le \frac{K{(r+m{d}_1)}^2}{4\alpha d_2^{2}r}-\frac{m}{d_2}(\frac{\beta d_1}{\alpha }u+d_{2}v),\hfill \end{array}$$

Then, by Lemma 3, we have

$$\frac{\beta d_1}{\alpha }u+d_{2}v\le \frac{K{(r+m{d}_1)}^2}{4\alpha d_{2}rm}.$$

Hence, for any $x\in \Omega$, we have $v\left(x\right)<\frac{K{(r+m{d}_1)}^2}{4\alpha d_2^{2}rm}:=V_K$. □

Theorem 7.

Suppose that $K\ne \frac{(\beta u-a\eta -m)\pm \sqrt{{(\beta u-a\eta -m)}^2-4\eta ma}}{2\eta }$ and the equilibrium point $E^{\ast }=(u^{\ast },v^{\ast })$ exists. Then there exists a constant $\underline{C}>0$ such that for any solution $\left(u\right(x),v(x\left)\right)$ of system (29), we always have

$$\underset{\overline{\Omega }}{min}u\left(x\right)\ge \underline{C},\underset{\overline{\Omega }}{min}v\left(x\right)\ge \underline{C}.$$

(31)

Proof. 

Denote

$$c_1\left(x\right)=d_1^{-1}(\frac{r}{1+bv}r(1-\frac{u}{K})-\frac{\alpha v}{a+u}),c_2\left(x\right)=d_2^{-1}(\frac{\beta u}{a+u}-m-\eta u).$$

Then, from (30), there exists a constant $\widehat{C}>0$ such that $\left|\right|c_1{\left|\right|}_{\infty }$, $\left|\right|c_2{\left|\right|}_{\infty }\le \widehat{C}$. Moreover, by Lemma 4, there exists a constant $C_2>0$ such that

$$\underset{\overline{\Omega }}{max}u\left(x\right)\le C_2\underset{\overline{\Omega }}{min}u\left(x\right),\underset{\overline{\Omega }}{max}v\left(x\right)\le C_2\underset{\overline{\Omega }}{min}v\left(x\right).$$

Now, we will prove (31) by contradiction. Assume that (31) is not true. There then exists a sequence of solution $\left\{(u_n,v_n)\right\}$ of system (29) satisfying

$$\underset{\overline{\Omega }}{max}{u}_n\to 0\mathrm{or}\underset{\overline{\Omega }}{max}{v}_n\to 0,$$

(32)

which implies that there exists a subsequence of $\left\{(u_n,v_n)\right\}$ converging to $(u_0,v_0)$ in $C^2\left(\overline{\Omega }\right)$. For simplicity, we still denote the subsequence by $\left\{(u_n,v_n)\right\}$. Clearly, we have $u_0\le K$ and moreover, from condition (32), we have either $u_0\equiv 0$ or $v_0\equiv 0$.

Note that integrating the second equation of (29) for $u_n$ and $v_n$, we obtain

$${\int }_{\Omega }{v}_n(\frac{\beta u_n}{a+u_n}-m-\eta u_n)\mathrm{d}x=0.$$

(33)

If $u_0\equiv 0$, then no matter whether $v_0\neg \equiv 0$ or $v_0\equiv 0$, we have $v_n>0$ uniformly. Therefore, when n is large enough, we have

$$v_n(\frac{\beta u_n}{a+u_n}-m-\eta u_n)\to -m{v}_n<0,$$

which contradicts Equation (33).

If $u_0\neg \equiv 0$, then $v_0\equiv 0$. Moreover, it is easy to see that $u_0$ satisfies Equation (29) which implies that $u_n\to u_0\equiv K$ as $n\to \infty$ and therefore we have

$$\frac{\beta u_n}{a+u_n}-m-\eta u_n\to \frac{\beta K}{a+K}-m-\eta K.$$

Since $K\ne \frac{(\beta u-a\eta -m)\pm \sqrt{{(\beta u-a\eta -m)}^2-4\eta ma}}{2\eta }$ and the positive equilibrium point $(u^{\ast },v^{\ast })$ exists, we have then

$$v_n(\frac{\beta u_n}{a+u_n}-m-\eta u_n)\to v_n(\frac{\beta K}{a+K}-m-\eta K),$$

which contradicts Equation (33) and completes the proof. □

4.2. Non-Existence of Non-Constant Positive Solutions

Theorem 8.

(i) 

If $\frac{\beta K}{a+K}-m-\eta \underline{C}<0$, then system (29) possesses only one nonzero solution $E_1=(K,0)$;

(ii) 

Given any system parameters r, K, α, β, d, m, and η, there exists always a constant $d^{\ast }>0$ such that if $d_1>d^{\ast }$ and $d_2>d^{\ast }$, then system (29) possesses the unique positive solution $E^{\ast }=(u^{\ast },v^{\ast })$.

Proof. 

(i) Assume that there exists a nonconstant positive solution $\left(u\right(x),v(x\left)\right)$ for system (29), then by Theorem 6, we have $u\left(x\right)\le K$, for $x\in \Omega$. Applying the condition $\frac{\beta K}{a+K}-m-\eta \underline{C}$$<0$, we have then

$$-d_2\Delta v=v(-m+\frac{\beta u}{a+u}-\eta u)\le v(\frac{\beta K}{a+K}-m-\eta \underline{C})\le 0,x\in \Omega$$

which leads to $v\equiv 0$ in $\Omega$ by Lemma 3. This gives a contradiction.

(ii) Assume that $\left(u\right(x),v(x\left)\right)$ is any nonnegative solution of (29) and define

$$\begin{array}{cc}\hfill & \overline{u}={|\Omega |}^{-1}{\int }_{\Omega }u\left(x\right)\mathrm{d}x\ge 0,\overline{v}={|\Omega |}^{-1}{\int }_{\Omega }v\left(x\right)\mathrm{d}x\ge 0.\hfill \end{array}$$

Multiplying $u-\overline{u}$ on both sides of the first equation of (29), we obtain

$$\begin{array}{cc}\hfill d_1{\int }_{\Omega }{|\nabla (u-\overline{u})|}^2\mathrm{d}x=& {\int }_{\Omega }(u-\overline{u})\left(\frac{ru}{1+bv}(1-\frac{u}{K})-\frac{r\overline{u}}{1+b\overline{v}}(1-\frac{\overline{u}}{K})-\frac{\alpha uv}{a+u}+\frac{\alpha \overline{u}\overline{v}}{a+\overline{u}}\right)\mathrm{d}x\hfill \\ \hfill =& {\int }_{\Omega }(\frac{r}{1+bv}+\frac{\alpha a\overline{v}}{(a+u)(a+\overline{u})}){(u-\overline{u})}^2\mathrm{d}x-{\int }_{\Omega }(\frac{\alpha u}{a+u}\hfill \\ \hfill & +\frac{rb\overline{u}}{(1+bv)(1+b\overline{v})})(u-\overline{u})(v-\overline{v})\mathrm{d}x\hfill \\ \hfill & +{\int }_{\Omega }(\frac{\beta u}{a+u}-\eta \overline{u}-m){(v-\overline{v})}^2\mathrm{d}x.\hfill \end{array}$$

(34)

In the same way, we obtain

$$\begin{array}{cc}\hfill d_2{\int }_{\Omega }{|\nabla (v-\overline{v})|}^2\mathrm{d}x=& {\int }_{\Omega }{(v-\overline{v})}^2(\frac{\beta u}{a+u}-\eta \overline{u}-m)\overline{v})\mathrm{d}x\hfill \\ \hfill & +{\int }_{\Omega }(\frac{\beta a\overline{v}}{(a+u)(a+\overline{u})}-\eta v)(u-\overline{u})(v-\overline{v})\mathrm{d}x.\hfill \end{array}$$

(35)

Let $A,B$ be two constants defined by

$$A=r+\frac{2\alpha +\beta }{2a}\overline{C}+\frac{\alpha }{2}+\frac{rb}{2}\overline{C},B=\frac{\beta }{a}-\eta \underline{C}-m+\frac{\beta }{2a}+\frac{\alpha }{2}+\frac{rb}{2}\overline{C}.$$

Then by using the $Poincar\acute{e}$ inequality and (34) and (35), a direct calculation gives

$$\begin{array}{cc}\hfill & d_1{\mu }_1{\int }_{\Omega }{\left(u-\overline{u}\right)}^2\mathrm{d}x+d_2{\mu }_1{\int }_{\Omega }{\left(v-\overline{v}\right)}^2\mathrm{d}x\hfill \\ & \le d_1{\int }_{\Omega }|\nabla (u-\overline{u}){|}^2\mathrm{d}x+d_2{\int }_{\Omega }{|\nabla (v-\overline{v})|}^2\mathrm{d}x\hfill \\ & \le (r+\frac{2\alpha +\beta }{2a}\overline{C}+\frac{\alpha }{2}+\frac{rb}{2}\overline{C}){\int }_{\Omega }{(u-\overline{u})}^2\mathrm{d}x\hfill \\ \hfill & +(\frac{\beta }{a}-\eta \underline{C}-m+\frac{\beta }{2a}+\frac{\alpha }{2}+\frac{rb}{2}\overline{C}){\int }_{\Omega }{(v-\overline{v})}^2\mathrm{d}x\hfill \\ \hfill & =A{\int }_{\Omega }{(u-\overline{u})}^2\mathrm{d}x+B{\int }_{\Omega }{(v-\overline{v})}^2\mathrm{d}x,\hfill \end{array}$$

(36)

Denote $d^{\ast }=\frac{1}{{\mu }_1}max\{A,B\}$, then (36) implies that if $min\{d_1,d_2\}>d^{\ast }$, then we have $\nabla (u-\overline{u})=\nabla (v-\overline{v})=0$, which implies that $\left(u\right(x),v(x\left)\right)$ is a constant solution. □

4.3. Existence of Non-Constant Positive Solution

Let $\mathbf{w}={(u,v)}^T$ and $\mathbf{W}=\{\mathbf{w}\in {\left[C^2\left(\overline{\Omega }\right)\right]}^2|{\partial }_n\mathbf{w}=0,\mathrm{on}\partial \Omega \}$. Denote

$$\begin{array}{cc}\hfill & {\mathbf{W}}^{+}=\{\mathbf{w}\in \mathbf{W}|u>0,v>0,\mathrm{on}\overline{\Omega }\}.\hfill \\ \hfill & \mathbf{F}\left(\mathbf{w}\right)={\left(\frac{ru}{1+bv}(1-\frac{u}{K})-\frac{\alpha uv}{a+u},\frac{\beta uv}{a+u}-mv-\eta uv\right)}^T.\hfill \end{array}$$

Then Equation (29) can be expressed as

$$-D\Delta \mathbf{w}=\mathbf{F}\left(\mathbf{w}\right)\mathrm{in}\frac{\partial \mathbf{w}}{\partial n}=0\mathrm{on}\Omega ,$$

(37)

with $D=\mathrm{diag}(d_1,d_2)$. Let I be the identity operator and clearly Equation (29) is equivalent to

$$\mathcal{F}\left(\mathbf{w}\right)=\mathbf{w}-{(I-\Delta )}^{-1}(D^{-1}\mathbf{F}\left(\mathbf{w}\right)+\mathbf{w})=0.$$

(38)

Therefore, for any solution ${\mathbf{w}}^{\ast }=(u^{\ast },v^{\ast })$, we have

$${\mathcal{F}}_{\mathbf{w}}\left({\mathbf{w}}^{\ast }\right)=I-{(I-\Delta )}^{-1}(D^{-1}{\mathbf{F}}_{\mathbf{w}}\left({\mathbf{w}}^{\ast }\right)+I).$$

Recall that ${\mu }_j$, for $j>0$, and eigenvalues of $\Delta$ are defined in (14). Define the matrix $B_j$ as

$$B_j=I-\frac{1}{1+{\mu }_j}[D^{-1}{\mathbf{F}}_{\mathbf{w}}\left({\mathbf{w}}^{\ast }\right)+I]=\frac{1}{1+{\mu }_j}[{\mu }_{j}I-D^{-1}{\mathbf{F}}_{\mathbf{w}}\left({\mathbf{w}}^{\ast }\right)].$$

Notice that $\mathcal{F}$ and $B_j$, for any $j>0$, have the same eigenvalues. When every eigenvalues of ${\mathcal{F}}_{\mathbf{w}}\left({\mathbf{w}}^{\ast }\right)$ is nonzero, let $\xi$ denote the number of its negative eigenvalues and define

$$\mathrm{index}(\mathcal{F}(·),{\mathbf{w}}^{\ast })={(-1)}^{\xi }.$$

Denote

$$\begin{array}{cc}\hfill & H(d_1,d_2;\mu )=\det [\mu I-D^{-1}{\mathbf{F}}_{\mathbf{w}}\left({\mathbf{w}}^{\ast }\right)]=\frac{1}{d_1{d}_2}\det [\mu D-{\mathbf{F}}_{\mathbf{w}}\left({\mathbf{w}}^{\ast }\right)],\hfill \\ \hfill & ={\mu }^2-\frac{B_1(u^{\ast },v^{\ast })}{d_1}\mu -\frac{B_2(u^{\ast },v^{\ast })B_3(u^{\ast },v^{\ast })}{d_1{d}_2},\hfill \end{array}$$

(39)

where

$$\begin{array}{cc}\hfill & B_1(u^{\ast },v^{\ast })=\frac{r}{1+b{v}^{\ast }}-\frac{2r{u}^{\ast }}{(1+b{v}^{\ast })K}-\frac{\alpha a{v}^{\ast }}{{(a+u^{\ast })}^2},\hfill \\ \hfill & B_2(u^{\ast },v^{\ast })=-\frac{br{u}^{\ast }(1-\frac{u^{\ast }}{K})}{{(1+b{v}^{\ast })}^2}-\frac{\alpha u^{\ast }}{a+u^{\ast }},\hfill \\ \hfill & B_3(u^{\ast },v^{\ast })=\frac{a\beta v^{\ast }}{{(a+u^{\ast })}^2}-\eta v^{\ast }.\hfill \end{array}$$

Obviously, if $d_2{B}_1^2(u^{\ast },v^{\ast })+4d_1{B}_2(u^{\ast },v^{\ast })B_3(u^{\ast },v^{\ast })>0$ and $B_1(u^{\ast },v^{\ast })>0$, then there exist two positive solutions for $H(d_1,d_2;\mu )=0$ as follows

$$\begin{array}{cc}\hfill & {\mu }^{-}=\frac{d_2{B}_1(u^{\ast },v^{\ast })-\sqrt{d_2^2{B}_1^2(u^{\ast },v^{\ast })+4d_1{d}_2{B}_2(u^{\ast },v^{\ast })B_3(u^{\ast },v^{\ast }))}}{2d_1{d}_2},\hfill \\ \hfill & {\mu }^{+}=\frac{d_2{B}_1(u^{\ast },v^{\ast })+\sqrt{d_2^2{B}_1^2(u^{\ast },v^{\ast })+4d_1{d}_2{B}_2(u^{\ast },v^{\ast })B_3(u^{\ast },v^{\ast }))}}{2d_1{d}_2}.\hfill \end{array}$$

A direct calculation gives

$$\underset{d_2\to \infty }{lim}{\mu }^{-}(d_1,d_2)=0,\underset{d_2\to \infty }{lim}{\mu }^{+}(d_1,d_2)=\frac{B_1(u^{\ast },v^{\ast })}{d_1}.$$

Lemma 5.

If $H(d_1,d_2;{\mu }_i)\ne 0$$(i\ge 0)$, then $\mathrm{index}(\mathcal{F}(·),{\mathbf{w}}^{\ast })={(-1)}^{\xi }$, where

$$\xi =\sum _{H(d_i,d_2;{\mu }_i)<0}m\left({\mu }_i\right),$$

with $m\left({\mu }_i\right)$ being the algebraic multiplicity of ${\mu }_i$.

Theorem 9.

Assume that $m>m_0$ and

$$\frac{d_2}{d_1}>\frac{-4B_2(u^{\ast },v^{\ast })B_3(u^{\ast },v^{\ast })}{B_1^2(u^{\ast },v^{\ast })}.$$

Furthermore, if there exists $i,j\in N$, such that $0\le {\mu }_j<{\mu }^{-}<{\mu }_{j+1}\le {\mu }_i<{\mu }^{+}<{\mu }_{i+1}$ and $\sum _{k=j+1}^{i}m\left({\mu }_k\right)$ is odd, then Equation (29) has at least one non-constant positive solution.

Proof. 

Recall that $d^{\ast }$ is a constant defined in Theorem 8. Consider a mapping $\mathcal{A}:[0,1]$$\to {\mathbf{W}}^{+}$ by

$$\mathcal{A}(\mathbf{w},t)\triangleq {(-\Delta +I)}^{-1}\left(\begin{array}{c}u+\left(\frac{1-t}{d^{\ast }}+\frac{t}{d_1}\right)f(u,v)\hfill \\ v+\left(\frac{1-t}{d^{\ast }}+\frac{t}{d_2}\right)g(u,v)\hfill \end{array}\right),$$

which indicates that $E^{\ast }$ is the unique fixed point of $\mathcal{A}(·,0)$. On the one hand, it is easy to see that

$$\mathrm{deg}(I-\mathcal{A}(·,0),\Lambda ,0)=\mathrm{index}(I-\mathcal{A}(·,0),E^{\ast })=1.$$

(40)

Notice that any solution of (29) is, in fact, a fixed point of $\mathcal{A}(·,1)$. Denote $\mathcal{F}=I-\mathcal{A}(·,1)$, and then if ${\mathbf{w}}^{\ast }$ is the unique solution of (29), then we have

$$\mathrm{deg}(I-\mathcal{A}(·,1),\Lambda ,0)=\mathrm{index}(\mathcal{F},{\mathbf{w}}^{\ast })={(-1)}^{\sum _{k=j+1}^{i}m\left({\mu }_k\right)}=-1.$$

(41)

□

5. Numerical Simulation

5.1. The Stability of the Equilibrium $E_1=(K,0)$

Let $\Omega =(0,\pi )$ and the parameters of system (1) be chosen as $r=0.5,K=4,$$b=0.4,m=0.25,\alpha =0.5$, $\beta =1$, $a=0.5$, $\eta =0.5$. Obviously, the parameters satisfy the condition $\frac{K\beta }{K+a}-K\eta -m<0$. Therefore, $E_1=(K,0)$ is locally asymptotically stable (see Figure 1).

5.2. Hopf Bifurcation and Chaos Induced by Delay

The parameters are chosen as $r=0.9,b=10,K=11,a=5,\alpha =0.9,\beta =0.3,$$m=0.1,\eta =0.01$. Then the positive equilibrium points can be obtained as $E_1^{\ast }=(5,0.6902)$ and $E_2^{\ast }=(10,0.3226)$, and $E_2$ is always unstable. By a direct calculation, we know that the parameters satisfy the condition of Theorem 3. Therefore, when $\tau =0$, $E_1^{\ast }$ is locally asymptotically stable (Figure 2). In addition, we find that $E_1=(K,0)$ is also locally asymptotically stable. Figure 3 shows that $E_1^{\ast }$ and $E_1=(11,0)$ are both locally asymptotically stable, which means that under some initial condition predators may be extinct.

According to Equations (24) and (28), we can obtain that the critical value ${\tau }^{\ast }=17.3$. Thus, according to Theorem 5, $E_1^{\ast }$ is locally asymptotically stable (Figure 4) when $\tau \in [0,17.3)$ and is unstable when $\tau >17.3$.

With the further increase in delay, we find that delay can even lead to chaos phenomena. Figure 5 and Figure 6 show that the system switches between period and chaos, where we use spatially averaged population density.

Now, stipulate the parameters as $r=0.9,b=50,K=2,a=0.5,\alpha =0.6,$$\beta =0.3,m=0.2,\eta =0.01$ and in this case that system (1) possesses a unique positive equilibrium $E^{\ast }=(1.2056,0.1329)$. Solving Equation (23), we obtain two positive solutions ${\omega }_1=0.0526$ and ${\omega }_2=0.0220$. Moreover, the critical values of $\tau$ can be obtained, from (25), as

$$\begin{array}{ccc}{\tau }_1^j& =& \hfill 25.9,145.5,265.0,384.6,\dots ,\\ {\tau }_2^j& =& \hfill 123.4,408.4,693.5,978.6,\dots .\end{array}$$

In addition, from (26), we obtain

$${\left.\frac{\mathrm{d}\mathrm{Re}\lambda \left(\tau \right)}{\mathrm{d}\tau }\right|}_{\tau ={\tau }_1^{\left(j\right)},\lambda =i{\omega }_1}=147.1872>0,{\left.\frac{\mathrm{d}\mathrm{Re}\lambda \left(\tau \right)}{\mathrm{d}\tau }\right|}_{\tau ={\tau }_2^{\left(j\right)},\lambda =i{\omega }_2}=-14.3336<0,$$

which implies that eigenvalues cross the imaginary axis from left to right at $\tau ={\tau }_1^j$, and from right to left at $\tau ={\tau }_2^j$. Figure 7 shows the simulations for different values of $\tau$, which indicates that when $\tau \in [0,{\tau }_1^1)\bigcup ({\tau }_2^1,{\tau }_1^2)$, $E^{\ast }$ is asymptotically stable, and when $\tau \in ({\tau }_1^1,{\tau }_2^1)\bigcup ({\tau }_1^2,+\infty )$, $E^{\ast }$ is unstable.

5.3. The Effect of the Cost of Fear

The parameters are chosen as $r=0.9,K=11,a=5,\alpha =0.9,\beta =0.3,m=0.1,\eta =0.01$, $d_1=0.002$, $d_2=2$, and $\Omega =(0,\pi )$. Figure 8 shows the change of critical value ${\tau }_0$ when the cost of fear b varies in $[2,50]$, which indicates that the stability of $E^{\ast }$ can be changed by the cost of fear and moreover, the Hopf bifurcation can also be induced by it.

5.4. The Effect of Anti-Predator Parameter $\eta$

Choose the parameters as $r=0.9,b=10,K=11,a=5,\alpha =0.9,\beta =0.3,m=0.1$, $d_1=0.002$, $d_2=2$, and $\Omega =(0,\pi )$, and let the parameter vary in $[0,0.01]$. Figure 9 shows that with the increase of the parameter $\eta$, the stability region increases first and then decreases, which means that the stability of $E^{\ast }$ can be changed by the parameter $\eta$ and moreover, the Hopf bifurcation can also be induced by it. In addition, Figure 10 shows that with the increase of $\eta$, the biomass of the prey increases, and the biomass of the predator increases first and then decreases.

6. Conclusions

In this paper, we took into account the fact that prey population growth rates naturally slow down when predator apprehension rises. However, the fear effect does not immediately slow the expansion of a prey species; rather, there is a lag time before the effect is felt. We have taken into account the possibility of a lag time between the moment a prey animal perceives the fear of a predator through chemical/vocal cues and the time that fear slows its rate of growth. To make our model more accurate, we included the delay phenomena. The primary objective of this research is to investigate the effect of postponing the cost of fear on the dynamics of a population. To get there, we looked at a diffusive predator–prey model with and without delay analytically and numerically, and found that the latter was more advanced.

First, we obtained the existence of the solution of system (1), and the stability of the equilibria, including the extinct equilibrium, the predator-free equilibrium, and the positive equilibrium.

Then we investigated the Hopf bifurcation induced by the time delay, and gave some theoretical results about the non-existence and existence of the non-constant solution of the proposed model.

We finally completed our investigation into the delayed model and uncovered the complex dynamics of the system through numerical analysis. While the non-delay system displays stable coexistence, an internal equilibrium point’s stability can be changed from steady focus to limit cycle oscillation and back again when time is added or subtracted. The system enters the chaotic phase after a large number of transitions between stable focus and limit cycle oscillation. For this reason, in a straightforward two-dimensional predator–prey model with fear effect, a larger value of delay can lead to anarchy. In addition, we further explored the effect of the cost of fear; the stability of the inter equilibrium can be changed by the cost of fear, moreover, it can also cause Hopf bifurcation.

Recently, many experts have studied real-world problems such as applications of the extended rational sine–cosine and sinh–cosh techniques to some nonlinear complex models arising in mathematical physics [29]. However, application of this model in this paper will require real data and newly developed analytical methods to manage a prey–predator system.