Dynamics in a Predator–Prey Model with Cooperative Hunting and Allee Effect

This paper deals with a diffusive predator–prey model with two delays. First, we consider the local bifurcation and global dynamical behavior of the kinetic system, which is a predator–prey model with cooperative hunting and Allee effect. For the model with weak cooperation, we prove the existence of limit cycle, and a loop of heteroclinic orbits connecting two equilibria at a threshold of conversion rate p=p#, by investigating stable and unstable manifolds of saddles. When p>p#, both species go extinct, and when p<p#, there is a separatrix. The species with initial population above the separatrix finally become extinct, and the species with initial population below it can be coexisting, oscillating sustainably, or surviving of the prey only. In the case with strong cooperation, we exhibit the complex dynamics of system, including limit cycle, loop of heteroclinic orbits among three equilibria, and homoclinic cycle with the aid of theoretical analysis or numerical simulation. There may be three stable states coexisting: extinction state, coexistence or sustained oscillation, and the survival of the prey only, and the attraction basin of each state is obtained in the phase plane. Moreover, we find diffusion may induce Turing instability and Turing–Hopf bifurcation, leaving the system with spatially inhomogeneous distribution of the species, coexistence of two different spatial-temporal oscillations. Finally, we consider Hopf and double Hopf bifurcations of the diffusive system induced by two delays: mature delay of the prey and gestation delay of the predator. Normal form analysis indicates that two spatially homogeneous periodic oscillations may coexist by increasing both delays.


Introduction
A predator-prey model can be described by the following equations [1]: where u and v stand for the densities of prey and predator, respectively. f (u) represents the per capita prey growth rate in the absence of predators, and G(u, v) is the functional response characterizing predation. d describes the rate of biomass conversion from predation, and m is the death rate of predator. Many kinds of functional response have been proposed [2]. Among them, Holling type I, II, and III [3][4][5][6] are discussed widely, which are prey-dependent functional response. Cooperative behavior within a species is very common in nature, of which cooperative hunting is the most widely distributed form in animals. In mammals, the most famous examples include wolves [7], African wild dogs [8], lions [9], and chimpanzees [10], who cooperate for hunting preys. Cooperative behaviors are also widespread in other living where r 1 is the per capita intrinsic growth rate of the prey, K 1 is the environmental carrying capacity for the prey, a 1 (a 1 < K 1 ) is the Allee threshold of the prey population, b 1 is the attack rate per predator and prey, c 1 measures the degree of cooperation of predator, p 1 is the prey conversion to predator, and m 1 is the per capita death rate of predator. They discussed the existence and stability of the equilibria and the occurrence of Hopf bifurcation. They also devised a best strategy for culling the predator by maximizing the prey population and minimizing the predators along with the costs associated with the control. In fact, the preys and predators distribute inhomogeneously in different locations, and the spatial diffusion plays an important part in the process of population evolution.
It has been widely accepted that time delays have very complex effect on the dynamics of a system, for example, some delays can destroy the stability of equilibria and induce various oscillations and periodic solutions (see, for example, in [35][36][37][38] and references cited therein). There are time delays in almost every process of population interaction, so it is more realistic to introduce time delays when we model the interaction of predator and prey. For example, after consuming the prey, the corresponding biomass cannot be added to the predator population immediately. Thus, the reproduction of predator after predating the prey will be mediated by some discrete time lag (say τ 2 ) required for gestation of predator, which can be regarded as gestation delay. Therefore, the growth rate of predator depends on the number of individuals present at time t − τ 2 . There is an extensive literature about the studies of the dynamics of predator-prey models with the effect of gestation delay of the predator [35][36][37][38]. There is a time lag (say τ 1 ) in the growth to maturity of the prey, and such a maturation delay of the prey has also been considered [39][40][41][42][43]. However, to the best of our knowledge, there is very little literatures on delayed predator-prey system with Allee effect [44].
Taking into account the diffusion and delays in system (2), we consider the following diffusive system with two delays: Here, τ 1 is the time delay due to the maturation of the prey, τ 2 is time delay due to gestation of the predator, and we define τ = max{τ 1 , τ 2 }. d 1 , d 2 > 0 are the diffusion coefficients characterizing the rates of the spatial dispersion of the prey and predator population, respectively.
Systems with multiple delays have attracted much attention [45][46][47][48][49][50][51]. Generally, delay may induce Hopf bifurcation, and if Hopf bifurcation curves intersect, double Hopf bifurcation may arise. To figure out the effect of delay on the dynamics of systems, Hopf bifurcation and double Hopf bifurcation induced by delay have been investigated [52][53][54][55][56]. However, in most of literatures we mentioned above, the system is reduced into a system with one delay. In fact, systems with multiple delays conform to reality better than one single delay. Recently, the research on the dynamics and bifurcation analysis of system with two simultaneously varying delays are of great interest to scholars. In [49,57], double Hopf bifurcation induced by two delays is studied by different methods. Through the analysis of double Hopf bifurcation, we can classify the topological structures of various bifurcating solutions. By the classification, the dynamics in the neighborhood of the double Hopf bifurcation point in system can be obtained completely.
For the convenience of discussion and reducing the number of parameters, we nondimensionalize the model by letting Dropping the hats for simplicity of notations, we obtain the simplified dimensionless system where r = K 1 r 1 , a = a 1 K 1 , c = c 1 b 2 1 , and p = b 1 p 1 K 1 m . Recalling the definition of Allee effect, one always has 0 < a < 1. In fact, for fixed K 1 , r 1 , a 1 , b 1 , and m, we have c and p directly proportional to c 1 and p 1 , respectively. Thus, we refer to c and p as the degree of cooperative hunting of the predator and the conversion rate from the prey to the predator, respectively, in the context. By the explanation above, we make the following assumption always: (H) r, a, c, m, p and c are all positive, and 0 < a < 1.
The Allee effect can drive the extinction of low-density species, some of which are preyed by predators that engage in cooperative hunting. This paper investigate the dynamics of system (4) with cooperative hunting and Allee effect from the point of view bifurcation analysis, which can provide useful insights into the intersection of the predator and prey, and help us figure out how to prevent such a prey population from extinction. We find out that the number of constant steady states of (4) is highly dependent on the strength of cooperation in predator c. When c < 1 r(1−a) , there is one or no positive steady states, while when c > 1 r(1−a) , there is no, one or two positive steady states. We say the cooperation in the former case is weak, and that in the latter is strong. First of all, we figure out the global dynamical behavior of the kinetic system through the analysis of local bifurcation and the stable and unstable manifolds, taking p as a bifurcation parameter. The existence of some connecting orbits is obtained, thus, a partition of the phase space is given, i.e., the attraction basin of each stability state is obtained. For the model with weak cooperation, we give the condition of the stability of the positive steady state and the occurrence of Hopf bifurcation. Through detailed analysis on the stable and unstable manifolds of saddles, we prove theoretically there is a loop of heteroclinic orbits when p = p # . When p > p # , the extinction steady state E 0 (0, 0) is globally asymptotically stable. When p < p # , there is a separatrix, which is the stable manifold of E a . The orbits above the separatrix converge to E 0 , and the orbits below it may converge to E 1 , E * or a limit cycle. Thus, there are three types of bistability: bistability between two boundary equilibria, between a boundary equilibrium and a interior equilibrium, and between a boundary equilibrium and a limit cycle. For the model with strong cooperation, it turns out there are many complex phenomena, including limit cycle, homoclinic cycle, and a loop of heteroclinic orbits between two or among three equilibria. The model may display tristable dynamics, including extinction state, coexistence or sustained oscillation, and the survival of the prey only.
Second, we discuss the diffusive system (4) with no delays. Taking diffusive coefficients as parameters, we give the condition of Turing instability and Turing-Hopf bifurcation in both weak and strong cooperation cases. We discuss the effect of diffusion on the dynamics, and illustrate complex spatio-temporal dynamical behavior, including the existence of spatially inhomogeneous steady state, coexistence of two spatially inhomogeneous periodic solutions.
Third, taking two time delays as bifurcation parameters, we study the effect of two simultaneous varying delays on the dynamics of system (4). After some detailed analysis on the characteristic equation with three exponential terms, the stability switching curves are obtained, on which there is a pair of purely imaginary eigenvalues. An intersection of the stability switching curves may be a double Hopf bifurcation point. We derive the normal form on the center manifold near the double Hopf bifurcation point, and give the unfoldings near the critical points, which indicate that two spatially homogeneous periodic oscillations may coexist by increasing both delays.
The scheme of this paper is arranged as follows. In Section 2, we investigate the existence of constant steady states for system (4). In Section 3, we analyze the local bifurcation and the global dynamics of the kinetic system (7) with weak cooperation and strong cooperation, respectively. In Section 4, we consider the effect of diffusion on system (4) through investigating Turing instability and Turing-Hopf bifurcation. In Section 5, we consider the diffusive system with two delays, and Hopf bifurcation and double Hopf bifurcation induced by two delays are analyzed on the center manifold via normal form approach. In Section 6, we carry out some numerical simulations to support our theoretical analysis and illustrate complex dynamical phenomenon.

