Turing–Hopf bifurcation in the predator–prey model with cross-diffusion considering two different prey behaviours’ transition

In this paper, we study the Turing–Hopf bifurcation in the predator–prey model with cross-diffusion considering the individual behaviour and herd behaviour transition of prey population subject to homogeneous Neumann boundary condition. Firstly, we study the non-negativity and boundedness of solutions corresponding to the temporal model, spatiotemporal model and the existence and priori boundedness of solutions corresponding to the spatiotemporal model without cross-diffusion. Then by analysing the eigenvalues of characteristic equation associated with the linearized system at the positive constant equilibrium point, we investigate the stability and instability of the corresponding spatiotemporal model. Moreover, by calculating and analysing the normal form on the centre manifold associated with the Turing–Hopf bifurcation, we investigate the dynamical classification near the Turing–Hopf bifurcation point in detail. At last, some numerical simulations results are given to support our analytic results.


Introduction
Since the groundbreaking works of Lotka [1] and Volterra [2], the predator-prey model is used to describe the dynamical interaction between two species and has been widely researched by many scholars in the fields of biology and mathematics [3,4]. It is well known that there are many functional response functions which are used to describe the interactional effect between the prey and predator species. For instance, the Holling I-III types [5,6], ratio-dependent type [7,8], Beddington-DeAngelis type [9], the different types with Allee effect [10][11][12] and so on.
Notice the fact that in natural ecosystem, many species may gather together and form herds to search for food resources or to defend the predators, which means that all members of a group do not interact at a time. This behaviour is often called as the herd behaviour. Recently, the authors in [13] have proposed a more realistic predator-prey model to describe this behaviour which can be written as where u(t) and v(t) stand for the prey and predator densities at time t, respectively. Furthermore, β γ is the death rate of the predator in the absence of prey, and γ is the conversion or consumption rate of prey to predator. The basic idea is that the prey population often gather together in huge herds with the strongest individuals on the border and the weakest will stay in the middle region of the enclosed region. Therefore, the prey which gathers together in the boundary region will be hunted by the predator. That is to say that, the prey population that occupies the outermost positions in the herds will be hunted by the predator. In recent years, many predator-prey models with herd behaviour without diffusion term have been studied. Saha et al. [14] have studied a predator-prey model with herd behaviour and disease in prey incorporating prey refuge, which is an eco-epidemiological model of an infected predator-prey system. They assume that the susceptible prey shows the herd behaviour. The conditions for which the equilibrium point changes their stability and also the conditions for occurring Hopf bifurcation have been analysed. Manna et al. [15] have studied the dynamics of a predator-prey model with Allee effect and herd behaviour, where both predator and prey show herd behaviour. Here, the authors used the linear functional response function. A steady-state analysis has been performed, and some conditions for Hopf bifurcation are derived. Maiti et al. [16] have studied the dynamical behaviours of a predator-prey system with herd behaviour and the Holling type II functional response function, where both the predator and prey show herd behaviours. The positivity, boundedness and stability of equilibrium point are discussed.
From [17][18][19][20], we see that various reaction-diffusion predator-prey models have been extensively studied in the last decades. By considering the fact that in real natural environment, apart from the natural dispersive force of movement of an individual which is usually referred as the self-diffusion, there exists a mutual interference between individuals and is usually referred as the cross-diffusion, see [21][22][23] for details. Recently, cross-diffusion term has appeared in different fields, such as the population dynamics and ecology [24][25][26][27] and chemical reactions [28][29][30][31]. Furthermore, cross-diffusion is taken into consideration in predator-prey model which is used to measure the situation that the prey keeps away from the predator and conversely. Notice also that in recent years, many researchers have shifted from the study of the formation of stationary spatial patterns induced by Turing instability to the study of the formation of spatiotemporal patterns. For example, by combining with the model (1) and considering the self-diffusion and cross-diffusion, Tang et al. [32] have proposed a predator-prey model with herd behaviour and cross-diffusion, and they study the spatiotemporal dynamics near Turing-Hopf bifurcation point of the proposed model. Faria [33] computed the normal forms on centre manifolds (or other invariant manifolds) for partial functional differential equations (PFDEs) near the positive constant equilibrium point. Especially, when a Hopf bifurcation occurs, Faria [33] used the obtained normal forms to study the qualitative behaviours of solutions on those manifolds. In [34], Song et al. have proposed a rigorous procedure for calculating the normal form for the codimensiontwo Turing-Hopf bifurcation in the general reactiondiffusion system and to investigate the dynamical classification near the Turing-Hopf bifurcation point in detail.
Notice the fact that if the herd is too small, it may not be possible for the herd to form an appropriate group defence. In other words, if the herd is too small, the boundary of the herd may consist of the total number of the population. Thus, De Assis et al. [35] have proposed a modified predator-prey model with herd behaviour by considering large prey populations which display herd behaviour, and small prey populations which display individual behaviour. From a mathematical point of view, the interaction term a 1 uv can be used to describe a small number of prey populations, and a 2 √ uv can be used to describe a large prey populations. Thus, the response function of their model is described by a piecewise function where h 1 represents a critical threshold of group size for effective defence. By considering that in real life, a smooth transition between ineffective and effective group defence is expected, at least for some species, and the function g (u) is discontinuous, the authors in [35] have proposed the following response function where a and h 2 are parameters for which the biological interpretation will be explained after some calculations. When u → 0 and u → ∞, the authors in [35] have created a smooth transition function f (u) that approximates g(u) by letting g(u)/ f (u) goes to one. Through simple calculations, the authors in [35] obtained that a = a 2 and h 2 = (a 2 /a 1 ) 2 . Furthermore, it follows that h 1 = (a 2 /a 1 ) 2 from the continuity of g(u), i.e. the critical threshold h 1 of group size for effective defence and the threshold h 2 is consistent. Thus, f (u) can be seen as an approximation of g(u) with a smooth transition for the individual behaviour a 1 u and the group defence a 2 √ u. For convenience, the authors in [35] let h 1 = h 2 = h. Then, the functional response function of their proposed model can be written as Here, by setting a 1 = 1 and a 2 = √ 20, then h 1 = h 2 = h = 20, we plot the graphs of f (u) and g(u) in Fig. 1. The coordinate of point P is (20,20).
Thus, based on the above facts, the authors in [35] have proposed a modified predator-prey model with herd behaviour described by the following ordinary differential equations where r is the intrinsic growth rates of the prey, K is the carrying capacity of the prey, a is the maximum value of prey consumed by per predator per unit time, h is a threshold for the transition between herd grouping and solitary behaviour, m is the death rate of the predator in the absence of prey, and e is the conversion or consumption rate of prey to predator.
To the best of our knowledge, there are no results on spatiotemporal dynamics near Turing-Hopf bifurcation point of the model (2). Therefore, by assuming that the prey and predator populations are in an isolate patch and neglecting the impact of migration, including immigration and emigration, and introducing the spatial diffusion into the model (2), we consider the following modified predator-prey model with herd behaviour and cross-diffusion subject to homogeneous Neumann boundary condition for x ∈ Ω := (0, π) with ∈ R + . That is where u(x, t) and v(x, t) describe the prey and predator densities at the spatial location x and at time t, respectively, the nonnegative constants d 11 and d 22 are the self-diffusion coefficients of the prey and the predator, respectively, and the nonnegative constants d 12 , d 21 are the cross-diffusion coefficients, which describe the respective population fluxes of preys and predators resulting from the presence of the other species, respectively, Δu = ∂ 2 u/∂ x 2 , Δv = ∂ 2 v/∂ x 2 , and n is the outward unit normal vector at the smooth boundary ∂Ω. Moreover, φ(x) and ψ(x) are the initial functions. Furthermore, it is necessary to assume that d 11 d 22 > d 12 d 21 , which indicates that the flux of the respective densities in the spatial domain depends more strongly on their own density than on the other [36]. Especially, we point out that the condition d 11 d 22 > d 12 d 21 can also be obtained in Sect. 3 which is one of the necessary conditions for the occurrence of Turing instability for model (4). In this paper, by calculating the normal form for the codimension-two Turing-Hopf bifurcation in the model (4), we investigate the dynamical classification near the Turing-Hopf bifurcation point. It is necessary to point out the fact that, although we use the method of computing normal form which presented in [34], since we choose d 11 as the bifurcation parameter and the existence of cross-diffusion, the procedure of computing B 11 and B 13 needs to be deduced again, see "Appendix" for details.
The remaining part of this article is organized as follows. In Sect. 2, we show the non-negativity and boundedness of solutions corresponding to the temporal model and the spatiotemporal model, respectively. Furthermore, the existence and priori boundedness of solutions corresponding to the spatiotemporal model without cross-diffusion are also researched. Section 3 is devoted to the stability analysis of the proposed spatiotemporal model, including the stability analysis for the case without self-diffusion and cross-diffusion and the stability analysis for the case with self-diffusion and cross-diffusion. Furthermore, we plot the bifurcation diagram for model (4) in d 11 − δ plane which can be found in Sect. 3. In Sect. 4, we compute the normal form on the centre manifold for Turing-Hopf bifurcation corresponding to the model (4) by using the method in [34]. Some numerical simulations are given to support the theoretical results in Sect. 5. Finally, we conclude this paper in Sect. 6.

