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.


Introduction
The interaction between predator and prey has become increasingly complex due to longterm 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 CONTACT Xin-You Meng xymeng@lut.edu.cn;Yan Feng 1742734928@qq.comestablished the first predator-prey model with prey refuge: 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: where mx (m ∈ [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), a three-component model composed of single and double prey species is proposed in the reference [40]: ( 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] where N and P are the densities of the prey and the predator respectively.(N, P) is Holling-type I function, i.e. (N, P) = λN, here λ > 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) by replacing the constant attack rate λ with a density-dependent term: (N, P) = (λ + aP)N, ( 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: 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) 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.
The first general model describing IGP was exploited by Holt and Polis [22], which takes the following form: where R(t), N(t) and P(t) represent the densities of the shared resource, IG prey and IG predator, respectively.ρ 2 (R, N, P)R and ρ 3 (R, N, P)N are functional responses of the IG predator to the shared resource and IG prey, respectively, ρ 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ϕ(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 Hollingtype 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: 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.⊂ R N (N ≤ 3) is a fixed spatial bounded domain with smooth boundary ∂ and ν is the outward unit normal vector of boundary ∂ .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; α 1 , α 2 b and α 3 are maximum numbers of different prey that can be eaten by the respective predator per unit time; β 1 , β 2 and β 3 are a conversion rate of different prey captured by respective predator.m 1 ∈ [0, 1) and m 2 ∈ [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).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 U = (X, Y, Z) T , D = diag(d 11 , d 22 , d 33 ) and Then model (7) can be rewritten as 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) are discussed.In Section 3, the diffusion-driven Turing instability of model ( 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.

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).The ODE model is as follows:

The existence and local stability of equilibria
It is obvious that (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 and (H 2 ) : ; (4) the semi-trivial equilibrium in the absence of IG prey E 3 (X 3 , 0, Z 3 ) if assumption (H 3 ) holds, where and (H 3 ) : Table 1.The number of positive roots of Equation (10) The sign of coefficients Number of sign changes Number of positive roots (5) the co-exist equilibrium E * (X * , Y * , Z * ) is expressed as follows: and Z * determined by the following equation: where , 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) is shown in Table 1.
Based on the above analysis, it can be seen that the existence of the co-exist equilibrium of model ( 9) is very complicated.Therefore, we only discuss the second and fourth cases in this paper, that is, model ( 9) has only one co-exist equilibrium if the assumption (H 4 ) : Next, we will use the standard linearization method to discuss the local stability of each equilibrium.The Jacobian matrix of the model ( 9) at any non-negative equilibrium (X, Y, Z) is given by where Lemma 2.1: The trivial equilibrium E 0 (0, 0, 0) of model ( 9) is always unstable.

Proof:
The characteristic equation of model ( 9) at the trivial equilibrium E 0 is Thus we can get that the corresponding eigenvalues are So the trivial equilibrium E 0 is unstable.
Proof: From matrix (11), the characteristic equation of model (9) at equilibrium E 2 takes the form where Therefore, λ 1 = B 33 , λ 2 and λ 3 are the roots of the following equation: where Hence, the semi-trivial equilibrium E 2 is locally asymptotically stable if (H 5 ) is true.
Proof: From the Jacobian matrix (11), the characteristic equation of model ( 9) at the semitrivial equilibrium E 3 takes where Further, λ 1 = C 22 , λ 2 and λ 3 are the roots of following equation: 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 * (X * , Y * , Z * ).According to the matrix (11), the characteristic equation of model ( 9) at E * is where Lemma 2.5: All roots of Equation (12) Further, the co-existence equilibrium E * loses stability when G 5 G 6 − G 7 passes through 0. That is, the Hopf bifurcation occurs at the co-existence equilibrium E * 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, where Without loss of generality, let m * 1 be the positive root of Equation (13).Then Equation ( 12) has a negative root denoted λ 1 and a pair of complex conjugates roots which are denoted Substituting λ 2 = u(m 1 ) + iv(m 1 ) into Equation ( 12) and separating the real and imaginary parts gives Because u(m 1 ) = 0, from Equation ( 15) we have Then we substitute (16) in Equation ( 14) and get So Equation ( 17) has a root u(m * 1 ) if and only if )) = 0. Furthermore, since 8u 2 + 8uG 5 + 2(G 2  5 + G 6 ) = 0 has no real root, u = 0 is the only root.Thereby, when m 1 = m * 1 , the roots of ( 12) are λ 1 = −G 5 , λ 2,3 = ±i √ G 6 .We will prove that du dm 1 | m 1 =m * 1 = 0. Differentiating Equation ( 17) with respect to m 1 gives When 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) at E * are functions of m 1 and G i (m 1 )(i = 5, 6, 7) at m 1 = m * 1 satisfy the following conditions: Thus Hopf bifurcation occurs at the co-existence equilibrium E * when m 1 = m * 1 .

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) by using the centre manifold theory and normal form method.
For the sake of convenience, we translate the equilibrium where Let p 1 and p 2 be the eigenvectors corresponding to eigenvalues λ 1 = −G 5 and λ 2,3 = ±i √ G 6 ±iω of J U (U * ) at m 1 = m * 1 , respectively.We can obtain p 1 and p 2 , where , e 32 = 0, e 33 = 1.
Then we obtain where By the transformation (X, Y, Z) T = E(x, y, z) T , system (19) becomes where Thus the system (20) can be written as By the algorithm proposed by Hassard et al. [17], we get and Then, we can compute the following values: 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.