The Existence of Constant Steady States
First, we consider the existence of steady states of model (4). As (4) is two-dimensional, it is hard to consider non-constant steady states. Thus, we only discuss the existence of constant steady states.
System (4) has three boundary equilibria: (i) E 0 (0, 0), which means the extinction of both the species; (ii) E a (a, 0), which is induced by Allee effect; and (iii) E 1 (1, 0), which means the extinction of the predator and the survival of the prey, achieving its carrying capacity. Now, we discuss the existence of interior equilibria. Obviously, any interior equilibrium has components u > 0 and v > 0 satisfying i.e., In fact, r(1 − u)(u − a) = (1 + cv)v is an ellipse, sitting in the first and forth quadrant, and intersecting the horizontal axis at a and 1. pu(1 + cv) = 1 intersects the horizontal axis at 1 p , and is decreasing and concave in (0, 1 p ). Denote the u-nullcline r(1 − u)(u − a) = (1 + cv)v as v = f (u), and the v-nullcline pu(1 + cv) = 1 as v = g(u). The existence of interior equilibria may have the following different cases. Proposition 1. (i) When p ≥ 1 a , system (4) has no interior equilibria. (ii) When 1 < p < 1 a , system (4) has a unique positive constant equilibrium E * (u * , v * ) with a < u * < 1 p . (iii) Suppose that p < 1. If c < 1 r(1−a) , system (4) has no interior equilibria. If c > 1 r(1−a) , there exists a p SN < 1, such that system (4) has (a) two interior equilibria E * and E * R when p SN < p < 1; (b) one interior equilibrium E * * R when p = p SN ; and (c) no interior equilibria when p < p SN .
Proof. In fact, when a < 1 p < 1, two nullclines have two intersections: E * in the first quadrant and E * R in the forth quadrant. Thus, there is a unique interior equilibrium E * (u * , v * ) with u * and v * satisfying (5). When 1 p decreases to a, E * collides with E a , and E * R is still in the forth quadrant. When 1 p < a, E * moves into the forth quadrant, leaving no interior equilibria. When 1 p increases to 1, v-nullcline and u-nullcline intersect at E 1 (1, 0), with the slopes of tangents at E 1 (1, 0) being − 1 c and −r(1 − a), respectively. If − 1 c < −r(1 − a), E * collides with E 1 when p = 1, E * R is in the forth quadrant, and there is no interior equilibria (see Figure 1a). Moreover, when p < 1, there is no interior equilibria. If − 1 c > −r(1 − a), E * R (sitting in the forth quadrant when 1 < p < 1 a ) collides with E 1 when p = 1, and there is a unique interior equilibrium E * in the first quadrant (see Figure 1b). When p < 1 and close to 1, E * R moves into the first quadrant, and there are two interior equilibria E * and E * R . If p < 1 goes on to decrease, there exists a p SN , such that the u-nullcline v = f (u) and v-nullcline v = g(u) tangent at E * * R , where E * and E * R collide, and there is a unique interior equilibrium E * * R . When p < p SN , there is no interior equilibria.

Remark 1.
In the absence of cooperative hunting within the predator, i.e., c = 0, system (4) always has a unique interior equilibrium ). The predator and prey coexist if and only if 1 < p < 1 a , and the predator is distinct if p < 1. However, from proposition 1, the predator and prey still coexist when p SN < p < 1 with strong cooperation, which means that strong cooperation is beneficial for the survival of predator.
Proposition 1 tells us that the strength of cooperation c = 1 r(1−a) is a threshold for the system (4). To be specific, if c < 1 r(1−a) , system (4) has at most one coexisting equilibrium, while if c > 1 r(1−a) , system (4) may have two coexisting equilibria. Thus, in the following discussion, when c < 1 r(1−a) , we say the cooperation among predators is weak; when c > 1 r(1−a) , the cooperation is strong.

Bifurcation and Global Dynamics of Kinetic System
Before the discussion of the stability of equilibria of system (4), we should first figure out the stability of the equilibria of the kinetic system We investigate the local bifurcation and global dynamics when the cooperation is weak and strong, i.e., when c < 1 r(1−a) and c > 1 r(1−a) .

Bifurcation and Global Dynamics of Kinetic System with Weak Cooperative Hunting
From proposition 1, if c < 1 r(1−a) , system (4) has a unique interior equilibrium. In this section, taking p as bifurcation parameter, we first investigate the Hopf bifurcation near the unique interior equilibrium. Through studying the stable manifold and unstable manifold of saddles, we prove the existence of loop of heteroclinic orbits, and get the global dynamics of system (7).
3.1.1. Stability of All Equilibria and Hopf Bifurcation at E * The Jacobian matrices of the function on the right-hand of system (7) around E 0 , E a , and E 1 are , respectively, from which we can easily get the local stability of the boundary equilibria.
Lemma 1. For system (7), is an unstable node if p > 1 a , and it is a saddle if p < 1 a ; and (iii) E 1 (1, 0) is a stable node if p < 1, and it is a saddle if p > 1.
For the interior equilibrium E * (u * , v * ), the corresponding Jacobian matrix is and thus Thus E * may be node or focus, and its stability depends on the trace Thus, we have the following conclusions on the stability of E * . Theorem 1. If c < 1 r(1−a) , then there exists a unique p H ∈ (1, 2 a+1 ) such that the unique interior equilibrium E * is locally asymptotically stable when 1 < p < p H , and unstable when p H < p < 1 a . Moreover, system (7) undergoes a Hopf bifurcation at E * when p = p H .

