Dynamical behaviour of an intraguild predator–prey model with prey refuge and hunting cooperation

ABSTRACT

An intraguild predator–prey model including prey refuge and hunting cooperation is investigated in this paper. First, for the corresponding ordinary differential equation model, the existence and stability of all equilibria are given, and the existence of Hopf bifurcation, direction and stability of bifurcating periodic solutions are investigated. Then, for partial differential equation model, the diffusion-driven Turing instability is obtained. What is more, the existence and non-existence of the non-constant positive steady state of the reaction–diffusion model are established by the Leray–Schauder degree theory and some priori estimates. Next, some numerical simulations are performed to support analytical results. The results showed that prey refuge can change the stability of model and even have a stabilizing effect on this model, meanwhile the hunting cooperation can make such model without diffusion unstable, but make such model with diffusion stable. Lastly, a brief conclusion is concluded in the last section.

KEYWORDS:

• Intraguild predation
• prey refuge
• hunting cooperation
• Leray–Schauder degree theory

MATHEMATICS SUBJECT CLASSIFICATIONS:

• 34D20
• 35B35
• 92B05

1. Introduction

The interaction between predator and prey has become increasingly complex due to long-term convolution. To avoid being hunted by predators and prevent their own species from extinction, the prey gradually adopts the collective defence, escape, hiding and other methods to reduce their mortality rate. One of two common ways is the prey hiding from predators through rapid escape, the other is that the prey can avoid predators by establishing various refuges (i.e. spatial refuges) to provide guarantee for population reproduction [20, 21, 32, 51]. Although, the refuge reduces the mortality rate of the prey, the birth rate of the prey will also be reduced due to the lack of resources and harsh environment in refuges. Therefore, the refuge has both positive and negative effects on the prey [20, 25].

In 1946, Crombic did one experiment on two species of beetles, in which the larvae of prey species were placed in glass tubes to avoid being eaten by predators, then obtained the coexistence of the two species. By experiment, he not only verified the existence of natural refuges for the prey but also introduced them into Lotka–Volterra model. Thus Crombic established the first predator–prey model with prey refuge: $\{\begin{array}{rl}& \frac{\mathrm{d}x}{\mathrm{d}t}=ax-b(x-m)y,\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=c(x-m)y-\mathrm{d}y,\end{array}$where m is a fixed number of the prey being protected by the refuge. To easily establish the predator–prey system with prey refuge in mathematics, Taylor [44] divided prey refuges into two types: one is that the prey enters the refuge in a fixed number [6, 23], and the other is that the prey enters the refuge in a fixed proportion [10]. In addition, Ruxton [37] derived a continuous-time model under the assumption that the refuge term is proportional to the density of predator and showed that prey refuge term has a stabilizing effect. Many empirical and theoretical scholars have investigated the consequences of prey refuge and concluded that prey refuge has a stabilizing effect on the considered predator–prey interaction, and increasing refuge can prevent the prey from dying out [9, 11–14, 24, 30, 50].

Recently, Kar [24] investigated a prey–predator model incorporating prey–refuge and independent harvesting: (1) $\{\begin{array}{rl}& \frac{\mathrm{d}x}{\mathrm{d}t}=\alpha x(1-\frac{x}{k})-\frac{\beta (1-m)xy}{1+a(1-m)x}-q_1{E}_{1}x,\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=-\gamma y+\frac{c\beta (1-m)xy}{1+a(1-m)x}-q_2{E}_{2}y,\end{array}$(1) where mx ($m\in [0,1)$) is a refuge protecting of the prey and $(1-m)x$ is the prey available to the predator. Using the harvesting efforts as controls, the optimal harvest policy is obtained by Pontryagin's maximal principle [24].

Based on the model (1(1) $\{\begin{array}{rl}& \frac{\mathrm{d}x}{\mathrm{d}t}=\alpha x(1-\frac{x}{k})-\frac{\beta (1-m)xy}{1+a(1-m)x}-q_1{E}_{1}x,\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=-\gamma y+\frac{c\beta (1-m)xy}{1+a(1-m)x}-q_2{E}_{2}y,\end{array}$(1) ), a three-component model composed of single and double prey species is proposed in the reference [40]: (2) $\{\begin{array}{rl}& \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=\alpha x_1(1-\frac{x_1}{k})-\frac{{\beta }_1(1-m_1)x_1{x}_2}{1+a_1(1-m_1)x_1}-\frac{{\beta }_2(1-m_2)x_1{x}_3}{1+a_2(1-m_2)x_1},\\ & \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=-d_1{x}_2+\frac{c_1{\beta }_1(1-m_1)x_1{x}_2}{1+a_1(1-m_1)x_1}-{\sigma }_1{x}_2{x}_3,\\ & \frac{\mathrm{d}{x}_3}{\mathrm{d}t}=-d_2{x}_3+\frac{c_2{\beta }_2(1-m_2)x_1{x}_3}{1+a_2(1-m_2)x_1}-{\sigma }_2{x}_2{x}_3.\end{array}$(2) Here Holling type II function contains a constant proportion of prey shelters. The authors [40] also considered competition among predators for food (prey) and shelters, and obtained the threshold values of some parameters, the feasibility and stability conditions of the equilibrium, and the ranges of important parameters allowing different types of bifurcations in such model.

Prey species seek the shelter, which makes it difficult for predator species to find food. So a series of adaptive characteristics for predators have been produced in the evolution process, such as sharp teeth, claws, fangs and other tools, using bait pursuit, collective hunting and other ways. For example, carnivores [29], lions [35], wolves [41] and African wild dogs [5] often capture and kill their prey together. Ants [34], spiders [47] and birds [18] also seek out and attack their prey together. These behaviours have caused widespread attention.

Cosner [4] first considered hunting cooperation in predator–prey model and proposed a functional response for the predators who hunt for food in a spatially linear formation and aggregate when they encounter a large group of prey. Berec [2] discussed in detail how cooperative foraging did not always help to sustain the population but can have some destabilizing effect on the underlying systems too. In 2017, Alves and Hilker [1] investigated a prey–predator model with hunting cooperation in a Holling type I function: (3) $\{\begin{array}{rl}& \frac{\mathrm{d}N}{\mathrm{d}t}=rN(1-\frac{N}{K})-\mathrm{\Phi }(N,P)P,\\ & \frac{\mathrm{d}P}{\mathrm{d}t}=e\mathrm{\Phi }(N,P)P-mP,\end{array}$(3) where N and P are the densities of the prey and the predator respectively. $\mathrm{\Phi }(N,P)$ is Holling-type I function, i.e. $\mathrm{\Phi }(N,P)=\lambda N$, here $\lambda >0$ is the attack rate per predator and prey. In the case of hunting cooperation, the function depends on both prey and predator densities. They assumed that cooperative predators benefit from their behaviour, the success of attacks on prey increases with predator density. They [1] represented this assumption in the model (3(3) $\{\begin{array}{rl}& \frac{\mathrm{d}N}{\mathrm{d}t}=rN(1-\frac{N}{K})-\mathrm{\Phi }(N,P)P,\\ & \frac{\mathrm{d}P}{\mathrm{d}t}=e\mathrm{\Phi }(N,P)P-mP,\end{array}$(3) ) by replacing the constant attack rate λ with a density-dependent term: (4) $\mathrm{\Phi }(N,P)=(\lambda +aP)N,$(4) where a>0 is the predator cooperation in hunting and aP is the cooperation term. They showed that hunting cooperation is beneficial to the predator population by increasing the attack rate, identified a scenario in which the hunting cooperation produces Allee effects in predators and allowed the latter to persist when the prey does not sustain them in the absence of hunting cooperation. However, hunting cooperation can turn detrimental to predators when prey density drastically decreases due to the increased predation pressure, which decreases inhalation of the predator. Hunting cooperation can also destabilize the system and promote a sudden collapse of the predator [8, 31].

In addition, Sapkota et al. [39] proposed a three-species food chain model with hunting cooperation among the middle predator as follows: (5) $\{\begin{array}{rl}& \frac{\mathrm{d}{u}_1}{\mathrm{d}t}=u_1(a_1-b_1{u}_1-\frac{{\omega }_0(1+\alpha u_2)u_2}{(1+\alpha u_2)u_1+D_0}),\\ & \frac{\mathrm{d}{u}_2}{\mathrm{d}t}=u_2(-a_2+\frac{{\omega }_0(1+\alpha u_2)u_1}{(1+\alpha u_2)u_1+D_1}-\frac{{\omega }_2{u}_3}{u_2+D_2}),\\ & \frac{\mathrm{d}{u}_3}{\mathrm{d}t}=u_3^2(a_3-\frac{{\omega }_3}{u_2+D_3}),\end{array}$(5) where $u_1,u_2$ and $u_3$ denote the prey, middle predator and top predator respectively. The top predator preys on the middle predator and the middle predator preys on the prey. The hunting cooperation among the middle predator affects interestingly on the numbers of the predator and prey. The linear stability of model (5(5) $\{\begin{array}{rl}& \frac{\mathrm{d}{u}_1}{\mathrm{d}t}=u_1(a_1-b_1{u}_1-\frac{{\omega }_0(1+\alpha u_2)u_2}{(1+\alpha u_2)u_1+D_0}),\\ & \frac{\mathrm{d}{u}_2}{\mathrm{d}t}=u_2(-a_2+\frac{{\omega }_0(1+\alpha u_2)u_1}{(1+\alpha u_2)u_1+D_1}-\frac{{\omega }_2{u}_3}{u_2+D_2}),\\ & \frac{\mathrm{d}{u}_3}{\mathrm{d}t}=u_3^2(a_3-\frac{{\omega }_3}{u_2+D_3}),\end{array}$(5) ) and the two-parameter numerical analysis on hunting cooperation are conducted [39].

Under the same or similar nutrient level in an ecosystem, the relationship between two populations may be both predation and competition [3, 49]. Consequently, Polis, Myers and Holt [36] first proposed the concept of intraguild predation (IGP) system. Intraguild predation is the mixture of competition and predation. That is, IGP involves three species: shared resource, IG prey and IG predator (see Figure 1) and describes an interaction in which two species compete for shared resources and also eat each other. The IG prey feeds on only the shared resource while the IG predator feeds on both the IG prey and the shared resource.

Figure 1. The schematic diagram of intraguild prey–predator model.

The first general model describing IGP was exploited by Holt and Polis [22], which takes the following form: (6) $\{\begin{array}{rlr}& \frac{\mathrm{d}R}{\mathrm{d}t}& =R(\phi (R)-{\rho }_1(R,N,P)N-{\rho }_2(R,N,P)P),\\ & \frac{\mathrm{d}N}{\mathrm{d}t}& =N(e_1{\rho }_1(R,N,P)R-{\rho }_3(R,N,P)P-m_1),\\ & \frac{\mathrm{d}P}{\mathrm{d}t}& =P(e_2{\rho }_2(R,N,P)R+e_3{\rho }_3(R,N,P)N-m_2),\end{array}$(6) where $R(t),N(t)$ and $P(t)$ represent the densities of the shared resource, IG prey and IG predator, respectively. ${\rho }_2(R,N,P)R$ and ${\rho }_3(R,N,P)N$ are functional responses of the IG predator to the shared resource and IG prey, respectively, ${\rho }_1(R,N,P)R$ is the functional response of the IG prey to the shared resource. $m_1$ and $m_2$ are density-independent death rates. The parameters $e_1$ and $e_2$ convert resource consumption into reproduction for the IG prey and IG predator, respectively; $e_3$ is the benefit enjoyed by the IG predator from its consumption of the IG prey; $R\phi (R)$ is recruitment of the shared resource.

In reality, according to Fick's law, species have spatial heterogeneity, and individuals move to areas with lower density to increase the chance of survival. Based on this fact, the importance of the spatial model has been recognized by scholars in various fields and has been continuously investigated in depth. In recent decades, many experiments showed that the spatial models can generate the rich spatiotemporal dynamics (see [7, 15, 16, 26, 33, 38, 42, 43, 45, 48, 52, 53] and references therein). Specially, Tiwari et al. [45] discussed a diffusive predator–prey system including mutually interfering predators, nonlinear harvesting in predators with Crowley–Martin function. They [45] not only obtained the local and global stability of interior equilibrium and the existence and nonexistence of non-constant steady state solutions but also gave the conditions for Turing instabilities of the diffusive system. Han et al. [16] proved the well-posedness of the temporal and spatiotemporal models. For the temporal model, they [16] analysed its temporal dynamics. For the spatiotemporal model, they [16] not only investigated its permanence, stability of non-negative constant steady states and Turing instability, but also gave the existence and non-existence of non-constant positive steady states by Leray Schauder degree theory.

Therefore, based on the above research, it has rich research significance by considering the intraguild predation relationship from a biological point of view. Because the relationship has both predation and competitive, and many scholars have done much work on the intraguild predation model in recent years. Thus we will also consider the intraguild predation model. Both the IG prey and the IG predator feed on shared resources, thus we will consider the effect of refuge on the shared resource to protect the population from extinction. Of course, the IG prey and IG predator can promote them to capture food with greater efficiency through cooperation, thus the hunting cooperation is also considered. The movement among the three populations makes it more realistic for the model to introduce self-diffusion. As a result, taking into account of prey refuge in a Holling-type I functional response of the IG prey to the shared resource and in a Holling-type II functional response of the IG predator to the shared resource, hunting cooperation in a Holling-type I functional response of the IG predator to the IG prey, we will investigate an intraguild predation model as follows: (7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }t}& =d_{11}\mathrm{\Delta }X+rX-a{X}^2-d_{1}X-{\alpha }_1(1-m_1)XY\\ & -\frac{{\alpha }_2(1-m_2)XZ}{1+b(1-m_2)X},x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }t}& =d_{22}\mathrm{\Delta }Y+{\beta }_1(1-m_1)XY-({\alpha }_3+cZ)YZ-d_{2}Y,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }t}& =d_{33}\mathrm{\Delta }Z+\frac{{\beta }_2(1-m_2)XZ}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)YZ-d_{3}Z,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }\nu }& =\frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }\nu }=0,x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ X(x,0)& =X_0(x)\ge 0,Y(x,0)=Y_0(x)\ge 0,Z(x,0)=Z_0(x)\ge 0,x\in \mathrm{\Omega },\end{array}\end{array}$(7) where $X(x,t),Y(x,t)$ and $Z(x,t)$ represent the density of the shared resource, IG prey and IG predator at the position x and time t, respectively. $\mathrm{\Omega }\subset {\mathbb{R}}^N(N\le 3)$ is a fixed spatial bounded domain with smooth boundary $\mathrm{\partial }\mathrm{\Omega }$ and ν is the outward unit normal vector of boundary $\mathrm{\partial }\mathrm{\Omega }$. Generally, $d_{11},d_{22}$ and $d_{33}$ respectively are the self-diffusion rates of the three species; r and a respectively denote the birth rate and intraspecific competition rate of the shared resource. $d_1,d_2$ and $d_3$ are natural death rates of the three species; ${\alpha }_1,\frac{{\alpha }_2}{b}$ and ${\alpha }_3$ are maximum numbers of different prey that can be eaten by the respective predator per unit time; ${\beta }_1,{\beta }_2$ and ${\beta }_3$ are a conversion rate of different prey captured by respective predator. $m_1\in [0,1)$ and $m_2\in [0,1)$ are fixed proportion of the shared resource by refuge, the remaining of the shared resource $(1-m_1)X$ and $(1-m_2)X$ are available to IG prey and IG predator, respectively. c is the same meaning as a of (4(4) $\mathrm{\Phi }(N,P)=(\lambda +aP)N,$(4) ). All the parameters are the positive constants.

To investigate the properties of the corresponding ordinary differential model and the influence of diffusion in pattern formation as well as the influence of refuge effect and hunting cooperation, we denote $\mathbf{U}=(X,Y,Z{)}^T,\mathbf{D}=\mathrm{diag}(d_{11},d_{22},d_{33})$ and $\mathbf{J}(\mathbf{U})=\left(\begin{array}{c}{J}_1(\mathbf{U})\\ J_2(\mathbf{U})\\ J_3(\mathbf{U})\end{array}\right)=\left(\begin{array}{c}rX-a{X}^2-d_{1}X-{\alpha }_1(1-m_1)XY-{\displaystyle \frac{{\alpha }_2(1-m_2)XZ}{1+b(1-m_2)X}}\\ {\beta }_1(1-m_1)XY-({\alpha }_3+cZ)YZ-d_{2}Y\\ {\displaystyle \frac{{\beta }_2(1-m_2)XZ}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)YZ-d_{3}Z}\end{array}\right).$Then model (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }t}& =d_{11}\mathrm{\Delta }X+rX-a{X}^2-d_{1}X-{\alpha }_1(1-m_1)XY\\ & -\frac{{\alpha }_2(1-m_2)XZ}{1+b(1-m_2)X},x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }t}& =d_{22}\mathrm{\Delta }Y+{\beta }_1(1-m_1)XY-({\alpha }_3+cZ)YZ-d_{2}Y,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }t}& =d_{33}\mathrm{\Delta }Z+\frac{{\beta }_2(1-m_2)XZ}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)YZ-d_{3}Z,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }\nu }& =\frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }\nu }=0,x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ X(x,0)& =X_0(x)\ge 0,Y(x,0)=Y_0(x)\ge 0,Z(x,0)=Z_0(x)\ge 0,x\in \mathrm{\Omega },\end{array}\end{array}$(7) ) can be rewritten as (8) $\{\begin{array}{ll}{\mathbf{U}}_t-\mathbf{D}\mathrm{\Delta }\mathbf{U}=\mathbf{J}(\mathbf{U}),& x\in \mathrm{\Omega },t>0,\\ {\displaystyle \frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ \mathbf{U}(x,0)\ge 0,& x\in \mathrm{\Omega }.\end{array}$(8) The remaining of this paper is organized as follows. In Section 2, the existence and stability of equilibria, the existence, direction and stability of the Hopf bifurcation of model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) are discussed. In Section 3, the diffusion-driven Turing instability of model (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }t}& =d_{11}\mathrm{\Delta }X+rX-a{X}^2-d_{1}X-{\alpha }_1(1-m_1)XY\\ & -\frac{{\alpha }_2(1-m_2)XZ}{1+b(1-m_2)X},x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }t}& =d_{22}\mathrm{\Delta }Y+{\beta }_1(1-m_1)XY-({\alpha }_3+cZ)YZ-d_{2}Y,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }t}& =d_{33}\mathrm{\Delta }Z+\frac{{\beta }_2(1-m_2)XZ}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)YZ-d_{3}Z,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }\nu }& =\frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }\nu }=0,x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ X(x,0)& =X_0(x)\ge 0,Y(x,0)=Y_0(x)\ge 0,Z(x,0)=Z_0(x)\ge 0,x\in \mathrm{\Omega },\end{array}\end{array}$(7) ) is investigated, and the existence and non-existence of non-constant positive steady states are discussed by establishing positive upper and lower bounds and using Leray–Schauder degree theory. Numerical simulations are employed to confirm the theoretical results in Section 4. In Section 5, some discussions are given.

2. Stability and Hopf bifurcation of ODE model