Non-negativity and boundedness of solutions corresponding to the temporal model
With a non-dimensionalized change of variables and let then model (3) becomes the following non-dimensiona lized model where h = h/K represents the critical threshold of group size for effective defence in non-dimensional scaling. Since we assumed that there is an observable group defence effect, it is reasonable to take 0 < h < K , and hence 0 < h < 1.
The temporal model associated with the model (4) is For the temporal model (5), we have the following results.
Theorem 1 Suppose that γ > 0, β > 0, 0 < h < 1 and consider the system given by (5) and its trajectories u(t), v(t); if the initial value u(0) ≥ 0, v(0) ≥ 0, then the solutions u(t) and v(t) are always non-negative. Furthermore, for any solution (u(t), v(t)) of system (5), we have Proof In fact, the system (5) is composed of the following three subsystems and Since it is easy to verify that the solution of (7) will approach the origin (0, 0) along the v-axis as t → ∞ and the solution of (refh8) will tend to (1, 0) along the u-axis as t → ∞, we only need to consider the solutions of (6). Assume X (t) = (u(t), v(t)) is a solution of system (6). If it does not remain in the first quadrant, then the solution either hits the u-axis or the v-axis in finite time. Thus, by combining with the analysis results of the subsystems (7) and (refh8), we can obtain that any solution (u(t), v(t)) of system (6) with nonnegative initial value (u(0), v(0)) will remain positive for all t > 0, or will approach either the original or (1, 0) along the axes. Next, the boundedness of solutions is confirmed. Let This concludes the proof.
where |Ω| is the length of the bounded domain Ω.
Proof Denote then by integrating on both sides of each equation on the region Ω in system (4), we have Let , then by combining with Eq. (9) and noticing the homogeneous Neumann boundary defined in system (4), we can obtain from which we obtain that This proof is completed.

