Turing-Hopf bifurcation analysis in a diffusive Gierer-Meinhardt model

The reaction-diffusion Gierer-Meinhardt system in one dimensional bounded domain is considered in the present paper. The Hopf bifurcation is investigated, which is found to be degenerate. With the aid of Maple, the normal form associated with the degenerate Hopf bifurcation is obtained to determinate the existence of Bautin bifurcation. We get the universal unfolding for the Bautin bifurcation so that we can identify the stability of periodic solutions. Then, the existence of the codimension-two Turing-Hopf bifurcation is further investigated. To research the spatiotemporal dynamics of the model near the Turing-Hopf bifurcation point, the method of the multiple time scale analysis is adopted to derive the amplitude equations. It is noted that the Gierer-Meinhardt model may show the spatial, temporal or the spatiotemporal patterns, such as the nonconstant steady state, spatially homogeneous periodic solutions and the spatially inhomogeneous periodic solutions. Finally, some numerical simulations are presented to demonstrate the applicability of the theoretical results.


Introduction
As we know, the development of an organism is a complicated phenomenon involving a set of more elementary processes, one of which in embryology and regeneration is the formation of a spatial pattern of tissues. For a long time many researchers tried to find the underlying mechanisms and explain the formation of organs. In 1950s, the British mathematician Turing [1] explained the mechanism of pattern formation of chemical or morphogen concentrations, and showed that the coupled reaction diffusion equations could cause the spatial patterns due to the effect of diffusion. In 1972, an activatorinhibitor reaction-diffusion system was proposed by Gierer and Meinhardt [2] during their study of relatively simple molecular mechanisms based on auto-and cross catalysis. In 1974, they obtained some sufficient conditions to ensure the formation of spatial patterns in [3]. Since then, the Gierer-Meinhardt model was frequently utilized to present the formation of morphogenesis. Generally, the Gierer-Meinhardt model is described as follows ∂V(x,t) ∂t (1 .1) where m, n, k and l are non-negative integers. Now consider the Gierer-Meinhardt model with m = l = 2, n = 1 and k = 0. Then, the above model reduces to ∂V(x,t) ∂t where U(x, t) and V(x, t) denote the concentrations of activator and inhibitor, respectively; ρ 0 ρ is the exogenous source of activator; θ 1 and θ 2 are the degradation coefficients of the activator and inhibitor; cρ and c 1 ρ 1 are cross-reaction coefficients; D 1 and D 2 are the diffusion coefficients of activator and inhibitor, respectively; cρ U 2 (x,t) V(x,t) denotes the interaction of two reactants, U 2 (x, t) in the numerator stands for the self-catalysis of activator and V(x, t) in the denominator represents the effect of activator by inhibitor; c 1 ρ 1 U 2 (x, t) is the effect of autocatalysis of activator on inhibitor.
For the Gierer-Meinhardt system (1.2) with the saturating term, Song et al. [4] identified the parameter region where possible Turing instability could happen and obtained the amplitude equations at the critical value of bifurcation through the multiple scale method. Then they found that some patterns, such as spot-like pattern, stripe-like pattern and the coexistence pattern, could appear. Chen et al. [5] worked on the Gierer-Meinhardt system with a saturation in the activator production and investigated the stability of the unique positive constant steady state solution, Hopf bifurcations and steady state bifurcations. They obtained a global bifurcation diagram of non-trivial periodic orbits and steady state solutions. For the Gierer-Meinhardt system (1.2) without the saturating term, some related works have been reported. Wu et al. [6] investigated the effects of diffusion on the stability of equilibrium and the bifurcated limit cycle from Hopf bifurcation. Moreover, with some conditions, the diffusion can drive the Turing instability and those diffusion-driven instabilities would lead to the occurrence of various patterns. In [7], the authors considered the Gierer-Meinhardt model without the saturation term and investigated the Turing instability of the positive equilibrium, the existence of the Turing-Hopf and spatial resonance bifurcation. Ruan [8] considered the Gierer-Meinhard model of morphogenesis, showing that the homogeneous equilibrium solution and the homogeneous periodic solution could turn to be diffusively unstable, if the diffusion coefficients of the two substances were chosen suitably. Liu et al. [9] studied the Hopf bifurcation and steady state bifurcation of a reaction-diffusion Gierer-Meinhardt model of morphogenesis subject to Dirichlet fixed boundary conditions in a one-dimensional spatial domain, obtained the spatially nonhomogenous periodic orbits and nonconstant positive steady state solutions.
In [10], the activator-inhibitor model with different sources was considered with the exogenous sources of activator being zero, that is to say, this activation-inhibition model is formed in a sealed circumstance. By the non-dimension scaling transformation [10], then the Gierer-Meinhardt system could be simplified into Based on multiple scale analysis and the obtained amplitude equations, it was showed that the model could admit spotted patterns, stripe patterns and the coexistence of spotted and stripe patterns. For the model (1.3) with a saturation in the activator production, Chen et al. [11] proved that the unique constant steady state solution was globally attractive by employing an upper and lower solution method. They showed that the parameter area for the formation of spatiotemporal patterns would be limited. Yang et al. [12] investigated the conditions of the stability of the unique positive equilibrium in a diffusive activator-inhibitor model, the existence of Hopf and steady state bifurcations. They obtained the normal forms corresponding to the Hopf bifurcation and steady state bifurcation. Results about other bifurcations for the Gierer-Meinhardt model with different saturation terms were also reported in [13,14]. However, the direction of Hopf bifurcation for the model (1.3) in the references was not reported in detail and only the spatial patterns induced by Turing instability were presented. So it is worth considering whether the exogenous sources of activators will cause the spatiotemporal patterns to happen. We will find that the Hopf bifurcation in (1.3) is degenerate, and the bifurcation analysis could still be processed.
The normal form theory and center manifold reduction [15] are commonly used to determine the direction of the Hopf bifurcation and the stability of the bifurcated periodic solution, see, for example, [16,17]. Since the Hopf bifurcation is degenerate in the present paper, we will use the formal procedure in Maple to get the unfolding of Hopf bifurcation. The normal form associated with the degenerated Hopf bifurcation could be established, then by analyzing the obtained normal form, one can analyze the degenerate Hopf bifurcation, and the stability of bifurcated periodic solutions. Moreover, to obtain richer and more complex spatiotemporal dynamical behaviors of the Gierer-Meinhardt system, we will find that the spatiotemporal dynamics will be exhibited in the model near the Turing-Hopf bifurcation point by the multiple time scale analysis.
To this end, a diffusion activator-inhibitor model with Neumann boundary conditions and positive initial conditions can be formulated as follows where D 1 and D 2 are diffusion coefficients of activator u(x, t) and inhibitor v(x, t) at position x and time t, respectively; ∆ = ∂ 2 ∂x 2 denotes the Laplacian operator; Ω = (0, π) is a bounded domain in R with the boundary ∂Ω. Parameters D 1 , D 2 , g and e are positive. In the present paper, we will use the standard multiple scale analysis [18] to investigate the dynamical behaviors of Gierer-Meinhardt model.
The rest of this paper is organized as follows. In Section 2, the ordinary differential equations (ODEs) of model (1.4) will be first focused on. The normal form associated with the Hopf bifurcation will be given by using the symbolic language Maple to determine the existence of the Bautin bifurcation. Moreover, we will get the universal unfolding of the ODEs of model (1.4) for the Bautin bifurcation, so that we could identify the bifurcation. Then for the partial differential equations (PDEs) of model (1.4), by choosing the parameters g and D 2 as the Turing and Hopf bifurcation parameters, the codimension-two Turing-Hopf bifurcation point in the (D 2 , g) plane will be located. The existence of the Turing-Hopf bifurcation will be discussed in detail. In Section 3, the multiple time scale analysis is adopted to get the amplitude equations near the Turing-Hopf bifurcation point. In Section 4, numerical simulations are carried out to illustrate the theoretical analysis. Finally, some discussions and conclusions are made in Section 5.

