Dynamical analysis of a delayed diffusive predator–prey model with schooling behaviour and Allee effect

Abstract

In this paper, a delayed diffusive predator–prey model with schooling behaviour and Allee effect is investigated. The existence and local stability of equilibria of model without time delay and diffusion are given. Regarding the conversion rate as bifurcation parameter, Hopf bifurcation of diffusive system without time delay is obtained. In addition, the local stability of the coexistent equilibrium and existence of Hopf bifurcation of system with time delay are discussed. Moreover, the properties of Hopf bifurcation are studied based on the centre manifold and normal form theory for partial functional differential equations. Finally, some numerical simulations are also carried out to confirm our theoretical results.

Keywords:

• Schooling behaviour
• time delay
• Hopf bifurcation
• Predator–prey
• Allee effect

1. Introduction

In natural ecosystems, many species live together in groups. Although selfish behaviour can sometimes have an evolutionary advantage, all individuals face the same public goods as pointed out in reference [31]. There are benefits from living in groups, which, for example, increase the chance of finding mates and foraging effectively for food, reduce the risk of capture and weather and so on.

This type of herd among different species takes different forms due to the species evolution, such as forms of line, square, oval and amoeboid. Recently, the predator–prey model in which one population gets together in group while the other exhibits individualistic has been investigated by many authors [6,9,11,17,26,29,33,34,40,44–46]. Braza [6] analysed the dynamics of the square root system with Holling Type-II response functions. Zhu et al. [46] considered two-dimensional predator–prey model with the square root and delay. Moreno and Plataniab [26] presented a good method of flexibility and provided a high analytical tractability corresponding to the interest rate level by using harmonic oscillators. Tang et al. [32] chose diffusion coefficients and delay as the bifurcation parameters, and investigated the dynamics of predator–prey model with group behaviour and hyperbolic mortality.

In fact, many fish species live in a form of herd. The way of herd life is divided mainly into two types: shoaling type and schooling type, where shoaling behaviour is that members of a group swim alone in the manner of connecting and schooling behaviour is that all members swim in the same direction and same speed. Both behaviours depend on selecting mate or other facts [27]. Many biologists have found that large numbers of fishes exhibit schooling behaviour. For example, Brattey [5] found adult code in dense schools. Mitsunaga et al. [25] studied schooling behaviour of juvenile tuna. Of course, some researchers [7,16,24] proposed that small fishes such as herring, capelin and sardine show spectacular schooling behaviour. Thus it is important to consider the schooling behaviour in the predator–prey system for above fishes. In the predator–prey system of fishes, the consideration of schooling behaviour can better explain their mechanism. Ajraldi et al. [1] proposed a new predator–prey model which the individuals of one population gather together in herds, while the other one shows a more individualistic behaviour, (1) $\begin{array}{l}\frac{du}{dt}=ru\left(1-\frac{u}{K}\right)-a\sqrt{u}v,\\ \frac{dv}{dt}=a\gamma \sqrt{u}v-\text{β}v,\end{array}$(1) where $u\left(t\right)$ and $v\left(t\right)$ stand for the prey density and predator density at the time t, respectively. r denotes the prey net reproduction rate and K denotes their carrying capacity. a is the predators' hunting rate, γ is the conversion or consumption rate of prey to predator, and β is the death rate of the predator in the absence of prey.

What's more, it is well known that most countries of the world achieve the economic benefits of natural resource management by restoring and maintaining the ecosystem's health, the productivity and biodiversity and the overall quality of life in a manner that integrates social and economic objectives. Of course, it also satisfies the need for humans to benefit from natural resources. The prey–predator systems with harvesting have been widely studied by a large number of scholars [13,18–23]. Meng et al. [19] studied a predator–prey system with harvesting prey and disease in prey species and found that the optimal harvesting effort is closely related to the incubation period of the infectious disease, and the maximum value of the optimal harvesting decreases with the increase of the time delay. Although some works have already been done to understand the predator–prey system considering fishes with schooling behaviour by authors [8,30], the study of predator–prey system in which both predator and prey take into account schooling behaviour is rare so far. Manna et al. [15] considered fishes containing both predator and prey with schooling behaviour and proposed the model as follows: (2) $\begin{array}{l}\frac{dx}{dt}=rx\left(1-\frac{x}{k}\right)-b\sqrt{x}\sqrt{y}-q_{1}Ex,\\ \frac{dy}{dt}=c\sqrt{x}\sqrt{y}-\mathrm{d}y-q_{2}Ey,\end{array}$(2) where $x\left(t\right)$ and $y\left(t\right)$ denote the density of prey fishes and predator fishes at time t, respectively. It is assumed that the biomass conversion of predator is expressed by square root functional response which both predator and prey with schooling behaviour.

Recently, the Allee effect has been introduced in the growth of population widely [3,4,20,28,35,36,42]. Allee [2] proposed the Allee effect, which mainly divided into two classes: strong and weak [35]. Strong Allee is a phenomenon that population will extinct when the density below the fixed threshold, and the weak Allee is the case that the growth rate of population reduces still keeps positive at low population density [4]. Lewis and Kareiva [12] considered the strong Allee effect in the prey growth rate $M(1-M)(M-\text{β}),$where β quantifies the intensity of strong Allee effect with $0<\text{β}<1$.

In recent years, diffusion has attracted a large attention in the biological models [11,14,37,41]. Since fishes distribute in ocean inhomogeneous in different space at one time t. Thus diffusion should be taken into the model. In a real ecological environment, the time delay incorporating into the resource limitation of the prey logistic equation also should be considered. Taking the effect of the time delay in resource limitation of the prey logistic equation into predator–prey model can make the model be realistic. Generally speaking, there are some economic benefits for some fish in aquatic food chain. Some fish is harvested for human food. Thus harvesting is an essential element in our life. On the one hand, the harvesting is one of the economic sources of human economy. On the other hand, it is one method to control the balance of ecosystem with harvesting. The main aim of this paper intends to see whether the time delay can exert an impact on the delayed-diffusion model with schooling behaviour, Allee effect and linear predator harvesting.

Combining the work of references [1,15] with Allee effect, the delayed-diffusive model is proposed as follows: (3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) where $M(x,t)$ and $N(x,t)$ denote the density of prey fishes and predator fishes at time t and position x, respectively. $d_1$ and $d_2$ are the diffusion constants for the prey fishes and predator fishes. α and γ are the prey fishes uptake rate and conversion rate of which required for growth of the predator fishes, respectively. δ is the mortality rate of the predator fishes. τ is the time delay incorporated into the resource limitation of the prey logistic equation. The constant q is the catchability coefficient of the number of the predator fishes, and E is the harvesting effort.

The organization of this paper is as follows. In Section 2, the existence and local stability of the equilibria of system (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) without delay and diffusion are discussed. In Section 3, the local stability of the coexistent equilibrium and the existence of Hopf bifurcation of system (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) without time delay are investigated. In Section 4, Hopf bifurcation and its properties of system (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) with time delay based on normal form theory and centre manifold theorem for partial functional differential equation are also discussed. In Section 5, numerical simulations are given to illustrate the results. Some discussions and conclusions are included in the end.

2. System without diffusion and delay

System (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) without diffusion and time delay becomes the form of an ordinary differential system (4) $\begin{array}{l}\dot{M}=M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \dot{N}=\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\end{array}$(4) with initial conditions $M\left(0\right)=\varphi \left(0\right)\ge 0,N\left(0\right)=\psi \left(0\right)\ge 0.$

2.1. The existence of all equilibria

In this section, we will discuss the existence of equilibria of system (4(4) $\begin{array}{l}\dot{M}=M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \dot{N}=\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\end{array}$(4) ).

It is easy to obtain that system (4(4) $\begin{array}{l}\dot{M}=M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \dot{N}=\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\end{array}$(4) ) has four equilibriums, the extinction equilibrium $S_0(0,0)$, the predator extinction equilibria $S_1(1,0)$ and $S_{\text{β}}(\text{β},0)$, and the coexistence equilibrium $S^{\ast }(M^{\ast },N^{\ast })$.

For the coexistence equilibrium $S^{\ast }(M^{\ast },N^{\ast })$, the $M^{\ast }$ and $N^{\ast }$ satisfy the following equations: (5) $\begin{array}{l}\sqrt{M}(1-M)(M-\text{β})-\alpha \sqrt{N}=0,\\ \gamma \sqrt{M}-\delta \sqrt{N}-qE\sqrt{N}=0.\end{array}$(5) From the second equation of (5(5) $\begin{array}{l}\sqrt{M}(1-M)(M-\text{β})-\alpha \sqrt{N}=0,\\ \gamma \sqrt{M}-\delta \sqrt{N}-qE\sqrt{N}=0.\end{array}$(5) ), we get (6) $N={\left(\frac{\gamma }{\delta +qE}\right)}^{2}M,$(6) and put (6(6) $N={\left(\frac{\gamma }{\delta +qE}\right)}^{2}M,$(6) ) into the first equation of (5(5) $\begin{array}{l}\sqrt{M}(1-M)(M-\text{β})-\alpha \sqrt{N}=0,\\ \gamma \sqrt{M}-\delta \sqrt{N}-qE\sqrt{N}=0.\end{array}$(5) ), which leads to the equation (7) $M^2-(1+\text{β})M+\text{β}+\frac{\alpha \gamma }{\delta +qE}=0.$(7) If $\mathrm{\Delta }=(1-\text{β}{)}^2-\frac{4\alpha \gamma }{\delta +qE}=0$, that is $\gamma =\frac{1}{4\alpha }(1-\text{β}{)}^2(\delta +qE)$, then Equation (7(7) $M^2-(1+\text{β})M+\text{β}+\frac{\alpha \gamma }{\delta +qE}=0.$(7) ) has only one positive real root $M_1^{\ast }=\frac{1+\text{β}}{2}.$If $\mathrm{\Delta }>0$, that is $\gamma <\frac{1}{4\alpha }(1-\text{β}{)}^2(\delta +qE)$, then Equation (7(7) $M^2-(1+\text{β})M+\text{β}+\frac{\alpha \gamma }{\delta +qE}=0.$(7) ) has two different positive real roots $M_2^{\ast }=\frac{1+\text{β}+\sqrt{\mathrm{\Delta }}}{2},M_3^{\ast }=\frac{1+\text{β}-\sqrt{\mathrm{\Delta }}}{2}.$From (6(6) $N={\left(\frac{\gamma }{\delta +qE}\right)}^{2}M,$(6) ), we have $N_1^{\ast }=\frac{(1+\text{β}){\gamma }^2}{2(\delta +qE{)}^2},N_2^{\ast }=\frac{(1+\text{β}+\sqrt{\mathrm{\Delta }}){\gamma }^2}{2(\delta +qE{)}^2},N_3^{\ast }=\frac{(1+\text{β}-\sqrt{\mathrm{\Delta }}){\gamma }^2}{2(\delta +qE{)}^2}.$From above analysis, we obtain the following results.

