Bifurcations and Turing patterns in a diffusive Gierer–Meinhardt model

. In this paper, the Hopf bifurcations and Turing bifurcations of the Gierer– Meinhardt activator-inhibitor model are studied. The very interesting and complex spatially periodic solutions and patterns induced by bifurcations are analyzed from both theoretical and numerical aspects respectively. Firstly, the conditions for the existence of Hopf bifurcation and Turing bifurcation are established in turn. Then, the Turing instability region caused by diffusion is obtained. In addition, to uncover the diffusion mechanics of Turing patterns, the dynamic behaviors are studied near the Turing bifurcation by using weakly nonlinear analysis techniques, and the type of spatial pattern was predicted by the amplitude equation. And our results show that the spatial patterns in the Turing instability region change from the spot, spot-stripe to stripe in order. Finally, the results of the analysis are verified by numerical simulations.


Introduction
In general, reaction-diffusion systems [4,14,15] are used to describe models in which the concentration of one or more substances diffuses in space and is affected by the diffusion and inter-conversion of substances.In 1952, A. M. Turing [23] mathematically proposed the conclusion that the homogeneous steady state in a reaction-diffusion system becomes destabilized under certain conditions, that is, the initial steady-state solution of the reaction-diffusion system becomes unstable due to the introduction of a diffusion term.This instability caused by diffusion is often referred to as Turing instability.Thereafter, Turing instability has received a great amount of attention from a wide range of scholars and has become a typical problem in the formation of spatio-temporal patterns [1,7,9,12,16,18,21,26].The various results of pattern formation in the reaction-diffusion system are specified as follows.The Turing-Murray principle was proposed by James Murray [16], which investigated the reaction-diffusion systems of animal bodies and tails and their Turing instability.Schepers and Markus [21] demonstrated that cellular automata can produce Turing patterns in the activator-inhibitor system that is qualitatively consistent with various experiments in chemistry.A diffusion model with a Degn-Harrison reaction scheme is considered by Li et al. [12], and the local and global structure of the steady-state bifurcation is established by the technique of spatial decomposition and implicit function theorem.These works demonstrated that Turing patterns can emerge in a number of ecological and chemical systems.
To uncover the diffusion mechanism of Turing patterns and to examine the actual format of Turing patterns in the real world, we will select the activator-inhibitor model [6] proposed by Gierer and Meinhardt to study the typologies of Turing patterns.The activator-inhibitor model shows that two substances can resist each other's action, and can also be used to depict the formation of polar structures, animal structures, and periodic structures (dots on animals).In recent decades, a large literature has been devoted to the study of this system, as seen in [2,11,13,20,25] and the references therein, which can be written as where (i) u and v represent the concentration of activator and inhibitor respectively, D u and D v are their corresponding diffusion constants, and ∂u ∂t means the change in the concentration of the activator per unit of time.
(ii) ρ u > 0, ρ v > 0 represent the baseline yield of the activator and the inhibitor, separately, and µ u , µ v are the decay rate.
For the Gierer-Meinhardt system (1.1),Ruan [20] demonstrated that diffusion can cause homogeneous equilibrium solutions and homogeneous periodic solutions to become unstable.Liu et al. [13] investigated the multiple bifurcation analysis and spatiotemporal patterns in the one-dimensional Gierer-Meinhardt model.Wu et al. [25] performed a Hopf bifurcation analysis of this diffusion model and studied the direction and stability of Hopf bifurcation by standard central manifold theorem.Stability and Hopf bifurcation analysis on a simplified Gierer-Meinhardt model were studied by Asheghi [2], and the direction of the Hopf bifurcation was obtained by the normal form theory.The investigation conducted by Li et al. [11] pertained to the analysis of Turing patterns observed in a broad-spectrum Gierer-Meinhardt model of morphogenesis.In the particular case, when ρ u = ρ v = 0, a simple scale transformation model is as follows where (i) u and v stand for u(x, y, t) and v(x, y, t), (x, y) ∈ Ω ⊂ R 2 , β denotes the decay rate of the activator.
None of the above-mentioned literature deals with the formation of Turing patterns on two-dimensional space in the Gierer-Meinhardt model.However, for chemical systems, patterns on a two-dimensional plane will be more realistic, more intuitive, and abundant than those on a one-dimensional plane [17,24].For the one-dimensional space, only spot patterns and strip patterns exist.However, in two dimensions, not only spots and strips but also patterns such as spot-strip coexistence and maze shapes may appear.To more clearly understand the mechanisms of pattern formation in Gierer-Meinhardt model, we will study the spatio-temporal evolution pattern of the system (1.2) in two dimensions space.
In this paper, the dynamical behaviors of the system (1.2) are studied by using the decay rate of activator β as a bifurcation parameter.The existing conditions of the Hopf bifurcations and the Turing bifurcations are established in turn.The very interesting and complex patterns (spot patterns, spot-stripe coexistence patterns, and stripe patterns) induced by the Turing bifurcation are analyzed from both theoretical and numerical aspects by a multi-scale method [3,5,27].And our results show that the decay rate of the activator β can affect the dynamical behavior of the system (1.2).The system will occur Turing instability when the decay rate β is within a certain region, the impact of diffusion on the system will be diminished as the decay rate β increases.
The layout of this paper is organized as follows.In Section 2, the conditions for the existence of Hopf bifurcation and the Turing instability with spatial inhomogeneity are discussed analytically.In Section 3, the amplitude equation near the instability threshold is derived using weakly nonlinear analysis, and different solutions to the amplitude equation and its stability are investigated.And the correctness of the theoretical part of the analysis is verified by numerical simulations in space.In Section 4, finally, some conclusions and discussions are given.