Analysis of the corresponding ODE system
In this subsection, we consider the ODE system with respect to (1.4) It is not difficult to find that system (2.1) has one organism-free equilibrium: E 0 := (0, 0) and a positive equilibrium E * := (u * , v * ) = ( g e , g e ), normally we are interested in the positive equilibrium. The Jacobian matrix of (2.1) evaluated at the positive equilibrium E * is Therefore the corresponding characteristic equation for J(E * ) is the eigenvalues are where tr(J(E * )) = 1 − g, detJ(E * ) = g > 0.
Theorem 2.1. The positive equilibrium E * of system (2.1) is locally asymptotically stable if and is unstable if Moreover, the Hopf bifurcation will occur at the point E * in system (2.1) when g = 1.
Proof. If g > 1 holds, the eigenvalues have negative real parts, thus the equilibrium E * of system (2.1) is locally asymptotically stable; if g < 1 holds, then the eigenvalues have positive real parts, thus the equilibrium E * of system (2.1) is unstable; if g = 1 holds, then tr(J(E * )) = 1−g = 0, so the eigenvalues are a pair of conjugate pure imaginary roots λ 1,2 = ±iω * = ±i √ g = ±i. Note that then we choose g as the Hopf bifurcation parameter, the Hopf bifurcation may occur in system (2.1) on the critical line L 1 : g = 1.
In the following, we will employ the method in [15] to obtain the direction of Hopf bifurcation. Put u = u − g e andv = v − g e into system (2.1), then we have when g = 1, we have λ 1,2 = ±iω * = ±i √ g, and (ω * ) 2 = g = 1 > 0. Let T be the matrix that transforms the Jacobian matrix into a standard form Under the transformation , Then the stability of Hopf bifurcation in system (2.1) at E * is determined by the sign of the following quantity where all partial derivatives are calculated at the bifurcation point (u, v, g) = (0, 0, 1).
By calculation, we find that This implies that Hopf bifurcation is degenerate [19,20]. Next, the technique in [21] will be employed to analyze the degenerate bifurcation. First the normal form of system (2.1) under the Hopf bifurcation condition should be presented. In [21], the author mainly introduced how to use a perturbation technique to find a unique normal form for a given set of differential equations, the procedure is formulated as the Maple code. Now assume g = 1, that is to say, the Hopf bifurcation occurs in model (2.1), then system (2.1) becomes Since the linear part of the system is not the standard form, therefore from the above transformation, then we have Before using Maple, we need to determine the data in the input file. The input file contains the data of the input functions f i , the number of the real eigenvalues M1, the number of the pairs of complex conjugate eigenvalues M2, and the order of the normal forms. Therefore, for system (2.6), M1 = 0 is the number of the non-zero real eigenvalues, and M2 = 0 is the number of the pairs of complex conjugate eigenvalues. The system (2.4) is dimension N = 2 + 0 + 2 × 0 = 2, indicating that the system has only a pair of purely imaginary eigenvalues(±iω * = ±i). Order = 5 indicates that the procedure will be executed up to sixth order approximations. The input file is represented as that in Table 1, which is shown below.
From [22], we have the universal unfolding of system (2.1) for Bautin bifurcation is Therefore, from the universal unfolding (2.13) we can study the distribution of limit cycles in system (2.1).
Remark 2.1. Since σ = 0 at the parameter value g = 1, system (2.1) undergoes the Buatin bifurcation. The bifurcation diagram of Bautin bifurcation can be found in some literatures [23,24]. System (2.1) always has an equilibrium, may have no, one or two limit cycles. Here, one of two limit cycles is stable and another is unstable.
Now the numerical simulations are carried out to illustrate the Hopf bifurcation. The results are shown in Figure 1A and 1B. Take parameters as g = 1.5, e = 1 in Figure 1A, parameters are g = 1, e = 1 in Figure 1B.  Figure 1. The equilibrium E * is stable and the stable periodic solution occur in system (2.1).
Here (A) for the stable equilibrium and (B) for Hopf bifurcation.