Theorem 2.1

1. If $0<\text{β}<1$, then the extinction equilibrium $S_0(0,0)$, the predator extinction equilibrium $S_1(1,0)$ and $S_{\text{β}}(\text{β},0)$ always exist.

2. If $0<\text{β}<1$, then the coexistence equilibrium $S_1^{\ast }(M_1^{\ast },N_1^{\ast })$ exists when $\gamma =\frac{1}{4\alpha }(1-\text{β}{)}^2(\delta +qE)$; the coexistence equilibria $S_2^{\ast }(M_2^{\ast },N_2^{\ast })$ and $S_3^{\ast }(M_3^{\ast },N_3^{\ast })$ exist when $\gamma <\frac{1}{4\alpha }(1-\text{β}{)}^2(\delta +qE)$.

In order to better analyse the equilibria of system (4(4) $\begin{array}{l}\dot{M}=M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \dot{N}=\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\end{array}$(4) ), we show the existence of non-negative equilibria in Figure 1. The values of parameter except β are $\alpha =0.9,q=0.3,E=15,\gamma =0.8,\delta =0.1$. Those figures are the possible number of equilibria under five different cases. The manganese violet lines represent the curves of implicit function $M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N}=0$, and the green lines represent the curves of implicit function $\gamma \sqrt{M}\sqrt{N}-\delta N-qEN=0$. (a) When $\text{β}=0.22$, the two curves do not intersect at any point in the interior of the feasible domain (the green curves is below the manganese violet lines), which means that there is no positive equilibrium; (b) when $\text{β}=0.21$, there are two equal equilibria, that is there is only one equilibrium; (c) when $\text{β}=0.15$, the both curves cross twice, which means that there are two equilibriums; (d) when $\text{β}=0.1$ and E = 480, there are two boundary equilibriums; (e) when $\text{β}=0.01$, there are two positive equilibriums and one extinction equilibrium.

Figure 1. The possible number of equilibria of system (4(4) $\begin{array}{l}\dot{M}=M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \dot{N}=\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\end{array}$(4) ) with different values of parameter β: (a) $\text{β}=0.22$, (b) $\text{β}=0.221$, (c) $\text{β}=0.15$, (d) $\text{β}=0.1$, (e) $\text{β}=0.01$.

2.2. The local stability of equilibria

In this section, we discuss the local stability of system (4(4) $\begin{array}{l}\dot{M}=M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \dot{N}=\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\end{array}$(4) ) at equilibria.

Because system (4(4) $\begin{array}{l}\dot{M}=M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \dot{N}=\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\end{array}$(4) ) is not linearizable about $S_0(0,0)$, $S_1(1,0)$ and $S_{\text{β}}(\text{β},0)$, the local stability of them cannot be studied. In the following, we discuss the local stability of coexistence equilibria of system (4(4) $\begin{array}{l}\dot{M}=M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \dot{N}=\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\end{array}$(4) ) when $0<\text{β}<1$.

Let $\overline{S}=(\overline{M},\overline{N})$ be the arbitrary coexistence equilibrium, we have a linear system of system (4(4) $\begin{array}{l}\dot{M}=M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \dot{N}=\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\end{array}$(4) ) at $\overline{S}$ as follows: (8) $\dot{u}\left(t\right)=\overline{A}u\left(t\right),$(8) where $u\left(t\right)=\left(M\right(t),N(t){)}^T$, $\overline{A}=\left(\begin{array}{cc}-3{\overline{M}}^2+2(\text{β}+1)\overline{M}-\text{β}-\frac{\alpha \sqrt{\overline{N}}}{2\sqrt{\overline{M}}}& -\frac{\alpha \sqrt{\overline{M}}}{2\sqrt{\overline{N}}}\\ \frac{\gamma \sqrt{\overline{N}}}{2\sqrt{\overline{M}}}& \frac{\gamma \sqrt{\overline{M}}}{2\sqrt{\overline{N}}}-\delta -qE\end{array}\right)\doteq \left(\begin{array}{cc}{\overline{A}}_{11}& {\overline{A}}_{12}\\ {\overline{A}}_{21}& {\overline{A}}_{22}\end{array}\right).$The characteristic equation of system (8(8) $\dot{u}\left(t\right)=\overline{A}u\left(t\right),$(8) ) is (9) ${\lambda }^2-({\overline{A}}_{11}+{\overline{A}}_{22})\lambda +{\overline{A}}_{11}{\overline{A}}_{22}-{\overline{A}}_{12}{\overline{A}}_{21}=0.$(9) According to the characteristic equation (9(9) ${\lambda }^2-({\overline{A}}_{11}+{\overline{A}}_{22})\lambda +{\overline{A}}_{11}{\overline{A}}_{22}-{\overline{A}}_{12}{\overline{A}}_{21}=0.$(9) ) and the coexistence equilibrium $\overline{S}=S_1^{\ast }(M_1^{\ast },N_1^{\ast })$, we have that $\text{β}=1-2\sqrt{\frac{\alpha \gamma }{\delta +qE}}$, and $\begin{array}{l}{\lambda }_1+{\lambda }_2={\overline{A}}_{11}+{\overline{A}}_{22}=\frac{\alpha \gamma }{2(\delta +qE)}-\frac{\delta +qE}{2},\\ {\lambda }_1{\lambda }_2={\overline{A}}_{11}{\overline{A}}_{22}-{\overline{A}}_{12}{\overline{A}}_{21}=\left[\frac{1}{4}(\text{β}-1{)}^2-\frac{\alpha \gamma }{2(\delta +qE)}\right]\left(-\frac{\delta +qE}{2}\right)+\frac{\alpha \gamma }{4}.\end{array}$We know that $\gamma >\frac{1}{\alpha }(\delta +qE)$ from $\text{β}=1-2\sqrt{\frac{\alpha \gamma }{\delta +qE}}$ and $0<\text{β}<1$. Therefore, we have ${\lambda }_1=0$, ${\lambda }_2=\frac{\alpha \gamma }{2(\delta +qE)}-\frac{\delta +qE}{2}$ for $\gamma >\frac{1}{\alpha }(\delta +qE)$, that is the positive equilibrium $S_1^{\ast }(M_1^{\ast },N_1^{\ast })$ is unstable.

Similarly, we have that $S_2^{\ast }$ and $S_3^{\ast }$ exist if $\gamma <\frac{1}{4\alpha }(1-\text{β}{)}^2(\delta +qE)$ from Theorem 2.1. We have that at the positive equilibrium $\overline{S}=S_2^{\ast }(M_2^{\ast },N_2^{\ast })$, $\begin{array}{l}{\lambda }_1+{\lambda }_2=-\frac{1}{2}(\text{β}-1{)}^2-\frac{1}{2}(\text{β}+1)\sqrt{\mathrm{\Delta }}+\frac{5\alpha \gamma }{2(\delta +qE)}-\frac{\delta +qE}{2},\\ {\lambda }_1{\lambda }_2=\left[-\frac{1}{2}(\text{β}-1{)}^2-\frac{1}{2}(\text{β}+1)\sqrt{\mathrm{\Delta }}+\frac{5\alpha \gamma }{2(\delta +qE)}\right]\left(-\frac{\delta +qE}{2}\right)+\frac{\alpha \gamma }{4}.\end{array}$If $\gamma <\frac{1}{5\alpha }(1-\text{β}{)}^2(\delta +qE)$, then we have $-\frac{1}{2}(\text{β}-1{)}^2-\frac{1}{2}(\text{β}+1)\sqrt{\mathrm{\Delta }}+\frac{5\alpha \gamma }{2(\delta +qE)}<0$ which leads to that ${\lambda }_1+{\lambda }_2<0$ and ${\lambda }_1{\lambda }_2>0$. Therefore, when $\gamma <\frac{1}{5\alpha }(1-\text{β}{)}^2(\delta +qE)$, Equation (9(9) ${\lambda }^2-({\overline{A}}_{11}+{\overline{A}}_{22})\lambda +{\overline{A}}_{11}{\overline{A}}_{22}-{\overline{A}}_{12}{\overline{A}}_{21}=0.$(9) ) has two negative real part roots, that is the positive equilibrium $S_2^{\ast }(M_2^{\ast },N_2^{\ast })$ is locally asymptotically stable.

Similarly, when $\overline{S}=S_3^{\ast }(M_3^{\ast },N_3^{\ast })$, Equation (9(9) ${\lambda }^2-({\overline{A}}_{11}+{\overline{A}}_{22})\lambda +{\overline{A}}_{11}{\overline{A}}_{22}-{\overline{A}}_{12}{\overline{A}}_{21}=0.$(9) ) has two negative real part roots if $\frac{1}{\alpha }(\delta +qE{)}^2<\gamma <\frac{1}{4\alpha }(1-\text{β}{)}^2(\delta +qE)$, that is the positive equilibrium $S_3^{\ast }(M_3^{\ast },N_3^{\ast })$ is locally asymptotically stable.

