Abstract: In this paper, we consider the dynamics of delayed Gierer–Meinhardt system, which is used as a classic example to explain the mechanism of pattern formation. The conditions for the occurrence of Turing, Hopf and Turing–Hopf bifurcation are established by analyzing the characteristic equation. For Turing–Hopf bifurcation, we derive the truncated third-order normal form based on the work of Jiang et al. [11], which is topologically equivalent to the original equation, and theoretically reveal system exhibits abundant spatial, temporal and spatiotemporal patterns, such as semistable spatially inhomogeneous periodic solutions, as well as tristable patterns of a pair of spatially inhomogeneous steady states and a spatially homogeneous periodic solution coexisting. Especially, we theoretically explain the phenomenon that time delay inhibits the formation of heterogeneous steady patterns, found by S. Lee, E. Gaffney and N. Monk [The influence of gene expression time delays on Gierer–Meinhardt pattern formation systems, Bull. Math. Biol., 72(8):2139–2160, 2010.]
Keywords: Gierer–Meinhardt system, delay, Turing–Hopf bifurcation, normal form, spatiotempo- ral patterns.
Article
Turing–Hopf bifurcation and spatiotemporal patterns in a Gierer–Meinhardt system with gene expression delay
Zhao Shuangrui wanghb@hit.edu.cn
Wang Hongbin wanghb@hit.edu.cn
Jiang Weihua jiangwh@hit.edu.cn
Recepción: 10 Febrero 2020
Revisado: 22 Septiembre 2020
Publicación: 01 Mayo 2021
In developmental biology, embryonic development is mediated by morphogens. It is a sig- nal molecule that determines the location, differentiation and fate of many surrounding cells [9]. In [26], Turing showed that two diffusible morphogens could instigate diffusion- driven symmetry breaking and bifurcation. Diffusion can destroy the stability of spatial homogeneous steady state, that is, the stability process can evolve into an instability with diffusion effect. Since then, a large amount of literature research on Turing instability mechanism has emerged in developmental biology. A fundamental property for insta- bilities in such systems is short-range activation and long-range inhibition [13]. In de- velopmental biology, the Gierer–Meinhardt model is used as a classic example to explain the mechanism of pattern formation. The so-called Gierer–Meinhardt model is [6]
(1)where u(x, t) and v(x, t) are concentrations of activator and inhibitor at (x, t), respec- tively, and 
. are all positive constants.
Numerous studies have been done for model (1). The detailed stability and diffusion driven instability, i.e., Turing instability of such system are provided in [21, 28]. These re- sults further theoretically verify Turing’s idea. The sufficient conditions for the occurrence of Hopf bifurcation is performed in [18] showing the existence of spatially homogeneous periodic solutions. Recently, Yang et al. [29] have investigated the conditions for the existence of Turing–Hopf bifurcation to reveal the system exhibits various spatiotemporal patterns.
Sakuma et al. [22] have emphasized the role of gene expression in morphogenesis utilizing in situ hybridization, which records mRNA levels as illustrated in developmental self-organisation via Nodal and Lefty gene products in zebrafish mesodermal induction. The timescales of transcriptional and translational are estimated to be in the range of 10 minutes to several hours [16]. However, previous studies usually ignore the role of gene expression delays caused by transcription and translation in kinetics. Therefore, a diffusive Gierer–Meinhardt system with time delay is proposed in [15] as follows
(2)The two activators molecules can reversibly bind to a receptor and eventually induce the production of additional activator (or inhibitor) molecules, but this production has been delayed for period of time 
. For simplicity, we assume throughout time delay of both gene expression events is same and constant [15].
A large number of literature studies show that dynamics depends crucially on the time delay parameter. Time delay will destroy the stability of the steady state and lead to a temporally periodic solution, that is, the Hopf bifurcation occurs [3, 8, 19, 20]. In [2], Chen et al. have studied the Hopf bifurcation analysis for model (2). It is proved that system possesses temporal pattern. Lee et al. [15] have explored the influence of gene expression time delays on pattern formation, what is more, proved that the delayed Gierer–Meinhardt model exhibits abundant Turing pattern through numerical simulation.
Furthermore, in this paper, we will focus on the joint effect of diffusion and time delay on the pattern for model (2). It is noted that Jiang et al. [12] have mainly considered the impact of diffusion and delay on the Schnakenberg model, which arises from simple chemical reaction systems with limit cycle behavior [23], and found that the system exhibits affluent dynamic behavior (the unfoldings for normal forms at Turing–Hopf sin- gularity are case Ia and case III according to [7]). We observe that the Gierer–Meinhardt model has more complex nonlinear terms and based on the work of Lee et al. [15] in which revealed model exhibits numerous interesting Turing pattern by numerical simulation. Therefore, these findings have inspired us to study the Gierer–Meinhardt model to explore how diffusion and delay essentially affect the formation of pattern and whether the model can produce more complex Turning patterns comparing with the Schnakenberg model.
In this paper, we will build the existence of Turing, Hopf and Turing–Hopf bifurcation for the delayed Gierer–Meinhardt model. Firstly, we get the Turing bifurcation curve by analyzing the characteristic equation, which is continuous and piecewise smooth, and system undergoes Turing–Turing bifurcation at the nonsmooth points. Especially, we obtain the spatial inhomogeneous steady state in multifarious profiles, which depends on the wave number. This provides a theoretical explanation for the existence of spatial inhomogeneous periodic solutions with high frequency at low diffusion rates. Next, we take time delay as the bifurcation parameter proving that the model will undergo Hopf bifurcation at critical values. It is worth mentioning that characteristic equation contains a second-order transcendental term, which results in solving the parameter values for the purely imaginary eigenvalues i.is reduced to an eighth-order polynomial of ., hence the critical parameter values can hardly be explicitly solved. Finally, we will mainly focus on the interaction of diffusion and delay on the model from the perspective of Turing–Hopf by method of combining the central manifold theorem and the normal form theory [4,5,10,27]. It is worth noting that codimension-2 Turing–Hopf bifurcation is usu- ally to be applied to explain spatiotemporal phenomena in chemical reaction, predator– prey models, developmental biology, etc. [1, 12, 24, 25, 29, 30]. Jiang et al. [11] have derived the formulas of calculating normal forms for a general delayed reaction diffusion equation with Neumann boundary condition, which can greatly simplify the complexity of calculation. By employing these formulas we theoretically prove the existence of various spatiotemporal patterns instead of computational simulations [14, 15], such as semistable spatially inhomogeneous periodic solutions, as well as tristability of a pair of spatially inhomogeneous steady states and a spatially homogeneous periodic solution coexisting, in addition, quantitatively give the specific existence region of various forms of solutions near the Turing–Hopf singularity.
This paper is organized as follows. In Section 2, by analyzing characteristic equations
at the positive constant steady state we have established the conditions for Turing, Hopf and Turing–Hopf bifurcation. In Section 3, normal forms truncated to order 3 of the delayed Gierer–Meinhardt systems in the neighborhood of Turing–Hopf singularity are derived by applying normal form method [4, 27] and generic formulas evolved in [11]. In Section 4, we analyze the reduced Gierer–Meinhardt systems with gene expression delay and present that the systems exhibits various interesting spatial, temporal and spa- tiotemporal patterns. Moreover, numerical simulations are shown to illustrate the previous theoretical results. Finally, there is a brief conclusion in Section 5. Throughout the paper, N is the set of all positive integers, and 
represents the set of all nonnegative integers.
In this section, we consider the Turing bifurcation and Hopf bifurcation for system (2) with the homogeneous Neumann boundary condition.
For the sake of convenience, by applying the following scalings [15]:
where Tsis a arbitrary timescale. Dropping the tilde, system (2) becomes the following nondimensionalized system:
(3)where 
Obviously, there is a unique positive equilibrium
. Linearizing system (3) at (u*, v*), we obtain
(4)Let µk, k N., be the eigenvalue of Laplace operator ∆ with Neumann boundary condition in one dimensional spatial domain (0, π). Then µk= K2 and the characteristic equation of (4) is
(5)Where
(6)In particular, for
. = 0, (5) turns into
(7)Where
Throughout this paper, we assume that
(N0) 1 > 2q/(p + 1) − q > 0.
By (N0) we know that all eigenvalues of (7) with .k = 0 have negative real parts.
For any k ∈ N, define
Obviously, DET. k = 0 whenever 
. The following lemma gives the proper- ties of 
.
Suppose that (N0) holds. Then we have
(i) For any fixed
reaches its maximum εmaxat extreme point 
,
and ε*(k, D). is monotonically decreasing .increasing. in D for D > Dm(k) (Dk < D < Dm(k)).
(ii) For any k ∈ N, the equation
has a unique root 
for D, which is given by
Moreover,
In order to understand the properties of e*(k, D) more intuitively, we give the graph of e*(k, D) for different k in Fig. 1.
In the following, we define
where 
From the following analysis we will note that e*(D) is actually the Turing bifurcation curve and 
is the Turing–Turing bifurcation points, which are plotted in Fig. 1. The properties of e*(D) are similar to the Turing bifurcation curve in delayed Schnakenberg systems, in [12], it gives a detailed explanation so we can omit here.
Assume that (N0) holds and D > γ. Then
(i) If 
for some 
then . is a simple root of (5)with k = k1, and all the other roots of (5)have strictly negative parts for τ = 0. Furthermore, Let λ = .(k, τ, ε) be the root of (5)with k = k1. such that λ(k1, τ, ε.(D)) = 0. Then d.(k1, τ, ε)/de|e=e∗(D) < 0.
(ii) If 
, then . is a simple root of (5)for both k.
and k. + 1.
Proof of Lemma 2. (i) DET. = 0 if and only 
So, for any k ∈N, y = 0 is always a root of (5) with such . when 
Then by the properties of .e*(D) and e*(k, D), when D E(Dk1, k1+1, Dk1−1, k1 ) for some k1 and 
is the root of (5) with k = k1. Moreover, 
is simple, by D>γ we can obtain that
By (N0) we know that TRk < 0 for all k ∈ N. and DETk > 0 for all k ∈ N, k/= k1. So all the other roots of (5) for T = 0 has negative real parts when 
Differentiating (5) with respect to . and due to D > γ, we can deduce that
(ii) The proof is similar to proof of (i), which is omitted here.
From the above analysis we can get the following important conclusions about Turing bifurcation of system (3).
Assume that (N0) holds and D > γ. Then
(i) For 
then (u*, v*) of system (3)is asymptotically stable for τ = 0, and if 0 < ε < ε*(D), then (u*, v*) is unstable.
(ii) For D ∈ (Dk, k+1, Dk−1, k), system (3)will undergoes k-mode Turing bifurca- tion at 
(iii) When D = Dk, k+1, (k, k + 1)-mode Turing–Turing bifurcation occurs at ε = e∗(Dk, k+1).
The definition of .-mode Turing bifurcation, which occurs in the last theorem, as well as ..-mode Hopf bifurcation and (k1, k2)-mode Turing–Hopf bifurcation, which will occur in the following parts, are given on page 6 of [11]. A (k, k + 1)-mode Turing– Turing bifurcation can be defined in the same manner as (k1, k2)-mode Turing–Hopf. The reader can refer to [12], so, it will not be repeated here.
In this section, we study the Hopf bifurcation in the case of 
We will employ the method proposed in [2] and [17] to analyze the distribution of characteristic roots of (5).
For some 
satisfy D.(iwk, τ, ε) = 0, then we have
If w.τ/2 /= π/2 + jπ, j ∈ Z, then let θk = tan(wkτ/2), and we have e−iwkτ = (1 − iθk).(1 + iθk). Separating the real and imaginary parts, we have
(8)Denote
And
We define
If D(wk) /= 0, then we can solve from (8) that
(9)and from (9) we find that wk satisfied
(10)If D(wk) = 0, in order to make sure the solvability of (8) for ., then we have
and hence wk satisfies (10) in this case as well. Simplifying (10), we conclude that wksatisfies a polynomial equation with degree 8:
(11)Where
and
is a positive root of
(12)If w.τ/2 = π/2 + jπ, j ∈ Z, then 
and hence
is still a positive root of (12). From the above analysis we have the following lemma.
For 
satisfies 
If we can give the conditions under which the converse of Lemma 3 is true, then we have obtained the purely imaginary roots of the characteristic equation (5). In [2], Chen et al. provided an effective method to solve this problem. Denote
Therefore, we obtain the following lemma.
(See [2].) For k ∈ No, there is wk > 0 satisfying 
then ±iωk are a simple pair of purely imaginary roots of (5)when
where 
Moreover, 
, there exists 
which is the unique root of (5)for 
for some small 
in addition,
For 
(in the sense that arctan 
the same result.
In order to discuss the existence of the positive root of the (12), we have adopted the method in [17].
For h(zk) of (12), we have 
Set
(13)Let 
Then (13) becomes
Where
Define
(14)
(15)h(zk) = 0 has at least one positive root if and only if one of the following conditions is satisfied, where rk, qk, hk are defined in (6), and 
(1 ≤ . ≤ 3) are defined as in (14)and (15), respectively:
(H1) 
.
(H2) 
(H3) 
there exists at least one 
such that 
By Routh–Hurwitz stability criterion it is easy to prove that there exists a 
such that (12) has no positive roots for k > Ko. In other words, (12) exists positive roots only possible for a finite number of 0 ≤ k ≤ Ko.
Denote
and one of conditions (H1), (H2) and (H3) holds .
Suppose K is not empty, and without loss of generality, we assume that it has four positive roots denoted by 
Then (11) has four positive roots say 
Let
Then ±iwk,lis a pair of purely imaginary roots of (5) with 
,
be the root of (5) near 
satisfying
Clearly, the sequence 
is increasing in j, and
Thus, we can define
(16)From the above analysis we arrive at the following conclusion on the Hopf bifurcation of system (3).
Assume that (N0) and ε > ε*(D) are satisfied. τ*, θ*, ω*, k2and l0are defined as in (16).
(i) If none of (H1)–(H3) in Lemma . is satisfied, then all the roots of (5)have neg- ative real parts for all τ ≥ 0. Therefore, equilibrium point 
of system (2)is asymptotically stable;
(ii) If one of (H1)–(H3) in Lemma . is satisfied, then (12)has at least one positive root, all the roots of (5)have negative real parts when τ ∈ [0, τ*). Moreover, 
and 
then when 
, all the roots of (5)have negative real parts except a pair of simple purely imaginary roots ±iw*, system (2)undergoes k2-mode Hopf bifurcation near 
and for 
with some small ϵ > 0, (5)has exactly one pair of conjugate complex roots with positive real parts.
Based on the analysis above, we have obtained the following Turing–Hopf bifurcation theorem.
For system (3), assume (N0) holds, D > γ and one of conditions (H1)–(H3) in Lemma . holds. Given k1 ∈ N0, k2 ∈ K . Then the constant steady state 
is locally asymptotically stable when ε > ε*(D) and 0 ≤ τ < τ*. Moveover, for D ∈ (Dk1,k1+1, Dk1−1, k1 ), system (3)undergoes (k1, k2)-mode Turing–Hopf bifurcation near 
at 
.
In the following, we consider the spatiotemporal pattern of system (3) induced by Turing–Hopf bifurcation. We apply the method in [11] to calculate the normal forms of Turing–Hopf bifurcation for (τ, ε) near the bifurcation point (.T*, ε*). First, normalize the delay T in system (3) by time scaling 
and translate 
into origin. Then system (3) can be translated into
(17)First of all, we define the following real-value Hilbert space
where H2(0, lπ) is a standard Sobolev space, and the corresponding complexification is 
with the complex valued .. inner product
Let 
denote the phase space with the sup norm. We write 
for 
According to [11], system (17) can be rewritten as
(18)where 
diag 
is a bounded linear operator 
function and satisfies 
.
From (18) we can obtain that
(19)with 
and 
are the second and third Fréchet derivative of
of . (α, ϕ) at . = 0, respectively.
By (2.4), (2.6), (2.7), (2.9) and (2.10) of [11] we can get that
(20)with
For 
in a small neighbourhood of (0, 0), it follows from [11] that the normal forms of (18) for Ω = (0, π) up to the third order are
(21)According to [11], the coefficients in (21) can be calculated by the following lemma in the case of 
(which is referred to Turing–Hopf bifurcation of Hopf– Pitchfork type).
In the following, we will give examples to illustrate the various spatiotemporal patterns of system (3) with a fixed set of parameter values.
Referring to [15], choose 
diffusion coefficient D = 0.5. We can obtain
Then 1-mode Turing bifurcation of (3) will occur at e*(0.5) = 0.0824. Through numerical simulation, we can know that .w* = 0.0721 and T* = 2.6753.
Furthermore, under the above given parameters, normal form for (1, 0)-mode Turing– Hopf bifurcation truncated to order 3 is
In addition, the equivalent planar system is
(22)According to [7], the unfolding for (22) is case Ib. The bifurcation set for system (3) in ε, τ -plane is shown in Fig. 2(a), where the critical bifurcation lines are
Therefore, the τ, ε-plane is divided into six regions by critical bifurcation curves
which are denoted by 
in different region
the dynamical behaviors of (22) can be described by corresponding phase portraits; see Fig. 2(b).
For system (3)with p = 0.1, q = 0.8, γ = 0.1 and D = 0.5, when param- eters .τ, ε. are chosen near Turing–Hopf bifurcation point
= (0.0824, 2.6753), system (3)exhibits the following complex dynamics:
(i) When (τ, ε) ∈ D., the steady state
of system (3)is asym.totically stable. Otherwise, the steady state
is unstable, while (τ, ε) ∈. D1.
(ii) When (τ, ε) ∈ D2, system (3)possesses a pair of stable spatially inhomogeneous steady states through Turing bifurcation at
when (τ, ε) crosses L1, which reflects that system exhibits spatial patterns and bistability.
(iii) When 
there are a pair of stable spatially inhomogeneous steady states and an unstable spatially homogeneous periodic solution for system (3). Moreover, there exist a pair of heteroclinic orbits, which connect the spatially homogeneous periodic solution to spatially inhomogeneous steady states, respec- tively. That is, the system exhibits semistable patterns and bistability.
(iv) When 
there are two unstable spatially inhomogeneous periodic solution, which are bifurcated from the spatially homogeneous periodic orbit on ., and a pair of spatially inhomogeneous steady states, and a spatially homogeneous periodic solution for system (3)still remains. Moreover, there are four heteroclinic orbits connecting spatially inhomogeneous periodic solution to the spatially homogeneous periodic solution and spatially inhomogeneous steady states, respectively. Therefore, system shows semistable spatiotemporal patterns and tristability; see Fig. 3.
(v) When
system (3)possesses a stable spatially homogeneous periodic solution and a pair of unstable spatially inhomogeneous steady states. Moreover, there are a pair of heteroclinic orbits connecting spatially inhomogeneous steady states to the spatially homogeneous periodic solution, respectively, which reflects that system exhibits semistable patterns and bistability.
(vi) When
system (3)has a stable spatially homogeneous periodic solution, which represents that system exhibits temporal patterns
In this part, we select
Then 
and T*. = 2.6753. According to [7], the unfolding for the cylindrical coordinate equation of (21) is case III. When (τ, ε) = (2.7253, 0.0445), there are two spatially inhomogeneous periodic solutions; see Fig. 4. If (τ, ε) = (2.6253, 0.0705), two spatially inhomogeneous steady states coexist; see Fig. 5. It is easy to observe that the patterns are different from the last case.
In this paper, we have builded the existence conditions of Turing, Hopf and Turing–Hopf bifurcation for Gierer–Meinhardt system with gene expression delay. The first Turing bifurcation curve 
is given, which is continuous and piecewise smooth, and at nonsmooth points system (3), has undergone Turing–Turing bifurcation. Besides, it has found that diffusion can induce Turing instability, resulting in spatially nonhomogeneous steady states.
Using . as bifurcation parameter, we have further investigated the Hopf bifurcation for system (3) near
. Based on the method of [2] and [17], we have overcome the complexity of solving the purely imaginary roots for a second-order transcendental polynomial and given the critical values . for the occurrence of Hopf bifurcation at which spatially homogeneous periodic solution will be bifurcated from
.
In order to explore the joint effect of diffusion and time delay, we further investigate the Turing–Hopf bifurcation. Utilizing the general formula established in [11], we have derived the normal form for system (3) near the Turing–Hopf singularity. It is theoreti- cally proved that system exhibits spatial, temporal and spatiotemporal patterns, such as semistable spatially inhomogeneous periodic solutions, as well as tristable phenomena of a pair of stable spatially inhomogeneous steady states and a spatially homogeneous peri- odic solution coexisting, which is not found in the another representative Turing model, Schnakenberg system [12]. Significantly, the morphogenetic system can delineate the patterns of animal epidermis, as a consequence, various spatiotemporal patterns provide an explanation for the appearance of animal epidermal patterns. Our research provide an analytical means for understanding the behavior of delayed biological self-organizing systems and the parameter space for when bifurcations occur. The analysis have shown that Turing–Hopf bifurcation provides an explanation for the formation of spatiotemporal pattern. The numerical results well confirmed the frontal theoretical analyses.
Lee et al. have found that the inclusion of time delay induces the generation of oscilla- tion patterns [15, Fig. 5 A, C or Fig. 5 D, F or Fig. 6 D, E], which can be explained by Hopf bifurcation taking time delay as the parameter. Furthermore, time delay will greatly post- pone the formation of the pattern [15, Fig. 5 A, B], which is also confirmed in our research. By analyzing the phase portraits near the Turing–Hopf singularity of the system
we find that the introduction of time delay result in there are semistable patterns from spatially homogeneous periodic solution (or spatially nonhomogeneous periodic solution bifurcated from spatially homogeneous periodic solution) to spatially nonhomogeneous steady state, which attract its nearby solutions, therefore, the occurrence of the semistable patterns will postpone the time to reach the heterogeneous steady pattern. In addition, it is known from Theorem 1 that Turing–Turing bifurcation will occur in model (3). Such kind of degenerate bifurcation is usually used to explain the formation of the superposition of spatial patterns.
