Bifurcation and Patterns Analysis for a Spatiotemporal Discrete Gierer-Meinhardt System

: The Gierer-Meinhardt system is one of the prototypical pattern formation models. The bifurcation and pattern dynamics of a spatiotemporal discrete Gierer-Meinhardt system are investigated via the couple map lattice model (CML) method in this paper. The linear stability of the ﬁxed points to such spatiotemporal discrete system is analyzed by stability theory. By using the bifurcation theory, the center manifold theory and the Turing instability theory, the Turing instability conditions in ﬂip bifurcation and Neimark–Sacker bifurcation are considered, respectively. To illustrate the above theoretical results, numerical simulations are carried out, such as bifurcation diagram, maximum Lyapunov exponents, phase orbits, and pattern formations.


Introduction
In Turing's pioneering work, the diffusion terms on pattern formation for reaction diffusion systems are the main factor [1].In recent decades, Turing instability has been investigated in biology, physics, chemistry, embryogenesis, and others [2][3][4][5].
Many biological, chemical and physical phenomena in nature can be described by the reaction diffusion equation, such as patterns and wave speeds.To describe the spatial patterns formation for the tissue structures in embryology and regeneration, Gierer and Meinhardt proposed several reaction diffusion equations-molecular models in 1972 [6].The Gierer-Meinhardt system is expressed as the following form ∂a(t,x,y) ∂t = ρ 0 ρ + cρ a r (t,x,y) h s − µa(t, x, y) + d 1 ∆a(t, x, y), ∂h(t,x,y) ∂t = c ρ a T (t,x,y) where a(t, x, y) and h(t, x, y) are the density of the activator and inhibitor at time t > 0 and spatial location (x, y) respectively.d 1 and d 2 are the constant diffusion parameters to the activator and inhibitor, respectively; ρρ 0 is the source concentration for the activator, here ρ and ρ are constants; ρ is the one for the inhibitor; the activator and the inhibitor are removed by the first order kinetics µa and vh, respectively.Lots of works on the stability and bifurcation problems of stationary solutions and simulation research have been performed to study the dynamical behaviors of these systems [7][8][9][10][11][12][13][14][15], references therein.Ward and Wei [7] studied the stability and oscillatory instability of symmetric k-spike equilibrium solutions to the Gierer-Meinhardt reaction-diffusion system in a one-dimensional spatial domain for various ranges of the reaction-time constant and the diffusivity of the inhibitor field dynamics.Wei and Winter [8] constructed solutions with a single interior condensation point for the two-dimensional Gierer-Meinhardt system with strong coupling.Ruan [9] investigated the instability of equilibrium points and the periodic solutions under diffusive effects, which were stable without diffusion via the perturbation method for such model with common sources.When r = 2, s = 2, T = 1 and u = 0, Wang et al. [10] studied Hopf bifurcation and Turing instability to such a system.Turing instability to the semi-discrete Gierer-Meinhardt model was considered in [11], and Turing bifurcation and chaos for the spatiotemporal discrete Gierer-Meinhardt were studied in [12].Wu et al. [15] studied Turing instability and Hopf bifurcation to such a system, when r = 2, s = 1, T = 2 and u = 0.The influence of gene expression time delay on the patterns to Gierer-Meinhardt system was explored [16], and the influence of fractional Laplacian on the multi-bump solutions to Gierer-Meinhardt system was explored [17].Moreover, some properties of the development of retinotectal projections in amphibians and fish can be described by using the Gierer-Meinhardt system [18].Most applications of Gierer-Meinhardt system can be seen in Refs.[19][20][21].
Noticed that the coupled map lattices (CMLs) method can illustrate the complex dynamics of competitive system [2], physics system [22], population model [18], etc. Due to the efficiency of its numerical experiments, it is becoming an important branch of nonlinear dynamics.Some existing methods and theoretical results can be applied to study the dynamics of the CMLs model [23][24][25][26][27][28].In addition, the Neimark-Sacker bifurcation and Turing bifurcation analysis were investigated to spatiotemporal discrete nutrient-phytoplankton model with time delay in Ref. [29] by using such methods.We will investigate the model theoretically to determine the conditions for such bifurcations to a modified Gierer-Meinhardt system based on CMLs.Then, the influence of parameters on the patterns formation can be illustrated quantitatively.Hence, the CMLs model can give a better description and prediction of pattern formation, which has been applied to the phytoplankton-zooplankton model [27,28] and predator-prey model [30] for pattern formation.
In order to obtain the spatiotemporal discrete Gierer-Meinhardt system, the CMLs method apply to a modified Gierer-Meinhardt model in this paper.By stability and bifurcation analysis, there are many interesting dynamical behaviors which cannot be generalized by the corresponding continuous Gierer-Meinhardt system for the classical bifurcation analysis, such as defect patterns.In this paper, we find other mechanisms (for example, flip-Turing bifurcation, Neimark-Sacker-Turing bifurcation), except Turing instability mechanism.Based on the analysis of these mechanisms, the circle, spiral of spatial patterns are found.
The remainder of this paper is organized as follows.The Gierer-Meinhardt system with CMLs model and its stability analysis are developed in Section 2 .The detailed theoretical analysis of flip, Neimark-Sacker, Turing bifurcation and co-dimensional 2 bifurcation (such as Neimark-Sacker-Turing bifurcation) are carried out for the spatiotemporal discrete Gierer-Meinhardt system in Section 3. Numerical simulations are provided to illustrate the theoretical results in Section 3 and show the dynamical behaviors and spatial patterns in Section 4.