In order to analyse the local stability of the extinction equilibrium, we redefine the variables $M\left(t\right)$, $N\left(t\right)$ as $M\left(t\right)=X^2\left(t\right)$, $N\left(t\right)=Y^2\left(t\right)$. Then system (4(4) $\begin{array}{l}\dot{M}=M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \dot{N}=\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\end{array}$(4) ) transforms to the following form: (10) $\begin{array}{l}\dot{X}=-\frac{1}{2}[X^5-(1+\text{β})X^3+\text{β}X+\alpha Y],\\ \dot{Y}=\frac{1}{2}[\gamma X-(\delta +qE\left)Y\right].\end{array}$(10) Clearly, system (10(10) $\begin{array}{l}\dot{X}=-\frac{1}{2}[X^5-(1+\text{β})X^3+\text{β}X+\alpha Y],\\ \dot{Y}=\frac{1}{2}[\gamma X-(\delta +qE\left)Y\right].\end{array}$(10) ) has two equilibrium points, namely the extinction equilibrium ${\tilde{S}}_0(0,0)$and the coexistence equilibrium ${\tilde{S}}_{\ast }(X^{\ast },Y^{\ast })$. At the coexistence  equilibrium, we have $Y^{\ast }=\frac{\gamma }{\delta +qE}{X}^{\ast },$and X satisfies (11) $(X^{\ast }{)}^5-(1+\text{β}\left)\right(X^{\ast }{)}^3+\left(\text{β}+\frac{\gamma \alpha }{\delta +qE}\right)X^{\ast }=0.$(11) Discussing the local stability of system (4(4) $\begin{array}{l}\dot{M}=M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \dot{N}=\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\end{array}$(4) ) is equivalent to discussing the local stability of system (10(10) $\begin{array}{l}\dot{X}=-\frac{1}{2}[X^5-(1+\text{β})X^3+\text{β}X+\alpha Y],\\ \dot{Y}=\frac{1}{2}[\gamma X-(\delta +qE\left)Y\right].\end{array}$(10) ) at the extinction equilibrium ${\tilde{S}}_0(0,0)$ when $0<\text{β}<1$. Let $\tilde{S}=(\tilde{X},\tilde{Y})$ be the arbitrary equilibrium of system (10(10) $\begin{array}{l}\dot{X}=-\frac{1}{2}[X^5-(1+\text{β})X^3+\text{β}X+\alpha Y],\\ \dot{Y}=\frac{1}{2}[\gamma X-(\delta +qE\left)Y\right].\end{array}$(10) ), we have a linear system of system (10(10) $\begin{array}{l}\dot{X}=-\frac{1}{2}[X^5-(1+\text{β})X^3+\text{β}X+\alpha Y],\\ \dot{Y}=\frac{1}{2}[\gamma X-(\delta +qE\left)Y\right].\end{array}$(10) ) at $\tilde{S}$ as follows: (12) $\dot{v}\left(t\right)=\tilde{A}v\left(t\right),$(12) where $v\left(t\right)=\left(X\right(t),Y(t){)}^T$, $\tilde{A}=\left(\begin{array}{lr}-\frac{5}{2}{X}^4+\frac{3}{2}(1+\text{β})X^2-\frac{1}{2}\text{β}& -\frac{\alpha }{2}\\ \frac{\gamma }{2}& -\frac{1}{2}(\delta +qE)\end{array}\right)=\left(\begin{array}{cc}{\tilde{A}}_{11}& {\tilde{A}}_{12}\\ {\tilde{A}}_{21}& {\tilde{A}}_{22}\end{array}\right).$The characteristic equation of system (12(12) $\dot{v}\left(t\right)=\tilde{A}v\left(t\right),$(12) ) is (13) ${\lambda }^2-({\tilde{A}}_{11}+{\tilde{A}}_{22})\lambda +{\tilde{A}}_{11}{\tilde{A}}_{22}-{\tilde{A}}_{12}{\tilde{A}}_{21}=0.$(13) The characteristic equation (13(13) ${\lambda }^2-({\tilde{A}}_{11}+{\tilde{A}}_{22})\lambda +{\tilde{A}}_{11}{\tilde{A}}_{22}-{\tilde{A}}_{12}{\tilde{A}}_{21}=0.$(13) ) at ${\tilde{S}}_0(0,0)$ is (14) ${\lambda }^2+\frac{1}{2}(\text{β}+\delta +qE)\lambda +\frac{1}{4}\left[\text{β}\right(\delta +qE)+\alpha \gamma ]=0.$(14) From Equation (14(14) ${\lambda }^2+\frac{1}{2}(\text{β}+\delta +qE)\lambda +\frac{1}{4}\left[\text{β}\right(\delta +qE)+\alpha \gamma ]=0.$(14) ), we can clearly see that ${\lambda }_1+{\lambda }_2=-\frac{1}{2}(\text{β}+\delta +qE)<0$ and ${\lambda }_1{\lambda }_2=\frac{1}{4}\left[\text{β}\right(\delta +qE)+\alpha \gamma ]>0$. Therefore, Equation (14(14) ${\lambda }^2+\frac{1}{2}(\text{β}+\delta +qE)\lambda +\frac{1}{4}\left[\text{β}\right(\delta +qE)+\alpha \gamma ]=0.$(14) ) has two negative real part roots. Thus the extinction equilibrium ${\tilde{S}}_0(0,0)$ is always locally asymptotically stable.

The characteristic equation (13(13) ${\lambda }^2-({\tilde{A}}_{11}+{\tilde{A}}_{22})\lambda +{\tilde{A}}_{11}{\tilde{A}}_{22}-{\tilde{A}}_{12}{\tilde{A}}_{21}=0.$(13) ) at ${\tilde{S}}_{\ast }(X^{\ast },Y^{\ast })$ is (15) $\begin{array}{rl}& {\lambda }^2+\left[(1+\text{β})X_{\ast }^2-2\text{β}-\frac{5\alpha \gamma }{2(\delta +qE)}+\frac{1}{2}(\delta +qE)\right]\lambda \\ & +\frac{1}{2}(1+\text{β})(\delta +qE)X_{\ast }^2-(\delta +qE)-\alpha \gamma =0.\end{array}$(15) If $\gamma <\frac{2}{5\alpha }(\delta +qE)\left(\right(1+\text{β})X_{\ast }^2-2\text{β})$, then Equation (15(15) $\begin{array}{rl}& {\lambda }^2+\left[(1+\text{β})X_{\ast }^2-2\text{β}-\frac{5\alpha \gamma }{2(\delta +qE)}+\frac{1}{2}(\delta +qE)\right]\lambda \\ & +\frac{1}{2}(1+\text{β})(\delta +qE)X_{\ast }^2-(\delta +qE)-\alpha \gamma =0.\end{array}$(15) ) has two negative real part roots, that is the positive equilibrium ${\tilde{S}}_{\ast }(X^{\ast },Y^{\ast })$ is locally asymptotically stable.

By analysis above, we get the following results.

Theorem 2.2

1. If $0<\text{β}<1$, then the extinction equilibrium ${\tilde{S}}_0(0,0)$ is always locally asymptotically stable for all γ.

2. If $\gamma <\frac{2}{5\alpha }(\delta +qE)\left(\right(1+\text{β})X_{\ast }^2-2\text{β})$, then the positive equilibrium ${\tilde{S}}_{\ast }(X^{\ast },Y^{\ast })$ is locally asymptotically stable.

In order to better show the existence and the stability of the equilibria between system (4(4) $\begin{array}{l}\dot{M}=M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \dot{N}=\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\end{array}$(4) ) and (10(10) $\begin{array}{l}\dot{X}=-\frac{1}{2}[X^5-(1+\text{β})X^3+\text{β}X+\alpha Y],\\ \dot{Y}=\frac{1}{2}[\gamma X-(\delta +qE\left)Y\right].\end{array}$(10) ), we summarize those results shown in Tables 1–4 as follows.

Table 1. Feasibility condition for equilibria of system (4(4) $\begin{array}{l}\dot{M}=M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \dot{N}=\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\end{array}$(4) ).

Equilibrium	Interpretation	Feasibility condition
$S_0(0,0)$	system disappearance	always
$S_1(1,0)$	prey only	always
$S_{\text{β}}(\text{β},0)$	prey only	always
$S_1^{\ast }(M_1^{\ast },N_1^{\ast })$	whole system coexistence	$\gamma =\frac{1}{4\alpha }(1-\text{β}{)}^2(\delta +qE)$
$S_2^{\ast }(M_2^{\ast },N_2^{\ast })$	whole system coexistence	$\gamma <\frac{1}{4\alpha }(1-\text{β}{)}^2(\delta +qE)$
$S_3^{\ast }(M_3^{\ast },N_3^{\ast })$	whole system coexistence	$\gamma <\frac{1}{4\alpha }(1-\text{β}{)}^2(\delta +qE)$

Table 2. Stability condition for equilibria of system (4(4) $\begin{array}{l}\dot{M}=M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \dot{N}=\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\end{array}$(4) ).

Equilibrium	Interpretation	Stability condition
$S_0(0,0)$	system disappearance	–
$S_1(1,0)$	prey only	–
$S_{\text{β}}(\text{β},0)$	prey only	–
$S_1^{\ast }(M_1^{\ast },N_1^{\ast })$	whole system coexistence	$\gamma >\frac{1}{\alpha }(\delta +qE)$
$S_2^{\ast }(M_2^{\ast },N_2^{\ast })$	whole system coexistence	$\gamma <\frac{1}{5\alpha }(1-\text{β}{)}^2(\delta +qE)$
$S_3^{\ast }(M_3^{\ast },N_3^{\ast })$	whole system coexistence	$\frac{1}{\alpha }(\delta +qE{)}^2<\gamma <\frac{1}{4\alpha }(1-\text{β}{)}^2(\delta +qE)$

Table 3. Feasibility condition for equilibria of system (10(10) $\begin{array}{l}\dot{X}=-\frac{1}{2}[X^5-(1+\text{β})X^3+\text{β}X+\alpha Y],\\ \dot{Y}=\frac{1}{2}[\gamma X-(\delta +qE\left)Y\right].\end{array}$(10) ).

