1. Introduction

Competition and predation are two fundamental ecological relationships among species and have been widely studied [1]. Recently, it is been recognized that intraguild predation (IGP), which is a combination of competition and predation, has significant impacts on the distribution, abundance, persistence and evolution of the species involved [2]. As a result, growing attention has been paid to IGP models [3,4,5,6,7,8,9].

The general framework of IGP described below was established by Holt and Polis [5]

$\left\{\begin{array}{lll} \dot{R}(t) = R(\varphi(R)-\rho_1(R,N,P)N-\rho_2(R,N,P)P),\\ \dot{N}(t) = N(e_1\rho_1(R,N,P)R-\rho_3(R,N,P)P-m_1),\\ \dot{P}(t) = P(e_2\rho_2(R,N,P)R+e_3\rho_3(R,N,P)N-m_2),\\ \end{array}\right. \label{Eq1.1}$

(1)

where $R(t), N(t), P(t)$ represent the densities of basal resource, IG prey and IG predator, respectively. The quantities $\rho_2(R,N,P)R$ and $\rho_3(R,N,P)N$ are functional responses of the IG predator to the resource and IG prey, respectively; $\rho_1(R,N,P)R$ is the functional response of the IG prey to the basal resource; and $m_1$ and $m_2$ are density-independent morality rates. The parameters $e_1$ and $e_2$ are the conversion rates of resource consumption into reproduction for the IG prey and IG predator, respectively; the parameter $e_3$ denotes the conversion rate of the IG predator from its consumption of IG prey; $R\varphi(R)$ is recruitment of the basal resource.

Functional response describes how the consumption rate of individual consumers varies with respect to resource density and is often used to model predator-prey interactions. For IGP models, several functional response functions have been studied. For instance, Velazquez $et$ $al.$ [10] and Hsu $et$ $al.$ [11] investigated the case with a linear functional response. Abrams and Fung [12] considered Holling type-II functional response. Verdy and Amarasekare [13] and Freeze $et$ $al.$ [14] investigated Holling type-II and ratio-dependent functional responses, respectively. Kang and Wedekin [15] considered Holling-III functional response.

Note that the reproduction of predator following the consumption of prey is not instantaneous, but rather is mediated by some reaction-time lag required for gestation. Time delay plays an important role in ecology and it can induce very complex dynamical behaviors [16,17,18,19,20,21,22]. For IGP models, it has been shown that a time delay greatly impacts their dynamics [23,24]. In [24], Shu $et$ $al.$ investigated the complex dynamics of the following IGP model

$\label{Eq1.2} \left\{\begin{array}{lll} \dot{R}(t) = rR(t)(1-\frac{R(t)}{K})-c_1R(t)N(t)-c_2R(t)P(t),\\ \dot{N}(t) = e_1c_1R(t-\tau)N(t-\tau)-c_3N(t)P(t)-m_1N(t),\\ \dot{P}(t) = e_2c_2R(t)P(t)+e_3c_3N(t)P(t)-m_2P(t), \end{array}\right.$

(2)

where $r$ is the growth rate of $R$ in the absence of $N$ and $P$, K is the carrying capacity of resource. $c_1$ is the predation rate of IG prey for resource, $c_2$ is the predation rate of IG predator for resource, $c_3$ is the consumption rate of IG predator to IG prey and all other parameters have the same meanings as those in (1).

Note that for each species, individuals tend to migrate towards regions with lower population densities. Hence the species are distributed over space and interact with each other within their spatial domains. To take spatial effects into consideration, reaction diffusion equations become a natural choice [25,26,27,28,29,30,31,32,33,34,35,36]. In this work, we consider a reaction diffusion IGP model with delay and Beddington-DeAngelis functional response.

Suppose $\Omega\subset\mathbb{R}^n$ is a bounded domain with smooth boundary $\partial\Omega$. Let $R(t,x)$, $N(t,x)$, $P(t,x)$ represent the densities of basal resource, IG prey and IG predator at time $t$ and location $x$, respectively. The basal resource is assumed to grow logistically. We assume the basal resource is consumed by the IG prey at a rate $c_1R(t,x)N(t,x)$, and the IG prey is consumed by the IG predator is $c_3N(t,x)P(t,x)$ at time $t$ and location $x$. In this paper, we will assume the functional response takes the Beddington-DeAngelis (B-D) form, i.e., the consumption of the resource by the IG predator is characterized by $\frac{c_2 P(t,x)R(t,x)}{1+a_1R(t,x)+a_2P(t,x)}$. The reproduction of IG prey from consuming the basal resource is $e_1c_1R(t-\tau,x)N(t-\tau,x)$, where the time-lag parameter is introduced in a manner analogous to the treatment in [24]. We further assume the populations cannot cross the boundary of $\Omega$. Our model then reads as

$\label{Eq1.3} \left\{\begin{array}{llll} \frac{\partial R(t,x)}{\partial t} = \widetilde d_1\Delta R+R\left(r(1-\frac{R}{K})-c_1N-\frac{c_2 P}{1+a_1R+a_2P}\right), t>0,x\in\Omega\\ \frac{\partial N(t,x)}{\partial t} = \widetilde d_2\Delta N+e_1c_1 N(t-\tau,x)R(t-\tau,x)-c_3NP-m_1N, t>0, x\in\Omega,\\ \frac{\partial P(t,x)}{\partial t} = \widetilde d_3\Delta P+P\left(\frac{e_2c_2 R}{1+a_1R+a_2P}+e_3c_3N-m_2\right), t>0, x\in\Omega, \\ \frac{\partial R}{\partial \nu} = \frac{\partial N}{\partial \nu} = \frac{\partial P}{\partial \nu} = 0, t>0, x\in\partial\Omega,\\ R(t,x) = \widetilde \phi_1(t,x)\geq 0, (t,x)\in[-\tau,0]\times\Omega,\\ N(t,x) = \widetilde \phi_2(t,x)\geq 0, (t,x)\in[-\tau,0]\times\Omega,\\ P(t,x) = \widetilde \phi_3(t,x)\geq 0, (t,x)\in[-\tau,0]\times\Omega,\\ \end{array}\right.$

(3)

where $\widetilde d_1, \widetilde d_2, \widetilde d_3$ denote the diffusion coefficients of the three species, respectively; $\Delta$ is the Laplacian operator in the $n$ dimensional space, $\nu$ is the outward unit normal vector on $\partial\Omega,$ and the homogeneous Neumann boundary conditions reflect the situation where the population cannot across the boundary of $\Omega.$ The meanings and units of the parameters of model (3) are summarized in Table 1.

For rescalling, we let

$\begin{array}{l} {u_1}(t,x) =\frac{{R(t,x)}}{K},{u_2}(t,x) =\frac{{{c_1}N(t,x)}}{r},{u_3}(t,x) =\frac{{{c_2}P(t,x)}}{r},{\gamma _1} =\frac{{{m_1}}}{r},{\gamma _2} =\frac{{{m_2}}}{r},\\ {\beta _1} =\frac{{{e_1}{c_1}K}}{r},{\beta _2} =\frac{{{e_2}{c_2}K}}{r},\alpha =\frac{{{c_3}}}{{{c_2}}},\beta =\frac{{{e_3}{c_3}}}{{{c_1}}},b ={a_1}K,c =\frac{{{a_2}r}}{{{c_2}}},\\ {d_1} =\frac{{{{\tilde d}_1}}}{{r{L^2}}},{d_2} =\frac{{{{\tilde d}_2}}}{{r{L^2}}},{d_3} =\frac{{{{\tilde d}_3}}}{{r{L^2}}},\hat x =\frac{x}{L},\hat t =tr,\hat \tau =\tau r,\\ {\phi _1}(t,x) =\frac{{{{\tilde \phi }_1}(t,x)}}{K},{\phi _2}(t,x) =\frac{{{c_1}{{\tilde \phi }_2}(t,x)}}{r},{\phi _3}(t,x) =\frac{{{c_2}{{\tilde \phi }_3}(t,x)}}{r}. \end{array}$

Then model (1.3) becomes

$\label{Eq1.4} \left\{\begin{array}{llll} \frac{\partial u_1(t,x)}{\partial t} = d_1\Delta u_1(t,x)+f_1({\bf{u}}, {\bf{v}}),t>0,x\in\Omega,\\ \frac{\partial u_2(t,x)}{\partial t} = d_2\Delta u_2(t,x)+f_2({\bf{u}}, {\bf{v}}), t>0,x\in\Omega,\\ \frac{\partial u_3(t,x)}{\partial t} = d_3\Delta u_3(t,x)+f_3({\bf{u}}, {\bf{v}}), t>0,x\in\Omega,\\ \frac{\partial u_1(t,x)}{\partial \nu} = \frac{\partial u_2(t,x)}{\partial \nu} = \frac{\partial u_3(t,x)}{\partial \nu} = 0, t>0, x\in \partial\Omega,\\ u_1(t,x) = \phi_1(t,x)\geq 0, (t,x)\in[-\tau,0]\times\Omega,\\ u_2(t,x) = \phi_2(t,x)\geq 0,(t,x)\in[-\tau,0]\times\Omega,\\ u_3(t,x) = \phi_3(t,x)\geq 0,(t,x)\in[-\tau,0]\times\Omega,\\ \end{array}\right.$

(4)

where ${\bf{u}} = (u_1(t,x), u_2(t,x), u_3(t,x)), {\bf{v}} = (u_1(t-\tau,x), u_2(t-\tau,x))$ and

$f_1({\bf{u}},{\bf{v}}) = u_1(t,x)\left(1-u_1(t,x)-u_2(t,x)-\frac{u_3(t,x)}{1+bu_1(t,x)+cu_3(t,x)}\right),$

$f_2({\bf{u}}, {\bf{v}}) = \beta_1u_1(t-\tau,x)u_2(t-\tau,x)-\alpha u_2(t,x)u_3(t,x)-\gamma_1u_2(t,x),$

$f_3({\bf{u}}, {\bf{v}}) = u_3(t,x)\left(\frac{\beta_2u_1(t,x)}{1+bu_1(t,x)+cu_3(t,x)}+\beta u_2(t,x)-\gamma_2\right).$

Throughout the paper, we denote $\bar{\Omega} = \Omega\bigcup\partial\Omega$, $D_T = (0,T]\times\Omega, \bar{D}_T = [0,T]\times\bar{\Omega}$, $Q_0 = [-\tau,0]\times\Omega, \bar{Q}_0 = [-\tau,0]\times\bar{\Omega},$ $Q_T = [-\tau,T]\times\Omega, \bar{Q}_T = [-\tau,T]\times\bar{\Omega}.$ Denote by $C^\gamma(D_T)$ the space of Hölder continuous functions in $D_T$ with exponent $\gamma\in(0,1).$ The space of continuous functions in $\bar {D}_T$ is denoted by $C(\bar{D}_T)$. For vector-value functions we use the product spaces

${\cal C}({\bar D_T}) \equiv C({\bar D_T}) \times C({\bar D_T}) \times C({\bar D_T}),{{\cal C}^\gamma }({\bar D_T}) \equiv {C^\gamma }({\bar D_T}) \times {C^\gamma }({\bar D_T}) \times {C^\gamma }({\bar D_T}).$

Denote

$X =\left\{ {{{({u_1},{u_2},{u_3})}^T} \in {H^2}(\Omega ) \times {H^2}(\Omega ) \times {H^2}(\Omega ):\frac{{\partial {u_1}}}{{\partial \nu }} =\frac{{\partial {u_2}}}{{\partial \nu }} =\frac{{\partial {u_3}}}{{\partial \nu }} =0{\mkern 1mu}\;{\rm{on}}\;{\mkern 1mu} \partial \Omega } \right\}$

with the usual inner product $\langle\cdot,\cdot\rangle$.

The rest of the paper is organized as follows. In Section 2, we study the existence and uniqueness of solution of (4) and estimate the solution's priori bounds. In Section 3, we discuss the existence of nonnegative spatially homogeneous steady states. In Section 4, we carry out stability analysis and Hopf bifurcation analysis about the unique positive spatially homogeneous steady state of System (4). Numerical simulations are presented in Section 5 to illustrate the impacts of delay, diffusion and the functional response on the dynamics of our IGP model. We conclude this paper with a brief summary and discussion in Section 6.

---

2. Existence of solution of System (4) and priori bound estimation

Theorem 2.1. Consider System (4), we have the following conclusions.

(ⅰ) Given any initial condition $(\phi_1(t,x),\phi_2(t,x),\phi_3(t,x))\in \mathcal{C}^\gamma(Q_0)$ with

$0\leq \phi_i(t,x)\leq L_i, \;(t,x)\in Q_0,\; i = 1, 2, 3, \label{IC}$

(5)

where $L_i, i = 1, 2, 3$ are positive constants satisfying

$1\leq L_1\leq \frac{\gamma_1}{\beta_1},\; L_2<\frac{\gamma_2}{\beta}, \; L_3\geq\frac{1}{c}\left[\frac{\beta_2L_1}{\gamma_2-\beta L_2}-bL_1-1\right], \label{ICm}$

(6)

System (4) admits a unique solution $(u_1(t,x), u_2(t,x), u_3(t,x))$ satisfying

$0\leq u_i(t,x)\leq L_i,\quad \mathrm{ for }\; t>0,\; x\in\Omega,\;i =1, 2, 3.$

(ⅱ) For any solution $(u_1(t,x), u_2(t,x), u_3(t,x))$ of System (4), it holds true that

$\limsup\limits_{t\rightarrow\infty}u_1(t,x)\leq 1, \, \limsup\limits_{t\rightarrow \infty}\int_\Omega u_2(t,x)dx\leq J_1, \, \limsup\limits_{t\rightarrow \infty}\int_\Omega u_3(t,x)dx\leq J_2$

where $J_1 = (\frac{\beta_1}{4\gamma_1}+\beta_1)|\Omega|$, $J_2 = (\frac{\beta_2}{4\kappa}+\frac{\beta_1\beta}{\alpha\kappa}J_1+\beta_2)|\Omega|$ with $\kappa = \min\{\gamma_1, \gamma_2\}$.

Furthermore, if $d_1 = d_2 = d_3$ and $\tau = 0$, then for any $x\in\bar\Omega,$

$\begin{split} &\limsup\limits_{t\rightarrow\infty}u_2(t,x)\leq \frac{\beta_2\alpha}{4\kappa\beta}+\frac{\beta_1}{\kappa}(\frac{\beta_1}{4\gamma_1}+\beta_1)+\frac{\beta_2\alpha}{\beta},\\ &\limsup\limits_{t\rightarrow\infty}u_3(t,x)\leq\frac{\beta_2}{4\kappa}+\frac{\beta_1\beta}{\alpha\kappa}(\frac{\beta_1}{4\gamma_1}+\beta_1)+\beta_2. \end{split}$

Proof. Note that $f_i({\bf{u, v}})$ is mixed quasi-monotone in a subset $\Lambda\times \Lambda^\ast$ of $\mathbb{R}^3\times \mathbb{R}^2$ for each $i = 1, 2, 3$, we can apply [27,Therorem 2.2] to establish the existence and uniqueness of solutions to System (4). To this end, we first need to construct a pair of coupled upper and lower solutions of System (4), which we denote by $(\widetilde {u}_1, \widetilde {u}_2, \widetilde {u}_3)$ and $(\widehat{u}_1, \widehat{u}_2, \widehat{u}_3)$, respectively. In view of [27,Definition 2.1], the required upper and lower solutions should satisfy the boundary-initial inequalities and the following differential inequalities

$\label{Eq2.3} \left\{ \begin{array}{l} \frac{\partial\widetilde{u}_1}{\partial t}\geq d_1\Delta \widetilde{u}_1+\widetilde{u}_1(1-\widetilde{u}_1-\widehat{u}_2-\frac{\widehat{u}_3}{1+b\widetilde{u}_1+c\widehat{u}_3})\\ \frac{\partial\widetilde{u_2}}{\partial t}\geq d_2\Delta\widetilde{u}_2+\beta_1\widetilde{u}_1\widetilde{u}_2-\alpha\widehat{u}_3\widetilde{u}_2-\gamma_1\widetilde{u}_2,\\ \frac{\partial\widetilde{u}_3}{\partial t}\geq d_3\Delta\widetilde{u}_3+\frac{\beta_2\widetilde{u}_1\widetilde{u}_3}{1+b\widetilde{u}_1+c\widetilde{u}_3} +\beta\widetilde{u}_2\widetilde{u}_3-\gamma_2\widetilde{u}_3 \end{array} \right.$

(7)

and

$\label{Eq2.4} \left\{ \begin{array}{l} \frac{\partial\widehat{u}_1}{\partial t}\leq d_1\Delta \widehat{u}_1+\widehat{u}_1(1-\widehat{u}_1-\widetilde{u}_2-\frac{\widetilde{u}_3}{1+b\widehat{u}_1+c\widetilde{u}_3})\\ \frac{\partial\widehat{u}_2}{\partial t}\leq d_2\Delta\widehat{u}_2+\beta_1\widehat{u}_1\widehat{u}_2-\alpha\widetilde{u}_3\widehat{u}_2-\gamma_1\widehat{u}_2\\ \frac{\partial\widehat{u}_3}{\partial t}\leq d_3\Delta\widehat{u}_3+\frac{\beta_2\widehat{u}_1\widehat{u}_3}{1+b\widehat{u}_1+c\widehat{u}_3}+\beta\widehat{u}_2\widehat{u}_3-\gamma_2\widehat{u}_3. \end{array} \right.$

(8)

Take $(\widehat{u}_1, \widehat{u}_2, \widehat{u}_3) = (0,0,0)$ and $(\widetilde {u}_1, \widetilde {u}_2, \widetilde {u}_3) = (L_1, L_2, L_3)$. Clearly, $\frac{\partial \widetilde u_i}{\partial \nu} = 0\geq 0\geq 0 = \frac{\partial \widehat u_i}{\partial \nu}$. It follows from (5) and (6) that the pair $(\widehat{u}_1, \widehat{u}_2, \widehat{u}_3) = (0,0,0)$ and $(\widetilde {u}_1, \widetilde {u}_2, \widetilde {u}_3) = (L_1, L_2, L_3)$ are coupled upper and lower solutions of System (4). It is easy to check that $f_i({\bf{u, \bf v}})\,(i = 1,2,3)$ satisfy the Lipschitz condition for $0\leq u_i\leq \mathfrak{M}_i, 0\leq v_i\leq L_i, i = 1,2,3,$ and we denote the Lipschitz constants by $K_i, i = 1,2,3.$ By [27,Theorem 2.2], System (4) admits a unique global solution $(u_1(t,x), u_2(t,x), u_3(t,x))$, which satisfies

$(0,0,0)\leq (u_1(t,x), u_2(t,x), u_3(t,x))\leq (L_1, L_2, L_3), \quad \quad t\geq 0,\quad x\in\Omega.$

This completes the proof of (ⅰ).

Next we establish the priori bound of solutions to System (4). To estimate $u_1(t,x)$, we observe that $u_1(t,x)$ satisfies

$\left\{\begin{array}{ll} \frac{\partial u_1(t,x)}{\partial t}\leq d_1\Delta u_1(t,x)+u_1(t,x)(1-u_1(t,x)),& t>0, x\in\Omega,\\ \frac{\partial u_1(t,x)}{\partial n} = 0,& t>0, x\in\partial\Omega. \end{array}\right.$

It follows from the standard comparison principle [37,Lemma 3.4.2] of parabolic equations that $\limsup\limits_{t\rightarrow\infty}u_1(t,x)\leq 1$. Thus for any $\varepsilon>0,$ there exists a $T_1>0$ such that $u_1 (t,x)\leq 1+\varepsilon$ for $t\geq T_1.$

To estimate the priori bounds of $u_2(t,x)$ and $u_3(t,x)$, we set

$U_1(t) = \int_\Omega u_1(t,x)dx, U_2(t) = \int_\Omega u_2(t,x)dx, U_3(t) = \int_\Omega u_3(t,x)dx.$

Then

$\begin{split} \frac{dU_1(t)}{dt} = &\int_\Omega \frac{\partial u_1}{\partial t}dx = \int_\Omega d_1\Delta u_1dx+\int_\Omega [u_1(t,x)(1-u_1(t,x)-u_2(t,x)-\\ &\frac{u_3(t,x)}{1+bu_1(t,x)+cu_3(t,x)})]dx, \end{split}$

$\begin{split} \frac{dU_2(t)}{dt} = &\int_\Omega \frac{\partial u_2}{\partial t}dx = \int_\Omega d_2\Delta u_2dx+\int_\Omega [\beta_1u_1(t-\tau,x)u_2(t-\tau,x)-\\ &\alpha u_2(t,x)u_3(t,x)-\gamma_1u_2(t,x)]dx, \end{split}$

$\begin{array}{l} \frac{{d{U_3}(t)}}{{dt}} = \int_\Omega {\frac{{\partial {u_3}}}{{\partial t}}} dx = \int_\Omega {{d_3}} \Delta {u_3}dx + \\ \quad \quad \quad \quad \quad \quad \quad \quad \quad \int_\Omega {[{u_3}(t,x)} (\frac{{{\beta _2}{u_1}(t,x)}}{{1 + b{u_1}(t,x) + c{u_3}(t,x)}} + \beta {u_2}(t,x) - {\gamma _2})]dx. \end{array}$

From the Neumann boundary conditions, we further obtain

$\begin{array}{l} \frac{d(\beta_1 U_1(t)+U_2(t+\tau))}{dt}\\ = \beta_1\int_\Omega\frac{\partial u_1}{\partial t}dx+\int_\Omega\frac{\partial u_2(t+\tau,x)}{\partial t}dx\\ = \beta_1\int_\Omega(u_1-u_1^2)dx-\int_\Omega\frac{\beta_1u_1u_3}{1+bu_1+cu_3}dx\\ \;\;\;\;-\int_\Omega\alpha u_2(t+\tau,x)u_3(t+\tau,x)dx -\int_\Omega\gamma_1u_2(t+\tau,x)dx\\ \leq\frac{\beta_1}{4}|\Omega|+\gamma_1\beta_1(1+\varepsilon)|\Omega|-\gamma_1(\beta_1U_1(t)+U_2(t+\tau)), \quad t>T_1. \end{array}$

By the comparison principle, we have

$\limsup\limits_{t\rightarrow \infty}(\beta_1U_1(t)+U_2(t+\tau))\leq \frac{\beta_1}{4\gamma_1}|\Omega|+\beta_1|\Omega|\equiv J_1.$

Similarly, there exists $T_2>T_1$ such that $\int_\Omega u_2(t,x)dx = U_2(t)\leq J_1+\varepsilon$ for $t\geq T_2.$ Thus, for $t\geq T_2+\tau,$ we have

$\begin{array}{l} \quad \frac{{d({\beta _2}{U_1}(t) + \frac{\beta }{\alpha }{U_2}(t) + {U_3}(t))}}{{dt}}\\ = {\beta _2}\int_\Omega {({u_1} - u_1^2)} dx + \int_\Omega {\frac{{{\beta _1}\beta }}{\alpha }} {u_1}(t - \tau ,x){u_2}(t - \tau ,x)dx - \frac{{\beta {\gamma _1}}}{\alpha }\int_\Omega {{u_2}} dx - {\gamma _2}\int_\Omega {{u_3}} dx\\ \le \frac{{{\beta _2}}}{4}|\Omega | + \frac{{{\beta _1}\beta }}{\alpha }(1 + \varepsilon )\int_\Omega {{u_2}} (t - \tau ,x)dx - \frac{{\beta {\gamma _1}}}{\alpha }{U_2} - {\gamma _2}{U_3}\\ \le \frac{{{\beta _2}}}{4}|\Omega | + \frac{{{\beta _1}\beta }}{\alpha }(1 + \varepsilon )({J_1} + \varepsilon )|\Omega | + {\beta _2}\kappa (1 + \varepsilon )|\Omega | - \kappa ({\beta _2}{U_1} + \frac{\beta }{\alpha }{U_2} + {U_3}), \end{array}$

where $\kappa = \min\{\gamma_1, \gamma_2\}.$ This implies that

$\limsup\limits_{t\rightarrow \infty}(\beta_2U_1+\frac{\beta}{\alpha}U_2+U_3)\leq \frac{\beta_2}{4\kappa}|\Omega|+\frac{\beta_1\beta}{\alpha\kappa}J_1|\Omega|+\beta_2|\Omega|\equiv J_2.$

and hence $\limsup\limits_{t\rightarrow \infty}\int_\Omega u_3(t,x)dx = \limsup\limits_{t\rightarrow \infty}U_3(t)\leq J_2.$

For the case with $d_1 = d_2 = d_3 = d$ and $\tau = 0$, we can similarly show that for any $\varepsilon>0$, there exists $T_3>T_1$ such that $0\leq u_1(t,x)\leq 1+\varepsilon$ and $0\leq u_2(t,x)\leq \frac{\beta_1}{4\gamma_1}+\beta_1$ for $t>T_3$. Moreover, let $S(t,x) = \beta_2u_1(t,x)+\frac{\beta}{\alpha}u_2(t,x)+u_3(t,x)$. Then

$\left\{\begin{array}{ll} \frac{\partial S}{\partial t} = d\Delta S+\beta_2(u_1-u_1^2)-\beta_2u_1u_2+\frac{\beta\beta_1}{\alpha}u_1u_2-\frac{\beta\gamma_1}{\alpha}u_2-\gamma_2u_3, t>T_3, x\in \Omega,\\ \frac{\partial S}{\partial \nu} = 0, t>T_3, x\in\partial\Omega,\\ S(T_3,x) = \beta_2u_1(T_3,x)+\frac{\beta}{\alpha}u_2(T_3,x)+u_3(T_3,x), x\in\Omega. \end{array}\right.$

Thus for $t>T_3$, we have

$\begin{split} &\beta_2(u_1-u_1^2)-\beta_2u_1u_2+\frac{\beta\beta_1}{\alpha}u_1u_2-\frac{\beta\gamma_1}{\alpha}u_2-\gamma_2u_3\\ \leq &\frac{\beta_2}{4}+\frac{\beta_1\beta}{\alpha}(1+\varepsilon)(\frac{\beta_1}{4\gamma_1}+\beta_1+\varepsilon)+\beta_2\kappa(1+\varepsilon)-\kappa S. \end{split}$

Consider the system

$\left\{ {\begin{array}{*{20}{l}} {\frac{{\partial W}}{{\partial t}} = d\Delta W + \frac{{{\beta _2}}}{4} + \frac{{{\beta _1}\beta }}{\alpha }(1 + \varepsilon )(\frac{{{\beta _1}}}{{4{\gamma _1}}} + {\beta _1} + \varepsilon ) + {\beta _2}\kappa (1 + \varepsilon ) - \kappa W,t > {T_3},x \in \Omega }\\ {\frac{{\partial W}}{{\partial \nu }} = 0,\quad t > {T_3},x \in \partial \Omega ,}\\ {W({T_3},x) = {\beta _2}{u_1}({T_3},x) + \frac{\beta }{\alpha }{u_2}({T_3},x) + {u_3}({T_3},x),x \in \Omega .} \end{array}} \right.$

It follows from [37,Theorem 2.4.6] that the solution $W(t,x)$ satisfies

$\lim\limits_{t\rightarrow\infty}W(t,x) = \frac{\beta_2}{4\kappa}+\frac{\beta_1\beta}{\alpha\kappa} (1+\varepsilon)(\frac{\beta_1}{4\gamma_1}+\beta_1+\varepsilon)+\beta_2(1+\varepsilon).$

The comparison argument implies that

$\limsup\limits_{t\rightarrow \infty}u_2(t,x)\leq \limsup\limits_{t\rightarrow \infty}\frac{\alpha}{\beta}S(t,x) \leq \frac{\beta_2\alpha}{4\kappa\beta}+\frac{\beta_1}{\kappa}(\frac{\beta_1}{4\gamma_1}+\beta_1)+\frac{\beta_2\alpha}{\beta}$

and

$\limsup\limits_{t\rightarrow \infty}u_3(t,x)\leq \limsup\limits_{t\rightarrow \infty}S(t,x) \leq \frac{\beta_2}{4\kappa}+\frac{\beta_1\beta}{\alpha\kappa}(\frac{\beta_1}{4\gamma_1}+\beta_1)+\beta_2.$

This completes the proof.

---

3. Spatially homogeneous steady states of System (4)

Same as in [24], we denote by $\mathcal {R}_i = \frac{\beta_i}{\gamma_i} = \frac{e_ic_iK}{m_i}\,(i = 1,2)$ the reproduction numbers for the IG prey and IG predator, respectively. Consider (4), we easily have the following existence results on trivial and semi-trivial spatially homogeneous steady states.

Proposition 1. (ⅰ) The trivial steady state $E_0 = (0, 0, 0)$ always exists.

(ⅱ) There is a weakly semi-trivial steady state in the absence of IG Prey and IG Predator $E_1 = (1, 0, 0).$

(ⅲ) The IG Prey-only strong semi-trivial steady state $E_{10}: = (\frac{1}{\mathcal{R}_1}, 1-\frac{1}{\mathcal{R}_1}, 0)$ exists if and only if $\mathcal{R}_1>1.$

(ⅳ) The IG Predator-only strong semi-trivial steady state $E_{01}: = (\hat{u}_1, 0, \mathcal{R}_2(1-\hat{u}_1)\hat{u}_1)$ exists if and only if $\mathcal{R}_2>1+b$, where

$\hat{u}_1 = \frac{-(\mathcal{R}_2-c\mathcal{R}_2-b)+\sqrt{(\mathcal{R}_2-c\mathcal{R}_2-b)^2+4c\mathcal{R}_2}}{2c\mathcal{R}_2}.$

System (4) admits a positive steady state $E^\ast: = (u_1^\ast, u_2^\ast, u^\ast_3)$ if $E^\ast$ is a positive solution to the following three equations:

$\label{Eq3.1} 1-u_1-u_2-\frac{u_3}{1+bu_1+cu_3} = 0,$

(9)

$\label{Eq3.2} \beta_1u_1-\alpha u_3-\gamma_1 = 0,$

(10)

$\label{Eq3.3} \frac{\beta_2u_1}{1+bu_1+cu_3}+\beta u_2-\gamma_2 = 0.$

(11)

It follows from (10) that

$\label{Eq3.4} u_1^\ast = \frac{\alpha}{\beta_1}u_3^\ast+\frac{\gamma_1}{\beta_1}.$

(12)

Combining (9), (10) and (12), we obtain

$\label{Eq3.5} {u_3^\ast}^2+pu_3^\ast+q = 0,$

(13)

where

$\begin{split} &p = \frac{c\gamma_2\beta_1^2+b\alpha\beta_1(\gamma_2-\beta)+\alpha\beta_1(\beta-\beta_2)+c\beta\beta_1(\gamma_1-\beta_1)+\beta\beta_1^2+2b\alpha\beta\gamma_1}{\alpha\beta(b\alpha+c\beta_1)},\\ &q = \frac{\beta_1^2(\gamma_2-\beta)+b\beta_1\gamma_1(\gamma_2-\beta)+\beta_1\gamma_1(\beta-\beta_2)+\beta(b\gamma_1^2-\beta_1^2)}{\alpha\beta(b\alpha+c\beta_1)}. \end{split}$

For the distribution of roots of Eq. (13), we have the following results.

Lemma 3.1. (ⅰ) If $p<0, p^2-4q>0$ and $q>0$, then Eq. (13) has two positive roots $u_{3\pm}^\ast = \frac{-p\pm\sqrt{p^2-4q}}{2}.$

(ⅱ) If $p<0$ and $p^2-4q = 0$, then Eq. (13) has a unique positive root $u_{30}^\ast = \frac{-p}{2}.$

(ⅲ) If $q<0$, then Eq. (13) has a unique positive root $u_3^\ast = u_{3+}^\ast.$

Remark 1. Set ${\mathcal{R}}_3 = \frac{\beta}{\gamma_2}$, then $q<0$ provided that $\mathcal{R}_1>b$ and $\mathcal{R}_2>\mathcal{R}_3>1$.

According to Lemma 3.1, the following proposition is valid.

Proposition 2. (ⅰ) When $p^2-4q>0,p<0,q>0$ and

${\mathcal{R}_2}<\min\left\{\frac{(b\alpha+c\beta_1)(-p+\sqrt{p^2-4q})+2(b\gamma_1+\beta_1)}{\alpha(-p+\sqrt{p^2-4q})+2\gamma_1}, \frac{(b\alpha+c\beta_1)(-p-\sqrt{p^2-4q})+2(b\gamma_1+\beta_1)}{\alpha(-p-\sqrt{p^2-4q})+2\gamma_1}\right\}$

System (4) has two positive constant steady states $E_{+}^\ast = (u_{1+}^\ast, u_{2+}^\ast, u_{3+}^\ast)$ and

$\begin{array}{l} E_ - ^ * = (u_{1 - }^ * ,u_{2 - }^ * ,u_{3 - }^ * ),\,where\;u_{1 + }^ * = \frac{{\alpha u_{3 + }^ * + {\gamma _1}}}{{{\beta _1}}},u_{2 + }^ * = \frac{{{\gamma _2}}}{\beta } - \frac{{{\beta _2}}}{\beta } \times \\ \frac{{\alpha u_{3 + }^ * + {\gamma _1}}}{{(b\alpha + c{\beta _1})u_{3 + }^ * + b{\gamma _1} + {\beta _1}}},u_{3 + }^ * = \frac{{ - p + \sqrt {{p^2} - 4q} }}{2}\,\;and\;u_{1 - }^ * = \frac{{\alpha u_{3 - }^ * + {\gamma _1}}}{{{\beta _1}}},u_{2 - }^ * = \frac{{{\gamma _2}}}{\beta } - \frac{{{\beta _2}}}{\beta } \times \\ \frac{{\alpha u_{3 - }^ * + {\gamma _1}}}{{(b\alpha + c{\beta _1})u_{3 - }^ * + b{\gamma _1} + {\beta _1}}},u_{3 - }^ * = \frac{{ - p - \sqrt {{p^2} - 4q} }}{2}. \end{array}$

(ⅱ) When $p^2-4q = 0,p<0$ and ${\mathcal{R}_2}<\frac{-(b\alpha+c\beta_1)p+2(b\gamma_1+\beta_1)}{2\gamma_1-\alpha p}$, $E_{+}^\ast$ and $E_{-}^\ast$ merge, denoted by $E_0^\ast = (u_{10}^\ast, u_{20}^\ast, u_{30}^\ast)$, where $u_{10}^\ast = -\frac{p\alpha}{2\beta_1}+\frac{\gamma_1}{\beta_1}, u_{20}^\ast = \frac{\gamma_2}{\beta}-\frac{\beta_2}{\beta}\frac{2\gamma_1-p\alpha}{2(b\gamma_1+\beta_1)-(b\alpha+c\beta_1)p}, u_{30}^\ast = -\frac{p}{2}.$

(ⅲ) When $q<0$ and $\mathcal{R}_2<\frac{(b\alpha+c\beta_1)(-p+\sqrt{p^2-4q})+2(b\gamma_1+\beta_1)}{\alpha(-p+\sqrt{p^2-4q})+2\gamma_1}$, System (4) has only one positive constant steady state $E^\ast = (u_1^\ast, u_2^\ast, u_3^\ast),$ where $u_1^\ast = u_{1+}^\ast, u_2^\ast = u_{2+}^\ast, u_3^\ast = u_{3+}^\ast.$

Remark 2. There exist parameter values such that Proposition 3.4 holds. For example, choosing $\alpha = 0.68, \beta = 0.9, \beta_1 = 1.9, \beta_2 = 1.8, \gamma_1 = 0.2, \gamma_2 = 0.76,, b = 2.5, c = 10.$ A direct calculation yields only one positive steady state

$E^\ast = (0.1891, 0.7562, 0.2344).$

---

4. Dynamics of System (4)

Let $0 = \mu_1<\mu_2\leq\mu_3\leq\cdots$ be the eigenvalues of $-\Delta$ on $\Omega$ under no-flux boundary conditions, and $E(\mu_i)$ be the eigenspace corresponding to $\mu_i$ with multiplicity $m_i\geq 1, i\in\mathbb{N}\equiv\{1,2,\cdots\}.$ Set ${\bf{X}}_{ij}: = \{{\bf{e}}\cdot\phi_{ij}: {\bf{e}}\in\mathbb{R}^3\}$, where $\{\phi_{ij}\}$ is an orthonormal basis of $E(\mu_i)$ for $j = 1,2,\cdots,\mathrm{dim}E(\mu_i)$. For ${\bf{X}}: = \{{\bf{w}}\in\mathcal{C}^1(\bar\Omega):\frac{\partial w_1}{\partial \nu} = \frac{\partial w_2}{\partial \nu} = \frac{\partial w_3}{\partial \nu} = 0\,\mathrm{on}\, \partial\Omega\}$, we have the following lemma from [37].

Lemma 4.1.

${\bf{X}} = \mathop \oplus \limits_{i = 1}^\infty {{\bf{X}}_i},{\mkern 1mu} {\rm{where}}{\mkern 1mu} {{\bf{X}}_i} = \mathop \oplus \limits_{j = 1}^{{\rm{dim}}E({\mu _i})} {{\bf{X}}_{ij}}.$

---

4.1. Stability of $E_{0}$ and $E_{1}$

In the following, we consider the stability of $E_0$ and $E_1$.

Theorem 4.2.

(ⅰ) The trivial steady state $E_0$ is always unstable.

(ⅱ) The semi-trivial steady state $E_1$ is locally asymptotically stable if $\mathcal{R}_1<1$ and $\mathcal{R}_2<1+b$ and unstable if either $\mathcal{R}_1>1$ or $\mathcal{R}_2>1+b.$

Proof. Linearizing System (4) at a constant steady state ${\bf{u}}^\diamond = (u_1^\diamond, u_2^\diamond, u_3^\diamond)$ gives

$\label{Eq4.1} \frac{\partial {\bf{u}}}{\partial t} = (D\Delta+\hat J_1){\bf{u}}+\hat J_2{\bf{u_\tau}},$

(14)

where $D =$diag$\{d_1, d_2, d_3\}, {\bf{u}} = (u_1(t,x), u_2(t,x), u_3(t,x)), {\bf{u_\tau}} = (u_1(t-\tau,x), u_2(t-\tau,x), u_3(t-\tau,x)$ and

$\begin{split} &\hat J_1 = \left [ \begin{array}{ccc} 1-2u_1^\diamond-u_2^\diamond-\frac{u_3^\diamond+c{u_3^\diamond}^2}{(1+bu_1^\diamond+cu_3^\diamond)^2}&-u_1^\diamond&-\frac{(1+bu_1^\diamond)u_1^\diamond}{(1+bu_1^\diamond+cu_3^\diamond)^2}\\ 0&-\alpha u_3^\diamond-\gamma_1&-\alpha u_2^\diamond\\ \frac{\beta_2(1+cu_3^\diamond)u_3^\diamond}{(1+bu_1^\diamond+cu_3^\diamond)^2}&\beta u_3^\diamond&\frac{\beta_2u_1^\diamond(1+bu_1^\diamond)}{(1+bu_1^\diamond+cu_3^\diamond)^2}+\beta u_2^\diamond-\gamma_2 \end{array}\right],\\ &\hat J_2 = \left [\begin{array}{ccc} 0&0&0\\ \beta_1u_2^\diamond&\beta_1u_1^\diamond&0\\ 0&0&0 \end{array}\right]. \end{split}$

From Lemma 4.1, we know that the eigenvalues of the System (4) is confined on the subspace ${\bf{X}}_i$, and $\lambda$ is an eigenvalue of (14) on ${\bf{X_i}}$ if and only if it is an eigenvalue of the matrix $-\mu_iD+J^\ast$, where $J^\ast = \hat J_1+\hat J_2e^{-\lambda\tau}.$ Then the characteristic equation of (14) is

$\label{Eq4.2} \mathrm{det}(\lambda I_3+\mu_iD-J^\ast) = 0,$

(15)

where $I_3$ stands for the $3\times3$ identity matrix.

(ⅰ) If ${\bf{u}}^\diamond = E_0$, then we obtain the following characteristic equation

$(\lambda+d_1\mu_i-1)(\lambda+d_2\mu_i+\gamma_1)(\lambda+d_3\mu_i+\gamma_2) = 0,$

which gives three sets of eigenvalues, namely, $\lambda_{1i} = -d_1\mu_i+1, \lambda_{2i} = -d_2\mu_i-\gamma_1, \lambda_{3i} = -d_3\mu_i-\gamma_2, i = 1,2,\cdots.$ Clearly, for $i = 1$, $\lambda_{11} = 1>0.$ From [40,Corollary 1.11], the trivial steady state $E_0$ is unstable.

(ⅱ) If ${\bf{u}}^\diamond = E_1$, then we obtain the following characteristic equation

$(\lambda+d_1\mu_i+1)(\lambda+d_2\mu_i-\beta_1e^{-\lambda\tau}+\gamma_1)(\lambda+d_3\mu_i-\frac{\beta_2}{1+b}+\gamma_2) = 0,$

which gives the eigenvalues $\lambda_{1i} = -d_1\mu_i-1, \lambda_{3i} = -d_3\mu_i+\frac{\beta_2}{1+b}-\gamma_2$ and $\lambda_{2i}$ is determined by $\lambda_{2i}+d_2\mu_i+\gamma_1-\beta_1e^{-\lambda_{2i}\tau} = 0.$ Clearly, $\lambda_{1i}<0$ for all $\mu_i.$ If $\mathcal{R}_2<1+b$, we get $\lambda_{3i}<0$ for all $\mu_i$. If $\mathcal{R}_1<1$, then we have $\beta_1<\gamma_1+d_2\mu_i$ for all $\mu_i$. Thus it follows from [24,Lemma 6] that the eigenvalues $\lambda_{2i}$ have negative real parts for all $\mu_i.$ It follows from [40,Corollary 1.11] that the equilibrium $E_1$ is locally asymptotically stable for $\mathcal{R}_1<1$ and $\mathcal{R}_2<1+b$. If $\mathcal{R}_1>1$ then there exists $\mu_1 = 0$ such that $d_2\mu_1+\gamma_1<\beta_1,$ so it follows from the [24,Lemma 6] that at least one of the eigenvalues $\lambda_{2i}$ has a positive real part. If $\mathcal{R}_2>1+b$ then there exists $\mu_1 = 0$ such that $\lambda_{31}>0.$ From [40,Corollary 1.11], the steady state $E_1$ is unstable in either case.

Theorem 4.3. Suppose that $\mathcal{R}_1<1, \mathcal{R}_2<1+b$ and $1-\frac{\gamma_2}{\beta}>\frac{1}{c}$. Then $E_1$ is globally asymptotically stable.

Proof. Let $(\widehat{u}_1, \widehat{u}_2, \widehat{u}_3) = (\varepsilon,0,0)$ and $(\widetilde {u}_1, \widetilde {u}_2, \widetilde {u}_3) = \left( \mathfrak{M}_1, \mathfrak{M}_2, \mathfrak{M}_3\right)$, where $\mathfrak{M}_1 = 1+\varepsilon<\frac{\gamma_1}{\beta_1}, \mathfrak{M}_2 = \frac{\gamma_2}{\beta}-\varepsilon>0$, $\mathfrak{M}_3 = \frac{1}{c}(\frac{\beta_2\mathfrak{M}_1}{\gamma_2-\beta\mathfrak{M}_2}-b\mathfrak{M}_1-1) = \frac{1}{c}\left(\frac{\beta_2(1+\varepsilon)}{\beta\varepsilon}-b(1+\varepsilon)-1\right)$ $>0$ and $\varepsilon$ is arbitrary small positive constant.

We claim that $(\widehat{u}_1, \widehat{u}_2, \widehat{u}_3) = (\varepsilon,0,0)$ and $(\widetilde {u}_1, \widetilde {u}_2, \widetilde {u}_3) = \left( \mathfrak{M}_1, \mathfrak{M}_2, \mathfrak{M}_3\right)$ are also coupled upper and lower solutions of System (4). In fact, it follows from $1-\frac{\gamma_2}{\beta}>\frac{1}{c}$ that $1-\frac{\gamma_2}{\beta}>\frac{\mathfrak{M}_3}{1+c\mathfrak{M}_3}$, and thus we get $\varepsilon(1-\varepsilon-(\frac{\gamma_2}{\beta}-\varepsilon)-\frac{\mathfrak{M}_3}{1+b\varepsilon+c\mathfrak{M}_3})\geq 0.$ Hence, we obtain $\widehat{u}_1(1-\widehat{u}_1-\widetilde{u}_2-\frac{\widetilde u_3}{1+b\widehat{u}_1+c\widetilde{u}_3})\geq 0.$ It is easy to check that $(\widehat{u}_1, \widehat{u}_2, \widehat{u}_3) = (\varepsilon,0,0)$ and $(\widetilde {u}_1, \widetilde {u}_2, \widetilde {u}_3) = \left( \mathfrak{M}_1, \mathfrak{M}_2, \mathfrak{M}_3\right)$ satisfy the differential inequalities in (7) and (8). Thus the claim holds. Define the iterated sequences $(\bar u_1^{(m)}, \bar u_2^{(m)}, \bar u_3^{(m)})$ and $(\underline u_1^{(m)}, \underline u_2^{(m)}, \underline u_3^{(m)})$ as follows:

$\begin{equation*}\begin{split} &\bar u_1^{(m)}=\bar u_1^{(m-1)}+\frac{1}{K_1}\left[\bar u_1^{(m-1)}\left(1-\bar u_1^{(m-1)}-\underline u_2^{(m-1)}-\frac{\underline u_3^{(m-1)}}{1+b\bar u_1^{(m-1)}+c\underline u_3^{(m-1)}}\right)\right],\\ &\underline u_1^{(m)}=\underline u_1^{(m-1)}+\frac{1}{K_1}\left[\underline u_1^{(m-1)}\left(1-\underline u_1^{(m-1)}-\bar u_2^{(m-1)}-\frac{\bar u_3^{(m-1)}}{1+b\underline u_1^{(m-1)}+c\bar u_3^{(m-1)}}\right)\right],\,\\ &\bar u_2^{(m)}=\bar u_2^{(m-1)}+\frac{1}{K_2}\left[\bar u_2^{(m-1)}\left(\beta_1\bar u_1^{(m-1)}-\alpha\underline u_3^{(m-1)}-\gamma_1\right)\right],\\ &\underline u_2^{(m)}=\underline u_2^{(m-1)}+\frac{1}{K_2}\left[\underline u_2^{(m-1)}\left(\beta_1\underline u_1^{(m-1)}-\alpha\bar u_3^{(m-1)}-\gamma_1\right)\right],\\ &\bar u_3^{(m)}=\bar u_3^{(m-1)}+\frac{1}{K_3}\left[\bar u_3^{(m-1)}\left(\frac{\beta_2\bar u_1^{(m-1)}}{1+b\bar u_1^{(m-1)}+c\bar u_3^{(m-1)}}+\beta\bar u_2^{(m-1)}-\gamma_2\right)\right],\\ &\underline u_3^{(m)}=\underline u_3^{(m-1)}+\frac{1}{K_3}\left[\underline u_3^{(m-1)}\left(\frac{\beta_2\underline u_1^{(m-1)}}{1+b\underline u_1^{(m-1)}+c\underline u_3^{(m-1)}}+\beta\underline u_2^{(m-1)}-\gamma_2\right)\right], \end{split}\end{equation*}$

where $m = 1,2,\cdots, (\bar u_1^{(0)}, \bar u_2^{(0)}, \bar u_3^{(0)}) = (\mathfrak{M}_1, \mathfrak{M}_2, \mathfrak{M}_3)$, $(\underline u_1^{(0)}, \underline u_2^{(0)}, \underline u_3^{(0)}) = (\varepsilon, 0, 0)$ and $K_i, i = 1,2,3$ are the Lipschitz constants in Theorem 2.1. It is easy to see that ${\bf{f(u,v)}}\equiv (f_1({\bf{u, v}}), f_2({\bf{u, v}}), f_3({\bf{u, v}}))$ is a $C^1$ function of ${\bf{u, v}}$ and is mixed quasi-monotone in a subset $\Lambda\times\Lambda^\ast$ of $\mathbb{R}^3\times\mathbb{R}^2$. We can deduce from the induction method that

${\bf{\widehat u}}\leq {\bf{\underline u^{(m)}}}\leq{\bf{\underline u^{(m+1)}}} \leq {\bf\bar u^{(m+1)}}\leq {\bf{\bar u^{(m)}}}\leq {\bf\widetilde u}. \label{Eq4.3}$

(16)

It follows from (16) that the limits

$\begin{split} &\lim\limits_{m\rightarrow\infty}\bar u_1^{(m)} = \bar u_1, \lim\limits_{m\rightarrow\infty}\bar u_2^{(m)} = \bar u_2, \lim\limits_{m\rightarrow\infty}\bar u_3^{(m)} = \bar u_3,\\ &\lim\limits_{m\rightarrow\infty}\underline u_1^{(m)} = \underline u_1, \lim\limits_{m\rightarrow\infty}\underline u_2^{(m)} = \underline u_2, \lim\limits_{m\rightarrow\infty}\underline u_3^{(m)} = \underline u_3 \label{Eq4.4} \end{split}$

(17)

exist and satisfy the following equations

$\label{Eq4.5} f_1(\bar u_1, \underline u_2, \underline u_3) = 0, f_2(\bar u_1, \bar u_2, \underline u_3) = 0, f_3(\bar u_1, \bar u_2, \bar u_3) = 0,$

(18)

$\label{Eq4.6} f_1(\underline u_1, \bar u_2, \bar u_3) = 0, f_2(\underline u_1, \underline u_2, \bar u_3) = 0, f_3(\underline u_1, \underline u_2, \underline u_3) = 0,$

(19)

where

${\bf{\widehat u}} = (\widehat u_1, \widehat u_2, \widehat u_3), {\bf{\underline u^{(m)}}} = (\underline u_1^{(m)},\underline u_2^{(m)}, \underline u_3^{(m)}), {\bf{\underline u^{(m+1)}}} = (\underline u_1^{(m+1)},\underline u_2^{(m+1)}, \underline u_3^{(m+1)}),$

${\bf\bar u^{(m+1)}} = (\bar u_1^{(m+1)},\bar u_2^{(m+1)}, \bar u_3^{(m+1)}), {\bf{\bar u^{(m)}}} = (\bar u_1^{(m)},\bar u_2^{(m)}, \bar u_3^{(m)}), {\bf{\widetilde u}} = (\widetilde u_1, \widetilde u_2, \widetilde u_3).$

Since $\underline u_2^{(0)} = 0$ and $\underline u_3^{(0)} = 0,$ we get $\underline u_2^{(m)} = 0$ and $\underline u_3^{(m)} = 0$. It follows from (17) that $\underline u_2 = \underline u_3 = 0.$ Thus, it follows from the first equality of (18) that $\bar u_1 = 1$. Substituting $\bar u_1 = 1$ into the second equality of (18) and noting that $\mathcal{R}_1<1$ yields $\bar u_2 = 0$. In view of the third equality of (18), we have $\bar u_3\left(\frac{\beta_2}{1+b+c\bar u_3}-\gamma_2\right) = 0.$ Since $\mathcal{R}_2<1+b$, we obtain $\frac{\beta_2}{1+b+c\bar u_3}-\gamma_2<0$. This implies that $\bar u_3 = 0.$ Hence, it follows from the first equality of (19) that $\underline u_1 = 1.$ Therefore, we have $(\bar u_1, \bar u_2, \bar u_3) = (1, 0, 0) = (\underline u_1, \underline u_2, \underline u_3).$ It follows from [27,Theorem 3.2] that for any initial function ${\bf\phi} = (\phi_1(t,x), \phi_2(t,x), \phi_3(t,x))$ satisfying ${\bf{\widehat u}} \leq {\bf{\phi}} \leq {\bf{\widetilde u}}$ in $Q_0$, the solution ${\bf{u}}\equiv(u_1(t,x), u_2(t,x), u_3(t,x))$ of System (4) satisfies $\lim\limits_{t\rightarrow\infty}{\bf{u}} = (1,0,0).$ This completes the proof.

---

4.2. Stability of $E_{10}$ and $E_{01}$

Next, we consider the stability of the two strong semi-trivial steady states: $E_{10}$ and $E_{01}.$ For the IG prey-only strong semi-trivial steady state $E_{10}$, we have the following result.

Theorem 4.4. Consider System (4) with $\mathcal{R}_1>1.$

(ⅰ) If $-\frac{\beta_2}{\mathcal{R}_1+b}-\beta (1-\frac{1}{\mathcal{R}_1})+\gamma_2<0$, then $E_{10}$ is unstable.

(ⅱ) If $-\frac{\beta_2}{\mathcal{R}_1+b}-\beta (1-\frac{1}{\mathcal{R}_1})+\gamma_2>0$ and $1<\mathcal{R}_1\leq 3$, then $E_{10}$ is locally asymptotically stable for all $\tau\geq 0.$

(ⅲ) If $-\frac{\beta_2}{\mathcal{R}_1+b}-\beta (1-\frac{1}{\mathcal{R}_1})+\gamma_2>0$ and $\mathcal{R}_1>3$, then there exists $\tau_{00}>0$ such that $E_{10}$ is locally asymptotically stable for $\tau\in[0,\tau_{00})$ and is unstable for $\tau>\tau_{00}.$ Further there exists a sequence of delays, $\{\tau_i^j\}_{j = 0}^{+\infty}$ for $i = 1,2,\cdots, N_1,$ at which $E_{10}$ undergoes Hopf bifurcations. Here, $\tau_{00}$ and $\tau_i^j$ are given in the proof of this theorem.

Proof. For $E_{10}$, the characteristic equation is

$m_1(\lambda)[g_1(\lambda)+h_1(\lambda)e^{-\lambda\tau}] = 0, \label{Eq4.7}$

(20)

with $m_1(\lambda) = \lambda+d_3\mu_i-\frac{\beta_2}{\mathcal{R}_1+b}-\beta (1-\frac{1}{\mathcal{R}_1})+\gamma_2,$ $h_1(\lambda) = -(\lambda+d_1\mu_i+\frac{1}{\mathcal{R}_1})\frac{\beta_1}{\mathcal{R}_1}+\frac{\beta_1}{\mathcal{R}_1}(1-\frac{1}{\mathcal{R}_1}),$ and $g_1(\lambda) = (\lambda+d_1\mu_i+\frac{1}{\mathcal{R}_1})(\lambda+d_2\mu_i+\gamma_1).$ Note that $-\frac{\beta_2}{\mathcal{R}_1+b}-\beta (1-\frac{1}{\mathcal{R}_1})+\gamma_2<0$, there exists $\mu_1 = 0$ such that $-d_3\mu_1+\frac{\beta_2}{\mathcal{R}_1+b}+\beta (1-\frac{1}{\mathcal{R}_1})-\gamma_2>0$ holds. Therefore, $m_1(\lambda)$ has at least one zero with a positive real part and the characteristic Eq. (20) has at least one positive root with a positive real part. The proof of $(ⅰ)$ is complete.

Denote

$\bar{E}_1 = d_1\mu_i+\frac{1}{\mathcal{R}_1}>0, \bar L_1 = d_2\mu_italic>0, \quad \quad \bar J_1 = \gamma_1(1-\frac{1}{\mathcal{R}_1})>0.$

Then we have

$\label{Eq4.8} g_1(\lambda)+h_1(\lambda) = \lambda^2+(\bar E_1+\bar L_1)\lambda+\bar E_1\bar L_1+\bar J_1.$

(21)

Since $\bar E_1+\bar L_1>0$ and $\bar E_1\bar L_1+\bar J_1>0,$ it is easy to see that Eq. (21) has no positive zeros for all $\mu_i.$

Define

$\begin{split} \label{Eq4.9} G(u)\equiv&\mid g_1(i\sqrt{u})\mid^2-\mid h_1(i\sqrt u)\mid^2\\ = &u^2+(\bar L_1^2+2\bar L_1\gamma_1+\bar E_1^2)u+\bar E_1^2\bar L_1^2+2\bar E_1^2\bar L_1\gamma_1-\bar J_1^2+2\bar E_1\gamma_1\bar J_1. \end{split}$

(22)

Since $\bar L_1^2+2\bar L_1\gamma_1+\bar E_1^2$ is positive for all $\mu_i,$ if $\mathcal{F}(\mu_i)\equiv \bar E_1^2\bar L_1^2+2\bar E_1^2\bar L_1\gamma_1-\bar J_1^2+2\bar E_1\gamma_1\bar J_1>0$ then $G(u)$ has no positive zeros.

Next, we discuss the distribution of positive zeros of $G(u).$ Clearly $\mathcal{F}(0) = \frac{\beta_1^2}{\mathcal{R}_1^2}(1-\frac{1}{\mathcal{R}_1})(\frac{3}{\mathcal{R}_1}-1).$ If $1<\mathcal{R}_1\leq 3$, then $\mathcal{F}(0)\geq0.$ Since $\mathcal{F}(\mu_i) = (d_1\mu_i+\frac{1}{\mathcal{R}_1})^2d_2^2\mu_i^2+2(d_1\mu_i+\frac{1}{\mathcal{R}_1})^2d_2\mu_i\gamma_1+\frac{\beta_1}{\mathcal{R}_1}(1-\frac{1}{\mathcal{R}_1})[2(d_1\mu_i+\frac{1}{\mathcal{R}_1})\frac{\beta_1}{\mathcal{R}_1}-\frac{\beta_1}{\mathcal{R}_1}(1-\frac{1}{\mathcal{R}_1})]$ is increasing in $\mu_i,$ we obtain $\mathcal{F}(\mu_i)\geq 0$ for all $\mu_i.$ Thus, $G(u)$ has no positive zeros for all $\mu_i.$ It follows from [24,Lemma 11] that $E_{10}$ is locally asymptotically stable for all $\tau\geq 0.$ This completes the proof of (ⅱ).

If $\mathcal{R}_1>3$, then $\mathcal{F}(0)<0.$ Since $\mathcal{F}(\mu_i)$ is increasing in $\mu_i$ and $\lim\limits_{i\rightarrow\infty}\mathcal{F}(\mu_i) = \infty$, there exists a constant $n\in\mathbb{N}$ such that

$\mathcal{F}(\mu_i)\geq 0 \mbox{ for }\quad italic>N_1, \mbox{ and } \mathcal{F}(\mu_i)<0 \mbox{ for }\quad i\leq N_1.$

This implies that (22) has no positive root for $italic>N_1$, has $N_1$ positive roots for $i\leq N_1$, denoted by $u_1, u_2, \cdots,u_{N_1}$. From (22), we get $u_i = \omega_i^2 = \frac{-Tr_i+\sqrt{Tr_i^2-4\delta_i}}{2},$ where $Tr_i = \bar L_1^2+2\bar L_1\gamma_1+\bar E_1^2, \delta_i = \bar E_1^2\bar L_1^2+2\bar E_1^2\bar L_1\gamma_1-\bar J_1^2+2\bar E_1\bar J_1\gamma_1$. Similar to the argument of [24], we obtain

$\tau_{i}^j = \frac{1}{\omega_{i}}\left\{\arccos\left(\frac{\mathcal{B}_{1i}} {\sqrt{\mathcal{B}_{1i}^2+\mathcal{C}_{1i}^2}}\right)+2j\pi\right\}, \quad i = 1,2,\cdots,N_1, \quad j = 0,1,2,\cdots,$

where

$\begin{split} \mathcal{B}_{1i} = &(\frac{\beta_1}{\mathcal{R}_1})d_1\mu_i+\frac{\beta_1}{\mathcal{R}_1^2}-\frac{\beta_1}{\mathcal{R}_1} (1-\frac{1}{\mathcal{R}_1})((d_1\mu_i+\frac{1}{\mathcal{R}_1})(d_2\mu_i+\gamma_1)-\omega_i^2)+\\ &\frac{\beta_1}{\mathcal{R}_1}\omega_i^2(d_1\mu_i+d_2\mu_i+\frac{1}{\mathcal{R}_1}+\gamma_1), \end{split}$

$\begin{split} \mathcal{C}_{1i} = &-(\frac{\beta_1}{\mathcal{R}_1} d_1\mu_i+\frac{\beta_1}{\mathcal{R}_1^2}-\frac{\beta_1}{\mathcal{R}_1} (1-\frac{1}{\mathcal{R}_1}))(d_1\mu_i+d_2\mu_i+\frac{1}{\mathcal{R}_1}+\gamma_1)\omega_i+\\ &\frac{\beta_1}{\mathcal{R}_1}\omega_i((d_1\mu_i+\frac{1}{\mathcal{R}_1})(d_2\mu_i+\gamma_1)-\omega_i^2). \end{split}$

Denote

$\tau_{00} = \min\limits_{i = 1,2,\cdots,N_1}\{\tau_i^0\}.$

It follows easily from $-\frac{\beta_2}{\mathcal{R}_1+b}-\beta (1-\frac{1}{\mathcal{R}_1})+\gamma_2>0$ that $d_3\mu_i-\frac{\beta_2}{\mathcal{R}_1+b}-\beta (1-\frac{1}{\mathcal{R}_1})+\gamma_2>0$ for all $\mu_i$. Then all zeros of $m_1(\lambda)$ have negative real parts. Since all zeros of $g_1(\lambda)+h_1(\lambda)$ have negative real parts, we conclude that all zeros of Eq. (20) have negative real parts for $\tau = 0$. Since $\tau_{00}$ is the minimum value of $\tau$ so that Eq. (20) has purely imaginary roots, applying Lemma 1.1 of [39], we get $E_{10}$ is locally asymptotically stable for $\tau\in[0,\tau_{00})$ and is unstable for $\tau>\tau_{00}.$ From (22), we obtain $G^\prime(u_i)>0$ for $i = 1,2,\cdots, N_1$. Thus, $E_{10}$ undergoes a sequence of Hopf bifurcations as $\tau$ increases through $\tau_i^j$ for $i = 1,2,\cdots,N_1, j = 0,1,2,\cdots.$ This completes the proof of (ⅲ).

For the IG predator-only strong semi-trivial steady state $E_{01}$, we have the following result.

Theorem 4.5. Consider System (4) with $\mathcal{R}_2>1+b.$

(ⅰ) If $\alpha\mathcal{R}_2\hat u_1(1-\hat u_1)+\gamma_1<\beta_1 \hat u_1$, then $E_{01}$ is unstable.

(ⅱ) If $\alpha\mathcal{R}_2\hat u_1(1-\hat u_1)+\gamma_1>\beta_1 \hat u_1$ and $b<\hat u_1(\mathcal{R}_2+b)$, then $E_{01}$ is locally asymptotically stable for all $\tau\geq 0.$ Here $\hat u_1$ is defined in Proposition 3.1.

Proof. For ${\bf{u}}^\diamond = (\hat{u}_1, 0, \mathcal{R}_2(1-\hat{u}_1)\hat{u}_1)$, we get the characteristic equation

$\label{Eq4.10} m_2(\lambda)g_2(\lambda) = 0,$

(23)

with

$m_2(\lambda) = (\lambda+d_2\mu_i+\alpha\hat u_3+\gamma_1-\beta_1\hat u_1e^{-\lambda\tau}),$

and

$\begin{array}{l} {g_2}(\lambda ) = (\lambda + {d_1}{\mu _i} + {{\hat u}_1} - \frac{{b{{\hat u}_1}{{\hat u}_3}}}{{{{(1 + b{{\hat u}_1} + c{{\hat u}_3})}^2}}})(\lambda + {d_3}{\mu _i} + {\gamma _2} - \frac{{{\beta _2}{{\hat u}_1}(1 + b{{\hat u}_1})}}{{{{(1 + b{{\hat u}_1} + c{{\hat u}_3})}^2}}}) + \\ \quad \quad \quad \quad \frac{{{\beta _2}{{\hat u}_1}{{\hat u}_3}(1 + b{{\hat u}_1})(1 + c{{\hat u}_3})}}{{{{(1 + b{{\hat u}_1} + c{{\hat u}_3})}^4}}}, \end{array}$

where

$\hat{u}_1 = \frac{-(\mathcal{R}_2-c\mathcal{R}_2-b)+\sqrt{(\mathcal{R}_2-c\mathcal{R}_2-b)^2+4c\mathcal{R}_2}}{2c\mathcal{R}_2}, \hat{u}_3 = \mathcal{R}_2(1-\hat{u}_1)\hat{u}_1.$

Since $\alpha\mathcal{R}_2\hat u_1(1-\hat u_1)+\gamma_1<\beta_1 \hat u_1$ and $\hat u_3 = \mathcal{R}_2\hat u_1(1-\hat u_1)$ then there exists $\mu_1 = 0$ such that $d_2\mu_1+\alpha\hat u_3+\gamma_1<\beta_1\hat u_1$ holds. Hence, it follows from [24,Lemma 6] that $m_2(\lambda)$ has at least one zero with a positive real part. The proof of $(ⅰ)$ is complete.

Denote

$\begin{align*} &\bar E_2 = d_1\mu_i+\hat u_1-\frac{b\hat u_1\hat u_3}{(1+b\hat u_1+c\hat u_3)^2},\, \bar L_2 = d_3\mu_i+\gamma_2-\frac{\beta_2\hat u_1(1+b\hat u_1)}{(1+b\hat u_1+c\hat u_3)^2},\\ &\bar J_2 = \frac{\beta_2\hat u_1\hat u_3(1+b\hat u_1)(1+c\hat u_3)}{(1+b\hat u_1+c\hat u_3)^4}. \end{align*}$

Then we have

$g_2(\lambda) = \lambda^2+(\bar E_2+\bar L_2)\lambda+\bar E_2\bar L_2+\bar J_2.$

Since $\bar E_2+\bar L_2 = d_1\mu_i+d_3\mu_i+\hat u_1-\frac{b\hat u_1\hat u_3}{(1+b\hat u_1+c\hat u_3)^2}+\gamma_2-\frac{\beta_2\hat u_1(1+b\hat u_1)}{(1+b\hat u_1+c\hat u_3)^2}$ and $\bar E_2\bar L_2+\bar J_2 = (d_1\mu_i+\hat u_1-\frac{b\hat u_1\hat u_3}{(1+b\hat u_1+c\hat u_3)^2})(d_3\mu_i+\gamma_2-\frac{\beta_2\hat u_1(1+b\hat u_1)}{(1+b\hat u_1+c\hat u_3)^2})+\frac{\beta_2\hat u_1\hat u_3(1+b\hat u_1)(1+c\hat u_3)}{(1+b\hat u_1+c\hat u_3)^2}.$ Obviously, $\gamma_2 = \frac{\beta_2\hat u_1}{1+b\hat u_1+c\hat u_3}>\frac{\beta_2\hat u_1(1+b\hat u_1)}{(1+b\hat u_1+c\hat u_3)^2}$. Since $b<\hat u_1(\mathcal{R}_2+b)$, we have $b\hat u_3<\mathcal{R}_2^2\hat u_1^2$. Noting that $\mathcal{R}_2\hat u_1 = 1+b\hat u_1+c\hat u_3$, we obtain $\frac{b\hat u_3}{(1+b\hat u_1+c\hat u_3)^2}<1$. This implies $\hat u_1>\frac{b\hat u_1\hat u_3}{(1+b\hat u_1+c\hat u_3)^2}$. Thus $g_2(\lambda)$ has no nonnegative zeros for all $\tau\geq 0$. It follows from $\alpha\mathcal{R}_2\hat u_1(1-\hat u_1)+\gamma_1>\beta_1 \hat u_1$ that all roots of $m_2(\lambda)$ have negative real parts. Thus all zeros of Eq. (23) have negative real parts for all $\mu_i$, and hence $E_{01}$ is locally asymptotically stable. This completes the proof of (ⅱ).

---

4.3. Stability of the unique positive spatially homogenous steady state $E^\ast$

---

4.3.1. Stability and Hopf bifurcation

In this subsection, by taking $\tau$ as the bifurcation parameter, we investigate the stability and the Hopf bifurcation near the unique positive spatially homogeneous steady state $E^\ast.$ For this purpose, we always assume

$({H_1})\;q < 0\;{\rm{and}}\;{{\cal R}_2} < \frac{{(b\alpha + c{\beta _1})( - p + \sqrt {{p^2} - 4q} ) + 2(b{\gamma _1} + {\beta _1})}}{{\alpha ( - p + \sqrt {{p^2} - 4q} ) + 2{\gamma _1}}}.$

The above assumption guarantees the uniqueness of the positive spatially homogeneous steady state $E^*$. For the spatially homogeneous positive steady state $E^\ast = (u_1^\ast, u_2^\ast, u_3^\ast)$, the characteristic equation is given below:

$\label{Eq4.11} \lambda^3+b_{2i}\lambda^2+b_{1i}\lambda+b_{0i}+e^{-\lambda\tau}(c_{2}\lambda^2+c_{1i}\lambda+c_{0i}) = 0,$

(24)

where

$\begin{align*} b_{2i} = d_1\mu_i+d_2\mu_i+d_3\mu_i+u_1^\ast-bA_3+\beta_1u_1^\ast+c\beta_2A_3,\\ b_{1i} = (d_1\mu_i+u_1^\ast-bA_3)(d_2\mu_i+\beta_1u_1^\ast+d_3\mu_i+c\beta_2A_3)+\\ \quad \quad (d_2\mu_i+\beta_1u_1^\ast)(d_3\mu_i+c\beta_2A_3)+\alpha\beta u_2^\ast u_3^\ast+A_1A_2\beta_2u_1^\ast u_3^\ast, \end{align*}$

$\begin{array}{l} b_{0i}=(d_1\mu_i+u_1^\ast-bA_3)(d_2\mu_i+\beta_1u_1^\ast)(d_3\mu_i+c\beta_2A_3)+(d_1\mu_i+u_1^\ast-bA_3)\alpha\beta u_2^\ast u_3^\ast\\ \quad \quad -\alpha\beta_2A_2u_1^\ast u_2^\ast u_3^\ast+(d_2\mu_i+\beta_1u_1^\ast)A_1A_2\beta_2u_1^\ast u_3^\ast,\\ c_2=-\beta_1u_1^\ast,\\ c_{1i}=\beta_1u_1^\ast u_2^\ast-\beta_1u_1^\ast(d_1\mu_i+d_3\mu_i+u_1^\ast-bA_3+c\beta_2A_3),\\ c_{0i}=-d_1d_3\beta_1u_1^\ast\mu_i^2+d_3\mu_i\beta_1u_1^\ast u_2^\ast-d_1\mu_ic\beta_2A_3\beta_1u_1^\ast-d_3\mu_i(u_1^\ast-bA_3)\beta_1u_1^\ast+\\ \quad \quad (c\beta_2A_3+A_1\beta u_3^\ast)\beta_1u_1^\ast u_2^\ast-(u_1^\ast-bA_3)c\beta_2A_3\beta_1u_1^\ast-A_1A_2\beta_2 u_1^\ast u_3^\ast\beta_1 u_1^\ast,\\ A_1=\frac{1+bu_1^\ast}{(1+bu_1^\ast+cu_3^\ast)^2}, A_2=\frac{1+cu_3^\ast}{(1+bu_1^\ast+cu_3^\ast)^2}, A_3=\frac{u_1^\ast u_3^\ast}{(1+bu_1^\ast+cu_3^\ast)^2}. \end{array}$

Denote $M_1 = u_1^\ast-bA_3, M_2 = \beta_1u_1^\ast, M_3 = c\beta_2A_3, M_4 = \alpha\beta u_2^\ast u_3^\ast, M_5 = A_1A_2\beta_2u_1^\ast u_3^\ast$, $M_6 = \alpha\beta_2A_2u_1^\ast u_2^\ast u_3^\ast, M_7 = A_1\beta u_2^\ast u_3^\ast.$ Clearly, $M_italic>0$ for $i = 2,3,4,5,6,7.$

Furthermore, we assume that

$({H_2}){\mkern 1mu} \;{M_1} > 0.$

$({H_3}){\mkern 1mu} \;{M_6} - {M_2}{M_7} > 0.$

$({H_4})\;{M_1}{M_4} + {M_2}{M_3}u_2^ * + {M_2}{M_7} - {M_6} > 0.$

$(H_3^\prime )\;{M_6} - {M_2}{M_7} \le 0.$

$(H_4^\prime )\;{M_1}({M_1}{M_3} + {M_5} + u_2^ * {M_2}) + {M_3}({M_1}{M_3} + {M_4} + {M_5}) + {M_6} - {M_2}{M_7} > 0.$

Theorem 4.6. (ⅰ) Assume that $(H_1)-(H_4)$ hold. Then the spatially homogeneous positive steady state $E^\ast$ of System (4) with $\tau = 0$ is locally asymptotically stable.

(ⅱ) Assume that $(H_1)-(H_2)$ and $(H_3^\prime)-(H_4^\prime)$ hold. Then the spatially homogeneous positive steady state $E^\ast$ of System (4) with $\tau = 0$ is locally asymptotically stable.

Proof. When $\tau = 0,$ (24) reduces to the following equation

$\label{Eq4.12} \lambda^3+(b_{2i}+c_2)\lambda^2+(b_{1i}+c_{1i})\lambda+b_{0i}+c_{0i} = 0.$

(25)

Since $M_1>0$, a direct calculation yields

$b_{2i}+c_2 = (d_1+d_2+d_3)\mu_i+M_1+M_3>0.$

$\begin{array}{l} {b_{1i}} + {c_{1i}} = ({d_1}{d_2} + {d_1}{d_3} + {d_2}{d_3})\mu _i^2 + ({M_3}{d_1} + {M_3}{d_2} + {M_1}{d_2} + {M_1}{d_3}){\mu _i} + \\ {M_1}{M_3} + {M_4} + {M_5} + u_2^*{M_2} > 0. \end{array}$

$\begin{split} b_{0i}+c_{0i} = &d_1d_2d_3\mu_i^3+(M_3 d_1d_2+M_1 d_2d_3)\mu_i^2+\\ &(M_2u_2^\ast d_3+M_1M_3d_2+M_4d_1+M_5d_2)\mu_i+M_1M_4+M_2M_3u_2^\ast+\\ &M_2M_7-M_6. \end{split}$

$\begin{equation*}\begin{split} (b_{2i}+c_2)(b_{1i}+c_{1i})-(b_{0i}+c_{0i})=&d_1\mu_i(b_{1i}+c_{1i}-d_2d_3\mu_i^2-M_3d_2\mu_i-M_4)+\\ &d_3\mu_i(b_{1i}+c_{1i}-u_2^\ast M_2)+\\ &d_{2}\mu_i(b_{1i}+c_{1i}-M_1d_3\mu_i-M_1M_3-M_5)+\\ &M_1(b_{1i}+c_{1i}-M_4)+M_3(b_{1i}+c_{1i}-u_2^\ast M_2)+\\ &M_6-M_2M_7. \end{split}\end{equation*}$

Since ${b_{1i}} + {c_{1i}} - {d_2}{d_3}\mu _i^2 - {M_3}{d_2}{\mu _i} - {M_4} > 0,{b_{1i}} + {c_{1i}} - u_2^*{M_2} > 0,{b_{1i}} + {c_{1i}} -$${M_1}{d_3}{\mu _i} - {M_1}{M_3} - {M_5} > 0,$${b_{1i}} + {c_{1i}} - {M_4} > 0,{b_{1i}} + {c_{1i}} - u_2^ * {M_2} > 0$, it follows from $(H_3)$ that $(b_{2i}+c_2)(b_{1i}+c_{1i})-(b_{0i}+c_{0i})>0$ for all $\mu_i$.

It follows from $(H_2)$ and $(H_4)$ that $b_{0i}+c_{0i}>0$ for all $\mu_i$. Thus, by the Routh-Hurwitz stability criterion, all the roots of (25) have negative real parts. This completes the proof of part (ⅰ). Similarly, noting that $(b_{2i}+c_2)(b_{1i}+c_{1i})-(b_{0i}+c_{0i})$ is increasing in $\mu_i$ under $(H_2)$, we can prove part (ⅱ).

Next, we discuss the effect of the delay $\tau\neq 0$ on the stability of the positive steady state $E^\ast$. We first determine critical values of $\tau$ at which a pair of simple purely imaginary eigenvalues appears.

Let $\lambda = \mathrm{i}\omega\,(\omega>0)$ be a root of Eq. (24). Substituting $\lambda = \mathrm{i}\omega$ into Eq. (24) yields

$-\mathrm{i}\omega^3-b_{2i}\omega^2+\mathrm{i}b_{1i}\omega+b_{0i}+(-c_2\omega^2+\mathrm{i}\omega c_{1i}+c_{0i})e^{-\mathrm{i}\omega\tau} = 0,$

which implies that

$\label{Eq4.13} b_{2i}\omega^2-b_{0i} = (c_{0i}-c_2\omega^2)\cos\omega\tau+c_{1i}\omega\sin\omega\tau,$

(26)

$\label{Eq4.14} -\omega^3+b_{1i}\omega = (c_{0i}-c_2\omega^2)\sin\omega\tau-c_{1i}\omega\cos\omega\tau.$

(27)

It follows from (26) and (27) that

$\label{Eq4.15} \omega^6+(b^2_{2i}-2b_{1i}-c_2^2)\omega^4+(b_{1i}^2-2b_{0i}b_{2i}+2c_{0i}c_2-c_{1i}^2)\omega^2+b_{0i}^2-c_{0i}^2 = 0.$

(28)

Let $\omega^2 = s,$ and denote $p_i = b^2_{2i}-2b_{1i}-c_2^2,$ $q_i = b_{1i}^2-2b_{0i}b_{2i}+2c_{0i}c_2-c_{1i}^2, r_i = b_{0i}^2-c_{0i}^2.$ Then (28) is reduced to

$\label{Eq4.16} h(s)\equiv s^3+p_is^2+q_is+r_i = 0.$

(29)

From (24), we get

$\begin{array}{l} \quad \quad {p_i} = (d_1^2 + d_2^2 + d_3^2)\mu _i^2 + 2({M_1}{d_1} + {M_2}{d_2} + {M_3}{d_3}){\mu _i} + \\ \quad \quad \quad \quad M_1^2 + M_3^2 - 2{M_4} - 2{M_5},\\ {b_{0i}} - {c_{0i}} = {d_1}{d_2}{d_3}\mu _i^3 + ({d_1}{d_2}{M_3} + {d_2}{d_3}{M_1} + 2{d_1}{d_3}{M_2})\mu _i^2 + \\ \quad \quad \quad \quad (2{d_1}{M_2}{M_3} + {d_2}{M_1}{M_3} + 2{d_3}{M_1}{M_2} + {d_1}{M_4} + {d_2}{M_5} - {d_3}u_2^ * {M_2}){\mu _i}\\ \quad \quad \quad \quad + 2{M_1}{M_2}{M_3} + 2{M_2}{M_5} + {M_1}{M_4} - u_2^ * {M_2}{M_3} - {M_2}{M_7} - {M_6}. \end{array}$

In the following, we need to seek conditions required for Eq. (29) to have at least one positive root. For this purpose, we further make the following hypotheses: $(H_5)$$2M_1M_2M_3+2M_2M_5+M_1M_4-u_2^\ast M_2M_3-M_2M_7-M_6<0.$

$\begin{array}{l} ({H_6})\;2{d_1}{M_2}{M_3} + {d_2}{M_1}{M_3} + 2{d_3}{M_1}{M_2} + {d_1}{M_4} + {d_2}{M_5} - {d_3}u_2^ * {M_2} \ge 0. \end{array}$

Since $b_{0i}-c_{0i}$ is increasing in $\mu_i$ under $(H_2)$ and $(H_6)$, it follows from $(H_5)$ that there exists $N_2\in\mathbb{N}$ such that

$b_{0i}-c_{0i}<0 \quad \mathrm{for}\quad 1\leq i\leq N_2 \quad \mathrm{and} \quad b_{0i}-c_{0i}\geq 0\quad \mathrm{for}\quad i\geq N_2+1,\quad i\in\mathbb{N}.$

According to the above analysis, we have the following lemma.

Lemma 4.7.(ⅰ) Assume that $(H_2)$, $(H_4)-(H_6)$ hold. Then Eq. (29) has at least one positive root for each $i\in\{1,2,\cdots,N_2\}$.

(ⅱ) Assume that and $(H_2), (H_3^\prime)$ and $(H_5)-(H_6)$ hold. Then Eq. (29) has at least one positive root for each $i\in\{1,2,\cdots,N_2\}$.

Proof. It follows from $(H_2)$ and $(H_4)$ that $b_{0i}+c_{0i}>0$ for all $i\in\mathbb{N}$. Thus $r_i<0$ if and only if $b_{0i}-c_{0i}<0.$ It follows from $(H_2)$, $(H_5)$ and $(H_6)$ that there exists $\mu_i\,(i = 1,2,\cdots,N_2)$ such that $r_i<0.$ Since $\lim\limits_{s\rightarrow\infty}h(s) = \infty$ for fixed $\mu_i$, then (29) has at least one positive root for each $i\in\{1,2,\cdots,N_2\}$. This completes proof of part (ⅰ). Similarly, we can prove part (ⅱ).

Remark 3. From Lemma 4.2, without loss of generality, for each $i, 1\leq i\leq N_2$, we may assume that it has three positive roots, which are denoted by $s_{1,i}, s_{2,i}, s_{3,i}$, respectively. Then for each $i, 1\leq i\leq N_2$, (18) has three positive roots $\omega_{k,i} = \sqrt{s_{k,i}}, k = 1,2,3.$ By (26) and (27), we get

$\cos\omega_{k,i}\tau_{k,i} = \frac{(c_{1i}-b_{2i}c_2)\omega_{k,i}^4+(b_{2i}c_{0i}+b_{0i}c_2-c_{1i}b_{1i})\omega_{k,i}^2-b_{0i}c_{0i}}{(c_{0i}-c_2\omega_{k,i}^2)^2+c_{1i}^2\omega_{k,i}^2}.$

Let

$\label{Eq4.17} \tau_{k,j}^i = \frac{\left(\arccos\left(\frac{(c_{1i}-b_{2i}c_2)\omega_{k,i}^4+(b_{2i}c_{0i}+b_{0i}c_2-c_{1i}b_{1i})\omega_{k,i}^2-b_{0i}c_{0i}}{(c_{0i}-c_2\omega_{k,i}^2)^2+c_{1i}^2\omega_{k,i}^2}\right)+2j\pi\right)}{\omega_{k,i}}$

(30)

for $1\leq i\leq N_2, k = 1,2,3, j = 0,1,2,\cdots.$ Then $\pm\mathrm{i}\omega_{k,i}$ is a pair of purely imaginary roots of (28) with $\tau = \tau_{k,j}^i.$ Define

$\label{Eq4.18} \tau^\ast = \tau_{k_0,0}^{i_0} = \min\limits_{k = {1,2,3},\,i = {1,2,\cdots,N_2}}\tau_{k,0}^i, \quad \quad \omega^\ast = \omega_{k_0,i_0}.$

(31)

Lemma 4.8. Let $\lambda(\tau) = \alpha(\tau)\pm\mathrm{i}\beta(\tau)$ be the roots of Eq. (14) near $\tau = \tau^\ast$ satisfying $\alpha(\tau^\ast) = 0, \beta(\tau^\ast) = \omega^\ast.$ Suppose that $h^\prime((\omega^{\ast})^2)>0.$ Then $\pm\mathrm{i}\omega^\ast$ is a pair of simple purely imaginary roots of Eq. (24). Moreover, the following transversality condition holds:

$\mathrm{sign}\left \{\frac{d(\mathrm{Re}\lambda(\tau))}{d\tau}\right\}_{\tau = \tau^\ast,\lambda = \mathrm{i}\omega^\ast}>0.$

Proof. We denote $P(\lambda) = \lambda^3+b_{2i}\lambda^2+b_{1i}\lambda+b_{0i}, Q(\lambda) = c_{2}\lambda^2+c_{1i}\lambda+c_{0i}.$ Then (24) can be rewritten as

$\label{Eq4.19} P(\lambda)+Q(\lambda)e^{-\lambda\tau} = 0.$

(32)

It is easy to know (26) and (27) are equivalent to the following equations

$\begin{split} &\mathrm{Re}P(\mathrm{i}\omega) = -\mathrm{Re}Q(\mathrm{i}\omega)\cos\omega\tau-\mathrm{Im} Q(\mathrm{i}\omega)\sin\omega\tau,\\ &\mathrm{Im}P(\mathrm{i}\omega) = \mathrm{Re}Q(\mathrm{i}\omega)\sin\omega\tau-\mathrm{Im} Q(\mathrm{i}\omega)\cos\omega\tau. \end{split}$

Thus

$\label{Eq4.20} h(\omega^2) = (\mathrm{Re}P(\mathrm{i}\omega))^2+(\mathrm{Im}P(\mathrm{i}\omega))^2- ((\mathrm{Re}Q(\mathrm{i}\omega))^2+(\mathrm{Im}Q(\mathrm{i}\omega))^2).$

(33)

Differentiating both sides of (33) with respect to $\omega$ yields

$\label{Eq4.21} 2\omega h^\prime(\omega^2) = \mathrm{i}[P^\prime(\mathrm{i}\omega)\bar{P}(\mathrm{i}\omega)-\bar P^\prime(\mathrm{i}\omega)P(\mathrm{i}\omega)-Q^\prime(\mathrm{i}\omega)\bar{Q}(\mathrm{i}\omega)+\bar Q^\prime(\mathrm{i}\omega)Q(\mathrm{i}\omega)].$

(34)

Substituting $\mathrm{i}\omega^\ast$ into (32) yields

$\label{Eq4.22}|P(\mathrm{i}\omega^\ast)| = |Q(\mathrm{i}\omega^\ast)|.$

(35)

If $\mathrm{i}\omega^\ast$ is not simple, then $\mathrm{i}\omega^\ast$ must satisfy $P^\prime(\mathrm{i}\omega^\ast)+[Q^\prime(\mathrm{i}\omega^\ast)-\tau^\ast Q(\mathrm{i}\omega^\ast)]e^{-\mathrm{i}\omega^\ast\tau^\ast} = 0.$ Note that $e^{-\mathrm{i}\omega^\ast\tau^\ast} = -P(\mathrm{i}\omega^\ast)/Q(\mathrm{i}\omega^\ast),$ we obtain $\tau^\ast = \frac{Q^\prime(\mathrm{i}\omega^\ast)}{Q(\mathrm{i}\omega^\ast)}-\frac{P^\prime(\mathrm{i}\omega^\ast)}{P(\omega^\ast)}.$ Using (34) and (35), we have

$\begin{split} \mathrm{Im}\tau^\ast = &\mathrm{Im}\left[\frac{Q^\prime(\mathrm{i}\omega^\ast)}{Q(\mathrm{i}\omega^\ast)}-\frac{P^\prime(\mathrm{i}\omega^\ast)}{P(\mathrm{i}\omega^\ast)}\right] = \mathrm{Im}\left[\frac{Q^\prime(\mathrm{i}\omega^\ast)\bar Q(\mathrm{i}\omega^\ast)}{Q(\mathrm{i}\omega^\ast)\bar Q(\mathrm{i}\omega^\ast)}-\frac{P^\prime(\mathrm{i}\omega^\ast)\bar P(\mathrm{i}\omega^\ast)}{P(\mathrm{i}\omega^\ast)\bar P(\mathrm{i}\omega^\ast)}\right]\\ & = \frac{1}{Q(\mathrm{i}\omega^\ast)\bar Q(\mathrm{i}\omega^\ast)}\mathrm{Im}\left[Q^\prime(\mathrm{i}\omega^\ast)\bar Q(\mathrm{i}\omega^\ast)-P^\prime(\mathrm{i}\omega^\ast)\bar P(\mathrm{i}\omega^\ast)\right]\\ & = \frac{Q^\prime(\mathrm{i}\omega^\ast)\bar Q(\mathrm{i}\omega^\ast)-P^\prime(\mathrm{i}\omega^\ast)\bar P(\mathrm{i}\omega^\ast)-\bar{Q}^\prime(\mathrm{i}\omega^\ast)Q(\mathrm{i}\omega^\ast)+\bar{P}^\prime(\mathrm{i}\omega^\ast) P(\mathrm{i}\omega^\ast)}{2\mathrm{i}|Q(\mathrm{i}\omega^\ast)|^2}\\ & = \frac{\omega^\ast h^\prime((\omega^{\ast})^2)}{|Q(\mathrm{i}\omega^\ast)|^2}. \end{split}$

This is a contradiction. Thus $\pm\mathrm{i}\omega^\ast$ is a pair of simple purely imaginary roots of Eq. (24).

Since $\pm\mathrm{i}\omega^\ast$ are simple purely imaginary roots and

$P^\prime(\mathrm{i}\omega^\ast)+[Q^\prime(\mathrm{i}\omega^\ast)-\tau^\ast Q(\mathrm{i}\omega^\ast)]e^{-\mathrm{i}\omega^\ast\tau^\ast}\neq 0,$

we may consider $\lambda = \lambda(\tau)$ to be a differentiable function. Differentiating (22) with respect to $\tau$ yields

$\frac{d\lambda(\tau)}{d\tau} = \frac{\lambda Q(\lambda)}{P^\prime(\lambda)e^{\lambda\tau}+Q^\prime(\lambda)-\tau Q(\lambda)}.$

Using (32) again, we obtain

$\label{Eq4.23} \left(\frac{d\lambda(\tau)}{d\tau}\right)^{-1} = -\frac{P^\prime(\lambda)}{\lambda P(\lambda)}+\frac{Q^\prime(\lambda)}{\lambda Q(\lambda)}-\frac{\tau}{\lambda}.$

(36)

Thus, from (36), we have

$\begin{split} \mathrm{sign}\left \{\frac{d(\mathrm{Re}\lambda(\tau))}{d\tau}\right\}_{\tau = \tau^\ast,\lambda = \mathrm{i}\omega^\ast} = &\mathrm{sign}\left\{\mathrm{Re}\left[-\frac{P^\prime(\lambda)}{\lambda P(\lambda)}+\frac{Q^\prime(\lambda)}{\lambda Q(\lambda)}-\frac{\tau}{\lambda}\right]\right\}_{\tau = \tau^\ast,\lambda = \mathrm{i}\omega^\ast}\\ = &\mathrm{sign}\left\{\mathrm{Re}\bar\lambda\left[Q^\prime(\lambda)\bar{Q(\lambda)}-P^\prime(\lambda)\bar{P(\lambda)}\right]\right\}_{\lambda = \mathrm{i}\omega^\ast}\\ = &\mathrm{sign}\left\{(\omega^\ast)^2h^\prime((\omega^\ast)^2)\right\}>0. \end{split}$

This completes the proof.

When $\tau\in[0,\tau^\ast)$, we know that (24) has no roots on the imaginary axis. By the eigenvalue theory of [39], the sum of orders of the zeros of Eq. (24) for $\tau\in[0,\tau^\ast)$ is equal to Eq. (25). Then Eq. (24) only has negative real part roots for $\tau\in[0,\tau^\ast)$, which implies that $(u_1^\ast, u_2^\ast, u_3^\ast)$ is locally asymptotically stable for $\tau\in[0,\tau^\ast).$ Combining Theorem 4.6, Lemma 4.7 and 4.8, we arrive the following theorem.

Theorem 4.9. Assume that either $(H_1)-(H_6)$ or $(H_1),(H_2), (H_3^\prime),(H_4^\prime),(H_5), (H_6)$ hold. We have the following results

(ⅰ) The spatially homogeneous positive steady state $E^\ast = (u_1^\ast, u_2^\ast, u_3^\ast)$ of System (4) is locally asymptotically stable for $\tau\in [0,\tau^\ast)$ and unstable when $\tau>\tau^\ast$, where $\tau^\ast$ is difined in (31).

(ⅱ) Furthermore, suppose that $h^\prime((\omega^\ast)^2)>0.$ Then System (4) undergoes Hopf bifurcation at the positive steady state $E^\ast = (u_1^\ast, u_2^\ast, u_3^\ast)$ when $\tau = \tau_{k,j}^i$, i.e., a family of spatially periodic solutions bifurcate from $E^\ast = (u_1^\ast, u_2^\ast, u_3^\ast)$ when $\tau$ crosses through the critical values $\tau_{k,j}^i,$ where $\tau_{k,j}^i$ is defined in (30).

Remark 4. There exist parameter values such that hypotheses $(H_1)-(H_6)$ hold. For example, we choose parameter values $\alpha = 0.7,\beta = 0.9,{\beta _1} = 1.95,{\beta _2} = 1.8,$${\gamma _1} = 0.2,{\gamma _2} = 0.8,b = 2.5,c = 12,{d_1} = 0.25,{d_2} = 0.28,{d_3} = 0.2.$ Using any CAS, it is easy to check that hypotheses $(H_1)-(H_6)$ are satisfied.

---

4.3.2. Properties of Hopf bifurcations

In this section, we investigate the direction of Hopf bifurcation and the stability of bifurcated periodic solutions by using the normal form theory and center manifold reduction. For convenience, for fixed $k\in\{1,2,3\}, j = 0,1,2,\cdots,$ we denote $\tau_{k,j}^i (1\leq i\leq N_2)$ by $\hat\tau$, and denote $\omega_{k,i} (1\leq i\leq N_2)$ by $\hat\omega$.

Let $\tau = \hat\tau+\mu, \mu\in \mathbb{R}$ and $\widetilde u_1(t,\cdot) = u_1(\tau t,\cdot)-u_1^\ast, \widetilde u_2(t,\cdot) = u_2(\tau t,\cdot)-u_2^\ast, \widetilde u_3(t,\cdot) = u_3(\tau t, \cdot)-u_3^\ast$ and $\widetilde U(t) = (\widetilde u_1(t,\cdot), \widetilde u_2(t,\cdot),\widetilde u_3(t,\cdot))$ as $u_1-u_1^\ast, u_2-u_2^\ast, u_3-u_3^\ast$, then drop the tildes for simplicity. System (1.4) can be rewritten in the phase space $\mathcal{C} = C([-1,0], X)$ as

$\dot U(t) = \hat\tau D\triangle U(t)+L(\hat\tau)(U_t)+F(U_t,\mu), \label{Eq4.24}$

(37)

where $D = \mathrm{diag}\{d_1, d_2, d_3\}, L(\delta)(\cdot):\mathcal{C}\rightarrow X$ and $F: \mathcal{C}\times R\rightarrow X$ are given by

$L(\delta)\phi = \delta\left[ \begin{array}{c} (\frac{bu_1^\ast u_3^\ast}{(1+bu_1^\ast+cu_3^\ast)^2}-u_1^\ast)\phi_1(0)-u_1^\ast\phi_2(0)-\frac{(1+bu_1^\ast)u_1^\ast}{(1+bu_1^\ast+cu_3^\ast)^2}\phi_3(0)\\ \beta_1u_2^\ast\phi_1(-1)+\beta_1u_1^\ast\phi_2(-1)-\beta_1u_1^\ast\phi_2(0)-\alpha u_2^\ast\phi_3(0)\\ \frac{\beta_2(1+cu_3^\ast)u_3^\ast}{(1+bu_1^\ast+cu_3^\ast)^2}\phi_1(0)+\beta u_3^\ast\phi_2(0)-\frac{c\beta_2u_1^\ast u_3^\ast}{(1+bu_1^\ast+cu_3^\ast)^2}\phi_3(0) \end{array}\right]$

and

$F(\phi,\mu) = \mu D\Delta\phi(0)+L(\mu)\phi+f(\phi,\mu),$

where

$\begin{split} f(\phi,\mu) = &(\hat\tau+\mu)\times\\ &\left[\begin{array}{l} (\phi_1(0)+u_1^\ast)(-\phi_1(0)-\phi_2(0)-\frac{\phi_3(0)+u_3^\ast}{1+b(\phi_1(0)+u_1^\ast)+c(\phi_3(0)+u_3^\ast)}+\\ \frac{u_3^\ast}{1+bu_1^\ast+cu_3^\ast})-((\frac{bu_1^\ast u_3^\ast}{(1+bu_1^\ast+cu_3^\ast)^2}-u_1^\ast)\phi_1(0)-u_1^\ast\phi_2(0)-\\ \frac{(1+bu_1^\ast)u_1^\ast}{(1+bu_1^\ast+cu_3^\ast)^2}\phi_3(0))\\\\ \beta_1\phi_1(-1)\phi_2(-1)-\alpha \phi_2(0)\phi_3(0)\\\\ (\phi_3(0)+u_3^\ast)(\frac{\beta_2(\phi_1(0)+u_1^\ast)}{1+b(\phi_1(0)+u_1^\ast)+c(\phi_3(0)+u_3^\ast)}+\beta \phi_2(0)-\frac{\beta_2 u_1^\ast}{1+bu_1^\ast+cu_3^\ast})\\ -(\frac{\beta_2(1+cu_3^\ast)u_3^\ast}{(1+bu_1^\ast+cu_3^\ast)^2}\phi_1(0)+\beta u_3^\ast\phi_2(0)-\frac{c\beta_2u_1^\ast u_3^\ast}{(1+bu_1^\ast+cu_3^\ast)^2}\phi_3(0)) \end{array}\right] \end{split}$

for $\phi = (\phi_1, \phi_2, \phi_3)\in\mathcal{C}.$ Let $\mathcal{A}$ to be the infinitesimal generator of the semigroup induced by the solution of the linearized equation of (37)

$\dot U(t) = \hat\tau D\triangle U(t)+L(\hat\tau)(U_t).$

Thus Eq. (37) can be written in the following abstract form

$\frac{dU_t}{dt} = \mathcal{A}U_t+X_0F(U_t,\mu),$

where $X_0(\theta) = \left\{ \begin{array}{ll} 0,\theta\in[-1,0),\\ I, \theta = 0. \end{array}\right.$ Recall that the Banach space decomposition ${\bf{X}} = \bigoplus_{i = 0}^\infty{\bf{X}}_i.$ In view of the Resiz representation theorem, there exists a $3\times 3$ matrix function $\eta(\theta, \hat\tau)\,(-1\leq\theta\leq 0)$, whose entries are bounded variation such that

$-\hat\tau D\mu_i\phi(0)+L(\hat\tau)(\phi) = \int^0_{-1}d[\eta(\theta, \hat\tau)]\phi(\theta),$

for $\phi\in C([-1,0], R^3).$ We can choose

$\begin{equation*}\begin{split} \eta(\theta, \hat\tau)=& \hat\tau\left[ \begin{array}{ccc} -d_1\mu_i+\frac{bu_1^\ast u_3^\ast}{(1+bu_1^\ast+cu_3^\ast)^2}-u_1^\ast&-u_1^\ast&-\frac{(1+bu_1^\ast)u_1^\ast}{(1+bu_1^\ast+cu_3^\ast)^2}\\ 0&-d_2\mu_i-\beta_1u_1^\ast&-\alpha u_2^\ast\\ \frac{\beta_2(1+cu_3^\ast)u_3^\ast}{(1+bu_1^\ast+cu_3^\ast)^2}&\beta u_3^\ast&-d_3\mu_i-\frac{c\beta_2u_1^\ast u_3^\ast}{(1+bu_1^\ast+cu_3^\ast)^2} \end{array}\right]\\ &\times\delta(\theta)+\hat\tau\left[ \begin{array}{ccc} 0&0&0\\ -\beta_1u_2^\ast&-\beta_1 u_1^\ast&0\\ 0&0&0 \end{array}\right]\delta(\theta+1), \end{split}\end{equation*}$

where $\delta$ is a Dirac delta function.

Let us define $C^\ast = C([0,1], \mathbb{R}^{3\ast}),$ where $\mathbb{R}^{3\ast}$ is the $3-$dimensional vector space of row vectors, $\mathcal{A}^\ast$ with domain dense in $C^\ast$ and range in $C^\ast$. Let $\mathcal {P}$ and $\mathcal {P}^\ast$ be the center subspace, that is, the generalized eigenspace of $\mathcal {A}$ and $\mathcal {A}^\ast$ associated with $\Lambda^\ast = \{\mathrm{i}\hat\omega\hat\tau, -\mathrm{i}\hat\omega\hat\tau\}.$

For $\Phi\in C([-1,0], R^3), \Psi\in C([0,1], R^{3\ast})$, we define

$\mathcal{A}(\Phi(\theta)) = \left\{\begin{array}{ll} \frac{d\Phi(\theta)}{d\theta}, \quad \quad \quad \quad \ \theta\in[-1,0),\\ \int_{-1}^0[d\eta(\theta, \hat\tau)]\Phi(\theta), \quad \theta = 0, \end{array}\right.$

and

$\mathcal{A}^\ast(\Psi(s)) = \left\{\begin{array}{ll} -\frac{d\Psi(s)}{ds},&s\in(0,1],\\ \int_{-1}^0\Psi(-\theta)[d\eta(\theta, \hat\tau)],&s = 1. \end{array}\right.$

Then $\mathcal{A}^\ast$ is the formal adjoint of $\mathcal{A}$ under the bilinear pairing

$\begin{align*}\begin{split} \langle\psi, \phi\rangle = &\psi(0)\phi(0)-\int^0_{-1}\int^\theta_0\psi(\xi-\theta)d[\eta(\theta, \hat\tau)]\phi(\xi)d\xi\\ = &\bar\psi(0)\phi(0)+\hat\tau\int^0_{-1}\bar\psi(\xi+1) \left[ \begin{array}{ccc} 0&0&0\\ \beta_1u_2^\ast&\beta_1u_1^\ast&0\\ 0&0&0 \end{array}\right]\phi(\xi)d\xi, \end{split}\end{align*}$

for $\phi\in C([-1,0], R^3), \psi\in C([0,1], R^{3\ast}).$

In view of the definition of the two infinitesimal generators $\mathcal{A}$ and $\mathcal{A}^\ast$, we have the following conclusions.

Lemma 4.10. Let

$\begin{split} \eta_1 = &\frac{(d_3\mu_i+M_3)(G_2^2+H_2^2)-\beta u_3^\ast(G_1G_2+H_1H_2)}{u_3^\ast\beta_2A_2(G_2^2+H_2^2)}+\\ &\mathrm{i}\frac{(\hat\omega(G_2^2+H_2^2)-\beta u_3^\ast(G_2H_1-G_1H_2))}{u_3^\ast\beta_2A_2(G_2^2+H_2^2)}, \end{split}$

$\begin{split} \eta_1^\ast = &\frac{-\alpha u_2^\ast(G_3G_4+H_3H_4)-(d_3\mu_i+M_3)(G_4^2+H_4^2)}{u_1^\ast A_1(G_4^2+H_4^2)}-\\ &\mathrm{i}\frac{\alpha u_2^\ast(G_4H_3-G_3H_4)+\hat\omega(G_4^2+H_4^2)}{u_1^\ast A_1(G_4^2+H_4^2)}, \end{split}$

$\eta_2 = \frac{G_1G_2+H_1H_2+\mathrm{i}(G_2H_1-G_1H_2)}{G_2^2+H_2^2},\eta_2^\ast = \frac{G_3G_4+H_3H_4+\mathrm{i}(G_4H_3-G_3H_4)}{G_4^2+H_4^2},$

where $G_i, H_i\,(i = 1,2,3,4)$ are defined as follows

$G_1 = M_5-\hat\omega+(d_1\mu_i+M_1)(d_3\mu_i+M_3), \quad \quad \quad \ H_1 = \hat\omega(d_1\mu_i+M_1+d_3\mu_i+M_3),$

$G_2 = \beta u_3^\ast d_1\mu_i+\beta M_1u_3^\ast-\beta_2 u_1^\ast u_3^\ast A_2, \quad \quad \quad \quad uad \quad H_2 = \beta u_3^\ast\hat\omega,$

$G_3 = -M_5-(d_3\mu_i+M_3)(d_1\mu_i+M_1)+(\hat\omega)^2,\quad \quad H_3 = -\hat\omega(d_3\mu_i+M_3)-\hat\omega(d_1\mu_i+M_1),$

$G_4 = u_2^\ast A_1M_2\cos\hat\omega\hat\tau+\alpha u_2^\ast(d_1\mu_i+M_1), \quad \quad uad \quad H_4 = -u_2^\ast A_1M_2\sin\hat\omega\hat\tau+\alpha u_2^\ast\hat\omega.$

Then

$p(\theta) = e^{\mathrm{i}\hat\omega\hat\tau\theta}(\eta_1, \eta_2, 1)^T$ is the eigenfunction of $\mathcal{A}$ with respect to $\mathrm{i}\hat\omega\hat\tau$;

$p^\ast (s) = e^{-\mathrm{i}\hat\omega\hat\tau s}(\eta_1^\ast, \eta_2^\ast, 1)$ is the eigenfunction of $\mathcal{A}^\ast$ with respect to $\mathrm{i}\hat\omega\hat\tau.$

Proof. The proof is standard and we omit it here.

Clearly, from Lemma 4.10, we know the center subspace of Eq. (37) is

$\mathcal{P} = \mathrm{span}\{ p(\theta), \bar {p(\theta)}\}, \mathcal{P}^\ast = \mathrm{span}\{p^\ast(s), \bar {p^\ast (s)}\}.$

Then $\mathcal{C}$ can be decomposed as $\mathcal{C} = \mathcal{P}\oplus\mathcal{Q},$ where

$\mathcal{Q} = \{\psi\in\mathcal{C}: \langle\widehat\psi, \psi\rangle = 0 \quad \mathrm{for}\,\mathrm{all}\, \widehat\psi\in \mathcal{P}^\ast \}.$

It follows from Lemma 4.12 that

$\begin{split} \langle p^\ast(s), p(\theta)\rangle = &\bar{ p^\ast(0)}p(0)+\\ &\hat\tau\int_{-1}^0e^{-\mathrm{i}\hat\omega\hat\tau(\xi+1)}(\bar{\eta^\ast_1}, \bar{\eta^\ast_2}, 1)\left[ \begin{array}{ccc} 0&0&0\\ \beta_1u_2^\ast&\beta_1u_1^\ast&0\\ 0&0&0 \end{array}\right](\eta_1, \eta_2, 1)^Te^{\mathrm{i}\hat\omega\hat\tau\xi}d\xi\\ = &\bar{\eta^\ast_1}\eta_1+\bar{\eta^\ast_2}\eta_2+1+\beta_1u_2^\ast\hat\tau e^{-\mathrm{i}\hat\omega\hat\tau}\eta_1\bar{\eta_2^\ast}+\beta_1u_1^\ast\hat\tau e^{-\mathrm{i}\hat\omega\hat\tau}\eta_2\bar{\eta_2^\ast}. \end{split}$

Thus we choose $\mathfrak{D} = 1/(\bar{\eta^\ast_1}\eta_1+\bar{\eta^\ast_2}\eta_2+1+\beta_1u_2^\ast\hat\tau e^{-\mathrm{i}\hat\omega\hat\tau}\eta_1\bar{\eta_2^\ast}+\beta_1u_1^\ast\hat\tau e^{-\mathrm{i}\hat\omega\hat\tau}\eta_2\bar{\eta_2^\ast}).$ Let $\Phi = (p(\theta), \bar{p(\theta)}), \Psi = \mathfrak{D}(p^\ast(s), \bar {p^\ast(s)})^T\equiv (q(s), \bar {q(s)})^T$, then $\langle \Psi, \Phi\rangle = I,$ where $I$ is the identity matrix in $\mathbb{R}^{2\times2}.$

In what follows, in order to determine the bifurcation direction and stability, we compute the coordinates to describe the center manifold $\mathcal{C}_0$. As the formulas to be developed for the bifurcation direction and stability are all relative to $\mu = 0$ only, we set $\mu = 0$ in Eq. (4.24) and obtain a center manifold

$W(z, \bar z)(\theta) = W_{20}(\theta)\frac{z^2}{2}+W_{11}(\theta)z\bar z+W_{02}(\theta)\frac{\bar z^2}{2}+\cdots \label{Eq4.25}$

(38)

with the range in $\mathcal{Q}.$ The flow of Eq. (4.24) on the center manifold can be written as

$U_t = \Phi\cdot(z(t), \bar {z(t)})^T+W(z(t),\bar {z(t)}). \label{Eq4.26}$

(39)

Moreover, from (39), $z$ satisfies

$\begin{split} \dot z(t) = &\frac{d}{dt}\langle q(s), U_t\rangle = \langle q(s), \mathcal{A}U_t\rangle+\langle q(s), X_0F(U_t,0)\rangle\\ = &\mathrm{i}\hat\omega\hat\tau z+\bar {q(0)}f(W(z,\bar z)+2\mathrm{Re}\{zp(\theta)\},0)\\ = &\mathrm{i}\hat\omega\hat\tau z+g(z,\bar z), \label{Eq4.27} \end{split}$

(40)

where

$\begin{split} g(z, \bar z) = g_{20}\frac{z^2}{2}+g_{11}z\bar z+g_{02}\frac{\bar z^2}{2}+g_{21}\frac{z^2\bar z}{2}+\cdot\cdot\cdot. \end{split}$

By the Taylor expansion

$\begin{array}{l} \frac{{v + {v^ * }}}{{1 + b(u + {u^ * }) + c(v + {v^ * })}} = \frac{{{v^ * }}}{{1 + b{u^ * } + c{v^ * }}} - \frac{{b{v^ * }u}}{{{{(1 + b{u^ * } + c{v^ * })}^2}}} + \frac{{(1 + b{u^ * })v}}{{{{(1 + b{u^ * } + c{v^ * })}^2}}} + \\ {C_{11}}{u^2} + {C_{12}}uv + {C_{13}}{v^2} + {C_{14}}{u^2}v + {C_{15}}u{v^2} + \\ {C_{16}}{u^3} + {C_{17}}{v^3} + O(4),\\ \frac{{u + {u^ * }}}{{1 + b(u + {u^ * }) + c(v + {v^ * })}} = \frac{{{u^ * }}}{{1 + b{u^ * } + c{v^ * }}} + \frac{{(1 + c{v^ * })u}}{{{{(1 + b{u^ * } + c{v^ * })}^2}}} - \frac{{c{u^ * }v}}{{{{(1 + b{u^ * } + c{v^ * })}^2}}}\\ + {C_{21}}{u^2} + {C_{22}}uv + {C_{23}}{v^2} + {C_{24}}{u^2}v + {C_{25}}u{v^2} + \\ {C_{26}}{u^3} + {C_{27}}{v^3} + O(4), \end{array}$

(41)

where

$\begin{array}{l} {C_{11}} = \frac{{{b^2}{v^ * }}}{{{{(1 + b{u^ * } + c{v^ * })}^3}}},{C_{12}} = \frac{{bc{v^ * } - b(1 + b{u^ * })}}{{{{(1 + b{u^ * } + c{v^ * })}^3}}},{C_{13}} = - \frac{{c(1 + b{u^ * })}}{{{{(1 + b{u^ * } + c{v^ * })}^3}}},\\ {C_{14}} = \frac{{{b^2}(1 + b{u^ * } - 2c{v^ * })}}{{{{(1 + b{u^ * } + c{v^ * })}^4}}},{C_{15}} = \frac{{bc(2(1 + b{u^ * }) - c{v^ * })}}{{{{(1 + b{u^ * } + c{v^ * })}^4}}},{C_{16}} = \frac{{ - {b^3}{v^ * }}}{{{{(1 + b{u^ * } + c{v^ * })}^4}}}, \end{array}$

$\begin{array}{l} {C_{17}} = \frac{{{c^3}(1 + b{u^ * })}}{{{{(1 + b{u^ * } + c{v^ * })}^4}}},{C_{21}} = - \frac{{b(1 + c{v^ * })}}{{{{(1 + b{u^ * } + c{v^ * })}^3}}},{C_{22}} = \frac{{bc{u^ * } - c(1 + c{v^ * })}}{{{{(1 + b{u^ * } + c{v^ * })}^3}}},\\ {C_{23}} = \frac{{{c^2}{u^ * }}}{{{{(1 + b{u^ * } + c{v^ * })}^3}}},{C_{24}} = \frac{{bc(2(1 + c{v^ * }) - b{u^ * })}}{{{{(1 + b{u^ * } + c{v^ * })}^4}}},{C_{25}} = \frac{{{c^2}(1 + c{v^ * } - 2b{u^ * })}}{{{{(1 + b{u^ * } + c{v^ * })}^4}}},\\ {C_{26}} = \frac{{{b^3}(1 + c{v^ * })}}{{{{(1 + b{u^ * } + c{v^ * })}^4}}},{C_{27}} = - \frac{{{c^3}{u^ * }}}{{{{(1 + b{u^ * } + c{v^ * })}^4}}}. \end{array}$

Noting that (40), we get

$g(z,\bar z) = \bar {\mathfrak{D}}(\bar{\eta_1^\ast}, \bar{\eta_2^\ast}, 1)f(W(z,\bar z)+zp(\theta)+\bar z\bar {p(\theta)},0). \label{Eq4.29}$

(42)

Substituting (38) into (42) and combining (41) yield

$\begin{split} g_{20} = &2\bar{\mathfrak{D}}\hat\tau[\bar{\eta_1^\ast}\left(({bA_3}/{u_1^\ast}-C_{11}u_1^\ast-1)\eta_1^2-\eta_1\eta_2-C_{13}u_1^\ast-(A_1+C_{12}u_1^\ast)\eta_1\right) +\\ &\bar{{\eta_2^\ast}}(\beta_1\eta_1\eta_2e^{-2\mathrm{i}\hat\omega\hat\tau}-\alpha\eta_2) +\bar{{\eta_3^\ast}}(u_3^\ast\beta_2C_{21}\eta_1^2+C_{23}u_3^\ast\beta_2-M_3/u_3^\ast+\\ &(\beta_2A_2+\beta_2C_{22}u_3^\ast)\eta_1+\beta\eta_2)], \end{split}$

$\begin{split} g_{11} = &\bar{\mathfrak{D}}\hat\tau[\bar{\eta_1^\ast}(2({bA_3}/{u_1^\ast}-C_{11}u_1^\ast-1)\eta_1\bar{\eta_1}-(\eta_1\bar{\eta_2}+\eta_2\bar{\eta_1})-2C_{13}u_1^\ast-\\ &(A_1+C_{12}u_1^\ast)(\eta_1+\bar{\eta_1}))+\bar{\eta_2^\ast}(\beta_1(\eta_1\bar{\eta_2}+\bar{\eta_1}\eta_2)-\alpha(\eta_2+\bar{\eta_2}))+\\ &\bar{\eta_3^\ast}(2u_3^\ast\beta_2C_{21}\eta_1\bar{\eta_1}+2(C_{23}u_3^\ast\beta_2-M_3/u_3^\ast)+\\ &(\beta_2A_2+\beta_2C_{22}u_3^\ast)(\eta_1+\bar{\eta_1})+\beta(\eta_2+\bar{\eta_2}))], \end{split}$

$\begin{split} g_{21} = &2\bar{\mathfrak{D}}\hat\tau[\bar{\eta_1^\ast}(2({bA_3}/{u_1^\ast}-C_{11}u_1^\ast-1)(\eta_1W_{11}^{(1)}(0)+\frac{\bar{\eta_1}}{2}W_{20}^{(1)}(0)) -(\frac{W_{20}^{(1)}(0)}{2}\bar{\eta_2}+\\ &W_{11}^{(2)}(0)\eta_1+W_{11}^{(1)}(0)\eta_2+\frac{W_{20}^{(2)}(0)}{2}\bar{\eta_1})-2C_{13}u_1^\ast(W_{11}^{(3)}(0)+\frac{W_{20}^{(3)}(0)}{2})-\\ &(A_1+C_{12}u_1^\ast)(\frac{W_{20}^{(1)}(0)}{2}+W_{11}^{(1)}(0)+W_{11}^{(3)}(0)\eta_1+\frac{W_{20}^{(3)}(0)}{2}\bar{\eta_1})-\\ &3(C_{11}+C_{16}u_1^\ast)\eta_1\eta_1\bar{\eta_1}-(C_{12}+C_{14}u_1^\ast)(\eta_1\eta_1+\eta_1\bar{\eta_1}+\bar{\eta_1}\eta_1)-3C_{17}u_1^\ast-\\ &(C_{13}+C_{15}u_1^\ast)(2\eta_1+\bar{\eta_1}))+\bar{\eta_2^\ast}(\beta_1\eta_1W_{11}^{2}(-1)e^{-\mathrm{i}\hat\omega\hat\tau}+ \beta_1\bar{\eta_1}e^{\mathrm{i}\hat\omega\hat\tau}\frac{W_{20}^{2}(-1)}{2}+\\ &\beta_1\bar{\eta_2}e^{\mathrm{i}\hat\omega\hat\tau}\frac{W_{20}^{1}(-1)}{2}+\beta_1\eta_2W_{11}^{1}(-1)e^{-\mathrm{i}\hat\omega\hat\tau}-\alpha\eta_2W_{11}^{(3)}(0) -\alpha\bar{\eta_2}\frac{W_{20}^{(3)}(0)}{2}-\\ &\alpha\frac{W_{20}^{(2)}(0)}{2}-\alpha W_{11}^{(2)}(0))+\bar{\eta_3^\ast}(2u_3^\ast\beta_2C_{21}(W_{11}^{(1)}(0)\eta_1+\frac{\bar {\eta_1}}{2}W_{20}^{(1)}(0))+\\ &2(C_{23}u_3^\ast\beta_2-M_3/u_3^\ast)(W_{11}^{(3)}(0)+\frac{W_{20}^{(3)}(0)}{2})+(\beta_2A_2+\beta_2C_{22}u_3^\ast)(\frac{W_{20}^{(1)}(0)}{2}+\\ &\eta_1W_{11}^{(3)}(0)+\bar{\eta_1}\frac{W_{20}^{(3)}(0)}{2}+W_{11}^{(1)}(0))+\beta(\eta_2W_{11}^{(3)}(0)+\bar{\eta_2}\frac{W_{20}^{(3)}(0)}{2} + \end{split}$

$\begin{split}&\frac{W_{20}^{(2)}(0)}{2}+W_{11}^{(2)}(0))+3C_{26}u_3^\ast\beta_2\eta_1\eta_1\bar{\eta_1}+3(C_{23}\beta_2+u_3^\ast\beta_2C_{27})+\\ &(C_{21}\beta_2+u_3^\ast\beta_2C_{24})(\eta_1^2+2\eta_1\bar{\eta_1})+(C_{22}\beta_2+u_3^\ast\beta_2C_{25})(2\eta_1+\bar{\eta_1}))], \end{split}$

$g_{02} = 2\bar{g_{20}}\bar{\mathfrak{D}}/\mathfrak{D}.$

Since $W(z(t),\bar{z(t)})$ satisfies the following equation

$\begin{split} \dot W = &\mathcal{A}W+X_0f(\Phi\cdot(z, \bar z)^T+W(z,\bar z),0)-\Phi\Psi(0)f(\Phi\cdot(z,\bar z)^T+W(z,\bar z),0)\\ = &\mathcal{A}W+H_{20}\frac{z^2}{2}+H_{11}z\bar z+H_{02}\frac{\bar z^2}{2}+\cdot\cdot\cdot \label{Eq4.30} \end{split}$

(43)

we obtain

$\left\{\begin{array}{ll} (2\mathrm{i}\hat\omega\hat\tau-\mathcal{A})W_{20} = H_{20},\\ -\mathcal{A}W_{11} = H_{11},\\ (-2\mathrm{i}\hat\omega\hat\tau-\mathcal{A})W_{02} = H_{02}. \end{array}\right. \label{Eq4.31}$

(44)

Since $\mathcal{A}$ has only two eigenvalues $\pm \mathrm{i}\hat\omega\hat\tau$, then Eq. (44) has only a unique solution $W_{ij}.$

We first compute $H_{ij}(\theta), \theta\in[-1, 0].$ From (43), we know that for $-1\leq\theta<0,$

$H(z, \bar z) = -\Phi\Psi(0)f(\Phi\cdot(z,\bar z)^T+W(z,\bar z), 0).$

Therefore, by comparing the coefficients, and notice that

$H(z, \bar z)(0) = f(\Phi\cdot(z, \bar z)^T+W(z,\bar z), 0)-\Phi\Psi(0)f(\Phi\cdot(z,\bar z)^T+W(z,\bar z), 0),$

we obtain

$\begin{split} &H_{20}(\theta) = \\ &\left\{\begin{array}{ll} -(g_{20}p(\theta)+\bar{g_{02}}\bar{p(\theta)}),\quad \quad \theta\in[-1,0),\\ 2\hat\tau\left[\begin{array}{c} ({bA_3}/{u_1^\ast}-C_{11}u_1^\ast-1)\eta_1^2-\eta_1\eta_2-C_{13}u_1^\ast-(A_1+C_{12}u_1^\ast)\eta_1\\ \beta_1\eta_1\eta_2e^{-2\mathrm{i}\hat\omega\hat\tau}-\alpha\eta_2\\ u_3^\ast\beta_2C_{21}\eta_1^2+C_{23}u_3^\ast\beta_2-M_3/u_3^\ast+(\beta_2A_2+\beta_2C_{22}u_3^\ast)\eta_1+\beta\eta_2 \end{array}\right]\\ -[g_{20}p(0)+\bar{g_{02}}\bar {p(0)}], \quad \quad \theta = 0. \end{array}\right. \end{split}$

$\begin{split} &H_{11}(\theta) = \\ &\left\{\begin{array}{ll} -(g_{11}p(\theta)+\bar{g_{11}}\bar{p(\theta)}),\quad \quad \theta\in[-1,0),\\ 2\hat\tau\left[\begin{array}{c} (\frac{bA_3}{u_1^\ast}-C_{11}u_1^\ast-1)\eta_1\bar{\eta_1}-\mathrm{Re}\{\eta_1\bar{\eta_2}\}-C_{13}u_1^\ast -(A_1+C_{12}u_1^\ast)\mathrm{Re}\{\eta_1\}\\ \beta_1\mathrm{Re}\{\eta_1\bar{\eta_2}\}-\alpha\mathrm{Re}\{\eta_2\}\\ u_3^\ast\beta_2C_{21}\eta_1\bar{\eta_1}+C_{23}u_3^\ast\beta_2-\frac{M_3}{u_3^\ast}+(\beta_2A_2+\beta_2C_{22}u_3^\ast)\mathrm{Re}\{\eta_1\}+\beta\mathrm{Re}\{\eta_2\} \end{array}\right]\\ -[g_{11}p(0)+\bar{g_{11}}\bar{p(0)}], \quad \quad \theta = 0. \end{array}\right. \end{split}$

It follows from (44) and the definition of $\mathcal{A}$ that

$\begin{split} &\dot W_{20}(\theta) = 2\mathrm{i}\omega_n\hat\tau W_{20}(\theta)+[g_{20}p(\theta)+\bar{g_{02}}\bar p(\theta)], -1\leq\theta\leq 0,\\ &-\dot W_{11}(\theta) = -[g_{11}p(\theta)+\bar{g_{11}}\bar{p(\theta)}], -1\leq \theta\leq 0. \end{split}$

Noting that $p(\theta) = p(0)e^{\mathrm{i}\hat\omega\hat\tau\theta}, -1\leq\theta\leq 0,$ we have

$\begin{split} &W_{20}(\theta) = [\frac{\mathrm{i}g_{20}}{\hat\omega\hat\tau}p(\theta)+\frac{\mathrm{i}\bar{g_{02}}}{3\hat\omega\hat\tau}\bar{p(\theta)}] +E_1e^{2\mathrm{i}\hat\omega\hat\tau\theta},\\ &W_{11}(\theta) = [\frac{-\mathrm{i}g_{11}}{\hat\omega\hat\tau}p(\theta)+\frac{\mathrm{i}\bar{g_{11}}}{\hat\omega\hat\tau}\bar{p(\theta)}] +E_2. \label{Eq4.32} \end{split}$

(45)

Utilizing the definition of $\mathcal{A}$, (44) and (45) yields

$\begin{split} &E_1 = 2\left[\begin{array}{c} ({bA_3}/{u_1^\ast}-C_{11}u_1^\ast-1)\eta_1^2-\eta_1\eta_2-C_{13}u_1^\ast-(A_1+C_{12}u_1^\ast)\eta_1\\ \beta_1\eta_1\eta_2e^{-2\mathrm{i}\hat\omega\hat\tau}-\alpha\eta_2\\ u_3^\ast\beta_2C_{21}\eta_1^2+C_{23}u_3^\ast\beta_2-M_3/u_3^\ast+(\beta_2A_2+\beta_2C_{22}u_3^\ast)\eta_1+\beta\eta_2 \end{array}\right]\times\\ &\left[ \begin{array}{ccc} 2\mathrm{i}\hat\omega+d_1\mu_i+M_1&u_1^\ast&u_1^\ast A_1\\ -\beta_1u_2^\ast e^{-2\mathrm{i}\hat\omega\hat\tau}&2\mathrm{i}\hat\omega+d_2\mu_i+M_2(1-e^{-2\mathrm{i}\hat\omega\hat\tau})&\alpha u_2^\ast\\ -\beta_2A_2u_3^\ast &-\beta u_3^\ast&2\mathrm{i}\hat\omega+d_3\mu_i+M_3 \end{array}\right ]^{-1}\\ \end{split}$

and

$\begin{array}{l} {E_2} = 2{\left[ {\begin{array}{*{20}{c}} {{d_1}{\mu _i} + u_1^ * - b{A_3}}&{u_1^ * }&{u_1^ * {A_1}}\\ { - {\beta _1}u_2^ * }&{{d_2}{\mu _i}}&{\alpha u_2^ * }\\ { - {\beta _2}{A_2}u_3^ * }&{ - \beta u_3^ * }&{{d_3}{\mu _i} + c{\beta _2}{A_3}} \end{array}} \right]^{ - 1}} \times \\ \quad \left[ {\begin{array}{*{20}{c}} {(\frac{{b{A_3}}}{{u_1^ * }} - {C_{11}}u_1^ * - 1){\eta _1}\overline {{\eta _1}} - {\rm{Re}}\{ {\eta _1}\overline {{\eta _2}} \} - {C_{13}}u_1^ * - ({A_1} + {C_{12}}u_1^ * ){\rm{Re}}\{ {\eta _1}\} }\\ {{\beta _1}{\rm{Re}}\{ {\eta _1}\overline {{\eta _2}} \} - \alpha {\rm{Re}}\{ {\eta _2}\} }\\ {u_3^ * {\beta _2}{C_{21}}{\eta _1}\overline {{\eta _1}} + {C_{23}}u_3^ * {\beta _2} - \frac{{{M_3}}}{{u_3^ * }} + ({\beta _2}{A_2} + {\beta _2}{C_{22}}u_3^ * ){\rm{Re}}\{ {\eta _1}\} + \beta {\rm{Re}}\{ {\eta _2}\} } \end{array}} \right]. \end{array}$

Now, we can compute the following values

$\begin{split} &c_1(0) = \frac{\mathrm{i}}{2\hat\omega\hat\tau}(g_{20}g_{11}-2|g_{11}|^2-\frac{|g_{02}|^2}{3})+\frac{ g_{21}}{2}, \nu_2 = -\frac{\mathrm{Re}(c_1(0))}{\mathrm{Re}(\lambda^\prime(\hat\tau))},\\ &\beta_2 = 2\mathrm{Re}(c_1(0)), T_2 = -\frac{\mathrm{Im}(c_1(0))+\nu_2\mathrm{Im}(\lambda^\prime(\hat\tau))}{\hat\omega\hat\tau}, \end{split}$

which determine the properties of bifurcating periodic solutions at critical value $\hat\tau,$ that is, $\nu_2$ determines the directions of the Hopf bifurcation; $\beta_2$ determines the stability of the bifurcating periodic solutions; $T_2$ determines the period of bifurcating periodic solutions. Moreover, by Hassard [41], we have the following result.

Theorem 4.11. Assume that the conditions of Theorem 4.5 are satisfied, we have

(ⅰ) If $\nu_2>0\,(<0)$, then the direction of the Hopf bifurcation is forward (backward).

(ⅱ)If $\beta_2<0\,(>0)$, then the bifurcating periodic solutions are orbitally stable (unstable).

(ⅲ) If $T_2>0\,(<0)$, then the period of the bifurcating periodic solutions increases (decreases).

---

5. Numerical simulations

---

5.1. The spatiotemporal dynamics in one-dimensional space

In this subsection, we numerically explore the dynamic behavior of System (4) with one-dimensional space, namely $n=1, \Delta=\frac{\partial^2}{\partial x^2}$ and we take $\Omega=(0,\pi).$

For the choice of the parameter values in System (4), we refer to [10,11,12,15] and choose parameter values as follows

$({P_1}):\alpha = 0.7,\beta = 0.9,{\beta _1} = 1.95,{\beta _2} = 1.8,{\gamma _1} = 0.2,{\gamma _2} = 0.8,b = 2.5,c = 5.5,{d_1} =$$0.4,{d_2} = 0.3,{d_3} = 0.2$

and the initial conditions as

$(I{C_1}):{\phi _1}(t,x) = 0.1717 + 0.001\cos x,{\phi _2}(t,x) = 0.7509 + 0.001\cos x,{\phi _3}(t,x) =$$0.1926 + 0.001\cos x.$

With these parameter values, System (4) admits a unique positive spatially homogeneous steady state $E^\ast = (u_1^\ast, u_2^\ast, u_3^\ast)\approx (0.1717, 0.7509, 0.1926).$ It is easy to check that the hypotheses $(H_1)-(H_6)$ hold and $h^\prime((\omega^\ast)^2) = 0.0938>0$. By Theorem 4.9, local Hopf bifurcation occurs at $\tau^\ast\approx0.7859$. We use the forward Euler method to find numerical solutions to System (4) with $\tau = 0.7<\tau^*$ and $\tau = 1.2>\tau^*$, respectively. As illustrated in Figures 1 and 2, when $\tau = 0.7<\tau^*$, solutions of (4) approach the steady state $E^*$, while when $\tau>\tau^*$, sustained oscillations are observed. Calculations give $c_1(0) =-4.2942-30.9399\mathrm{i}, \nu_2 = 83.06, \beta_2 = -8.5884, T_2 = 185.7969.$ Thus the Hopf bifurcation is forward and the bifurcated periodic solutions from $E^\ast$ are stable and the period of bifurcated periodic solution increases in $\tau$ for $\tau>\tau^\ast.$

---

5.2. The spatiotemporal dynamics in two-dimensional space

Consider System (4) with $u_1=u_1(t,x,y), u_2=u_2(t,x,y), u_3=u_3(t,x,y)$ and $\Delta=\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}$. For this purpose, the domain of System (4) is confined to a fixed spatial domain $\Omega=[0, L]\times [0,L]\subset R^2$ with $L=400$. we solve System (4) on a grid with $400\times 400$ sites by a simple Euler method with a time step size of $\delta t=0.01$ and a space step size of $\delta x=\delta y=1$. The discretization of the Laplacian term takes the form

$\begin{split} \Delta u|_{(i,j)} = &\frac{1}{s^2}[A_l(i,j)u(i-1,j)+A_r(i,j)u(i+1,j)+\\ &A_d(i,j)u(i,j-1)+A_u(i,j)u(i,j+1)-4u(i,j)], \end{split}$

where $(i,j)$ denote the lattice sites and $s = 1$ is the lattice constant. The matrix elements of $A_l, A_r, A_d, A_u$ are unity except at the boundary. When $(i,j)$ lies on the left boundary, that is $i = 0,$ we define $A_l(i,j)u(i-1,j)\equiv u(i+1,j)$, which guarantees zero-flux of individuals in the left boundary. Similarly we define $A_r(i,j), A_d(i,j), A_u(i,j)$ such that the zero-flux boundary condition is satisfied.

With the given Neumann boundary conditions, the eigenvalues of $-\Delta$ on $\Omega$ are $\mu_i = \frac{\pi^2}{L^2}(n^2+m^2), n, m\in \mathbb{Z},$ where $\mathbb{Z}$ represents the integer set. In order to discuss the impacts of delay and diffusion on the dynamics of System (4), we will compare the temporal model (that is, System (4) without diffusion) with System (4). We take parameters as

$({P_2}):{\mkern 1mu} \alpha = 0.7,\beta = 0.9,{\beta _1} = 1.95,{\beta _2} = 1.85,{\gamma _1} = 0.2,{\gamma _2} = 0.8,b = 2.5,c$$= 5,{d_1} = 1,{d_2} = 2,{d_3} = 4.$

We consider two types of different initial conditions:

$(IC_2): \phi_1(t,x,y) = u_1^{\ast}+0.001\sin x\sin y, \phi_2(t,x,y) = u_2^{\ast}+0.001\sin x\sin y,$

$\phi_3(t,x,y) = u_3^{\ast}+0.001\sin x \sin y,$

and

$(IC_2^\prime):\left\{\begin{array}{lll} u_1(t,x,y) = 0.1754-\varepsilon_1(x-0.1y-225)(x-0.1y-675),\\ u_2(t,x,y) = 0.7419-\varepsilon_2(x-450)-\varepsilon_3(y-150),\\ u_3(t,x,y) = 0.2029-\varepsilon_4(x-350)-\varepsilon_5(y-200)\\ \end{array}\right.$

for $(t,x,y)\in[-\tau,0]\times\Omega.$ Here $\varepsilon_1 = 2\times 10^{-7}, \varepsilon_2 = 3\times 10^{-5}, \varepsilon_3 = 1.2\times 10^{-3}, \varepsilon_4 = 6\times 10^{-5}, \varepsilon_5 = 3\times 10^{-5}.$

Under $(P_2)$ and $(IC_2)$, it is easy to check that all conditions of Theorem 4.9 are satisfied and there is a unique positive steady state $E^{\ast}\approx(0.1754, 0.7419, 0.2029)$. The corresponding Hopf bifurcation value is computed as $\tau^1_{1,0}\approx 0.8028$.

Figure 3 depicts the population dynamics of the temporal model and the spatiotemporal model at $\tau = 1.$ Both the temporal model and the spatiotemporal model undergo periodic oscillations. However, the temporal model exhibits irregular transient oscillations initially.

If we increase $\tau$ to $\tau = 1.5$ and use initial condition $(IC_2')$, as shown in Figure 4, the temporal model still exhibits regular oscillations, while the spatiotemporal model exhibits irregular oscillations and the calculated Lyapunov exponent is $0.0011>0$ (By the method proposed in [42]), which indicates the occurrence of chaos.

Figure Figure 5 depicts the snapshots of the contour maps of specie $u_1$ for the temporal model and the spatiotemporal model at time $t = 5000$. The temporal model exhibits the spiral wave pattern. However, the spatiotemporal model presents the chaotic wave pattern. To have a better understanding on the evolution process of the spatiotemporal pattern, in Figure Figure 6, we present the snapshots of contour maps of the basal resource $u_1$ at t = 200,500,1000,1200,1500,2500, respectively. As pointed out in [43], the spirals are usually observed under suitable parametric conditions near Turing-Hopf bifurcation threshold. In addition, in the spatially extended system the existence of a stable limit cycle normally results in the formation of chaotic spatiotemporal patterns [44]. As can be seen from Figure Figure 6, as time $t$ increases, an chaotic wave spatial pattern is gradually formed starting from a regular spiral wave pattern.

---

5.2.1. The effect of delay

To explore the impact of delay, in Figure Figure 7, we take the snapshots of the contour maps of specie $u_1$ at time $t = 1500$ for several different values of $\tau$. As can be seen in Figure Figure 7, the time delay can lead to the formation of an irregular spatial pattern from a regular spiral pattern in the whole domain as the time delay increases and surpasses some critical value.

---

5.2.2. The effect of diffusion

As seen from Figure Figure 6, System (4) has a regular spiral wave pattern when $\tau = 1.5$, $d_1 = 1$, $d_2 = 2$ and $d_3 = 4$ at $t = 1500$. To numerically examine how the diffusion affects the pattern, we take the snapshots of contour maps of $u_1$ at $t = 1500$ with several different choices of the diffusion coefficients. As shown in Figure Figure 8, spiral wave pattern emerges firstly when $d_1 = d_2 = d_3 = 0$, then as the three diffusion coefficients change to $d_1 = d_2 = 0.01,d_3 = 0.04$, the spiral wave structure disappears around the center of the spirals wave, with the increase in these diffusion coefficients, it grows steadily, and eventually the chaotic wave pattern dominates the whole domain. Differing from the instability mechanism in Figure Figure 7, the spirals wave loses its stability due to the Doppler effect [45].

---

5.2.3. The impact of the prey saturation constant $b$

Figure 9 demonstrates how the prey saturation constant constant $b$ affects the pattern formation of the basal resource $u_1$ at time $t = 1500$ with $\tau = 1.5:$ when $b$ is small, we observe a pattern with stripes firstly; as the constant $b$ increases, the spiral wave pattern emerges, then it grows steadily, as $b$ goes beyond a certain value, chaotic wave pattern appears.

---

5.2.4. The impact of the predator interference constant $c$

To see how the predator interference constant $c$ influences the spatiotemporal pattern, we numerically simulation (4) with different values of $c$ and plot the snapshots of the contour maps of the basal resource $u_1$ at $t = 1500$ in Figure 10. It is illustrated in Figure Figure 10 that the predator interference constatn $c$ can also lead to the formation of chaotic wave spatial pattern, which can be preceded from the evolution of a regular spiral patterns as the predator interference constant $c$ decreases.

---

6. A summary and discussion

In this work, we have investigated the spatiotemporal dynamics of a diffusive IGP model with delay and the Beddington-DeAngelis functional response. we have established locally asymptotically stability results of the trivial, semi-trivial and strong semi-trivial steady states. In the case that there is a unique positive spatially homogeneous steady state $E^\ast$, we have carried out the Hopf bifurcation analysis. Unlike competition models with monotone response functions ([17]) where delay does not induce sustained oscillations, in our IGP models, delay promotes complex dynamics including bistability, and the emergence of spiral wave pattern and chaotic wave pattern.

Compared with the temporal model in [24], we also observe bistability is possible in System (4). In addition, the diffusion also has impacts on the formation of spatiotemporal patterns as it can change the distribution of characteristic roots of the corresponding characteristic equations, and hence has an important effect on the dynamics for the constant steady state of System (4). This has been illustrated via numerical simulations as well (See Figure 8). Moreover, we have observed that the functional response can also influence the formation of complex patterns. As demonstrated in Figures 9 and 10, the functional responses can also trigger the emergence of spiral wave pattern and chaotic wave spatial pattern.

---

Acknowledgments

The authors were very grateful to two anonymous reviewers' very helpful comments and suggestions. The project was partially supported by the National Natural Science Foundation of China (No.51479215, No. 11401577) and the Natural Sciences and Engineering Research Council of Canada (RGPIN-2015-05686).

---