In this section, we mainly consider the existence and stability of equilibria, the existence, direction and stability of the Hopf bifurcation of the ordinary differential equation (ODE) model corresponding to model (8(8) $\{\begin{array}{ll}{\mathbf{U}}_t-\mathbf{D}\mathrm{\Delta }\mathbf{U}=\mathbf{J}(\mathbf{U}),& x\in \mathrm{\Omega },t>0,\\ {\displaystyle \frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ \mathbf{U}(x,0)\ge 0,& x\in \mathrm{\Omega }.\end{array}$(8) ). The ODE model is as follows: (9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9)

2.1. The existence and local stability of equilibria

It is obvious that (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) have the following nonnegative constant steady states:

1. the trivial equilibrium $E_0(0,0,0)$;

2. the semi-trivial equilibrium in the absence of IG prey and IG predator $E_1(X_1,0,0)$ if condition $(H_1)$ holds, where $X_1=\frac{r-d_1}{a}$ and $(H_1):r>d_1$;

3. the semi-trivial equilibrium in the absence of IG predator $E_2(X_2,Y_2,0)$ if condition $(H_2)$ is true, where $X_2=\frac{d_2}{{\beta }_1(1-m_1)}$, $Y_2=\frac{{\beta }_1(r-d_1)(1-m_1)-a{d}_2}{{\alpha }_1{\beta }_1(1-m_1{)}^2}$ and $(H_2):{\beta }_1(r-d_1)(1-m_1)>a{d}_2$, i.e. $0<m_1<\frac{{\beta }_1(r-d_1)-a{d}_2}{{\beta }_1(r-d_1)}$;

4. the semi-trivial equilibrium in the absence of IG prey $E_3(X_3,0,Z_3)$ if assumption $(H_3)$ holds, where $X_3=\frac{d_3}{({\beta }_2-b{d}_3)(1-m_2)}$, $Z_3=\frac{{\beta }_2[(r-d_1)({\beta }_2-b{d}_3)(1-m_2)-a{d}_3]}{{\alpha }_2(1-m_2{)}^2({\beta }_2-b{d}_3{)}^2}$ and $(H_3):{\beta }_2>b{d}_3,(r-d_1)({\beta }_2-b{d}_3)(1-m_2)-a{d}_3>0$;

5. the co-exist equilibrium $E^{\ast }(X^{\ast },Y^{\ast },Z^{\ast })$ is expressed as follows: $X^{\ast }=\frac{c(Z^{\ast }{)}^2+{\alpha }_3{Z}^{\ast }+d_2}{{\beta }_1(1-m_1)},Y^{\ast }=\frac{(b{d}_3-{\beta }_2)(1-m_2)X^{\ast }+d_3}{{\beta }_3({\alpha }_3+c{Z}^{\ast })[1+b(1-m_2)X^{\ast }]}$

and $Z^{\ast }$ determined by the following equation: (10) $A_1{Z}^5+A_2{Z}^4+A_3{Z}^3+A_4{Z}^2+A_{5}Z+A_6=0,$(10) where $\begin{array}{rl}{A}_1& =ab{c}^3{\beta }_3(1-m_2),A_2=3ab{c}^2{\alpha }_3{\beta }_3(1-m_2),\\ A_3& =c{\beta }_3\{ab(1-m_2)(2c{d}_2+3{\alpha }_3^2)+c{\beta }_1(1-m_1)[a-b(1-m_2)(r-d_1)]\},\\ A_4& ={\beta }_3\{ab{\alpha }_3(1-m_2)(4c{d}_2+{\alpha }_3^2)+2c{\alpha }_3{\beta }_1(1-m_1)[a-b(1-m_2)(r-d_1)]\}\\ & -c{\beta }_1(1-m_1{)}^2(1-m_2)[{\alpha }_1({\beta }_2-b{d}_3)-{\alpha }_2{\beta }_1{\beta }_3],\\ A_5& =ab{d}_2{\alpha }_3^2{\beta }_3(1-m_2)+{\beta }_3(c{d}_2+{\alpha }_3^2)\{ab{d}_2(1-m_2)\\ & +{\beta }_1(1-m_1)[a-b(1-m_2)(r-d_1)]\}\\ & -{\beta }_1(1-m_1{)}^2\{c{\beta }_1{\beta }_3(r-d_1)+{\alpha }_3(1-m_2)[{\alpha }_1({\beta }_2-b{d}_3)-{\alpha }_2{\beta }_1{\beta }_3]\},\\ A_6& =d_2{\alpha }_3{\beta }_3\{ab{d}_2(1-m_2)+{\beta }_1(1-m_1)[a-b(1-m_2)(r-d_1)]\}\\ & -{\beta }_1(1-m_1{)}^2\{{\alpha }_3{\beta }_1{\beta }_3(r-d_1)+{\alpha }_1[d_2(1-m_2)({\beta }_2-b{d}_3)-d_3{\beta }_1(1-m_1)]\}.\end{array}$We can know that $A_1>0,A_2>0$ and $A_3>0$ can lead to $A_4>0$. From Descartes' Rule of Signs, it follows that the number of positive roots of Equation (10(10) $A_1{Z}^5+A_2{Z}^4+A_3{Z}^3+A_4{Z}^2+A_{5}Z+A_6=0,$(10) ) is shown in Table 1.

Table 1. The number of positive roots of Equation (10(10) $A_1{Z}^5+A_2{Z}^4+A_3{Z}^3+A_4{Z}^2+A_{5}Z+A_6=0,$(10) )

The sign of coefficients	Number of sign changes	Number of positive roots
$A_3>0,A_4>0,A_5>0,A_6>0$	0	0
$A_3>0,A_4>0,A_5>0,A_6<0$	1	1
$A_3>0,A_4>0,A_5<0,A_6>0$	2	2  or  0
$A_3>0,A_4>0,A_5<0,A_6<0$	1	1
$A_3<0,A_4>0,A_5>0,A_6>0$	2	2  or  0
$A_3<0,A_4>0,A_5>0,A_6<0$	3	3  or  1
$A_3<0,A_4>0,A_5<0,A_6>0$	4	4  or  2  or  0
$A_3<0,A_4>0,A_5<0,A_6<0$	3	3  or  1
$A_3<0,A_4<0,A_5>0,A_6>0$	2	2  or  0
$A_3<0,A_4<0,A_5>0,A_6<0$	3	3  or  1
$A_3<0,A_4<0,A_5<0,A_6>0$	2	2  or  0
$A_3<0,A_4<0,A_5<0,A_6<0$	1	1

Based on the above analysis, it can be seen that the existence of the co-exist equilibrium of model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) is very complicated. Therefore, we only discuss the second and fourth cases in this paper, that is, model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) has only one co-exist equilibrium if the assumption $(H_4):r>d_1,b{d}_3>{\beta }_2,A_3>0,A_4>0,A_6<0$ holds.

Next, we will use the standard linearization method to discuss the local stability of each equilibrium. The Jacobian matrix of the model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) at any non-negative equilibrium $(X,Y,Z)$ is given by (11) ${\mathbf{J}}_{\mathbf{U}}=\left(\begin{array}{ccc}{J}_{11}& J_{12}& J_{13}\\ J_{21}& J_{22}& J_{23}\\ J_{31}& J_{32}& J_{33}\end{array}\right),$(11) where $\begin{array}{rl}{J}_{11}& =r-2aX-d_1-{\alpha }_1(1-m_1)Y-\frac{{\alpha }_2(1-m_2)Z}{[1+b(1-m_2)X{]}^2},J_{12}=-{\alpha }_1(1-m_1)X,\\ J_{13}& =-\frac{{\alpha }_2(1-m_2)X}{1+b(1-m_2)X},J_{21}={\beta }_1(1-m_1)Y,\\ J_{22}& ={\beta }_1(1-m_1)X-({\alpha }_3+cZ)Z-d_2,J_{23}=-{\alpha }_{3}Y-2cYZ,\\ J_{31}& =\frac{{\beta }_2(1-m_2)Z}{[1+b(1-m_2)X{]}^2},J_{32}={\beta }_3({\alpha }_3+cZ)Z,\\ J_{33}& =\frac{{\beta }_2(1-m_2)X}{1+b(1-m_2)X}+{\beta }_3({\alpha }_{3}Y+2cYZ)-d_3.\end{array}$

Lemma 2.1

The trivial equilibrium $E_0(0,0,0)$ of model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) is always unstable.

Proof.

The characteristic equation of model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) at the trivial equilibrium $E_0$ is $[\lambda -(r-d_1)](\lambda +d_2)(\lambda +d_3)=0.$Thus we can get that the corresponding eigenvalues are ${\lambda }_1=r-d_1>0,{\lambda }_2=-d_2<0,{\lambda }_3=-d_3<0.$So the trivial equilibrium $E_0$ is unstable.

Lemma 2.2

If the condition $(H_1)$ and $a>max\{\frac{{\beta }_1(1-m_1)}{d_2},\frac{({\beta }_2-b{d}_3)(1-m_2)}{d_3}\}$ hold, then the semi-trivial equilibrium $E_1(X_1,0,0)$ is locally asymptotically stable.

Proof.

Through the Jacobian matrix (11(11) ${\mathbf{J}}_{\mathbf{U}}=\left(\begin{array}{ccc}{J}_{11}& J_{12}& J_{13}\\ J_{21}& J_{22}& J_{23}\\ J_{31}& J_{32}& J_{33}\end{array}\right),$(11) ), the characteristic equation of model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) at the equilibrium $E_1$ is $\begin{array}{rl}& [\lambda -(d_1-r)][\lambda -\frac{{\beta }_1(1-m_1)(r-d_1)-a{d}_2}{a}]\\ & [\lambda -\frac{({\beta }_2-b{d}_3)(1-m_2)(r-d_1)-a{d}_3}{a+b(1-m_2)(r-d_1)}]=0.\end{array}$Thus one eigenvalue is ${\lambda }_1=d_1-r<0$ and other eigenvalues ${\lambda }_2=\frac{{\beta }_1(1-m_1)(r-d_1)-a{d}_2}{a}<0$, ${\lambda }_3=-\frac{({\beta }_2-b{d}_3)(1-m_2)(r-d_1)-a{d}_3}{a+b(1-m_2)(r-d_1)}<0$ if and only if $a>max\{\frac{{\beta }_1(1-m_1)}{d_2},\frac{({\beta }_2-b{d}_3)(1-m_2)}{d_3}\}$. Thus the semi-trivial equilibrium $E_1$ is locally asymptotically stable.

Lemma 2.3

Under the case of conditions $(H_2)$ and $(H_5)$, the semi-trivial equilibrium $E_2(X_2,Y_2,0)$ is locally asymptotically stable, where $(H_5):a>\frac{{\alpha }_1{\beta }_1(-m_1{)}^2[d_2({\beta }_2-b{d}_3)(1-m_2)-d_3{\beta }_1(1-m_1)]}{d_2{\alpha }_3{\beta }_3[b{d}_2(1-m_2)+{\beta }_1(1-m_1)]}+\frac{{\alpha }_3{\beta }_1(1-m_1)(r-d_1)}{d_2{\alpha }_3}$.

Proof.

From matrix (11(11) ${\mathbf{J}}_{\mathbf{U}}=\left(\begin{array}{ccc}{J}_{11}& J_{12}& J_{13}\\ J_{21}& J_{22}& J_{23}\\ J_{31}& J_{32}& J_{33}\end{array}\right),$(11) ), the characteristic equation of model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) at equilibrium $E_2$ takes the form $(\lambda -B_{33})[{\lambda }^2-(B_{11}+B_{22})\lambda +B_{11}{B}_{22}-B_{12}{B}_{21}]=0,$where $\begin{array}{rl}{B}_{11}& =r-d_1-2a{X}_2-{\alpha }_1(1-m_1)Y_2,B_{12}=-{\alpha }_1(1-m_1)X_2,\\ B_{21}& ={\beta }_1(1-m_1)Y_2,\\ B_{22}& ={\beta }_1(1-m_1)X_2-d_2,B_{33}=\frac{{\beta }_2(1-m_2)X_2}{1+b(1-m_2)X_2}+{\alpha }_3{\beta }_3{Y}_2-d_3.\end{array}$Therefore, ${\lambda }_1=B_{33}$, ${\lambda }_2$ and ${\lambda }_3$ are the roots of the following equation: ${\lambda }^2+G_1\lambda +G_2=0,$where $G_1=-(B_{11}+B_{22})>0,G_2=B_{11}{B}_{22}-B_{12}{B}_{21}>0$. It has that ${\lambda }_2<0$ and ${\lambda }_3<0$. Hence, the semi-trivial equilibrium $E_2$ is locally asymptotically stable if $(H_5)$ is true.

Lemma 2.4

If conditions $(H_3)$ and $(H_6)$ hold, then the semi-trivial equilibrium $E_3(X_3,0,Z_3)$ is locally asymptotically stable, where $(H_6):a\ge \frac{{\overline{d}}_3[{\beta }_2(1-m_2)-1]}{d_3[2{\beta }_2(1-m_2)-1]}$ and ${\beta }_1<\frac{c{\beta }_2({\overline{d}}_3-a{d}_3)({\overline{d}}_1+{\overline{d}}_2)}{d_3{\alpha }_2^2(1-m_1)(1-m_2{)}^3({\beta }_2-b{d}_3{)}^3}$, here ${\overline{d}}_1=(1-m_2)({\beta }_2-b{d}_3)[{\alpha }_2{\alpha }_3(1-m_2)({\beta }_2-b{d}_3)+c{\beta }_2(r-d_1)]-ac{\alpha }_3{\beta }_2$, ${\overline{d}}_2=d_2{\alpha }_2^2(1-m_2{)}^4({\beta }_2-b{d}_3{)}^4$,${\overline{d}}_3=(r-d_1)(1-m_2)({\beta }_2-b{d}_3)$.

Proof.

From the Jacobian matrix (11(11) ${\mathbf{J}}_{\mathbf{U}}=\left(\begin{array}{ccc}{J}_{11}& J_{12}& J_{13}\\ J_{21}& J_{22}& J_{23}\\ J_{31}& J_{32}& J_{33}\end{array}\right),$(11) ), the characteristic equation of model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) at the semi-trivial equilibrium $E_3$ takes $(\lambda -C_{22})[{\lambda }^2-(C_{11}+C_{33})\lambda +C_{11}{C}_{33}-C_{13}{C}_{31}]=0,$where $\begin{array}{rl}{C}_{13}& =-\frac{{\alpha }_2(1-m_2)X_3}{1+b(1-m_2)X_3},C_{31}=\frac{{\beta }_2(1-m_2)Z_3}{[1+b(1-m_2)X_3{]}^2},\\ C_{33}& =\frac{{\beta }_2(1-m_2)X_3}{1+b(1-m_2)X_3}-d_3,\\ C_{11}& =r-d_1-2a{X}_3-\frac{{\alpha }_2(1-m_2)Z_3}{[1+b(1-m_2)X_3{]}^2},\\ C_{22}& ={\beta }_1(1-m_1)X_3-({\alpha }_3+c{Z}_3)Z_3-d_2.\end{array}$Further, ${\lambda }_1=C_{22}$, ${\lambda }_2$ and ${\lambda }_3$ are the roots of following equation: ${\lambda }^2+G_3\lambda +G_4=0,$where $G_3=-(C_{11}+C_{33}),G_4=C_{11}{C}_{33}-C_{13}{C}_{31}>0$. It follows that ${\lambda }_1<0,{\lambda }_2<0,{\lambda }_3<0$ if and only if $C_{22}<0,G_3\ge 0$. That is, when the assumption $(H_6)$ holds, the semi-trivial equilibrium $E_3$ is locally asymptotically stable.

Next, we analyse the local stability of the co-existence equilibrium $E^{\ast }(X^{\ast },Y^{\ast },Z^{\ast })$. According to the matrix (11(11) ${\mathbf{J}}_{\mathbf{U}}=\left(\begin{array}{ccc}{J}_{11}& J_{12}& J_{13}\\ J_{21}& J_{22}& J_{23}\\ J_{31}& J_{32}& J_{33}\end{array}\right),$(11) ), the characteristic equation of model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) at $E^{\ast }$ is (12) ${\lambda }^3+G_5{\lambda }^2+G_6\lambda +G_7=0,$(12) where $\begin{array}{rl}{G}_5& =-(D_{11}+D_{22}+D_{33}),G_6=D_{11}{D}_{33}-D_{13}{D}_{31}-D_{12}{D}_{21}-D_{23}{D}_{32},\\ G_7& =D_{11}{D}_{23}{D}_{32}+D_{12}{D}_{21}{D}_{33}-D_{21}{D}_{13}{D}_{32}-D_{12}{D}_{23}{D}_{31},\\ D_{11}& =\frac{b{\alpha }_2(1-m_2)Z^{\ast }}{[1+b(1-m_2)X^{\ast }{]}^2}-a{X}^{\ast },D_{12}=-{\alpha }_1(1-m_1)X^{\ast },D_{33}=c{\beta }_3{Y}^{\ast }{Z}^{\ast },\\ D_{21}& ={\beta }_1(1-m_1)Y^{\ast },D_{22}=0,D_{23}=-{\alpha }_3{Y}^{\ast }-2c{Y}^{\ast }{Z}^{\ast },\\ D_{31}& =\frac{{\beta }_2(1-m_2)Z^{\ast }}{[1+b(1-m_2)X^{\ast }{]}^2},D_{32}={\beta }_3({\alpha }_3+c{Z}^{\ast })Z^{\ast },D_{13}=-\frac{{\alpha }_2(1-m_2)X^{\ast }}{1+b(1-m_2)X^{\ast }}.\end{array}$Assuming that $(H_7):D_{11}<-D_{33},D_{11}{D}_{33}-D_{13}{D}_{31}<D_{12}{D}_{21}+D_{23}{D}_{32},D_{11}{D}_{23}{D}_{32}-D_{21}{D}_{13}{D}_{32}>D_{12}{D}_{23}{D}_{31}-D_{12}{D}_{21}{D}_{33}$ holds, we can conclude that $G_5>0,G_6>0$ and $G_7>0$.

Lemma 2.5

All roots of Equation (12(12) ${\lambda }^3+G_5{\lambda }^2+G_6\lambda +G_7=0,$(12) ) have negative real parts if and only if the conditions $(H_4),(H_7)$ and $G_5{G}_6-G_7>0$ are satisfied. Equation (12(12) ${\lambda }^3+G_5{\lambda }^2+G_6\lambda +G_7=0,$(12) ) has one negative root and a pair of complex roots with a positive real part if and only if the conditions $(H_4),(H_7)$ and $G_5{G}_6-G_7<0$ are true.

Based on Lemma 2.5, we obtain the following theorem.

Theorem 2.6

Assume that conditions $(H_4)$ and $(H_7)$ are satisfied, the co-existence equilibrium $E^{\ast }(X^{\ast },Y^{\ast },Z^{\ast })$ is locally asymptotically stable if $G_5{G}_6-G_7>0$, but is unstable if $G_5{G}_6-G_7<0$. Further, the co-existence equilibrium $E^{\ast }$ loses stability when $G_5{G}_6-G_7$ passes through 0. That is, the Hopf bifurcation occurs at the co-existence equilibrium $E^{\ast }$ when $G_5{G}_6-G_7=0$.

Next, we investigate the existence of the Hopf bifurcation by considering the prey refuge $m_1$ as a bifurcation parameter. After some computations, $G_5{G}_6-G_7=0$ reduces to (13) $G(m_1)=G_8{m}_1^2+G_9{m}_1+G_{10}=0,$(13) where $\begin{array}{rl}{G}_8& ={\alpha }_1{\beta }_1(\lambda -D_{33})X^{\ast }{Y}^{\ast },\\ G_9& ={\beta }_1{Y}^{\ast }{D}_{13}{D}_{32}-{\alpha }_1{X}^{\ast }{D}_{23}{D}_{31}-2{\alpha }_1{\beta }_1(\lambda -D_{33})X^{\ast }{Y}^{\ast },\\ G_{10}& ={\lambda }^3-(D_{11}+D_{33}){\lambda }^2+(D_{11}{D}_{33}-D_{13}{D}_{31}-D_{23}{D}_{32})\lambda +D_{11}{D}_{23}{D}_{32}\\ & +{\alpha }_1{\beta }_1(\lambda -D_{33})X^{\ast }{Y}^{\ast }-{\beta }_1{Y}^{\ast }{D}_{13}{D}_{32}+{\alpha }_1{X}^{\ast }{D}_{23}{D}_{31}.\end{array}$Without loss of generality, let $m_1^{\ast }$ be the positive root of Equation (13(13) $G(m_1)=G_8{m}_1^2+G_9{m}_1+G_{10}=0,$(13) ). Then Equation (12(12) ${\lambda }^3+G_5{\lambda }^2+G_6\lambda +G_7=0,$(12) ) has a negative root denoted ${\lambda }_1$ and a pair of complex conjugates roots which are denoted ${\lambda }_{2,3}=u(m_1)\pm iv(m_1)$ such that $u(m_1^{\ast })=0,v(m_1^{\ast })\ne 0$.

Substituting ${\lambda }_2=u(m_1)+iv(m_1)$ into Equation (12(12) ${\lambda }^3+G_5{\lambda }^2+G_6\lambda +G_7=0,$(12) ) and separating the real and imaginary parts gives (14) $\begin{array}{rl}& u^3-3u{v}^2+G_5(u^2-v^2)+G_{6}u+G_7=0,\end{array}$(14) (15) $\begin{array}{rl}& 3u^{2}v-v^3+2G_{5}uv+G_{6}v=0.\end{array}$(15) Because $u(m_1)\ne 0$, from Equation (15(15) $\begin{array}{rl}& 3u^{2}v-v^3+2G_{5}uv+G_{6}v=0.\end{array}$(15) ) we have (16) $v^2=3u^2+2G_{5}u+G_6.$(16) Then we substitute (16(16) $v^2=3u^2+2G_{5}u+G_6.$(16) ) in Equation (14(14) $\begin{array}{rl}& u^3-3u{v}^2+G_5(u^2-v^2)+G_{6}u+G_7=0,\end{array}$(14) ) and get (17) $8u^3+8u^2{G}_5+2(G_5^2+G_6)u+G_5{G}_6-G_7=0.$(17) So Equation (17(17) $8u^3+8u^2{G}_5+2(G_5^2+G_6)u+G_5{G}_6-G_7=0.$(17) ) has a root $u(m_1^{\ast })$ if and only if $G_5(u(m_1^{\ast }))G_6(u(m_1^{\ast }))-G_7(u(m_1^{\ast }))=0$. Furthermore, since $8u^2+8u{G}_5+2(G_5^2+G_6)=0$ has no real root, u = 0 is the only root. Thereby, when $m_1=m_1^{\ast }$, the roots of (12(12) ${\lambda }^3+G_5{\lambda }^2+G_6\lambda +G_7=0,$(12) ) are ${\lambda }_1=-G_5,{\lambda }_{2,3}=\pm i\sqrt{G_6}$.