Equilibrium	Interpretation	Feasibility condition
${\tilde{S}}_0(0,0)$	system disappearance	–
${\tilde{S}}_{\ast }(X^{\ast },Y^{\ast })$	whole system coexistence	satisfies (11(11) $(X^{\ast }{)}^5-(1+\text{β}\left)\right(X^{\ast }{)}^3+\left(\text{β}+\frac{\gamma \alpha }{\delta +qE}\right)X^{\ast }=0.$(11) )

Table 4. Stability condition for equilibria of system (10(10) $\begin{array}{l}\dot{X}=-\frac{1}{2}[X^5-(1+\text{β})X^3+\text{β}X+\alpha Y],\\ \dot{Y}=\frac{1}{2}[\gamma X-(\delta +qE\left)Y\right].\end{array}$(10) ).

Equilibrium	Interpretation	Stability condition
${\tilde{S}}_0(0,0)$	system disappearance	stable
${\tilde{S}}_{\ast }(X^{\ast },Y^{\ast })$	whole system coexistence	$\gamma <\frac{2}{5\alpha }(\delta +qE)\left(\right(1+\text{β})X_{\ast }^2-2\text{β})$

3. Bifurcation analysis of system without delay

In this section, we will analyse system (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) without delay and study the existence of Hopf bifurcation in one-dimensional space domain $\mathrm{\Omega }=(0,\pi )$. We will only discuss the dynamics of the coexistence equilibrium $S_2^{\ast }$. For convenience, we denote the coexistence equilibrium $S_2^{\ast }$= $S^{\ast }$.

When $\tau =0$, system (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) becomes following diffusion equations: (16) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t>0,\\ M(x,0)=\varphi \left(x\right)\ge 0,N(x,0)=\psi \left(x\right)\ge 0,x\in [0,\pi ].\end{array}$(16) We denote that $k^2(k\in {\mathbb{N}}_0)$ is  eigenvalue of $-{\mathrm{\partial }}^2/\mathrm{\partial }{x}^2$ with the Neumann boundary conditions in the $\mathrm{\Omega }=(0,\pi )$. In the following discussion, we take the conversion rate γ as bifurcation parameter. Then, the linear operator of system (16(16) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t>0,\\ M(x,0)=\varphi \left(x\right)\ge 0,N(x,0)=\psi \left(x\right)\ge 0,x\in [0,\pi ].\end{array}$(16) ) at $S^{\ast }$ is $L\left(\gamma \right)=\left(\begin{array}{cc}{d}_1\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}+L_{11}& L_{12}\\ L_{21}& d_2\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}+L_{22}\end{array}\right),$where $\begin{array}{l}{L}_{11}=-\frac{1}{2}(\text{β}-1{)}^2-\frac{1}{2}(\text{β}+1)\sqrt{\mathrm{\Delta }}+\frac{5\alpha \gamma }{2(\delta +qE)},L_{12}=-\frac{\alpha (\delta +qE)}{2\gamma },\\ L_{21}=\frac{{\gamma }^2}{2(\delta +qE)},L_{22}=-\frac{\delta +qE}{2}.\end{array}$Based on [10,43], the eigenvalue of $L\left(\gamma \right)$ is given by the sequence operator $L_k\left(\gamma \right)$, where $\begin{array}{l}{L}_k\left(\gamma \right)=\left(\begin{array}{cc}-d_1{k}^2+L_{11}& L_{12}\\ L_{21}& -d_2{k}^2+L_{22}\end{array}\right).\end{array}$The corresponding characteristic equation is (17) ${\lambda }^2-T_k\lambda +J_k=0,$(17) where (18) $\begin{array}{l}{T}_k=-k^2(d_1+d_2)+L_{11}+L_{22},\\ J_k=d_1{d}_2{k}^4-(d_1{L}_{22}+d_2{L}_{11})k^2+L_{11}{L}_{22}-L_{12}{L}_{21}.\end{array}$(18) From Theorem 2.1, we know that $T_k<0$ and $J_k>0$ when $\gamma <{\gamma }_1$, here ${\gamma }_1=\frac{1}{5\alpha }(1-\text{β}{)}^2(\delta +qE)$. Therefore, all eigenvalues of Equation (17(17) ${\lambda }^2-T_k\lambda +J_k=0,$(17) ) have negative real part roots. So the coexistence equilibrium $S^{\ast }$ is locally asymptotic stability.

Now, we consider the situation in the condition ${\gamma }_1<\gamma <{\gamma }_2$, here ${\gamma }_2=\frac{1}{4\alpha }(1-\text{β}{)}^2(\delta +qE)$, and investigate the existence of Hopf bifurcations of system (16(16) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t>0,\\ M(x,0)=\varphi \left(x\right)\ge 0,N(x,0)=\psi \left(x\right)\ge 0,x\in [0,\pi ].\end{array}$(16) ).

Based on [43], if there exists $j\in {\mathbb{N}}_0$ for certain critical value ${\gamma }_k^H$ such that $T_j\left({\gamma }_k^H\right)=0,J_j\left({\gamma }_k^H\right)>0,and{T}_i\left({\gamma }_k^H\right)\ne 0,J_i\left({\gamma }_k^H\right)\ne 0,fori\ne j,$and for a only pair of complex eigenvalue $\alpha \left(\gamma \right)\pm i\omega \left(\gamma \right)$ near imaginary axis, which leads to (19) ${\alpha }^{\prime }\left({\gamma }_k^H\right)\ne 0,$(19) then system (16(16) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t>0,\\ M(x,0)=\varphi \left(x\right)\ge 0,N(x,0)=\psi \left(x\right)\ge 0,x\in [0,\pi ].\end{array}$(16) ) undergoes Hopf bifurcation at the $S^{\ast }$ near the $\gamma ={\gamma }_k^H$. From Equation (18(18) $\begin{array}{l}{T}_k=-k^2(d_1+d_2)+L_{11}+L_{22},\\ J_k=d_1{d}_2{k}^4-(d_1{L}_{22}+d_2{L}_{11})k^2+L_{11}{L}_{22}-L_{12}{L}_{21}.\end{array}$(18) ), we get $T_j\left({\gamma }_k^H\right)=0$ is equivalent to ${\gamma }_k^H=\frac{1}{2\alpha }[-2k^2(d_1+d_2)+(\text{β}-1{)}^2-(\delta +qE\left)]\right(\delta +qE).$There are finite bifurcating points ${\gamma }_k^H\in ({\gamma }_1,{\gamma }_2)$. Namely, there exists a non-negative inter $N_k\in {\mathbb{N}}_0$ such that ${\gamma }_k^H$ is bifurcating points for $0\le k\le N_k$, otherwise, ${\gamma }_k^H$ is not bifurcating points for $k>N_k$, which is the largest integer. Thus we have that $k^2\le \frac{L_{11}+L_{12}}{d_1+d_1}$.

Lemma 3.1

There are $N_k+1$ bifurcating points ${\gamma }_k^H\in ({\gamma }_1,{\gamma }_2)$ satisfying $0<{\gamma }_{N_k}^H<\cdots {\gamma }_k^H\cdots <{\gamma }_1^H<{\gamma }_0^H,$for any $0\le k\le N_k$. It is clearly that $T_j\left({\gamma }_k^H\right)=0$ and $T_i\left({\gamma }_k^H\right)\ne 0$ for any $j\ne i$.

In order to illustrate that $J_j\left({\gamma }_k^H\right)>0$ is true for all $j\in {\mathbb{N}}_0$, we give one lemma as follows.

Lemma 3.2

Suppose that $0<\text{β}<1$, if the hypothesis $\left(H\right):\frac{L_{12}{L}_{21}-2L_{11}{L}_{22}-\sqrt{{\mathrm{\Delta }}_k}}{L_{22}^2}<\frac{d_1}{d_2}<\frac{L_{12}{L}_{21}-L_{11}{L}_{22}+\sqrt{{\mathrm{\Delta }}_k}}{L_{22}^2},$holds, where ${\mathrm{\Delta }}_k=(4L_{12}{L}_{21}-2L_{11}{L}_{22}{)}^2-L_{11}^2{L}_{22}^2>0$ for $\gamma <{\gamma }_2$. Then $J_j\left({\gamma }_k^H\right)>0$ holds for all $j\in {\mathbb{N}}_0$.

Now, we consider the transverse condition in (19(19) ${\alpha }^{\prime }\left({\gamma }_k^H\right)\ne 0,$(19) ).

Lemma 3.3

Suppose that the eigenvalues near the purely imaginary are $\alpha \left(\gamma \right)\pm i\omega \left(\gamma \right)$ at $\gamma ={\gamma }_k^H$ for $0\le k\le N_k$. Then, we have ${\alpha }^{\prime }\left({\gamma }_k^H\right)=\frac{T_k^{\prime }\left({\gamma }_k^H\right)}{2}=-\frac{\alpha }{2(\sigma +qE)}<0,\omega \left({\gamma }_k^H\right)=\sqrt{J_j\left({\gamma }_k^H\right)}>0.$

Combining with the above discussion, we obtain the following conclusion.

Theorem 3.4

The periodic solutions are spatially homogeneous bifurcating from $\gamma ={\gamma }_0^H$, but the periodic solutions are spatially non-homogeneous bifurcating from $\gamma ={\gamma }_{N_k}^H$ for $1\le k\le N_k$.

Remark 3.5

The method of discussing the positive equilibria $S_1^{\ast }(M_1^{\ast },N_1^{\ast })$ and $S_3^{\ast }(M_3^{\ast },N_3^{\ast })$ is similar to that of the positive equilibrium $S_2^{\ast }(M_2^{\ast },N_2^{\ast })$. Thus we omit the corresponding results.

4. Hopf bifurcation and its properties of system with delay

