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

1.   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 activator-inhibitor 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

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

where $U(x, t)$ and $V(x, t)$ denote the concentrations of activator and inhibitor, respectively; $\rho_{0}\rho$ is the exogenous source of activator; $\theta_{1}$ and $\theta_{2}$ are the degradation coefficients of the activator and inhibitor; $c\rho$ and $c_{1}\rho_{1}$ are cross-reaction coefficients; $D_1$ and $D_2$ are the diffusion coefficients of activator and inhibitor, respectively; $c\rho\frac{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}\rho_{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; $\Delta = \frac{\partial^2}{\partial x^2}$ denotes the Laplacian operator; $\Omega = (0, \pi)$ is a bounded domain in $\mathbb{R}$ with the boundary $\partial\Omega$. 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.

2.   Existence of the Turing-Hopf bifurcatin

2.1.   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^*) = (\frac{g}{e}, \frac{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 $\lambda_{1, 2} = \pm i\omega^{*} = \pm i\sqrt{g} = \pm 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 $\bar{u} = u-\frac{g}{e}$ and $\bar{v} = v-\frac{g}{e}$ into system (2.1), then we have

where

when $g = 1$, we have $\lambda_{1, 2} = \pm i\omega^{*} = \pm i\sqrt{g}$, and $(\omega^{*})^2 = g = 1 > 0$. Let $T$ be the matrix that transforms the Jacobian matrix into a standard form

Under the transformation

system (2.3) becomes

where

therefore

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

then

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

and (2.3) turns into

where

and

Since the linear part of the system is not the standard form, therefore from the above transformation, then we have

Let $u = x_{1}$, $v = x_{2}$ 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\times 0 = 2$, indicating that the system has only a pair of purely imaginary eigenvalues$(\pm i\omega^{*} = \pm 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.

Now, after executing the Maple source code $program 1$ which is mentioned in ^[21] and the input file (see Table 1), one has the output file, which includes $\mathscr{D}_{i}r$, $\mathscr{D}_{i}\phi$ and the solutions $x_{i}$. The results:

and

where the fifth and sixth order terms are omitted for simplicity. From Eqs (2.7) and (2.8), under $g = 1$ the normal form can be obtained

From Eq (2.10), we know that the the second Lyapunov coefficient associated with degenerated Hopf bifurcation $\sigma_{2} = \frac{e^{4}}{3}\neq 0$. So the Bautin bifurcation may arise from $E^{*}$ at $g = g* = 1$. Up to now, we need to further calculate the universal unfolding of system (2.1) for Bautin bifurcation so that we can identify the stability of periodic solutions. From the above discussion, we know that $g$ and $e$ are chosen as bifurcation parameters and two small perturbations $\tau = (\tau_{1}, \tau_{2})$ should be introduced to system (2.1), namely $g = 1+\tau_{1}$, $e = 1+\tau_{2}$. Then using the parameter transformation $g = 1+\tau_{1}$, $e = 1+\tau_{2}$ and the variable transformation

we can transform system (2.3) into the following system

Now, the Maple program can be applied to Eqs (2.11) to obtain the normal form. The normal form is obtained from the output in the following

where

From ^[22], we have the universal unfolding of system (2.1) for Bautin bifurcation is

where

Therefore, from the universal unfolding (2.13) we can study the distribution of limit cycles in system (2.1).

Remark 2.1. Since $\sigma = 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.

2.2.   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^* = (\frac{g}{e}, \frac{g}{e})$. So the characteristic equation for the positive equilibrium $E^*$ is

where

with

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_{2}^{n_{*}}, g^{*})$, where $D_{2}^{n}$, $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)\in R_{1}$ and unstable for $(D_{2}, g)\in R_{2}$. Here $R_{1}$ and $R_{2}$ are defined by

with

and

where

(2) The critical line of the Turing instability is defined by $L_{2}: g = g_{n_{*}}(D_{2}), D_{2} > 0$; the Turing instability occurs for $(D_{2}, g)\in\tilde{R_{2}}$, where

(3) System (1.4) undergoes the Turing-Hopf bifurcation at the point $(D_{2}^{n_{*}}, 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\in 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\in \left(0, \sqrt{\frac{1}{D_{1}}}\right)$. If $n\in \left(0, \sqrt{\frac{1}{D_{1}}}\right)$ exists which means $0 < D_{1} < 1$, then there is the Turing instability. Otherwise there is no Turing instability. Therefore when $n\in \left(0, \sqrt{\frac{1}{D_{1}}}\right)$ and $0 < D_{1} < 1$, the Turing instability exists.

The critical lines of Hopf and Turing bifurcations intersect at the point $P^{*} = (D_{2}^{n}, g^*)$ where

and when $n\in \left(0, \sqrt{\frac{1}{D_{1}}}\right)$, $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.

Case (1) $0 < D_{1}\leq\sqrt{2}-1$, then when $n\in\left(0, \left[\sqrt{\frac{\sqrt{2}-1}{D_{1}}}\right]\right]$, $k(n)$ is monotonically increasing; when $n\in\left(\left[\sqrt{\frac{\sqrt{2}-1}{D_{1}}}\right], \sqrt{\frac{1}{D_{1}}}\right)$, $k(n)$ is monotonically decreasing;

Case (2) $\sqrt{2}-1 < D_{1} < 1$, then when $n\in\left(0, \sqrt{\frac{1}{D_{1}}}\right)$, $k(n)$ is monotonically decreasing,

where $[\cdot]$ is the integer part function. Let $n_{1} = 1$ and $n_{2} = \left[\sqrt{\frac{1}{D_{1}}}\right]$.

Next, under Case (1), that is $0 < D_{1}\leq\sqrt{2}-1$. We will discuss the value of the critical wave number $n_{*}$ in the following cases.

① if $k(n_{1})\leq k(n_{2})$, then we choose $n_{*} = n_{1} = 1$;

② if $k(n_{1}) > k(n_{2})$, then we choose $n_{*} = n_{2} = \left[\sqrt{\frac{1}{D_{1}}}\right]$.

Under Case (2), we know that then when $n\in\left(0, \sqrt{\frac{1}{D_{1}}}\right)$, $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_{2}^{n_{*}}, g^{*})$, $\Delta_{0}(\lambda) = 0$ has a pair of purely imaginary roots $\pm i\sqrt{J_{0}}$, and $\Delta_{n_{*}}(\lambda) = 0$ has a root $\lambda = 0$ with a negative root $\lambda = -T_{n_{*}}(g)$. For the other wave number $n\neq0, n_{*}$, all roots of $\Delta_{n}(\lambda) = 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_{2}^{n_{*}}, g^{*})$, then the proof is completed.

