WOLFRAM|DEMONSTRATIONS PROJECT

Turing Pattern in a Reaction-Diffusion System

​
time
0
In the last few decades, developmental biologists have extensively used the reaction-diffusion model to explain pattern formation in living organisms. The original model was proposed by Alan Turing in 1952 [1]. The model is based on the idea that pattern formation results from two fundamental mechanisms: (1) coupled catalytic and autocatalytic reactions in a space element between two chemical species, an activator and an inhibitor, and (2) transfer of the interacting species to and from the neighboring space elements through a diffusional transport mechanism. Under appropriate reaction and diffusion conditions, a periodic pattern is formed from an initially homogeneous spatial distribution of activator and inhibitor [2, 3]. Examples of pattern formation can be found in biology, chemistry (the famous Belousov–Zhabotinskii reaction), physics, and mathematics [4, 5].
To illustrate the mechanism of pattern formation, consider the hypothetical activator-inhibitor reaction sequence: ​
(R1)
2A+C
→
3A
(R2)
2A+D
→
2A+B
(R3)
B+C
⇄
E
(R4)
A
→
F
(R5)
B
→
G
The species
D
,
E
,
F
and
G
are supposed to be sufficiently abundant in the reaction mixture that their respective concentrations
[D]
,
[E]
,
[F]
, and
[G]
can be considered constant. The activator
A
reacts with species
C
to produce more
A
by the autocatalytic reaction R1 (so that
A
is reactant, catalyst, and product), but it also promotes the production of the inhibitor
B
by the catalytic reaction R2, in which
A
is a catalyst. Both activator and inhibitor decay with time (reactions R4 and R5).
Because of the equilibrium reaction R3, the concentrations of
B
and
C
are such that
[B]∝1/[C]
. Denoting by
k
j
the reaction constant of reaction
j
, the reaction rates of
A
and
B
are:
r
A
=
2
k
1
[A]
[C]-
k
4
[A]=
-1
2

k
1
[A]
[B]
-
k
4
[A]
and
r
B
=
2
k
2
[A]
-
k
5
[B]
.
From the expression of the reaction rate of
A
,
r
A
, the species
B
inhibits the activator production, hence its name: the larger its concentration, the lower the production rate of
A
.
Introducing the variables
u∝[A]
and
v∝[B]
, the governing equations of the reaction-diffusion system in 2D can be written in the following nondimensional form:
∂u
∂t
=
A
2
∂
u
∂
2
x
+
2
∂
u
∂
2
y
+
2
u
/v-βu
,
∂v
∂t
=
B
2
∂
v
∂
2
x
+
2
∂
v
∂
2
y
+
2
u
-v
,
with
0≤x,y≤1.
Set
β=1
. For the formation of spatial patterns, the diffusion rates of activator and inhibitor should be very different: set the diffusion coefficients to be

A
=0.0005
2
cm
/s
and

B
=0.01 
2
cm
/s
.
For the numerical solution of the system of ODEs, assume (1) periodic boundary conditions as well as (2) the initial conditions:
u(x,y,t=0)=
ψ
uxy
,
v(x,y,t=0)=0.1+
ψ
vxy
,
where
 
ψ
uxy
and
 
ψ
vxy
are numbers taken randomly in the interval
[0,1]
.
The Chebyshev orthogonal collocation method applied with
N
p
=14
nodes in both spatial directions transforms the system of two coupled nonlinear PDEs into a system of 392 nonlinear coupled first-order ordinary differential equations. This system of ODEs is solved using the built-in Mathematica command NDSolve. The snapshots show the results for the 2D inhibitor concentration distribution
v(x,y)
at times
t=0
and at
t=10
. It is interesting to note, in particular, the emergence of Turing patterns similar to leopard spots in the concentration distribution at
t=10
.