In the following, we will study the effect of delay on system. By linearizing system (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) at interior $S^{\ast }(M^{\ast },N^{\ast })$, we obtain (20) $\left(\begin{array}{c}\frac{\mathrm{\partial }\left(M\right)}{\mathrm{\partial }\left(t\right)}\\ \frac{\mathrm{\partial }\left(N\right)}{\mathrm{\partial }\left(t\right)}\end{array}\right)=D\mathrm{\Delta }\left(\begin{array}{c}M\left(t\right)\\ N\left(t\right)\end{array}\right)+A\left(\begin{array}{c}M\left(t\right)\\ N\left(t\right)\end{array}\right)+B\left(\begin{array}{c}M(t-\tau )\\ N(t-\tau )\end{array}\right),$(20) with (21) $D=\left(\begin{array}{cc}{d}_1& 0\\ 0& d_2\end{array}\right),A=\left(\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right),B=\left(\begin{array}{cc}{B}_{11}& 0\\ 0& 0\end{array}\right),$(21) where $\begin{array}{rl}{A}_{11}& =\frac{1-\text{β}}{2}+\frac{\alpha \gamma }{2(\delta +qE)},A_{12}=-\frac{\alpha (\delta +qE)}{2\gamma },A_{21}=\frac{{\gamma }^2}{2(\delta +qE)},\\ A_{22}& =-\frac{\delta +qE}{2},B_{11}=-\frac{{\text{β}}^2-1}{2}.\end{array}$Obviously, the characteristic equation of the linearization of (20(20) $\left(\begin{array}{c}\frac{\mathrm{\partial }\left(M\right)}{\mathrm{\partial }\left(t\right)}\\ \frac{\mathrm{\partial }\left(N\right)}{\mathrm{\partial }\left(t\right)}\end{array}\right)=D\mathrm{\Delta }\left(\begin{array}{c}M\left(t\right)\\ N\left(t\right)\end{array}\right)+A\left(\begin{array}{c}M\left(t\right)\\ N\left(t\right)\end{array}\right)+B\left(\begin{array}{c}M(t-\tau )\\ N(t-\tau )\end{array}\right),$(20) ) is equivalent to the sequence of the transcendental equation (22) $det(\lambda I-D_k-A-B{e}^{-\lambda \tau })=0,$(22) where I is the $2\times 2$ identify matrix and $D_k=-k^{2}diag\{d_1,d_1\}$, $k\in {\mathbb{N}}_0=0,1,2,\dots$

According to Equation (22(22) $det(\lambda I-D_k-A-B{e}^{-\lambda \tau })=0,$(22) ), the characteristic equations for the positive equilibrium $S^{\ast }$ are the following sequence of quadratic transcendental equation: (23) ${\lambda }^2+T_1+T_2+T_3-(\lambda +T_2)B_{11}{e}^{-\lambda \tau }=0,$(23) where $T_1=d_1{k}^2-A_{11},T_2=d_2{k}^2-A_{22},T_3=T_1{T}_2-A_{12}{A}_{21}$.

4.1. Hopf bifurcation induced by delay

In this part, we will discuss the existence of Hopf bifurcation by considering time delay as bifurcation parameter. From Theorem 2.2, the positive equilibrium $S^{\ast }$ of system (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) is locally asymptotically stable when $\gamma <{\gamma }_1$. Now, we assume that $\lambda =i\omega (\omega >0)$ is a root of the characteristic equation (23(23) ${\lambda }^2+T_1+T_2+T_3-(\lambda +T_2)B_{11}{e}^{-\lambda \tau }=0,$(23) ) when $\gamma <{\gamma }_1$. Since the following equation about ω is (24) $-{\omega }^2+(T_1+T_2-B_{11}{e}^{-i\omega \tau })\lambda +T_3-T_2{B}_{11}{e}^{-i\omega \tau }=0.$(24) Separating real and imaginary parts, we get (25) $\begin{array}{l}{T}_2{B}_{11}\cos \left(\omega \tau \right)+B_{11}\omega \sin \left(\omega \tau \right)=-{\omega }^2+T_3,\\ B_{11}\omega \cos \left(\omega \tau \right)-T_2{B}_{11}\sin \left(\omega \tau \right)=\omega (T_1+T_2).\end{array}$(25) Then, by squaring two equations above and adding up them together, we obtain (26) ${\omega }^4+B_k{\omega }^2+C_k=0,$(26) where $B_k=(T_1+T_2{)}^2-2T_3-B_{11}^2$, $C_k=T_3^2-(T_2{B}_{11}{)}^2$.

In order to discuss the existence of the positive roots of Equation (23(23) ${\lambda }^2+T_1+T_2+T_3-(\lambda +T_2)B_{11}{e}^{-\lambda \tau }=0,$(23) ), we need discuss the distribution of the roots of Equation (26(26) ${\omega }^4+B_k{\omega }^2+C_k=0,$(26) ). Let ${\mathrm{\Delta }}_{\omega }=B_k^2-4C_k$. Some assumptions are given as follows:

1. (i) ${\mathrm{\Delta }}_{\omega }<0$; (ii) ${\mathrm{\Delta }}_{\omega }>0,B_k>0$ and $C_k\ge 0$; (iii) ${\mathrm{\Delta }}_{\omega }=0,B_k\ge 0$,

2. (i) $C_k\ge 0$; (ii) ${\mathrm{\Delta }}_{\omega }=0,B_k<0$,

3. ${\mathrm{\Delta }}_{\omega }>0$, $B_k<0$ and $C_k\ge 0$.

Lemma 4.1

1. If any one of the $\left(H_1\right)$ holds, then Equation (26(26) ${\omega }^4+B_k{\omega }^2+C_k=0,$(26) ) has no real positive root, thus Equation (23(23) ${\lambda }^2+T_1+T_2+T_3-(\lambda +T_2)B_{11}{e}^{-\lambda \tau }=0,$(23) ) has no purely imaginary roots with $\tau \ne 0$.

2. If any one of the $\left(H_2\right)$ holds, then Equation (26(26) ${\omega }^4+B_k{\omega }^2+C_k=0,$(26) ) has one real positive root, that is, Equation (23(23) ${\lambda }^2+T_1+T_2+T_3-(\lambda +T_2)B_{11}{e}^{-\lambda \tau }=0,$(23) ) has a pair of purely imaginary roots $\pm i{\omega }_k^{+}=\pm i\sqrt{\frac{-B_k+\sqrt{{\mathrm{\Delta }}_{\omega }}}{2}}$.

3. If $\left(H_3\right)$ holds, then Equation (26(26) ${\omega }^4+B_k{\omega }^2+C_k=0,$(26) ) has two real positive roots, thus Equation (23(23) ${\lambda }^2+T_1+T_2+T_3-(\lambda +T_2)B_{11}{e}^{-\lambda \tau }=0,$(23) ) has two pair of purely imaginary roots $\pm i{\omega }_k^{\pm }=\pm i\sqrt{\frac{-B_k\pm \sqrt{{\mathrm{\Delta }}_{\omega }}}{2}}$.

From (25(25) $\begin{array}{l}{T}_2{B}_{11}\cos \left(\omega \tau \right)+B_{11}\omega \sin \left(\omega \tau \right)=-{\omega }^2+T_3,\\ B_{11}\omega \cos \left(\omega \tau \right)-T_2{B}_{11}\sin \left(\omega \tau \right)=\omega (T_1+T_2).\end{array}$(25) ), we have (27) ${\tau }_{k,j}^{\pm }=\frac{1}{{\omega }_k^{\pm }}\arccos \frac{T_1({\omega }_k^{\pm }{)}^2+T_2{T}_3}{(T_2^2+{\omega }_k^2)B_{11}}+\frac{2j\pi }{{\omega }_k^{\pm }},j=0,1,2,\dots$(27) Let ${\lambda }_k\left(\tau \right)={\alpha }_k\left({\tau }_{k,j}\right)+i{\omega }_k^{\pm }\left(\tau \right)$ be the root of Equation (23(23) ${\lambda }^2+T_1+T_2+T_3-(\lambda +T_2)B_{11}{e}^{-\lambda \tau }=0,$(23) ) such that ${\alpha }_k\left({\tau }_{k,j}\right)=0$ and ${\omega }_k\left({\tau }_{k,j}\right)={\omega }_k$. Then we can get the following transversality conditions.

Lemma 4.2

Suppose $\gamma <{\gamma }_1$. If $\left(H_2\right)$ and $\left(H_3\right)$ hold, then we have

1. when ${\mathrm{\Delta }}_{\omega }=0$, then $\mathrm{R}\mathrm{e}\left(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\right){\mid }_{\tau ={\tau }_{k,j}^{\pm }}=0;$

2. when ${\mathrm{\Delta }}_{\omega }>0$, then $\mathrm{R}\mathrm{e}\left(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\right){\mid }_{\tau ={\tau }_{k,j}^{+}}>0$, $\mathrm{R}\mathrm{e}\left(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\right){\mid }_{\tau ={\tau }_{k,j}^{-}}<0$.

Proof.

Differentiating both sides of Equation (23(23) ${\lambda }^2+T_1+T_2+T_3-(\lambda +T_2)B_{11}{e}^{-\lambda \tau }=0,$(23) ) with corresponding to τ, we get ${\left(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\right)}^{-1}=-\frac{(2\lambda +T_1+T_2)e^{\lambda \tau }}{\lambda (\lambda +T_2)B_{11}}+\frac{1}{\lambda (\lambda +T_2)}-\frac{\tau }{\lambda }.$From (27(27) ${\tau }_{k,j}^{\pm }=\frac{1}{{\omega }_k^{\pm }}\arccos \frac{T_1({\omega }_k^{\pm }{)}^2+T_2{T}_3}{(T_2^2+{\omega }_k^2)B_{11}}+\frac{2j\pi }{{\omega }_k^{\pm }},j=0,1,2,\dots$(27) ), we have $\begin{array}{rl}& \mathrm{R}\mathrm{e}{\left[{\left(\frac{\mathrm{d}\lambda }{d\tau }\right)}^{-1}\right]}_{\lambda =\pm i{\omega }_k^{\pm }}\\ & =\mathrm{R}\mathrm{e}{\left[\left(-\frac{\left.2i\omega +T_1+T_2\right)e^{-i\omega \tau }}{i\omega B_{11}(i\omega +T_2)}-\frac{1}{i\omega (i\omega +T_2)}\right.\right]}_{\lambda =\pm i{\omega }_k^{\pm }}\\ & =\mathrm{R}\mathrm{e}{\left[-\frac{\left.2i\omega +T_1+T_2\right)(\cos \omega \tau +\sin \omega \tau )}{i\omega B_{11}(i\omega +T_2)}-\frac{1}{i\omega (i\omega +T_2)}\right]}_{\lambda =\pm i{\omega }_k^{\pm }}\\ & ={\left[\frac{2{\omega }^2+(T_1+T_2{)}^2-2T_3-B_{11}^2}{B_{11}^2({\omega }^2+T_2^2)}\right]}_{\lambda =\pm i{\omega }_k^{\pm }}\\ & ={\left[\frac{\pm \sqrt{{\mathrm{\Delta }}_{\omega }}}{B_{11}^2({\omega }^2+T_2^2)}\right]}_{\lambda =\pm i{\omega }_k^{\pm }}.\end{array}.$Thus when ${\mathrm{\Delta }}_{\omega }\ne 0$, we have $\mathrm{R}\mathrm{e}{\left[{\left(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\right)}^{-1}\right]}_{\lambda =i{\omega }_k^{+}}>0,\mathrm{R}\mathrm{e}{\left[{\left(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\right)}^{-1}\right]}_{\lambda =i{\omega }_k^{-}}<0.$The proof of this lemma is completed.