We will prove that $\frac{\mathrm{d}u}{\mathrm{d}{m}_1}{|}_{m_1=m_1^{\ast }}\ne 0$. Differentiating Equation (17(17) $8u^3+8u^2{G}_5+2(G_5^2+G_6)u+G_5{G}_6-G_7=0.$(17) ) with respect to $m_1$ gives (18) $\begin{array}{rl}& 24u^2\frac{\mathrm{d}u}{\mathrm{d}{m}_1}+16u{G}_5\frac{\mathrm{d}u}{\mathrm{d}{m}_1}+8u^2\frac{\mathrm{d}{G}_5}{\mathrm{d}{m}_1}+2\frac{\mathrm{d}u}{\mathrm{d}{m}_1}(G_5^2+G_6)\\ & +2(2G_5\frac{\mathrm{d}{G}_5}{\mathrm{d}{m}_1}+\frac{\mathrm{d}{G}_6}{\mathrm{d}{m}_1})u+\frac{\mathrm{d}(G_5{G}_6-G_7)}{\mathrm{d}{m}_1}=0.\end{array}$(18) When $m_1=m_1^{\ast }$ in (18(18) $\begin{array}{rl}& 24u^2\frac{\mathrm{d}u}{\mathrm{d}{m}_1}+16u{G}_5\frac{\mathrm{d}u}{\mathrm{d}{m}_1}+8u^2\frac{\mathrm{d}{G}_5}{\mathrm{d}{m}_1}+2\frac{\mathrm{d}u}{\mathrm{d}{m}_1}(G_5^2+G_6)\\ & +2(2G_5\frac{\mathrm{d}{G}_5}{\mathrm{d}{m}_1}+\frac{\mathrm{d}{G}_6}{\mathrm{d}{m}_1})u+\frac{\mathrm{d}(G_5{G}_6-G_7)}{\mathrm{d}{m}_1}=0.\end{array}$(18) ), we have that from $u(m_1^{\ast })=0$ $\frac{\mathrm{d}u}{\mathrm{d}{m}_1}{|}_{m_1=m_1^{\ast }}=-\frac{1}{2(G_5^2+G_6)}\frac{\mathrm{d}(G_5{G}_6-G_7)}{\mathrm{d}{m}_1}{|}_{m_1=m_1^{\ast }}\ne 0.$It shows that the hypotheses of Hopf theorem are satisfied. Thus we obtain the following result.

Theorem 2.7

Suppose that the coefficients of Equation (12(12) ${\lambda }^3+G_5{\lambda }^2+G_6\lambda +G_7=0,$(12) ) at $E^{\ast }$ are functions of $m_1$ and $G_i(m_1)(i=5,6,7)$ at $m_1=m_1^{\ast }$ satisfy the following conditions: $G_i(m_1^{\ast })>0,i=5,6,7,G_5{G}_6-G_7{|}_{m_1=m_1^{\ast }}=0,\frac{\mathrm{d}(G_5{G}_6-G_7)}{\mathrm{d}{m}_1}{|}_{m_1=m_1^{\ast }}\ne 0.$Thus Hopf bifurcation occurs at the co-existence equilibrium $E^{\ast }$ when $m_1=m_1^{\ast }$.

2.2. Direction of Hopf bifurcation and stability of the bifurcating periodic solutions

In this part, we will derive an explicit algorithm to determine the direction of Hopf bifurcation and stability of the bifurcating periodic solutions for model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) by using the centre manifold theory and normal form method.

For the sake of convenience, we translate the equilibrium ${\mathbf{U}}^{\mathbf{\ast }}=(X^{\ast },Y^{\ast },Z^{\ast })$ to the origin by the transformation $\overline{X}=X-X^{\ast },\overline{Y}=Y-Y^{\ast },\overline{Z}=Z-Z^{\ast }$. We denote $\overline{X},\overline{Y},\overline{Z}$ by X, Y, Z to calculate conveniently, so the model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) becomes (19) $\left(\begin{array}{c}{\displaystyle \frac{\mathrm{d}X}{\mathrm{d}t}}\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}\end{array}\right)={\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\mathbf{\ast }})\left(\begin{array}{c}X\\ Y\\ Z\end{array}\right)+\left(\begin{array}{c}{N}_1(\mathbf{U})\\ N_2(\mathbf{U})\\ N_3(\mathbf{U})\end{array}\right),$(19) where $\begin{array}{rl}{N}_1(\mathbf{U})& =(\frac{2b{\alpha }_2(1-m_2{)}^2{Z}^{\ast }}{[1+b(1-m_2)X^{\ast }{]}^3}-2a)X^2-\frac{{\alpha }_2(1-m_2)}{[1+b(1-m_2)X^{\ast }{]}^2}XZ\\ & -{\alpha }_1(1-m_1)XY-\frac{6b^2{\alpha }_2(1-m_2{)}^3{Z}^{\ast }}{[1+b(1-m_2)X^{\ast }{]}^4}{X}^3+\frac{b{\alpha }_2(1-m_2{)}^2}{[1+b(1-m_2)X^{\ast }{]}^3}{X}^{2}Z,\\ N_2(\mathbf{U})& =-2c{Y}^{\ast }{Z}^2+{\beta }_1(1-m_1)XY-({\alpha }_3+2c{Z}^{\ast })YZ,\\ N_3(\mathbf{U})& =\frac{6b^2{\beta }_2(1-m_2{)}^3{Z}^{\ast }}{[1+b(1-m_2)X^{\ast }{]}^4}{X}^3-\frac{b{\beta }_2(1-m_2{)}^2{X}^2(Z+2Z^{\ast })}{[1+b(1-m_2)X^{\ast }{]}^3}\\ & +2c{\beta }_3{Y}^{\ast }{Z}^2+\frac{{\beta }_2(1-m_2)}{[1+b(1-m_2)X^{\ast }{]}^2}XZ+({\alpha }_3{\beta }_3+2c{\beta }_3{Z}^{\ast })YZ.\end{array}$Let $p_1$ and $p_2$ be the eigenvectors corresponding to eigenvalues ${\lambda }_1=-G_5$ and ${\lambda }_{2,3}=\pm i\sqrt{G_6}\triangleq \pm i\omega$ of ${\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\mathbf{\ast }})$ at $m_1=m_1^{\ast }$, respectively. We can obtain $p_1$ and $p_2$, $p_1=\left(\begin{array}{c}{p}_{11}\\ p_{12}\\ p_{13}\end{array}\right),p_2=\left(\begin{array}{c}{p}_{21}\\ p_{22}\\ p_{23}\end{array}\right),$where $\begin{array}{rl}{p}_{11}& =\frac{-[(J_{11}{J}_{13}+J_{12}{J}_{23}){\omega }^2+J_{12}^2{J}_{21}{J}_{23}]+[J_{12}(J_{11}{J}_{23}-J_{13}{J}_{21})-J_{13}{\omega }^2]i\omega }{(J_{12}{J}_{21}+{\omega }^2{)}^2+J_{11}^2{\omega }^2},\\ p_{12}& =\frac{(J_{11}{J}_{23}-J_{13}{J}_{21})J_{12}{J}_{21}-J_{13}{J}_{21}{\omega }^2-p_{12}^{I}i\omega }{(J_{11}{J}_{21}+{\omega }^2{)}^2+J_{11}^2{\omega }^2},\\ p_{12}^I& =J_{23}(J_{12}{J}_{21}+{\omega }^2)+J_{11}(J_{11}{J}_{23}-J_{13}{J}_{21}),p_{13}=1,\\ p_{21}& =\frac{[(J_{11}{J}_{13}+J_{12}{J}_{23})G_5^2-J_{12}^2{J}_{21}{J}_{23}]-[J_{12}(J_{11}{J}_{23}-J_{13}{J}_{21})+J_{13}{G}_5^2]G_5}{(J_{12}{J}_{21}-G_5^2{)}^2-J_{11}^2{G}_5^2},\\ p_{22}& =\frac{(J_{11}{J}_{23}-J_{13}{J}_{21})J_{12}{J}_{21}+J_{13}{J}_{21}{G}_5^2+p_{22}^I{G}_5}{(J_{11}{J}_{21}-G_5^2{)}^2-J_{11}^2{G}_5^2},\\ p_{22}^I& =J_{23}(J_{12}{J}_{21}-G_5^2)+J_{11}(J_{11}{J}_{23}-J_{13}{J}_{21}),p_{23}=1.\end{array}$Denote the matrix $E=(\mathrm{Re}(p_1),-\mathrm{Im}(p_1),p_2)=\left(\begin{array}{ccc}{e}_{11}& e_{12}& e_{13}\\ e_{21}& e_{22}& e_{23}\\ e_{31}& e_{32}& e_{33}\end{array}\right),$with $\begin{array}{rl}{e}_{11}& =-\frac{(J_{11}{J}_{13}+J_{12}{J}_{23}){\omega }^2+J_{12}^2{J}_{21}{J}_{23}}{(J_{12}{J}_{21}+{\omega }^2{)}^2+J_{11}^2{\omega }^2},e_{12}=\frac{J_{13}{\omega }^2-J_{12}(J_{11}{J}_{23}-J_{13}{J}_{21})}{(J_{12}{J}_{21}+{\omega }^2{)}^2+J_{11}^2{\omega }^2},\\ e_{13}& =p_{21},e_{21}=\frac{(J_{11}{J}_{23}-J_{13}{J}_{21})J_{12}{J}_{21}-J_{13}{J}_{21}{\omega }^2}{(J_{11}{J}_{21}+{\omega }^2{)}^2+J_{11}^2{\omega }^2},e_{23}=p_{22},e_{31}=1,\\ e_{22}& =\frac{J_{23}(J_{12}{J}_{21}+{\omega }^2)+J_{11}(J_{11}{J}_{23}-J_{13}{J}_{21})}{(J_{11}{J}_{21}+{\omega }^2{)}^2+J_{11}^2{\omega }^2},e_{32}=0,e_{33}=1.\end{array}$Then we obtain $E^{-1}=\frac{1}{F}\left(\begin{array}{ccc}{f}_{11}& f_{12}& f_{13}\\ f_{21}& f_{22}& f_{23}\\ f_{31}& f_{32}& f_{33}\end{array}\right),$where $\begin{array}{rl}F& =e_{11}(e_{22}{e}_{33}-e_{23}{e}_{32})+e_{12}(e_{23}{e}_{31}-e_{21}{e}_{33})+e_{13}(e_{21}{e}_{32}-e_{22}{e}_{31}),\\ f_{11}& =e_{22}{e}_{33}-e_{23}{e}_{32},f_{12}=e_{13}{e}_{32}-e_{12}{e}_{33},f_{13}=e_{12}{e}_{23}-e_{13}{e}_{22},\\ f_{21}& =e_{23}{e}_{31}-e_{21}{e}_{33},f_{22}=e_{11}{e}_{33}-e_{13}{e}_{31},f_{23}=e_{13}{e}_{21}-e_{11}{e}_{23},\\ f_{31}& =e_{21}{e}_{32}-e_{22}{e}_{31},f_{32}=e_{12}{e}_{31}-e_{11}{e}_{32},f_{33}=e_{11}{e}_{22}-e_{12}{e}_{21}.\end{array}$By the transformation $(X,Y,Z{)}^T=E(x,y,z{)}^T$, system (19(19) $\left(\begin{array}{c}{\displaystyle \frac{\mathrm{d}X}{\mathrm{d}t}}\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}\end{array}\right)={\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\mathbf{\ast }})\left(\begin{array}{c}X\\ Y\\ Z\end{array}\right)+\left(\begin{array}{c}{N}_1(\mathbf{U})\\ N_2(\mathbf{U})\\ N_3(\mathbf{U})\end{array}\right),$(19) ) becomes (20) $\left(\begin{array}{c}{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}\\ {\displaystyle \frac{\mathrm{d}y}{\mathrm{d}t}}\\ {\displaystyle \frac{\mathrm{d}z}{\mathrm{d}t}}\end{array}\right)=E^{-1}\left(\begin{array}{c}{\displaystyle \frac{\mathrm{d}X}{\mathrm{d}t}}\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}\end{array}\right)=E^{-1}{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\mathbf{\ast }})E\left(\begin{array}{c}x\\ y\\ z\end{array}\right)+\left(\begin{array}{c}{n}_1(x,y,z)\\ n_2(x,y,z)\\ n_3(x,y,z)\end{array}\right),$(20) where $\begin{array}{rl}& \left(\begin{array}{c}{n}_1(x,y,z)\\ n_2(x,y,z)\\ n_3(x,y,z)\end{array}\right)\\ & =E^{-1}\left(\begin{array}{c}{N}_1(e_{11}X+e_{12}Y+e_{13}Z,e_{21}X+e_{22}Y+e_{23}Z,e_{31}X+e_{32}Y+e_{33}Z)\\ N_2(e_{11}X+e_{12}Y+e_{13}Z,e_{21}X+e_{22}Y+e_{23}Z,e_{31}X+e_{32}Y+e_{33}Z)\\ N_3(e_{11}X+e_{12}Y+e_{13}Z,e_{21}X+e_{22}Y+e_{23}Z,e_{31}X+e_{32}Y+e_{33}Z)\end{array}\right),\end{array}$and $E^{-1}{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\mathbf{\ast }})E=\left(\begin{array}{ccc}0& -\omega & 0\\ \omega & 0& 0\\ 0& 0& -G_5\end{array}\right).$Thus the system (20(20) $\left(\begin{array}{c}{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}\\ {\displaystyle \frac{\mathrm{d}y}{\mathrm{d}t}}\\ {\displaystyle \frac{\mathrm{d}z}{\mathrm{d}t}}\end{array}\right)=E^{-1}\left(\begin{array}{c}{\displaystyle \frac{\mathrm{d}X}{\mathrm{d}t}}\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}\end{array}\right)=E^{-1}{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\mathbf{\ast }})E\left(\begin{array}{c}x\\ y\\ z\end{array}\right)+\left(\begin{array}{c}{n}_1(x,y,z)\\ n_2(x,y,z)\\ n_3(x,y,z)\end{array}\right),$(20) ) can be written as (21) $\left(\begin{array}{c}{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}\\ {\displaystyle \frac{\mathrm{d}y}{\mathrm{d}t}}\\ {\displaystyle \frac{\mathrm{d}z}{\mathrm{d}t}}\end{array}\right)=\left(\begin{array}{ccc}0& -\omega & 0\\ \omega & 0& 0\\ 0& 0& -G_5\end{array}\right)\left(\begin{array}{c}x\\ y\\ z\end{array}\right)+\left(\begin{array}{c}{n}_1(x,y,z)\\ n_2(x,y,z)\\ n_3(x,y,z)\end{array}\right).$(21) By the algorithm proposed by Hassard et al. [17], we get $\begin{array}{rl}{g}_{11}& =\frac{1}{4}[\frac{{\mathrm{\partial }}^2{n}_1}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2{n}_1}{\mathrm{\partial }{y}^2}+i(\frac{{\mathrm{\partial }}^2{n}_2}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2{n}_2}{\mathrm{\partial }{y}^2})],\\ g_{02}& =\frac{1}{4}[\frac{{\mathrm{\partial }}^2{n}_1}{\mathrm{\partial }{x}^2}-\frac{{\mathrm{\partial }}^2{n}_1}{\mathrm{\partial }{y}^2}-2\frac{{\mathrm{\partial }}^2{n}_2}{\mathrm{\partial }x\mathrm{\partial }y}+i(\frac{{\mathrm{\partial }}^2{n}_2}{\mathrm{\partial }{x}^2}-\frac{{\mathrm{\partial }}^2{n}_2}{\mathrm{\partial }{y}^2}+2\frac{{\mathrm{\partial }}^2{n}_1}{\mathrm{\partial }x\mathrm{\partial }y})],\\ g_{20}& =\frac{1}{4}[\frac{{\mathrm{\partial }}^2{n}_1}{\mathrm{\partial }{x}^2}-\frac{{\mathrm{\partial }}^2{n}_1}{\mathrm{\partial }{y}^2}+2\frac{{\mathrm{\partial }}^2{n}_2}{\mathrm{\partial }x\mathrm{\partial }y}+i(\frac{{\mathrm{\partial }}^2{n}_2}{\mathrm{\partial }{x}^2}-\frac{{\mathrm{\partial }}^2{n}_2}{\mathrm{\partial }{y}^2}-2\frac{{\mathrm{\partial }}^2{n}_1}{\mathrm{\partial }x\mathrm{\partial }y})],\end{array}$and $\begin{array}{rl}{G}_{21}& =\frac{1}{8}[\frac{{\mathrm{\partial }}^3{n}_1}{\mathrm{\partial }{x}^3}+\frac{{\mathrm{\partial }}^3{n}_1}{\mathrm{\partial }x\mathrm{\partial }{y}^2}+\frac{{\mathrm{\partial }}^3{n}_2}{\mathrm{\partial }{x}^2\mathrm{\partial }y}+\frac{{\mathrm{\partial }}^3{n}_2}{\mathrm{\partial }{y}^3}\\ & +i(\frac{{\mathrm{\partial }}^3{n}_2}{\mathrm{\partial }{x}^3}+\frac{{\mathrm{\partial }}^3{n}_2}{\mathrm{\partial }x\mathrm{\partial }{y}^2}+\frac{{\mathrm{\partial }}^3{n}_2}{\mathrm{\partial }x\mathrm{\partial }{y}^2}-\frac{{\mathrm{\partial }}^3{n}_1}{\mathrm{\partial }{x}^2\mathrm{\partial }y}-\frac{{\mathrm{\partial }}^3{n}_1}{\mathrm{\partial }{y}^3})],\\ n_{11}& =\frac{1}{4}(\frac{{\mathrm{\partial }}^2{n}_3}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2{n}_3}{\mathrm{\partial }{y}^2}),n_{20}=\frac{1}{4}(\frac{{\mathrm{\partial }}^2{n}_3}{\mathrm{\partial }{x}^2}-\frac{{\mathrm{\partial }}^2{n}_3}{\mathrm{\partial }{y}^2}-2i\frac{{\mathrm{\partial }}^2{n}_3}{\mathrm{\partial }x\mathrm{\partial }y}).\end{array}$By solving the equations $\mu w_{11}=-n_{11}$ and $(\mu -2i\omega )w_{20}=-n_{20}$, we have $w_{11}=-\frac{n_{11}}{\mu },w_{20}=-\frac{n_{20}}{\mu -2i\omega }$. Moreover, we have $\begin{array}{rl}{G}_{110}& =\frac{1}{2}[\frac{{\mathrm{\partial }}^2{n}_1}{\mathrm{\partial }x\mathrm{\partial }z}+\frac{{\mathrm{\partial }}^2{n}_2}{\mathrm{\partial }y\mathrm{\partial }z}+i(\frac{{\mathrm{\partial }}^2{n}_2}{\mathrm{\partial }x\mathrm{\partial }z}-\frac{{\mathrm{\partial }}^2{n}_1}{\mathrm{\partial }y\mathrm{\partial }z})],\\ G_{101}& =\frac{1}{2}[\frac{{\mathrm{\partial }}^2{n}_1}{\mathrm{\partial }x\mathrm{\partial }z}-\frac{{\mathrm{\partial }}^2{n}_2}{\mathrm{\partial }y\mathrm{\partial }z}+i(\frac{{\mathrm{\partial }}^2{n}_1}{\mathrm{\partial }y\mathrm{\partial }z}+\frac{{\mathrm{\partial }}^2{n}_2}{\mathrm{\partial }x\mathrm{\partial }z})],\\ g_{21}& =G_{21}+2G_{110}{w}_{11}+G_{101}{w}_{20}.\end{array}$Then, we can compute the following values: (22) $\begin{array}{rl}{c}_1(0)& =\frac{i}{2\omega }(g_{20}{g}_{11}-2|g_{11}{|}^2-\frac{1}{3}|g_{02}{|}^2)+\frac{1}{2}{g}_{21},\\ {\nu }_1& =-\frac{\mathrm{Re}\{c_1(0)\}}{\mathrm{Re}\{{\lambda }^{\mathrm{\prime }}(m_1^{\ast })\}},\\ {\varrho }_1& =-\frac{\mathrm{Im}\{c_1(0)\}+{\nu }_1\mathrm{Im}\{{\lambda }^{\mathrm{\prime }}(m_1^{\ast })\}}{\omega },\\ {\varsigma }_1& =2\mathrm{Re}\{c_1(0)\},\end{array}$(22) which determine the direction of Hopf bifurcation and the stability of bifurcating periodic solutions. Based on the above analysis and the conclusions of Hassard et al. [17], we have the following results.

