# Analyzing Turing’s Systems via Dynamic Bifurcation Theory

Xige Yang and Dapeng Li

###### Abstract

In this paper, we introduce a novel approach to study reaction diffusion
systems – dynamic transition theory approach developed in [13]. This approach
generalizes Turing’s classical result (linear stability analysis) on pattern
formation and cast some new insights into Turing’s systems. Specifically, we
studied the Turing’s instability and dynamic transition phenomenon for a
Turing’s system, and expressions of the critical parameters
$L(D_{u},D_{v})$,and $D_{v}^{*}$ are derived. These two simple parameters are
sufficient to provide us enough information on the Turing’s instability result
as well as the dynamic transition behavior of the system. As an application,
based on the method we establish in this paper, we found that the Schnakenberg
system has two different transition types : single real eigenvalue transition
and double real eigenvalues transition. These transition types are interpreted
using phase diagrams.

###### Contents

  1. 1 Introduction
  2. 2 Statement of the problem and main results
    1. 2.1 Turing’s instability
  3. 3 Statement of the problem and main results
    1. 3.1 Turing’s instability
    2. 3.2 A glance into main results
  4. 4 Key results in dynamic transition theory
    1. 4.1 Dynamic transition theory – basic setup
    2. 4.2 Exchange of stability
    3. 4.3 The eigenvalue problem
    4. 4.4 Center manifold reduction
  5. 5 Stability and dynamic transition of Turing’s systems
    1. 5.1 Critical parameters of the Turing’s system
    2. 5.2 Geometric insights
    3. 5.3 Center manifold reduction of Turing’s system
  6. 6 Application to the Schnakenberg system
    1. 6.1 Mathematical form of the Schnakenberg system
    2. 6.2 A necessary and sufficient condition for Turing’s instability
    3. 6.3 The critical parameter of the Schnakenberg system
    4. 6.4 Phase Transition of the Schnakenberg system
      1. 6.4.1 Single real eigenvalue transition of the Schnakenberg system
      2. 6.4.2 The double real eigenvalues transition of the Schnakenberg system.
    5. 6.5 An illustrative example
  7. 7 Discussion

##  1 Introduction

Patterns are universal phenomena in physics, chemistry, biology, geography,
economics, and even sociology. In biology, there had been lots studies on
pattern formation since Alan Turing published his celebrated paper The
Chemical Basis of Morphogenesis, which put forth a mathematical model for
spatial pattern formation. Since then, the mathematical model, now known as a
reaction-diffusion system, states a stable equilibrium solution without
diffusion may become unstable because of diffusion.

One of the most popular methods to understand biological pattern formation is
by using reaction diffusion systems, first introduced by Alan Turing in 1952
[26]. In his paper, Turing introduced the concept of ”pre-pattern” as a
precursor of the real pattern we observe, and hypothesized that these spatial
pre-patterns are generated by biochemicals, which he called morphogen. He did
lots of experiments on his own, mixing morphogens in a well-stirred system,
and found that the uniform steady state is stable to small perturbations.
Moreover, he showed that two different morphogens in the same system can
produce unstable patterns, now known as diffusion-driven instability (DDI).
That is, under suitable choices of parameters, the homogenous steady state of
the system will lose its stability.

Turing’s work had a great influence on interdisciplinary subjects such as
mathematical biology, biophysics, non-equilibrium physical chemistry and
complexity science. For example, in 1960s, IIya Prigogine and his
collaborators followed up Turing’s work, formulated and analyzed a model for
the Belousov-Zhabotinsky reaction, which is found in the 1950s an now known as
a classical example of a self-organizing chemical reaction [28]. Fig. 1 is an
chemical experiment of Belousov-Zhabotinsky reaction, and Fig. 2 is a
numerical simulation of BZ reaction using at different times (pictures
generated by MATLAB). In Prigogine’s work, they proposed a mathematical model
called Brusselator model to explain how Turing’s Pattern is generated, and
their model also explained why the spontaneous creation of order (or
spontaneous symmetry breaking) is not forbidden by the Second Law of
Thermodynamics [20, 19]. This work This work earned him the 1977 Nobel Prize
in Chemistry.

Fig. 1 is an chemical experiment of Belousov-Zhabotinsky reaction, and Fig. 2
is a numerical simulation of BZ reaction using at different times (pictures
generated by MATLAB).

Figure 1: An experiment of BZ reaction. Rotational chemical wave is observed
in seconds, indicating chemical oscillation. Details of the experiment can be
found at http://www.pojman.com/NLCD-movies/NLCD-movies.html. Figure 2: A
numerical simulation of BZ reaction using the Brusselator model at different
time slice (time is nondimensionalized). Under a suitable set of parameters,
rotation chemical waves are reproduced. The simulation is achieved combining
method of line (MOL) with central difference spatial scheme.

Latter developments of Belousov-Zhabotinsky reaction are more or less
extensions of Prigogine’s pioneering works. For example, in Fields et. al
developed a model involving 11 reactions and 12 species [3], and can be
further reduced in a equation system of three species [11, 13]. In 1972,
Alfred Gierer and Hans Meinhardt introduced the concept of activator-inhibitor
system, a special kind of reaction diffusion system [4, 15]. After linear
analysis of the Brusselator model, they noticed that only activator-inhibitor
type system may generate Turing’s instability and pattern formation. Such
activator-inhibitor system with the Neumann boundary condition is given and
analysised in [23], where $u$ is the activator satisfying $f^{\prime}_{u}>0$,
and $v$ is the inhibitor satisfying $g^{\prime}_{v}<0$. Besides, we also refer
the readers to [21], which provided us with necessary conditions for
instability of a higher-dimensional Turing’s system.

Other researching directions on Turing’s pattern consider more on agent-based
models as well as stochastic mechanism. For example, [30] partitioned the
whole domain into small subdomains, and built up a master equation system
based on these subdomains to interpret how domain size can affect the behavior
of Turing’s pattern.[29] studied the robustness property of the pattern on
each subdomain by applying spectral method on external noise. Results of both
papers are validated by numerical methods, and summarized in a more general
review article [14].

On the other hand, the dynamic transition theory developed in [13], which is
established by Ma and Wang, is a powerful mathematical tool to study the
nonlinear dissipative system. Based on the dynamical phase transition theory,
Ma and Wang have studied more than twenty kinds of phase transition phenomena,
including the Taylor’s instability, Rayleigh-Bénard convection in the fluid
dynamics, and ENSO phenomenon in atmospheric circulation. Based on the theory,
we study the phase transition phenomena of the Schnakenberg system.

This paper is organized as follows. In Section 3, we review some previous
results and declare our main works. In Section 4, the dynamic transition
theory will be introduced. In Section 5, the general method of analyzing
Turing’s systems will be addressed, including the necessity and sufficiency
condition for Turing’s instability, the method to derive the critical
parameter, the classification of dynamic transitions, and center manifold
reduction result. In the Section 6, based on the method established in Section
5, we derive the critical parameter $D_{v}^{*}$ of the Schnakenberg system.
Due to the principle of exchanges of stabilities of this system, it tells us
in a clear way that under which condition the Schnakenberg system can generate
the Turing’s instability, which types of dynamic transitions the system may
contain, and the physical meanings of the dynamic transitions.

##  2 Statement of the problem and main results

###  2.1 Turing’s instability

The kind of instability driven by diffusion, which can be widely found in the
activator-inhibitor systems, is so-called Turing’s instability. Turing firstly
noticed that such activator-inhibitor system can generate a stationary
pattern, if the inhibitor diffuses faster than the activator. For over a half
century, Turing’s instability and pattern formation has been studied widely.
In this paper, we study the following system

$\displaystyle\begin{cases}\dfrac{\partial u}{\partial t}=D_{u}\Delta
u+f(u,v),\\\\[10.0pt] \dfrac{\partial v}{\partial t}=D_{v}\Delta
v+g(u,v).\end{cases}$ (2.1)

where $u$ is the concentration of a short-range autocatalytic substance (i.e.
activator), and $v$ is the long-range antagonist (i.e. inhibitor) of $u$.

We obtain the necessary and sufficient condition of Turing’s instability for
the system (3.1). The critical numbers $L_{1}$ and $D_{v}^{*}$ reflecting
Turing’s instability, which are determined by diffusion coefficient $D_{u}$
and $D_{v}$, are derived. These critical numbers $L_{1}$ and $D_{v}^{*}$ can
clearly tell us in which condition the system can generate the Turing’s
instability, and which type of transition it has. The study in this paper can
also point out that the phase transition of Turing’s instability has two
types-single real eigenvalue transition and double real eigenvalues
transition.

##  3 Statement of the problem and main results

###  3.1 Turing’s instability

The kind of instability driven by diffusion, which can be widely found in the
activator-inhibitor systems, is so-called Turing’s instability. Turing firstly
noticed that such activator-inhibitor system can generate a stationary
pattern, if the inhibitor diffuses faster than the activator. For over a half
century, Turing’s instability and pattern formation has been studied widely.
In this paper, we study the following system

$\displaystyle\begin{cases}\dfrac{\partial u}{\partial t}=D_{u}\Delta
u+f(u,v),\\\\[10.0pt] \dfrac{\partial v}{\partial t}=D_{v}\Delta
v+g(u,v).\end{cases}$ (3.1)

where $u$ is the concentration of a short-range autocatalytic substance (i.e.
activator), and $v$ is the long-range antagonist (i.e. inhibitor) of $u$.

We obtain the necessary and sufficient condition of Turing’s instability for
the system (3.1). The critical numbers $L_{1}$ and $D_{v}^{*}$ reflecting
Turing’s instability, which are determined by diffusion coefficient $D_{u}$
and $D_{v}$, are derived. These critical numbers $L_{1}$ and $D_{v}^{*}$ can
clearly tell us in which condition the system can generate the Turing’s
instability, and which type of transition it has. The study in this paper can
also point out that the phase transition of Turing’s instability has two
types-single real eigenvalue transition and double real eigenvalues
transition.

###  3.2 A glance into main results

In [21] studying Turning’s instabilities, the authors derived the following
necessary condition for Turing’s instability. In their context, they
considered the following n-dimensional reaction-diffusion system:

$\displaystyle\begin{aligned} &\dfrac{\partial u_{j}}{\partial t}=D_{j}\Delta
u_{j}+f_{j}(u),\\\
&u=(u_{1},u_{2},...u_{n}),x\in\mathbb{R},t\in[0,T),T>0\end{aligned}$

with initial condition

$\displaystyle u_{j}(x,0)=u_{j0}(x),\quad j=1,2,...,n.$ (3.2)

and the associated linearized equation of (3.2) is

$\displaystyle u_{t}=Au,\quad u(x,0)=u_{0}(x)$ (3.3)

where $A=\\{a_{ij}\\}_{1\leq i,j\leq n}$ then we have the following result
with respect to (3.3) as follows:

######  Theorem 3.1 (Satnoianu et. al 2000).

If the kinetic system (3.3) of the problem (3.2) is s-stable then no Turing
bifurcation is possible from the uniform steady state solution $u_{s}$ for any
$n\geq 1$.

In particular, when $n=2$, the necessary conditions for Turing’s instability
are:

$\displaystyle a_{11}+a_{22}<0,\quad a_{11}a_{22}-a_{12}a_{21}>0,$ (3.4)
$\displaystyle
D_{1}a_{11}+D_{2}a_{22}>\sqrt{D_{1}D_{2}(a_{11}a_{22}-a_{12}a_{21})}>0;$ (3.5)

which is the standard result of Turing’s instability [26, 16]. One of our main
contributions is to generalize the above result, and to give a sufficient and
necessary condition for Turing’s instability when $n=2$:

######  Theorem 3.2 (Theorem 5.1).

If

$det(A)>0,~{}~{}Tr(A)<0$ (3.6)

where $A$ is as in (4.2), and

$\displaystyle
L=-\frac{det(A)}{D_{u}D_{v}}+\frac{D^{2}}{4D_{u}^{2}D_{v}^{2}},$ (3.7)
$\displaystyle D=D_{v}a_{11}+D_{u}a_{22},$ (3.8) $\displaystyle
det(A)=a_{11}a_{22}-a_{12}a_{21}.$ (3.9)

then we get two conclusions:
(1) If $D\leq 0$ and $L\leq 0$, system (3.1) is Turing stable.
(2) If $L>0$, system (3.1) generates Turing’s instability if and only if there
exists an eigenvalue $\lambda_{k}$ of $-\Delta$ such that
$-\sqrt{L}<\lambda_{k}-\frac{D}{2D_{u}D_{v}}<\sqrt{L}$.

Besides, the critical parameter $D_{v}^{*}$ defined in (5.9) can be found as
follows

######  Theorem 3.3 (Theorem 5.2).

Let $D_{u}$ be fixed, and (5.5) holds true, then the critical value of the
transition for system (3.1) is $D_{v}^{*}$, i.e, the system (3.1) generates
the Turing’s instability if and only if $D_{v}>D_{v}^{*}$.