Linear stability and Turing instability
In this section, we consider the influence of the diffusion on the stability of the co-existence equilibrium E * 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 = μ 1 < μ 2 < μ 3 < • • • be the eigenvalues of the operator − on with the homogeneous Neumann boundary condition and E(μ i ) be the eigenspace corresponding to μ i in H 1 ( ).Let {φ ij : j = 1, 2, . . ., dim(E(μ i ))} be an orthonormal basis of E(μ i ), and Suppose that (H 4 ) − (H 7 ) hold, let L = D + J U (U * ).The linearization system of (8) at U * is U t = LU.For each i ≥ 1, X i is invariant under the operator L, and λ is an eigenvalue of L if and only if it is an eigenvalue of the matrix −μ i D + J U (U * ) for some i ≥ 1, in which there is an eigenvector in X i .The characteristic polynomial of −μ i D + J U (U * ) is given by where where In the following, we will prove that there exists a positive constant δ such that Let λ = μ i ρ, we have Furthermore

A-priori bounds for the positive steady state
The corresponding steady-state problem of ( 8) is In the following, the generic constants C 1 , C 2 , C * , C, 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 , α i , β 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).First, we need two results as follows.

Lemma 3.4 (Harnack Inequality [28]):
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( , k) and C( , k) such that, when d ii ≥ k, i = 1, 2, 3, the positive solution (X, Y, Z) to (26) Proof: First of all, we will prove that this estimates max Applying Lemma 3.3 to the first equation of ( 26), it yields X ≤ r a .Define A simple application of Lemma 3.3 to the third equation of ( 26) obtains max ¯ Z ≤ C 1 min ¯ Z for a positive constant C 1 = C 1 ( , k).This together with (29) Because of (28), by the regularity for elliptic equations we have that X, Y, Z belong to C 1+α ( ¯ ), and C 1+α ( ¯ ) norms of them depend on known parameters r, a, b, c, d i , α i , β i , m j , (i = 1, 2, 3, j = 1, 2).Using the regularity of elliptic equations again, the estimate ( 27) follows.

Theorem 3.6 (Lower Bound):
Let k be a fixed positive number.Then there exist a positive constant C = C( , k) such that, when d ii ≥ k, i = 1, 2, 3, the positive solution (X, Y, Z) to (26) Proof: We suppose that (30) does not hold.Then there exists a sequence {(d ) and the corresponding positive solutions (X j , Y j , Z j ) to ( 26) such that min ¯ X j → 0, or min ¯ Y j → 0, or min ¯ Z j → 0, as j → ∞.
According to (27), this implies that max As The standard regularity theorem for elliptic equations yields that there exists a subsequence of {(X j , Y j , Z j )} ∞ j=1 denoted by {(X j , Y j , Z j )} ∞ j=1 , and three non-negative functions 3 as j → ∞.By (31), we find that (32), then we have Next, we will consider three cases as follows.
Case (i).X ≡ 0. Since X j ≡ X, as j → ∞, then β 1 (1 − m 1 )X j − (α 3 + cZ j )Z j − d 2 < 0 on ¯ , for all j 1. Integrating the differential equation for Y j over by parts, we have which is a contradiction to the second equation of (32).Case (ii).Y ≡ 0. Since Y j ≡ Y, as j → ∞, then Z j < 0 on ¯ , for all j 1, and which is a contradiction to the third equation of (32).Case (iii).Z ≡ 0, X ≡ 0, Y ≡ 0 on ¯ , then the Hopf lemma gives X > 0, Y > 0.
(a) d11 = ∞.Under the circumstances, X satisfies When (1−m 1 )k 1 , we get Y ≡ 0, which is a contradiction. If By the third equation of ( 32), there exists xj ∈ such that i.e.Y j (x j ) = . Furthermore, we assume that xj → x as j → ∞.So Equation (34) implies . Then from the first equation of ( 33), we get which is a contradiction to Lemma 4.5 [7].
By Lemma 4.4 [7], (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].
Notice that Y ≡ k 2 if d22 = ∞.The third equation of (37) means which makes X to be a constant k 1 , a contradiction to Lemma 4.5 [7].Therefore, the proof of this theorem is completed.

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) 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.Proof: Assume that (X, Y, Z) is a positive solution to (26) and denote Multiplying the first equation of ( 26) and then integrating on by parts yields By the similar method, we can get From the above results and the ε-Young inequality, it follows that where ε is a small enough positive number.Using the Poincáre inequality, we get Hence, this proof is completed.

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), that is the stationary pattern can exist for the system (7) with self-diffusion.
Let X be the Hilbert space as in Section 3.1, and define Next, ( 26) can be written as follows: Thus U is a positive solution of (38) if and only if it meets where (I − ) −1 is the inverse of I − with the homogeneous Neumann boundary condition.As L(•) is a compact perturbation of the identity operator, for any B = B(C), the Leray-Schauder degree deg(L(•), 0, B) is well-defined if L(U) = 0 on ∂B.Furthermore, simple computation yields , where κ is the sum of algebraic multiplicity of negative eigenvalues of D U L(U * ).
We refer to the decomposition (23) in the following discussion of the eigenvalues of D U L(U * ).First of all, for each integer i ≥ 1, 1 ≤ j ≤ dimE(μ i ), X ij is invariant under D U L(U * ), and λ is an eigenvalue of D U L(U * ) if and only if λ(1 + μ i ) is an eigenvalue of the matrix μ i I − D −1 J U (U * ).Thus D U L(U * ) is invertible if and only if this matrix is non-singular.Let To calculate the index (L(•), U * ), we cite the following lemma.