Existence of the Turing-Hopf bifurcation
In this subsection, we will study the diffusive effects on the stability of the positive equilibrium and show the existence of the Turing-Hopf bifurcation. From the above discussion, system (1.4) has the only positive equilibrium E * = ( g e , g e ). So the characteristic equation for the positive equilibrium E * is ∆ n (λ) = λ 2 + T n λ + J n = 0, n ∈ N 0 {0, 1, 2, 3, · · · }, (2.14) where so the eigenvalues are In what follows, we mainly research the diffusion-driven instability and the Turing-Hopf bifurcation analysis of system (1.4).
Theorem 2.2. Assume that D 1 , D 2 > 0, there are a diffusion-driven Turing instability and the Turing-Hopf bifurcation of system (1.4) at P * = (D n * 2 , g * ), where D n 2 , g n (D 2 ) are defined by (2.15), (2.16) and n * is the critical wave number for the Turing-Hopf bifurcation defined below. Then we have the following stability and instability results for system (1.4).
(1) The positive equilibrium E * is locally asymptotically stable for (D 2 , g) ∈ R 1 and unstable for (D 2 , g) ∈ R 2 . Here R 1 and R 2 are defined by ,
(3) System (1.4) undergoes the Turing-Hopf bifurcation at the point (D n * 2 , g * ). Proof. From Theorem 2.1, the Hopf bifurcation occurs at the line defined by L 1 : g = 1. To make sure the occurrence of Turing bifurcation, it is necessary that g > 1 and there exists some n ∈ N 0 /{0} such that J n < 0, which means g < g n (D 2 ). Thus the critical parameter value of the Turing instability is defined by J n = 0, which leads to the Turing critical line It is noted that the line L 2 : g = g n (D 2 ) is linear and its slope k(n) is given by To make the lines L 1 and L 2 intersect, it is necessary that k(n) > 0 must be satisfied. That is to say, 1 − D 1 n 2 > 0, n ∈ 0, 1 D 1 . If n ∈ 0, 1 D 1 exists which means 0 < D 1 < 1, then there is the Turing instability. Otherwise there is no Turing instability. Therefore when n ∈ 0, 1 D 1 and 0 < D 1 < 1, the Turing instability exists.
The critical lines of Hopf and Turing bifurcations intersect at the point P and when n ∈ 0, 1 D 1 , 0 < D 1 < 1, the intersecting point P * exists. In the following, we are going to choose n which minimizes k(n) as the critical wave number for the Turing-Hopf bifurcation. Now we provide some instructions about the monotone intervals of k(n). Differentiating k(n) with respect to n, we have Then we get the following instructions.
is the integer part function. Let n 1 = 1 and n 2 = 1 D 1 . Next, under Case (1), that is 0 < D 1 ≤ √ 2−1. We will discuss the value of the critical wave number n * in the following cases.
1 if k(n 1 ) ≤ k(n 2 ), then we choose n * = n 1 = 1; 2 if k(n 1 ) > k(n 2 ), then we choose n * = n 2 = 1 D 1 . Under Case (2), we know that then when n ∈ 0, 1 D 1 , k(n) is monotonically decreasing, then we choose n * = n 2 as the critical wave number. In the rest, taking g as a parameter we consider the transversality condition below Through the above analysis, when (D 2 , g) = (D n * 2 , g * ), ∆ 0 (λ) = 0 has a pair of purely imaginary roots ±i √ J 0 , and ∆ n * (λ) = 0 has a root λ = 0 with a negative root λ = −T n * (g). For the other wave number n 0, n * , all roots of ∆ n (λ) = 0 have negative real parts. Together with the transversality conditions (2.2) and (2.19), then we can arrive at the conclusion that system (1.4) undergoes the Turing-Hopf bifurcation at (D 2 , g) = (D n * 2 , g * ), then the proof is completed.