From (27(27) ${\tau }_{k,j}^{\pm }=\frac{1}{{\omega }_k^{\pm }}\arccos \frac{T_1({\omega }_k^{\pm }{)}^2+T_2{T}_3}{(T_2^2+{\omega }_k^2)B_{11}}+\frac{2j\pi }{{\omega }_k^{\pm }},j=0,1,2,\dots$(27) ), we get ${\tau }_{k,0}^{\pm }<{\tau }_{k,j}^{\pm }.$ Define ${\tau }_{k,0}=\min \left\{{\tau }_{k,0}^{\pm }\mathrm{o}\mathrm{r}{\tau }_{k,0}^{+}\right\},0\le k\le N_k.$

Theorem 4.3

Suppose that $\gamma <{\gamma }_1$. The following results are true.

1. If the assumption $\left(H_1\right)$ holds, then $S^{\ast }$ is locally asymptotically stable whenever $\tau \ge 0$.

2. If one of the assumptions $\left(H_2\right)$ or $\left(H_3\right)$ holds, then $S^{\ast }$ is locally asymptotically stable whenever $\tau \in [0,{\tau }_{k,o})$, and it is unstable when $\tau >{\tau }_{k,0}$.

3. If one of the assumption $\left(H_2\right)$ or $\left(H_3\right)$ holds, The system (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) undergoes a Hopf bifurcation around $S^{\ast }$ for $\tau ={\tau }_{k,j}^{\pm },0\le k\le N_k,j=0,1,2,\dots$

4.2. Properties of Hopf bifurcation

In this section, we will consider the direction of the Hopf bifurcation, the stability and the period of bifurcating periodic solution of system (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) by applying the centre manifold and normal form theory [10,38]. In reminder of this paper, we only give the main conclusions.

Let $\tilde{M}(x,t)=M(x,\tau t)-M^{\ast },\tilde{N}(x,t)=N(x,\tau t)-N^{\ast }$, we drop the tilde for simplification of nation. Then we transform the system (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) into the following equations: (28) $\begin{array}{rl}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}& =\tau [d_1\mathrm{\Delta }M+(M+M^{\ast })(1-M(t-1)-M^{\ast })(M+M^{\ast }-\text{β})\\ & -\alpha \sqrt{M+M^{\ast }}\sqrt{N+N^{\ast }}],\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}& =\tau [d_2\mathrm{\Delta }N+\gamma \sqrt{M+M^{\ast }}\sqrt{N+N^{\ast }}-\delta (N+N^{\ast })-qE(N+N^{\ast })].\end{array}$(28) for $x\in (0,\pi )$, and t>0. Let $\tau =\tilde{\tau }+\sigma ,U_1\left(t\right)=M(\cdot ,t),U_2\left(t\right)=N(\cdot ,t),U\left(t\right)=(U_1,U_2{)}^T.$Then (28(28) $\begin{array}{rl}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}& =\tau [d_1\mathrm{\Delta }M+(M+M^{\ast })(1-M(t-1)-M^{\ast })(M+M^{\ast }-\text{β})\\ & -\alpha \sqrt{M+M^{\ast }}\sqrt{N+N^{\ast }}],\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}& =\tau [d_2\mathrm{\Delta }N+\gamma \sqrt{M+M^{\ast }}\sqrt{N+N^{\ast }}-\delta (N+N^{\ast })-qE(N+N^{\ast })].\end{array}$(28) ) is written as an abstract equation in the space $\mathcal{C}=C\left(\right[-1,0],X)$, here $X=\left\{\right(M,N)\in W^{2,2}(0,\pi )\mid \frac{\mathrm{\partial }M}{\mathrm{\partial }x}=\frac{\mathrm{\partial }N}{\mathrm{\partial }x}=0atx=0,\pi \}$, (29) $\frac{dU\left(t\right)}{dt}=\tau D\mathrm{\Delta }U\left(t\right)+L\left(\tau \right)U_t+F(U_t,\sigma ),$(29) where $L\left(\tau \right)(\cdot ):\mathcal{C}\to X$ and $F:\mathcal{C}\times R\to X$ are given by $\begin{array}{l}L\left(\tau \right)\left(\varphi \right)=\tau A\left(\begin{array}{c}{\varphi }_1\left(0\right)\\ {\varphi }_2\left(0\right)\end{array}\right)+\tau B\left(\begin{array}{c}{\varphi }_1(-1)\\ {\varphi }_2(-1)\end{array}\right),F(U_t,\sigma )=\left(\begin{array}{c}{f}_1(U_t,\sigma )\\ f_2(U_t,\sigma )\end{array}\right),\end{array}$here A, B, D defined as (21(21) $D=\left(\begin{array}{cc}{d}_1& 0\\ 0& d_2\end{array}\right),A=\left(\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right),B=\left(\begin{array}{cc}{B}_{11}& 0\\ 0& 0\end{array}\right),$(21) ), $\varphi =({\varphi }_1,{\varphi }_2{)}^T\in \mathcal{C}$ and $\begin{array}{rl}{f}_1(\varphi ,\sigma )& =\frac{1}{2}\left[f_{11}{\varphi }_1^2\left(0\right)+2f_{12}{\varphi }_1\left(0\right){\varphi }_1(-1)+2f_{13}{\varphi }_1\left(0\right){\varphi }_2\left(0\right)+f_{33}{\varphi }_2^2\left(0\right)\right]\\ & +\frac{1}{6}\left[f_{111}{\varphi }_1^3\left(0\right)+3f_{113}{\varphi }_1^2\left(0\right){\varphi }_2\left(0\right)+3f_{133}{\varphi }_1\left(0\right){\varphi }_2^2\left(0\right)+f_{333}{\varphi }_2^3\left(0\right)\right],\\ f_2(\varphi ,\sigma )& =\frac{1}{2}\left[g_{11}{\varphi }_1^2\left(0\right)+2g_{13}{\varphi }_1\left(0\right){\varphi }_2\left(0\right)+g_{33}{\varphi }_2^2\left(0\right)\right]\\ & +\frac{1}{6}\left[g_{111}{\varphi }_1^3\left(0\right)+3g_{113}{\varphi }_1^2\left(0\right){\varphi }_2\left(0\right)+3g_{133}{\varphi }_1\left(0\right){\varphi }_2^2\left(0\right)+g_{333}{\varphi }_2^3\left(0\right)\right],\end{array}$here $\begin{array}{rl}{f}_{11}& =\frac{1-{\text{β}}^2}{2}+\frac{\alpha \gamma }{2(\delta +qE)(1+\text{β})},f_{113}=\frac{\alpha }{2(1+\text{β}{)}^2}{\left(\frac{\gamma }{(\delta +qE)}\right)}^{-1},\\ g_{13}& =\frac{\delta +qE}{(1+\text{β}{)}^2},f_{12}=\frac{\text{β}-1}{2},g_{111}=-\frac{3{\gamma }^2}{2(\delta +qE)(1+\text{β}{)}^2},\\ g_{333}& =-\frac{3\gamma }{2(1+\text{β}{)}^2}{\left(\frac{\gamma }{\delta +qE}\right)}^{-5},f_{13}=-\frac{\alpha (\delta +qE)}{\gamma (1+\text{β}{)}^2},\\ f_{111}& =-\frac{3\alpha \gamma }{2(\delta +qE)(1+\text{β}{)}^2},f_{333}=-\frac{3\alpha }{2(1+\text{β}{)}^2}{\left(\frac{\gamma }{\delta +qE}\right)}^{-5},\\ f_{33}& =-\frac{\alpha }{2(1+\text{β})}{\left(\frac{\gamma }{\delta +qE}\right)}^{-3},f_{133}=\frac{\alpha }{2(1+\text{β}{)}^2}{\left(\frac{\gamma }{\delta +qE}\right)}^{-3},\\ g_{11}& =-\frac{{\gamma }^2}{2(\delta +qE)(1+\text{β})},g_{33}=\frac{\gamma }{2(1+\text{β})}{\left(\frac{\gamma }{\delta +qE}\right)}^{-3},\\ g_{113}& =-\frac{\delta +qE}{2(1+\text{β}{)}^2},g_{133}=-\frac{\gamma }{2(1+\text{β}{)}^2}{\left(\frac{\gamma }{\delta +qE}\right)}^{-3}.\end{array}$It is clear that the origin $(0,0)$ is an equilibrium of (28(28) $\begin{array}{rl}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}& =\tau [d_1\mathrm{\Delta }M+(M+M^{\ast })(1-M(t-1)-M^{\ast })(M+M^{\ast }-\text{β})\\ & -\alpha \sqrt{M+M^{\ast }}\sqrt{N+N^{\ast }}],\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}& =\tau [d_2\mathrm{\Delta }N+\gamma \sqrt{M+M^{\ast }}\sqrt{N+N^{\ast }}-\delta (N+N^{\ast })-qE(N+N^{\ast })].\end{array}$(28) ). When $\sigma =0$, the linear equation is (30) $\frac{dU\left(t\right)}{dt}=\tilde{\tau }D\mathrm{\Delta }U\left(t\right)+L\left(\tilde{\tau }\right)U_t,$(30) which has a pair of simple purely imaginary eigenvalues ${\mathrm{\Lambda }}_k=\{i{\omega }_k\tilde{\tau },-i{\omega }_k\tilde{\tau }\}$.