Remark 2.
When there is no cooperative hunting in predator, i.e., c = 0, the v-nullcline is u = 1 p , which is vertical. trJ E * | c=0 = − r p ( 2 p − a − 1), detJ E * | c=0 = mv * . The system undergoes a Hopf bifurcation near E * when 1 p = a+1 2 , which is the top of the u-nullcline v = r(1 − u)(u − a). In fact, it has been discussed widely that the stability of positive equilibrium can be stated graphically by the u-nullcline when the v-nullcline is vertical: (u * , v * ) is unstable if the v-nullclines intersect to the left of a local maximum of the u-nullcline, and stable if they intersect to the right [29,[58][59][60], and thus, Hopf bifurcation occurs at the "top of the hump" of u-nullcline. However, for system (7) in this paper, the v-nullcline is not vertical, and we can verify that the intersection of the u-nullcline and the v-nullcline corresponding to the Hopf bifurcation point p = p H is on the right of the "top of the hump" of u-nullcline. In fact, the "top of the hump" of u-nullcline v = f (u) achieves at u = a+1 2 , at which the corresponding value of p is denoted by p top . When p = p top , E * can be also proved to be unstable. Comparing the results in systems with and without cooperation, we can conclude that cooperative hunting is more likely to bring instability into E * .
In order to determine the direction of Hopf bifurcation and the stability of the Hopf bifurcating periodic solution, we need to calculate the normal form near the Hopf bifurcation point. Through direct calculation following the steps in [61], we can get the following truncated normal form: Recalling that α (p H ) > 0, the direction of Hopf bifurcation and stability of Hopf bifurcating periodic solution are determined by the first Lyapunov coefficient a(p H ), and we have the following conclusion. Theorem 2. If c < 1 r(1−a) , system (7) undergoes a Hopf bifurcation at E * when p = p H . (i) If a(p H ) < 0, the bifurcating periodic solution is orbitally asymptotically stable, and it is bifurcating from E * as p increases and passes p H . (ii) If a(p H ) > 0, the bifurcating periodic solution is unstable, and it is bifurcating from E * as p decreases and passes p H .

The Global Dynamics of Kinetic System with Weak Cooperative Hunting
From the previous discussion, we have known the stability of all equilibria when c < 1 r(1−a) . Motivated by the work of [29,30], we investigate the global dynamics of system (7). To this end, we should first study the properties of stable manifold of E a and unstable manifold of E 1 , denoted by Γ s p and Γ u p , respectively. For the region bounded by Γ s p and Γ u p , we prove the existence or nonexistence of periodic orbits. We leave them in Appendix A, and only list main conclusions on the existence of a loop of heteroclinic orbit and global dynamics in the following.
From the monotonicity of S and U given in Proposition A1, we have the following conclusion. Proposition 2. If c < 1 r(1−a) , then there exists a unique p # ∈ (1, 1 a ), such that Γ s p # = Γ u p # , forming a heteroclinic orbit from E 1 to E a .
In fact, there is another heteroclinic orbit from E a to E 1 , which is formed by the unstable manifold of E a and the stable manifold of E 1 on the u−axis. Thus, there is a loop of heteroclinic orbits from E 1 to E a , and then back to E 1 .
a , E 0 (0, 0) is globally asymptotically stable (see Figure 3a). (ii) There is an ε 1 > 0 such that if 1 a − ε 1 < p < 1 a , Γ s p connects E * to E a , and the extinction equilibrium E 0 (0, 0) is globally asymptotically stable (see Figure 3b). (iii) There is an ε 2 > 0 such that if 1 < p < 1 + ε 2 , then Γ u p connects E 1 to E * . The orbits through any point above Γ s p converge to E 0 , and the orbits through any point below Γ s p converge to E * (see Figure 3c). (iv) If p < 1, the orbits through any point above Γ s p converge to E 0 , and the orbits through any point below Γ s p converge to E 1 (see Figure 3d). We leave the proof in Appendix A.2. Theorem 3 provides a global description of the dynamical behaviour of system (7) for p ∈ R + \(1 + ε 2 , 1 a − ε 1 ) = R + \I. When p ∈ I, if every periodic orbit of system (7) is orbitally stable, thus there can be at most one such orbit. From Proposition A5, we have the following conclusion.
Theorem 4. Suppose that c < 1 r(1−a) , the first Lyapunov coefficient a(p H ) < 0, and every periodic orbit of system (7) is orbitally stable. Then, p H < p # . (i) If 1 < p < p H , Γ u p connects E 1 to E * . The orbits through any point above Γ s p converge to E 0 , and the orbits through any point below Γ s p converge to E * . (ii) If p H < p < p # , E * is a repellor, and there is a unique limit cycle under Γ s p . The orbits through any point below Γ s p converge to the limit cycle. (see Figure 4a,b).
(iii) If p = p # , Γ s p = Γ u p , there are two heteroclinic orbits forming a loop of heteroclinic orbits from E 1 to E a and back to E 1 (see Figure 4c). The orbits through any point exterior to the loop converge to E 0 , and the orbits through any point interior to the cycle converge to the loop. (iv) If p # < p < 1 a , Γ s p connects E * to E a , and the extinction equilibrium E 0 (0, 0) is globally asymptotically stable (see Figure 4d). Figure 4. Phase portrait of (7) for c < 1 r(1−a) when (a) p H < p < p # and close to p H ; (b) p H < p < p # and close to p # ; (c) p = p # ; (d) p # < p < 1 a .

Remark 3.
We leave the proof in Appendix A.3. From the former conclusion on the global dynamics of model (7) with weak cooperation, we find that no matter how we choose the value of p, both of the two species are extinct if the ratio of predator to prey is high. Now, we discuss the population behavior when the ratio of predator to prey is low. When p < 1, because of the low conversion rate, the predator is extinct, while the prey reach the capacity of the environment. When 1 < p < p H , the predator and prey coexist, and tend to a constant value. When p H < p < p # , the predator and prey coexist, but oscillate sustainably. The amplitude of oscillation increases as p increases, and the minimum of the predator population is close to zero when p is close to p # , which will increase of risk on the extinction of predator. When the value of p is too high, high growth rate of prey leads to over hunting, and finally, both of the two species are extinct.

Remark 4.
Similar to the case in the model with Allee effect, there is a unique positive equilibrium in the model without Allee effect, and the varying p makes E * unstable and induces Hopf bifurcation. However, there is an additional saddle point (a, 0) in the model with Allee effect, and the stable manifold of the saddle point separates the positive quadrant into two positive invariant regions. If the initial populations are above the stable manifold of (a, 0), extinction happens for both species in the model with Allee effect, while coexistence happens for the model without Allee effect. This extinction is mediated by too high predator density that drives the prey density below the Allee threshold and both species become extinct.

Bifurcation and Global Dynamics of Kinetic System with Strong Cooperative Hunting
By strong cooperation we mean c > 1 r(1−a) . In this section, we first analyze the stability of interior equilibrium, then we prove the existence of Hopf bifurcation and give the condition of the existence and nonexistence of loop of heteroclinic orbits. We also exhibit the complex dynamics of model (7), such as limit cycle, loop of heteroclinic orbits, and homoclinic cycle.

Stability and Hopf Bifurcation of Interior Equilibria
Now, we consider the stability of interior equilibria. For E * or E * R , the corresponding Jacobian matrix is and thus Figure 5). The component u * of E * is always satisfying Figure 1). It follows that Thus, E * may be node or focus, and its stability depends on the trace Similar as the proof in Theorem 1, solving for trJ E * = 0, when a+1 Figure 2). Recalling that ∂u * ∂p < 0, corresponding to u * H , we obtain a unique p H (p SN < p H < 2 a+1 ), such that trJ E * > 0 when p ∈ (p H , 1 a ), and trJ E * < 0 when p ∈ (p SN , p H ). Let µ = α(p) ± iω(p) be the roots of µ 2 − trJ E * µ + detJ E * = 0 when p is near p H . We can prove α (p) | p=p H > 0, in a similar way as in Theorem 1. Therefore, we can get the following theorem.
Theorem 5. When c > 1 r(1−a) , there exists a unique p H ∈ (p SN , 2 a+1 ) such that E * is locally asymptotically stable if p SN < p < p H , and unstable if p H < p < 1 a . Moreover, system (7) undergoes a Hopf bifurcation at E * when p = p H .