Amplitude equations for the Turing-Hopf bifurcation
In this subsection, the spatiotemporal dynamics near the Turing-Hopf bifurcation point (D 2 , g) = (D n * 2 , g * ) will be explored by using the multiple time scale method [18]. From Theorem 2.2, for given parameters g, D 1 and D 2 , we can obtain the critical wave number n * in terms of (2.17) and (2.19). In the following, the process to get the the amplitude equations will be divided into three steps.
Step 1. Putting u = u − u * , v = v − v * into system (1.4), then we can obtain (3.1) Let U = (u, v) T , then the system (3.1) can be transformed into and (3.5) Step 2. The solution of system (3.2) near the Turing-Hopf bifurcation point (D n * 2 , g * ) can be represented by where W 1 is the amplitude of Hopf mode, W 2 is the amplitude of Turing mode and n * represents the critical wave number. To obtain the solution, we need to express it by the original parameters. First of all, we expand the solution U in the following forms with a small parameter ε: In the same way, expand the Turing-Hopf bifurcation parameters D 2 , g and the nonlinear term N in the small neighborhood with the small parameter ε. Then we have Next, we expand the time scale of system (3.2). Let T 0 = t, T 2 = ε 2 t, and T 0 , T 2 can be regard as independent variables, then the form of the derivative with respect to time t can be represented as follows After that, substitute (3.3)-(3.5) and (3.7)-(3.9) into (3.2), then expanding equation with respect to different orders of ε, comparing with the coefficients of ε on the both sides of equation, we get collecting the coefficients of ε 2 on the both sides of equation, we have comparing with the coefficients of ε 3 on the both sides of equation, we obtain From Eq (3.10), the general solution can be represented as a 1 e iω * T 0 + W 2 (T 2 )a 2 e in * x + c.c., (3.13) where the c.c. represents conjugate terms. Substituting (3.13) into (3.10), one has For Eq (3.11), the right-hand side of (3.11) is h 2 , which has the formulas of u 2 1 , u 1 v 1 and v 2 1 . Since these formulas contain and the conjugate terms of the above terms. Thus the solution of Eq (3.11) can be expressed as (3.14) where b 1 , b 2 , b 3 , b 4 , b 5 and b 6 are undetermined coefficients. Putting (3.14) into (3.11), comparing with the coefficients of similar terms W 2 1 e 2iω * T 0 , |W 1 | 2 , W 2 2 e 2in * x , |W 2 | 2 , W 1 W 2 e iω * T 0 +in * x , W 1W2 e iω * T 0 −in * x , which are on both sides of the equation, then we get As for Eq (3.12), putting (3.13) and (3.14) into Eq (3.12), then we obtain where NST stands for the non-secular terms and what is more, we have where c 3 =2e(b 11 + 2Re{b 21 }), In order to get the unique solution of Eq (3.12), the formula on the right-hand side of (3.15) must satisfy the Fredholm solvability conditions, that is to say, the following equations must be satisfied and the inner product has the form (p * j ) T u j = 0, j = 1, 2. Thus, from Eq (3.16), we can get the amplitude equations as follows Step 3. Putting the transformations W 1 = ρ 1 e iθ , W 2 = ρ 2 into system (3.17), we obtain the equivalent amplitude equations of the Turing-Hopf bifurcation in the real coordinates: = Re{ϕ 1 }ρ 1 + Re{ϕ 2 }ρ 3 1 + Re{ϕ 3 }ρ 1 ρ 2 2 , ρ 2 = Re{ψ 1 }ρ 2 + Re{ψ 2 }ρ 3 2 + Re{ψ 3 }ρ 2 ρ 2 1 , θ = Im{ϕ 1 } + Im{ϕ 2 }ρ 2 1 + Im{ϕ 3 }ρ 2 2 . (3.18) where Re{·} and Im{·} represent the real part and the imaginary part of ·, respectively. Then after we truncate the third order terms and remove the azimuth terms, the amplitude equations of the Turing-Hopf bifurcation become (3.19)