Next, we consider the linear functional differential equation (31) $\frac{dZ\left(t\right)}{dt}=-\tilde{\tau }D{K}^2\mathrm{\Delta }Z\left(t\right)+L\left(\tilde{\tau }\right)Z_t.$(31) Applying the Reisz representation theorem [39], there exists $2\times 2$ matrix-valued function ${\zeta }_k(\theta ,\sigma )$ of bounded variation for $\theta \in [-1,0]$ such that $-\tilde{\tau }D{K}^2\varphi \left(0\right)+L_{\tilde{\tau }}\left(\varphi \right)={\int }_{-1}^{0}d{\zeta }_k(\theta ,\sigma )\varphi \left(\theta \right),for\varphi \in C\left(\right[-1,0],R^2).$In fact, we can take ${\zeta }_k(\theta ,\sigma )=\tilde{\tau }\left[A\delta \right(\theta )+B\delta (\theta +1\left)\right],$where $\delta \left(\theta \right)$ is a Dirac delta function.

Further, we suppose that $A_{\tilde{\tau }}$ is the infinitesimal generator corresponding to (31(31) $\frac{dZ\left(t\right)}{dt}=-\tilde{\tau }D{K}^2\mathrm{\Delta }Z\left(t\right)+L\left(\tilde{\tau }\right)Z_t.$(31) ) and $A^{\ast }$ is the formal adjoint of A under the bilinear paring (32) $\begin{array}{rl}〈\psi \left(s\right),\varphi \left(\theta \right)〉& =\overline{\psi }\left(0\right)\varphi \left(0\right)-{\int }_{-1}^0{\int }_0^{\theta }\overline{\psi }(\xi -\theta )\mathrm{d}{\zeta }_k\left(\theta \right)\varphi \left(\xi \right)d\xi \\ & =\overline{\psi }\left(0\right)\varphi \left(0\right)+\tilde{\tau }{\int }_{-1}^0\overline{\psi }(\xi +1)B\varphi \left(\xi \right)d\xi ,\end{array}$(32) for $\varphi \in C^1\left(\right[-1,0],R^2)$, $\psi \in C^1\left(\right[-1,0],R^2)$, where ${\zeta }_K\left(\theta \right)=\zeta (\theta ,0)$. From the above discussion, we can know that $\pm i{\omega }_k$ are the eigenvalues of $A_{\tilde{\tau }}$, thus $\pm i{\omega }_k$ are also the eigenvalues of $A^{\ast }$.