3.   Amplitude equations for the Turing-Hopf bifurcation

In this subsection, the spatiotemporal dynamics near the Turing-Hopf bifurcation point $(D_{2}, g) = (D_{2}^{n_{*}}, 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

Let $U = (u, v)^T$, then the system (3.1) can be transformed into

where

Set

and

Step 2. The solution of system (3.2) near the Turing-Hopf bifurcation point $(D_{2}^{n_{*}}, 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 $\varepsilon$:

In the same way, expand the Turing-Hopf bifurcation parameters $D_{2}$, $g$ and the nonlinear term $\mathbf{N}$ in the small neighborhood with the small parameter $\varepsilon$. Then we have

where

Next, we expand the time scale of system (3.2). Let $T_{0} = t$, $T_{2} = \varepsilon^{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 $\varepsilon$, comparing with the coefficients of $\varepsilon$ on the both sides of equation, we get

collecting the coefficients of $\varepsilon^{2}$ on the both sides of equation, we have

comparing with the coefficients of $\varepsilon^{3}$ on the both sides of equation, we obtain

From Eq (3.10), the general solution can be represented as

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 $\mathbf{h_{2}}$, which has the formulas of $u_{1}^{2}$, $u_{1}v_{1}$ and $v_{1}^{2}$. Since these formulas contain

and the conjugate terms of the above terms. Thus the solution of Eq (3.11) can be expressed as

where $\mathbf{b_{1}}$, $\mathbf{b_{2}}$, $\mathbf{b_{3}}$, $\mathbf{b_{4}}$, $\mathbf{b_{5}}$ and $\mathbf{b_{6}}$ are undetermined coefficients. Putting (3.14) into (3.11), comparing with the coefficients of similar terms $W_{1}^{2}e^{2i\omega^{*}T_{0}}$, $|W_{1}|^{2}$, $W_{2}^{2}e^{2in_{*}x}$, $|W_{2}|^{2}$, $W_{1}W_{2}e^{i\omega^{*}T_{0}+in_{*}x}$, $W_{1}\bar{W_{2}}e^{i\omega^{*}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 $\mathbf{NST}$ stands for the non-secular terms and

what is more, we have

where

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

where

and the inner product has the form $(\bar{\mathbf{p}}_{j}^{*})^{T}\mathbf{u}_{j} = 0$, $j = 1, 2$.

Thus, from Eq (3.16), we can get the amplitude equations as follows

where

Step 3. Putting the transformations $W_{1} = \rho_{1}e^{i\theta}$, $W_{2} = \rho_{2}$ into system (3.17), we obtain the equivalent amplitude equations of the Turing-Hopf bifurcation in the real coordinates:

Let $\xi_{1} = \mathbf{Re}\{\varphi_{1}\}$, $\xi_{2} = \mathbf{Re}\{\psi_{1}\}$, $\xi_{3} = \mathbf{Im}\{\varphi_{1}\}$, $A_{1} = \mathbf{Re}\{\varphi_{2}\}$, $A_{2} = \mathbf{Re}\{\varphi_{3}\}$, $A_{3} = \mathbf{Re}\{\psi_{2}\}$, $A_{4} = \mathbf{Re}\{\psi_{3}\}$, $A_{5} = \mathbf{Im}\{\varphi_{2}\}$, $A_{6} = \mathbf{Im}\{\varphi_{3}\}$, where $\mathbf{Re}\{\cdot\}$ and $\mathbf{Im}\{\cdot\}$ represent the real part and the imaginary part of $\cdot$, respectively. Then after we truncate the third order terms and remove the azimuth terms, the amplitude equations of the Turing-Hopf bifurcation become

4.   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_{2}^{n_{*}}, 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_{2}^{n_{*}} = 4$ and

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.

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_{2}^{n_{*}}, g^{*})$.

Next, we will calculate the equilibria of system (4.1). Furthermore, note that $\rho_{1} > 0$, and $\rho_{2}$ is an arbitrary real number. Hence, the different constant solutions of system (4.1) are:

Therefore, the critical bifurcation lines are performed as follows

We know that the boundary equilibria $Q_{1}$ and $Q_{2}^{\pm}$ are bifurcated from the origin on the critical line $T_{1}$ and $T_{2}$, respectively. The two internal equilibrium points $Q_{3}^{\pm}$ are bifurcated from the boundary equilibria $Q_{1}$ and $Q_{2}^{\pm}$ 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 $\mathscr{R}_{1}-\mathscr{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}^{\pm}$ and $Q_{3}^{\pm}$ 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_{2}^{n_{*}}, 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 $\mathscr{R}_{1}$, system (4.1) only has one equilibrium: $Q_{0}$. The equilibrium $Q_{0}$ is asymptotically stable. That means system (1.4) only has a stable constant positive equilibrium $E^*$. Then, we choose $(d_{2}, g_{2}) = (0.1, 0.5)\in \mathscr{R}_{1}$ and the initial value $u(x, 0) = 1.5-0.001\cos(x)$, $v(x, 0) = 1.5-0.001\cos(x)$, as shown in Figure 4.

In region $\mathscr{R}_{2}$, system (4.1) has five equilibria: $Q_{0}$, $Q_{2}^{\pm}$ and $Q_{3}^{\pm}$. The equilibria $Q_{0}$ and $Q_{3}^{\pm}$ are unstable, $Q_{2}^{\pm}$ is stable. That means the original model (1.4) has an unstable constant steady state solution, two unstable spatially inhomogeneous periodic solutions and two stable nonconstant steady state solutions. Then, we choose $(d_{2}, g_{2}) = (0.5, 0.01)\in \mathscr{R}_{2}$ and the initial value $u(x, 0) = 1.01-0.51\cos(x)$, $v(x, 0) = 1.01-0.51\cos(x)$, as shown in Figure 5.

In region $\mathscr{R}_{3}$, system (4.1) has three equilibria: $Q_{0}$, $Q_{2}^{\pm}$. The equilibria $Q_{0}$ and $Q_{2}^{\pm}$ are stable. That means the original model (1.4) has a stable constant steady state solution and two stable nonconstant steady state solutions. Then, we choose $(d_{2}, g_{2}) = (0.1, 0.0015)\in \mathscr{R}_{2}$ and the initial value $u(x, 0) = 1.0015-0.001\cos(x)$, $v(x, 0) = 1.0015-0.001\cos(x)$, as shown in Figure 6.

In region $\mathscr{R}_{4}$, system (4.1) has four equilibria: $Q_{0}$, $Q_{1}$ and $Q_{2}^{\pm}$. The equilibria $Q_{0}$ and $Q_{1}$ are unstable, $Q_{2}^{\pm}$ is stable. That means the original model (1.4) has a unstable constant steady state solution, a unstable spatially homogeneous periodic solution and two stable nonconstant steady state solutions. Then, we choose $(d_{2}, g_{2}) = (-0.01, -0.03)\in \mathscr{R}_{4}$ and the initial value $u(x, 0) = 0.97-0.001\cos(x)$, $v(x, 0) = 0.97-0.1\cos(x)$, as shown in Figure 7.

In region $\mathscr{R}_{5}$, system (4.1) also has four equilibria: $Q_{0}$, $Q_{1}$ and $Q_{2}^{\pm}$. The equilibrium $Q_{2}^{\pm}$ is stable, and the equilibria $Q_{0}$ and $Q_{1}$ are unstable. That means the original model (1.4) has two stable nonconstant steady state solutions, a unstable constant steady state solution and a unstable spatially homogeneous periodic solution. Then, we choose $(d_{2}, g_{2}) = (-0.01, -0.00244)\in \mathscr{R}_{5}$ and the initial value $u(x, 0) = 0.99756-0.0091\cos(x)$, $v(x, 0) = 0.99756-0.1\cos(x)$, as shown in Figure 8.

In region $\mathscr{R}_{6}$, system (4.1) also has four equilibria: $Q_{0}$, $Q_{1}$ and $Q_{3}^{\pm}$. The equilibrium $Q_{0}$ and $Q_{1}$ are stable, and $Q_{3}^{\pm}$ 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)\in \mathscr{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.

5.   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.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Nos. 11971032, 62073114).

Conflict of interest

The authors declare that there are no conflicts of interest.