Numerical simulations
In this section, some numerical simulations are carried out to verify the above theoretical analysis. In the pervious sections, the amplitude equations of system (1.4) are obtained near the Turing-Hopf bifurcation. However, what kind of dynamical behaviors will appear in the system near the Turing-Hopf point (D n * 2 , g * )? Then what we are going to focus on is what kinds of dynamical patterns of system (1.4) near the bifurcation point will occur. We mainly study the amplitude Eq (3.19). First of all, we set D 1 = 0.6, e = 1, g = 1, it follows from Theorem 2.2 that n * = 1, D n * 2 = 4 and L 1 : g = 1, L 2 : g = g n (D 2 ) = 0.25D 2 , D 2 > 0.
The straight line L 1 intersects with L 2 at the point P * (4, 1) in the D 2 − g plane, system (1.4) will undergoes the Turing-Hopf bifurcation near the positive equilibrium E * at the point P * . The stable region for the positive equilibrium E * in the D 2 − g plane is shown in Figure 2. Figure 2. Bifurcation diagram near the Turing-Hopf bifurcation point P * (D n * 2 , g * ) in the plane of d 2 − g 2 .
According to the calculation procedure in Section 2, by using numerical simulation, we have Then, the amplitude equations truncated to the third order terms are where d 2 and g 2 are perturbation parameters for the Turing-Hopf point P * (D n * 2 , g * ). Next, we will calculate the equilibria of system (4.1). Furthermore, note that ρ 1 > 0, and ρ 2 is an arbitrary real number. Hence, the different constant solutions of system (4.1) are: , f or all d 2 , g 2 , three boundary equilibria : Therefore, the critical bifurcation lines are performed as follows T 1 : g 2 = 0; T 2 : g 2 = 0.25d 2 ; T 3 : g 2 = 0.2439d 2 , d 2 < 0; We know that the boundary equilibria Q 1 and Q ± 2 are bifurcated from the origin on the critical line T 1 and T 2 , respectively. The two internal equilibrium points Q ± 3 are bifurcated from the boundary equilibria Q 1 and Q ± 2 on the critical line T 3 and T 4 , respectively. These four straight lines T 1 , T 2 , T 3 and T 4 divide the parameter plane d 2 − g 2 into six regions marked by R 1 − R 6 , and the parameter plane d 2 − g 2 can be shown in Figure 3. It is worth noting that the equilibria Q 0 , Q 1 , Q ± 2 and Q ± 3 of the amplitude equation (4.1) correspond to the positive constant solution, the spatially homogeneous periodic solution, the nonconstant steady state and spatially inhomogeneous periodic solution of system (1.4). Therefore, the dynamics of the system (1.4) near the Turing-Hopf bifurcation point p * (D n * 2 , g * ) in the parameter plane d 2 − g 2 can be identified on the basis of the dynamics of the normal form system (4.1). In the following, we shall introduce the dynamics of system (1.4) when the values of the parameters are chosen differently in these six regions.
In region R 5 , system (4.1) also has four equilibria: Q 0 , Q 1 and Q ± 2 . The equilibrium Q ± 2 is stable, and the equilibria Q 0 and Q 1 are unstable. That means the original model   In region R 6 , system (4.1) also has four equilibria: Q 0 , Q 1 and Q ± 3 . The equilibrium Q 0 and Q 1 are stable, and Q ± 3 is unstable. That means the original model (1.4) has a stable constant steady state solution, a stable spatially homogeneous periodic solution and two unstable spatially inhomogeneous periodic solutions. Then, we choose (d 2 , g 2 ) = (−0.42, −0.02) ∈ R 6 and the initial value u(x, 0) = 0.98 − 0.001 cos(x), v(x, 0) = 0.98 − 0.001 cos(x), as shown in Figure 9.