Using a new technique called dynamic bifurcation theory as in Section 5.2., we
will also provide a geometric explanation for Turing’s instability.

Specifically, we give an application of the above result to study solution
behaviors of Schnakenberg reaction diffusion system [22]. This simple system
can produce a great amount of different dynamic patterns with modifications of
several simple parameters. Nevertheless, we propose the following necessary
and sufficient condition for Turing’s instability using center manifold
reduction.

######  Theorem 3.4 (Theorem 6.1).

Assume $\frac{b-a}{a+b}-(a+b)^{2}<0$
(1) if $-\frac{\sqrt{4d}}{r(a+b)}<d\frac{b-a}{a+b}-(a+b)^{2}<0$, i.e.

$\displaystyle\max{(d,1)}(b-a)<(a+b)^{3}<(d+\frac{\sqrt{4d}}{r})(b-a)$ (3.10)

hen system (6.4) is Turing stable.
(2) if $d\frac{b-a}{a+b}-(a+b)^{2}>\frac{\sqrt{4d}}{r(a+b)}$ or
$d\frac{b-a}{a+b}-(a+b)^{2}<-\frac{\sqrt{4d}}{r(a+b)}$

$\displaystyle b-a<(a+b)^{3}<(d-\frac{\sqrt{4d}}{r})(b-a)$ (3.11)
$\displaystyle\text{or~{}}\max{(1,d+\frac{\sqrt{4d}}{r})}(b-a)<(a+b)^{3}$
(3.12)

then system (6.4) generates Tuing’s instability if and only if there exists an
eigenvalue $\lambda_{k_{0}}$ of $-\Delta$ such that
$\beta_{-}<\lambda_{k_{0}}<\beta_{+}$.

Further more, we discussed the case of single and double eigenvalue
bifurcation respectively. For example, in the case of double eigenvalue
bifurcation, we are able to get the explicit expression for the system above

######  Theorem 3.5 (Theorem 6.5).

The system (6.4) has a stable steady state $(0,0)$ for $d<d_{0}$, and
bifurcates a new steady state $(x_{0},y_{0})$ for $d>d_{0}$, if and only if
the following conditions hold true

$\displaystyle A_{1}+A_{4}<0,A_{1}A_{4}-A_{2}A_{3}>0,$ (3.13)

where

$\displaystyle A_{1}=\beta_{\lambda_{i_{1}}1}+2a_{20}x_{0}+a_{11}y_{0},$
$\displaystyle A_{2}=2y_{0}a_{02}+a_{11}y_{0},$ $\displaystyle
A_{3}=2b_{02}x_{0}^{2}+b_{11}y_{0},$ $\displaystyle
A_{4}=\beta_{\lambda_{i_{2}}1}+2b_{20}y_{0}+b_{11}x_{0}.$

and there is a stable attractor $u^{d}$ of the system bifurcated from
$d>d_{0}$ as follows

$\displaystyle u^{d}=x_{0}\xi e_{\lambda_{i_{1}}}+y_{0}\eta
e_{\lambda_{i_{2}}}$

besides, we are also able to provide a geometric interpretation of the above
bifurcation results, using phase plane diagram.

##  4 Key results in dynamic transition theory

###  4.1 Dynamic transition theory – basic setup

Transitions are found throughout our everyday lives. Before studying details
of transition problems, we need a good understanding about the nature world in
advance. The laws of nature are usually represented by differential equations,
which can be regarded as dynamical systems – both finite and infinite-
dimensional. In this section, we briefly introduce the key ingredients of
dynamic transition theory developed by Ma and Wang [13].

We start with the reaction-diffusion system proposed in section 3.1, with
Dirichlet/Neumann boundary condition. Without loss of generality, suppose the
steady state $(u_{0},v_{0})$ =(0,0), and take Taylor expansion of (3.1) at
$(0,0)$ as follows:

$\displaystyle\begin{cases}\dfrac{\partial u}{\partial t}=D_{u}\Delta
u+a_{11}u+a_{12}v+g_{1}(u,v),\\\\[10.0pt] \dfrac{\partial v}{\partial
t}=D_{v}\Delta v+a_{21}u+a_{22}v+g_{2}(u,v),\end{cases}$ (4.1)

where

$A:=\left(\begin{array}[]{cc}a_{11}&a_{12}\\\
a_{21}&a_{22}\end{array}\right)=\left(\begin{array}[]{cc}f_{u}(0,0)&f_{v}(0,0)\\\
g_{u}(0,0)&g_{v}(0,0)\end{array}\right).$ (4.2)

Let

$\displaystyle H=L^{2}(\Omega,\mathbb{R}^{2}),$ (4.3) $\displaystyle
H_{1}=\\{(u,v)\in H^{2}(\Omega,\mathbb{R}^{2}):\frac{\partial u}{\partial
n}=0~{}\text{on}~{}\partial\Omega\\}.$ (4.4)

Define operator $L_{\mu}=A_{\mu}+B_{\mu}$ and $G_{\mu}$ as follows

$\displaystyle L_{\mu},G_{\mu}:H_{1}\rightarrow H,$ (4.5) $\displaystyle
A_{\mu}u=(D_{u}\Delta u,D_{v}\Delta v),$ (4.6)

$B_{\mu}u=\left(\begin{array}[]{cc}a_{11}&a_{12}\\\
a_{21}&a_{22}\end{array}\right)u,$ (4.7) $G_{\mu}(u)=(g_{1}(u),g_{2}(u)).$
(4.8)

So $L_{\mu}$ and $G_{\mu}$ are respectively the linear and nonlinear part of
(4.1). Therefore, system (4.1) can be rewritten as

$\displaystyle\left\\{\begin{array}[]{ll}\dfrac{\partial u}{\partial
t}=L_{\mu}u+G_{\mu}(u),\\\ u(0)=\varphi,\end{array}\right.$ (4.11)

where $u=(u,v)$ ,
$\mu=(D_{u},D_{v},r)\in\mathbb{R}_{3}^{+}=\\{(x_{1},x_{2},x_{3})\in\mathbb{R}^{3}:x_{i}\geq
0,i=1,2,3\\}$ and $G_{\mu}(u)=o(\|u\|_{X_{\alpha}})$.

In the next section, we will introduce our main result. We start with solving
eigenvalue problem of the system (4.11).

###  4.2 Exchange of stability

Bifurcation theory in ODE is already well-developed. It is well-know in the
context of ODE bifurcation theory that a bifurcation happens when the max real
part of the eigenvalues of the right hand side Jacobian matrix changes sign
when some parameter $\mu$ passes a critical value $\mu_{0}$.

A natural question arises: how can we extend the ODE bifurcation theory to PDE
case? Without too many doubts, the answer is still related to max real part of
the eigenvalues. However, we need to handle with the potential difficulty that
PDE problems are of infinitely dimensional.

Luckily enough, as we will see in Section 4.3, under certain conditions, a
large group of infinite dimensional operators have nice spectral properties.
Therefore, for this kind of systems, we can introduce the Exchange of
stability property (first coined by Davis [2] and formally explored by Ma and
Wang [10, 13]),.

######  Definition 4.1 (Principle of Exchange of Stability, PES).

Consider the system (4.11) let $\\{\beta_{i}\in\mathbb{C}|i=1,2,....\\}$ be
all eigenvalues of $L_{\mu}$ (counting multiplicities), and suppose that they
satisfy

$\displaystyle
Re~{}\beta_{i}(\mu)\begin{cases}<0\quad\textbf{if~{}}\mu<\mu_{0},\\\
=0\quad\textbf{if~{}}\mu=\mu_{0},\\\
>0\quad\textbf{if~{}}\mu>\mu_{0},\end{cases}$ $\displaystyle 1\leq i\leq m$
(4.12) $\displaystyle Re~{}\beta_{i}(\mu)\neq 0\hskip 80.0pt$ $\displaystyle
j\geq m+1$ (4.13)

for some $\mu_{0}\in\mathbb{R}^{1}$. Then we say the system (4.11) satisfies
Principle of Exchange of Stability or PES for short.

Definition. 4.1 provided us a natural way to divide the eigenvalues of (4.11)
into two different ways. More specifically, let
$\\{e_{1}(\mu),e_{2}(\mu),...,e_{m}(\mu)\\}$ and
$\\{e^{*}_{1}(\mu),e^{*}_{2}(\mu),...,e^{*}_{m}(\mu)\\}$ be the eigenvectors
of $L_{\mu}$ and its conjugate operator $L^{*}_{\mu}$ corresponding to the
eigenvalues $\beta_{i}(\mu),1\leq i\leq m$ respectively. According to Theorem
3.4 in [9] (to my knowledge, this had been the first time that the theorem
ever appeared)

######  Theorem 4.2 (Spectral decomposition of a linear completely continuous
field).

Let $L=A+B:X_{1}\to X$ be a linear completely continuous field, then we have
following results:

  1. 1.

if $\\{\beta_{k}|k\geq 1\\}\in\mathbb{C}$ are eigenvalues of $L$ and
$\\{e_{1},e_{2},...,e_{m}\\}$, $\\{e^{*}_{1},e^{*}_{2},...,e^{*}_{m}\\}$ be
the corresponding eigenvectors of $L_{\mu}$ and its conjugate operator
$L^{*}_{\mu}$, then

$\displaystyle\braket{e_{i},e^{*}_{j}}=\delta_{ij}$ (4.14)

  2. 2.

$X$ can be decomposed into the following direct sum

$\displaystyle\begin{cases}X=\bar{E_{1}}\oplus E_{2},\\\
E_{1}=\text{span~{}}\\{e_{k}|k\geq 1\\}\cap X,\\\ E_{2}=\\{\mathbf{v}\in
X|\braket{\mathbf{v},e_{k}^{*}}_{X}=0\forall k\geq 1\\};\end{cases}$ (4.15)

  3. 3.

$E_{1}$ and $E_{2}$ are invariant spaces of $L^{-1}$ and

$\displaystyle\lim_{n\to\infty}||L^{-1}\mathbf{v}||_{X}^{1/n}=0,\forall\mathbf{v}\in
E_{2};$ (4.16)

  4. 4.

Let $\\{\gamma_{1},\gamma_{2},...,\gamma_{k}\\}\subset\mathbb{C}$ be
eigenvalues of $L^{-1}$ (counting multiplicity) in the order
$|\gamma_{1}|\geq|\gamma_{2}|\geq...\geq|\gamma_{k}|$,
$\\{f^{*}_{1},f^{*}_{2},...,f^{*}_{k}\\}\subset{X\otimes\mathbb{C}}\ $ be
corresponding eigenvalues ($X\otimes\mathbb{C}$ is the complexification of
$X$) of $(L^{-1})^{*}=L^{*-1}$, and let
$E_{k}^{*}=\text{span}\\{f^{*}_{1},f^{*}_{2},...,f^{*}_{k}\\}$
($E_{k}^{*}=\emptyset$ if $k=0$). If

$\displaystyle\rho_{k+1}:=\sup_{u\in
X\otimes\mathbb{C}}\lim_{n\to\infty}|\braket{L^{-n}u,u}^{1/n}>0|$ (4.17)

then there is an eigenvalue $gamma_{k+1}\in\mathbb{C}$ of $L^{-1}$ with
$|\gamma_{k+1}|=\rho_{k+1}$ and $|\gamma_{k+1}|\leq|\gamma_{k}|$.

######  Sketch of the proof.

Conclusion 1 and 2 are analogs to Jordan’s Decomposition of a finite dimension
matrix (Theorem 3.3 of [9]); Conclusion 3 is a direct outcome of orthogonality
between the two spaces $E_{1}$ and $E_{2}$; Finally, combining Conclusion 1-3
together, we can decompose any vector $\bullet\in X\otimes\mathbb{C}$ in a
proper way, yielding Conclusion 4. Details of the proof can be found in [9].

Applying Theorem (4.2) to (4.11), we see that

$\displaystyle X_{1}=E_{1}\oplus E_{2},\quad X=\bar{E_{1}}\oplus E_{2},$
(4.18)

where

$\displaystyle\begin{cases}E_{1}=\text{span~{}}\\{e_{k}|k\geq 1\\}\cap X,\\\
\bar{E_{1}}=\text{~{}closure~{}of~{}}E_{1},\\\ E_{2}=\\{\mathbf{v}\in
X|\braket{\mathbf{v},e_{k}^{*}}_{X}\forall k\geq 1\\}.\end{cases}$ (4.19)

Another important element in dynamic bifurcation theory is how to classify
transition types. This can be done by making use of Principle 4.1 (PES). The
following theorem is a basic principle of transitions from equilibrium states.
It provides sufficient conditions and a basic classification for transitions
of nonlinear dissipative systems (see Theorem 2.1.3 of [13]).

######  Theorem 4.3 (Classification of transition types).