Theorem 2.8

Suppose that assumption $(H_4)$ satisfies, Hopf bifurcation occurs at the co-existence equilibrium $E^{\ast }$ of model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) when $m_1=m_1^{\ast }$. If ${\nu }_1>0(<0)$, then the direction of Hopf bifurcation is supercritical (subcritical); if ${\varrho }_1>0(<0)$, then the period of bifurcating periodic solutions increases (decreases); if ${\varsigma }_1<0(>0)$, then the bifurcating periodic solutions are stable (unstable).

3. Dynamics of the reaction–diffusion model

3.1. Linear stability and Turing instability

In this section, we consider the influence of the diffusion on the stability of the co-existence equilibrium $E^{\ast }$ and investigate the diffusion-induced spatially inhomogeneous steady state. In [46], Turing concluded that the reaction–diffusion model may exhibit spatial patterns under the fact that the equilibrium is linearly stable in the absence of diffusion, but becomes linearly unstable in the presence of diffusion. Such instability is called a Turing instability or diffusion-driven instability.

Let $0={\mu }_1<{\mu }_2<{\mu }_3<\cdots$ be the eigenvalues of the operator $-\mathrm{\Delta }$ on Ω with the homogeneous Neumann boundary condition and $E({\mu }_i)$ be the eigenspace corresponding to ${\mu }_i$ in $H^1(\mathrm{\Omega })$. Let $\{{\varphi }_{ij}:j=1,2,\dots ,\dim (E({\mu }_i))\}$ be an orthonormal basis of $E({\mu }_i)$, and $\mathbf{X}=[H^1(\mathrm{\Omega }){]}^3$, ${\mathbf{X}}_{ij}=\{\mathbf{c}{\varphi }_{ij}:\mathbf{c}\in {\mathbb{R}}^3\}$. Then (23) $\mathbf{X}=\underset{i=1}{\overset{+\mathrm{\infty }}{⨁}}{\mathbf{X}}_i,\mathrm{and}{\mathbf{X}}_i=\underset{j=1}{\overset{\dim E({\mu }_i)}{⨁}}{\mathbf{X}}_{ij}.$(23) Suppose that $(H_4)-(H_7)$ hold, let $\mathcal{L}=\mathbf{D}\mathrm{\Delta }+{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })$. The linearization system of (8(8) $\{\begin{array}{ll}{\mathbf{U}}_t-\mathbf{D}\mathrm{\Delta }\mathbf{U}=\mathbf{J}(\mathbf{U}),& x\in \mathrm{\Omega },t>0,\\ {\displaystyle \frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ \mathbf{U}(x,0)\ge 0,& x\in \mathrm{\Omega }.\end{array}$(8) ) at ${\mathbf{U}}^{\ast }$ is ${\mathbf{U}}_t=\mathcal{L}\mathbf{U}$. For each $i\ge 1,{\mathbf{X}}_i$ is invariant under the operator $\mathcal{L}$, and λ is an eigenvalue of $\mathcal{L}$ if and only if it is an eigenvalue of the matrix $-{\mu }_i\mathbf{D}+{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })$ for some $i\ge 1$, in which there is an eigenvector in ${\mathbf{X}}_i$. The characteristic polynomial of $-{\mu }_i\mathbf{D}+{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })$ is given by ${\phi }_i(\lambda )={\lambda }^3+H_{1i}{\lambda }^2+H_{2i}\lambda +H_{3i},$where $\begin{array}{rl}{H}_{1i}& =(d_{11}+d_{22}+d_{33}){\mu }_i+G_5>0,\\ H_{2i}& =(d_{11}{d}_{22}+d_{11}{d}_{33}+d_{22}{d}_{33}){\mu }_i^2-[D_{33}{d}_{11}+(D_{11}+D_{33})d_{22}+D_{11}{d}_{33}]{\mu }_i+G_6,\\ H_{3i}& =d_{11}{d}_{22}{d}_{33}{\mu }_i^3-(D_{33}{d}_{11}+D_{11}{d}_{33})d_{22}{\mu }_i^2-[D_{23}{D}_{32}{d}_{11}+D_{12}{D}_{21}{d}_{33}\\ & +(D_{13}{D}_{31}-D_{11}{D}_{33})d_{22}]{\mu }_i+G_7,\end{array}$here $D_{ij}$ and $G_i$ are given in (12(12) ${\lambda }^3+G_5{\lambda }^2+G_6\lambda +G_7=0,$(12) ). If $d_{33}\ge d_{11}$, then it follows that $H_{1i},H_{2i},H_{3i}>0$. Furthermore, we have (24) $H_{1i}{H}_{2i}-H_{3i}=G_{11}{\mu }_i^3+G_{12}{\mu }_i^2+G_{13}{\mu }_i+G_5{G}_6-G_7,$(24) where $\begin{array}{rl}{G}_{11}& =(d_{11}+d_{22}+d_{33})(d_{11}{d}_{22}+d_{11}{d}_{33}+d_{22}{d}_{33})-d_{11}{d}_{22}{d}_{33}>0,\\ G_{12}& =(d_{11}{d}_{33}+d_{22}^2+2d_{11}{d}_{22}+2d_{22}{d}_{33})G_5-(d_{11}+d_{33})(D_{33}{d}_{11}+D_{11}{d}_{33})>0,\\ G_{13}& =(D_{11}+D_{33})[D_{33}{d}_{11}+D_{11}{d}_{33}+(D_{11}+D_{33})d_{22}]+G_6(d_{11}+d_{33})\\ & +D_{23}{D}_{32}{d}_{11}+D_{12}{D}_{21}{d}_{33}-(D_{12}{D}_{21}+D_{23}{D}_{32})d_{22}.\end{array}$$G_{13}>0$ if $d_{22}\ge d_{11},d_{33}$. Recalling that $G_5{G}_6-G_7>0$, we conclude that $H_{1i}{H}_{2i}-H_{3i}>0$ for all $i\ge 0$. From Routh–Hurwitz criterion, for each $i\ge 1$, ${\lambda }_{i,1},{\lambda }_{i,2},{\lambda }_{i,3}$ of ${\phi }_i=0$ have negative real parts.

In the following, we will prove that there exists a positive constant δ such that (25) $\mathrm{Re}\{{\lambda }_{i,1}\},\mathrm{Re}\{{\lambda }_{i,2}\},\mathrm{Re}\{{\lambda }_{i,3}\}\le -\delta ,\mathrm{\forall }i\ge 1.$(25) Let $\lambda ={\mu }_i\rho$, we have ${\phi }_i(\lambda )={\mu }_i^3{\rho }^3+H_{1i}{\mu }_i^2{\rho }^2+H_{2i}{\mu }_i\rho +H_{3i}\triangleq {\tilde{\phi }}_i(\rho ).$Furthermore, it follows that $\begin{array}{rl}\underset{i\to \mathrm{\infty }}{lim}\frac{{\tilde{\phi }}_i(\rho )}{{\mu }_i^3}& ={\rho }^3+(d_{11}+d_{22}+d_{33}){\rho }^2+(d_{11}{d}_{22}+d_{11}{d}_{33}+d_{22}{d}_{33})\rho +d_{11}{d}_{22}{d}_{33}\\ & \triangleq \tilde{\phi }(\rho ).\end{array}$Applying the Routh–Hurwitz criterion, it follows that the three roots ${\rho }_1,{\rho }_2,{\rho }_3$ of $\tilde{\phi }(\rho )=0$ have negative real parts. Thus there exists a positive constant $\tilde{\delta }$ such that $\mathrm{Re}\{{\rho }_1\},\mathrm{Re}\{{\rho }_2\},\mathrm{Re}\{{\rho }_3\}\le -\tilde{\delta }$. By continuity, we can find that there exists $i_0$ such that the three roots ${\rho }_{i,1},{\rho }_{i,2},{\rho }_{i,3}$ of ${\tilde{\phi }}_i(\rho )=0$ satisfy $\mathrm{Re}\{{\rho }_{i,1}\},\mathrm{Re}\{{\rho }_{i,2}\},\mathrm{Re}\{{\rho }_{i,3}\}\le -\frac{\tilde{\delta }}{2},\mathrm{\forall }i\ge i_0$. In turn, $\mathrm{Re}\{{\lambda }_{i,1}\},\mathrm{Re}\{{\lambda }_{i,2}\},\mathrm{Re}\{{\lambda }_{i,3}\}\le -\frac{{\mu }_i\tilde{\delta }}{2}\le -\frac{\tilde{\delta }}{2},\mathrm{\forall }i\ge i_0.$Let $\underset{1\le i\le i_0}{max}\{\mathrm{Re}\{{\lambda }_{i,1}\},\mathrm{Re}\{{\lambda }_{i,2}\},\mathrm{Re}\{{\lambda }_{i,3}\}\}=-\overline{\delta }$, then $\overline{\delta }>0$, and Equation (25(25) $\mathrm{Re}\{{\lambda }_{i,1}\},\mathrm{Re}\{{\lambda }_{i,2}\},\mathrm{Re}\{{\lambda }_{i,3}\}\le -\delta ,\mathrm{\forall }i\ge 1.$(25) ) holds for $\delta =min\{\overline{\delta },\frac{\tilde{\delta }}{2}\}$. Hence, the spectrum of $\mathcal{L}$ lies in $\mathrm{Re}\lambda \le -\delta$. By Theorem 5.1.1 in reference [19], we have the local stability of ${\mathbf{U}}^{\mathbf{\ast }}$.

Theorem 3.1