The existence and priori boundedness of solutions corresponding to the spatiotemporal model without cross-diffusion
In this section, we give out a sufficient condition for the existence of a positive solution of system (4) without cross-diffusion. Meanwhile, we derive a priori boundedness of the solution. When d 12 = d 21 = 0, system (4) becomes and Ω ⊂ R is a bounded domain with smooth boundary. Then we have the following results: and Ω represents the closure of Ω; are the lower and upper solutions of system (10), respectively. According to Theorem 3.3 in Section 8.3 of Chapter 8 in [37], we know that (10) has a unique globally This completes the proof of conclusion (i).
From the above discussion, we can obtain that 0 One can see that u * (t) → 1 as t → ∞. Thus, for any ε > 0, there exists t 0 > 0 such that u(x, t) ≤ 1 + ε for t > t 0 and x ∈ Ω, which implies that lim By using the homogeneous Neumann boundary condition, we can obtain It is well known that the solution N (t) of In terms of comparison principle and using (11) we obtain that, for This completes the proof of conclusion (ii).

Stability analysis of the corresponding spatiotemporal model
The system (4) has two boundary equilibrium points (0, 0) and (1, 0). Moreover, when 0 < δ < (1/ √ h + 1), the system (4) has a unique positive constant equilib- Defining a real-valued Hilbert space with the inner product defined by It is well known that the eigenvalue problem has eigenvalues λ k = −k 2 / 2 with corresponding normalized eigenfunctions where k ∈ N 0 = N ∪ {0} is often called wave number, and N 0 = N∪{0} is the set of all non-negative integers, N = {1, 2, ...} represents the set of all positive integers. Thus, the eigenfunctions of dΔ on X corresponding to its eigenvalues have the form where a k , b k ∈ R and γ k (x) is defined by Eq. (14). Let be an eigenfunction of dΔ + A corresponding to an eigenvalue λ, that is Then, we have which follows that the eigenvalues of dΔ + A are given by the eigenvalues of L k with k ∈ N 0 = N∪{0}. Notice that the solutions of Eq. (12) can be obtained by using the eigenvalues and eigenvectors of the matrix dΔ+ A. Then by the Fourier expansion which can be seen as a nontrivial solutions of Eq. (12), where and by substituting (15) into (12), we can obtain the following sequence of characteristic equations Δ k , i.e. where By combining with Eqs. (16) and (17), and noticing that the necessary condition for the occurrence of Turing instability is d 11 According to [32], we know that if there exist a nonnegative integer k 1 and a positive integer k 2 = k 1 such that Δ k 1 = 0 has a pair of purely imaginary roots and Δ k 2 = 0 has a simple zero root, and no other roots of (16) has a zero real part, and the corresponding transversality conditions hold, then we call the bifurcation in this case as a Turing-Hopf bifurcation.