Consider the system (4.11) with , if it satisfies PES then it always undergoes
a dynamic transition from $(\mathbf{u}_{0},\mu)=(0,\mu_{0})$ (w.l.o.g, we can
set $0$ to be its steady state), and there is a neighborhood $U\subset X$ of
$\mathbf{u}=0$ such that the transition is one of the following three types:

  1. 1.

Continuous Transition: there exists an open and dense set
$\tilde{U}_{\mu}\subset U$ such that for every $\phi\in\tilde{U}_{\mu}$ , the
solution $\mathbf{u}_{\mu}$ of (4.11) satisfies

$\displaystyle\lim_{\mu\to\mu_{0}}\limsup_{t\to\infty}||\mathbf{u}_{\mu}(t,\phi)||=0.$

  2. 2.

Jump Transition: for every $\mu_{0}<\mu<\mu_{0}\epsilon$ with some
$\epsilon>0$, there is an open and dense set $\tilde{U}_{\mu}\subset U$ and a
number $\delta>0$ independent of $\mu$ such that for any
$\phi\in\tilde{U}_{\mu}$,

$\displaystyle\limsup_{t\to\infty}||\mathbf{u}_{\mu}(t,\phi)||\geq\delta>0,$

  3. 3.

Mixed Transition: for every $\mu_{0}<\mu<\mu_{0}\epsilon$ with some
$\epsilon>0$, $U$ can be decomposed into two open (not necessarily connected)
sets $U_{1}^{\mu}$ and $U_{2}^{\mu}$:
$\bar{U}=\bar{U}^{\mu}_{1}\cup\bar{U}^{\mu}_{2}$, $\o=U^{\mu}_{1}\cap
U^{\mu}_{2}$ such that

$\displaystyle\lim_{\mu\to\mu_{0}}\limsup_{t\to\infty}||\mathbf{u}_{\mu}(t,\phi)||=0,$
$\displaystyle\forall\phi\in U_{1}^{\mu}$
$\displaystyle\limsup_{t\to\infty}||\mathbf{u}_{\mu}(t,\phi)||\geq\delta>0$
$\displaystyle\forall\phi\in U_{1}^{\mu},$

here $U_{1}^{\mu}$ and $U_{2}^{\mu}$ are called metastable domains.

Figure 3: An illustration of different transition types. (a): continuous; (b):
jump and (c): mixed.

###  4.3 The eigenvalue problem

Since the operator $L_{\mu}$ is defined on a Hilbert space, and its range is
on another Hilbert space, we need to introduce how to define and solve
eigenvalues problem of a infinite dimensional space.

Denote $\mu_{k}$ and $e_{k}$ to be eigenvalues and eigenvectors of $-\Delta$
on $H$ respectively, due to the spectral decomposition theory [9],
$\\{e_{k}\\}$ is a base of $H_{1}$, where $e_{k}$ and $\mu_{k}$ satisfy the
following equations

$\displaystyle\begin{cases}-\Delta e_{k}=\mu_{k}e_{k},\\\\[5.0pt]
e_{k}=0\text{~{}or~{}}\dfrac{\partial e_{k}}{\partial
t}=0~{}on~{}\partial\Omega,\quad k\geq 0\end{cases}$ (4.20)

and $0=\mu_{0}<\mu_{1}<\mu_{2}<\ldots$ . By the eigenvalue $\\{\mu_{k}\\}$ of
$-\Delta$, we can get the eigenvalues $\beta_{k}^{(i)}$ and the eigenvectors
of $L_{\mu}$ satisfying the following equations