Suppose that $(H_4),(H_7),G_5{G}_6-G_7>0$ and $d_{22}\ge d_{33}\ge d_{11}$ are true, the co-existence equilibrium $E^{\ast }(X^{\ast },Y^{\ast },Z^{\ast })$ of model (8(8) $\{\begin{array}{ll}{\mathbf{U}}_t-\mathbf{D}\mathrm{\Delta }\mathbf{U}=\mathbf{J}(\mathbf{U}),& x\in \mathrm{\Omega },t>0,\\ {\displaystyle \frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ \mathbf{U}(x,0)\ge 0,& x\in \mathrm{\Omega }.\end{array}$(8) ) is local asymptotically stable if $m_1\in (m_1^{\ast },1)$, but is unstable if $m_1\in (0,m_1^{\ast })$.

The Routh–Hurwitz criterion implies that $E^{\ast }$ is unstable in model (8(8) $\{\begin{array}{ll}{\mathbf{U}}_t-\mathbf{D}\mathrm{\Delta }\mathbf{U}=\mathbf{J}(\mathbf{U}),& x\in \mathrm{\Omega },t>0,\\ {\displaystyle \frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ \mathbf{U}(x,0)\ge 0,& x\in \mathrm{\Omega }.\end{array}$(8) ) if $H_{3i}<0$ or $H_{1i}{H}_{2i}-H_{3i}<0$. Consequently, to discuss whether $E^{\ast }$ is stable (i.e. $m_1\in (m_1^{\ast },1)$) in model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) but unstable in partial differential equation (PDE) model (8(8) $\{\begin{array}{ll}{\mathbf{U}}_t-\mathbf{D}\mathrm{\Delta }\mathbf{U}=\mathbf{J}(\mathbf{U}),& x\in \mathrm{\Omega },t>0,\\ {\displaystyle \frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ \mathbf{U}(x,0)\ge 0,& x\in \mathrm{\Omega }.\end{array}$(8) ), we consider the following limits: $\begin{array}{rl}{L}_1& :=\underset{d_{11}\to \mathrm{\infty }}{lim}\frac{H_{2i}}{d_{11}}=(d_{22}+d_{33}){\mu }_i^2-D_{33}{\mu }_i,\\ L_2& :=\underset{d_{11}\to \mathrm{\infty }}{lim}\frac{H_{3i}}{d_{11}}=d_{22}{d}_{33}{\mu }_i^3-D_{33}{d}_{22}{\mu }_i^2-D_{23}{D}_{32}{d}_{11}{\mu }_i,\\ L_3& :=\underset{d_{11}\to \mathrm{\infty }}{lim}\frac{H_{1i}{H}_{2i}-H_{3i}}{d_{11}}=(d_{22}+d_{33}){\mu }_i^3-D_{33}{\mu }_i^2,\end{array}$where $D_{ij}(i,j=1,2,3)$ is given in (12(12) ${\lambda }^3+G_5{\lambda }^2+G_6\lambda +G_7=0,$(12) ). Suppose that ${\mu }_i>0$, we get $L_2>0$ if $d_{33}>\frac{c^2{\beta }_3{Y}^{\ast }{Z}^{\ast }}{4({\alpha }_3+2c{Z}^{\ast })({\alpha }_3+c{Z}^{\ast })}{d}_{22}$. $L_1,L_3<0$ if $(d_{22}+d_{33}){\mu }_i-D_{33}<0$, i.e. $m_1^{\ast }<m_1<{\tilde{m}}_1:=\{m_1|(d_{22}+d_{33}){\mu }_i-c{\beta }_3{Y}^{\ast }(m_1)Z^{\ast }(m_1)=0\}.$ Combining with $d_{11}$ large enough, we have $H_{1i}{H}_{2i}-H_{3i}<0$. Therefore, diffusion-driven Turing instability occurs when $d_{11}$ is large enough. In addition, we define ${\overline{m}}_1:=\{m_1|H_{1i}{H}_{2i}-H_{3i}=0\}$ and ${\hat{m}}_1=min\{{\tilde{m}}_1,{\overline{m}}_1\}$.

Theorem 3.2

If the assumptions $(H_4),(H_7)$ and $G_5{G}_6-G_7>0$ hold, then there exists an unbounded region $\mathrm{△}:=\{(d_{11},d_{22},d_{33})|d_{11}\gg d_{33}>\frac{c^2{\beta }_3{Y}^{\ast }{Z}^{\ast }}{4({\alpha }_3+2c{Z}^{\ast })({\alpha }_3+c{Z}^{\ast })}{d}_{22}\},$such that for any $(d_{11},d_{22},d_{33})\in \mathrm{△}$, $E^{\ast }$ is local asymptotically stable when $m_1\in ({\hat{m}}_1,1)$ in reaction–diffusion model (8(8) $\{\begin{array}{ll}{\mathbf{U}}_t-\mathbf{D}\mathrm{\Delta }\mathbf{U}=\mathbf{J}(\mathbf{U}),& x\in \mathrm{\Omega },t>0,\\ {\displaystyle \frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ \mathbf{U}(x,0)\ge 0,& x\in \mathrm{\Omega }.\end{array}$(8) ), but is unstable when $m_1\in (m_1^{\ast },{\hat{m}}_1)$. That is, Turing instability occurs.

3.2. A-priori bounds for the positive steady state

The corresponding steady-state problem of (8(8) $\{\begin{array}{ll}{\mathbf{U}}_t-\mathbf{D}\mathrm{\Delta }\mathbf{U}=\mathbf{J}(\mathbf{U}),& x\in \mathrm{\Omega },t>0,\\ {\displaystyle \frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ \mathbf{U}(x,0)\ge 0,& x\in \mathrm{\Omega }.\end{array}$(8) ) is (26) $\{\begin{array}{ll}-d_{11}\mathrm{\Delta }X=J_1(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{22}\mathrm{\Delta }Y=J_2(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{33}\mathrm{\Delta }Z=J_3(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(26) In the following, the generic constants $C_1,C_2,C_{\ast },\underset{\_}{C},\overline{C},$ etc., will depend on the domain Ω and the dimension N. However, when Ω and N are fixed, we will not mention this dependence explicitly. For convenience, we denote the constants $r,a,b,c,d_i,{\alpha }_i,{\beta }_i,m_j,(i=1,2,3,j=1,2)$ collectively by Λ. The main purpose of this section is to give a-priori positive upper and lower bounds for the positive solutions of (26(26) $\{\begin{array}{ll}-d_{11}\mathrm{\Delta }X=J_1(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{22}\mathrm{\Delta }Y=J_2(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{33}\mathrm{\Delta }Z=J_3(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(26) ). First, we need two results as follows.

Lemma 3.3

Maximum Principle [27]

Let $g(x,\omega )\in C(\mathrm{\Omega }\times {\mathbb{R}}^1)$ and $b_j(x)\in C(\overline{\mathrm{\Omega }}),j=1,2,\dots ,N$.

(1)

If $\omega \in C^2(\mathrm{\Omega })\cap C^1(\overline{\mathrm{\Omega }})$ satisfies $\{\begin{array}{ll}{\displaystyle -\mathrm{\Delta }\omega (x)\le \sum _{j=1}^N{b}_i(x){\omega }_{x_j}+g(x,\omega (x))}& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }\omega }{\mathrm{\partial }\nu }\le 0}& x\in \mathrm{\partial }\mathrm{\Omega },\end{array}$and $\omega (x_0)=\underset{\overline{\mathrm{\Omega }}}{max}\omega (x)$, then $g(x_0,\omega (x_0))\ge 0$.

(2)

If $\omega \in C^2(\mathrm{\Omega })\cap C^1(\overline{\mathrm{\Omega }})$ satisfies $\{\begin{array}{ll}{\displaystyle -\mathrm{\Delta }\omega (x)\ge \sum _{j=1}^N{b}_i(x){\omega }_{x_j}+g(x,\omega (x))}& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }\omega }{\mathrm{\partial }\nu }\ge 0}& x\in \mathrm{\partial }\mathrm{\Omega },\end{array}$and $\omega (x_0)=\underset{\overline{\mathrm{\Omega }}}{min}\omega (x)$, then $g(x_0,\omega (x_0))\le 0$.

Lemma 3.4

Harnack Inequality [28]

Let $\omega \in C^2(\mathrm{\Omega })\cap C^1(\overline{\mathrm{\Omega }})$ be a positive solution to $-\mathrm{\Delta }\omega (x)=c(x)\omega (x)$ with $c\in C(\overline{\mathrm{\Omega }})$, subject to homogeneous Neumann boundary condition. Then there exists a positive constant $C=C_{\ast }(N,\mathrm{\Omega },\Vert c{\Vert }_{\mathrm{\infty }})$ such that $\underset{\overline{\mathrm{\Omega }}}{max}\omega \le C_{\ast }\underset{\overline{\mathrm{\Omega }}}{min}\omega .$

Consequently, the results of upper and lower bounds can be given as follows:

Theorem 3.5

Upper Bound

Let k be a fixed positive number. Then there exist positive constants $C(\mathrm{\Lambda },k)$ and $\overline{C}(\mathrm{\Lambda },k)$ such that, when $d_{ii}\ge k,i=1,2,3,$ the positive solution $(X,Y,Z)$ to (26(26) $\{\begin{array}{ll}-d_{11}\mathrm{\Delta }X=J_1(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{22}\mathrm{\Delta }Y=J_2(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{33}\mathrm{\Delta }Z=J_3(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(26) ) satisfies (27) $\begin{array}{rl}& \underset{\overline{\mathrm{\Omega }}}{max}X\le C(\mathrm{\Lambda },k)\underset{\overline{\mathrm{\Omega }}}{min}X,\underset{\overline{\mathrm{\Omega }}}{max}Y\le C(\mathrm{\Lambda },k)\underset{\overline{\mathrm{\Omega }}}{min}Y,\underset{\overline{\mathrm{\Omega }}}{max}Z\le C(\mathrm{\Lambda },k)\underset{\overline{\mathrm{\Omega }}}{min}Z,\\ & \Vert X{\Vert }_{C^{2+\alpha }(\overline{\mathrm{\Omega }})}\le \overline{C}(\mathrm{\Lambda },k),\Vert Y{\Vert }_{C^{2+\alpha }(\overline{\mathrm{\Omega }})}\le \overline{C}(\mathrm{\Lambda },k),\Vert Z{\Vert }_{C^{2+\alpha }(\overline{\mathrm{\Omega }})}\le \overline{C}(\mathrm{\Lambda },k).\end{array}$(27)

Proof.

First of all, we will prove that this estimates (28) $\underset{\overline{\mathrm{\Omega }}}{max}X\le \overline{C}(\mathrm{\Lambda },k),\underset{\overline{\mathrm{\Omega }}}{max}Y\le \overline{C}(\mathrm{\Lambda },k),\underset{\overline{\mathrm{\Omega }}}{max}Z\le \overline{C}(\mathrm{\Lambda },k).$(28) Applying Lemma 3.3 to the first equation of (26(26) $\{\begin{array}{ll}-d_{11}\mathrm{\Delta }X=J_1(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{22}\mathrm{\Delta }Y=J_2(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{33}\mathrm{\Delta }Z=J_3(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(26) ), it yields $X\le \frac{r}{a}$. Define $K(x)=d_{11}X+\xi d_{22}Y+\eta d_{33}Z$, where $\xi =\frac{{\alpha }_1}{{\beta }_1},\eta =\frac{{\alpha }_2}{{\beta }_2}$, then we get $\{\begin{array}{ll}-\mathrm{\Delta }K=(r-d_1)X-a{X}^2-(\xi +\eta {\beta }_3)({\alpha }_3+cZ)YZ-\xi d_{2}Y-\eta d_{3}Z,& \mathrm{in}\mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }K}{\mathrm{\partial }\nu }=0,}& \mathrm{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$Let $x_0\in \mathrm{\Omega }$ such that $K(x_0)=\underset{\overline{\mathrm{\Omega }}}{max}K$. It follows from Lemma 3.3 that (29) $\begin{array}{rl}& (\xi +\eta {\beta }_3)({\alpha }_3+cZ(x_0))Y(x_0)Z(x_0)+\xi d_{2}Y(x_0)+\eta d_{3}Z(x_0)\\ & \le (r-d_1)X(x_0)-a{X}^2(x_0)\le \frac{(r-d_1{)}^2}{4a}.\end{array}$(29) A simple application of Lemma 3.3 to the third equation of (26(26) $\{\begin{array}{ll}-d_{11}\mathrm{\Delta }X=J_1(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{22}\mathrm{\Delta }Y=J_2(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{33}\mathrm{\Delta }Z=J_3(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(26) ) obtains $\underset{\overline{\mathrm{\Omega }}}{max}Z\le C_1\underset{\overline{\mathrm{\Omega }}}{min}Z$ for a positive constant $C_1=C_1(\mathrm{\Lambda },k)$. This together with (29(29) $\begin{array}{rl}& (\xi +\eta {\beta }_3)({\alpha }_3+cZ(x_0))Y(x_0)Z(x_0)+\xi d_{2}Y(x_0)+\eta d_{3}Z(x_0)\\ & \le (r-d_1)X(x_0)-a{X}^2(x_0)\le \frac{(r-d_1{)}^2}{4a}.\end{array}$(29) ) gets $\underset{\overline{\mathrm{\Omega }}}{max}Z\le C_2$ for a positive constant $C_2=C_2(\mathrm{\Lambda },k)$. Similarly, we can prove that $\underset{\overline{\mathrm{\Omega }}}{max}Y\le C_3\underset{\overline{\mathrm{\Omega }}}{min}Y,\underset{\overline{\mathrm{\Omega }}}{max}Y\le C_4$ and $\underset{\overline{\mathrm{\Omega }}}{max}X\le C_5\underset{\overline{\mathrm{\Omega }}}{min}X$ for some positive constants $C_3=C_3(\mathrm{\Lambda },k),C_4=C_4(\mathrm{\Lambda },k)$ and $C_5=C_5(\mathrm{\Lambda },k)$.

Define $C(\mathrm{\Lambda },d)=max\{C_1,C_3,C_5\}$ and $\overline{C}(\mathrm{\Lambda },d)=max\{\frac{r}{a},C_2,C_4\}$. Because of (28(28) $\underset{\overline{\mathrm{\Omega }}}{max}X\le \overline{C}(\mathrm{\Lambda },k),\underset{\overline{\mathrm{\Omega }}}{max}Y\le \overline{C}(\mathrm{\Lambda },k),\underset{\overline{\mathrm{\Omega }}}{max}Z\le \overline{C}(\mathrm{\Lambda },k).$(28) ), by the regularity for elliptic equations we have that X, Y, Z belong to $C^{1+\alpha }(\overline{\mathrm{\Omega }})$, and $C^{1+\alpha }(\overline{\mathrm{\Omega }})$ norms of them depend on known parameters $r,a,b,c,d_i,{\alpha }_i,{\beta }_i,m_j,(i=1,2,3,j=1,2)$. Using the regularity of elliptic equations again, the estimate (27(27) $\begin{array}{rl}& \underset{\overline{\mathrm{\Omega }}}{max}X\le C(\mathrm{\Lambda },k)\underset{\overline{\mathrm{\Omega }}}{min}X,\underset{\overline{\mathrm{\Omega }}}{max}Y\le C(\mathrm{\Lambda },k)\underset{\overline{\mathrm{\Omega }}}{min}Y,\underset{\overline{\mathrm{\Omega }}}{max}Z\le C(\mathrm{\Lambda },k)\underset{\overline{\mathrm{\Omega }}}{min}Z,\\ & \Vert X{\Vert }_{C^{2+\alpha }(\overline{\mathrm{\Omega }})}\le \overline{C}(\mathrm{\Lambda },k),\Vert Y{\Vert }_{C^{2+\alpha }(\overline{\mathrm{\Omega }})}\le \overline{C}(\mathrm{\Lambda },k),\Vert Z{\Vert }_{C^{2+\alpha }(\overline{\mathrm{\Omega }})}\le \overline{C}(\mathrm{\Lambda },k).\end{array}$(27) ) follows.

Theorem 3.6

Lower Bound

Let k be a fixed positive number. Then there exist a positive constant $\underset{\_}{C}=\underset{\_}{C}(\mathrm{\Lambda },k)$ such that, when $d_{ii}\ge k,i=1,2,3,$ the positive solution $(X,Y,Z)$ to (26(26) $\{\begin{array}{ll}-d_{11}\mathrm{\Delta }X=J_1(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{22}\mathrm{\Delta }Y=J_2(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{33}\mathrm{\Delta }Z=J_3(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(26) ) satisfies (30) $\underset{\overline{\mathrm{\Omega }}}{min}X\ge \underset{\_}{C},\underset{\overline{\mathrm{\Omega }}}{min}Y\ge \underset{\_}{C},\underset{\overline{\mathrm{\Omega }}}{min}Z\ge \underset{\_}{C}.$(30)

Proof.

We suppose that (30(30) $\underset{\overline{\mathrm{\Omega }}}{min}X\ge \underset{\_}{C},\underset{\overline{\mathrm{\Omega }}}{min}Y\ge \underset{\_}{C},\underset{\overline{\mathrm{\Omega }}}{min}Z\ge \underset{\_}{C}.$(30) ) does not hold. Then there exists a sequence $\{(d_{11}^j,d_{22}^j,d_{33}^j){\}}_{j=1}^{\mathrm{\infty }}$ with $d_{11}^j,d_{22}^j,d_{33}^j\in [d,\mathrm{\infty })\times [d,\mathrm{\infty })\times [d,\mathrm{\infty })$ and the corresponding positive solutions $(X_j,Y_j,Z_j)$ to (26(26) $\{\begin{array}{ll}-d_{11}\mathrm{\Delta }X=J_1(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{22}\mathrm{\Delta }Y=J_2(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{33}\mathrm{\Delta }Z=J_3(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(26) ) such that $\underset{\overline{\mathrm{\Omega }}}{min}{X}_j\to 0,\mathrm{or}\underset{\overline{\mathrm{\Omega }}}{min}{Y}_j\to 0,\mathrm{or}\underset{\overline{\mathrm{\Omega }}}{min}{Z}_j\to 0,\mathrm{as}j\to \mathrm{\infty }.$According to (27(27) $\begin{array}{rl}& \underset{\overline{\mathrm{\Omega }}}{max}X\le C(\mathrm{\Lambda },k)\underset{\overline{\mathrm{\Omega }}}{min}X,\underset{\overline{\mathrm{\Omega }}}{max}Y\le C(\mathrm{\Lambda },k)\underset{\overline{\mathrm{\Omega }}}{min}Y,\underset{\overline{\mathrm{\Omega }}}{max}Z\le C(\mathrm{\Lambda },k)\underset{\overline{\mathrm{\Omega }}}{min}Z,\\ & \Vert X{\Vert }_{C^{2+\alpha }(\overline{\mathrm{\Omega }})}\le \overline{C}(\mathrm{\Lambda },k),\Vert Y{\Vert }_{C^{2+\alpha }(\overline{\mathrm{\Omega }})}\le \overline{C}(\mathrm{\Lambda },k),\Vert Z{\Vert }_{C^{2+\alpha }(\overline{\mathrm{\Omega }})}\le \overline{C}(\mathrm{\Lambda },k).\end{array}$(27) ), this implies that (31) $\underset{\overline{\mathrm{\Omega }}}{max}{X}_j\to 0,\mathrm{or}\underset{\overline{\mathrm{\Omega }}}{max}{Y}_j\to 0,\mathrm{or}\underset{\overline{\mathrm{\Omega }}}{max}{Z}_j\to 0,\mathrm{as}j\to \mathrm{\infty }.$(31) As $d_{11}^j\ge d,d_{22}^j\ge d,d_{33}^j\ge d$, we assume $(d_{11}^j,d_{22}^j,d_{33}^j)\to ({\hat{d}}_{11},{\hat{d}}_{22},{\hat{d}}_{33})\in [d,\mathrm{\infty })\times [d,\mathrm{\infty })\times [d,\mathrm{\infty })$ as $j\to \mathrm{\infty }$. By the maximum principle, we get $X_j\le \frac{r}{a}$. Integrating by parts yields (32) $\{\begin{array}{rl}& {\int }_{\mathrm{\Omega }}{X}_j[r-a{X}_j-d_1-{\alpha }_1(1-m_1)Y_j-\frac{{\alpha }_2(1-m_2)Z_j}{1+b(1-m_2)X_j}]\mathrm{d}x=0,\\ & {\int }_{\mathrm{\Omega }}{Y}_j[{\beta }_1(1-m_1)X_j-({\alpha }_3+c{Z}_j)Z_j-d_2]\mathrm{d}x=0,\\ & {\int }_{\mathrm{\Omega }}{Z}_j[\frac{{\beta }_2(1-m_2)X_j}{1+b(1-m_2)X_j}+{\beta }_3({\alpha }_3+c{Z}_j)Y_j-d_3]\mathrm{d}x=0.\end{array}$(32) The standard regularity theorem for elliptic equations yields that there exists a subsequence of $\{(X_j,Y_j,Z_j){\}}_{j=1}^{\mathrm{\infty }}$ denoted by $\{(X_j,Y_j,Z_j){\}}_{j=1}^{\mathrm{\infty }}$, and three non-negative functions $X,Y,Z\in C^2(\overline{\mathrm{\Omega }})$, such that $(X_j,Y_j,Z_j)\to (X,Y,Z)$ in $[C^2(\overline{\mathrm{\Omega }}){]}^3$ as $j\to \mathrm{\infty }$. By (31(31) $\underset{\overline{\mathrm{\Omega }}}{max}{X}_j\to 0,\mathrm{or}\underset{\overline{\mathrm{\Omega }}}{max}{Y}_j\to 0,\mathrm{or}\underset{\overline{\mathrm{\Omega }}}{max}{Z}_j\to 0,\mathrm{as}j\to \mathrm{\infty }.$(31) ), we find that $X\equiv 0$ or $Y\equiv 0$ or $Z\equiv 0$. Let $j\to \mathrm{\infty }$ in (32(32) $\{\begin{array}{rl}& {\int }_{\mathrm{\Omega }}{X}_j[r-a{X}_j-d_1-{\alpha }_1(1-m_1)Y_j-\frac{{\alpha }_2(1-m_2)Z_j}{1+b(1-m_2)X_j}]\mathrm{d}x=0,\\ & {\int }_{\mathrm{\Omega }}{Y}_j[{\beta }_1(1-m_1)X_j-({\alpha }_3+c{Z}_j)Z_j-d_2]\mathrm{d}x=0,\\ & {\int }_{\mathrm{\Omega }}{Z}_j[\frac{{\beta }_2(1-m_2)X_j}{1+b(1-m_2)X_j}+{\beta }_3({\alpha }_3+c{Z}_j)Y_j-d_3]\mathrm{d}x=0.\end{array}$(32) ), then we have (33) $\{\begin{array}{rl}& {\int }_{\mathrm{\Omega }}X[r-aX-d_1-{\alpha }_1(1-m_1)Y-\frac{{\alpha }_2(1-m_2)Z}{1+b(1-m_2)X}]\mathrm{d}x=0,\\ & {\int }_{\mathrm{\Omega }}Y[{\beta }_1(1-m_1)X-({\alpha }_3+cZ)Z-d_2]\mathrm{d}x=0,\\ & {\int }_{\mathrm{\Omega }}Z[\frac{{\beta }_2(1-m_2)X}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)Y-d_3]\mathrm{d}x=0.\end{array}$(33) Next, we will consider three cases as follows.

Case (i). $X\equiv 0$.

Since $X_j\equiv X$, as $j\to \mathrm{\infty }$, then ${\beta }_1(1-m_1)X_j-({\alpha }_3+c{Z}_j)Z_j-d_2<0$ on $\overline{\mathrm{\Omega }}$, for all $j\gg 1$. Integrating the differential equation for $Y_j$ over Ω by parts, we have $\begin{array}{rl}-d_{22}^j{\int }_{\mathrm{\Omega }}\mathrm{\Delta }{Y}_j\mathrm{d}x& ={\int }_{\mathrm{\Omega }}{Y}_j[{\beta }_1(1-m_1)X_j-({\alpha }_3+c{Z}_j)Z_j-d_2]\mathrm{d}x<0,\\ & \text{for all}j\gg 1,\end{array}$which is a contradiction to the second equation of (32(32) $\{\begin{array}{rl}& {\int }_{\mathrm{\Omega }}{X}_j[r-a{X}_j-d_1-{\alpha }_1(1-m_1)Y_j-\frac{{\alpha }_2(1-m_2)Z_j}{1+b(1-m_2)X_j}]\mathrm{d}x=0,\\ & {\int }_{\mathrm{\Omega }}{Y}_j[{\beta }_1(1-m_1)X_j-({\alpha }_3+c{Z}_j)Z_j-d_2]\mathrm{d}x=0,\\ & {\int }_{\mathrm{\Omega }}{Z}_j[\frac{{\beta }_2(1-m_2)X_j}{1+b(1-m_2)X_j}+{\beta }_3({\alpha }_3+c{Z}_j)Y_j-d_3]\mathrm{d}x=0.\end{array}$(32) ).

Case (ii). $Y\equiv 0$.

Since $Y_j\equiv Y$, as $j\to \mathrm{\infty }$, then $\frac{{\beta }_2(1-m_2)X_j}{1+b(1-m_2)X_j}+{\beta }_3({\alpha }_3+c{Z}_j)Y_j-d_3={\beta }_3({\alpha }_3+c{Z}_j)Y_j+\frac{({\beta }_2-b{d}_3)(1-m_2)X_j-d_3}{1+b(1-m_2)X_j}{Z}_j<0$ on $\overline{\mathrm{\Omega }}$, for all $j\gg 1$, and $\begin{array}{rl}-d_{33}^j{\int }_{\mathrm{\Omega }}\mathrm{\Delta }{Z}_j\mathrm{d}x& ={\int }_{\mathrm{\Omega }}{Z}_j[\frac{{\beta }_2(1-m_2)X_j}{1+b(1-m_2)X_j}+{\beta }_3({\alpha }_3+c{Z}_j)Y_j-d_3]\mathrm{d}x<0,\\ & \text{for all}j\gg 1,\end{array}$which is a contradiction to the third equation of (32(32) $\{\begin{array}{rl}& {\int }_{\mathrm{\Omega }}{X}_j[r-a{X}_j-d_1-{\alpha }_1(1-m_1)Y_j-\frac{{\alpha }_2(1-m_2)Z_j}{1+b(1-m_2)X_j}]\mathrm{d}x=0,\\ & {\int }_{\mathrm{\Omega }}{Y}_j[{\beta }_1(1-m_1)X_j-({\alpha }_3+c{Z}_j)Z_j-d_2]\mathrm{d}x=0,\\ & {\int }_{\mathrm{\Omega }}{Z}_j[\frac{{\beta }_2(1-m_2)X_j}{1+b(1-m_2)X_j}+{\beta }_3({\alpha }_3+c{Z}_j)Y_j-d_3]\mathrm{d}x=0.\end{array}$(32) ).

Case (iii). $Z\equiv 0,X\not\equiv 0,Y\not\equiv 0$ on $\overline{\mathrm{\Omega }}$, then the Hopf lemma gives X>0, Y >0.

(a) ${\hat{d}}_{11}=\mathrm{\infty }$. Under the circumstances, X satisfies $\{\begin{array}{ll}-\mathrm{\Delta }X=0,& X\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }=0,}& X\in \mathrm{\partial }\mathrm{\Omega },\end{array}$which implies that $X=k_1$ is a positive constant. If ${\hat{d}}_{22}<\mathrm{\infty }$, then Y meets $\{\begin{array}{ll}-{\hat{d}}_{22}\mathrm{\Delta }Y=Y({\beta }_1(1-m_1)k_1-d_2),& Y\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=0,}& Y\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$When ${\beta }_1\ne \frac{d_2}{(1-m_1)k_1}$, we get $Y\equiv 0$, which is a contradiction.

If ${\beta }_1=\frac{d_2}{(1-m_1)k_1}$ or ${\hat{d}}_{22}=\mathrm{\infty }$, then $Y=k_2$ is a constant. By the third equation of (32(32) $\{\begin{array}{rl}& {\int }_{\mathrm{\Omega }}{X}_j[r-a{X}_j-d_1-{\alpha }_1(1-m_1)Y_j-\frac{{\alpha }_2(1-m_2)Z_j}{1+b(1-m_2)X_j}]\mathrm{d}x=0,\\ & {\int }_{\mathrm{\Omega }}{Y}_j[{\beta }_1(1-m_1)X_j-({\alpha }_3+c{Z}_j)Z_j-d_2]\mathrm{d}x=0,\\ & {\int }_{\mathrm{\Omega }}{Z}_j[\frac{{\beta }_2(1-m_2)X_j}{1+b(1-m_2)X_j}+{\beta }_3({\alpha }_3+c{Z}_j)Y_j-d_3]\mathrm{d}x=0.\end{array}$(32) ), there exists ${\check{x}}_j\in \mathrm{\Omega }$ such that (34) $\frac{{\beta }_2(1-m_2)X_j({\check{x}}_j)}{1+b(1-m_2)X_j({\check{x}}_j)}+{\beta }_3({\alpha }_3+c{Z}_j({\check{x}}_j))Y_j({\check{x}}_j)-d_3=0,$(34) i.e. $Y_j({\check{x}}_j)=\frac{d_3+(b{d}_3-{\beta }_2)(1-m_2)X_j({\check{x}}_j)}{{\beta }_3({\alpha }_3+c{Z}_j({\check{x}}_j))[1+b(1-m_2)X_j({\check{x}}_j)]}$. Furthermore, we assume that ${\check{x}}_j\to \check{x}$ as $j\to \mathrm{\infty }$. So Equation (34(34) $\frac{{\beta }_2(1-m_2)X_j({\check{x}}_j)}{1+b(1-m_2)X_j({\check{x}}_j)}+{\beta }_3({\alpha }_3+c{Z}_j({\check{x}}_j))Y_j({\check{x}}_j)-d_3=0,$(34) ) implies $Y(\check{x})=\frac{d_3+(b{d}_3-{\beta }_2)(1-m_2)k_1}{{\alpha }_3{\beta }_3[1+b(1-m_2)k_1]}$. Then from the first equation of (33(33) $\{\begin{array}{rl}& {\int }_{\mathrm{\Omega }}X[r-aX-d_1-{\alpha }_1(1-m_1)Y-\frac{{\alpha }_2(1-m_2)Z}{1+b(1-m_2)X}]\mathrm{d}x=0,\\ & {\int }_{\mathrm{\Omega }}Y[{\beta }_1(1-m_1)X-({\alpha }_3+cZ)Z-d_2]\mathrm{d}x=0,\\ & {\int }_{\mathrm{\Omega }}Z[\frac{{\beta }_2(1-m_2)X}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)Y-d_3]\mathrm{d}x=0.\end{array}$(33) ), we get $k_1=\frac{r-d_1-{\alpha }_1(1-m_1)k_2}{a},$which is a contradiction to Lemma 4.5 [7].

(b) ${\hat{d}}_{11}<\mathrm{\infty }$. If ${\hat{d}}_{22}<\mathrm{\infty }$, then X and Y satisfy (35) $\{\begin{array}{ll}-{\hat{d}}_{11}\mathrm{\Delta }X=X[r-aX-d_1-{\alpha }_1(1-m_1)Y],& X,Y\in \mathrm{\Omega },\\ -{\hat{d}}_{22}\mathrm{\Delta }Y=Y[{\beta }_1(1-m_1)X-d_2],& X,Y\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }={\displaystyle \frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=0,}}& X,Y\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(35) By Lemma 4.4 [7], (35(35) $\{\begin{array}{ll}-{\hat{d}}_{11}\mathrm{\Delta }X=X[r-aX-d_1-{\alpha }_1(1-m_1)Y],& X,Y\in \mathrm{\Omega },\\ -{\hat{d}}_{22}\mathrm{\Delta }Y=Y[{\beta }_1(1-m_1)X-d_2],& X,Y\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }={\displaystyle \frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=0,}}& X,Y\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(35) ) has no non-constant positive solution. Hence $(X,Y,Z{)}^T$ is a constant vector and Z = 0, but it is a contradiction to Lemma 4.5 [7].

If ${\hat{d}}_{22}=\mathrm{\infty }$, set ${\hat{Z}}_j=\frac{Z_j}{\Vert Z_j{\Vert }_{\mathrm{\infty }}}$, then $X_j,Y_j,{\hat{Z}}_j$ satisfy (36) $\{\begin{array}{rl}& -d_{11}^j\mathrm{\Delta }{X}_j=X_j[r-a{X}_j-d_1-{\alpha }_1(1-m_1)Y_j-\frac{{\alpha }_2(1-m_2){\hat{Z}}_j\Vert Z_j{\Vert }_{\mathrm{\infty }}}{1+b(1-m_2)X_j}],\\ & -d_{22}^j\mathrm{\Delta }{Y}_j=Y_j[{\beta }_1(1-m_1)X_j-({\alpha }_3+c{\hat{Z}}_j\Vert Z_j{\Vert }_{\mathrm{\infty }}){\hat{Z}}_j\Vert Z_j{\Vert }_{\mathrm{\infty }}-d_2],\\ & -d_{33}^j\mathrm{\Delta }{\hat{Z}}_j={\hat{Z}}_j[\frac{{\beta }_2(1-m_2)X_j}{1+b(1-m_2)X_j}+{\beta }_3({\alpha }_3+c{\hat{Z}}_j\Vert Z_j{\Vert }_{\mathrm{\infty }})Y_j-d_3],\\ & \frac{\mathrm{\partial }{X}_j}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }{Y}_j}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }{\hat{Z}}_j}{\mathrm{\partial }\nu }=0.\end{array}$(36) Similarly, there exists a subsequence $\{({\hat{X}}_j,{\hat{Y}}_j,{\hat{Z}}_j){\}}_{j=1}^{\mathrm{\infty }}$, denoted by $\{({\hat{X}}_j,{\hat{Y}}_j,{\hat{Z}}_j){\}}_{j=1}^{\mathrm{\infty }}$, and three non-negative functions $(\hat{X},\hat{Y},\hat{Z})\in [C^2(\overline{\mathrm{\Omega }}){]}^2$, such that $({\hat{X}}_j,{\hat{Y}}_j,{\hat{Z}}_j)\to (\hat{X},\hat{Y},\hat{Z})$ in $[C^2(\overline{\mathrm{\Omega }}){]}^2$. Due to $\Vert {\hat{Z}}_j{\Vert }_{\mathrm{\infty }}=1$, we obtain that $\Vert \hat{Z}{\Vert }_{\mathrm{\infty }}=1$. In addition, we notice that $\underset{j\to \mathrm{\infty }}{lim}\Vert Z_j{\Vert }_{\mathrm{\infty }}=0$. Therefore, let $j\to \mathrm{\infty }$ in (36(36) $\{\begin{array}{rl}& -d_{11}^j\mathrm{\Delta }{X}_j=X_j[r-a{X}_j-d_1-{\alpha }_1(1-m_1)Y_j-\frac{{\alpha }_2(1-m_2){\hat{Z}}_j\Vert Z_j{\Vert }_{\mathrm{\infty }}}{1+b(1-m_2)X_j}],\\ & -d_{22}^j\mathrm{\Delta }{Y}_j=Y_j[{\beta }_1(1-m_1)X_j-({\alpha }_3+c{\hat{Z}}_j\Vert Z_j{\Vert }_{\mathrm{\infty }}){\hat{Z}}_j\Vert Z_j{\Vert }_{\mathrm{\infty }}-d_2],\\ & -d_{33}^j\mathrm{\Delta }{\hat{Z}}_j={\hat{Z}}_j[\frac{{\beta }_2(1-m_2)X_j}{1+b(1-m_2)X_j}+{\beta }_3({\alpha }_3+c{\hat{Z}}_j\Vert Z_j{\Vert }_{\mathrm{\infty }})Y_j-d_3],\\ & \frac{\mathrm{\partial }{X}_j}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }{Y}_j}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }{\hat{Z}}_j}{\mathrm{\partial }\nu }=0.\end{array}$(36) ), then we get (37) $\{\begin{array}{rl}& -{\hat{d}}_{11}\mathrm{\Delta }X=X[r-aX-d_1-{\alpha }_1(1-m_1)Y],\\ & -\mathrm{\Delta }Y=0,\\ & -{\hat{d}}_{33}\mathrm{\Delta }\hat{Z}=\hat{Z}[\frac{{\beta }_2(1-m_2)X}{1+b(1-m_2)X}+{\alpha }_3{\beta }_{3}Y-d_3],\\ & \frac{\mathrm{\partial }{X}_j}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }{Y}_j}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }{\hat{Z}}_j}{\mathrm{\partial }\nu }=0.\end{array}$(37) If ${\hat{d}}_{33}=\mathrm{\infty }$, then $\hat{Z}$ satisfies $-\mathrm{\Delta }\hat{Z}=0\hat{Z}\in \mathrm{\Omega },\frac{\mathrm{\partial }\hat{Z}}{\mathrm{\partial }\nu }=0,\hat{Z}\in \mathrm{\partial }\mathrm{\Omega }$, and $\Vert \hat{Z}{\Vert }_{\mathrm{\infty }}=1$, which means that $\hat{Z}=1$. If ${\hat{d}}_{33}<\mathrm{\infty }$, an application of strong maximum principle and the Hopf boundary lemma to the third equation of (37(37) $\{\begin{array}{rl}& -{\hat{d}}_{11}\mathrm{\Delta }X=X[r-aX-d_1-{\alpha }_1(1-m_1)Y],\\ & -\mathrm{\Delta }Y=0,\\ & -{\hat{d}}_{33}\mathrm{\Delta }\hat{Z}=\hat{Z}[\frac{{\beta }_2(1-m_2)X}{1+b(1-m_2)X}+{\alpha }_3{\beta }_{3}Y-d_3],\\ & \frac{\mathrm{\partial }{X}_j}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }{Y}_j}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }{\hat{Z}}_j}{\mathrm{\partial }\nu }=0.\end{array}$(37) ) gives $\hat{Z}>0$ on $\overline{\mathrm{\Omega }}$. The above analysis shows that $\hat{Z}>0$ on $\overline{\mathrm{\Omega }}$ for both cases ${\hat{d}}_{33}=\mathrm{\infty }$ and ${\hat{d}}_{33}<\mathrm{\infty }$.

Notice that $Y\equiv k_2$ if ${\hat{d}}_{22}=\mathrm{\infty }$. The third equation of (37(37) $\{\begin{array}{rl}& -{\hat{d}}_{11}\mathrm{\Delta }X=X[r-aX-d_1-{\alpha }_1(1-m_1)Y],\\ & -\mathrm{\Delta }Y=0,\\ & -{\hat{d}}_{33}\mathrm{\Delta }\hat{Z}=\hat{Z}[\frac{{\beta }_2(1-m_2)X}{1+b(1-m_2)X}+{\alpha }_3{\beta }_{3}Y-d_3],\\ & \frac{\mathrm{\partial }{X}_j}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }{Y}_j}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }{\hat{Z}}_j}{\mathrm{\partial }\nu }=0.\end{array}$(37) ) means $k_2=\frac{d_3+(b{d}_3-{\beta }_2)k_1}{{\alpha }_3{\beta }_3[1+b(1-m_2)k_1]}$. Then X satisfies a logistic-type elliptic equation $-d\mathrm{\Delta }u=u[r-au-d_1-{\alpha }_1(1-m_1)k_2]$, which makes X to be a constant $k_1$, a contradiction to Lemma 4.5 [7]. Therefore, the proof of this theorem is completed.

3.3. Non-existence of non-constant positive steady states

In this section, we apply the energy method to prove the non-existence of the non-constant positive steady states of (26(26) $\{\begin{array}{ll}-d_{11}\mathrm{\Delta }X=J_1(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{22}\mathrm{\Delta }Y=J_2(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{33}\mathrm{\Delta }Z=J_3(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(26) ) by the effect of large diffusivity of the three species, which manifests that the three species can disperse fast. This means that the three populations will be spatially homogeneously distributed.

Theorem 3.7

If there exist three constants $d_{11}^{\ast },d_{22}^{\ast },d_{33}^{\ast }$ such that $d_{11}>d_{11}^{\ast },d_{22}>d_{22}^{\ast }$ and $d_{33}>d_{33}^{\ast }$, then the problem (26(26) $\{\begin{array}{ll}-d_{11}\mathrm{\Delta }X=J_1(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{22}\mathrm{\Delta }Y=J_2(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{33}\mathrm{\Delta }Z=J_3(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(26) ) has no non-constant solution.

Proof.

Assume that $(X,Y,Z)$ is a positive solution to (26(26) $\{\begin{array}{ll}-d_{11}\mathrm{\Delta }X=J_1(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{22}\mathrm{\Delta }Y=J_2(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{33}\mathrm{\Delta }Z=J_3(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(26) ) and denote $\overline{X}=\frac{1}{|\mathrm{\Omega }|}{\int }_{\mathrm{\Omega }}Xdx$, $\overline{Y}=\frac{1}{|\mathrm{\Omega }|}{\int }_{\mathrm{\Omega }}Ydx$ and $\overline{Z}=\frac{1}{|\mathrm{\Omega }|}{\int }_{\mathrm{\Omega }}Zdx$. Multiplying the first equation of (26(26) $\{\begin{array}{ll}-d_{11}\mathrm{\Delta }X=J_1(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{22}\mathrm{\Delta }Y=J_2(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{33}\mathrm{\Delta }Z=J_3(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(26) ) and then integrating on Ω by parts yields $\begin{array}{rl}& d_{11}{\int }_{\mathrm{\Omega }}|\mathrm{\nabla }(X-\overline{X}){|}^2\mathrm{d}x\\ & ={\int }_{\mathrm{\Omega }}\{[(r-d_1)-a(X+\overline{X})](X-\overline{X}{)}^2-{\alpha }_1(1-m_1)(XY-\overline{X}\overline{Y})(X-\overline{X})\\ & -\frac{{\alpha }_2(1-m_2)(XZ-\overline{X}\overline{Z})(X-\overline{X})+b{\alpha }_2(1-m_2{)}^{2}X\overline{X}(X-\overline{X})(Z-\overline{Z})}{[1+b(1-m_2)X][1+b(1-m_2)\overline{X}]}\}\mathrm{d}x\\ & \le (r+{\alpha }_1{C}_4+{\alpha }_2{C}_2){\int }_{\mathrm{\Omega }}(X-\overline{X}{)}^2\mathrm{d}x+\frac{r{\alpha }_1}{a}{\int }_{\mathrm{\Omega }}|X-\overline{X}||Y-\overline{Y}|\mathrm{d}x\\ & +\frac{r{\alpha }_2(a+br)}{a^2}{\int }_{\mathrm{\Omega }}|X-\overline{X}||Z-\overline{Z}|\mathrm{d}x.\end{array}$By the similar method, we can get $\begin{array}{rl}{d}_{22}{\int }_{\mathrm{\Omega }}|\mathrm{\nabla }(Y-\overline{Y}){|}^2\mathrm{d}x& \le (\frac{r{\beta }_2}{a}+d_2+{\alpha }_3{C}_2+c{C}_2^2){\int }_{\mathrm{\Omega }}(Y-\overline{Y}{)}^2\mathrm{d}x\\ & +{\beta }_1{C}_4{\int }_{\mathrm{\Omega }}|X-\overline{X}||Y-\overline{Y}|\mathrm{d}x\\ & +{\alpha }_3{C}_4{\int }_{\mathrm{\Omega }}|Y-\overline{Y}||Z-\overline{Z}|\mathrm{d}x,\\ d_{33}{\int }_{\mathrm{\Omega }}|\mathrm{\nabla }(Z-\overline{Z}){|}^2\mathrm{d}x& \le (\frac{r{\beta }_2(a+br)}{a^2}+d_3+{\alpha }_3{\beta }_3{C}_4){\int }_{\mathrm{\Omega }}(Z-\overline{Z}{)}^2\mathrm{d}x\\ & +{\beta }_2{C}_2{\int }_{\mathrm{\Omega }}|X-\overline{X}||Z-\overline{Z}|\mathrm{d}x\\ & +({\alpha }_3{\beta }_3{C}_2+c{\beta }_3{C}_2^2){\int }_{\mathrm{\Omega }}|Y-\overline{Y}||Z-\overline{Z}|\mathrm{d}x.\end{array}$From the above results and the ε-Young inequality, it follows that $\begin{array}{rl}& d_{11}{\int }_{\mathrm{\Omega }}|\mathrm{\nabla }(X-\overline{X}){|}^2\mathrm{d}x+d_{22}{\int }_{\mathrm{\Omega }}|\mathrm{\nabla }(Y-\overline{Y}){|}^2\mathrm{d}x+d_{33}{\int }_{\mathrm{\Omega }}|\mathrm{\nabla }(Z-\overline{Z}){|}^2\mathrm{d}x\\ & \le (r+{\alpha }_1{C}_4+{\alpha }_2{C}_2){\int }_{\mathrm{\Omega }}(X-\overline{X}{)}^2\mathrm{d}x+(\frac{r{\alpha }_1}{a}+{\beta }_1{C}_4){\int }_{\mathrm{\Omega }}|X-\overline{X}||Y-\overline{Y}|\mathrm{d}x\\ & +(\frac{r{\alpha }_2(a+br)}{a^2}+{\beta }_2{C}_2){\int }_{\mathrm{\Omega }}|X-\overline{X}||Z-\overline{Z}|\mathrm{d}x\\ & +({\alpha }_3{C}_4+{\alpha }_3{\beta }_3{C}_2+c{\beta }_3{C}_2^2){\int }_{\mathrm{\Omega }}|Y-\overline{Y}||Z-\overline{Z}|\mathrm{d}x\\ & +(\frac{r{\beta }_2}{a}+d_2+{\alpha }_3{C}_2+c{C}_2^2){\int }_{\mathrm{\Omega }}(Y-\overline{Y}{)}^2\mathrm{d}x\\ & +(\frac{r{\beta }_2(a+br)}{a^2}+d_3+{\alpha }_3{\beta }_3{C}_4){\int }_{\mathrm{\Omega }}(Z-\overline{Z}{)}^2\mathrm{d}x\\ & :={\rho }_1{\int }_{\mathrm{\Omega }}(X-\overline{X}{)}^2\mathrm{d}x+{\rho }_2{\int }_{\mathrm{\Omega }}|X-\overline{X}||Y-\overline{Y}|\mathrm{d}x+{\rho }_3{\int }_{\mathrm{\Omega }}|X-\overline{X}||Z-\overline{Z}|\mathrm{d}x\\ & +{\rho }_4{\int }_{\mathrm{\Omega }}|Y-\overline{Y}||Z-\overline{Z}|\mathrm{d}x+{\rho }_5{\int }_{\mathrm{\Omega }}(Y-\overline{Y}{)}^2\mathrm{d}x+{\rho }_6{\int }_{\mathrm{\Omega }}(Z-\overline{Z}{)}^2\mathrm{d}x\\ & \le {\int }_{\mathrm{\Omega }}({\rho }_1+\frac{\epsilon {\rho }_2}{2}+\frac{{\rho }_3}{2\epsilon })(X-\overline{X}{)}^2\mathrm{d}x+{\int }_{\mathrm{\Omega }}({\rho }_5+\frac{{\rho }_2}{2\epsilon }+\frac{{\rho }_4}{2\epsilon })(Y-\overline{Y}{)}^2\mathrm{d}x\\ & +{\int }_{\mathrm{\Omega }}({\rho }_6+\frac{{\rho }_3\epsilon }{2}+\frac{{\rho }_4\epsilon }{2})(Z-\overline{Z}{)}^2\mathrm{d}x,\end{array}$where ε is a small enough positive number. Using the Poinc$\acute{a}$re inequality, we get $\begin{array}{rl}& {\int }_{\mathrm{\Omega }}{d}_{11}\mu |\mathrm{\nabla }(X-\overline{X}){|}^2\mathrm{d}x+{\int }_{\mathrm{\Omega }}{d}_{22}\mu |\mathrm{\nabla }(Y-\overline{Y}){|}^2\mathrm{d}x+{\int }_{\mathrm{\Omega }}{d}_{33}\mu |\mathrm{\nabla }(Z-\overline{Z}){|}^2\mathrm{d}x\\ & \le {\int }_{\mathrm{\Omega }}{d}_{11}|\mathrm{\nabla }(X-\overline{X}){|}^2\mathrm{d}x+{\int }_{\mathrm{\Omega }}{d}_{22}|\mathrm{\nabla }(Y-\overline{Y}){|}^2\mathrm{d}x+{\int }_{\mathrm{\Omega }}{d}_{33}|\mathrm{\nabla }(Z-\overline{Z}){|}^2\mathrm{d}x.\end{array}$Let $\begin{array}{rl}{d}_{11}^{\ast }& =\frac{1}{\mu }({\rho }_1+\frac{\epsilon {\rho }_2}{2}+\frac{{\rho }_3}{2\epsilon }),d_{22}^{\ast }=\frac{1}{\mu }({\rho }_5+\frac{{\rho }_2}{2\epsilon }+\frac{{\rho }_4}{2\epsilon }),\\ d_{33}^{\ast }& =\frac{1}{\mu }({\rho }_6+\frac{{\rho }_3\epsilon }{2}+\frac{{\rho }_4\epsilon }{2}).\end{array}$If $d_{11}>d_{11}^{\ast },d_{22}>d_{22}^{\ast }$ and $d_{33}>d_{33}^{\ast }$, then $X=\overline{X}=$ constant, $Y=\overline{Y}=$ constant and $Z=\overline{Z}=$ constant. Hence, this proof is completed.

3.4. Existence of non-constant positive steady states

The previous part analysed the case that a non-constant positive steady states would not occur. In this part, by using the Leray–Schauder degree theory, we will discuss the existence of non-constant positive steady states of (26(26) $\{\begin{array}{ll}-d_{11}\mathrm{\Delta }X=J_1(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{22}\mathrm{\Delta }Y=J_2(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{33}\mathrm{\Delta }Z=J_3(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(26) ), that is the stationary pattern can exist for the system (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }t}& =d_{11}\mathrm{\Delta }X+rX-a{X}^2-d_{1}X-{\alpha }_1(1-m_1)XY\\ & -\frac{{\alpha }_2(1-m_2)XZ}{1+b(1-m_2)X},x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }t}& =d_{22}\mathrm{\Delta }Y+{\beta }_1(1-m_1)XY-({\alpha }_3+cZ)YZ-d_{2}Y,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }t}& =d_{33}\mathrm{\Delta }Z+\frac{{\beta }_2(1-m_2)XZ}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)YZ-d_{3}Z,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }\nu }& =\frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }\nu }=0,x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ X(x,0)& =X_0(x)\ge 0,Y(x,0)=Y_0(x)\ge 0,Z(x,0)=Z_0(x)\ge 0,x\in \mathrm{\Omega },\end{array}\end{array}$(7) ) with self-diffusion.

Let $\mathbf{X}$ be the Hilbert space as in Section 3.1, and define $\begin{array}{rl}\mathbf{W}& =\{\mathbf{U}\in [C^1(\overline{\mathrm{\Omega }}){]}^3,{\mathrm{\partial }}_{\nu }\mathbf{U}=0\mathrm{on}\mathrm{\partial }\mathrm{\Omega }\},{\mathbf{W}}^{\mathbf{+}}=\{\mathbf{U}\in \mathbf{W},X,Y,Z>0\mathrm{on}\overline{\mathrm{\Omega }}\},\\ B(C)& =\{X,Y,Z\in \mathbf{W}|C^{-1}<X,Y,Z<C\mathrm{on}\overline{\mathrm{\Omega }}\},C>0.\end{array}$Next, (26(26) $\{\begin{array}{ll}-d_{11}\mathrm{\Delta }X=J_1(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{22}\mathrm{\Delta }Y=J_2(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{33}\mathrm{\Delta }Z=J_3(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(26) ) can be written as follows: (38) $\{\begin{array}{ll}-\mathbf{D}\mathrm{\Delta }\mathbf{U}=\mathbf{J}(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(38) Thus $\mathbf{U}$ is a positive solution of (38(38) $\{\begin{array}{ll}-\mathbf{D}\mathrm{\Delta }\mathbf{U}=\mathbf{J}(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(38) ) if and only if it meets $L(d_{11},d_{22},d_{33};\mathbf{U}):=\mathbf{U}-(\mathbf{I}-\mathrm{\Delta }{)}^{-1}\{{\mathbf{D}}^{-1}\mathbf{J}(\mathbf{U})+\mathbf{U}\}=0\mathrm{in}{\mathbf{W}}^{\mathbf{+}},$where $(\mathbf{I}-\mathrm{\Delta }{)}^{-1}$ is the inverse of $\mathbf{I}-\mathrm{\Delta }$ with the homogeneous Neumann boundary condition. As $L(\cdot )$ is a compact perturbation of the identity operator, for any $B=B(C)$, the Leray–Schauder degree $\mathrm{deg}(L(\cdot ),0,B)$ is well-defined if $L(\mathbf{U})\ne 0$ on $\mathrm{\partial }B$. Furthermore, simple computation yields $D_{\mathbf{U}}L({\mathbf{U}}^{\ast })=\mathbf{I}-(\mathbf{I}-\mathrm{\Delta }{)}^{-1}\{{\mathbf{D}}^{-1}{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })+\mathbf{I}\}.$If $D_{\mathbf{U}}L({\mathbf{U}}^{\ast })$ is invertible, the index of L at ${\mathbf{U}}^{\ast }$ is defined as index $(L(\cdot ),{\mathbf{U}}^{\ast })=(-1{)}^{\kappa }$, where κ is the sum of algebraic multiplicity of negative eigenvalues of $D_{\mathbf{U}}L({\mathbf{U}}^{\ast })$.

We refer to the decomposition (23(23) $\mathbf{X}=\underset{i=1}{\overset{+\mathrm{\infty }}{⨁}}{\mathbf{X}}_i,\mathrm{and}{\mathbf{X}}_i=\underset{j=1}{\overset{\dim E({\mu }_i)}{⨁}}{\mathbf{X}}_{ij}.$(23) ) in the following discussion of the eigenvalues of $D_{\mathbf{U}}L({\mathbf{U}}^{\ast })$. First of all, for each integer $i\ge 1,1\le j\le \dim E({\mu }_i)$, ${\mathbf{X}}_{ij}$ is invariant under $D_{\mathbf{U}}L({\mathbf{U}}^{\ast })$, and λ is an eigenvalue of $D_{\mathbf{U}}L({\mathbf{U}}^{\ast })$ if and only if $\lambda (1+{\mu }_i)$ is an eigenvalue of the matrix ${\mu }_i\mathbf{I}-{\mathbf{D}}^{-1}{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })$. Thus $D_{\mathbf{U}}L({\mathbf{U}}^{\ast })$ is invertible if and only if this matrix is non-singular. Let $M(\mu ):=\det \{\mu \mathbf{I}-{\mathbf{D}}^{-1}{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })\}=\frac{1}{d_{11}{d}_{22}{d}_{33}}\det \{\mu \mathbf{D}-{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })\}.$To calculate the index $(L(\cdot ),{\mathbf{U}}^{\mathbf{\ast }})$, we cite the following lemma.

Lemma 3.8

[7]

Suppose that the matrix ${\mu }_i\mathbf{I}-{\mathbf{D}}^{-1}{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })$ is non-singular for all $i\ge 1$. Then $\mathrm{index}(L(\cdot ),{\mathbf{U}}^{\mathbf{\ast }})=(-1{)}^{\kappa },$where $\kappa =\sum _{i\ge 1,M({\mu }_i)<0}\dim E({\mu }_i)$.

The direct computation gives (39) $\begin{array}{rl}\det \{\mu \mathbf{D}-{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })\}& =W_1(d_{22}){\mu }^3+W_2(d_{22}){\mu }^2+W_3(d_{22})\mu -\det \{{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })\}\\ & :=W(d_{22};\mu ),\end{array}$(39) where $\begin{array}{rl}& W_1(d_{22})=d_{11}{d}_{22}{d}_{33},W_2(d_{22})=-(d_{11}{D}_{33}+d_{33}{D}_{11})d_{22},\\ & W_3(d_{22})=(D_{11}{D}_{33}-D_{13}{D}_{31})d_{22}-D_{23}{D}_{32}{d}_{11}-D_{12}{D}_{21}{d}_{33},\end{array}$here $D_{ij}$ are given in (12(12) ${\lambda }^3+G_5{\lambda }^2+G_6\lambda +G_7=0,$(12) ).

We find the dependence of W on $d_{22}$. Let ${\tilde{\mu }}_1(d_{22}),{\tilde{\mu }}_2(d_{22})$ and ${\tilde{\mu }}_3(d_{22})$ be the three roots of $W(d_{22};\mu )=0$ with $\mathrm{Re}\{{\tilde{\mu }}_1(d_{22})\}\le \mathrm{Re}\{{\tilde{\mu }}_2(d_{22})\}\le \mathrm{Re}\{{\tilde{\mu }}_3(d_{22})\}$. It is known that ${\tilde{\mu }}_1(d_{22}){\tilde{\mu }}_2(d_{22}){\tilde{\mu }}_3(d_{22})=\det \{{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })\}$ from (39(39) $\begin{array}{rl}\det \{\mu \mathbf{D}-{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })\}& =W_1(d_{22}){\mu }^3+W_2(d_{22}){\mu }^2+W_3(d_{22})\mu -\det \{{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })\}\\ & :=W(d_{22};\mu ),\end{array}$(39) ), and according to the previous analyses, $\det \{{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })\}<0$ can be obtained. Thus, due to $W_1(d_{22})>0$, one of ${\tilde{\mu }}_1(d_{22}),{\tilde{\mu }}_2(d_{22}),{\tilde{\mu }}_3(d_{22})$ is real and negative, and the product of the other two is positive.

Consider the following limits: $\begin{array}{rl}& \underset{d_{22}\to \mathrm{\infty }}{lim}\frac{W_1(d_{22})}{d_{22}}=d_{11}{d}_{33},\underset{d_{22}\to \mathrm{\infty }}{lim}\frac{W_2(d_{22})}{d_{22}}=-(D_{33}{d}_{11}+D_{11}{d}_{33}),\\ & \underset{d_{22}\to \mathrm{\infty }}{lim}\frac{W_3(d_{22})}{d_{22}}=D_{11}{D}_{33}-D_{13}{D}_{31},\\ & \underset{d_{22}\to \mathrm{\infty }}{lim}\frac{W(d_{22})}{d_{22}}=\mu [d_{11}{d}_{33}{\mu }^2-(D_{33}{d}_{11}+D_{11}{d}_{33})\mu +(D_{11}{D}_{33}-D_{13}{D}_{31})].\end{array}$If $(H_8):D_{11}{D}_{33}-D_{13}{D}_{31}<0$ holds, then the following lemma is true.

Lemma 3.9

Assume that $(H_8)$ holds, there exists a positive constant d such that when $d_{22}>d$, three roots ${\tilde{\mu }}_1(d_{22})$, ${\tilde{\mu }}_2(d_{22}),{\tilde{\mu }}_3(d_{22})$ of $W(d_{22};\mu )=0$ are all real and meet $\begin{array}{rl}& \underset{d_{22}\to \mathrm{\infty }}{lim}\frac{W_1(d_{22})}{d_{22}}=\frac{D_{33}{d}_{11}+D_{11}{d}_{33}-\sqrt{(D_{33}{d}_{11}-D_{11}{d}_{33}{)}^2+4d_{11}{d}_{33}{D}_{13}{D}_{31}}}{2d_{11}{d}_{33}}\triangleq \hat{\mu },\\ & \underset{d_{22}\to \mathrm{\infty }}{lim}\frac{W_2(d_{22})}{d_{22}}=0,\\ & \underset{d_{22}\to \mathrm{\infty }}{lim}\frac{W_3(d_{22})}{d_{22}}=\frac{D_{33}{d}_{11}+D_{11}{d}_{33}+\sqrt{(D_{33}{d}_{11}-D_{11}{d}_{33}{)}^2+4d_{11}{d}_{33}{D}_{13}{D}_{31}}}{2d_{11}{d}_{33}}\triangleq \check{\mu }.\end{array}$Furthermore, it obtains (40) $\{\begin{array}{rl}& -\mathrm{\infty }<{\tilde{\mu }}_1(d_{22})<0<{\tilde{\mu }}_2(d_{22})<{\tilde{\mu }}_3(d_{22}),\\ & W(d_{22};\mu )<0,\mathrm{if}\mu \in (-\mathrm{\infty },{\tilde{\mu }}_1(d_{22}))\cup ({\tilde{\mu }}_2(d_{22}),{\tilde{\mu }}_3(d_{22})),\\ & W(d_{22};\mu )>0,\mathrm{if}\mu \in ({\tilde{\mu }}_1(d_{22}),{\tilde{\mu }}_2(d_{22}))\cup ({\tilde{\mu }}_3(d_{22}),+\mathrm{\infty }).\end{array}$(40)

By determining the range of μ such that $M(\mu )<0$, we have the existence of non-constant steady states of (26(26) $\{\begin{array}{ll}-d_{11}\mathrm{\Delta }X=J_1(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{22}\mathrm{\Delta }Y=J_2(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{33}\mathrm{\Delta }Z=J_3(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(26) ).

Theorem 3.10

Suppose that $d_{11},d_{33}$ and Λ are fixed and the condition $(H_8)$ is true. If $\check{\mu }\in ({\mu }_n,{\mu }_{n+1})$ for some $n\ge 2$ and the $\sum _{i=2}^n\dim E({\mu }_i)$ is odd, then there exists a positive constant d such that (26(26) $\{\begin{array}{ll}-d_{11}\mathrm{\Delta }X=J_1(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{22}\mathrm{\Delta }Y=J_2(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{33}\mathrm{\Delta }Z=J_3(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(26) ) has at least one non-constant positive solution for $d_{22}>d$.

Proof.

By Lemma 3.9 and $\check{\mu }$, there exists a positive constant d such that (40(40) $\{\begin{array}{rl}& -\mathrm{\infty }<{\tilde{\mu }}_1(d_{22})<0<{\tilde{\mu }}_2(d_{22})<{\tilde{\mu }}_3(d_{22}),\\ & W(d_{22};\mu )<0,\mathrm{if}\mu \in (-\mathrm{\infty },{\tilde{\mu }}_1(d_{22}))\cup ({\tilde{\mu }}_2(d_{22}),{\tilde{\mu }}_3(d_{22})),\\ & W(d_{22};\mu )>0,\mathrm{if}\mu \in ({\tilde{\mu }}_1(d_{22}),{\tilde{\mu }}_2(d_{22}))\cup ({\tilde{\mu }}_3(d_{22}),+\mathrm{\infty }).\end{array}$(40) ) satisfies if $d_{22}>d$, and (41) ${\tilde{\mu }}_1(d_{22})<0={\mu }_1<{\tilde{\mu }}_2(d_{22})<{\mu }_2,{\tilde{\mu }}_3(d_{22})\in ({\mu }_n,{\mu }_{n+1}).$(41) By using the homotopic invariance of the topological degree, we prove that (26(26) $\{\begin{array}{ll}-d_{11}\mathrm{\Delta }X=J_1(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{22}\mathrm{\Delta }Y=J_2(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{33}\mathrm{\Delta }Z=J_3(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(26) ) has at least one non-constant positive solution for any $d_{22}>d$.

On the contrary, suppose that the conclusion is not true for some $d_{22}={\tilde{d}}_{22}>d$. Let $d_{22}={\tilde{d}}_{22}$. By Theorem 3.7, we can get that $d_{11}^{\ast }=\frac{{\rho }_1}{\mu }$, $d_{22}^{\ast }=\frac{{\rho }_5}{\mu }$, $d_{33}^{\ast }=\frac{{\rho }_6}{\mu }$. Then fix ${\hat{d}}_{11}>d_{11}^{\ast }$, ${\hat{d}}_{22}>d_{22}^{\ast }$, ${\hat{d}}_{33}=d_{33}^{\ast }+d_{33}$. Define $\mathbf{D}(t)=\mathrm{diag}(d_{11}(t),d_{22}(t),d_{33}(t))$ with $d_{ii}(t)=t{d}_{ii}+(1-t){\hat{d}}_{ii},i=1,2,3$ and $t\in [0,1]$, and consider the following elliptic equation problem: (42) $\{\begin{array}{ll}-\mathbf{D}(t)\mathrm{\Delta }\mathbf{U}=\mathbf{J}(\mathbf{U})& \mathrm{in}\mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }\nu }=0}& \mathrm{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(42) $\mathbf{U}$ is a non-constant positive solution of (26(26) $\{\begin{array}{ll}-d_{11}\mathrm{\Delta }X=J_1(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{22}\mathrm{\Delta }Y=J_2(\mathbf{U}),& x\in \mathrm{\Omega },\\ -d_{33}\mathrm{\Delta }Z=J_3(\mathbf{U}),& x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }X}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Y}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z}{\mathrm{\partial }\nu }=0,}& x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$(26) ) if and only if it is a solution of (42(42) $\{\begin{array}{ll}-\mathbf{D}(t)\mathrm{\Delta }\mathbf{U}=\mathbf{J}(\mathbf{U})& \mathrm{in}\mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }\nu }=0}& \mathrm{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(42) ) with t = 1. It is obvious that ${\mathbf{U}}^{\ast }$ is the unique constant positive solution of (42(42) $\{\begin{array}{ll}-\mathbf{D}(t)\mathrm{\Delta }\mathbf{U}=\mathbf{J}(\mathbf{U})& \mathrm{in}\mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }\nu }=0}& \mathrm{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(42) ) for any $0\le t\le 1$. For any $t\in [0,1]$, $\mathbf{U}$ is a non-constant positive solution of (42(42) $\{\begin{array}{ll}-\mathbf{D}(t)\mathrm{\Delta }\mathbf{U}=\mathbf{J}(\mathbf{U})& \mathrm{in}\mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }\nu }=0}& \mathrm{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(42) ) if and only if $L(t;\mathbf{U}):=\mathbf{U}-(\mathbf{I}-\mathrm{\Delta }{)}^{-1}\{{\mathbf{D}}^{-1}(t)\mathbf{J}(\mathbf{U})+\mathbf{U}\}=0\mathrm{in}{\mathbf{W}}^{\mathbf{+}}.$Obviously, $L(1;\mathbf{U})=L(\mathbf{U})$. According to Theorem 3.7, $L(0;\mathbf{U})=0$ has only the positive solution ${\mathbf{U}}^{\ast }$ in ${\mathbf{W}}^{\mathbf{+}}$. Then, we obtain $\begin{array}{rl}& D_{\mathbf{U}}L(0;{\mathbf{U}}^{\ast })=\mathbf{I}-(\mathbf{I}-\mathrm{\Delta }{)}^{-1}\{{\hat{\mathbf{D}}}^{-1}{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })+\mathbf{I}\},\\ & D_{\mathbf{U}}L(1;{\mathbf{U}}^{\ast })=\mathbf{I}-(\mathbf{I}-\mathrm{\Delta }{)}^{-1}\{{\mathbf{D}}^{-1}{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })+\mathbf{I}\}=D_{\mathbf{U}}L({\mathbf{U}}^{\ast }),\end{array}$where $\hat{\mathbf{D}}=\mathrm{diag}({\hat{d}}_{11},{\hat{d}}_{22},{\hat{d}}_{33})$. From (39(39) $\begin{array}{rl}\det \{\mu \mathbf{D}-{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })\}& =W_1(d_{22}){\mu }^3+W_2(d_{22}){\mu }^2+W_3(d_{22})\mu -\det \{{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })\}\\ & :=W(d_{22};\mu ),\end{array}$(39) ), we know that $M(\mu )=\frac{1}{d_{11}{d}_{22}{d}_{33}}W(d_{22};\mu ).$By (40(40) $\{\begin{array}{rl}& -\mathrm{\infty }<{\tilde{\mu }}_1(d_{22})<0<{\tilde{\mu }}_2(d_{22})<{\tilde{\mu }}_3(d_{22}),\\ & W(d_{22};\mu )<0,\mathrm{if}\mu \in (-\mathrm{\infty },{\tilde{\mu }}_1(d_{22}))\cup ({\tilde{\mu }}_2(d_{22}),{\tilde{\mu }}_3(d_{22})),\\ & W(d_{22};\mu )>0,\mathrm{if}\mu \in ({\tilde{\mu }}_1(d_{22}),{\tilde{\mu }}_2(d_{22}))\cup ({\tilde{\mu }}_3(d_{22}),+\mathrm{\infty }).\end{array}$(40) ) and (41(41) ${\tilde{\mu }}_1(d_{22})<0={\mu }_1<{\tilde{\mu }}_2(d_{22})<{\mu }_2,{\tilde{\mu }}_3(d_{22})\in ({\mu }_n,{\mu }_{n+1}).$(41) ), it follows that $\{\begin{array}{rl}& M({\mu }_1)=M(0)>0,\\ & M({\mu }_i)<0,2\le i\le n,\\ & M({\mu }_i)>0,i\ge n+1.\end{array}$Hence, zero is not an eigenvalue of the matrix ${\mu }_i\mathbf{I}-{\mathbf{D}}^{-1}{\mathbf{J}}_{\mathbf{U}}({\mathbf{U}}^{\ast })$ for all $i\ge 1$, and $\sum _{i\ge 1,M({\mu }_i)<0}\dim E({\mu }_i)=\sum _{i=2}^n\dim E({\mu }_i)={\kappa }_n$, which is odd. It follows from Lemma 3.8 that (43) $\mathrm{index}(L(1;\cdot ),{\mathbf{U}}^{\ast })=(-1{)}^{\kappa }=(-1{)}^{{\kappa }_n}=-1,\mathrm{index}(L(0;\cdot ),{\mathbf{U}}^{\ast })=(-1{)}^0=1.$(43) Now, by Theorems 3.6 and 3.5, there exists a positive constant C such that, for all $0\le t\le 1$, the positive solution of (42(42) $\{\begin{array}{ll}-\mathbf{D}(t)\mathrm{\Delta }\mathbf{U}=\mathbf{J}(\mathbf{U})& \mathrm{in}\mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }\nu }=0}& \mathrm{on}\mathrm{\partial }\mathrm{\Omega }.\end{array}$(42) ) satisfies $\frac{1}{C}<X,Y,Z<C$. Hence $L(t;\mathbf{U})\ne 0$ on $\mathrm{\partial }B(C)$ for all $0\le t\le 1$. By the homotopy invariance of the topological degree, we have that (44) $\mathrm{deg}(L(1;\cdot ),0,B(C))=\mathrm{deg}(L(0;\cdot ),0,B(C)).$(44) Nevertheless, both equation $L(0;\mathbf{U})=0$ and $L(1;\mathbf{U})=0$ have the unique positive solution ${\mathbf{U}}^{\ast }$ in $B(C)$. From (43(43) $\mathrm{index}(L(1;\cdot ),{\mathbf{U}}^{\ast })=(-1{)}^{\kappa }=(-1{)}^{{\kappa }_n}=-1,\mathrm{index}(L(0;\cdot ),{\mathbf{U}}^{\ast })=(-1{)}^0=1.$(43) ), we have that $\begin{array}{rl}& \mathrm{deg}(L(1;\cdot ),0,B(C))=\mathrm{index}(L(1;\cdot ),{\mathbf{U}}^{\ast })=-1,\\ & \mathrm{deg}(L(0;\cdot ),0,B(C))=\mathrm{index}(L(0;\cdot ),{\mathbf{U}}^{\ast })=1,\end{array}$which contradict (44(44) $\mathrm{deg}(L(1;\cdot ),0,B(C))=\mathrm{deg}(L(0;\cdot ),0,B(C)).$(44) ). Thus we finish the proof of this theorem.

4. Numerical simulation

In this section, with the help of Matlab software, numerical simulations will be given to illustrate our analytic results obtained in previous sections. All values of parameters in model (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }t}& =d_{11}\mathrm{\Delta }X+rX-a{X}^2-d_{1}X-{\alpha }_1(1-m_1)XY\\ & -\frac{{\alpha }_2(1-m_2)XZ}{1+b(1-m_2)X},x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }t}& =d_{22}\mathrm{\Delta }Y+{\beta }_1(1-m_1)XY-({\alpha }_3+cZ)YZ-d_{2}Y,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }t}& =d_{33}\mathrm{\Delta }Z+\frac{{\beta }_2(1-m_2)XZ}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)YZ-d_{3}Z,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }\nu }& =\frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }\nu }=0,x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ X(x,0)& =X_0(x)\ge 0,Y(x,0)=Y_0(x)\ge 0,Z(x,0)=Z_0(x)\ge 0,x\in \mathrm{\Omega },\end{array}\end{array}$(7) ) are given according to the meaning of parameters. What's more, because model (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }t}& =d_{11}\mathrm{\Delta }X+rX-a{X}^2-d_{1}X-{\alpha }_1(1-m_1)XY\\ & -\frac{{\alpha }_2(1-m_2)XZ}{1+b(1-m_2)X},x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }t}& =d_{22}\mathrm{\Delta }Y+{\beta }_1(1-m_1)XY-({\alpha }_3+cZ)YZ-d_{2}Y,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }t}& =d_{33}\mathrm{\Delta }Z+\frac{{\beta }_2(1-m_2)XZ}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)YZ-d_{3}Z,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }\nu }& =\frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }\nu }=0,x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ X(x,0)& =X_0(x)\ge 0,Y(x,0)=Y_0(x)\ge 0,Z(x,0)=Z_0(x)\ge 0,x\in \mathrm{\Omega },\end{array}\end{array}$(7) ) is not considered for the specific species a reasonable set of parameters are chosen as follows: $r=0.7,a=0.5,d_1=0.12,{\alpha }_1=1.9,{\alpha }_2=2.5,b=10,{\beta }_1=2.2,{\alpha }_3=2.57,d_2=0.2,{\beta }_2=1.1,{\beta }_3=2,d_3=0.5,m_1=0.7,m_2=0.4,c=\mathrm{0.1.}$ Since $X(t),Y(t)$ and $Z(t)$ represent the density of the shared resource, IG prey and IG predator at the time t, respectively, the initial values of system (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) are assumed to be $X(0)=0.41,Y(0)=0.21,Z(0)=0.2$.

First of all, we mainly illustrate the stability and Hopf bifurcation for non-spatial model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) by using the ode45 function of Matlab software. According to Theorem 2.6, we obtain that model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) has the positive equilibrium $E^{\ast }=(0.9947,0.0784,0.1764)$ and is locally asymptotically stable (see Figure 2). Subsequently, we choose $m_1$ of model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) as a bifurcation parameter. Based on the parameter values given above, we can get $m_1^{\ast }=0.59$. Furthermore, when $m_1=0.7>m_1^{\ast }$, $E^{\ast }$ is locally asymptotically stable (see Figure 3a), but $E^{\ast }$ is unstable when $m_1=0.54<m_1^{\ast }$ (see Figure 3b). That is, Hopf bifurcation takes place when $m_1=m_1^{\ast }$. That is, when we take measures to protect the shared resource, the dynamics of system is determined by the strength of protection. If the strength of protection is weak, periodic fluctuations occur in the system. Inverse, the system will keep well stability.

Figure 2. The co-existence equilibrium $E^{\ast }$ of model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) is locally asymptotically stable: (a) shared resource, (b) IG prey, (c) IG predator. (a) $X(t)$, (b) $Y(t)$ and (c) $Z(t)$.

Figure 3. Local dynamics of model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) at the co-existence equilibrium $E^{\ast }$: (a) $m_1=0.7>m_1^{\ast }=0.59$, (b) $m_1=0.54<m_1^{\ast }=0.59$.

Next, by choosing the $m_1=0,m_2=0,c=0$ and fixing all values of the other parameters, we obtain that the dynamical behaviour of the model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) at the positive equilibrium $E^{\ast }$ without prey refuge and hunting cooperation is shown in Figure 4. From Figure 4, the shared resource, IG prey and IG predator take periodically if IG predators lack of hunting cooperation and the shared resource is not protected.

Figure 4. The dynamical behaviour of the model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) at $E^{\ast }$ when $m_1=0,m_2=0,c=0$.

Based on Theorem 2.6, the dynamic behaviour of the system at the positive equilibrium $E^{\ast }$ without the effect of hunting cooperation is seen in Figure 5. Further, it can be concluded that the stability of model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) at $E^{\ast }$ is better with the increase of prey refuge.

Figure 5. The dynamical behaviour of model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) at $E^{\ast }$ when c = 0: (a) shared resource, (b) IG prey, (c) IG predator. (a) $X(t)$, (b) $Y(t)$ and (c) $Z(t)$.

In addition, according to Lemma 2.5, the dynamical behaviour of model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) at the positive equilibrium $E^{\ast }$ without prey refuge is seen in Figure 6. From Figure 6, the shared resource, IG prey and IG predator take periodically if the prey refuge of the shared resource is not considered. That is, we take measures to protect the resource, which is good for the long-term stability of the system.

Figure 6. The dynamics of model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) at $E^{\ast }$ when $m_1=0,m_2=0$: (a) shared resource, (b) IG prey, (c) IG predator. (a) $X(t)$, (b) $Y(t)$ and (c) $Z(t)$.

While the effect of prey refuge is considered, i.e. $m_1=0.7,m_2=0.4$, the model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) becomes stable when c is very small (see Figure 7), but becomes unstable when it gradually increases (see Figure 8). From the viewpoint of ecology, the hunting cooperation of IG predator may change the stability of system. That is, the hunting cooperation may be beneficial for IG predator. However, when the intensity of cooperation is gradually increasing, it may be bad for the system, which shows that the effect of hunting cooperation makes the model always unstable.

Figure 7. The model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) with prey refuge becomes stable when hunting cooperation c is very small: (a) shared resource, (b) IG prey, (c) IG predator. (a) $X(t)$. (b) $Y(t)$ and (c) $Z(t)$.

Figure 8. The model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) with prey refuge becomes unstable when hunting cooperation c gradually increases: (a) shared resource, (b) IG prey, (c) IG predator. (a) $X(t)$, (b) $Y(t)$ and (c) $Z(t)$.

Finally, we will check the spatial dynamics of model (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }t}& =d_{11}\mathrm{\Delta }X+rX-a{X}^2-d_{1}X-{\alpha }_1(1-m_1)XY\\ & -\frac{{\alpha }_2(1-m_2)XZ}{1+b(1-m_2)X},x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }t}& =d_{22}\mathrm{\Delta }Y+{\beta }_1(1-m_1)XY-({\alpha }_3+cZ)YZ-d_{2}Y,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }t}& =d_{33}\mathrm{\Delta }Z+\frac{{\beta }_2(1-m_2)XZ}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)YZ-d_{3}Z,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }\nu }& =\frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }\nu }=0,x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ X(x,0)& =X_0(x)\ge 0,Y(x,0)=Y_0(x)\ge 0,Z(x,0)=Z_0(x)\ge 0,x\in \mathrm{\Omega },\end{array}\end{array}$(7) ) by considering the diffusion. We will use pdepe function of Matlab software to solve PDE model and show intuitive results. We assume the spatial region Ω as an one-dimensional domain $[0,10\pi ]$ for the convenience of simulation. We keep all values of all parameters with the same of the above and choose the values of diffusions $d_{11}=10,d_{22}=20,d_{33}=15$ satisfying one condition of Theorem 3.1. Under this case, we can obtained the critical value $m_1^{\ast }=0.59$. Thus, according to Theorem 3.1, the co-existence equilibrium $E^{\ast }$ is stable when $m_1=0.7\in (m_1^{\ast },1)$, which is shown in Figure 9, but is unstable when $m_1=0.5\in (0,m_1^{\ast })$, which is shown in Figure 10. That is, it also can be concluded that the stability of model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) at $E^{\ast }$ is better with the increase of prey refuge.