Stability analysis for the case without self-diffusion and cross-diffusion
In the case of without self-diffusion or cross-diffusion (d 11 = d 12 = d 21 = d 22 = 0), the original system (4) becomes the following ordinary differential equation Clearly, the original system (4) and (19) have the same positive constant equilibrium point E * (u * , v * ). By a simple linear analysis, we can obtain the following result.

Theorem 4 For the positive constant equilibrium point
we have the following results on its stability.

Stability analysis for the case with self-diffusion and cross-diffusion
In this section, we investigate the effect of crossdiffusion on the positive constant equilibrium point E * (u * , v * ) of the original system (4). We have the following Theorem 5. Then (ii) For (d 11 , δ) = d * 11 , δ * , the equation Δ 0 = 0 has a pair of purely imaginary roots ±iω c and Δ k * = 0 has a simple zero root, and for the system (4), there are no other roots with zero real parts, where Proof The Turing bifurcation curve L k defined by (21) is followed from D k = 0, i.e.
By setting then there exists a k * > k * such that d * Thus, the Hopf bifurcation curve H 0 intersects with the Turing bifurcation curve L k at which is called as the Turing-Hopf bifurcation point. Next, we continue to verify the transversality condition. Taking d 11 as a bifurcation parameter and letting λ(d 11 ) be the root of Eq. (16) near d 11 = d * 11 satisfying λ d * 11 = 0. Differentiating the two sides of the characteristic equation (16) with respect to d 11 , we obtain By combining with Eq. (17), we can obtain Moreover, when δ > δ * , from the discussion in Sect.

