Dynamics of a Leslie–Gower predator–prey system with cross-diffusion

A Leslie–Gower predator–prey system with cross-diffusion subject to Neumann boundary conditions is considered. The global existence and boundedness of solutions are shown. Some sufficient conditions ensuring the existence of nonconstant solutions are obtained by means of the Leray–Schauder degree theory. The local and global stability of the positive constant steady-state solution are investigated via eigenvalue analysis and Lyapunov procedure. Based on center manifold reduction and normal form theory, Hopf bifurcation direction and the stability of bifurcating timeperiodic solutions are investigated and a normal form of Bogdanov–Takens bifurcation is determined as well.


Introduction
In ecological systems, the interaction of predator and prey has abundant dynamical features although the investigations on predator-prey models has improved and lasted for several decades, which are based on the pioneering works of Lotka and Volterra [34]. Moreover, more realistic models are proposed in view of laboratory experiments and observations. Leslie and Gower [17] first proposed the following predator-prey model where u(t) and v(t) represent the densities of prey and predators at time t, respectively; the parameters a 1 , b 1 c 1 and d 1 are positive constants; the term d 1 v/u is called the Leslie-Gower terms, which measures the loss in the predator population due to rarity of its favorite food. System (1.1) is regarded as a prototypical predator-prey system in the ecological studies. But the interaction terms in (1.1) are unbounded, which is not reasonable in the real world. By using Holling type II functional response [13] in both prey and predator interaction terms, a Leslie-Gower predator-prey system with saturated functional responses is obtained and takes the form (see [4]): (1.2) The model (1.2) is based on the biological fact that if the predator v is more capable of switching from its favorite food (the prey u) to other food options, then it has better ability to survive when the prey population is low; here a 1 and a 2 are the growth rates per capita of prey u and predator v, respectively; b 1 measures the strength of intraspecific competition among individuals of species u, and it is related to the carrying capacity of the prey; c 1 is the maximum value of the per capita reduction rate of u due to v, and c 2 is the maximum growth per capita of v due to predation of u; k 1 and k 2 measure the extent to which environment provides protection to prey u and predator v, respectively. Non-monotonic responses appear at the microbial level; when the nutrient concentration reaches at a high level an inhibitory effect of the specific growth rate can occur [3,6]. This may frequently be noticed when micro-organisms are used for waste decomposition or for water purification. Andrews [3] suggested a response function p(u) = mu k 1 +k 2 u+u 2 , known as Monod-Haldane response function, to model such an inhibitory effect at high concentrations. In particular, Sokol and Howell [31] derived a simplified Monod-Haldane type p(u) = mu k 1 +u 2 . A Leslie-Gower predator-prey system with a Monod-Haldane functional response takes the form: (1. 3) In mathematical ecology, population may be distributed non-homogeneously, and the predators and preys naturally develop strategies for survival. Thus, we may introduce diffusive structure, which can be illustrated as different concentration levels of predators and preys causing different movements. Diffusion means the movement of individuals from a higher to a lower concentration region, while cross diffusion implies the population fluxes of one species owing to the presence of the other species. In this paper, our concern is the following system with cross-diffusion rates in Ω × (0, ∞), in Ω × (0, ∞), ∂u ∂n u(x, 0) = ϕ(x) ≥ 0, v(x, 0) = ψ(x) ≥ 0, in Ω, (1.4) whose corresponding ordinary differential equations (ODEs) is (1.3) with all the parameters b 1 , m, k 1 and k 2 equal to 1. Here ∆ denotes the Laplacian operator on R N (N ≥ 1), Ω is a connected bounded open domain in R N , with a smooth boundary ∂Ω, n is the outward unit normal vector on ∂Ω. The homogeneous Neumann boundary condition means that the two species have zero flux across the boundary ∂Ω. The diffusion terms d j , j = 1, 2 stand for natural dispersive force of movement of an individual, while β describes the mutual interferences between individuals and is usually referred as the cross-diffusion pressure measuring the situation that the prey keeps away from the predator; a and b are the growth rates per capita of prey u and predator v. The parameters a, b, d 1 and d 2 are positive constants and β is non-negative constant. In some cases, the quantity v is not influenced by any cross diffusion in the sense that the coefficient β in the second equation of (1.4) vanishes, that is, we ignore the population migration of predators due to the presence of preys. In this situation, Li et al. [20] considered the following reaction-diffusion system in the one-dimensional space domain Ω = (0, π): in Ω × (0, ∞), in Ω × (0, ∞), ∂u ∂n = ∂v ∂n = 0 on ∂Ω × (0, ∞), (1.5) where d is the relative diffusion rate of predator v when the diffusion rate of prey is rescaled to 1. Li et al. [20] studied the Hopf bifurcation and steady-state bifurcation by taking d as the bifurcation parameter and described both the global structure of the steady-state bifurcation from simple eigenvalues and the local structure of the steady-state bifurcation from double eigenvalues by using space decomposition and the implicit function theorem. The presence of the cross-diffusion term causes more abundant dynamic behaviors. For example, the effect of cross diffusion on dynamics of predator-prey models has been studied in [5,7,22,24,30,35,37,43,44,47]. The relevant discussion is a bit difficult and requires more techniques than for models without cross-diffusion. In [5,24,43,44], the researchers mainly obtained the non-existence and existence of non-constant positive steady-states (patterns) and showed cross diffusion can create non-constant steady states. Gambino et al. [7] analyzed the linear stability of the positive equilibrium of a competitive Lotka-Volterra system, and showed the cross-diffusion is the key mechanism for the formation of spatial patterns through Turing bifurcation. Liu et al. [22] not only obtained the global existence result of solutions under an appropriate parameter condition, but also gave explicit parameter ranges of the existence of non-constant positive steady-states.
For system (1.4), we first discuss the influence of the cross-diffusion coefficient β on the global existence of the solution. As far as global existence is concerned, many researchers have some relevant works, for example, [22,26,33,41]. Wu et al. [41] and Tao [33] analyzed the predator-prey model with prey-taxis and discussed the effect of the prey-taxis term on the global existence of solutions of the system. Mu et al. [26] studied the global existence of classical solutions to a parabolic-parabolic chemotaxis system, but there are strict restrictions on functions in the system. Liu et al. [22] investigated the global existence of solutions of a parabolic-elliptic two-species competition model with cross diffusion.
Next, for a predator-prey system, what we are interested in is whether the various species can exist and takes the form of non-constant time-independent positive solutions. In [5,8,24,25,43,44], the authors have established the existence of stationary patterns in some predatorprey models in the presence of self-diffusion and cross-diffusion. Our results are a little different from theirs. We not only prove the existence of non-constant solution of system (1.4) when the cross-diffusion β is sufficiently large, but also we find infinitely many intervals of d 1 > 0 near zero such that (1.4) admits at least one nonconstant solution if d 1 belongs to such intervals. Moreover, researchers have paid more attention to Hopf bifurcation and steady state bifurcation (cf. [9,10,15,18,19,36,42,46]), and investigated some predator-prey models without cross diffusion term. Only a few works [23,45] have concentrated on the Bogdanov-Takens bifurcation phenomena of diffusive predator-prey systems with delay effect. In this paper, we study the Bogdanov-Takens bifurcation by regarding the cross-diffusion term β as one of bifurcation parameters.
The organization of the remaining part of the paper is as follows. In Section 2 we prove the global existence and boundedness results of solutions to (1.4) and in Section 3 we obtain a priori bounds of nonnegative steady state solutions. In Section 4 we deal with the nonexistence of non-constant positive steady states for sufficient large diffusion coefficient and consider the existence of non-constant positive steady states for a small range of diffusion coefficient and sufficient large cross-diffusion coefficient by using the Leray-Schauder degree theory. Section 5 is devoted to the local and global stability of homogeneous steady states. Center manifold reduction and normal form theory are employed in Section 6 not only to discuss the existence of Hopf bifurcation but also to determine the Hopf bifurcation direction and the stability of bifurcating time-periodic solutions. In Section 7 we observe that system (1.4) exhibits Bogdanov-Takens bifurcation phenomena. Finally in Section 8, some conclusions are presented and numerical simulations are carried out to illustrate some previous theoretical results.
For convenience, we introduce the following notations. Let H k (Ω) (k ≥ 0) be the Sobolev space of the L 2 -functions f defined on Ω whose derivatives f (n) (n = 1, . . . , k) also belong to L 2 (Ω). Denote the spaces X = {φ ∈ H 2 (Ω)| ∂φ ∂n = 0 on ∂Ω} and Y = L 2 (Ω). For a space Z, we also define the complexification of Z to be Z C Z ⊕ iZ = {x 1 + ix 2 |x 1 , x 2 ∈ Z}. Define an inner product on the complex-valued Hilbert space Y 2 C by