The Global Dynamics of Kinetic System with Strong Cooperative Hunting
When c > 1 r(1−a) , system (7) may have one or two interior equilibria, which brings a few complications to the dynamics of model (7). There are three different cases of global dynamics when c takes different sizes. In this section, we exhibit the dynamics combining theoretical analysis and numerical simulations.
First, we consider the properties of the stable manifold of E a , unstable manifold of E 1 , as well as stable and unstable manifold of E * R , which are in Appendix B. Then, we can determine the global dynamics of model (7) when p is chosen in the following ranges. As in Theorem 3, we can prove the following conclusion.
a , E 0 (0, 0) is globally asymptotically stable (see Figure 6a). (ii) If p < 1 a and near 1 a , Γ s p (E a ) connects E * to E a , and the extinction equilibrium E 0 (0, 0) is globally asymptotically stable (see Figure 6b). (iii) If p < p SN , the orbits through any point above Γ s p (E a ) converge to E 0 , and the orbits through any point below Γ s p (E a ) converge to E 1 (see Figure 6c).
The remaining question is how are the dynamics of system (7) for p chosen in the rest ranges, i.e., p SN < p < 1 a . As larger c means less steep of the v-nullcline, the critical points (for example, p H , p # ) may appear in different ranges, and there are the following three different cases.
Case 1 Choose suitable c such that c > 1 r(1−a) , U 1 (1 + ) < S a (1 + ), and p H > 1. From Theorem A1, there is a loop of heteroclinic orbits when p = p # . We have the following conclusion.
Theorem 7. Suppose that c > 1 r(1−a) , U 1 (1 + ) < S a (1 + ), and p H > 1. Assume that the first Lyapunov coefficient a(p H ) < 0, and every periodic orbit of system (7) is orbitally stable, then The orbits through any point above Γ s p (E a ) converge to E 0 , and the orbits through any point inside the stable manifold Γ s p (E * R ) converge to E * . The orbits through any point below Γ s p (E a ) and exterior to Γ s p (E * R ) converge to E 1 (see Figure 7a).
The orbits through any point above Γ s p (E a ) converge to E 0 , and the orbits through any point below Γ s p (E a ) converge to E * . (iii) When p H < p < p # , E * is a repellor, and there is a unique limit cycle under Γ s p (E a ). The orbits through any point below Γ s p (E a ) converge to the limit cycle (see Figure 7b,c).
, and there are two heteroclinic orbits forming a loop of heteroclinic orbits between E 1 and E a (see Figure 7d). The orbits through any point exterior to the cycle converge to E 0 , and the orbits through any point interior to the cycle converge to the cycle.
(v) If p # < p < 1 a , Γ s p connects E * to E a , and the extinction equilibrium E 0 (0, 0) is globally asymptotically stable (phase portrait is similar as Figure 6b).
As E * is stable when p < p H , then E * is stable for all p SN < p < 1. If there is a periodic orbit below Γ s p (E a ), it must encircle E * , and it is unstable from inside, which contradicts with our assumption of the orbital stability of periodic orbit. It means that there is no periodic orbit under Γ s p (E a ) for all p < 1. From Proposition A9, E * is the ω−limit set of Ω 22 , and Γ s p (E * R ) must be above Figure 7a). Thus, the orbits through any point inside the stable manifold Γ s Figure 7. Phase portrait of system (7) for case 1 when (a) p > p SN and close to p SN , (b) p > p H and close to p H , (c) p H < p < p # and close to p # , and (d) p = p # .
The orbits through any point above Γ s p (E a ) converge to E 0 , and the orbits through any point inside the stable manifold Γ s p (E * R ) converge to E * . The orbits through any point below Γ s p (E a ) and exterior to Γ s p (E * R ) converge to E 1 (phase portrait is similar as Figure 7a). (ii) When p H < p < 1, there is a unique limit cycle inside the stable manifold Γ s p (E * R ) of E * R (see Figure 8a). The orbits through any point interior to Γ s p (E * R ) converge to the limit cycle. (iii) When 1 < p < p # , there is a unique limit cycle under Γ s p . The orbits through any point below Γ s p converge to the limit cycle (see Figure 8b).
, there are two heteroclinic orbits forming a loop of heteroclinic orbits from E 1 to E a and back to E 1 (see Figure 8c). The orbits through any point exterior to the cycle converge to E 0 , and the orbits through any point interior to the cycle converge to the cycle.
(v) If p # < p < 1 a , Γ s p connects E * to E a , and the extinction equilibrium E 0 (0, 0) is globally asymptotically stable (phase portrait is similar as Figure 6b).

Proof.
We can prove (i) similar as in Theorem 7. (ii) As a(p H ) < 0 and α (p H ) > 0, there is a stable periodic orbit bifurcating from Hopf bifurcation when p > p H . As in the proof of Theorem 7, we can prove that Γ s p (E * R ) must be above Γ u p (E * R ). From Proposition A9, the unique limit cycle is the ω−limit set of Ω 22 . The proof of (iii)-(v) is similar as that of Theorem 4.  (7) for case 2 when (a) p H < p < 1 and close to p H ; (b) p H < p < p # and close to p # ; It is easy to prove that when 1 < p < 1 a , E 0 is globally asymptotically stable (see Figure 9e).

Remark 5. In this case, the location of
We can observe the following dynamics by numerical simulations.
(i) Loop of heteroclinic orbits. The downward branch of the unstable manifold of E * R connects E * R to E 1 , and the upward connects E * R to E a , which collides with the the stable manifold of E a . The upward and downward unstable manifold of E * R , together with the unstable manifold of E a on the u− axis, forms a loop of heteroclinic orbits among E * R , E a and E 1 when p = p # (see Figure 9d). (ii) Homoclinic cycle. The upward unstable manifold and the left stable manifold of E * R collide, which forms a homoclinic cycle. Denote the parameter p as p = p hom (see Figure 9c). (iii) Limit cycle induced by Hopf bifurcation (see Figure 9a,b).

Remark 6.
There are two positive equilibria E * and E * R in both of the model with or without Allee effect, where E * R is always a saddle. In the model without Allee effect, the stable manifold of E * R separates the positive quadrant into two positive invariant regions. Two species coexist if the initial populations are above the stable manifold. However, there is an additional stable manifold of saddle point (a, 0) in the model with Allee effect, and two stable manifolds of the two saddle points separate the positive quadrant into three positive invariant regions. Two species only coexist between the two stable manifolds, and both species go extinct if the initial populations is above the stable manifold of (a, 0). This means that the region of coexistence of two species is reduced due to the Allee effect. Figure 9. Phase portrait of system (7) for case 3 when (a) p > p H and near p H , there is a limit cycle inside the stable manifold of E * R ; (b) p < p hom and near p hom , there is a limit cycle inside the stable manifold of E * R ; (c) p = p hom , the upward unstable manifold and the left stable manifold of E * R collide, which forms a homoclinic cycle; (d) the upward unstable manifold of E * R , the unstable manifold of E a on the u−axis, and the downward unstable manifold of E * R forms a loop of heteroclinic orbits among E * R , E a , and E 1 when p = p # ; and (e) 1 < p < 1 a .

Diffusion-Driven Turing Instability and Turing-Hopf Bifurcation
In the previous section, we have discussed both local bifurcation and global dynamics of kinetic model (7). In this section, we consider the effect of the diffusion on the stability of the interior constant steady states; thus, we set τ 1 = τ 2 = 0 in system (4). If an equilibrium is linearly stable in the absence of diffusion, and it becomes unstable in the presence of diffusion, we call such an instability Turing instability. As E * R is always a saddle if it exists, thus Turing instability can only happen near E * in both cases: weak cooperation and strong cooperation. We first investigate the existence of Turing instability, then we consider Turing-Hopf bifurcation near E * .
For Neumann boundary condition, we define the real-valued Sobolev space The linearization of system (4) with τ 1 = τ 2 = 0 at the constant steady state E * (u * , v * ) is given by where It is well known that the eigenvalues of D∆ on X are −d 1 where From the previous section, in the absence of diffusion, i.e., a) ). If there exists an n ∈ N, such that ∆ n = 0 has roots with positive real part when p < p H , Turing instability occurs. As trJ E * < 0 when p < p H , and detJ E * > 0 always satisfies, from (15), we have T n < 0 for any n when p < p H . Therefore, when p < p H , the signs of the real parts of roots of (14) are determined by the signs of J n .
In fact, the curve of J n = 0 is a hyperbola on d 1 − d 2 plane, whose horizontal and vertical asymptotes are We only need to consider the right branch of the hyperbola, which intersects with the d 1 −axis at d 1 = detJ E * l 2 mpcu * v * n 2 . Thus, with respect to n, both the horizontal asymptote of the right branch of the hyperbola and the intersection with the d 1 −axis are decreasing. Denote Noting that u * > a+1 On the d 1 − d 2 plane, we call the boundary curve of stable region of steady state E * the Turing bifurcation curve l T , which is formed by a sequence of curve segments l T n (n = 1, 2, ...) where d 1,n is the intersection of d 2 = d T 2 (n + 1, d 1 ) and d 2 = d T 2 (n, d 1 ), and d 1,0 can be infinite.
When p = p H , T n < 0 for all n ∈ N, similar as the analysis in the case of p < p H , we can also get a Turing bifurcation curve, and we have the following conclusion.

Remark 7.
Assume, in the p − d 2 plane, d 2 = d T 2 (n, d 1 ) and p = p H intersect at a point TH, then the characteristic equation has a pair of purely imaginary roots and a zero root at TH. According to the general theory of [62], a Turing-Hopf bifurcation may appear. However, rigorous derivation for the normal form is impossible due to the implicit form of E * in this paper. We shall give some numerical illustrations to show the dynamics of the system near the Turing-Hopf bifurcation point.

The Dynamics of Diffusive System with Two Delays
In the previous section, we have found the conditions of stability and Turing instability of constant steady state E * of system (4) when τ 1 = τ 2 = 0. We know that if p < p H and d 2 > d T 2 (n, d 1 ) for all n ∈ N, the steady state E * is locally asymptotically stable. In this section, we investigate the diffusive system (4) under the assumption p < p H and d 2 > d T 2 (n, d 1 ), and focus on the effect of two delays on the dynamics of the diffusive system near E * .