Normal form for Turing-Hopf bifurcation in the reaction-diffusion system (4)
Notice that the method of computing normal form of general reaction-diffusion system which presented in [34] can be used to obtain the third-order truncated normal form of system (4) by a slight modification. Denote the Turing-Hopf bifurcation point as d * 11 , δ * . Introduce the perturbation parameters μ 1 and μ 2 by setting d 11 = d * 11 + μ 1 and δ = δ * + μ 2 such that (μ 1 , μ 2 ) = (0, 0) is the Turing-Hopf bifurcation point in the perturbation plane of μ 1 and μ 2 . Then the system (4) becomes Making the change of variables by the translation u = u − u * and v = v − v * , and dropping the bars, then the system (28) is transformed into For the system (29), when μ 1 = μ 2 = 0, Δ 0 (λ) = 0 has a pair of purely imaginary roots ±iω c , Δ k * (λ) = 0 has a simple zero root λ = 0 and a negative real root λ = −T k * , and if k = 0, k * , all roots of Δ k (λ) = 0 have negative real parts.
By the real-valued Hilbert space X which is defined in Sect. 3, the system (29) can be written as the following abstract ordinary differential equation (ODE), that Here, a 11 (δ * + μ 2 ), a 21 (δ * + μ 2 ) indicating a 11 and a 21 are dependent on δ * + μ 2 , respectively. Moreover, we have For the formal Taylor expansions of L and by Eq. (31), we can obtain where I 2 is a 2 × 2 identity matrix, then by a straightforward calculation, we can obtain that ξ 0 ∈ C 2 and ξ k * ∈ R 2 are the eigenvectors associated with the eigenvalues iω c and 0, respectively, and η 0 ∈ C 2 and η k * ∈ R 2 are the corresponding adjoint eigenvectors, where such that Here, col(.) represents the column vector. Furthermore, for vectors ϕ, ψ ∈ R 2 , we define their scalar product as ψ T , ϕ = ψ T ϕ. Notice that the phase space X can be decomposed as where dim P = 3 and π : X → P is the projection defined by According to (33), U ∈ X can be decomposed as where w = col (w 1 , w 2 ) and By letting Equation (34) can be rewritten as For the simplicity of notation, we denote and if let L U := dΔU + L 0 (U ) and denote by L 1 the restriction of L to Q, then system (30) is equivalent to a system of abstract ordinary differential equations (ODEs) in R 3 ×Q, with finite-and infinite-dimensional variables also separated in the linear term. That is, H (z, w, μ), According to [34], by a recursive transformation, the authors obtain that the normal form for Turing-Hopf bifurcation in system (30) reads aṡ The normal form (35) can be written in real coordinates v through the change of variables 3 , and then changing to cylindrical coordinates by v 1 = ρ cos Θ, v 2 = ρ sin Θ, v 3 = ς , we obtain, truncating at third-order terms and removing the azimuthal terṁ where According to Section 8.6 in [38], we know that the third-order truncated normal form (36) is exactly the same to the third-order truncated normal form of the four-dimensional smooth system depending on two parameters with Hopf-Hopf bifurcation. For the socalled simple case [38], i.e. κ 11 κ 22 > 0, the dynamics of system (4) near the bifurcation value is topologically equivalent to that of normal form (36). However, for the "difficult" cases, the original system (4) is never topologically equivalent to the truncated normal form (36). In general, for the "difficult" cases, in order to obtain the whole dynamics of the original (4), five-order or higher-order normal form needs to be calculated.
With the help of MATLAB software, the explicit values of the coefficients B 11 , B 21 , B 13 , B 23 , B 210 , B 102 ,  B 111 and B 003 can be obtained for the fixed parameters. Notice that the formulas which are used to calculate the above coefficients are rather complicated and we leave them in "Appendix".