Global existence and boundedness
In this section, we employ the method in [40] to obtain the global existence and boundedness of solutions of model (1.4). We need to establish some priori estimates. It is clear that the local existence of solutions to (1.4) was established by Amann [1]. This result can be summarized as follows.
Lemma 2.1. For each fixed p > N, assume that the initial data (ϕ, ψ) ∈ (W 1,p (Ω)) 2 satisfies ϕ ≥ 0 and ψ ≥ 0, then there exists a positive constant T max (the maximal existence time) such that (ϕ, ψ) determines a unique nonnegative classical solution (u( Proof. (i) The local existence of the solution to (1.4) follows from [1]. Denote by T max the maximal existence time of the solution. Next, we shall prove (2.1). On account of (1.4), we know that v(x, t) satisfies Clearly, v ≡ 0 is a sub-solution to problem (2.2). Hence, we can apply the maximum principle for parabolic equations to obtain that v(x, t) ≥ 0. Similarly, we can obtain u(x, t) ≥ 0. Also from (1.4) and v ≥ 0, we obtain that in Ω × (0, ∞), Then from comparison principle of parabolic equations, it is easy to verify u(x, t) ≤ c, where c is given in (2.1). This completes the proof of Lemma 2.1.
The above lemma means that, in the space W 1,p (Ω) each pair of the initial values ϕ and ψ can determine a unique nonnegative classical solution (u(x, t), v(x, t)), which is twice continuously differentiable with respect to x ∈ Ω and continuously differentiable with respect to t ∈ [0, T max ). Moreover, u(·, t), v(·, t) ∈ W 1,p (Ω) can be regarded as a continuous mapping with respect to t ∈ [0, T max ).
According to Amann's results [2], we need to establish the L ∞ bound of (u, v) in order to show its global existence. Based on Lemma 2.1, it is enough to establish the L ∞ bound of v(x, t). Firstly, we shall show that the solution v(x, t) is bounded in L 1 (Ω). In the proof, we need to use the following elementary inequality [39].
Then we haveU In view of Lemma 2.1, we know that 0 ≤ u(x, t) ≤ c for all (x, t) ∈Ω × [0, T max ) and hence that which, together with the Hölder inequality, implies thaṫ with r = max{a, b}. It follows from Lemma 2.2 that there exists a positive constant M such that U(t) + V(t) ≤ M for all t ∈ (0, T max ), and hence that there exists a positive constant C 0 such that (2.3) holds. The proof is completed.
Secondly, we will establish L p estimates for v(x, t) by using a weight function φ(u) similar to that in [32,38,41]. We now present some basic inequalities which will be used in the sequel (see [14,27]). In several places we shall need the following Poincaré's inequality: with arbitrary p > 1 and q > 0. Also, an essential role will be played by Gagliardo-Nirenberg interpolation inequality and consider a weight function Denote φ(u(x, t)) by φ(u), then we have It follows from system (1.4) that In virtue of (2.7), we know that φ (u), φ(u) > 0. Combining with v(x, t) ≥ 0, it is easy to see that Furthermore, using Young's inequality, we obtain (2.9) Similar to the above, we obtain (2.10) (2.11) Substituting (2.9), (2.10) and (2.11) into (2.8), we have Clearly, By a direct calculation, we obtain for 0 ≤ u ≤ c, where β and α satisfy (2.4) and (2.5) respectively, and Therefore, (2.14) It follows from (2.14) that (2.12) is simplified to be where C 1 = (bp + 2αac 2 )/p. By the Gagliardo-Nirenberg and Poincaré's inequality and (2.7) and (2.3), we have (2.17) Hence from (2.15) and (2.17) we obtain for all t ∈ (0, T max ), where 1 η > 1. By using Lemma 2.2 and (2.7), we conclude that there exists which is the desired result.
Finally, we establish the L ∞ bound of v(x, t) using Lemma 2.4.
It follows from Lemmas 2.4 and 2.1 that there exists a positive constant A 1 such that In virtue of (2.18) and the first equation of system (1.4) and the L p -estimate for parabolic equations, we obtain u(·, t) W 2 p (Ω) ≤ A 1 for all t ∈ (0, T max ). (2.19) This, together with the Sobolev embedding theorem (see [16]), yields We now turn to the second equation of (1.4), which can be rewritten as the non-divergence form: ∂v In virtue of Lemmas 2.4 and 2.1 and (2.19), we have Using (2.21), (2.20) and (2.22) and the L p -estimate for parabolic equations, we have v(·, t) W 2 p (Ω) ≤ A 1 for all t ∈ (0, T max ). Again, taking p to be sufficiently large and combing with the Sobolev embedding theorem (see [16] Hence, this proof is completed.
Obviously, from Lemmas 2.1 and 2.5 and [2], we conclude that Notice that in the proof of Lemma 2.1, for any positive constant ε 0 , there exists Hence we can replace c by a + ε 0 for t ∈ (t 1 , ∞). Similarly in Lemma 2.3, C 0 can be chosen to be independent of (ϕ, ψ). So Ω u(x, t) ≤ C 0 for t ∈ (t 2 , ∞) with t 2 > t 1 . Again in the proof of Lemmas 2.4 and 2.5, we can also replace c by a + ε 0 and then we can find where E and A are independent of (ϕ, ψ). In view of (2.23) and (2.24), there exists a constant where M 1 is independent of (ϕ, ψ). Therefore, we have the following theorem.

A priori estimates
Steady-state solutions of (1.4) satisfy the following system: It is easy to see that system (1.4) has a positive constant steady-state solution Next, we study the asymptotic behavior of positive solutions of (3.1) as d 1 is small or β is sufficiently large. For the first step of the asymptotic analysis, we derive a priori positive upper and lower bounds for positive solutions to (3.1).
for all x ∈Ω.
Proof. Let x 0 ∈Ω be a maximum point of u, i.e., u(x 0 ) = max x∈Ω u(x). Then by using the maximum principle [24] to the first equation of (3.1), one has a − u( By setting w = (d 2 + βu)v, we can reduce the second equation of (3.1) with the boundary condition to To obtain the lower bound for v, we define w(y 0 ) = minΩ w(x). Similarly, applying the maximum principle [24] to (3.2) yields v(y 0 ) ≥ b(1 + u 2 (y 0 )). According to the definition of w, we obtain Now, denote u(y 1 ) = minΩ u(x) for some y 1 ∈Ω. It follows from the maximum principle [24] that This, together with (3.3) and (3.4), implies that If a > κ then u(x) ≥ a − κ for all x ∈ Ω and hence the proof of Lemma 3.1 is completed.
In what follows, we shall show that u(x) ≥Č in the case where a ≤ κ. Let From Harnack's inequality (see [21]), there exists a positive constant C * such that Hence, it remains to prove that there is a positive constant ε such that maxΩ u( From the Sobolev embedding theorem and elliptic estimates, there exists a subsequence By the property of solutions of the logistic equation and min v n Then by dividing the first equation of (3.1) by u n L ∞ , we know that {ũ n } forms a sequence of positive solutions of Note that ũ n L ∞ = 1 for n ∈ N, then it follows from the elliptic regularity theory and the Sobolev embedding theorem that there exists a nonnegative functionũ ∞ ∈ C 1 (Ω) such that lim n→∞ũn =ũ ∞ in C 1 (Ω). This, combining with ũ ∞ L ∞ = 1, yields thatũ ∞ > 0. On the other hand, by integrating the first equation in (3.5) over Ω, we observe that Let n → ∞, and note thatũ ∞ > 0, u ∞ = 0 and v ∞ = b, then we have a = b, which contradicts our assumption. Therefore, we complete the proof of Lemma 3.1.

Existence/nonexistence of nonconstant solutions
Throughout the remaining part of this paper, we always assume that a > b.
Proof. For fixed a, b, β and Ω, Lemma 3.1 and standard regularity arguments tell that {(u n , v n )} ∞ n=1 has a convergent subsequence, which we still denote by {(u n , v n )} ∞ n=1 . According to the argument by Lou and Ni [24], we can obtain a positive constant K, which is independent of n, such that for n ∈ N. Together with Lemma 3.1, we can find a constantū ∈ [0, a] such that lim n→∞ u n =ū uniformly inΩ. Lemma 3.1 and the standard L p -estimate for elliptic equations mean that both {u n } ∞ n=1 and {v n } ∞ n=1 are uniformly bounded in W 2,p (Ω). Thus, the usual compactness argument implies lim n→∞ u n =ū in C 1 (Ω), (4.1) passing to subsequence. We can similarly get a nonnegative functionv such that passing to a subsequence. By setting n → ∞ in the weak form of the second equation of (3.1) and using the elliptic regularity theory, we know thatv satisfies Sinceū ∈ [0, a] is constant, the well-known property of the logistic equation implies thatv is also constant and satisfiesv Integrating the first equation of (3.1) yields   (d 2 , β, a, b, Ω) satisfying a > b, there exists a large positive constant D such that (3.1) with d 1 ≥ D has no nonconstant solutions.
Proof. Assume that (u, v) is a non-negative solution of (3.1) and denotē Then, multiplying u −ū the first equation in (3.1) by and integrating over Ω yield where the last inequality comes from Lemma 3.1. Recall the Poincaré-Wirtinger inequality where λ 1 is the least positive eigenvalue of −∆ with homogeneous Neumann boundary condition on ∂Ω. Then it follows from (4.5) that Similarly, multiplying by v −v the second equation of (3.1) and integrating the resulting expression lead us to By Lemma 3.1 and Young's inequality, for any ε > 0, one can find a positive constant K such that where κ is the positive number given in Lemma 3.1. Then the Poincaré-Wirtinger inequality implies Thus, when d 1 > 0 is large, (4.6) and (4.8) enable us to find a positive constant K 1 such that which implies that u is a constant if d 1 is large enough. Combining with (4.8) and (4.9), we deduce that (u, v) is a constant solution if d 1 > 0 is sufficiently large. Then the proof of Theorem 4.2 is completed. Recall that −∆ under Neumann boundary condition has eigenvalues 0 = λ 0 < λ 1 < · · · < λ n < · · · with lim n→∞ λ n = +∞. Let S i be the eigenspace associated with λ i with multiplicity n i . Let φ ij , 1 ≤ j ≤ n i , be the normalized eigenfunctions corresponding to λ i . Then the set (4.10) Next, we study the linearization of (3.1) at (u * , v * ), where e = (u * , v * ) is the unique positive constant solution of (1.4). Let Φ(U) = (d 1 u, d 2 v + βuv) T and in Ω, where C is a positive constant whose existence is guaranteed by Lemma 3.1. Note that the derivative Φ U (U) of Φ(U) with respect to U is a positive operator for all non-negative U, then Φ −1 U (U) exists and is a positive operator as well. Hence, U is a positive solution to (4.11) if and only if is a compact perturbation of the identity operator, the Leray-Schauder degree deg(F(·), 0, B) is well-defined if F(U) = 0 on ∂B. Note that D U F(e) = I − (I − ∆) −1 {Φ −1 U (e)G U (e) + I}, and that the index of F at e is defined as index(F(·), e) = (−1) γ provided that D U F(e) is invertible, where γ = ∑ m µ and m µ is the multiplicity of any negative eigenvalue µ of D U F(e); see [28] for more details.
We now consider the eigenvalues of D U F(e). First, for every integer i ≥ 0 and 1 ≤ j ≤ dimS i , X ij is invariant under D U F(e), and µ is an eigenvalue of D U F(e) on X ij if and only if it is an eigenvalue of the matrix Thus, D U F(e) is invertible if and only if, for all i ≥ 0, the above matrix is nonsingular. To calculate γ, we first define If H(λ i ) = 0, then for each 1 ≤ j ≤ dimS i , the number of negative eigenvalues of D U F(e) on X ij is odd if and only if H(λ i ) < 0. In conclusion, we have the following lemma (see [29]), which gives the explicit formula of calculating the index. To facilitate our computation of index(F(·), e), we will consider the sign of H(λ i ). Notice that our aim is to investigate the effect of the cross-diffusion coefficient β and diffusion coefficient d 1 on the existence of stationary patterns. Then we will concentrate on the dependence of H(λ i ) on β and d 1 . Note that Then, we have then H(λ) = 0 has two roots λ = λ + (β, d 1 ) and We first consider the dependence of H(λ) on β. When β is large enough, we have Λ > 0 and the two roots of H(λ) satisfy lim (4.14) Thus, we have the following existence result about the non-constant steady state solution: Theorem 4.5. Assume that a > b, b > (1+θ 2 )θ 1+3θ 2 andλ ∈ (λ n , λ n+1 ) for some n ≥ 1 and ∑ n i=1 n i (λ i ) is odd, then there exists a positive number β * such that system (3.1) with β ≥ β * has at least one non-constant positive solution.
Next we consider the dependence of H(λ) on d 1 . From the previous analysis, it follows that the roots of H(λ) = 0 are all negative if 2θb 1+θ 2 + βb d 2 +βθ − 1 < 0 and Λ(β, d 1 ) > 0. So, in this case, we can't obtain the existence of non-constant positive solutions by using the method of degree theory.

Stability of the positive constant solution
In this section, we firstly analyze the stability of the positive constant steady-state solution by eigenvalue analysis. And then, we will investigate the global stability of the positive constant steady-state solution. To investigate the local dynamical behavior of system (1.4) near the positive constant solution e, we need to consider the linearized operator L α 1 ,β of (1.4) with respect to (u, v) at (u * , v * ). Note that i=0 is a complete orthogonal base of X. Substituting them into the characteristic equation yields To investigate the stability of the positive steady-state solution, it suffices to study the characteristic equation detM(σ, α 1 , β, λ i ) = 0, that is, It is easy to know that two solutions of equation (5.1) have negative real parts if T i (α 1 , β) < 0 and D i (α 1 , β) > 0 for all i ≥ 0. Thus, we have the following results.
Lemma 5.1. If a > b, then all eigenvalues of L α 1 ,β have negative real parts, or equivalently, the homogenous steady-state e = (θ, b(1 + θ 2 )) is locally asymptotically stable, provided that one of the following conditions is satisfied: Now, we consider the global stability of e.
Proof. We discuss the global stability of e by Lyapunov method. Define Then where It is easy to see that 4d 1 . Hence, we have I 1 ≤ 0. Further computation gives Clearly, if 2b(a + c) < 1. Therefore, we have I 2 < 0. It follows from the above arguments that if the conditions of Lemma 5.2 are satisfied, then L (t) < 0 along all trajectories in the first quadrant except (u * , v * ). Therefore e = (u * , v * ) is globally asymptotically stable.

Hopf bifurcation
This section is devoted to the Hopf bifurcation at the nontrivial steady-state solution e = (u * , v * ) T of (1.4) with a > b.
To be more precise, as a pair of simple complex conjugate eigenvalues of the linearization around e = (u * , v * ) T cross the imaginary axis of the complex plane, the nontrivial steady-state solution e = (u * , v * ) T of (1.4) loses stability and a branch of small-amplitude limit cycles emerges from e = (u * , v * ) T . Throughout this section, we always assume that (H1) a > b, λ i is a simple eigenvalues of the linear operator −∆ subject to the homogeneous boundary condition ∂ ∂n u on ∂Ω, ϕ i is the eigenvector associated with λ i satisfying In what follows, by choosing the cross-diffusion coefficient β as the bifurcation parameter, we shall analyze the occurrence of Hopf bifurcation, the Hopf bifurcation direction and the stability of bifurcating time-periodic solutions. It follows from [11,12] that system (1.4) with a > b undergoes Hopf bifurcation near β = β i at the nontrivial steady-state solution e = (u * , v * ) T , where β i ∈ (0, ∞) satisfies Note that T i (α 1 , β) is monotone decreasing with respect to β, then it is easy to see that T i (α 1 , ·) has exactly one zero Next, we only need to verify D j (α 1 , β i ) = 0 for all j = i. Obviously, Therefore, we have D j (α 1 , β i ) < 0 for all j = i if < 0 and D j (α 1 , β i ) = 0 for all j = i if > 0 and Therefore, we shall consider Hopf bifurcation under the following assumptions: (H3) Either < 0 or > 0 and α 1 = α ± 1 .
For convenience, we call a Hopf bifurcation forward if there exist periodic solutions when parameter value β > β i ; and backward if β < β i . Under assumptions (H1), (H2) and (H3), L α 1 ,β i has exactly one pair of purely imaginary eigenvalues ±iω i with associated eigenvectors q i andq i , where ω i = D i (α 1 , β i ), q i = ρ i ϕ i , and the nonzero vector . Moreover, as β varies such that T i (α 1 , β) decreases and passes through 0, σ(β) varies from a complex number with a positive real part to a purely imaginary number and then to a complex number with a negative real part. This implies that a codimension one Hopf bifurcation for (1.4) occurs at β = β i . Namely, in every neighborhood of (U, β) = (e, β i ) there is a unique branch of time-periodic spatially non-homogeneous solutions U β (t, x), which tends to e as β → β i . The period T β of U β (t, x) satisfies that T β → 2π/ω i as β → β i . Under assumptions (H1), (H2) and (H3), −iω i is also an eigenvalue of L * α 1 ,β i with an associated eigenvector Next, we consider the bifurcation direction and stability of the bifurcating periodic solutions at β = β i according to [11,12]. Denote by G 2 = (G 2 1 , G 2 2 ) T , and G 3 = (G 3 1 , G 3 2 ) T the second-and third-order Fréchet derivatives of ∆Φ(U) + G(U) with respect to U at e = (u * , v * ), respectively. A straightforward computation yields for all ξ = (ξ 1 , ξ 2 ) T , (ζ 1 , ζ 2 ) T and ς = (ς 1 , ς 2 ) T ∈ X. It is well known that the following quantity determines the direction and stability of bifurcating periodic orbits U β (t, x) (see [11,12]) where and Therefore, we obtain the following result. Theorem 6.1. In addition to assumptions (H1), (H2) and (H3), a Hopf bifurcation for (1.4) occurs at β = β i if a > b. Namely, when a > b, in a neighborhood of (U, β) = (e, β i ) there is a branch of periodic solutions U β (x, t) satisfying U β (x, t) → e as β → β i . The period T β of U β (x, t) satisfies that T β → 2π/ω * as β → β i . Moreover, the bifurcation is backward (respectively, forward) if Re(Υ i ) < 0 (respectively, > 0).
Obviously, in Theorem 6.1, if λ i is not the principal eigenvalue of the linear operator −∆ subject to the homogeneous boundary condition ∂ ∂n u = 0 on ∂Ω, then the Hopf bifurcating periodic solutions U β (x, t) is spatially nonhomogeneous and unstable. However, if λ i is the principal eigenvalue λ 0 = 0, then the associated eigenvector ϕ 0 can be a positive constant function on Ω. In this case, assumption (H1) is obviously satisfied and α 1 − b is sufficiently close to zero. Hence, we can regard b as a bifurcation parameter. Obviously, we have T 0 (α 1 , β) = 0, θ 2 −1 and one of the following conditions is satisfied It is easy to evaluate σ(b) at b = b * to get Reσ (b * ) = θ 2 −1 2(1+θ 2 ) > 0. Thus, it remains to calculate the direction of Hopf bifurcation and the stability of bifurcating periodic orbits bifurcating from (U, b) = (e, b * ). In virtue of (6.1), we have

Corollary 6.2.
Under one of conditions (A1)-(A3), if a > b then in every neighborhood of (U, b) = (e, b * ) there is a branch of spatially homogeneous periodic solutions U b (x, t) satisfying U b (x, t) → e as b → b * and the Hopf bifurcation is forward (respectively, backward) and the bifurcation periodic solutions are orbitally asymptotically stable (respectively, unstable) if 2 − 3θ 2 + 6θ 4 − θ 6 < 0 (respectively,

Bogdanov-Takens bifurcation
Apart from the occurrence of Hopf bifurcation discussed so far, codimension 2 bifurcation such as Bogdanov-Takens bifurcation is also possible in system (1.4). In order to discuss codimension 2 bifurcation, in addition to taking β as a bifurcation parameter, we need another parameter. It is easy to see that T i (α 1 , β) depends on (β, α 1 θ, b) and D i (α 1 , β) on (β, b, α 1 , α 2 θ). More precisely, α 1 depends on θ and b, α 2 on θ. For convenience, we choose α 1 and β as bifurcation parameters. In this section, we investigate the Bogdanov-Takens bifurcation at the nontrivial steady-state solution e = (u * , v * ) T of (1.4) under the condition (H1) and the following assumption That is, the Bogdanov-Takens bifurcation is a bifurcation in a two-parameter family of system (1.4) at which e = (u * , v * ) T has a zero eigenvalue of geometric multiplicity one and algebraic where H(z 1 , z 2 , y) = F(e + Φz + y, α 1 , β) − Ψ, F(e + Φz + y, α 1 , β) .
According to [11], the normal form of Bogdanov-Takens bifurcation under conditions (H1) and (H4) is given by where For convenience, we denotẽ It is easy to see that detG > 0. Therefore, we have the following conclusion.

Conclusions and numerical simulations
In this paper, we have shown that all solutions of system (1.4) exist globally and are uniformly bounded if β satisfies (2.4). But we don't know whether the solution of system (1.4) can blow up in a finite time or exists globally if (2.4) does not hold. This is a problem filled with challenge. Next, this paper presents the existence of the non-constant positive steady states of system (1.4). In view of Theorems 4.5 and 4.6, we see that system (1.4) has a nonconstant positive steady state if either the diffusion coefficient d 1 is small or the cross-diffusion coefficient β is large. This implies the predator and prey species may coexist in the interacting habit nonuniformly if the predator disperses quickly from a high density of prey to a low density one, or the prey move slowly from a higher to a lower concentration region. Our theoretical analysis shows that the cross-diffusion phenomenon has the potential to play an important role in the coexistence information. From the biological point of view, our analysis gives a theoretical support for studying coexistence phenomena of reaction-diffusion systems with cross-diffusion. Sections 6 and 7 show that system (1.4) is capable of producing much more abundant dynamics than the corresponding ODEs. For example, system (1.4) may have multiple bifurcation under certain conditions, and both Hopf bifurcation and homoclinic bifurcation are possible. According to Section 7, we know that the ODEs associated with (1.4) (i.e., with d 1 = d 2 = β = 0) cannot show Bogdanov-Takens singularity, but system (1.4) can show Bogdanov-Takens singularity (see Theorem 7.1); this indicates that diffusion plays a fundamental role in producing a rich dynamics and even Bogdanov-Takens bifurcation phenomena. Meanwhile, the existence and properties of the spatially nonhomogeneous Hopf bifurcation of system (1.4) (i.e., λ i = λ 0 ) are established in Theorem 6.1, and Corollary 6.2 is devoted to spatially homogeneous Hopf bifurcation of system (1.4) (i.e., λ i = λ 0 ).
Finally, we present some numerical simulations to support and supplement our analytic results given in the previous sections. For the spatially homogeneous model (1.4), it follows from Corollary 6.2 that e is locally asymptotically stable if a > b, b < θ(1+θ 2 ) θ 2 −1 and one of conditions (A1)-(A3) is satisfied, and is unstable if θ > 1 and b > θ(1+θ 2 ) θ 2 −1 . In addition, when b passes through θ(1+θ 2 ) θ 2 −1 from the left of θ(1+θ 2 ) θ 2 −1 , e will lose its stability and a family of periodic solutions bifurcate from the interior equilibrium e. It also follows from Corollary 6.2 that the direction of Hopf bifurcation is forward and the bifurcating periodic solutions are asymptotically stable if 2 − 3θ 2 + 6θ 4 − θ 6 < 0. For system (1.4) with Ω = (0, 2π) and initial values u(x, t) = cos 2 x π and v(x, t) = cos 2 x π , if we fix θ = 2.4, then the critical point θ(1+θ 2 ) θ 2 −1 = 3.4084 and 2 − 3θ 2 + 6θ 4 − θ 6 = −7.3174 < 0. Next, we can choose the following three sets of parameter values, which satisfy conditions (A1)-(A3) respectively: Obviously, if the values of a and b are fixed, the mathematical phenomena described by the above three sets of parameters are quite similar. Without loss of generality, we illustrate our analytical results by numerical simulations only under the condition (P1). If a = 5.55 and b = 3.15 < 3.4084, then the positive constant solution e is locally asymptotically stable (see Figure 8.1). Choose a = 5.81, b = 3.41 > 3.4084, then we see that a limit cycle arises out of Hopf bifurcation around e (see Figure 8.2). Lemma 5.2 tells us the positive constant solution e is globally asymptotically stable if a > b, 2b(a + c) < 1 and β ≤ 2 Here, for system (1.4) with Ω = (0, 2π), we choose d 1 = d 2 = 1, a = 1.5, b = 0.15 β = 2.38, and initial values u(x, t) = cos 2 x π , v(x, t) = cos 2 x π , then c = a = 1.5, β < 2 c d 1 d 2 θ b(1+θ 2 ) = 2.3809 and 2b(a + c) = 0.9. Thus, as depicted in Figure 8.3, the positive constant solution e is globally asymptotically stable. Nevertheless, we do not know whether the conclusion of Lemma 5.2 holds true if the value b does not satisfy 2b(a + c) < 1. Therefore, we just find a sufficient condition ensuring the global asymptotical stability of the positive steady-state solution e. However, when a > b, we can get the critical value b * of the parameter b by fixing the value of θ. It follows from Corollary 6.2 that the positive steady-state solution e is locally asymptotically stable if b < b * and will lose its stability when b passes b * from the left of b * .