$\displaystyle\left\\{\begin{array}[]{ll}D_{u}\Delta(\varepsilon_{1i}^{k}e_{k})+a_{11}\varepsilon_{1i}^{k}e_{k}+a_{12}\varepsilon_{2i}^{k}e_{k}=\beta_{k}^{(i)}\varepsilon_{1i}^{k}e_{k}\\\\[5.0pt]
D_{v}\Delta(\varepsilon_{2}^{k}e_{k})+a_{21}\varepsilon_{1i}^{k}e_{k}+a_{22}\varepsilon_{2i}^{k}e_{k}=\beta_{k}^{(i)}\varepsilon_{2i}^{k}e_{k}\end{array}\right.$

where $k=0,1,2,3,4,\ldots$, and $i=1,2$.

Denote

$M_{\mu_{k}}=\left(\begin{array}[]{cc}a_{11}-\mu_{k}D_{u}&a_{12}\\\
a_{21}&a_{22}-\mu_{k}D_{v}\\\ \end{array}\right),$ (4.22)

then all the eigenvalues $\beta_{k}^{(i)}$ of $L_{\mu}$ are all eigenvalues of
$M_{\mu_{k}}$. It is easy to see that $\beta_{k}^{(i)}(i=1,2)$ are the two
solutions of the following equation

$\displaystyle\mu^{2}-Tr(M_{\mu_{k}})\mu+det(M_{\mu_{k}})=0,$ (4.23)

and

$\displaystyle\beta_{k}^{(2)}=\frac{Tr(M_{\mu_{k}})+\sqrt{Tr(M_{\mu_{k}})^{2}-4det(M_{\mu_{k}})}}{2},$
$\displaystyle\beta_{k}^{(1)}=\frac{Tr(M_{\mu_{k}})-\sqrt{Tr(M_{\mu_{k}})^{2}-4det(M_{\mu_{k}})}}{2},$

where

$\displaystyle Tr(M_{\mu_{k}})=a_{11}+a_{22}-\mu_{k}(D_{u}+D_{v}),$ (4.24)
$\displaystyle
det(M_{\mu_{k}})=D_{u}D_{v}\mu_{k}^{2}-(a_{22}D_{u}+a_{11}D_{v})\mu_{k}+det(A).$
(4.25)

###  4.4 Center manifold reduction

In physical science, it is often crucial to determine the asymptotic behavior
of a system at the critical threshold. For this purpose and for determining
the structure of the local attractor representing the transition states, the
most natural approach is to project the underlying system to the space
generated by the most unstable modes, preserving the dynamic transition
properties. This is achieved with center manifold reduction [25, 27, 6, 13].

Let $X_{1}$ and $X$ be two Banach spaces with $X_{1}\subset X$ a dense and
compact inclusion. Consider the following one-parameter nonlinear evolution
equation again (in the remaining of this chapter, we do not distinguish
$G(\cdot,\mu)$ and $G_{\mu}(\cdot)$):

$\displaystyle\begin{aligned} &\dfrac{\partial\mathbf{u}}{\partial
t}=L_{\mu}\mathbf{u}+G(\bullet,\mu),\\\
&\mathbf{u}(0,t)=\mathbf{u}_{0}(t),\quad\mathbf{x}\in\Omega,\mu\in\mathbb{R}^{1}\end{aligned}$
(4.26)

where $L_{\mu}=A_{\mu}+B_{\mu}:X_{1}\to X$ is a linear completely continuous
fields, (i.e. $A$ is a linear homeomorphism and $B_{\mu}$ a linearly compact
operator) depending continuously on $\mu$. Suppose that the system (4.26)
satisfies PES 4.1, by the Spectral Decomposition Theorem 4.2, $L_{\mu}$ can be
decomposed into $L_{\mu}=L^{(1)}_{\mu}+L^{(2)}_{\mu}$ such that for every
$\mu$ sufficiently close to $\mu_{0}$,

$\displaystyle\begin{aligned}
L^{(1)}_{\mu}=L_{\mu}|_{E_{1}}:E_{1}\to\tilde{E}_{1},\\\
L^{(2)}_{\mu}=L_{\mu}|_{E_{2}}:E_{2}\to\tilde{E}_{2},\end{aligned}$ (4.27)

where $\tilde{E}_{i}$ is the completion of $E_{i}\subset X_{1}$ in $X$, the
eigenvalues of $L^{(1)}_{\mu}$ have nonnegative real parts and $L^{(2)}_{\mu}$
have negative real parts at $\mu=\mu_{0}$. Therefore, (4.26) can be written as

$\displaystyle\begin{aligned} \dfrac{\partial\mathbf{u}_{1}}{\partial
t}=L^{(1)}_{\mu}\mathbf{u}_{1}+G_{1}(\mathbf{u}_{1},\mathbf{u}_{2},\mu),\\\
\dfrac{\partial\mathbf{u}_{2}}{\partial
t}=L^{(2)}_{\mu}\mathbf{u}_{2}+G_{2}(\mathbf{u}_{1},\mathbf{u}_{2},\mu).\end{aligned}$
(4.28)

where $\mathbf{u}=\mathbf{u}_{1}+\mathbf{u}_{2}$, $\mathbf{u}_{1}\in E_{1}$,
$\mathbf{u}_{2}\in E_{2}$,
$G_{i}(\mathbf{u}_{1},\mathbf{u}_{2},\mu)=P_{i}G(\mathbf{u},\mu)$ and
$P_{i}:X\to E_{i}$ are canonical projections.

Clearly, in finite dimensional case, $L_{\mu}$ is just a $n\times n$ matrix,
denote it by $A$. So (4.28) can be expressed as

$\displaystyle\begin{aligned}
\frac{d\mathbf{u}_{1}}{dt}=A^{(1)}_{\mu}\mathbf{u}_{1}+G_{1}(\mathbf{u}_{1},\mathbf{u}_{2},\mu),\\\
\frac{d\mathbf{u}_{2}}{dt}=A^{(2)}_{\mu}\mathbf{u}_{2}+G_{2}(\mathbf{u}_{1},\mathbf{u}_{2},\mu).\end{aligned}$
(4.29)

The followings are the well-known center manifold theorems. Theorem 4.4 can be
found in standard books for bifurcation theory, e.g. [27, 6]. Theorem 4.5 can
be found in [25, 9].

######  Theorem 4.4 (Finite dimensional case).

Suppose that all the eigenvalues of A have non-negative real parts, and all
the eigenvalues of B have negative (or positive) real parts. Then, for the
system with the condition that all eigenvalues of $A^{(1)}$ is non-negative
(resp. non-positive) and eigenvalues of $A^{(2)}$ negative (resp. positive),
then there exists a $C^{r}$ function

$\displaystyle\Phi(\cdot,\mu):\Omega\to\mathbb{R}^{n-m};\Omega\subset\mathbb{R}^{m}\text{~{}a~{}neighborhood~{}of~{}}\mathbf{x}=0,$

such that $\Phi(\mathbf{x},\mu)$ is continuous on $\mu$ and

  1. 1.

$\Phi(0,\mu)=0,\dfrac{\partial\Phi}{\partial x}(0,\mu)=0$;

  2. 2.

the set

$\displaystyle
M_{\mu}:=\\{(\mathbf{x},\mathbf{y})|\mathbf{x}\in\Omega\subset\mathbb{R}^{m},\mathbf{y}=\Phi(\mathbf{x},\mu)\\}$

called the (local) center manifold, is a local invariant manifold of (4.4);

  3. 3.

if $M_{\mu}$ is positive invariant (or negative invariant), namely
$\mathbf{z}(t,\phi)\in M_{\mu}(\mathbf{z}(-t,\phi)\in M_{\mu}),\forall t\geq
0$provided $\phi\in M_{\mu}$, then $M_{\mu}$ is an attracting set of (4.4)(or
a repelling set) , i.e. there is a neighborhood $U\subset\mathbb{R}^{n}$ of
$M_{\mu}$, as $\phi\in U$, we have

$\displaystyle\lim_{t\to\infty(\text{~{}or~{}}\infty)}dist(\mathbf{z}(t,\phi),M_{\mu})=0$

where $\mathbf{z}(t,\phi)=\\{\mathbf{x}(t,\phi),\mathbf{y}(t,\phi)\\}$ is the
solution of (5.2.3) with the initial condition $z(0,\phi)=\phi$.

######  Theorem 4.5 (Infinite dimensional case).

Suppose (4.26)-(4.28), and assume $X_{1}$ and $X$ are Hilbert spaces, then
there exists a neighborhood of $\mu$ given by $|\mu-\mu_{0}|<\delta$ for some
$\delta>0$, a neighborhood $O_{\mu}\subset E^{\mu}_{1}$ of $x=0$, and a
$C^{1}$ function $\Phi(\cdot,X):O_{\mu}\to E^{\mu}_{2}(\alpha)$ depending
continuously on $\mu$, where $\bar{E}^{\mu}_{2}(\alpha)$ is the completion of
$E^{\mu}_{2}(\alpha)$ in the $X_{\alpha}$-norm, with $0<\alpha<1$ such that

  1. 1.

$\Phi(0,\mu)=0,\dfrac{\partial\Phi}{\partial x}(0,\mu)=0$;

  2. 2.

the set

$\displaystyle M_{\mu}:=\\{(\mathbf{x},\mathbf{y})|\mathbf{x}\in
O_{\mu},\mathbf{y}=\Phi(\mathbf{x},\mu)\in E^{\mu}_{2}(\alpha)\\}$

called the (local) center manifold, is a local invariant manifold of (4.26);

  3. 3.

$(\mathbf{x}_{\mu}(t),\mathbf{y}_{\mu}(t))$ is a solution to (4.26), then
there is a $\beta_{\mu}>0$ and $C_{\mu}>0$ with $k_{\mu}$ depending on
$(\mathbf{x}_{\mu}(0),\mathbf{y}_{\mu}(0))$ such that

$\displaystyle||\mathbf{y}_{\mu}(t)-\Phi(\mathbf{x}_{\mu}(t),\mu)||_{X}\leq
C_{\mu}\exp{(-\beta_{\mu}t)}$ (4.30)

######  Remark (1).

It is noticeable that, Theorem 4.5 only works when $X$ and $X_{1}$ are Hilbert
spaces, due to their spectral properties. By the spectral decomposition
theorem 4.2, however, we can extend the result of Theorem 4.5 to any totally
continuous field $L$ in Banach spaces. Details can be found in Ma and Wang’s
book [9].

######  Remark (2).

Both Theorem 4.4 and 4.5 ascertained the existence of center manifold
functions for finite and infinite dimensional dynamical systems. However, they
do not provide a explicit way for constructing the central manifold. A
systematic construction can be found in Section 3 of [9] and Appendix A of
[13], which is skipped in this context.

While Theorems 4.4 and 4.5 show the significance of a center manifold
function, they do not tell how to explicitly find these functions. In general,
due to our knowledge, the only systematic way of calculating center manifold
functions are to assume polynomial structures of $Gs^{\prime}$ in (4.26), or
to use Taylor’s expansion near some steady state and critical parameter. For
detailed calculation, please refer to Chapter 18 of [27] for finite
dimensional cases, and Section 3.2 of [9] or Appendix A of [13]. For example,
if $G$ in (4.26) has the following form

$\displaystyle G(\mathbf{u},\mu)=\sum_{n=k}^{\infty}G_{n}(\mathbf{u},\mu),$
(4.31)

for some $k\geq 2$, where $G_{n}:\prod_{i=1}^{n}X_{1}\to X$ is an $n$-multiple
linear mapping, and
$G_{n}(\mathbf{u},\mu):=G_{n}(\mathbf{u},...,\mathbf{u},\mu)$. Then we have

######  Theorem 4.6 (Theorem 3.8 of [9]).

If $L_{\mu}$ is a sectorial operator and $G$ satisfies (4.31), then the center
manifold function $\Phi(\mathbf{x},\mu)$ in Theorem 4.5 can be expressed as

$\displaystyle\Phi(\mathbf{x},\mu)=(-L^{(2)}_{\mu})^{-1}P_{2}G_{k}(\mathbf{x},\mu)+O(|Re~{}\beta(\mu)|\cdot||\mathbf{x}||^{k})+o(||\mathbf{x}||^{k})$
(4.32)

where $P_{2}:X_{1}\to\tilde{E}_{2}$ is the canonical projection,
$\mathbf{x}\in E_{1}$ and $\beta(\mu)=(\beta_{1}(\mu),...,\beta_{m}(\mu))$ are
the eigenvalues of $L^{(1)}_{\mu}$.

By (4.32) or equation (A.1.14) of [13], it can be calculated that

$\displaystyle\begin{aligned}
&\Phi(\mathbf{x},\mu)=\sum_{j=m+1}^{\infty}\Phi_{j}(\mathbf{x},\mu)e_{j}\\\
&\Phi_{j}(\mathbf{x},\mu)=-\frac{1}{\beta_{j}(\mu)<e_{j},e_{j}^{*}>_{X,X^{*}}}<G_{j}(\mathbf{x},\mu),e^{*}_{j}>_{X,X^{*}}+o(k)\end{aligned}$
(4.33)

where $\mathbf{x}=\sum_{i=1}^{m}x_{i}e_{i}$ and
$\beta(\mu)=(\beta_{1}(\mu),...,\beta_{m}(\mu))$, $<\cdot,\cdot>_{X,X^{*}}$
denotes the dual product between the function space $X$ and its dual $X^{*}$.

##  5 Stability and dynamic transition of Turing’s systems

In this section, we study the stability and dynamic transition behavior of the
Turing’s system (3.1).

###  5.1 Critical parameters of the Turing’s system

The following theorem is the necessary and sufficient condition for Turing’s
instability of system (3.1). Let

$\displaystyle
L=-\frac{det(A)}{D_{u}D_{v}}+\frac{D^{2}}{4D_{u}^{2}D_{v}^{2}},$ (5.1)

and assume that

$det(A)>0,~{}~{}Tr(A)<0$ (5.2)

where $A$ is as in (4.2), and

$\displaystyle D=D_{v}a_{11}+D_{u}a_{22},$ (5.3) $\displaystyle
det(A)=a_{11}a_{22}-a_{12}a_{21}.$ (5.4)

######  Theorem 5.1.

Under the assumption (5.2), we obtain the following two conclusions: (1) If
$D\leq 0$ and $L\leq 0$, system (3.1 ) is Turing stable. (2) If $L>0$, system
(3.1) generates Turing’s instability if and only if there exists an eigenvalue
$\lambda_{k}$ of $-\Delta$ such that
$-\sqrt{L}<\lambda_{k}-\frac{D}{2D_{u}D_{v}}<\sqrt{L}$.

###### Proof.

First, we prove the first assertion. Let

$\displaystyle\begin{aligned}
h(\lambda,L)&=D_{u}D_{v}\lambda^{2}-(D_{u}a_{22}+D_{v}a_{11})\lambda+a_{11}a_{22}-a_{12}a_{21}\\\
&=D_{u}D_{v}[(\lambda-\frac{D}{2D_{u}D_{v}})^{2}-L]\end{aligned}$

If $D\leq 0$ and $L\leq 0$, combining with (4.24 ) and (5.1), it is easy to
check that $det(M_{\lambda_{k}})=h(\lambda_{k},L)>0$ and
$tr(M_{\lambda_{k}})>0$ holds true for all the eigenvalues $\\{\lambda_{k}\\}$
of the operator $-\Delta$, that is, all eigenvalues of $L_{\lambda}$ have
negative real part, which means that system (3.1)is Turing stable.

Secondly, we prove the second assertion.

Sufficiency. The condition $Tr(A)=a_{11}+a_{22}<0$ and
$det(A)=a_{11}a_{22}-a_{12}a_{21}>0$ mean that $(0,0)$ is a stable steady
state of the system (3.1) without diffusion. Let the solutions of
$h(\lambda,L)=0$ be $k_{\pm}=\frac{D}{2D_{u}D_{v}}\pm\sqrt{L}$. We can deduce
from $k_{-}<\lambda_{k}<k_{+}$ that $det(M_{\lambda_{k}})<0$ , that is, there
exists eigenvalue $\beta_{k}^{i_{0}}$ of $L_{\lambda}$ such that
$\beta_{k}^{i_{0}}>0$. Then $(0,0)$ is unstable steady state of system (3.1),
that is, system (3.1)generates Turing’s instability.

Necessity. Based on the definition of Turing’s instability, $(0,0)$ should be
the stable steady state of (3.1) without diffusion, that is, $a_{11}+a_{22}<0$
and $a_{11}a_{22}-a_{12}a_{21}>0$ hold true. $(0,0)$ is not the stable steady
state of (3.1), which means that there exists $\beta_{k}^{i_{0}}>0$ of
$L_{\lambda}$, such that $h(\lambda_{k},L)=det(M_{\lambda_{k}})<0$. In another
word, there is a eigenvalue $\lambda_{k}$ of the operator $-\Delta$ satisfying
$k_{+}<\lambda_{k}<k_{-}$, where $k_{\pm}=\frac{D}{2D_{u}D_{v}}\pm\sqrt{L}$.
The proof is complete. ∎

######  Remark 5.1.

It is also important to see that, when $h(\lambda,L)<0$, then the system (3.1)
is Turing unstable. This provide us a way to determine the transition
parameter of (3.1) as will be shown in Remark 5.3.

######  Remark 5.2.

$D>0$ and $L>0$ means $\frac{D_{v}}{D_{u}}>1$, which is a necessary condition
of the Turing’s instability. It means that the inhibitor diffuses faster than
the activator if Turing’s instability is achieved.

The eigenvalues of the operator $-\Delta$ with Dirichlet/Neumann boundary
condition in the case that $\Omega=(0,s)$ is as follows

$\displaystyle\\{\frac{\pi^{2}m^{2}}{s^{2}}:m\in N_{+}\\},$

then we can get the following corollary.

######  Corollary 1.

Under the assumptions in Theorem 5.1, assume that $D>0$ and $L>0$, and taking
$\Omega=(0,s)$, then system (3.1) generates Turing’s instability if and only
if there exists a integer $m\in N_{+}$ such that
$-\sqrt{L}<\frac{m^{2}\pi^{2}}{s^{2}}-\frac{D}{2D_{u}D_{v}}<\sqrt{L}$.

Assume that there exists some $i\geq 0$ and $D_{v}^{0}$ such that

$\displaystyle\lambda_{i}<k_{c}^{0}<\lambda_{i+1},~{}~{}detM_{\lambda_{i}}>0,~{}~{}detM_{\lambda_{i+1}}>0,$
(5.5)

where $\lambda_{i}$ and $\lambda_{i+1}$ are the eigenvalues of the operator
$-\Delta$, $M_{\lambda_{i}}(k=i)$ is as (5.5), and

$k_{c}^{0}=\frac{a_{11}D_{v}^{0}+a_{22}D_{u}}{2D_{u}D_{v}^{0}}.$ (5.6)

where $D_{v}^{0}$ is defined in (5.20). It is easy to check that, $k_{c}^{0}$
is the minimal point of the polynomial $h(\lambda,L)$.

Further denote that

$\displaystyle
D_{v}^{\lambda_{i}}=\frac{a_{22}D_{u}\lambda_{i}-det(A)}{D_{u}\lambda_{i}^{2}-a_{11}},$
(5.7) $\displaystyle
D_{v}^{\lambda_{i+1}}=\frac{a_{22}D_{u}\lambda_{i+1}-det(A)}{D_{u}\lambda_{i+1}^{2}-a_{11}},$
(5.8) $\displaystyle
D_{v}^{*}=\min\\{D_{v}^{\lambda_{i}},D_{v}^{\lambda_{i+1}}\\}.$ (5.9)

Note that since $\lambda_{i},i\geq 1$ are bounded below, it is easy to see
that $D_{v}^{*}$ is well defined. Theorem 5.1 can then be improved as follows

######  Theorem 5.2.

Suppose the conditions in Theorem 5.1 and (5.5) holds true, let $D_{u}$ be
fixed, then the critical value of the transition for system (3.1) is
$D_{v}^{*}$, i.e, the system (3.1) generates the Turing’s instability if and
only if $D_{v}>D_{v}^{*}$.

###### Proof.

Consider the following polynomial

$\displaystyle
h(\lambda)=D_{u}D_{v}\lambda^{2}-(a_{11}D_{v}+a_{22}D_{u})\lambda+det(A),$
(5.10)

and let $\lambda=\lambda_{i}$ or $\lambda=\lambda_{i+1}$ . By solving the
following equations

$\displaystyle
h(\lambda_{k})=D_{u}D_{v}\lambda_{i}^{2}-(a_{11}D_{v}+a_{22}D_{u})\lambda_{i}+det(A)=0,$
(5.11) $\displaystyle
h(\lambda_{i+1})=D_{u}D_{v}\lambda_{i+1}^{2}-(a_{11}D_{v}+a_{22}D_{u})\lambda_{i+1}+det(A)=0,$
(5.12)

we get

$\displaystyle
D_{v}^{\lambda_{i}}=\frac{a_{22}D_{u}\lambda_{i}-det(A)}{D_{u}\lambda_{i}^{2}-a_{11}},$
(5.13) $\displaystyle
D_{v}^{\lambda_{i+1}}=\frac{a_{22}D_{u}\lambda_{k+1}-det(A)}{D_{u}\lambda_{i+1}^{2}-a_{11}}.$
(5.14)

Besides,

$\displaystyle\min_{\lambda\in
R}\\{h(\lambda)\\}=det(A)-\frac{(a_{22}D_{u}+a_{11}D_{v})^{2}}{4D_{u}D_{v}},$
(5.15)

which is bounded below in variable $D_{v}$ for fixed $D_{u}$, so we can take

$D_{v}^{(i)}=\inf\\{D_{v}^{\lambda_{i}},D_{v}^{\lambda_{i+1}}\\}.$ (5.16)

Based on (4.25) and (5.5), obviously, $D_{v}^{*}$ is the critical point. ∎

######  Remark 5.3.

Here we give the method to find the $D_{v}^{0}$ and $k_{c}^{0}$. Let

$\displaystyle\inf_{\lambda\in
R}\\{h(\lambda)\\}=det(A)-\frac{(a_{22}D_{u}+a_{11}D_{v}^{0})^{2}}{4D_{u}D_{v}^{0}}=0,$
(5.17)

that is,

$\displaystyle
a_{11}^{2}(D_{v}^{0})^{2}-(4D_{u}det(A)-2a_{11}a_{22}D_{u})D_{v}^{0}+a_{22}^{2}D_{u}^{2}=0,$
(5.18)

Then we choose the larger positive root, and get

$\displaystyle k_{c}^{0}=\frac{a_{11}D_{v}^{0}+a_{22}D_{u}}{2D_{u}D_{v}^{0}},$
(5.19) $\displaystyle
D_{v}^{0}=\frac{q+\sqrt{q^{2}-4a_{11}^{2}a_{22}^{2}D_{u}^{2}}}{2a_{11}^{2}},$
(5.20)

where

$q=4D_{u}det(A)-2a_{11}a_{22}D_{u}.$

######  Remark 5.4.

Without loss of generality, let $D_{*}=D_{v}^{(1)}$, then the $PES$ condition
is shown as follows:

If $D_{v}^{\lambda_{i}}\neq D_{v}^{\lambda_{i+1}}$, then

$\displaystyle\beta_{i}^{(2)}\left\\{\begin{array}[]{ll}<0,~{}if~{}D_{v}<D_{*},\\\
=0,~{}if~{}D_{v}=D_{*},\\\ >0,~{}if~{}D_{v}>D_{*}.\end{array}\right.$ (5.24)
$\displaystyle\beta_{k}^{(j)}(D_{*})<0,~{}for~{}all~{}(k,j)\neq(i,2).$ (5.25)

If $D_{v}^{\lambda_{i}}=D_{v}^{\lambda_{i+1}}$, then

$\displaystyle\beta_{i}^{(2)}=\beta_{i+1}^{(2)}\left\\{\begin{array}[]{ll}<0,~{}if~{}D_{v}<D_{*},\\\
=0,~{}if~{}D_{v}=D_{*},\\\ >0,~{}if~{}D_{v}>D_{*}.\end{array}\right.$ (5.29)
$\displaystyle\beta_{k}^{(j)}(D_{*})<0,~{}for~{}all~{}(k,j)\neq(i,2)~{}and~{}(i+1,2).$
(5.30)

###  5.2 Geometric insights

In this section, we give the geometrical explanation for Turing’s instability
of system (3.1), which can help us understand the process of Turing’s losing
stability. Let

$\displaystyle k_{c}=\frac{D}{2D_{u}D_{v}},$ (5.31) $\displaystyle
h(\lambda,L)=D_{u}D_{v}\lambda^{2}-D\lambda+det(A)$
$\displaystyle=D_{u}D_{v}[(\lambda-k_{c})^{2}-L],$ (5.32) $\displaystyle
L_{1}=\min\\{L_{\lambda_{i}},L_{\lambda_{i+1}}\\},$ (5.33)

where

$\displaystyle L_{\lambda_{s}}=(k_{c}-\lambda_{i})^{2},s=i,i+1,$ (5.34)
$\displaystyle\lambda_{i}<k_{c}<\lambda_{i+1},~{}~{}detM_{\lambda_{i}}>0,~{}~{}detM_{\lambda_{i+1}}>0,$
(5.35)

$L$ and $D$ are as in (5.1) and (5.3). If $k_{c}=\frac{D}{2D_{u}D_{v}}$ is a
constant, $\lambda_{i}<k_{c}<\lambda_{i+1}$, based on (4.25) and (5.5),
obviously, $L_{1}$ is the critical parameter reflecting the Turing’s
instability and transition.

CASE 1:$L_{\lambda_{i}}\neq L_{\lambda_{i+1}}$.

Figure 4: A representative graph for Case 1. In this case,
$(\lambda_{i}-k_{c})^{2}=(\lambda_{i+1}-k_{c})^{2}$.

Without loss of generality, taking $L_{1}=L_{\lambda_{k}}$, which means that
$det(M_{\lambda_{k}})<0$ for $L>L_{1}$, in particular,

$\displaystyle\beta_{i}^{(2)}\left\\{\begin{array}[]{ll}<0,~{}if~{}L<L_{1},\\\
=0,~{}if~{}L=L_{1},\\\ >0,~{}if~{}L>L_{1}.\end{array}\right.$ (5.39)
$\displaystyle\beta_{k}^{(j)}(L_{1})<0,~{}for~{}all(k,j)\neq(i,2),$ (5.40)

The transition of case 1 is shown in Fig. 4.

The curves crossing the fixed point $(0,detA)$ in figure 4 is determined by
(5.2). Based on Theorem 5.1, system (3.1) generates Turing’s instability if
and only if there is a spectral point of $-\Delta$ falling into between the
two intersection point of $h(\lambda,L)$ and k-axis. Obviously, Fig. 4 shows
that the spectral point $\lambda_{i}$ of $-\Delta$ is exactly a intersection
point of $h(\lambda,L)$ and k-axis.

CASE 2:$L_{\lambda_{i}}=L_{\lambda_{i+1}}$.

Figure 5: A representative graph for Case 2. In this case,
$(\lambda_{i}-k_{c})^{2}\neq(\lambda_{i+1}-k_{c})^{2}$.

In the same way, we can get

$\displaystyle\beta_{j}^{(2)}\left\\{\begin{array}[]{ll}<0,ifL<L_{1},\\\
=0,ifL=L_{1},\\\ >0,ifL>L_{1},\end{array}\right.~{}for~{}j=i~{}and~{}i+1,$
(5.44) $\displaystyle\beta_{k}^{(h)}(L_{1})<0,for~{}all~{}(k,h)\neq(j,2),$
(5.45)

The transition of this case is shown in Fig. 5, and the
$\lambda_{i_{2}}=\lambda_{k}=\lambda_{i_{1}+1}$ are exactly the two
intersection points of $h(\lambda,L)$ and $\lambda$-axis.

In fact, the curve determined by $h(\lambda,L)$ move down as the increasing of
the parameter $L(D_{u},D_{v})$, and $L_{1}$ is the critical value of the
$L(D_{u},D_{v})$ determined by $h(\lambda,L)$ touching the eigenvalue
$\lambda_{i}$ or $\lambda_{i+1}$.

The two cases mean that Turing’s instability and phase transition of system (
3.1 ) has only two types, the single real eigenvalue transition for the case
that $L_{\lambda_{i}}\neq L_{\lambda_{i+1}}$, and the double real eigenvalue
transition for the case that $L_{\lambda_{i}}=L_{\lambda_{i+1}}$.

###  5.3 Center manifold reduction of Turing’s system

In addition to Section 4.4, center manifold reduction is also a basic tool to
calculate the bifurcated solution in dynamic transition theory, which was
established in [17]. Here we show how this tool can help to derive the
bifurcated solution.

Let $L_{1}=L_{\lambda_{i}}$, the $\xi$ and $\xi_{*}$ be the eigenvectors of
matrix $M_{\lambda_{i}}$ and $M_{\lambda_{i}}^{*}$ respectively, and the
center manifold function is shown as the follows:

$\displaystyle\Phi=\Phi(x):\text{span}{e_{\lambda_{i}}}\rightarrow(\text{span}{e_{\lambda_{i}}})^{\perp}$
(5.46)

For CASE 1, the center manifold reduction system is shown as follows

$\displaystyle\frac{dx}{dt}=\beta_{i}^{(2)}x+\frac{<G(x\xi
e_{\lambda_{i}}+\Phi(x)),\xi^{*}e_{\lambda_{i}}>_{H}}{<\xi
e_{\lambda_{i}},\xi^{*}e_{\lambda_{i}}>_{H}},$ (5.47)

where $\beta_{i}^{(2)}$ is as in (5.39), $e_{\lambda_{i}}$ is a eigenvector of
$-\Delta$ corresponding to $\lambda_{i}$, $\xi$ and $\xi^{*}$ are the
eigenvectors of the matrix $M_{\lambda_{i}}$ and $M_{\lambda_{i}}^{*}$
respectively.

For CASE 2, the center manifold reduction system is shown as follows

$\displaystyle\frac{dx}{dt}=\beta_{i}^{(2)}x+\frac{<G(x\xi
e_{\lambda_{i}}+y\eta
e_{\lambda_{i+1}}+\Phi(x,y)),\xi^{*}e_{\lambda_{i}}>_{H}}{<\xi
e_{\lambda_{i}},\xi^{*}e_{\lambda_{i}}>_{H}},$ (5.48)
$\displaystyle\frac{dy}{dt}=\beta_{i+1}^{(2)}y+\frac{<G(x\xi
e_{\lambda_{i}}+y\eta
e_{\lambda_{i+1}}+\Phi(x,y)),\eta^{*}e_{\lambda_{i+1}}>_{H}}{<\eta
e_{\lambda_{i+1}},\eta^{*}e_{\lambda_{i+1}}>_{H}},$ (5.49)

where $\beta_{i}^{(2)}$ and $\beta_{i+1}^{(2)}$ are as (5.44),
$e_{\lambda_{i}}$ and $e_{\lambda_{i+1}}$ are the eigenvectors of $-\Delta$
corresponding to $\lambda_{i}$ and $\lambda_{i+1}$ respectively, $\xi$ and
$\xi^{*}$ are the eigenvectors of the matrix $M_{\lambda_{i}}$ and
$M_{\lambda_{i}}^{*}$ respectively, $\eta$ and $\eta^{*}$ are the eigenvectors
(a $2\times 1$ vector of functions) of the matrix $M_{\lambda_{i+1}}$ and
$M_{\lambda_{i+1}}^{*}$ respectively, and
$\Phi(x,y):span\\{e_{\lambda_{i}},e_{\lambda_{i+1}}\\}\rightarrow(span\\{e_{\lambda_{i}},e_{\lambda_{i+1}}\\})^{\perp}$
is the center manifold function.

Therefore, we basically reduced an infinite dimensional system (3.1) to an one
or two-dimensional dynamical system, depending on the $i-$th eigenspace of
(3.1). (5.47)-(5.49) total determines all the bifurcated solutions of (3.1) by
expressing the bifurcated solutions locally using a center manifold function
defined in Theorem 4.5.

##  6 Application to the Schnakenberg system

The Schnakenberg system is a well-studied reaction-diffusion systems. It is a
classical example of non-equilibrium thermodynamics resulting in the
establishment of a nonlinear chemical oscillator [22]. It has also been used
to model the spatial distribution of a morphogen, e.g., the distribution of
calcium in the tips and whorl in Acetabularia [5]. As reviewed at the
beginning of this paper, morphogen-based mechanisms have been widely proposed
for tissue patterning, but only recently have there been sufficient
experimental data and adequate modeling for us to begin to understand how
various morphogens interact with cells and emergent patterns [7, 1].

Denote $X$, $A$, $B$ and $Y$ to be four different chemicals, Schnakenberg
considered the following chemical reaction

$\displaystyle X\xrightleftharpoons[k_{1}]{k_{-1}}A,\quad
B\overset{k_{2}}{\rightarrow}Y,\quad 2X+Y\overset{k_{3}}{\rightarrow}X$ (6.1)

###  6.1 Mathematical form of the Schnakenberg system

If concentrations $A$ and $B$ are approximately constants (e.g. $A$ and $B$
are abundant in the system), after proper nondimensionalization and impose
Dirichlet/Neumann boundary condition on a bounded domain
$\Omega\subset\mathbb{R}^{n}$, then the mathematical form of the Schnakenberg
model is given by

$\displaystyle\left\\{\begin{array}[]{ll}\frac{\partial u}{\partial t}=\Delta
u+r(a-u+u^{2}v),\\\ \frac{\partial v}{\partial t}=d\Delta
v+r(b-u^{2}v).\end{array}\right.$ (6.4)

where $u$ and $v$ are all the concentrations of the chemicals, $u$ is
activator, and $v$ is inhibitor. Obviously, the steady state of (6.4) is as
follows

$(u_{0},v_{0})=(a+b,\frac{b}{(a+b)^{2}}).$ (6.5)

Make the following substitution

$u=a+b+u_{1},~{}v=\frac{b}{(a+b)^{2}}+v_{1},$

then equation (6.4) can be rewritten as

$\displaystyle\left\\{\begin{array}[]{ll}\frac{\partial u_{1}}{\partial
t}=\Delta
u_{1}+r(\frac{b-a}{a+b}u_{1}+(a+b)^{2}v_{1}+2(a+b)u_{1}v_{1}-\frac{b}{(a+b)^{2}}u_{1}^{2}+u_{1}^{2}v_{1}),\\\
\frac{\partial v_{1}}{\partial t}=d\Delta
v_{1}-r(\frac{2b}{a+b}u_{1}+(a+b)^{2}v_{1}+2(a+b)u_{1}v_{1}-\frac{b}{(a+b)^{2}}u_{1}^{2}+u_{1}^{2}v_{1}).\end{array}\right.$
(6.8)

By (4.5)-(4.27), equation (6.8) is equivalents to following operator equation

$\displaystyle\left\\{\begin{array}[]{ll}\frac{\partial w}{\partial
t}=L_{\lambda}w+G(w,\lambda),\\\ v(0)=\varphi,\end{array}\right.$ (6.11)

where

$\displaystyle
L_{\lambda}=A_{\lambda}+B_{\lambda},G(w,\lambda):H_{1}\rightarrow H,$ (6.12)
$\displaystyle w=(u_{1},v_{1})^{T},$ (6.13) $\displaystyle
B_{\lambda}u=(\Delta u_{1},d\Delta v_{1}),$ (6.14) $\displaystyle
A_{\lambda}v=\left(\begin{array}[]{cc}r\frac{b-a}{a+b}&r(a+b)^{2}\\\
-r\frac{2b}{a+b}&-r(a+b)^{2}\end{array}\right)w,$ (6.17) $\displaystyle
G(w,\lambda)=(g_{1}(w,r),g_{2}(w,r))$ (6.18) $\displaystyle
g_{1}(w,r)=r(2(a+b)u_{1}v_{1}-\frac{b}{(a+b)^{2}}u_{1}^{2}+u_{1}^{2}v_{1}),$
(6.19) $\displaystyle
g_{2}(w,r)=-r(2(a+b)u_{1}v_{1}-\frac{b}{(a+b)^{2}}u_{1}^{2}+u_{1}^{2}v_{1}),$
(6.20) $\displaystyle\lambda=(1,\alpha,r)\in
R_{3}^{+}=\\{(x_{1},x_{2},x_{3})\in R^{3}:x_{i}\geq 0,i=1,2,3\\}.$

and also

$\displaystyle a_{11}=r\frac{b-a}{a+b},a_{12}=r(a+b)^{2},$ (6.21)
$\displaystyle a_{21}=-r\frac{2b}{a+b},a_{22}=-r(a+b)^{2}.$ (6.22)

###  6.2 A necessary and sufficient condition for Turing’s instability

For system (6.4), based on (5.1) and (5.3), we have

$\displaystyle D=rd\frac{b-a}{a+b}-r(a+b)^{2},$ (6.23) $\displaystyle
det(A)=r^{2}(a+b)^{2},$ (6.24) $\displaystyle
L=-\frac{r^{2}(a+b)^{2}}{d}+\frac{D^{2}}{4d^{2}},$ (6.25)

Let

$\beta_{\pm}=\frac{D\pm\sqrt{D^{2}-4d~{}det(A)}}{2d},$ (6.26)

based on Theorem (5.1) , we get the following conclusion

######  Theorem 6.1.

Assume $\frac{b-a}{a+b}-(a+b)^{2}<0$
(1) if $-\frac{\sqrt{4d}}{r(a+b)}<d\frac{b-a}{a+b}-(a+b)^{2}<0$, i.e.

$\displaystyle\max{(d,1)}(b-a)<(a+b)^{3}<(d+\frac{\sqrt{4d}}{r})(b-a)$ (6.27)

hen system (6.4) is Turing stable.
(2) if $d\frac{b-a}{a+b}-(a+b)^{2}>\frac{\sqrt{4d}}{r(a+b)}$ or
$d\frac{b-a}{a+b}-(a+b)^{2}<-\frac{\sqrt{4d}}{r(a+b)}$, i.e.

$\displaystyle(b-a)<(a+b)^{3}<(d-\frac{\sqrt{4d}}{r})(b-a)$ (6.28)
$\displaystyle\text{or~{}}\max{(1,d+\frac{\sqrt{4d}}{r})}(b-a)<(a+b)^{3}$
(6.29)

then system (6.4) generates Tuing’s instability if and only if there exists an
eigenvalue $\lambda_{k_{0}}$ of $-\Delta$ such that
$\beta_{-}<\lambda_{k_{0}}<\beta_{+}$.

###  6.3 The critical parameter of the Schnakenberg system

In fact, $d$ is an adjustable parameter for Schnakenberg system. That means
that the critical parameter is determined by $d$ . In the following, we will
give the critical parameter $d_{0}$.

Due to the method introduced in Section 4.4, by directed calculation we get

$\displaystyle
h(\lambda)=d\lambda^{2}-r(\frac{b-a}{b+a}d-(a+b)^{2})\lambda+r^{2}(a+b)^{2},$
(6.30) $\displaystyle k_{c}=\frac{r(b-a)d-r(a+b)^{3}}{2d(a+b)}.$ (6.31)

Let

$D_{v}^{0}=\frac{(4b^{2}+4ab)(a+b)^{2}+(a+b)^{2}\sqrt{(4b^{2}+4ab)^{2}+a^{2}-b^{2}}}{(b-a)^{2}},$
(6.32) $\displaystyle
k_{c}^{0}=r\frac{(b-a)D_{v}^{0}-(b+a)^{3}}{2D_{v}^{0}(b+a)},$ (6.33)
$\displaystyle
d_{0}^{\lambda_{k}}=\frac{-r^{2}(a+b)^{3}-r(a+b)^{3}\lambda_{k}}{(a+b)\lambda_{k}^{2}-r(b-a)},$
(6.34) $\displaystyle
d_{0}^{\lambda_{k+1}}=\frac{-r^{2}(a+b)^{3}-r(a+b)^{3}\lambda_{k+1}}{(a+b)\lambda_{k+1}^{2}-r(b-a)},$
(6.35) $\displaystyle
d_{0}=\min\\{d_{0}^{\lambda_{k}},d_{0}^{\lambda_{k+1}}\\},$ (6.36)

where $\lambda_{k}$ and $\lambda_{k+1}$ are the eigenvalues of the operator
$-\Delta$ such that

$\displaystyle\lambda_{k}<k_{c}^{0}<\lambda_{k+1},~{}h(\lambda_{k})>0,~{}h(\lambda_{k+1})>0,$
(6.37)

Based on Theorem (5.2) in section 2, then we get the following corollary.

######  Corollary 2.

Let $d_{0}$ be as in (6.36), then $d_{0}$ is the critical parameter reflecting
Turing’s instability and phase transition, i.e, if $d>d_{0}$, then Turing’s
instability appears and dynamic transition occurs.

###  6.4 Phase Transition of the Schnakenberg system

####  6.4.1 Single real eigenvalue transition of the Schnakenberg system

Based on the dynamic bifurcation theory in [17], without loss of generality,
assume that $\lambda_{i_{1}}=\lambda_{1}$, then the center manifold reduction
equation for system (6.4) is given by

$\frac{dx}{dt}=\beta_{\lambda_{1}}^{(2)}x+\frac{1}{\langle\xi
e_{\lambda_{1}},\xi^{*}e_{\lambda_{1}}\rangle}_{H}\langle G(x\xi
e_{\lambda_{1}}+\Phi(x)),\xi^{*}e_{\lambda_{1}}\rangle_{H}$ (6.38)

where

$M_{\lambda_{1}}\xi=0,M_{\lambda_{1}}^{*}\xi^{*}=0,$ (6.39)
$M_{\lambda_{1}}=\left(\begin{array}[]{cc}r\frac{b-a}{a+b}-\lambda_{1}^{2}&r(a+b)^{2}\\\
-r\frac{2b}{a+b}&-r(a+b)^{2}-\lambda_{1}^{2}\end{array}\right),M_{\lambda_{1}}^{*}=M_{\lambda_{1}}^{T}.$
(6.40)
$\xi=(\xi_{1},\xi_{2})=(a+b,\frac{(a+b)\lambda_{1}^{2}-(b-a)r}{r(a+b)^{2}}),$
(6.41)
$\xi^{*}=(\xi_{1}^{*},\xi_{2}^{*})=(\frac{2br(a+b)}{(b-a)r-(a+b)\lambda_{1}^{2}},a+b).$
(6.42)

Note

$h=\bigg{(}\frac{2br(a+b)^{2}}{(b-a)r-(a+b)\lambda_{1}^{2}}+\frac{(a+b)\lambda_{1}^{2}-(b-a)r}{r(a+b)}\bigg{)}\int_{\Omega}e_{\lambda_{1}}^{2}dx,$

thus we get the second order term and the third order term as follows.

$\displaystyle G_{2}=r(-Hu^{2}+Muv,Hu^{2}-Muv),$ (6.43) $\displaystyle
G_{3}=r(u^{2}v,-u^{2}v),$ (6.44) $\displaystyle H=\frac{b}{(a+b)^{2}},$ (6.45)
$\displaystyle M=2(a+b).$ (6.46)

We can also get

$\displaystyle\frac{1}{h}\langle
G_{2}(x\xi_{1}e_{\lambda_{1}},x\xi_{2}e_{\lambda_{1}}),\xi^{*}e_{\lambda_{1}}\rangle_{H}$
$\displaystyle=\frac{r}{h}\bigg{(}(\xi_{2}^{*}-\xi_{1}^{*})H\xi_{1}^{2}\int_{\Omega}e_{\lambda_{1}}^{3}dx+(\xi_{2}^{*}-\xi_{1}^{*})\xi_{1}\xi_{2}M\int_{\Omega}e_{\lambda_{1}}^{3}dx\bigg{)}x^{2}$
$\displaystyle=Px^{2}.$ (6.47)

Hence, the center manifold reduction equation (6.38) can be written as

$\frac{dx}{dt}=\beta_{\lambda_{1}}^{(2)}x+Px^{2}.$ (6.48)

Suppose

$\int_{\Omega}e_{\lambda_{1}}^{3}dx\neq 0.$ (6.49)

We have the following Theorem.

######  Theorem 6.2.

Let $\lambda_{1}^{2}\neq r$ and $\int e_{\lambda_{1}}^{3}dx\neq 0$, then the
system (6.4) has a transition at $d=d_{0}$, which is mixed transition. In
particular, the system bifurcates on each side of $d=d_{0}$ to a unique branch
$w^{d}$ of steady state solutions, such that the following assertions hold
true:
(1) When $d<d_{0}$, the bifurcated solution $w^{d}$ is a saddle, and the
stable manifold $\Gamma_{d}^{1}$ of $w^{d}$ separates the space H into two
disjoint open sets $U^{d}$ and$V^{d}$, such that $v=0\in V^{d}$ is an
attractor, and the orbits of (3.24 ) in $u^{d}$are far from $v=0$.
(2) When $d>d_{0}$, the stable manifold $\Gamma_{d}^{0}$ of $v=0$ separates
the neighborhood O of $u=0$ into two disjoint open sets $O_{1}^{d}$ and
$O_{2}^{d}$, such that the transition is jump in $O_{1}^{0}$, and is
continuous in $O_{2}^{d}$. The bifurcated solution $w^{d}\in O_{2}^{d}$is an
attractor such that for any $\phi\in O_{2}^{d}$,we have

$\lim_{t\rightarrow\infty}\|w(t,\phi)-w^{d}\|_{H}=0$ (6.50)

where $w(t,\phi)$ is the solution of (6.4) with $w(0,\phi)=\phi.$
(3) The bifurcated solution $w^{d}$ can be expressed as

$w^{d}=-\frac{\beta_{\lambda_{1}}^{(2)}}{P}\xi e_{\lambda_{1}}(x)$ (6.51)

$\xi$ and $\beta_{\lambda_{1}}^{(2)}$ are shown above.

###### Proof.

The reduction system (6.38) equivalents to

$\frac{dx}{dt}=\beta_{\lambda_{1}}^{(2)}x+Px^{2}.$ (6.52)

Due to the PES condition as follows:

$\displaystyle\beta_{\lambda_{1}}^{(2)}\left\\{\begin{array}[]{ll}<0,~{}if~{}d<d_{0};\\\
=0,~{}if~{}d=d_{0};\\\ >0,~{}if~{}d>d_{0};\end{array}\right.$ (6.56)
$\displaystyle\beta_{k}^{(j)}(d_{0})<0,~{}for~{}all(k,j)\neq(\lambda_{1}2).$
(6.57)

Therefore, all the above conclusions hold true.

The transition behavior of 6.56 is shown in Fig. 6.

$o$$v^{d}$$U_{2}^{d}$$U_{1}^{d}$$d<d_{0}$

$o$$v^{d}$$O_{1}^{d}$$O_{2}^{d}$$d>d_{0}$

Figure 6: Transition diagram for Theorem 6.2.

∎

If (6.49) is not true, i.e,

$\int_{\Omega}e_{\lambda_{1}}^{3}dx=0.$ (6.58)

We introduce the following parameter

$Q=\frac{1}{h}(\xi_{1}^{*}-\xi_{2}^{*})Mr\xi_{1}\int_{\Omega}\psi_{2}e_{\lambda_{1}}^{2}{\rm
d}x-2(\xi_{1}^{*}-\xi_{2}^{*})rH\xi_{1}\int_{\Omega}\psi_{1}e_{\lambda_{1}}^{2}{\rm
d}x$ (6.59)

Where$\psi=(\psi_{1},\psi_{2})$ satisfied the following equation

$\displaystyle\Delta\psi_{1}+r\frac{b-a}{b+a}\psi_{1}+r(a+b)^{2}\psi_{2}=-rH\xi_{1}^{2}e_{\lambda_{1}}^{2}+rM\xi_{1}\xi_{2}e_{\lambda_{1}}^{2},$
(6.60)
$\displaystyle\Delta\psi_{1}-r\frac{2b}{b+a}\psi_{1}-r(a+b)^{2}\psi_{2}=rH\xi_{1}^{2}e_{\lambda_{1}}^{2}-rM\xi_{1}\xi_{2}e_{\lambda_{1}}^{2}.$
(6.61)

By the Fredholm Alternative Theorem, under the condition (6.58), the equations
(6.60)-(6.61) have a unique solutions. Moreover, we have the following
conclusion.

######  Theorem 6.3.

Let $\lambda_{1}^{2}\neq r$, $\int e_{\lambda_{1}}^{3}dx=0$, and $Q$ be the
number given by (6.59), then the system (6.4) has a transition at $d=d_{0}$,
and the transition of (6.4) at $d=d_{0}$ is continuous if $Q<0$, jump if
$Q>0$. The following assertions hold true:
(1)If $Q>0$,(6.4) has no bifurcation when $d>d_{0}$, and has exact two
bifurcated solutions $w_{+}^{d}$and $w_{-}^{d}$ which are saddles when
$d<d_{0}$. Moreover, the stable manifolds $\Gamma_{+}^{d}$ and
$\Gamma_{-}^{d}$ of the two bifurcated solutions divide H into three disjoint
open sets $U_{+}^{d}$, $U_{0}^{d}$, $U_{-}^{d}$such that $v=0\in U_{0}^{d}$is
an attractor, and the orbits in $U_{\pm}^{d}$ are far from $v=0$.
(2)If $Q<0$,(6.4) has no bifurcation when $d<d_{0}$, and has exact two
bifurcated solutions $w_{+}^{d}$ and $w_{-}^{d}$ when $d>d_{0}$ , which are
attractors. In addition, there is a neighborhood $O\subseteq H$of $w=0$, such
that the stable manifold $\gamma$ of $w=0$ divides $O$ into two disjoint open
sets $U_{+}^{d}$ and $U_{-}^{d}$ such that $w_{+}^{d}\in U_{+}^{d}$,
$w_{-}^{d}\in U_{-}^{d}$, and$w_{\pm}^{d}$ attracts $U_{\pm}^{d}$;
(3)The bifurcated solutions $w_{\pm}^{d}$ can be expressed as

$w_{\pm}^{d}=\pm(-\frac{\beta_{\lambda_{1}}^{(2)}}{Q})^{\frac{1}{2}}\xi
e_{\lambda_{1}}(x)$ (6.62)

$\xi$ and $\beta_{\lambda_{1}}^{(2)}$ is shown above.

###### Proof.

To prove this result, we need to calculate the center manifold function
$\Phi(x)$. By the procedures in Section 3 of [9] and Appendix A of [13],
$\Phi(x)$ satisfies:

$L_{\lambda}\Phi=-P_{2}G(x\xi e_{\lambda_{1}})$ (6.63)

where $P_{2}:H\rightarrow E_{2}$ is the canonical projection, $L_{\lambda}$ is
as in (6.12), $e_{\lambda_{1}}$ and $e_{\lambda_{1}}^{*}$ are given by (6.58),
and

$E_{2}=\\{v\in H:(v,e_{\lambda_{1}}^{*})=0\\}$ (6.64)

we can see that

$\displaystyle G(x\xi
e_{\lambda_{1}})=(-rH\xi_{1}^{2}e_{\lambda_{1}}^{2}x^{2}+rM\xi_{1}\xi_{2}e_{\lambda_{1}}^{2}x^{2},$
(6.65) $\displaystyle
rH\xi_{1}^{2}e_{\lambda_{1}}^{2}x^{2}-rM\xi_{1}\xi_{2}e_{\lambda_{1}}^{2}x^{2}).$
(6.66)

Let

$\Phi=\psi x^{2}+0(2),\psi=(\psi_{1},\psi_{2})\in H.$ (6.67)

By (6.58), $(e_{\lambda_{1}}^{2},-e_{\lambda_{1}}^{2})\in E_{2}$. Hence, it
follows from (6.63) and (6.67) that

$L_{\lambda}\psi=(-r^{2}H\xi_{1}^{2}x^{2}+rM\xi_{1}\xi_{2}x^{2})(e_{\lambda_{1}}^{2},-e_{\lambda_{1}}^{2})$
(6.68)

By direct calculation, we obtain the following

$\displaystyle<G(x\xi e_{\lambda_{1}}+\psi x_{2}),\xi^{*}e_{\lambda_{1}}>_{H}$
(6.69)
$\displaystyle=((Mr\xi_{1}\psi_{2}e_{\lambda_{1}}-2rH\xi_{1}\psi_{1})x^{3}+0(3)$
(6.70)
$\displaystyle-(Mr\xi_{1}\psi_{2}e_{\lambda_{1}}-2rH\xi_{1}\psi_{1})x^{3}+0(3))=Qx^{3},$
(6.71)

where

$Q=\frac{1}{h}(\xi_{1}^{*}-\xi_{2}^{*})Mr\xi_{1}\int_{\Omega}\psi_{2}e_{\lambda_{1}}^{2}{\rm
d}x-2(\xi_{1}^{*}-\xi_{2}^{*})rH\xi_{1}\int_{\Omega}\psi_{1}e_{\lambda_{1}}^{2}{\rm
d}x$ (6.72)

In the end, we get the center manifold reduction system as follows:

$\frac{dx}{dt}=\beta_{\lambda_{1}}^{(2)}x+Qx^{3}.$ (6.73)

Whose steady state solutions are

$x_{\pm}=\pm(-\frac{\beta_{\lambda_{1}}^{(2)}}{Q})^{\frac{1}{2}}$ (6.74)

Therefore, following procedures found in Section 3 of [9] and Appendix A of
[13], all the above conclusions hold true.

Fig. 7 and 8 are corresponding transition diagrams for the cases when $Q>0$
and $Q<0$ respectively.

$o$$v_{+}^{d}$$v_{-}^{d}$$U_{0}^{d}$$U_{+}^{d}$$U_{-}^{d}$$d<d_{0}$

$o$$d>d_{0}$

Figure 7: Transition diagram when $Q>0$.

$o$$d<d_{0}$

$o$$v_{+}^{d}$$v_{-}^{d}$$U_{0}^{d}$$U_{+}^{d}$$U_{-}^{d}$$d>d_{0}$

Figure 8: Transition diagram when $Q<0$.

∎

If $\lambda_{1}^{2}=r$, for the case that $P=Q=0$, we have the following
conclusion.

######  Theorem 6.4.

The attractor $v=0$ is stable if $d<d_{0}$, and unstable if $d<d_{0}$ for
system (6.4), there is no solution bifurcated from the critical parameter
$d=d_{0}$.

###### Proof.

The condition $\lambda_{1}^{2}=r$ is equavalent to $\xi_{1}^{*}=\xi_{2}^{*}$,
so $P=Q=0$ . The conclusion can be simply obtained from the following
reduction system

$\frac{dx}{dt}=\beta_{\lambda_{1}}^{(2)}x$ (6.75)

The proof is complete. ∎

####  6.4.2 The double real eigenvalues transition of the Schnakenberg system.

For the PES (6.56), there exists real double eigenvalues transition for
$d>d_{0}$. let $i_{2}=i_{1}+1$ and $i=i_{1}$, then the dynamic transition
behavior of the system can be dictated by the center manifold reduction
equations for this case as follows.

$\displaystyle\left\\{\begin{array}[]{ll}\frac{dx}{dt}=\beta_{\lambda_{i_{1}}}^{2}x+\frac{1}{h_{1}}\langle
G_{2}(x\xi e_{\lambda_{i_{1}}}+y\eta
e_{\lambda_{i_{2}}}),\xi^{*}e_{\lambda_{i_{1}}}\rangle_{H},\\\
\frac{dy}{dt}=\beta_{\lambda_{i_{2}}}^{2}y+\frac{1}{h_{2}}\langle G_{2}(x\xi
e_{\lambda_{i_{1}}}+y\eta
e_{\lambda_{i_{2}}}),\eta^{*}e_{\lambda_{i_{2}}}\rangle_{H},\end{array}\right.$
(6.78)

where

$\xi=(\xi_{1},\xi_{2})=(a+b,\frac{(a+b)\lambda_{i_{1}}^{2}-(b-a)r}{r(a+b)^{2}})^{T},$
(6.79)
$\eta=(\eta_{1},\eta_{2})=(a+b,\frac{(a+b)\lambda_{i_{2}}^{2}-(b-a)r}{r(a+b)^{2}})^{T},$
(6.80)
$\xi=(\xi_{1}^{*},\xi_{2}^{*})=(\frac{2br(a+b)}{(b-a)r-(a+b)\lambda_{i_{1}}^{2}},a+b)^{T},$
(6.81)
$\eta^{*}=(\eta_{1}^{*},\eta_{2}^{*})=(\frac{2br(a+b)}{(b-a)r-(a+b)\lambda_{i_{2}}^{2}},a+b)^{T},$
(6.82)
$h_{1}=\bigg{(}\frac{2br(a+b)^{2}}{(b-a)r-(a+b)\lambda_{i_{1}}}+\frac{(a+b)\lambda_{i_{1}}-(b-a)r}{r(a+b)}\bigg{)}\int_{\Omega}e_{\lambda_{i_{1}}}^{2}dx,$
(6.83)
$h_{2}=\bigg{(}\frac{2br(a+b)^{2}}{(b-a)r-(a+b)\lambda_{i_{2}}}+\frac{(a+b)\lambda_{i_{2}}-(b-a)r}{r(a+b)}\bigg{)}\int_{\Omega}e_{\lambda_{i_{2}}}^{2}dx.$
(6.84)

Let

$\displaystyle s_{1}=-Hr\xi_{1}^{2}+Mr\xi_{1}\xi_{2},$ (6.85) $\displaystyle
s_{2}=-2Hr\xi_{1}\eta_{1}+Mr(\xi_{1}\eta_{2}+\xi_{2}\eta_{1}),$ (6.86)
$\displaystyle s_{3}=-2Hr\eta_{1}^{2}+Mr\eta_{1}\eta_{2},$ (6.87)

where

$H=\frac{b}{(a+b)^{2}},~{}~{}~{}~{}~{}M=2(a+b)$ (6.88)

then,

$\displaystyle G_{2}(x\xi e_{\lambda_{i_{1}}}+y\eta
e_{\lambda_{i_{2}}})=(M_{1},-M_{1}),$ (6.89) $\displaystyle
M_{1}=s_{1}e_{\lambda_{i_{1}}}^{2}x^{2}+s_{2}e_{\lambda_{i_{1}}}e_{\lambda_{i_{2}}}xy+s_{3}e_{\lambda_{i_{2}}}^{2}y^{2}.$
(6.90)

Denote

$\displaystyle
a_{20}=\frac{s_{1}(\xi_{1}^{*}-\xi_{2}^{*})}{h_{1}}\int_{\Omega}e_{\lambda_{i_{1}}}^{3}dx,$
(6.91) $\displaystyle
a_{11}=\frac{s_{2}(\xi_{1}^{*}-\xi_{2}^{*})}{h_{1}}\int_{\Omega}e_{\lambda_{i_{1}}}^{2}e_{\lambda_{i_{2}}}dx,$
(6.92) $\displaystyle
a_{02}=\frac{s_{3}(\xi_{1}^{*}-\xi_{2}^{*})}{h_{1}}\int_{\Omega}e_{\lambda_{i_{2}}}^{2}e_{\lambda_{i_{1}}}dx,$
(6.93) $\displaystyle
b_{20}=\frac{s_{1}(\eta_{1}^{*}-\eta_{2}^{*})}{h_{2}}\int_{\Omega}e_{\lambda_{i_{2}}}^{3}dx,$
(6.94) $\displaystyle
b_{11}=\frac{s_{2}(\eta_{1}^{*}-\eta_{2}^{*})}{h_{2}}\int_{\Omega}e_{\lambda_{i_{2}}}^{2}e_{\lambda_{i_{1}}}dx,$
(6.95) $\displaystyle
b_{02}=\frac{s_{3}(\eta_{1}^{*}-\eta_{2}^{*})}{h_{2}}\int_{\Omega}e_{\lambda_{i_{1}}}^{2}e_{\lambda_{i_{2}}}dx,$
(6.96)

then, (6.4) can be written as

$\displaystyle\left\\{\begin{array}[]{ll}\frac{dx}{dt}=\beta_{\lambda_{i_{1}}}^{2}x+a_{20}x^{2}+a_{11}xy+a_{02}y^{2},\\\
\frac{dy}{dt}=\beta_{\lambda_{i_{2}}}^{2}y+b_{20}x^{2}+b_{11}xy+b_{02}y^{2}.\end{array}\right.$
(6.99)

Denote by

$\displaystyle A_{1}=\beta_{\lambda_{i_{1}}1}+2a_{20}x_{0}+a_{11}y_{0},$
(6.100) $\displaystyle A_{2}=2y_{0}a_{02}+a_{11}y_{0},$ (6.101) $\displaystyle
A_{3}=2b_{02}x_{0}^{2}+b_{11}y_{0},$ (6.102) $\displaystyle
A_{4}=\beta_{\lambda_{i_{2}}1}+2b_{20}y_{0}+b_{11}x_{0}.$ (6.103)

Hence, the transition for (6.4) is equivalent to the transition for (6.99).
The transition theorem is stated as follows.

######  Theorem 6.5.

The system (6.4) has a stable steady state $(0,0)$ for $d<d_{0}$, but
bifurcates a new steady state $(x_{0},y_{0})$ for $d>d_{0}$, if and only if
the following conditions hold true

$\displaystyle A_{1}+A_{4}<0,A_{1}A_{4}-A_{2}A_{3}>0,$ (6.104)

and there is a stable attractor $u^{d}$ of the system bifurcated from
$d>d_{0}$ as follows

$u^{d}=x_{0}\xi e_{\lambda_{i_{1}}}+y_{0}\eta e_{\lambda_{i_{2}}}$

###### Proof.

The following matrix is the Jacobian of the system at $(x_{0},y_{0})$

$\left(\begin{array}[]{cc}A_{1}&A_{2}\\\ A_{3}&A_{4}\end{array}\right)$
(6.105)

The condition (6.104) tell us that the $(x_{0},y_{0})$ is a stable attractor,
which can be written as $u^{d}=x_{0}\xi e_{\lambda_{i_{1}}}+y_{0}\eta
e_{\lambda_{i_{2}}}$. ∎

###  6.5 An illustrative example

Taking $\Omega$ to be the rectangle $[0,10]\times[0,5]$, if we consider
Dirichlet boundary condition, then we have

$\displaystyle e_{m,n}(x)=cos(\frac{m\pi x}{10})cos(\frac{n\pi y}{5}),$
(6.106)
$\displaystyle\lambda_{m,n}(x)=(\frac{m\pi}{10})^{2}+(\frac{n\pi}{5})^{2},$
(6.107)

and $m,n=0,1,2,3,4\ldots$.

By directly calculating, we can get

$\displaystyle\begin{aligned} &\lambda_{0,1}=0.3948,\lambda_{1,0}=0.0987,\\\
&\lambda_{2,0}=0.3948,\lambda_{1,1}=0.4935,\\\
&\lambda_{1,2}=1.6778,\lambda_{0,2}=1.5791,\end{aligned}$ (6.108)

Figure 9: Plot of the critical parameter $d_{0}^{\lambda}$ as a function of
$\lambda$ when (a) $a=0.1,b=0.5,r=1$, (b) $a=1,b=0.5,r=1$. The red stars
correspond to the eigenvalues of $-\Delta$ in $\Omega$ as calculated in
(6.108)

For system (6.4), we consider two sets of parameters differing only in $a$,
specifically

(a) $\displaystyle a=0.1,b=0.5,r=1,$ (6.109) (b) $\displaystyle
a=1,b=0.5,r=1.$ (6.110)

then for case (a), we get (calculated using MATLAB)

$\displaystyle D_{v}^{0}=5.1648,k_{c}^{0}=0.2985$ (6.111)

since $\lambda_{i}$ in Theorem 5.2 is $\lambda_{1,0}=0.0987$ and
$\lambda_{0,1}=0.3948$, for case (a), $\lambda_{1,0}<k_{c}^{0}<\lambda_{0,1}$
is satisfied and the critical parameter $d_{0}=$ is given by

$\displaystyle
d_{0}=\min(d_{0}^{\lambda_{1,0}},d_{0}^{\lambda_{0,1}})=0.6157.$ (6.112)

Therefore, Turing’s instability and the transition of the system (6.4) with
coefficients (6.109) occurs when diffusion rate of $v$ is greater than
$d_{0}=0.6157$. The numerical simulation of this case is illustrated in Fig.
10 below.

Figure 10: Numerical simulation of the Schnakenberg system for case (a) when
$d=1$ under Dirichlet boundary condition $(u,v)=(u_{0},v_{0})$ on
$\partial\Omega$. The images shows temporal snapshots of spatial distributions
of chemical $V$. In this case, the homogeneous steady state losses its
stability and generates dynamics patterns (periodic solutions).

For case (b), however, we get

$\displaystyle D_{v}^{0}=55.1025,k_{c}^{0}=-0.1871$ (6.113)

since $k_{c}^{0}<0<\inf_{k\geq 0}{(\lambda_{k})}$, the assumption (6.37) does
not apply. However, because $D=-2.5833,L=-0.5816<0$, the system must be
stable. The numerical simulation of this case is illustrated in Fig. 11 below.

Figure 11: Numerical simulation of the Schnakenberg system for case (b) when
$d=1$ under Dirichlet boundary condition $(u,v)=(u_{0},v_{0})$ on
$\partial\Omega$ and was purturbed. The images shows temporal snapshots of
spatial distributions of chemical $V$. In this case, no Turing’s instability
occurs.

##  7 Discussion

In this paper, we have studied the stability and dynamic transition property
of Turing’s systems, using the new dynamical transition theory developed in
[9, 13]. Turing’s systems are famous for the so-called diffusion-driven
instabilities, resulting emerging patterns as a result of bifurcations from
the uniform steady state. Our analysis showed that only two bifurcation
parameters are enough to dictate all the dynamics of the Turing’s system near
a steady state, as shown in Section 5. Besides, with the help of center
manifold reduction procedure, we are able to locally derive all the dynamical
behavior of a infinite dimensional system near a steady state.

Turing’s instability is a widely studied topic in nonlinear science, lots of
works had been done with respect to this. Most related works, however, mainly
address on analyzing asymptotic behavior or numerical simulations, and mostly
focus on a specific Turing’s system, e.g. the Brusselator model, the Gray-
Scott system, the Schnakenberg system etc., as introduced in Section 1. Our
work, on the other hand, started from the general model (3.1) and its
analytical properties. From these analytical properties, by calculating
critical parameters and bifurcated solutions, we found that Turing’s
instability can be recognized as a critical phenomena, and consequently can be
studied in a similar fashion as other physical phenomena like
superconductivity, crystallization etc [12, 13]. Dynamic transition theory is
not only a fundamental method of studying critical phenomena, but also a
powerful tool of studying dynamics of a PDE system.

The method introduced in this paper is systematic, and can be applied to
various different kinds of two-component reaction diffusion systems (Turing’s
systems). In the future, we propose to build up a general method to study more
complicated phase transition behavior of a reaction diffusion system, for
example, second order phase transition and its biological and chemical
implications. In general, chaotic behavior will take place for a physical
system underwent second and higher order phase transitions under modern
classification [8, 24, 18, 13], and fruitful critical phenomena structures are
found in numerous experiments for each order of transition. In the case of
finite dimension, this kind of problems has been well studied from
mathematical point of view [27, 6]. Since most physical, chemical and
biological phenomena are modeled by PDEs, accompanying with infinite dimension
operators, constructing a systematic theory of studying second and higher
order phase transition on infinite dimensional dynamical systems is of great
importance for linking experimental results and mathematical theory, and that
will be our next goal.

## References

  * [1] AAM Arafa, SZ Rida, and H Mohamed.  Approximate analytical solutions of schnakenberg systems by homotopy analysis method.  Applied Mathematical Modelling, 36(10):4789–4796, 2012.
  * [2] Stephen H. Davis.  On the principle of exchange of stabilities.  Proc. R. Soc. Lond. A, 310(1502):341–358, 1969.
  * [3] Richard J. Field and Richard M. Noyes.  Oscillations in chemical systems. iv. limit cycle behavior in a model of a real chemical reaction.  The Journal of Chemical Physics, 60(5):1877–1884, 1974.
  * [4] A. Gierer and H. Meinhardt.  A theory of biological pattern formation.  Kybernetik, 12(1):30–39, Dec 1972.
  * [5] Brian C Goodwin and LEH Trainor.  Tip and whorl morphogenesis in acetabularia by calcium-regulated strain fields.  Journal of theoretical biology, 117(1):79–106, 1985.
  * [6] John Guckenheimer and Philip Holmes.  Nonlinear oscillations, dynamical systems, and bifurcations of vector fields, volume 42.  Springer Science & Business Media, 2013.
  * [7] JB Gurdon and P-Y Bourillot.  Morphogen gradient interpretation.  Nature, 413(6858):797, 2001.
  * [8] John Michael Kosterlitz and David James Thouless.  Ordering, metastability and phase transitions in two-dimensional systems.  Journal of Physics C: Solid State Physics, 6(7):1181, 1973.
  * [9] Tian Ma and Shouhong Wang.  Bifurcation Theory and Applications, volume 53.  World Scientific Series on Nonlinear Science Series A, 2005.
  * [10] Tian Ma and Shouhong Wang.  Dynamic bifurcation of nonlinear evolution equations.  Chinese Annals of Mathematics, 26(02):185–206, 2005.
  * [11] Tian Ma and Shouhong Wang.  Phase transitions for belousov-zhabotinsky reactions.  Mathematical Methods in the Applied Sciences, 34(11):1381–1397, July 2011.
  * [12] Tian Ma and Shouhong Wang.  Mathematical Principles of Theoretical Physics.  Science Press, Beijing, 524pp, 2015.
  * [13] Tian Ma and Shouhong Wang.  Phase Transition Dynamics.  Springer, New York, NY, 2015.
  * [14] Philip K. Maini., Thomas E. Woolley, Ruth E. Baker, Eamonn A. Gaffney, and S. Seirin Lee.  Turing’s model for biological pattern formation and the robustness problem.  Interface Focus, 2(4):487–496, August 2012.
  * [15] Hans Meinhardt.  Generation of biological patterns and form: some physical, mathematical, and logical aspects.  Prog Biophys Mol Biol, 37(1):1–47, 1981.
  * [16] James D Murray.  Mathematical biology. II Spatial models and biomedical applications $\\{$Interdisciplinary Applied Mathematics V. 18$\\}$.  Springer-Verlag New York Incorporated New York, 2001.
  * [17] Wei-Ming Ni.  Diffusion, cross-diffusion, and their spike-layer steady states.  Notices of the AMS, 45(1):9–18, 1998.
  * [18] Michael I Ojovan.  Ordering and structural changes at the glass–liquid transition.  Journal of Non-Crystalline Solids, 382:79–86, 2013.
  * [19] Alexander Pechenkin.  B p belousov and his reaction.  Journal of Biosciences, 34(3):365–371, Sep 2009.
  * [20] I. Prigogine and G. Nicolis.  Self-Organisation in Nonequilibrium Systems: Towards A Dynamics of Complexity, pages 3–12.  Springer Netherlands, Dordrecht, 1985.
  * [21] Razvan A. Satnoianu, Michael Menzinger, and Philip K. Maini.  Turing instabilities in general systems.  Journal of Mathematical Biology, 41(6):493–512, Dec 2000.
  * [22] J Schnakenberg.  Simple chemical reaction systems with limit cycle behaviour.  Journal of theoretical biology, 81(3):389–400, 1979.
  * [23] Junping Shi.  Bifurcation in infinite dimensional spaces and applications in spatiotemporal biological and chemical models.  Frontiers of Mathematics in China, 4(3):407–424, Sep 2009.
  * [24] Ricard V Solé, Susanna C Manrubia, Bartolo Luque, Jordi Delgado, and Jordi Bascompte.  Phase transitions and complex systems: Simple, nonlinear models capture complex systems at the edge of chaos.  Complexity, 1(4):13–26, 1996.
  * [25] Roger Temam.  Infinite-Dimensional Dynamical Systems in Mechanics and Physics.  Applied Mathematical Sciences. Springer, 1997.
  * [26] Alan Turing.  The chemical basis of morphogenesis.  Philosophical Transactions of the Royal Society of London B, 237(641):37–72, 1952.
  * [27] Stephen Wiggins.  Introduction to applied nonlinear dynamical systems and chaos, volume 2.  Springer Science & Business Media, 2003.
  * [28] Arthur T Winfree.  The prehistory of the belousov-zhabotinsky oscillator.  Journal of Chemical Education, 61(8):661, 1984.
  * [29] Thomas E. Woolley, Ruth E. Baker, Eamonn A. Gaffney, and Philip K. Maini.  Power spectra methods for a stochastic description of diffusion on deterministically growing domains.  Phys. Rev. E, 84:021915, Aug 2011.
  * [30] Thomas E. Woolley, Ruth E. Baker, Eamonn A. Gaffney, and Philip K. Maini.  Stochastic reaction and diffusion on growing domains: Understanding the breakdown of robust pattern formation.  Phys. Rev. E, 84:046216, Oct 2011.