Figure 9. The model (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }t}& =d_{11}\mathrm{\Delta }X+rX-a{X}^2-d_{1}X-{\alpha }_1(1-m_1)XY\\ & -\frac{{\alpha }_2(1-m_2)XZ}{1+b(1-m_2)X},x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }t}& =d_{22}\mathrm{\Delta }Y+{\beta }_1(1-m_1)XY-({\alpha }_3+cZ)YZ-d_{2}Y,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }t}& =d_{33}\mathrm{\Delta }Z+\frac{{\beta }_2(1-m_2)XZ}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)YZ-d_{3}Z,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }\nu }& =\frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }\nu }=0,x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ X(x,0)& =X_0(x)\ge 0,Y(x,0)=Y_0(x)\ge 0,Z(x,0)=Z_0(x)\ge 0,x\in \mathrm{\Omega },\end{array}\end{array}$(7) ) is stable when $m_1=0.7\in (m_1^{\ast },1)$: (a) shared resource, (b)IG prey, (c) IG predator. (a) $X(t)$, (b) $Y(t)$ and (c) $Z(t)$.

Figure 10. Unstable behaviour of model (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }t}& =d_{11}\mathrm{\Delta }X+rX-a{X}^2-d_{1}X-{\alpha }_1(1-m_1)XY\\ & -\frac{{\alpha }_2(1-m_2)XZ}{1+b(1-m_2)X},x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }t}& =d_{22}\mathrm{\Delta }Y+{\beta }_1(1-m_1)XY-({\alpha }_3+cZ)YZ-d_{2}Y,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }t}& =d_{33}\mathrm{\Delta }Z+\frac{{\beta }_2(1-m_2)XZ}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)YZ-d_{3}Z,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }\nu }& =\frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }\nu }=0,x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ X(x,0)& =X_0(x)\ge 0,Y(x,0)=Y_0(x)\ge 0,Z(x,0)=Z_0(x)\ge 0,x\in \mathrm{\Omega },\end{array}\end{array}$(7) ) when $m_1=0.5\in (0,m_1^{\ast })$: (a) shared resource, (b) IG prey, (c) IG predator. (a) $X(t)$. (b) $Y(t)$ and (c) $Z(t)$.