Numerical simulations
In this section, the third-order truncated normal form is obtained under the given parameters and some numerical simulations are made to support the results of our theoretical analysis. More precisely, some numerical simulations about the temporal patterns, spatial patterns and spatiotemporal patterns are given.
Let Ω = (0, π) and γ = 3, h = 0.1, d 12 = 1, d 21 = 10, d 22 = 15, then the system (4) at least undergoes Turing-Hopf bifurcation at the point (1.1467, 0.4404). By using the above-given parameters, the normal form truncated to the third-order terms isρ where μ 1 and μ 2 are perturbation parameters for the Turing-Hopf bifurcation point (1.1467, 0.4404). Notice that ρ > 0 and ς is arbitrary real number. System (37) has a zero equilibrium point A 0 (0, 0) for any μ 1 , μ 2 ∈ R, three boundary equilibrium points for 0.9290μ 1 − 2.5716μ 2 > 0, and two interior equilibrium points Table 1 The twelve unfoldings of (36) Table 2 The corresponding relationship between the equilibrium points of (37) and the solutions of original system (4) Equilibrium points of (37) Solutions of the original system (4) A 0 Positive constant steady state Spatially homogeneous periodic solution Two spatially inhomogeneous steady states with cos x-like shape in space Two spatially inhomogeneous periodic solutions with cos x-like shape in space These four lines divide the μ 1 − μ 2 parameter plane into six regions marked as D j , j = 1, 2, . . . , 6.
Based on Section 7.5 in [39], by the different signs of d, b, c, d − bc in Table 1, the system (36) has twelve distinct types of unfoldings, which are twelve essentially distinct types of phase portraits and bifurcation diagrams. More precisely, for the system (37), we have d = +1, b = 1.3593 > 0, c = 1.3860 > 0, d −bc = −0.8840 < 0. That is, the unfolding of the planner system (37) corresponding to the Case Ib in Table 1. Thus, the phase portraits and bifurcation diagrams corresponding to Case Ib can be given out, see Fig. 7.5.2 in Section 7.5 in [39] for details.
The dynamics of the original reaction-diffusion system (4) can be determined by the third-order truncated normal form (37) near the neighbourhood of the Turing-Hopf bifurcation point. According to [34], the  (1.1467, 0.4404) corresponding relationships between the equilibrium points of plane system (37) and the solutions of original system (4) are shown in Table 2.
Furthermore, by the defined critical bifurcation lines, the bifurcation diagram in the μ 1 − μ 2 parameter plane is shown in Fig. 3. The linearized equation of the system (37) at each equilibrium point is d dt  where More precisely, the coefficient matrices of linearized equation (38) at equilibrium points A 0 , respectively, where By combining with the bifurcation diagram in Fig.  3 and the linear stability theory, we can analyse the sign of the eigenvalues corresponding to the characteristic equations; thus, the stability and instability of each equilibrium point in regions D 1 -D 6 can be obtained.
In region D 1 , the third-order truncated normal form (37) has only one equilibrium point A 0 and it is asymptotically stable. This implies that the positive constant steady state E * of the original system (4) is asymptotically stable, as shown in Fig. 4 for (μ 1 , μ 2 ) = (0.60, 0.25) and the initial value u(x, 0) = 0.5615 − 0.25 cos(x), v(x, 0) = 0.7386 + 0.25 cos x.
In region D 3 , the third-order truncated normal form (37) has four equilibrium points: A 0 , A 1 , A + 2 and A − 2 . The equilibrium points A 0 , A + 2 and A − 2 are unstable and the equilibrium point A 1 is asymptotically stable. Thus, the original system (4) has an unstable positive constant steady state, two unstable spatially inhomogeneous steady states like cos x-shape in space, and a stable spatially homogeneous periodic solution. By choosing (μ 1 , μ 2 ) = (−0.16, −0.08) and the initial value u(x, 0) = 0.1961 + 0.2 cos x, v(x, 0) = 0.4729+0.2 cos x, the dynamics of the original system (4) evolves from the spatially inhomogeneous steady states to the spatially homogeneous periodic solution, as shown in Fig. 6.
In region D 4 , the third-order truncated normal form (37) has six equilibrium points: A 0 , A 1 , A ± 2 and A ± 3 . The equilibrium point A 0 is unstable, the equilibrium points A 1 , A ± 2 and A ±  to stable spatially homogeneous periodic solution, as shown in Fig. 8.
In region D 5 , the third-order truncated normal form has four equilibrium points: A 0 , A 1 , A + 2 and A − 2 . The equilibrium point A 0 is unstable and the equilibrium points A 1 , A + 2 and A − 2 are asymptotically stable. Thus, the original system (4) has an unstable positive constant steady state, a stable spatially homogeneous periodic solution and two stable spatially inhomogeneous steady states like cos x-shape in space. Taking the parameter (μ 1 , μ 2 ) = (0.17, −0.03) and the initial value u(x, 0) = 0.2389 − 0.15 cos x, v(x, 0) = 0.5455 − 0.15 cos x close to the unstable positive constant steady state, the dynamics of the original system (4) evolves from unstable positive constant steady state to stable spatially homogeneous periodic solution, as shown in Fig. 9a and b. Taking the parameter (μ 1 , μ 2 ) = (0.31, −0.05) and the initial value u(x, 0) = 0.2213 − 0.2 cos x, v(x, 0) = 0.5169 + 0.2 cos x close to the unstable positive constant steady state, the dynamics of the original system (4) evolves from unstable positive constant steady state to stable spatially inhomogeneous steady states, as shown in Fig.  9c and d.
In region D 6 , the third-order truncated normal form (37) has three equilibrium points: A 0 , A + 2 and A − 2 . The equilibrium point A 0 is unstable and the equilibrium points A + 2 and A − 2 are asymptotically stable. Thus, the original system (4) has an unstable positive con-stant steady state and two stable non-constant steady states like cos x-shape in space. For the fixed parameter (μ 1 , μ 2 ) = (0.75, 0.1675) and choosing different initial values, the original system (4) can converge to one of these two stable non-constant steady states, as shown in Fig. 10a