Turing instability and bifurcation analysis
In this section, the conditions for the existence of Hopf bifurcation and the Turing instability are discussed.
The local system corresponding to the diffusion system (1.2) is with a unique positive equilibrium The Jacobian matrix computed.The Jacobi matrix taken at the positive equilibrium E * is and the characteristic equation is as follows with Theorem 2.1.For the local system (2.1), when 0 < β < 1, the positive equilibrium E * is locally asymptotically stable, and the system (2.1) undergoes the Hopf bifurcation at β = 1.
Next, we study the diffusion-driven Turing instability of the diffusion system (1.2) under the basic assumption that the constant equilibrium E * (u * , v * ) of the system (1.2) is asymptotically stable (0 < β < 1).
In order to study the linear stability of the constant equilibrium E * (u * , v * ) of (1.2), we need to study the distribution of the roots of the characteristic equation of (1.2).The linearization of Equation (1.2) at the constant equilibrium point Assume the solution of (2.4) is that where k denotes the wave number with the expression k = (k x , k y ), and satisfies k = |k|.r is the spatial vector in two dimensions whose expression is r = (x, y).We can get the corresponding characteristic matrix is The characteristic equation is where Under Theorem 2.1, we have 0 < β < 1, thus for any positive natural number k, there always exist T k < 0. Then the instability condition of the positive equilibrium point E * (u * , v * ) of the system (1.2) should be that: existing a k > 0 make D k < 0. In other words, when D k < 0 (k > 0) is satisfied, there exists a diffusion-driven Turing instability.Since β > 0, the sufficient condition for D k < 0 is that the following two conditions H 1 and H 2 hold and Consider D k as a quadratic function of k 2 , the function D k can obtain the minimum value If H 1 and H 2 hold, then minD k T < 0, which indicates the occurrence of Turing instability.
In the following, we choose β as the parameter to study the conditions that make H 1 and H 2 hold.Regarding the Turing instability of the system (1.2), we obtain the following results.
Theorem 2.2.Assume that the positive equilibrium point E * of the corresponding local system (2.1) is stable, which is given by Theorem 2.1.For the reaction-diffusion system (1.2) T , 1) and Turing bifurcation occurs at where Proof.(I) From Theorem 2.1, we know that the positive equilibrium point E * is stable for 0 < β < 1.Therefore, when σ 1 ≥ σ 2 , we have σ 1 σ 2 ≥ 1 > β, hence H 1 is not satisfied.The conclusion (I) is proved.
(II) Under the conditions of Theorem 2.1, it is easy to get H 1 equivalent to β * < β < 1, where and H 2 is equivalent to the following condition obviously, Q 1 > 0. This means that h(β) = 0 has two positive roots, which are denoted as β (1) T and β (2) and h(β) > 0 if only and if 0 < β < β (1) T and β > β (2) T .In addition, we can get hence, we have the following inequality, T . (2.12) Therefore, H 1 , H 2 are both satisfied for β (2) T , H 1 , H 2 .Then we can conclude that Turing instability occurs only in the region β (2) T > 1, the positive equilibrium point E * is unstable, hence, there is no Turing instability.The conclusion (i) in (II) is proved.
To support the previous theoretical analysis, taking σ 1 = 0.3, σ 2 = 5, we can obtain β (2) T = 0.3497.According to Theorem 2.2, we know that Turing instability occurs for β ∈ (β T , 1).Therefore, to investigate the Turing pattern formation of system (1.2), we need to ensure that the control parameter β ∈ (0.3497, 1).By increasing the value of parameter β in (0.3497, 1), we can obtain the relationship between Re(λ) and k 2 (see T , which implies that there is no Turing instability.Therefore, T is the necessary condition for Turing instability to occur.In the following, we consider the Hopf bifurcation of the system (1. According to [8], n-mode Hopf bifurcation means that the characteristic equation (2.6) has a pair of purely imaginary roots, while the other roots have non-zero real parts and satisfy the corresponding transversal conditions.Theorem 2.3.Suppose one of the following conditions holds: The system (1.2) occurs 0-mode Hopf bifurcation at β = β H 0 = 1, where the characteristic F 0 (λ) = 0 have a pair of purely imaginary roots and other roots of the characteristic F k (λ) = 0 (k > 0) have negative real parts.Where β * and β (2) T are defined in (2.8) and (2.11).Proof.Since dT 0 dβ = 1 2 , then T 0 = 0 has a unique root β = β H 0 = 1, and obviously the transversal conditions satisfied.Moreover, T k < 0 (k ≥ 1) and T usable D k > 0 always satisfied.Thus the system (1.2) occurs 0-mode Hopf bifurcation.
In the next, we find the spatially inhomogeneous Hopf bifurcation for n ∈ N. Define which is the root of Tn There are the following conclusions.Theorem 2.4.Suppose one of the following conditions holds: The system (1.2) undergoes a n-mode Hopf bifurcation around E * (u * , v * ) at β H n for n ∈ N, where the characteristic equation (2.6) has a pair of purely imaginary roots, while all the other roots of F j (λ) = 0 (j ̸ = n l ) have non-zero real parts.Where β * and β ( T are defined in (2.8) and (2.11).Proof.To find the spatially inhomogeneous Hopf bifurcation points for n ∈ N, we have to seek the roots of ( n l ) 2 (σ 1 + σ 2 ) + 1 = β.Since   Remark 2.6.In Figure 2.2, B and C denote the Turing-Hopf bifurcation points corresponding to the (k T , 0)-mode and (k T , 1)-mode, respectively.Point B is located at the coordinates (2.33,1), while point C is located at (2.24, 1.041).To investigate the dynamical behaviors that may occur near these points, we performed numerical simulations.Notably, in the vicinity of point B and C, we observe spatially homogeneous periodic solutions, non-constant steady-state solutions and spatially homogeneous quasi-periodic solutions.These observations are visually depicted in Figure 2.5.These results provide valuable insights into the behavior of the system near the Turing-Hopf bifurcation point.
This section focuses on the stability, Hopf bifurcation, and Turing instability regions of the diffusive Gierer-Meinhardt activator-inhibitor system (1.2) and obtains the conditions for the occurrence of Turing bifurcation, 0-mode Hopf bifurcation, k-mode Hopf bifurcation.As it is known, pattern formation can be induced by Turing instability.To uncover the diffusion mechanics of Turing patterns, this paper requires us to investigate and analyze the dynamic behavior of the Turing bifurcation.To solve this problem, we will employ the amplitude equation as an effective tool.In the next section, we will consider the amplitude equation of the system (1.2).

The amplitude equation and pattern formation 3.1 The amplitude equation of Turing bifurcation
In this subsection, in order to reveal the effect of diffusion on Turing patterns, the amplitude equation of the system (1.2) near the Turing bifurcation β = β (2) T will be deduced by weakly nonlinear analysis [3,5,27].To begin with, we consider the third order polynomial system of the system (1.2), which can be expressed as where and Applying perturbation techniques to the system (3.1), a small parameter ε is introduced near the critical value β (2) T of the Turing bifurcation and satisfies the following form Meanwhile, the linear operator I can be decomposed into where In all cases, the initial values for u and v are given by (u 0 , v 0 ) = (0.9 + 0.01 cos(2x), 0.9 + 0.01 cos(2x)).
In addition, relating the variable U to the parameter ε can be written as Substituting (3.2) and (3.4) into system (3.1),we obtain the following equation where particularly, Accordingly, multiple time scales are introduced and the derivatives with respect to t are converted to O ε 2 : O ε 3 : where with s 22 and s 32 can be obtained by replacing f by g in s 21 and s 31 , and Firstly, we discuss the first order of ε, while (u 1 , v 1 ) T is the linear combinations that belong to the eigenvectors corresponding to zero eigenvalues.The general solution of equation (3.9) can be composed in the following form For convenience, we define It is clear that (ϕ, 1) T is a zero eigenvector of C k , and by simple calculation, we can obtain . Using the Fredholm solvability condition for (3.10), the zero eigenvectors of the adjoint operator I * T of I T is orthogonal to (3.10) right-hand side, and the eigenvector corresponding to the zero eigenvalues of . Using the Fredholm solvability condition to (3.10) By moving the term, we get the following formula Using the orthogonality condition for (3.10), we can obtain the following equations where Suppose that the solution of (3.10) has the following form where c.c represents the complex conjugate of all the preceding terms.Substituting (3.17) into (3.10),we can derive that Using the Fredholm solvability condition to (3.10), After simplification, we can obtain the following equations where T ) 2 , The solution of the reaction-diffusion system (1.2) at the Turing instability critical point has the following form Combining (3.4), (3.11), (3.17) and (3.20), the amplitude Z j can be transformed into the following form Z j = εM j + ε 2 V j + o(ε 3 ).Determined by the expressions of Z j and Eqs.(3.7), (3.11), (3.16) and (3.19) we can obtain the equation for the amplitude corresponding to Z 1 as follows where Analogously, we can derive two other amplitude equations Using the polar coordinate transform Z j = ρ j exp(iφ j ) (j = 1, 2, 3), where ρ = Z j , and φ j is the polar angle.Then substituting (3.21) into (3.22), the system (3.22)becomes where From the first equation of the system (3.23),there are only two conditions to consider: θ = 0 or π.The system (3.23) is stable for θ = 0, d > 0 and θ = π, d < 0. Hence, the system (3.23) can be reduced to the following form (3.24) As the results in the [17,22,24], by studying the existence and stability of the equilibrium points of the amplitude system (3.24),we know that the amplitude system (3.24) has five types of steady-state solutions with the following conclusions: (1) The amplitude system (3.24) has an equilibrium E 1 = (0, 0, 0), which is stable for µ < µ 2 and unstable for µ > µ 2 ; (2) When µw 1 > 0, the amplitude system (3.24) has an equilibrium E 2 = µ w 1 , 0, 0 , which is stable for µ > µ 3 with w 2 > w 1 > 0; (3) When w 1 + 2w 2 > 0, µ 1 < µ < 0 or w 1 + 2w 2 < 0, µ < 0, the system (3.24) has an equilibrium E (0) , which is always unstable; (4) When w 1 + 2w 2 > 0, µ 1 < µ, the system (3.24) has an equilibrium E (5) When w 2 > w 1 > 0, µ > µ 3 or w 1 < w 2 < 0, µ < µ 3 , the system (3.24) has an equilibrium

which is always unstable;
where By Theorem 2.2, Turing instability occurs at β ∈ (β T , 1) for the system (1.2), that is µ = 3 does not exist.According to the results in [17,24], the existence and stability of the equilibria of the amplitude system (3.24)correspond to the type of spatial patterns of the original system (1.2).E 1 and E (π) 3 correspond to the spot patterns, E 2 and E 4 correspond to the the stripe patterns and the mixed patterns, respectively.In addition, it is easy to know from the above discussion that µ 1 < µ 2 < µ 3 < µ 4 .Consequently, one obtains the following results: (1) The system (3.24)only has a equilibrium E (π) 3 for µ 2 < µ < µ 3 , which is stable, therefore, the system (1.2) only appear spot patterns; (2) When β crosses a critical value so that µ 3 < µ < µ 4 , the system (3.24) has two equilibria E 2 and E (π) 3 , correspondingly, the system (1.2) can occurs mixed patterns; (3) When µ 4 < µ, the system (3.24)only has a equilibrium E 2 , and then, stripe patterns will appear in the system (1.2).
Therefore, we are able to establish a connection between the initial reaction-diffusion equation and the amplitude equation presented in Table 3.1.This linkage not only sheds light on the underlying mechanisms of these mathematical models, but also provides a valuable theoretical framework for further research in this field of study [17].
In this subsection, we derive the amplitude equation (3.24) of the system (1.2) using the weakly nonlinear analysis method and obtain the conditions for the appearance of different Turing patterns.In the next subsection, we will verify theoretical analysis by numerical simulation.
Amplitude system (3.24)The original system (1.2) Mixed pattern Table 3.1:The correspondence between the amplitude system and the original system.

Numerical simulations of pattern formation
In this subsection, we will perform numerical simulations to verify the last part of the theoretical analysis.Taking the Parameters σ 1 = 0.5, σ 2 = 3.6, then we have T = 0.8095.
According to Theorem 2.2, Turing pattern will appear when β ∈ (0.8095, 1).Then we choose β = 0.99, and with simple calculations, the following results can be obtained and E (π) 3 = (1.2414,1.2414, 1.2414) represent a specific range of conditions that correspond to the fourth steady-state solution of the amplitude equation (as defined in (4)).Based on our previous analysis, the appearance of spot patterns in the reaction-diffusion system (1.2) is expected under these conditions (see Figure 3.1).Therefore, we can conclude that the formation of spot patterns in the system is likely to occur under the specified parameter values.
And thus obtain µ > µ 4 .The system (1.2) exhibits stripe patterns (see Figure 3.3), as predicted by previous theoretical findings, when the following conditions are met: ρ 1 = 0.2121, µ > µ 3 , w 2 > w 1 > 0, µw 1 = 0.0556 > 0, and E 2 = (0.2121, 0, 0).The corresponding steady-state solution of the amplitude equation is denoted as (2).From the above analysis, Table 3.2 was derived.effect on the system.Then the conditions for the Hopf bifurcation as well as the Turing bifurcation are established theoretically, and the effects of parameter β on the Hopf bifurcation and Turing bifurcation are discussed numerically.It is shown that under certain conditions, a diffusion-driven Turing instability occurs at the positive equilibrium point E * .For a fixed σ 1 , the Turing instability region in the β − σ 2 plane is surrounded by the Hopf bifurcation curve and the Turing bifurcation curve (see Figure 2.2).It can be concluded that there is no Turing instability for the higher decay rate of the activator.
For studying and analyzing the dynamic behavior near the Turing bifurcation, the corresponding amplitude equations are driven for the system (1.2) near the Turing bifurcation point by the weakly nonlinear analysis method, which can be used to predict the stability of the spatial pattern and its type.Based on theoretical analysis, the system will appear with spot patterns, mixed patterns, and stripe patterns, which can be verified by numerical simulations in the subsection 3.2.The results show that, with β as the adjustment parameter, the spatial patterns in the Turing instability region change from the spot patterns, and spotstripe coexistence patterns to stripe patterns in order.These spatial patterns can not only simulate and explain the chemical oscillations between activator concentrations and inhibitor concentrations in a better way but they can also be applied to medical tests [10].
Figure 2.1(a)) and the relationship between D k and k 2 (see Figure 2.1(b)), where Re(λ) is the real part of λ.From Figure 2.1(a) and Figure 2.1(b), it is easy to see that Re(λ) < 0 and D k > 0 always hold for β < β (2)

Figure 2 . 1 :
Figure 2.1: (a): the graph of the dispersion relation with respect to k 2 for different β; (b): the graph of D k (β) with respect to k 2 for different β.

2 )TRemark 2 . 5 .
Tn l = 0 has a unique root β = β H n for n ∈ N, and obviously the corresponding transversal conditions satisfied.Moreover, it is easy to get that Tn l is monotonically decreasing with respect to n, therefore Tj l(β H n ) > 0 for j < n and Tj l (β H n ) < 0 for j > n.By the proof of Theorem 2.3, we know that D k > 0 for one of the conditions in (I) or (II) holds.Thus the system undergoes n-mode Hopf bifurcation atβ H n .In addition, to more intuitively understand Theorem Theorem 2.2-Theorem 2.4, takingσ 1 = 0.4,we plot the stability regions and the existing region of Turing instability in σ 2 − β plane, as shown in Figure 2.2.According to Theorem 2.1-Theorem 2.4, in D 1 , the positive equilibrium E * is unstable and occurs Turing instability, and β = β (represents Turing bifurcation curve.In D 2 , the positive equilibrium E * is unstable but not occurs Turing instability.In D 3 and D 4 the positive equilibrium E * is asymptotically stable.Moreover, we set σ 1 = 0.4, σ 2 = 0.2, then the 0-mode Hopf bifurcation will occurs at β = β H 0 = 1.Taking β = 0.99 < β H 0 , the system (1.2) can occur the spatially homogeneous periodic solutions (as shown in Figure 2.3).We set σ 1 = 0.4, σ 2 = 3, n = 1, l = 8, thus β H 1 = 1.0531.And the 1-mode Hopf bifurcation will occurs at β = β H 1 .Taking β = 1.01 < β H 1 , the system (1.2) can appear the spatially inhomogeneous periodic solution (as shown in Figure 2.4).When β ∈ (β * , β (2) T ), at least one eigenvalue of D k has positive real part, then the Hopf bifurcating periodic solutions are always unstable.Particularly, for 0-mode Hopf bifurcation, bifurcating periodic solutions are unstable in the interval AB in Figure 2.2.

Figure 2 . 2 :
Figure 2.2: When σ 1 = 0.4, the Turing bifurcation curve and Hopf bifurcation curve in σ 2 − β plane.D 1 is the Turing instability region, D 2 denotes unstable regions in which do not occurs Turing unstable, D 3 and D 4 are both stable regions.And B represents (k T , 0)-mode Turing-Hopf bifurcation point, C stands for (k T , 1)-mode Turing-Hopf bifurcation point.