To check Turing instability, we choose that the values are $r=0.7,a=0.38,d_1=0.2,{\alpha }_1=2.4,{\alpha }_2=2.5,m_2=0.4,b=5,{\beta }_1=2.4,{\alpha }_3=2.7,c=0.1,d_2=0.3,{\beta }_2=0.8,{\beta }_3=1.5,d_3=0.3,d_{11}=1,d_{22}=0.05,d_{33}=0.01$. It is obtained that $m_1^{\ast }=0.502,{\hat{m}}_1=0.59$. According to Theorem 3.2, when $m_1=0.8\in ({\hat{m}}_1,1)$, the positive equilibrium $E^{\ast }$ is local asymptotically stable as shown in Figure 11, but Turing instability appears as shown in Figure 12 when $m_1=0.51\in (m_1^{\ast },{\hat{m}}_1)$. Generally speaking, the diffusion often makes the Turing instability take place. However, the prey refuge for IG prey also change the stability of the model with and without diffusion.

Figure 11. Model (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }t}& =d_{11}\mathrm{\Delta }X+rX-a{X}^2-d_{1}X-{\alpha }_1(1-m_1)XY\\ & -\frac{{\alpha }_2(1-m_2)XZ}{1+b(1-m_2)X},x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }t}& =d_{22}\mathrm{\Delta }Y+{\beta }_1(1-m_1)XY-({\alpha }_3+cZ)YZ-d_{2}Y,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }t}& =d_{33}\mathrm{\Delta }Z+\frac{{\beta }_2(1-m_2)XZ}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)YZ-d_{3}Z,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }\nu }& =\frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }\nu }=0,x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ X(x,0)& =X_0(x)\ge 0,Y(x,0)=Y_0(x)\ge 0,Z(x,0)=Z_0(x)\ge 0,x\in \mathrm{\Omega },\end{array}\end{array}$(7) ) is local asymptotically stable when $d_{11}=1$, $d_{22}=0.05$, $d_{33}=0.01$ and $m_1=0.8\in ({\hat{m}}_1,1)$: (a) shared resource, (b)IG prey, (c) IG predator. (a) $X(t)$, (b) $Y(t)$ and (c) $Z(t)$.