Conclusion and discussion
In this paper, a predator-prey model with crossdiffusion considering the prey individual behaviour and herd behaviour transition with homogeneous Neumann boundary condition is investigated. We first show that the non-negativity and boundedness of solutions corresponding to the model without self-diffusion and cross-diffusion and the model with self-diffusion and cross-diffusion. Then we show the existence and priori boundedness of solutions corresponding to the spatiotemporal model without cross-diffusion. In order to classify the possible dynamical classification near the Turing-Hopf bifurcation point, by using the method of computing the normal form presented in [34], the third-order truncated normal form (37) is given. By the obtained third-order truncated normal form (37), we obtain a zero equilibrium point corresponding to the positive constant steady state of the original sys- When (μ 1 , μ 2 ) = (0.31, −0.05) lies in region D 5 and let d 11 = 0.9, δ = 0.4004, the positive constant steady state E * (0.2213, 0.5169) is unstable and there is a heteroclinic orbit connecting the unstable positive constant steady state to stable spatially inhomogeneous steady states. The initial value is u(x, 0) = 0.2213 − 0.2 cos x, v(x, 0) = 0.5169 + 0.2 cos x tem (4), three boundary equilibrium points A 1 and A ± 2 in which A 1 corresponding to the spatially homogeneous periodic solution of the original system (4) and A ± 2 corresponding to the two spatially inhomogeneous steady states with cos x-like shape in space of the original system (4). Furthermore, two interior equilibrium points A ± 3 is also obtained corresponding to the two spatially inhomogeneous periodic solutions with cos xlike shape in space of the original system (4). Moreover, we obtain four inequality which is used to ensure the existence of these different types of equilibrium points. Notice that the four inequality can be seen as the four critical bifurcation lines; thus, by the defined critical bifurcation lines, the bifurcation diagram in the μ 1 -μ 2 parameters plane which includes six different regions is shown in Fig. 3.
By the numerical simulations, the rich dynamics such as the positive constant steady state, the spatially homogeneous periodic solution, spatially inhomogeneous steady states and spatially inhomogeneous periodic solutions have been found, which can be seen in Figs. 4, 5, 6, 7, 8, 9 and 10. Especially, we would like to mention that the interaction between the Turing bifurcation curve and the Hopf bifurcation curve may leads to the emergence of the spatially inhomogeneous periodic solutions, see Fig. 7a and b for details. We have to point out the fact that the method of computing normal form developing in this paper can also be used to the case without cross-diffusion by a slight modification. Furthermore, in order to the simplicity of notation for further normal form computation and for the convenience of carrying out the numerical simulations, we let Ω := (0, π). However, the spatial Ω can also be taken (0, π) with ∈ R + . Notice also that the general open interval ( a, b) can be transformed to (0, π) by a translation and rescaling.
In Sect. 1, we have pointed out that the procedure of computing B 11 and B 13 needs to be deduced again. In the following, the detailed derivation process of B 11 and B 13 is given. Following Section 3.1.1 in [34], by considering the formal Taylor expansion where F j (Φz x + w, μ), j ≥ 2 is the j-th Fréchet derivative of F(Φz x + w, μ). For simplicity, we set , α ∈ C.