Hopf and Double Hopf Bifurcation Induced by Two Delays
In this subsection, we investigate the existence of Hopf bifurcation induced by two delays by the method of stability switching curves given in [63], and give the condition of the existence of double Hopf bifurcation. The linearization of system (4) at the steady state E * (u * , v * ) is given by where u(x, t) and v(x, t) satisfy the homogeneous Neumann boundary condition.
From Wu [64], the corresponding characteristic equation of (18) is where I is a 2 × 2 identity matrix, M n = −n 2 /l 2 diag(d 1 , d 2 ), n ∈ N 0 = {0, 1, 2, ...}. Equation (20) can be written in the following form: where P 0,n (λ) = (λ + d 1 The characteristic equation with the form of (21) has been investigated by Lin and Wang [63]. They derived an explicit expression for the stability switching curves, on which there is a pair of purely imaginary roots for (21). Moreover, they gave a criterion to determine the crossing directions, i.e., on which side of the stability switching curve there are two more characteristic roots with positive real parts. Using this method, we can find all the stability switching curves in the (τ 1 , τ 2 ) plane and determine the crossing directions. We leave the details in Appendices C and D.
Moreover, we have the following Hopf bifurcation theorem with two parameters.

Theorem 10.
For each j ∈ {1, 2, · · · , N}, T j n , defined by (A9), is a Hopf bifurcation curve in the following sense: for any p ∈ T j n and for any smooth curve Γ intersecting with T j n transversely at p, we define the tangent of Γ at p by − → l . If ∂Reλ ∂ − → l | p = 0, and the other eigenvalues of (21) at p have non-zero real parts, then system (4) undergoes a Hopf bifurcation at p when parameters (τ 1 , τ 2 ) cross T j n at p along Γ.
The theorem can be proved by a similar method in [49].

Normal Form on the Center Manifold for Double Hopf Bifurcation
In this subsection, applying the normal form method of partial functional differential equations [65], we derive normal form of double Hopf bifurcation taking two delays as bifurcation parameters. Then, we can classify the topological structures of various bifurcating solutions, and get the dynamics in the neighborhood of the double Hopf bifurcation point in system (4).
For the real-valued Sobolev space X defined in (12), define the complexification space of the real-valued Hilbert space X by with the general complex-value L 2 inner product Let C := C([−1, 0], X C ) denote the phase space with the sup norm, and write u t ∈ C for u t (θ) = u(t + θ), −1 ≤ θ ≤ 0.
Without loss of generality, we assume τ 1 < τ 2 in this section. Denote the double Hopf bifurcation point by 1 and σ 2 = τ 2 − τ * 2 as two bifurcation parameters. Drop the bars, and denote U(t) = (u(t), v(t)) T , then system (4) can be written as where with A, B, C and D being defined in (19), and Consider the linearized system of (23) We know that the normalized eigenfunctions of D∆ on X corresponding to the eigenvalues −d 1 n 2 l 2 and −d 2 Write system (23) as where System (25) can be rewritten as an abstract ordinary differential equation on C [65] where A is the infinitesimal generator of the C 0 -semigroup of solution maps of the linear equation (23), defined by with dom(A) = {ϕ ∈ C :φ ∈ C , ϕ(0) ∈ dom(∆)}, and X 0 is given by X 0 (θ) = 0 for θ ∈ [−1, 0) and X 0 (0) = I. Clearly, A : C 1 0 ∩ C → C . Then, on B n , the linear equation is equivalent to the retarded functional differential equation on C 2 : By the Riesz representation theorem, there exists a matrix whose components are bounded variation functions η k ∈ BV([−1, 0], R 2×2 ) such that Let A k (k = 1, 2) denote the infinitesimal generator of the semigroup generated by (28), and A * k denotes the formal adjoint of A k under the bilinear form The calculations of normal form are very long, so we leave them in supplementary materials. Based on the derivation in the Supplementary Materials, the normal form truncated to the third order on the center manifold for double Hopf bifurcation is obtained. Making the polar coordinate transformation, then we obtain the following system corresponding to the truncated normal formρ 1 where ν 1 = 1 (ReB 11 σ 1 + ReB 21 σ 2 ), ν 2 = 1 (ReB 13 σ 1 + ReB 23 σ 2 ),

Numerical Simulations
In this section, we carry out some numerical simulations for the kinetic system (7), diffusive system without delays, and the diffusive system (4) with two delays, respectively. These numerical simulations confirm and illustrate our main theoretical results.

Numerical Simulations for Kinetic System
In this subsection, we carry out some numerical simulations for system (7). We fix the parameters as a = 0.23, r = 1.1, m = 0.31.

Numerical Simulations for Kinetic System with Weak Cooperative Hunting
For system (7), choose c = 0.25 such that c < 1 r(1−a) . To show the complex dynamics of system (7), different values of p are chosen, listed in Table 1 and 2. The corresponding phase portraits for different values of p have been illustrated in Figures 3 and 4, which is drawn by pplane8 [66].        Using the method in [67], we can get Hopf bifurcation point p H = 1.5432, and the first Lyapunov coefficient is −1.1211, which means that the bifurcating periodic solution is asymptotically stable, and it is bifurcating from E * as p increases past p H from Theorem 2. We can draw the bifurcation diagram in Figure 10. If p < p H = 1.5432, E * is asymptotically stable (the black curve in Figure 10a). The increase of p induces instability to E * . When p > 1.5432, E * becomes unstable (the red curve in Figure 10a), and a unique limit cycle appears (marked by LPC). The period of limit cycle is increasing as p increases from p H , which tends to infinite as p → 1.6491 (see Figure 10b). When p = p # = 1.6491, the limit cycle coincides with a loop of heteroclinic orbits. When p > p # = 1.6491, the loop of heteroclinic orbits breaks and the limit cycle disappears. The corresponding phase portrait has been drawn in Figure 4, and the results are consistent with Theorem 4. Moreover, the amplitude of oscillation increases as p increases, and the minimum of the predator population is close to zero when p is close to p # = 1.6491, which will increase the risk on the extinction of predator.

Numerical Simulations for Kinetic System with Strong Cooperative Hunting
In Section 3.2, we discuss three different cases of dynamics when c is chosen as different values. To show case 1, we fix c = 3 and vary p as a bifurcation parameter. We get p H = 1.1024, and the first Lyapunov coefficient is −2.3280. Choosing different values of p listed in Table 3, we have drawn the corresponding phase portraits for different values of p in Figure 7.   Figure 9.  To show case 2, choosing c = 5 and varying p, we can get p H = 0.9608, and the first Lyapunov coefficient is −2.9553. Choosing p as listed in Table 2, the corresponding phase portraits have been illustrated in Figure 8.
To show case 3, choosing c = 8 and varying p, we can get p SN = 0.7665, p H = 0.8306, and the first Lyapunov coefficient is −3.732. We can draw the bifurcation diagram in Figure 11. When p < 0.7665 (i.e., on the left of LP), there is no interior equilibria. When p > 0.7665 (i.e., on the right of LP), there are two equilibria E * (the black curve) and E * R (the blue dotted curve). If 0.7665 < p < 0.8306, E * is asymptotically stable (the blue curve). The increase of p induces instability to E * . When p > 0.8306, E * becomes unstable (the red curve), and a unique limit cycle appears (marked by LPC in Figure 11a). When p increasing from p H , the period of limit cycle is increasing, and it tends to infinite as p → 0.8919 (see Figure 11b). Moreover, the amplitude of oscillation increases as p increases, and the amplitude line touches the curve of v−value of E * R . It means that there is a cycle goes through E * R , which forms a homoclinic cycle. The limit cycle coincides with a homoclinic cycle when p = p hom = 0.8919. When p > 0.8919, the homoclinic cycle breaks, and the limit cycle disappears. Choosing p as listed in Table 3, the corresponding phase portraits have been illustrated in Figure 9.
(a) (b) Figure 11. The bifurcation diagram with a = 0.23, m = 0.31, r = 1.1, and c = 8, and the conversion rate p as the bifurcation parameter. (a) The black curve is the stable steady state E * , the red curve is the unstable steady state E * , and the dotted blue curve is the unstable steady state E * R . The blue curves denote the v−amplitude of the limit cycle, which begin at the Hopf bifurcation point p = p H = 0.8306. (b) The period of the limit cycles.

Numerical Simulations for Diffusive System without Delays
In this section, we carry out some numerical simulations for diffusive system (4) with τ 1 = τ 2 = 0. We fix the parameters as To illustrate the dynamics of (4) with no delays in the case of weak cooperation, we choose c = 0.25 such that c < 1 r(1−a) . Let p = 1.4 < p H = 1.5432, from the previous section, E * (0.6882, 0.1516) is stable when d 1 = d 2 = 0. On d 1 − d 2 plane, we can draw curves determined by (16) for n ∈ N, and we only illustrate four curves for n = 2, 3, 4, 5. Correspondingly, from (17), we can get four curve segments, forming partial Turing bifurcation curve l T (see Figure 12a). From Theorem 9, if we choose d 1 and d 2 in the region above the curve, the constant steady state E * is stable. If we choose d 1 = 1.496 and d 2 = 0.000688 in the region under l T , the steady state E * of the diffusive system is Turing unstable, and there is a spatially inhomogeneous steady state (see Figure 13).
If we fix d 1 = 1.496, from (16), we can draw Turing bifurcation curve on p − d 2 plane (see Figure 12b), and d 2 = d T 2 (5, d 1 ) intersects line p = p H at TH (1.5432, 0.0009). Choosing p = 1.5432, d 2 = 0.00081 near TH, two stable spatially inhomogeneous periodic solutions coexist when we choose two different initial values (see Figure 14). To illustrate the dynamics of (4) without delays in the case of strong cooperation, we choose c = 8 such that c > 1 r(1−a) . Let p = 0.8 < p H = 0.8306. From (16) and (17), we can draw partial Turing bifurcation curve l T on d 1 − d 2 plane (see Figure 15a). The constant steady state E * (0.7392, 0.0864) is stable when we choose d 1 and d 2 in the region above l T . If we choose d 1 = 0.13 and d 2 = 0.006 in the region below l T , from Theorem 9, the steady state E * of the diffusive system is Turing unstable. If we fix d 1 = 0.13, from (16), we can draw Turing bifurcation curve on p − d 2 plane (see Figure 15b), and d 2 = d T 2 (5, d 1 ) intersects line p = p H at TH (0.8306, 0.009). Choosing p = 0.8306, d 2 = 0.0065 near TH, we can illustrate two stable spatially inhomogeneous periodic solutions coexisting when we choose two different initial values (see Figure 16).

Numerical Simulations for Diffusive System with Two Delays
In this section, we carry out some numerical simulations for diffusive system (4). Fix which is the same as in (31). Fix c = 0.25 such that c < 1 r(1−a) , and choose p = 1.2 < p H , d 1 = 0.3, and d 2 = 0.4, such that d 2 > d T 2 (1, d 1 ) for all n ∈ N. From Theorem 9, one can get the unique constant steady state E * (0.80945, 0.11797) is locally asymptotically stable. Now, we illustrate the effect of two delays on the dynamics (4). Following the steps in Appendices C and D, we can draw all the stability switching curves on (τ 1 , τ 2 ) plane, and decide the crossing direction, which are shown in Figures A3-A7 in Appendix E. Combining all the stability switching curves together, and zooming in the part of (τ 1 , 20], we get the Hopf bifurcation curves shown in Figure 17a. Consider the bottom left region bounded by left-most curve of T 1 0 and the lowest curve of T 2 0 (see Figure 17a), which are both part of T 0 . The crossing directions of the two switching curves (the black line and blue line) are all pointing outside of the region. From Theorem 9, under the assumption p < p H and d 2 > d T 2 (n, d 1 ), E * is locally asymptotically stable. From Theorem 10, when (τ 1 , τ 2 ) cross from the left to the right of the stability switching curve τ 2(1) 0 (or cross from the bottom to the top of τ 1(1) 0 ), the characteristic equation has two more characteristic roots with positive real parts from Lemma A1, E * becomes unstable, and periodic solution occurs due to Hopf bifurcation. Similarly, when (τ 1 , τ 2 ) cross from the left to the right of stability switching curve τ 1(1) 1 (or from the left to the right of stability switching curve τ 1(1) 2 ), the characteristic equation has two more characteristic roots. Then, E * is still unstable, and system (4) undergoes Hopf bifurcation and periodic solution arises. However, the periodic solution is unstable, as there is unstable manifold. It means that stability switching curves are the boundary of stable region for E * . In the stable region in Figure 17a, the predator and the prey coexist, and when (τ 1 , τ 2 ) cross the stability switching curves  Figure 17a). The positive equilibrium E * is unstable when τ 2 ∈ [0, 2.31) and there are periodic solutions; E * becomes locally asymptotically stable and periodic solutions disappear when τ 2 ∈ (2.31, 10.41), and E * loses its stability when τ 2 > 10. 41 and periodic solution appears again. Multiple times stability switching means that the stability of E * can switch with the varying of parameters (τ 1 , τ 2 ), and periodic solutions may appear or disappear during the stability switching of E * .  Figure 17b. In region D 1 , the positive equilibrium is asymptotically stable. In region D 2 or D 6 , there is a stable periodic solution. When the parameters cross into the region D4, there are two stable homogeneous periodic solutions coexisting in D4.