Figure 12. Turing instability occurs in model (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }t}& =d_{11}\mathrm{\Delta }X+rX-a{X}^2-d_{1}X-{\alpha }_1(1-m_1)XY\\ & -\frac{{\alpha }_2(1-m_2)XZ}{1+b(1-m_2)X},x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }t}& =d_{22}\mathrm{\Delta }Y+{\beta }_1(1-m_1)XY-({\alpha }_3+cZ)YZ-d_{2}Y,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }t}& =d_{33}\mathrm{\Delta }Z+\frac{{\beta }_2(1-m_2)XZ}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)YZ-d_{3}Z,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }\nu }& =\frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }\nu }=0,x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ X(x,0)& =X_0(x)\ge 0,Y(x,0)=Y_0(x)\ge 0,Z(x,0)=Z_0(x)\ge 0,x\in \mathrm{\Omega },\end{array}\end{array}$(7) ) when $d_{11}=1$, $d_{22}=0.05$, $d_{33}=0.01$, and $m_1=0.51\in (m_1^{\ast },{\hat{m}}_1)$: (a) shared resource, (b) IG prey, (c) IG predator. (a) $X(t)$, (b) $Y(t)$ and (c) $Z(t)$.

Whether the hunting cooperation in PDE model has the same role in that ODE model? We choose that the new values of all parameters are as follows: $r=0.7,d_1=0.12,{\alpha }_1=1.9,{\alpha }_2=2.5,m_1=0.2,m_2=0.1,b=10,{\beta }_1=2.2,{\alpha }_3=2.57,c=0.1,d_2=0.2,{\beta }_2=1.1,{\beta }_3=2,d_3=\mathrm{0.5.}$ When $d_{11}=10,d_{22}=20$ and $d_{33}=15$, the effect of hunting cooperation on the model (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }t}& =d_{11}\mathrm{\Delta }X+rX-a{X}^2-d_{1}X-{\alpha }_1(1-m_1)XY\\ & -\frac{{\alpha }_2(1-m_2)XZ}{1+b(1-m_2)X},x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }t}& =d_{22}\mathrm{\Delta }Y+{\beta }_1(1-m_1)XY-({\alpha }_3+cZ)YZ-d_{2}Y,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }t}& =d_{33}\mathrm{\Delta }Z+\frac{{\beta }_2(1-m_2)XZ}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)YZ-d_{3}Z,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }\nu }& =\frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }\nu }=0,x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ X(x,0)& =X_0(x)\ge 0,Y(x,0)=Y_0(x)\ge 0,Z(x,0)=Z_0(x)\ge 0,x\in \mathrm{\Omega },\end{array}\end{array}$(7) ) is shown in Figure 13. From Figure 13, we know that the model (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }t}& =d_{11}\mathrm{\Delta }X+rX-a{X}^2-d_{1}X-{\alpha }_1(1-m_1)XY\\ & -\frac{{\alpha }_2(1-m_2)XZ}{1+b(1-m_2)X},x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }t}& =d_{22}\mathrm{\Delta }Y+{\beta }_1(1-m_1)XY-({\alpha }_3+cZ)YZ-d_{2}Y,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }t}& =d_{33}\mathrm{\Delta }Z+\frac{{\beta }_2(1-m_2)XZ}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)YZ-d_{3}Z,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }\nu }& =\frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }\nu }=0,x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ X(x,0)& =X_0(x)\ge 0,Y(x,0)=Y_0(x)\ge 0,Z(x,0)=Z_0(x)\ge 0,x\in \mathrm{\Omega },\end{array}\end{array}$(7) ) becomes more stable when the effect of hunting cooperation increases. That is, under the case diffusion, the hunting cooperation of IG predator can make the system stale.

Figure 13. Effect of hunting cooperation on model (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }t}& =d_{11}\mathrm{\Delta }X+rX-a{X}^2-d_{1}X-{\alpha }_1(1-m_1)XY\\ & -\frac{{\alpha }_2(1-m_2)XZ}{1+b(1-m_2)X},x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }t}& =d_{22}\mathrm{\Delta }Y+{\beta }_1(1-m_1)XY-({\alpha }_3+cZ)YZ-d_{2}Y,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }t}& =d_{33}\mathrm{\Delta }Z+\frac{{\beta }_2(1-m_2)XZ}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)YZ-d_{3}Z,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }\nu }& =\frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }\nu }=0,x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ X(x,0)& =X_0(x)\ge 0,Y(x,0)=Y_0(x)\ge 0,Z(x,0)=Z_0(x)\ge 0,x\in \mathrm{\Omega },\end{array}\end{array}$(7) ) when $m_1=0.2,m_2=0.1$, $d_{11}=10$, $d_{22}=20$ and $d_{33}=15$: (a) c = 0.45; (b) c = 1; (c) c = 2.5.

5. Conclusion

In this paper, we have investigated an intraguild prey–predator model with prey refuge and hunting cooperation. By choosing prey refuge $m_1$ as the Hopf bifurcation parameter, we obtained the condition of the existence of the Hopf bifurcation, and showed the direction of Hopf bifurcation and stability of the bifurcating periodic solutions by using the normal form theory and the centre manifold theorem. According to a lot of numerical simulations, we found that only hunting cooperation c always makes model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) unstable, but model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ) becomes stable with the increase of prey refuge effect $m_1$. Moreover, our research also explains the reasons for prey refuge and hunting cooperation in real life from the perspective of Hopf bifurcation, which effectively verifies the reliability of the theoretical results. From an ecological point of view, hunting cooperation makes the ODE model with intraguild predation unstable, which increases the mortality rate of the three populations, and the possibility of extinction occurs in severe cases, while the refuge effect is within a certain range which makes the model stable. That is, maintaining the survival of various species and ensuring that the population does not become extinct. In addition, we study the diffusive model (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }t}& =d_{11}\mathrm{\Delta }X+rX-a{X}^2-d_{1}X-{\alpha }_1(1-m_1)XY\\ & -\frac{{\alpha }_2(1-m_2)XZ}{1+b(1-m_2)X},x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }t}& =d_{22}\mathrm{\Delta }Y+{\beta }_1(1-m_1)XY-({\alpha }_3+cZ)YZ-d_{2}Y,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }t}& =d_{33}\mathrm{\Delta }Z+\frac{{\beta }_2(1-m_2)XZ}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)YZ-d_{3}Z,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }\nu }& =\frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }\nu }=0,x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ X(x,0)& =X_0(x)\ge 0,Y(x,0)=Y_0(x)\ge 0,Z(x,0)=Z_0(x)\ge 0,x\in \mathrm{\Omega },\end{array}\end{array}$(7) ), which has the same positive equilibrium with the ODE model (9(9) $\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{J}\mathbf{(}\mathbf{U}\mathbf{)}.$(9) ). It is shown that diffusion-driven Turing instability occurs under certain conditions. Finally, we also discuss the existence and non-existence of non-constant positive steady states of the system (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }t}& =d_{11}\mathrm{\Delta }X+rX-a{X}^2-d_{1}X-{\alpha }_1(1-m_1)XY\\ & -\frac{{\alpha }_2(1-m_2)XZ}{1+b(1-m_2)X},x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }t}& =d_{22}\mathrm{\Delta }Y+{\beta }_1(1-m_1)XY-({\alpha }_3+cZ)YZ-d_{2}Y,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }t}& =d_{33}\mathrm{\Delta }Z+\frac{{\beta }_2(1-m_2)XZ}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ)YZ-d_{3}Z,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }\nu }& =\frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }\nu }=0,x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ X(x,0)& =X_0(x)\ge 0,Y(x,0)=Y_0(x)\ge 0,Z(x,0)=Z_0(x)\ge 0,x\in \mathrm{\Omega },\end{array}\end{array}$(7) ) by establishing positive upper and lower bounds by using Leray–Schauder degree theory. Biologically, the refuge effect causes Turing instability within a certain range, which shows the population inhomogeneously distributed within the spatial domain and exhibits various patterns. As the hunting cooperative effect increases, the model becomes more stable. If two factors are properly introduced together, then the system will become more stable, and the population is homogeneously distributed in space. When the shared resource or IG prey populations diffuse rapidly, they will break the original stable state and are spatially inhomogeneously distributed. While the three populations diffuse rapidly, they will not change the stable state, that is, the population is spatially inhomogeneously distributed.

From the viewpoint of human profit, the humans will capture some corresponding organisms to obtain certain economic benefits and the degree of economic impact on people often depends on the implementation of the harvest in nature. In addition, some behaviours of predator or prey species cannot be displayed immediately in ecology, that is, there is a time delay. To make the model more practical, we will consider time delay and linear harvesting into our model, $\{\begin{array}{rl}\frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }t}& =d_{11}\mathrm{\Delta }X+rX-a{X}^2-d_{1}X-{\alpha }_1(1-m_1)XY\\ & -\frac{{\alpha }_2(1-m_2)XZ}{1+b(1-m_2)X},x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }t}& =d_{22}\mathrm{\Delta }Y+{\beta }_1(1-m_1)XY-({\alpha }_3+cZ)YZ-d_{2}Y-q_{1}EY,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }t}& =d_{33}\mathrm{\Delta }Z+\frac{{\beta }_2(1-m_2)XZ}{1+b(1-m_2)X}+{\beta }_3({\alpha }_3+cZ(t-\tau ))Y(t-\tau )Z(t-\tau )\\ & -d_{3}Z-q_{2}EZ,x\in \mathrm{\Omega },t>0,\\ \frac{\mathrm{\partial }X(x,t)}{\mathrm{\partial }\nu }& =\frac{\mathrm{\partial }Y(x,t)}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }Z(x,t)}{\mathrm{\partial }\nu }=0,x\in \mathrm{\partial }\mathrm{\Omega },t>0,\\ X(x,0)& =X_0(x)\ge 0,Y(x,0)=Y_0(x)\ge 0,Z(x,0)=Z_0(x)\ge 0,x\in \mathrm{\Omega },\end{array}$where E is the harvesting effort, $q_1$ and $q_2$ are the catchability coefficient of IG prey and IG predator, respectively. τ is the time delay due to the gestation of the IG predator. We leave this work in the future.

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 Nos. 12161054, 11661050 and 11861044), the National Natural Science Foundation of Gansu Province (Grant No. 20JR10RA156), the Doctoral Foundation of Lanzhou University of Technology and the HongLiu First-class Disciplines Development Program of Lanzhou University of Technology.