Stability Analysis
In [15], the authors studied the continuous modified Gierer-Meinhardt system which has the following form ∂a(t,x,y) ∂t where ∆ = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 in this paper.(x, y) is the spatial vector in two-dimensional space, which shows the position of a(t, x, y) or h(t, x, y).By discretizing the model ( 2), the CMLs model is developed as follows.One considers n × n lattices in a two-dimensional square domain.Each lattice is a site (i, j), i, j = 1, 2, ..., n includes two numbers which are the density of activator a(i, j, t) and the density of inhibitor h(i, j, t) at time t = T t + T 0 , where T 0 is initial time.Due to interactions and migration between two species, the density of two will vary with time.When discrete step increases from t to t + 1, the CMLs dynamics of the activator and inhibitor consists of two stages; one is the "reaction" stage and the other is "dispersal" stage .The dispersal stage can be obtained by the spatiotemporal discrete of ( 2) where δ and τ are the space interval and time interval for discretizing Equation (2), respectively.The discrete forms of the Laplacian operator ∆ d which can be shown in the following According to [31], the reaction terms can be obtained by discretizing the time terms of Equation ( 3) where Assume that all the parameters are non-negative and a (i,j,t) and h (i,j,t) are non-negative.The periodic boundary conditions to the CMLs models are considered in this paper.
For all i, j, t satisfy Based on the above analysis, the homogeneous dynamics can be determined by Then, the Equation ( 8) can be rewritten into the form of maps equation Hence, the homogeneous dynamics of equations can be directly analyzed by maps Equation (9).
The fixed points of maps Equation ( 9) are the solutions of the following system Clearly, system (10) has a unique fixed point (a . The Jacobian matrix associated to point (a * , h * ) is defined by According to [32], the fixed point is stable, if the two eigenvalues of |p(τ)| < 1 + q(τ), the fixed point (a * , h * ) is stable.By solving the above inequalities, we have the following result.
Proposition 1.The fixed point (a * , h * ) is stable, if one of the following conditions is satisfied. ( (H4) : In addition, if (H1) or (H2) holds, the fixed point (a * , h * ) is a stable node.If one of the conditions (H3)-(H5) holds, it is a stable focus.Here,

Bifurcation Analysis
In this section, taking τ as the critical bifurcation parameter, we will investigate flip, Naimark-Sacker, Turing bifurcation of system (3).

Flip Bifurcation
As known, if the fixed point (a * , h * ) loses its stability and undergoes a flip bifurcation, then the period-2 points are bifurcated from the fixed point.At the critical value of flip bifurcation, the two eigenvalues of J are λ 1 (τ) = −1 and |λ 2 (τ)| = 1.Hence, the bifurcation τ satisfies the following conditions: The center manifold reduction can be applied to determine the stability of the bifurcated periodic-2 points.Taking τ as an independent, the center manifold reduction should be applied variable into the maps (9), and let w = a − a * , z = h − h * and τ = τ − τ f ; hence, the maps ( 9) is turned into the following form where Moreover, O(4) represents high order (≥4) in the variables (w, z, τ), and where λ 2 = 1 − p(τ f ), and the maps (15) can be transformed into where The center manifold where H( ã, τ) = e 1 p2 + e 2 p τ + e 3 τ2 + O(3).Using the invariance of the center manifold, taking h = H( ã, τ) into the maps, and completing the coefficients of e j , j = 1, 2, 3, then , The maps ( 16) reduced to the central manifold W c (0, 0, 0) are the following form.
The occurrence of flip bifurcation for maps (17) also requires the flip bifurcation theorem in [33].

Neimark-Sacker Bifurcation
The Neimark-Sacker bifurcation occurs at the fixed point (a * , h * ) of maps (9); there exists a pair of conjugate complex eigenvalues of the maps (9); in addition, the modules of two eigenvalues are both 1, which means p(τ) − 4q(τ) < 0 and q(τ) = 1.Hence In addition, the Neimark-Sacker bifurcation theorem [32,33] requires the tranversality condition to be non-zero; in fact where The corresponding two eigenvalues are where here, In addition, the Neimark-Sacker bifurcation [33] requires the discriminatory quantity σ satisfy where with which are the second and third order partial derivatives of G 1 (a * , h * ) and G 2 (a * , h * ) at (0, 0).

Turing Bifurcation
Assume that one of (H1)-(H2) holds.For Turing bifurcation, the eigenvalues for the discrete Laplacian operator ∆ d should be studied in this part.Hence, the characteristic equation for ∆ d is defined by with periodic boundary conditions In [31], one knows that the discrete Laplacian operator ∆ d has the following forms eigenvalues Let ã(i, j, t) = a(i, j, t) − a * and h(i, j, t) = h(i, j, t) − h * , and notice that ∆ d ã(i, j, t) = ∆ d a(i, j, t) and ∆ d h(i, j, t) = ∆ d h(i, j, t), and taking this transformation into Equations ( 4)-( 7), one has ãi,j,t+1 hi,j,t+1 = a 10 a 01 b 10 b 01 ãi,j,t hi,j,t where a 10 , a 01 , b 10 , b 01 are defined in above analysis.If the perturbations are small, the higher order terms can be omitted.Assume that X ij kl as the eigenfunction of the eigenvalues λ kl , and multiplying the both sides of Equation ( 28) by X ij kl , one gets Summing Equation ( 29) for all of i and j, one obtains that Hence, the eigenvalues of the Jacobian of equations of system (31) are where In addition, one defines Z m (τ) = 1 is the threshold condition for the occurrence of Turing bifurcation.In addition, if P T (k, l, τ) 2 − 4Q T (k, l, τ) > 0 holds, the critical value τ T satisfies n max k=1,l=1 Theorem 3. Assume that one of (H1)-(H5) holds, and τ is close to τ T .If Z m (τ) > 1, the homogeneous steady state of model ( 3) with periodic conditions undergo Turing instability, and Turing patterns will occur.If Z m (τ) < 1, the homogeneous steady state is still stable; no Turing patterns will occur.

Numerical Simulations
In this section, we will give some examples to illustrate the results in Section 3. We can find the bifurcations and chaos, as well as many different Turing patterns in this part by numerical simulations.  1 and 2, respectively.Besides, the phase orbits for flip bifurcation and Neimark-Sacker bifurcation are shown in Figures 3 and 4, respectively.
From Figure 3, we can see that there exists stable periodic-2, 4, 8, 10 points as τ great than τ f sightly, as shown in Figure 3a-c,e, respectively.We can also find chaos (Figure 3d,f) as τ increases furthermore.
From Figure 4, we can see that there exists a stable limit cycle as τ great than τ N sightly, as shown in Figure 4a.We can also find chaos (Figure 4c) as τ increases furthermore.
During the route to chaos, we find that there exists periodic-11 windows, as shown in Figure 4b.

The Dynamics Behaviors for Spatially Heterogenous State
In this part, we show the spatiotemporally dynamics of Turing instability for flip bifurcation and Neimark-Sacker bifurcation.In order to ensure that a and h are nonnegative, we need In this sequel, we show the patterns for Neimark-Sacker Turing instability.Let τ = 4.11, the pattern induced by invariant circles, be shown in Figure 7(a1-a3) for t = 1000, t=5000, and t = 10,000, respectively.As τ increases to 4.3223 continuously, the circular pattern induced by periodic-11 orbit is shown in Figure 7(b1-b3) for t = 1000, t = 5000, and t = 10,000, respectively, which is more twisted than the top patterns of Figure 7.As τ reaches to 4.328 eventually, the pattern induced by the homogeneous chaotic oscillating states [27,31] are more twisted than the former two, as shown in Figure 7(c1-c3) for t = 1000, t = 5000, and t = 10,000, respectively,

Conclusions and Future Direction
The flip, Neimark-Sacker and Turing bifurcations of a spatiotemporal discrete Gierer-Meinhardt system are investigated in this paper.In addition, we illustrate the patterns induced by the flip-Turing and Neimark-Sacker Turing instability.Compared to the previous works [3,34], we found that the flip-Turing patterns and Neimark-Sacker-Turing patterns are similar with the patterns induced by the real Ginzburg-Landau equation which emerges as the amplitude equation near a Hopf instability to a continuous reactiondiffusion system, such as the defect turbulence.In fact, the coupled map lattice system is a dynamic system that discretizes time and space but its state variables still remain continuous.Hence, what is the relationship between the patterns of a continuous reaction-diffusion system and a spatiotemporal discrete one?It is worth investigating this phenomenon in further work.

Figure 5 .Figure 6 .
Figure 5. (a) Z m − τ diagram for Turing instability in flip bifurcation; (b) Z m − τ diagram for Turing instability in Neimark-Sacker bifurcation.In this sequel, we show the patterns for flip-Turing instability.Let τ = 3.6; the pattern induced by periodic-2 points[27,31] is shown in Figure6a, which is formed by two alteration states.Similarly, taking τ = 3.8 and τ = 3.96, the patterns induced by periodic-4 and -8 points are shown in Figure6b,c, respectively.Besides, on the paths to chaos, there exists a pattern induced by periodic-10 points, which is shown in Figure6e.Taking τ = 4.02 and τ = 4.18, the patterns induced by chaotic attractors are shown in Figure6d,f, respectively.It is clear to observe that the pattern in Figure6fis more fragmented than that one in Figure6d, which means that the self-organized symmetric patterns are broken and spatial chaotic characteristics are shown.