Concluding Remarks
We investigate the dynamics of a diffusive predator-prey model with two delays. First, we discuss the existence of constant coexisting steady states of system (4). It turns out that the strength of cooperation c = 1 r(1−a) is a threshold: if c < 1 r(1−a) , i.e., the cooperation is weak, there is at most one coexisting steady state, and if c > 1 r(1−a) , i.e., the cooperation is strong, there are at most two coexisting steady states.
To figure out the dynamics of system (4), we first carry out the study on the kinetic predator-prey model (7) with Allee effect in prey and cooperative hunting in predators. Compared with the work of Jang et al. [17], which mainly gives the local bifurcation analysis, we mainly focus on the global dynamics. We start with the kinetic system in the case of weak cooperation. Taking p, the conversion rate, as the bifurcation parameter, we prove that there is a loop of heteroclinic orbits between E 1 and E a formed by the intersection of the stable manifold of E a and the unstable manifold of E 1 when p = p # (p # > 1). It shows that p # is an important threshold, which distinguishes the globally stability for E 0 and the existence of a separatrix curve. To be specific, when p > p # , E 0 is globally stable, which means the extinction of both species, and we call the phenomenon overexploitation. When p < p # , there is a separatrix curve Γ s p , which is determined by the stable manifold of E a . High initial predator density will lead to extinction for both species, and conversely low initial predator density will approach either a steady state or a periodic oscillation.
When p < 1, E 1 (1, 0) is locally stable, which means that the prey achieves its carrying capacity and the predator is extinct. When 1 < p < p H , the stable state becomes E * , which implies the coexistence of both species. When p crosses p H , the stable state switches into a limit cycle. The period of the limit cycle increases with p, and tends to ∞ as p → p # . The existence of a limit cycle with a large amplitude increases the risk of the extinction for the predator.
In the case of strong cooperation, three different cases are considered. We prove the existence of a loop of heteroclinic orbits between E 1 and E a when p = p # in the first and second cases and nonexistence of such a loop of heteroclinic orbits in the third case. We also exhibit the existence of loop of heteroclinic orbits among E * R , E a , and E 1 formed by the unstable manifold of E * R and the unstable manifold E a when p = p # , and a homoclinic cycle when p = p hom by numerical simulations in case 3. When p > p # in case 1 and case 2 (p > p # in case 3), both species are distinct. When 1 < p < p # in case 1 and case 2, there is a separatrix curve Γ s p (E a ), and there are two stable states coexisting, which are extinction for both species and coexistence or oscillation. When p < p # (or p < p # ) and p SN < p < 1, there are two separatrix curves, the stable manifold of E a and the stable manifold of E * R , separating the first quadrant into three parts. If the initial population is above the stable manifold of E a , both species are extinct; if the initial population is inside the stable manifold of E * R , the stable state can be E * or limit cycle; if the initial population is below the stable manifold of E a and outside the stable manifold of the stable manifold of E * R , the prey achieves its carrying capacity and the predator is extinct. When p < p SN , there is one separatrix curve Γ s p (E a ), and there are two stable states coexisting. The species with initial population above the separatrix finally become extinct, however only the prey survive when initial population is below it.
There is a unique equilibrium for weak cooperation when 1 < p < 1 a . It is shown that there are at most two stable states coexist for different initial values separated by the stable manifold of E a . There are two interior equilibria E * and E * R in the case of strong cooperation when p SN < p < 1. The equilibrium E * R brings one more separatrix curve than the case of weak cooperation, thus there may be three stable states coexisting.
Second, we consider the diffusive system without delays, and focus on the effect of diffusion on the dynamics. Taking diffusion coefficients d 1 and d 2 as bifurcation parameters, we give the conditions of existence of Turing instability and Turing-Hopf bifurcation for the system with both weak and strong cooperation. We illustrate complex dynamics of the diffusive system, including the existence of spatially inhomogeneous steady state, coexistence of two spatially inhomogeneous periodic solutions.
Finally, we discuss the joint effect of two delays on the dynamical behavior of the diffusive system. Applying the method of stability switching curves, we find all the stability switching curves at which the characteristic roots are purely imaginary. Combining the geometrical and analytic method, we can decide the crossing direction of the characteristic roots as long as we confirm the positive direction for stability switching curves. Then, we can get the condition of existence of Hopf bifurcation. By searching the intersection of stability switching curves near the stable region, we get the double Hopf bifurcation point. Through the calculation of normal form of the system, we get the corresponding unfolding system and the bifurcation set, from which we can figure out the complete dynamics near the bifurcation point on the (τ 1 , τ 2 ) parametric plane. We theoretically prove and illustrate that near the bifurcation point, there are the phenomena of the stability of positive equilibrium: stable periodic solutions and the coexistence of two periodic solutions. From double Hopf bifurcation analysis, we can gain an insight into the local dynamics near the double Hopf bifurcation point. However, numerical simulations tells us that the occurrence of complex phenomena near bifurcation point is sensitive to the choice of initial values. When we choose the initial values with high ratio of predator to prey, both species are extinct. Thus, the complex phenomena near the double Hopf bifurcation probably accompany with an extinction steady state, which is consistent with the results in the kinetic system. Therefore, the analysis of global dynamics of kinetic system can help us comprehend these bistable dynamics for the diffusive system.