Let $p\left(\theta \right)=(1,p_1{)}^T{e}^{i{\omega }_k\tilde{\tau }\theta }$ and $p^{\ast }\left(\theta \right)=B_1(1,p_2{)}^T{e}^{-i{\omega }_k\tilde{\tau }\theta }$ be eigenvectors of $A_k$ and $A^{\ast }$ respecting to the eigenvalues $i{\omega }_k\tilde{\tau }$ and $-i{\omega }_k\tilde{\tau }$, respectively. We obtain $p_1=\frac{A_{21}}{i{\omega }_k-A_{22}},p_2=-\frac{A_{12}}{i{\omega }_k+A_{22}}.$We can choose the value of $B_1$ as $\overline{B_1}=[1+p_1\overline{p_2}+\tilde{\tau }{B}_{11}{e}^{-i{\omega }_k\tilde{\tau }}{]}^{-1},$such that $〈p^{\ast },p〉=1$. Thus we can computer the following quantities: $\begin{array}{rl}{g}_{20}& =\left\{\begin{array}{ll}0,& k=1,2,\dots ,N_k,\\ \tilde{\tau }{B}_1[f_{11}+2f_{13}{p}_1+f_{33}{p}_1^2+2f_{12}{p}_1{e}^{-i{\omega }_k\tilde{\tau }}\\ +\overline{p_2}(g_{11}+2g_{13}{p}_1+g_{33}{p}_1^2)],& k=0,\end{array}\right.\\ g_{11}& =\left\{\begin{array}{ll}0,& k=1,2,\dots ,N_k,\\ \tilde{\tau }{B}_1[f_{11}+f_{13}(p_1+\overline{p_1})+f_{33}{p}_1\overline{p_1}\\ +f_{12}{p}_1{e}^{-i{\omega }_k\tilde{\tau }}+f_{12}{p}_1{e}^{i{\omega }_k\tilde{\tau }}& k=0,\\ +\overline{p_2}(g_{11}+g_{13}({\overline{p}}_1+p_1)+g_{33}{p}_1{\overline{p}}_1)],\end{array}\right.\\ g_{02}& =\left\{\begin{array}{ll}0,& k=1,2,\dots ,N_k,\\ \tilde{\tau }{B}_1[f_{11}+2f_{13}{\overline{p}}_1+f_{33}{\overline{p}}_1^2+2f_{12}{\overline{p}}_1{e}^{-i{\omega }_k\tilde{\tau }}\\ +\overline{p_2}(g_{11}+2g_{13}{\overline{p}}_1+g_{33}{\overline{p}}_1^2)],& k=0,\end{array}\right.\\ g_{21}& =\frac{\tilde{\tau }\overline{B}}{\pi }\left[(f_{11}+g_{11}+(f_{13}+g_{13})p_1+f_{12}{\overline{p}}_1{e}^{-i{\omega }_k\tilde{\tau }})2W_{11}^{\left(0\right)}\left(0\right)\right.\\ & +(f_{13}+g_{13}+(f_{33}+g_{33}\left)p_1\right)2W_{11}^{\left(2\right)}\left(0\right)+\left(f_{11}+g_{11}+(f_{13}+g_{13}){\overline{p}}_1\right.\\ & \left.+f_{12}{\overline{p}}_1{e}^{-i{\omega }_k\tilde{\tau }}\right)W_{20}^{\left(1\right)}\left(0\right)+(f_{13}+g_{13}+(f_{33}+g_{33}\left){\overline{p}}_1\right)W_{20}^{\left(2\right)}\left(0\right)\\ & \left.+f_{12}{W}_{20}^{\left(1\right)}(-1)+f_{12}{W}_{11}^{\left(1\right)}(-1)\right]{\int }_0^{\pi }{\cos }^2\left(kx\right)dx\\ & +\frac{\tilde{\tau }\overline{B}}{\pi }\left[\left(\frac{1}{4}(f_{111}+g_{111})+\frac{1}{12}(f_{113}+g_{113})({\overline{p}}_1+2p_1)\right.\right.\\ & \left.+\frac{1}{12}(f_{133}+g_{133})(2{\overline{p}}_1+p_1)+\frac{1}{4}(f_{333}+g_{333}){\overline{p}}_1{p}_1^2\right]{\int }_0^{\pi }{\cos }^4\left(kx\right)dx.\end{array}$where $k=0,1,2,\dots ,N_k$, and $\begin{array}{rl}{W}_{20}\left(\theta \right)& =\left(\frac{i{g}_{20}p\left(\theta \right)}{{\omega }_k\tilde{\tau }}+\frac{i{\overline{g}}_{20}\overline{p}\left(\theta \right)}{3{\omega }_k\tilde{\tau }}\right)\cdot f_k+E_1{e}^{2i{\omega }_k\tilde{\tau }\theta },\\ W_{11}\left(\theta \right)& =\left(\frac{-i{g}_{11}p\left(\theta \right)}{{\omega }_k\tilde{\tau }}+\frac{i{\overline{g}}_{11}\overline{p}\left(\theta \right)}{3{\omega }_k\tilde{\tau }}\right)\cdot f_k+E_2,\end{array}$the definition of $f_k$ is the same to that in the reference [38], and $\begin{array}{rl}{E}_1& =\tilde{\tau }E{\chi }_1{\cos }^2\left(kx\right),E_2=\tilde{\tau }{E}^{\ast }{\chi }_2{\cos }^2\left(kx\right),\\ E_1& ={\left(\begin{array}{cc}i{\omega }_k+d_1{K}^2-A_{11}-B_{11}{e}^{-2i{\omega }_k\tilde{\tau }}& -A_{12}\\ -A_{21}& i{\omega }_k+d_2{K}^2-A_{22}\end{array}\right)}^{-1},\\ E^{\ast }& ={\left(\begin{array}{cc}{d}_1{K}^2-A_{11}-B_{11}& -A_{12}\\ -A_{21}& d_2{K}^2-A_{22}\end{array}\right)}^{-1},\\ {\chi }_1& =\left(\begin{array}{c}{f}_{11}+2f_{13}{p}_1+f_{33}{p}_1^2+2f_{12}{p}_1{e}^{-i{\omega }_k\tilde{\tau }}\\ g_{11}+2g_{13}{p}_1+g_{33}{p}_1^2\end{array}\right),\\ {\chi }_2& =\left(\begin{array}{c}{f}_{11}+f_{13}(p_1+\overline{p_1})+f_{33}{p}_1\overline{p_1}+f_{12}{p}_1{e}^{-i{\omega }_k\tilde{\tau }}+f_{12}{p}_1{e}^{i{\omega }_k\tilde{\tau }}\\ g_{11}+g_{13}({\overline{p}}_1+p_1)+g_{33}{p}_1{\overline{p}}_1\end{array}\right).\end{array}$Based on above analysis, the following values can be given as (33) $\begin{array}{rl}{\mathcal{C}}_1\left(0\right)& =\frac{i}{2{\omega }_k\tilde{\tau }}\left(g_{20}{g}_{11}-2|g_{11}{|}^2-\frac{|g_{02}{|}^2}{3}\right)+\frac{g_{21}}{2},\\ {\mu }_2& =-\frac{Re\left\{{\mathcal{C}}_1\left(0\right)\right\}}{Re\left\{{\lambda }^{\prime }\left({\tau }_0\right)\right\}},\\ {\text{β}}_2& =2Re\left\{{\mathcal{C}}_1\left(0\right)\right\},\\ T_2& =-\frac{Im\left\{{\mathcal{C}}_1\left(0\right)\right\}+{\mu }_{2}Im\left\{{\lambda }^{\prime }\left({\tau }_0\right)\right\}}{{\omega }_k\tilde{\tau }}.\end{array}$(33) By computing (33(33) $\begin{array}{rl}{\mathcal{C}}_1\left(0\right)& =\frac{i}{2{\omega }_k\tilde{\tau }}\left(g_{20}{g}_{11}-2|g_{11}{|}^2-\frac{|g_{02}{|}^2}{3}\right)+\frac{g_{21}}{2},\\ {\mu }_2& =-\frac{Re\left\{{\mathcal{C}}_1\left(0\right)\right\}}{Re\left\{{\lambda }^{\prime }\left({\tau }_0\right)\right\}},\\ {\text{β}}_2& =2Re\left\{{\mathcal{C}}_1\left(0\right)\right\},\\ T_2& =-\frac{Im\left\{{\mathcal{C}}_1\left(0\right)\right\}+{\mu }_{2}Im\left\{{\lambda }^{\prime }\left({\tau }_0\right)\right\}}{{\omega }_k\tilde{\tau }}.\end{array}$(33) ), we get the following theorem.

Theorem 4.4

Properties of bifurcating periodic solutions are determined as follows.

1. If ${\mu }_2>0({\mu }_2<0)$, then Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for $\tau >\tilde{\tau }(\tau <\tilde{\tau })$.

2. If ${\text{β}}_2<0({\text{β}}_2>0)$, then bifurcating periodic solutions on the centre manifold remain stable (unstable).

3. If $T_2>0(T_2<0)$, then the period of the periodic solution increases (decreases).

5. Numerical simulation

In this section, we carry out numerical simulations to verify our theoretical results by the use of software MATLAB. Let $\mathrm{\Omega }=[0,2\pi ](i.e.l=2)$, and most values of parameters in system (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) are mentioned in reference [1], $\text{β}=0.01,\alpha =0.9,\gamma =0.8,\delta =0.1,q=0.3,E=15,d_1=d_2=\mathrm{0.5.}$We obtain that the coexistent equilibrium of system with $\tau =0$ is $S^{\ast }=(0.802,0.243)$, $T_k<0$ and $J_k>0$ when $\gamma <{\gamma }_1$. Therefore, the coexistent equilibrium $S^{\ast }$ is locally asymptotically stable (see Figure  2).

Figure 2. System (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) at the coexistent equilibrium $S^{\ast }$ is locally asymptotically stable: (a) $M(x,t)$; (b) $N(x,t)$.

From Figure 2, we know that $S^{\ast }=(M^{\ast },N^{\ast })$ is locally asymptotically stable when $\tau =0$. As the time increases, the prey fishes and predator fishes population fluctuate at first and then tend to be stable.

By only calculation, we obtain ${\tau }_{k,0}={\tau }_{0,0}=2.03,forall0\le k\le N_k.$Thus the critical value is ${\tau }_{\ast }={\tau }_{0,0}=2.03$. By Theorem 3.4, system (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) is the locally asymptotically stable when $0<\tau <{\tau }_{\ast }$. The stability of system (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) is shown in Figure 3,

Figure 3. System (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) at the coexistent equilibrium $S^{\ast }$ is locally asymptotically stable when $\tau =40<{\tau }_{\ast }$: (a) $M(x,t)$; (b) $N(x,t)$.

From Figure 3, we know that $S^{\ast }=(M^{\ast },N^{\ast })$ is still asymptotically stable when $0<\tau <{\tau }_{\ast }$. As the time increases, the prey fishes and predator fishes population first fluctuate and even tend to be stable.

Further, according to (33(33) $\begin{array}{rl}{\mathcal{C}}_1\left(0\right)& =\frac{i}{2{\omega }_k\tilde{\tau }}\left(g_{20}{g}_{11}-2|g_{11}{|}^2-\frac{|g_{02}{|}^2}{3}\right)+\frac{g_{21}}{2},\\ {\mu }_2& =-\frac{Re\left\{{\mathcal{C}}_1\left(0\right)\right\}}{Re\left\{{\lambda }^{\prime }\left({\tau }_0\right)\right\}},\\ {\text{β}}_2& =2Re\left\{{\mathcal{C}}_1\left(0\right)\right\},\\ T_2& =-\frac{Im\left\{{\mathcal{C}}_1\left(0\right)\right\}+{\mu }_{2}Im\left\{{\lambda }^{\prime }\left({\tau }_0\right)\right\}}{{\omega }_k\tilde{\tau }}.\end{array}$(33) ), we obtain the values as follows: ${\mathcal{C}}_1\left(0\right)\approx -0.0491-0.1306i,{\mu }_2\approx 0.4931,{\text{β}}_2\approx -0.0982,T_2\approx \mathrm{0.2133.}$Therefore, the periodic solutions is stable (see Figure 4), Hopf-bifurcating is supercritical, and the period of the periodic solution increases.

Figure 4. The periodic solution is stable when $\tau =2.3>{\tau }_{\ast }$. (a) $M(x,t)$; (b) $N(x,t)$.

From Figure 4, we chose $\tau =2.3$. That is, when time delay τ increases across its critical value ${\tau }_{\ast }\approx 2.03$, a spatially periodic solution occurs near the positive equilibrium $S^{\ast }=(M^{\ast },N^{\ast })$. However, the bifurcating periodic solution bifurcating from the critical value ${\tau }_{\ast }$ may be stable on the whole phase space. A regular oscillation of the prey fishes and predator fishes population are showed.

In addition, although the local stability of the predator extinction equilibria is not given in the previous part, we can show their stability under some proper values of all parameters by numerical simulations. For example, we consider system (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) with the following values: $\text{β}=0.1,\alpha =0.9,\gamma =0.2,\delta =0.1,q=0.23,E=480,d_1=d_2=1.$The predator extinction equilibrium $S_1(1,0)$ is always locally asymptotically stable, which is shown in Figure 5.

Figure 5. The predator extinction equilibrium $S_1(1,0)$ is locally asymptotically stable: (a) $M(x,t)$; (b) $N(x,t)$.

From Figure 5, we can see that the number of prey fishes reaches its maximum, but the number of predator tends to zero.

Similarly, we consider system (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) with the following values: $\text{β}=0.9,\alpha =0.3,\gamma =0.2,\delta =0.1,q=0.01,E=10,d_1=d_2=1.$From simple computation, we obtain that $\gamma >\frac{1}{\alpha }(\delta +qE{)}^2\approx 0.13$. The predator extinction equilibrium $S_{\text{β}}(0.9,0)$ is locally asymptotically stable, which is shown in Figure 6.

Figure 6. The predator extinction equilibrium $S_{\text{β}}(0.9,0)$ is locally asymptotically stable: (a) $M(x,t)$; (b) $N(x,t)$.

From Figure 6, we can see that the number of prey fishes reaches its maximum, but the number of predator goes to zero.

6. Conclusion and discussion

Ajraldi et al. [1] explained the phenomenon that the individuals of one population gather together in herds, while the other shows a more individualistic behaviour. Manna et al. [15] proposed the predator–prey model considering both predator and prey with schooling behaviour. In this paper, the dynamic of a delayed diffusive predator–prey model with schooling behaviour, Allee effect and harvesting was analysed. Considering the advantage on the edge of capture of marginal predator, therefore, the relationship of them was expressed by the square root function. Meanwhile, Allee effect had been introduced in the prey intrinsic growth rate by considering the mate limitation and completion. Furthermore, we consider that fishing is an essential element in our life. We obtained that the stability of ordinary differential system around equilibria, as shown in Figures 2,5 and 6. From Theorem 3.4, by regarding conversion rate γ as bifurcation parameter, Hopf bifurcation of the diffusive system without time delay occurred. What's more, system (3(3) $\begin{array}{l}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}=d_1\mathrm{\Delta }M+M(1-M(t-\tau \left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}=d_2\mathrm{\Delta }N+\gamma \sqrt{M}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)=M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)=\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)\in [0,\pi ]\times [-\tau ,0],\end{array}$(3) ) stayed stability when the time delay $\tau <{\tau }_{\ast }$, as shown in Figure 3, but it lost its stability around $S^{\ast }$. Hopf-bifurcation occurred when time delay crossed its corresponding critical value, as shown in Figure 4. Based on normal form theory and centre manifold theorem for partial differential equations, we obtained the properties of Hopf bifurcation.

In order to make the model more practical, the time delay in maturation of predator should be considered into predator–prey model with schooling behaviour (34) $\begin{array}{rl}\frac{\mathrm{\partial }M}{\mathrm{\partial }t}& =d_1\mathrm{\Delta }M+M(1-M(t-{\tau }_1\left)\right)(M-\text{β})-\alpha \sqrt{M}\sqrt{N},\\ \frac{\mathrm{\partial }N}{\mathrm{\partial }t}& =d_2\mathrm{\Delta }N+\gamma \sqrt{M(t-{\tau }_2)}\sqrt{N}-\delta N-qEN,\\ M_x(0,t)& =M_x(\pi ,t)=N_x(0,t)=N_x(\pi ,t)=0,t\ge 0,\\ M(x,t)& =\varphi (x,t)\ge 0,N(x,t)=\psi (x,t)\ge 0,\\ (x,t)& \in [0,\pi ]\times [-\tau ,0],\end{array}$(34) where ${\tau }_1$ and ${\tau }_2$ are the time delay of the resource limitation of the prey logistic equation and the time delay in maturation of predator, respectively. We leave this work in the future.

Acknowledgments

The authors are grateful to the anonymous reviewers for their constructive comments and important suggestions.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work is supported by the National Natural Science Foundation of China (Grant No. 11661050,11860014), and the HongLiu First-class Disciplines Development Program of Lanzhou University of Technology.