Lemma 3.8 ([7]): Suppose that the matrix μ
The direct computation gives where Furthermore, it obtains By determining the range of μ such that M(μ) < 0, we have the existence of non-constant steady states of (26).
U is a non-constant positive solution of (26) if and only if it is a solution of (42) with t = 1.
It is obvious that U * is the unique constant positive solution of ( 42) for any 0 ≤ t ≤ 1.For any t ∈ [0, 1], U is a non-constant positive solution of (42) if and only if Obviously, L(1; U) = L(U).According to Theorem 3.7, L(0; U) = 0 has only the positive solution U * in W + .Then, we obtain where D = diag( d11 , d22 , d33 ).From (39), we know that By (40) and (41), it follows that Hence, zero is not an eigenvalue of the matrix μ i I − D −1 J U (U * ) for all i ≥ 1, and which contradict (44).Thus we finish the proof of this theorem.

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) are given according to the meaning of parameters.What's more, because model ( 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, α 1 = 1.9, 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) are assumed to be X(0 First of all, we mainly illustrate the stability and Hopf bifurcation for non-spatial model (9) by using the ode45 function of Matlab software.According to Theorem 2.6, we obtain that model ( 9) has the positive equilibrium E * = (0.9947, 0.0784, 0.1764) and is locally asymptotically stable (see Figure 2).Subsequently, we choose m 1 of model ( 9) as a bifurcation parameter.Based on the parameter values given above, we can get m * 1 = 0.59.Furthermore, when m 1 = 0.7 > m * 1 , E * is locally asymptotically stable (see Figure 3a), but E * is unstable when m 1 = 0.54 < m * 1 (see Figure 3b).That is, Hopf bifurcation takes place when m 1 = m * 1 .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.
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) at the positive equilibrium   E * 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.Based on Theorem 2.6, the dynamic behaviour of the system at the positive equilibrium E * without the effect of hunting cooperation is seen in Figure 5. Further, it can be concluded that the stability of model ( 9) at E * is better with the increase of prey refuge.
In addition, according to Lemma 2.5, the dynamical behaviour of model ( 9) at the positive equilibrium E * 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.
While the effect of prey refuge is considered, i.e. m 1 = 0.7, m 2 = 0.4, the model (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.
Finally, we will check the spatial dynamics of model ( 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π ] for the convenience of simulation.We keep all values of all parameters with the same of the above and choose the values of diffusions    which is shown in Figure 9, but is unstable when m 1 = 0.5 ∈ (0, m * 1 ), which is shown in Figure 10.That is, it also can be concluded that the stability of model ( 9) at E * is better with the increase of prey refuge.
To check Turing instability, we choose that the values are r = 0.    7) is shown in Figure 13.From Figure 13, we know that the model ( 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.

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) unstable, but model ( 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), which has the same positive equilibrium with the ODE model (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) 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, ∂X(x, t) ∂ν = ∂Y(x, t) ∂ν = ∂Z(x, t) ∂ν = 0, x ∈ ∂ , t > 0, X(x, 0) = X 0 (x) ≥ 0, Y(x, 0) = Y 0 (x) ≥ 0, Z(x, 0) = Z 0 (x) ≥ 0, x ∈ , 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.
investigated a prey-predator model with hunting cooperation in a Holling type I function: N, P)P − mP,
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) 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.