Conflicts of Interest:
The authors declare no conflict of interest.
p H < p < 1 a and near p H , E * is an unstable spiral. Thus, the set K = {p ∈ (1, 1 a ) : S(p) > v * } is a nonempty open set containing (1, p H + ε). Now we show that S(p) is monotonically decreasing on any component of K. Let p 1 < p 2 be two points in some interval in K. Denote the v-nullcline for p 1 and p 2 as v = g 1 (u) and v = g 2 (u). From the expression of v = g(u), we know that the curve of v = g 1 (u) is on the right of v = g 2 (u). The stable manifold of E a (a, 0) for p 1 and p 2 are graphs of function v 1 (u) and v 2 (u) defined for u > a, respectively. The corresponding eigenvectors at E a are Let u be the smallest such value. Then, we must have v 2 (u) ≥ v 1 (u) ≥ 0. However, which is a contradiction. Thus, the intersection of v 2 (u) and g 2 (u), S(p 2 ), is below the intersection of v 1 (u) and g 2 (u). For v 1 (u) between g 2 (u) and g 1 (u), from the vector field on the left of g 1 (u), the intersection of v 1 (u) and g 1 (u), S(p 1 ), is higher than that of v 1 (u) and g 2 (u). Thus, S(p 2 ) < S(p 1 ). Notice that the argument above shows that if p 2 ∈ K and S(p 2 ) > v * , then any p ∈ (1, p 2 ) also belongs to K and S(p) > S(p 2 ).
We can claim now K = (1, p a ) for some p a ∈ (p H , 1 a ). For p ∈ (1, p a ), S(p) > v * , and for p ∈ [p a , 1 a ), S(p) = v * . It remains to prove S(p) is decreasing when p ∈ [p a , 1 a ). From the vector field in (7), Γ s p moves towards the upper right, backward, before it meets the v-nullcline. Then, for p = p top , we have S(p) > v * , and thus p top < p a . It is obvious that v * is decreasing with respect to p for p ∈ (p top , 1 a ), and thus so does for p ∈ [p a , 1 a ). Therefore, S(p) is decreasing for both (1, p a ) and [p a , 1 a ). Obviously, p # is a threshold value of the property of Γ s p and Γ u p . When p < p # , Γ s p is above Γ u p , and when p > p # , Γ s p is below Γ u p . Let Ω denotes the bounded open subset of the positive quadrant, bounded by Γ s p , v-nullcline between Γ s p and Γ u p , Γ u p , and the segment from E 1 to E a on the u−axis (see Figure A1).  Figure A1. The dashed curve is the u-nullcline, and the dotted curve is the v-nullcline. Different positions of the stable manifold Γ s p of E a and the unstable manifold Γ u p of E 1 for (a) 1 < p < p # and (b) p # < p < 1 a .
Proposition A2. Suppose that c < 1 r(1−a) . (i) If 1 < p < p # , then S(p) > U(p). Moreover, all orbits in the positive quadrant above Γ s p converge to E 0 , and all orbits below Γ s p have their ω−limit sets in Ω, which is a positive invariant set.
(ii) If p # < p < 1 a , then S(p) < U(p), and Γ u p tends to E 0 . Moreover, all orbits in the positive quadrant above Γ u p converge to E 0 , and all orbits below Γ u p have their α−limit sets in Ω, which is a negative invariant set.
Proof. We prove the case for 1 < p < p # , and the other case can be proved similarly. Propositions A1 and 2 directly lead to S(p) > U(p). From Proposition A1, Γ s p enters E a from the region above the u-nullcline v = f (u), and meets the v-nullcline at (u s p , S(p)). On the right of (u s p , S(p)), Γ s p is still above v = f (u), since Γ u p is above the u-nullcline and S(p) > U(p). Moreover, Γ s p can not tend to E 1 . In fact, from the direction of the vector field in system (7), Γ s p goes to the lower right as t → −∞. Thus, Γ s p divides the first quadrant into two regions. Noticing that the first quadrant is invariant, we can get the results from Poncaré-Bendixson theorem.
Proposition A3. Suppose that c < 1 r(1−a) . (i) If 1 < p < p # , then either E * is stable or Ω contains a periodic orbits which is stable from the outside (both may be true). (ii) If p # < p < 1 a , then either E * is unstable or Ω contains a periodic orbits which is unstable from the outside (both may be true).
From Proposition A3, in order to figure out the dynamics in Ω in detail, we have to consider the existence and nonexistence of periodic orbits. In Section 3.1.1, we have discussed the existence of periodic orbits bifurcating from Hopf bifurcation. Now, we consider the nonexistence of periodic orbits.
there is no periodic orbits, and Γ s p connects E * to E a , which is a heteroclinic orbit. (ii) There is an ε 2 > 0 such that if 1 < p < 1 + ε 2 , there is no periodic orbits, and Γ u p connects E 1 to E * .
Proof. We prove the case for p < 1 a and near 1 a , and the latter case can be proved similarly. If S(p) = v * , Γ s p connects E * to E a , and there is no periodic orbits. When S(p) > v * , denote the intersection of Γ s p and v-nullcline as P(u s p , S(p)). It is easy to confirm that there is a p 1 such that S(p) = f (u p 1 ), where f (u p 1 ) is the component of the intersection of u-nullcline and v-nullcline corresponding to p 1 (denoted by P 1 (u p 1 , f (u p 1 )). Consider the region with vertices E a , P, P 1 , Q 1 ( 1 p 1 , 0), which is a negative invariant region. P 1 is well defined when p is close to 1 a as S(p) < f (p top ). As E * is the unique equilibrium in this region, thus if there are periodic orbits, they must encircle E * and lie wholly in this region. However, the divergence of the vector field of (7) is positive for p → 1 a − , as it is ra(1 − a) + m(pa − 1) at E a , which is positive. From Bendixson's criterion, there is no periodic orbits in the region. Therefore, according to Poincaré Bendixson theorem, E * is the α−limit set of Γ s p .
Theorem 3 provides a global description of the dynamical behaviour of system (7) for p ∈ R + \(1 + ε 2 , 1 a − ε 1 ) = R + \I. When p ∈ I, we can go further with the method in [30].
Proposition A5. Suppose that c < 1 r(1−a) . (i) If p H < p # , then (a) for p H < p < p # , there is at least one periodic solution in Ω which is stable from the outside and one (perhaps the same) stable from inside; (b) if a(p H ) > 0, then there is an ε > 0 such that if p H − ε < p < p H , there are at least two distinct periodic solutions, the inner of which is unstable while the outer is stable from the outside.
(ii) If p # < p H , then (a) for p # < p < p H , there is at least one periodic solution in Ω which is unstable from the outside and one (not necessarily distinct) unstable from inside; (b) if a(p H ) < 0, then there is an ε > 0 such that if p H < p < p H + ε, there are at least two distinct periodic solutions, the inner is stable while the outer is unstable from the outside.
Proof. We prove (i), and (ii) can be proved analogously. If p H < p < p # , E * is a repellor.
As Ω is a positive invariant region, (i 1 ) follows from the Poincaré-Bendixson theorem. If a(p H ) > 0, an unstable periodic solution bifurcates from E * as p decreases past p H .
Again Ω is a positive invariant region, from Poincaré-Bendixon theorem, there is a second periodic orbit exterior to the unstable bifurcating periodic orbit.
Appendix A.2. The Proof of Theorem 3 Proof. (i) If p ≥ 1 a , there is no interior equilibrium, then there is no periodic orbit in the first quadrant. Thus, every orbit converges to a boundary equilibrium. We have known that E 0 is a stable node, E 1 is a saddle, and E a is a unstable node if p > 1 a and a non-hyperbolic repellor if p = 1 a . Noting that the first quadrant is positive invariant, then E 0 is globally asymptotically stable.
(ii) and (iii) follows from Propositions A2 and A4 and the fact that E * is a sink when p > 1 and close to 1 and a source when p < 1 a and close to 1 a . (iv) If p decreases to p = 1, E * collides with E 1 . A transcritical bifurcation involving E * and E 1 occurs. E 1 changes its stability from a saddle point to a stable node. If p < 1, there is no interior equilibrium, and there is no periodic orbit in the first quadrant. E 0 and E 1 are both stable node. E a is a saddle, and its stable manifold Γ s p divides the first quadrant into two regions. The region above Γ s p is the attractive basin of E 0 , and the region below Γ s p is the attractive basin of E 1 .

Appendix A.3. Proof of Theorem 4
Proof. If a(p H ) < 0, from Theorem 2, stable periodic orbits bifurcating from Hopf bifurcation appear when p > p H . According to Proposition A5 (ii), if p # < p H , there will be at least two distinct periodic orbits, which contradicts our assumption about the uniqueness of periodic orbit. Thus, we have p H ≤ p # . As E * is asymptotically stable for any 1 < p < p H , if there is a periodic orbit, it must be unstable from inside, which contradicts our assumption about the orbital stability of periodic orbit. Therefore, there is no periodic orbit when 1 < p < p H , and the conclusion (i) follows from Proposition A2 (i). For the case of p # < p < 1 a , any orbit must encircle E * and lie wholly in Ω, which is negatively invariant by Proposition A2 (ii). Thus, if there is a periodic orbit, it must be unstable from outside, which contradicts our assumption about the orbital stability of periodic orbit. Therefore, the conclusion (iv) follows from Proposition A2 (ii). Moreover, together with the fact that Hopf bifurcating periodic orbit appears when p > p H , the nonexistence of periodic orbits for p # < p < 1 a implies that p H is strictly less that p # . If p H < p < p # , E * is a repellor, thus it follows from Proposition A3 (i) that there is a periodic orbit for all p H < p < p # . Noticing that we have proved there is no cycle for p > p # , according to the work of [68][69][70], the period of the unique cycle must tend to infinity as p → p # . Proposition A2 (i) completes the proof of (ii).
When p = p # , then p > p H , and E * is a repellor, from Poincaré-Bendixson theorem, the loop of heteroclicic orbits is the ω-limit set of orbits through any point in its interior.

Appendix B. Some Results and Proofs for the Case of Strong Cooperation
First, we consider the properties of the stable manifold of E a and unstable manifold of E 1 , denoted by Γ s p (E a ) and Γ u p (E 1 ), respectively.
Proposition A6. Suppose that c > 1 r(1−a) , and p SN < p < 1 a . (i) The orbit Γ s p (E a ) meets the v-nullcline v = g(u) at a point (u s p , S a (p)), where S a (p) ≥ v * := f (u * ), and S a (p) is a monotone decreasing function for p ∈ (p SN , 1 a ).
Proof. For the negative eigenvalue λ 1 of the Jacobian of E * R , the corresponding eigenvector is (1, mpv(1+cv) −mpcuv+λ 1 ). Noticing that the tangent vector of v = g(u) at E * R is mpv(1+cv) −mpcuv , the left part of Γ s p (E * R ) near E * R is below the v-nullcline, and the right part of Γ s p (E * R ) near E * R is above the v-nullcline. Obviously, the right part of Γ s p (E * R ) enters E * R from the lower right. From the vector field for (7), before the left part of Γ s p (E * R ) meets the v-nullcline, the curve under the u-nullcline directs lower right; the curve above the u-nullcline directs lower left. As it cannot cross the stable manifold Γ s p (E a ), thus, it is bounded before it meets the v-nullcline. Thus, the left part of Γ s p (E * R ) must meet the v-nullcline at a point, denoted by (u s p , S R (p)). Obviously, S R (p) ≥ v * .
(ii) can be proved similarly. Let . Similar as Proposition A7, we have the following conclusion.
Proposition A9. (i) If S R (p) > U R (p). All orbits inside the region boundary with Γ s p (E * R ) have their ω−limit sets in Ω 22 , which is a positive invariant set. (ii) If S R (p) < U R (p). All orbits below the left Γ u p (E * R ) have their α−limit sets in Ω 22 , which is a negative invariant set.

Appendix C. Stability Switching Curves
Applying the method proposed in Lin and Wang [63] to find the stability switching curves, we need to verify the assumptions (i)-(iv) in [63].
Lemma A1. As the delays (τ 1 , τ 2 ) vary continuously in R 2 + , the number of zeros (counting multiplicity) of D n (λ; τ 1 , τ 2 ) on C + can change only if a zero appears on or cross the imaginary axis.
Thus, for the stability switching curves corresponding to Ω j,n , T +j j 1 ,j 2 ,n is connected to T −j j 1 +δ a 1 ,j 2 −δ a 2 ,n at one end a j,n , and connected to T Definition A1. Any (τ 1 , τ 2 ) ∈ T n is called a crossing point, which makes D n (λ; τ 1 , τ 2 ) = 0 have at least one purely imaginary root iω with ω belongs to the crossing set Ω n . The set T n is called stability switching curves.

Appendix D. Crossing Directions
When (τ 1 , τ 2 ) varies and crosses the stability switching curves from one side to the other, the number of characteristic roots with positive real part may increase. We call it the crossing direction of stability switching curves.
In order to describe the crossing direction clearly, we need to specify the positive direction for stability switching curves T ±j j 1 ,j 2 ,n . We call the direction of the curve corresponding to increasing ω ∈ Ω j,n the positive direction, i.e., from (τ ±j 1,j 1 ,n (a j,n ), τ ∓j 2,j 2 ,n (a j,n )) to (τ ±j 1,j 1 ,n (b j,n ), τ ∓j 2,j 2 ,n (b j,n )). From the fact that T +j j 1 ,j 2 ,n is connected to T −j j 1 +δ a 1 ,j 2 −δ a To make our expression clear, we draw a schematic diagram of a part of stability switching curves corresponding to Ω j,n = [a j,n , b j,n ] (see Figure A2). The solid curve stands for T +j j 1 ,j 2 ,n , with two ends A(τ +j 1,j 1 ,n (a j,n ), τ −j 2,j 2 ,n (a j,n )), and B(τ +j 1,j 1 ,n (b j,n ), τ −j 2,j 2 ,n (b j,n )). The dashed curve denotes T −j j 1 +δ a 1 ,j 2 −δ a 2 ,n , which is connected to T +j j 1 ,j 2 ,n at A corresponding to a j,n , with the positive direction from A to C.
Thus, we can get a further conclusion which is more intuitive.
Remark A1. The region on the right of T +j j 1 ,j 2 ,n has two more characteristic roots with positive real parts, and the region on the right of T −j j 1 ,j 2 ,n has two fewer characteristic roots with positive real parts.

Appendix E. The Figures of Stability Switching Curves
Following the steps in Appendices C and D, we can draw all the stability curves on (τ 1 , τ 2 ) plane and decide the crossing direction.