Conclusions and discussions
In this paper, the Gierer-Meinhardt model with different sources with Neumann boundary conditions has been considered. Firstly, we focus on the ODEs system of the Gierer-Meinhardt model (1.4), we investigate the condition of the existence for Hopf bifurcation, and find that Hopf bifurcation is degenerate. Therefore, we further study the order of degenerate Hopf bifurcation. With the aid of the symbolic language Maple, we get the normal form associated degenerate Hopf bifurcation, we find that the system (2.1) undergoes a Bautin bifurcation. Furthermore, we get the the universal unfolding of system (2.1) for Bautin bifurcation so that we can identify the stability of periodic solutions. Next we mainly consider the diffusive Gierer-Meinhardt model (1.4), the Turing instability and the existence of the Turing-Hopf bifurcations are investigated. Then, by employing the multiple time scale technique, we obtain the amplitude equations for the Turing-Hopf bifurcation. Furthermore, by analyzing the obtained amplitude equations, the spatiotemporal patterns and their stability are classified. From numerical simulations, we find that the original Gierer-Meinhardt model may exhibit the nonconstant steady state, spatially homogeneous periodic solution and the spatially inhomogeneous periodic solution. The numerical results well confirmed the theoretical analysis we have obtained above. In our work, we choose D 2 and g as the bifurcate parameters. That is to say, the formation of spatiotemporal patterns is related to the concentrations of inhibitor, and the concentrations of inhibitor has a direct relationship with the exogenous sources. When the exogenous sources of activator are taken into account, the same spatiotemporal patterns have been found in [7]. We find that the absence of the exogenous sources of activator can not change the shape of spatiotemporal patterns, but the exogenous sources would affect the speed of patterns formation. Then, how about the spatial resonance in the Gierer-Meinhardt model? The spatial resonance in the Gierer-Meinhardt model need to be discussed in the future work.