Journal ArticleDOI

Periodic pattern formation in reaction-diffusion systems: an introduction for numerical simulation.

01 Sep 2004-Anatomical Science International (Blackwell Science Pty)-Vol. 79, Iss: 3, pp 112-123
TL;DR: In the present review, a detailed description of the definition of the Turing reaction-diffusion model is provided that is comprehensible without a special mathematical background, then a method for reproducing numerical calculations with Microsoft Excel is illustrated.
Abstract: The aim of the present review is to provide a comprehensive explanation of Turing reaction-diffusion systems in sufficient detail to allow readers to perform numerical calculations themselves. The reaction-diffusion model is widely studied in the field of mathematical biology, serves as a powerful paradigm model for selforganization and is beginning to be applied to actual experimental systems in developmental biology. Despite the increase in current interest, the model is not well understood among experimental biologists, partly because appropriate introductory texts are lacking. In the present review, we provide a detailed description of the definition of the Turing reaction-diffusion model that is comprehensible without a special mathematical background, then illustrate a method for reproducing numerical calculations with Microsoft Excel. We then show some examples of the patterns generated by the model. Finally, we discuss future prospects for the interdisciplinary field of research involving mathematical approaches in developmental biology.

Introduction

• Mathematical modeling, numerical simulation, pattern formation, Turing, also known as Key words.
• As the activator peaks grow, inhibitor peaks should also grow in response because the activator promotes the production of inhibitor.

What is a differential equation?

• These equations are called differential equations (more specifically, partial differential equations) and describe the rate of change of an internal state of a certain physical system in both time and space.
• To fully specify the problem requires three factors: (i) initial conditions; (ii) governing equations; and (iii) boundary conditions.
• The governing equations define the rules on how these values will change in time and space.
• Boundary conditions define how the system behaves at its boundary; for example, the system may be confined within a certain domain, so there would be no flux out of the boundaries.
• The authors will explain what these terms actually mean by using the simplest example.

Initial conditions

• In that case, the authors have 20 pieces of tissue (1/0.05 = 20) that contain both activator and inhibitor, so they define the initial distribution of activator and inhibitor by assigning 20 random numbers to each tissue piece.
• This set of random numbers represents the small noise that should exist in the actual system and, as the authors will see later, any set of random numbers as the initial condition will eventually lead to a similar periodic pattern in this simulation.

Discretization

• Because it is quite difficult to think directly about the time-course of spatial distribution of the molecules, at first the authors divide this rod-like structure into small pieces that have horizontal length dx (dx << 1, meaning dx is much smaller than 1) and suppose the spatial distribution of activator and inhibitor in these small pieces to be homogeneous.
• Then, the authors will think about the concentration change of activator and inhibitor in these small elements.
• There are two factors that affect the concentration of activator and inhibitor in these small pieces: (i) the interaction of activator and inhibitor within each element; and (ii) the transfer of activator and inhibitor between each element and its two nearest neighbors.
• Basically, time should be continuous, but the authors consider updating the system in discrete time steps, dt.

Reaction term

• At first the authors suppose the spatial distributions of activator and inhibitor at time m × dt are known and think about what will happen in each small tissue element after the short time period dt.
• Then, the authors will consider the events that occur solely inside each tissue piece.
• In the activator–inhibitor type reaction–diffusion model (Fig. 3), the activator promotes its own production and promotes the production of inhibitor and the inhibitor inhibits the activator production and decays with time.
• For mathematical simplicity, the authors allow negative values for p and q and set their initial values to 0.
• The dt term arises from observing that the actual increase in chemical concentration is obtained by multiplying the net rate of production by time.

Diffusion term

• Next, the authors will consider the interaction between a tissue element and its two nearest neighboring tissue elements during the time interval (m × dt, (m + 1) × dt).
• In the biological context, there are many ways to transmit signals spatially but, for simplicity, here the authors consider that both activator and inhibitor diffuse passively between tissue elements.
• Therefore, the amount of activator that is transferred from the right neighboring piece (arrow in Fig. 4) is: [4] where dp represents the diffusion coefficient of the activator.
• Similarly, the concentration change of inhibitor is: [10 ] in which dq is the diffusion coefficient of the inhibitor.
• This part, which represents the effect of diffusion, is called the ‘diffusion term’ in the reaction–diffusion system.

Boundary conditions

• As the authors have seen above, the tissues at the boundary should be treated separately.
• By evaluating these cells, the authors can obtain the complete spatial distribution of activator and inhibitor at time dt.

Zero-flux boundary condition

• The authors assume that the boundary is impermeable; that is, no material is transferred across it.
• This condition is also called the Newmann boundary condition.

Transformation to continuous equation

• The authors can obtain the continuous differential equations corresponding to the discrete equations above by making dt and dx infinitely small.
• At first, the authors define the concentration of activator and inhibitor as u(x, t) and v(x, t), which are now continuous functions (x and t are now real values instead of integers).
• This is a ‘partial derivative’ and is defined as a derivative of a function of several variables when all but one variable (the variable of interest) are held fixed during the differentiation.
• The authors describe the continuous equations only for deciphering equation 1, but actually an analytical treatment of the system is possible for the continuous case (see Murray, 2003).
• Numerical calculation of the reaction– diffusion system Parameters and equations used in the numerical calculations.

Governing equation

• First, the authors will calculate the concentration of activator (p) and inhibitor (q) at a certain tissue piece after time dt has passed.
• Let us calculate the activator and inhibitor concentrations for a certain column (C) and write the values to C4 and C5.
• The authors can obtain values in other cells by applying the same equation (Fig. 6).
• Fortunately, Excel automatically converts these equations to appropriate forms by simply copying and pasting cells C4 and C5 to other cells.

Numerical calculation

• Fortunately Excel again automatically converts the equations to appropriate forms by simply copying and pasting the whole row to the rows below.
• By doing this 100–200 times repeatedly, you can observe that the concentrations gradually form a periodic structure, as shown in Fig.
• Numerical calculation with Mathematica Obviously, it is too labor consuming to undertake the above procedure, so the authors generally use Mathematica (Wolfram Research, Champaign, IL, USA) to calculate the results.
• The details of the program can be provided electronically.
• (All the calculations in the present paper are performed by Mathematica and the source code (with additional instructions on linear stability analysis) is freely available on request from the authors.).

Properties of Turing reaction–diffusion systems

• Relationship between domain size and number of structures.
• One characteristic of reaction–diffusion systems is that they tend to form structures of similar size , so if the domain size is changed the number of structures should change, not the size of each structure.

Changing initial conditions

• One characteristic of the Turing reaction–diffusion system is that it has the ability to form de novo stable periodic patterns; that is, without any prepattern.
• As you can see in Fig. 10, if the initial condition is not homogeneous, pattern formation can occur sequentially.
• The final periodic structure is more or less the same.
• Detailed analysis on the pattern appearance speed has been recently undertaken in Miura and Maini (2004).

Changing diffusion coefficients

• As you can see from the above discussion, the reaction–diffusion system has the ability to make a periodic pattern from an almost homogeneous initial state and the wavelength of the pattern is decided by the parameters in the reaction and diffusion terms.
• This is experimentally assayed by Miura and Shiota (2000a) in a limb bud mesenchyme cell culture system.
• Actually, the result is easy to understand without numerical calculation.
• The result is confirmed by numerical calculation (Fig. 12).

Cross-type reaction–diffusion model

• The activator–inhibitor scheme is well known among developmental biologists, but actually there is another type of reaction–diffusion model that has the ability to generate periodic pattern.
• This system, usually called the substrate-depletion system, also consists of two hypothetical molecules, the substrate and the enzyme.
• Numerical simulations are shown in Fig. 13.
• In the activator–inhibitor system, the inhibitor should be expressed at chondrogenic sites and have an ability to inhibit chondrogenesis.

Reaggregated tissue experiment

• As the authors saw in the previous subsection, one important characteristic of the Turing reaction–diffusion system is an ability to form a periodic pattern from various initial states and the process is quite robust.
• Here, the authors try to emulate the most extreme case, where the tissue is dissociated into single cells and reaggregated.
• In contrast, if the number of structures is quite stable, the authors have to couple the reaction–diffusion system with other mechanisms to explain the stability.
• The reliability of the pattern formation mechanism was first discussed by Bard and Lauder (1974) and a possible scenario for robust pattern formation is proposed by Crampin et al. (1999).
• The authors introduce the latter in the following subsection.

Growing domain

• In actual biological systems (especially embryonic tissue), the size and shape of the pattern formation field is usually not constant, but grows.
• The Turing reaction–diffusion system has the property that it retains the periodic structure of fixed wavelength, so if the tissue grows the authors can expect that the pattern will change according to the growth.
• Next, the authors slightly modify the reaction term of the system (Fig. 16b).
• In the previous simulation, additional peaks are inserted between pre-existing peaks, but this time each peak is split into two peaks.
• This distinction of peak increase (called ‘mode doubling’) is investigated analytically by Crampin et al. (2002).

Two dimensions: stripe–spot selection

• The authors can easily modify the above program for simulation on a two-dimensional spatial domain (a flat surface).
• One finds stripes, spots or more complex patterns.
• One can see, in Fig. 17, a stripe- and spot-like pattern.the authors.
• Stripe–spot selection has been studied in a special case (Ermentrout, 1991; Lyons & Harrison, 1992), but the mechanism, in general, is not fully understood.
• Numerical calculations on two-dimensional curved surfaces (Varea et al., 1999) and three dimensions (Leppänen et al., 2002) have been done recently and some interesting features are observed concerning the connectivity of the periodic pattern.

Future prospects

• There are several groups of researchers who are trying to understand biological pattern formation using mathematical models.
• The field is expanding rapidly and a considerable number of researchers (although much less than in developmental biology) is involved.
• As far as the authors know, there are only two to three groups in the world that can deal with this kind of problem both theoretically and experimentally.
• It seems that the two cultural differences described above can be an energy barrier for experimental people to enter this field.

Acknowledgments

• The authors thank Professor Stuart Newman (New York Medical College) and Dr Chad Perlin (Department of Human Anatomy and Genetics, University of Oxford) for their critical reading of the manuscript.
• This work was supported by the Japan Society for the Promotion of Science.

Did you find this useful? Give us your feedback

Content maybe subject to copyright    Report

Anatomical Science International (2004) 79, 112–123
Blackwell Publishing, Ltd.
Review Article
Periodic pattern formation in reactiondiffusion systems:
An introduction for numerical simulation
Takashi Miura* and Philip K. Maini
* Department of Human Anatomy and Genetics and
Centre for Mathematical Biology, Mathematical Institute, University of Oxford, UK
Abstract
The aim of the present review is to provide a comprehensive explanation of Turing reaction–diffusion systems
in sufﬁcient detail to allow readers to perform numerical calculations themselves. The reaction–diffusion
model is widely studied in the ﬁeld of mathematical biology, serves as a powerful paradigm model for self-
organization and is beginning to be applied to actual experimental systems in developmental biology. Despite
the increase in current interest, the model is not well understood among experimental biologists, partly
because appropriate introductory texts are lacking. In the present review, we provide a detailed description
of the deﬁnition of the Turing reaction–diffusion model that is comprehensible without a special mathemat-
ical background, then illustrate a method for reproducing numerical calculations with Microsoft Excel. We
then show some examples of the patterns generated by the model. Finally, we discuss future prospects for
the interdisciplinary ﬁeld of research involving mathematical approaches in developmental biology.
Key words:
mathematical modeling, numerical simulation, pattern formation, Turing.
Introduction: Periodic pattern formation
in biological systems and the Turing
reactiondiffusion model
It has been reported that, in many cases during
development, periodic patterns emerge. Examples
include the skin pigment pattern in zebra (Bard,
1981; Murray, 2003), angelﬁsh (Kondo & Asai, 1995;
Shoji et al., 2003), zebraﬁsh (Asai et al., 1999) and
sea shells (Meinhardt, 1995), feather follicle forma-
tion (Jung et al., 1998), tooth development (Salazar-
Ciudad & Jernvall, 2002) and digit formation during
limb development (Newman & Frisch, 1979; Maini &
Solursh, 1991; Dowine and Newman, 1994, 1995;
Miura and Shiota, 2000a, 2000b; Miura et al., 2000;
Moftah et al., 2002). Certain aspects of these pattern
formation processes cannot be easily explained
simply by the combination of morphogen gradients.
Recently, some developmental biologists have
started using the Turing reaction–diffusion model for
this type of pattern formation. This model was origin-
ally proposed by the British mathematician Alan
Turing (Turing, 1952) and a signiﬁcant amount of
work has been done using this idea in the ﬁeld of
mathematical biology (for reviews, see Bard, 1990;
Meinhardt, 1995; Murray, 2003). This model hypo-
thesizes the existence of two molecules, an activator
and an inhibitor, and, if they interact with each other
in a speciﬁc manner (see below), a periodic pattern
is formed from a homogeneous initial spatial distribu-
tion of activator and inhibitor. The qualitative expla-
nation of pattern formation is as follows: because
the activator has an ability to enhance its own pro-
duction, any small peak of activator in the initial
distribution is ampliﬁed. As the activator peaks grow,
inhibitor peaks should also grow in response because
the activator promotes the production of inhibitor.
Inhibitor peaks should be less steep than activator
peaks owing to the assumption that the inhibitor has a
larger diffusion coefﬁcient, which results in the inhibi-
tion of new activator peak formation near pre-existing
peaks. This results in a periodic pattern of activator
and inhibitor peaks (Meinhardt, 1995; Kondo, 2002).
This explanation itself has become well known to
developmental biologists and some standard devel-
opmental biology textbooks have started to take up
this topic (see Wolpert, 1998; Gilbert, 2000). Most
experimental biologists can follow the description
of the logic, but very few actually use the model
because it is necessary to at least reproduce the
numerical simulation results to do any serious scien-
tiﬁc research using the model. One obstacle for this
is that there is virtually no introductory review that is
suitable for this purpose. Therefore, the aim of this
Correspondence: Takashi Miura, Department of Anatomy and
Developmental Biology, Kyoto University Graduate School of
Medicine, Yoshida Konoe-Chou, Sakyo-ku, 606-8501, Japan.
Email: miura-takashi@umin.ac.jp
This paper is based on the winning article of the 2000 Incentive
Prize of the Japanese Association of Anatomists.
Received 8 March 2004; accepted 14 April 2004.

Tu r ing reactiondiffusion systems 113
text is to provide enough information for the reader
to perform actual numerical calculations on reaction–
diffusion models. We will provide an actual program
written by Microsoft Excel (Redmond, WA, USA) and
will provide a Mathematica program on request. The
(2003; Volume II, Chapter 2), which contains the
basic mathematical explanation of why periodic pat-
terns are formed in the model.
Model equations
What is a differential equation?
Most of the aforementioned texts that deal with
reaction–diffusion systems are aimed at readers with
a mathematical background, so they immediately
plunge into the mantra-like string:
u = f(u, v) + d
u
u
v = g(u, v) + d
v
v [1]
which is very puzzling for most experimental biolog-
ists at ﬁrst glance. These equations are called differ-
ential equations (more speciﬁcally, partial differential
equations) and describe the rate of change of an
internal state of a certain physical system in both
time and space. To fully specify the problem requires
three factors: (i) initial conditions; (ii) governing equa-
tions; and (iii) boundary conditions.
Initial conditions mean the values of the given sys-
tem at the very beginning. The governing equations
deﬁne the rules on how these values will change in
time and space. Boundary conditions deﬁne how the
system behaves at its boundary; for example, the
system may be conﬁned within a certain domain, so
there would be no ﬂux out of the boundaries. We will
explain what these terms actually mean by using
the simplest example. To facilitate understanding, the
deﬁnitions of each variable are given in the Appendix.
Initial conditions
For simplicity, we will think about a quasi one-dimen-
sional rod-like embryonic tissue (Fig. 1). We deﬁne
horizontal length as 1 and vertical length as dy (<< 1,
much smaller than 1) and suppose the distribution
of molecules in the vertical direction is negligible.
Cells inside this tissue produce two diffusible sign-
aling molecules, the activator and the inhibitor, and
these molecules control the production (or degrada-
tion) of both molecules and diffuse to neighboring
cells. The distribution can be homogeneous or can
have some prepattern, as shown in later sections,
depending on the actual experimental situation.
Discretization
Because it is quite difﬁcult to think directly about the
time-course of spatial distribution of the molecules,
at ﬁrst we divide this rod-like structure into small
pieces that have horizontal length dx (dx << 1, mean-
ing dx is much smaller than 1) and suppose the
spatial distribution of activator and inhibitor in these
small pieces to be homogeneous. Then, we will
think about the concentration change of activator
and inhibitor in these small elements. There are two
factors that affect the concentration of activator and
inhibitor in these small pieces: (i) the interaction of
activator and inhibitor within each element; and (ii)
the transfer of activator and inhibitor between each
element and its two nearest neighbors.
We consider the sum of these two factors as the
concentration change of these molecules within a
speciﬁc element in a short period of time.
Basically, time should be continuous, but we con-
sider updating the system in discrete time steps, dt.
So, at time step m (where m is a positive integer),
a time of m × dt has actually passed (see Fig. 2). We
can deﬁne the concentration of activator molecule in
the n th tissue element from the left-hand boundary
at time m × dt as p (n, m) and the concentration of
inhibitor as q(n, m).
Reaction term
At ﬁrst we suppose the spatial distributions of activator
and inhibitor at time m × dt are known and think about
what will happen in each small tissue element after
the short time period dt. Then, we will consider the
events that occur solely inside each tissue piece. There
are two types of Turing reaction–diffusion model that
are known to generate spatial patterns: (i) the activator
inhibitor type; and (ii) the substrate–depletion type.
In the activator–inhibitor type reaction–diffusion model
(Fig. 3), the activator promotes its own production and
promotes the production of inhibitor and the inhibitor
inhibits the activator production and decays with time.
So, for example, if we set the rate of change of
activator and inhibitor as f(p, q) and g(p, q), respect-
ively, and consider:
f (p, q) = 0.6p q
g (p, q) = 1.5p 2q [2]
Figure 1. Deﬁnition of initial conditions and discretization in
space.

114 T. Miura and P.K. Maini
we then have a system in which p is the activator
and q is the inhibitor. These are called the ‘reaction
terms’.
Actual biologically relevant functions will have
saturation limits for the concentration of activator
and inhibitor, so the reaction terms will have more
complex form. For mathematical simplicity, we allow
negative values for p and q and set their initial values
to 0. If we set actual concentrations of activator and
inhibitor as P and Q, they will then be transformed
to p and q, with p = P P
0
and q = Q Q
0
, where P
0
and Q
0
are certain positive constants.
We can calculate the amount of change of activ-
ator and inhibitor during the time interval spanning
(m × dt, (m + 1) × dt) as:
f (p(m, n), q (m, n)) × dt
g(p (m, n), q(m, n)) × dt [3]
respectively. The dt term arises from observing that the
actual increase in chemical concentration is obtained
by multiplying the net rate of production by time.
Diffusion term
Next, we will consider the interaction between a
tissue element and its two nearest neighboring tissue
elements during the time interval (m × dt, (m + 1) × dt).
In the biological context, there are many ways to
transmit signals spatially but, for simplicity, here we
consider that both activator and inhibitor diffuse
passively between tissue elements.
If we suppose the concentration of activator in the
nth tissue element is p(n, m), then the concentration
of activator in the tissue element to the right will be
p(n + 1, m). In the case of simple diffusion, the amount
of molecule transferred from element n to n + 1 is
proportional to the concentration gradient (p(n + 1, m)
p(n, m))/dx and transverse length of the element dy
(Fick’s Law). Therefore, the amount of activator that
is transferred from the right neighboring piece (arrow
in Fig. 4) is:
[4]
where d
p
represents the diffusion coefﬁcient of the
activator. The concentration change induced by this
transfer can be obtained by dividing this value by
the area of the tissue element (dx × dy).
[5]
[6]
Figure 2. Discretization in space and time.
Figure 3. Schematic representation of the reaction term.
Figure 4. Transfer of activator or inhibitor between neighboring
tissue pieces.
d
pn m pnm
dx
dy dt
p
×
+−
××
(( , ) (, ))
1
d
pn m pnm
dx
dy dt dx dy
p
×
+−
×× ×
(( , ) (, ))
/( )
1
=
+−
×
(( , ) (, ))
dpn m pnm
dx
dt
p
1
2

Tu r ing reactiondiffusion systems 115
Similarly, we can obtain the amount of activator trans-
ferred from the left neighboring element as:
[7]
and the concentration change induced by the trans-
fer as:
[8]
Taken together, the change of activator concentration
in element n between time m × dt and (m + 1) × dt is:
[9]
Similarly, the concentration change of inhibitor is:
[10
]
in which d
q
is the diffusion coefﬁcient of the inhibitor.
This part, which represents the effect of diffusion,
is called the ‘diffusion term’ in the reactiondiffusion
system.
Governing equations
Taking both reaction and diffusion into consideration,
the changes of concentrations of activator and inhib-
itor during the time interval (m × dt, (m + 1) × dt) are:
p(n, m + 1) p(n, m) = (f(p(n, m), q(n, m)) +
d
p
(p(n + 1, m) + p(n 1, m) 2 × p(n, m))/dx
2
) × dt
q(n, m + 1) q(n, m) = (g(p(n, m), q(n, m)) +
d
q
(q(n + 1, m) + q (n 1, m) 2 × q(n, m))/dx
2
) × dt
[11]
From these, we can calculate the concentrations of
activator and inhibitor at time (m + 1) × dt as:
p(n, m + 1) = p(n, m) + (f(p(n, m), q(n, m)) +
d
p
(p(n + 1, m) + p(n 1, m) 2 × p(n, m))/dx
2
) × dt
q(n, m + 1) = q(n, m) + (g(p(n, m), q(n, m)) +
d
q
(q(n + 1, m) + q(n 1, m) 2 × q(n, m))/dx
2
) × dt
[12]
Therefore, if we know the initial spatial distribution of
activator and inhibitor, we can obtain the concentra-
tion distribution of activator and inhibitor at arbitrary
time by repeatedly applying this equation at each time
step dt.
Boundary conditions
Not all tissue pieces obey the rule described above.
In the left-most and right-most tissue we have to deﬁne
special conditions. There are several ways to do this,
depending on the biological situation. If the total
number of tissue elements is N
total
, then p(1, m) and
p(N
total
, m) have only one neighbor and we cannot
treat them as above; that is, equations 12 only hold
for integer n with n greater than 1 and less than
N
total
.
This condition poses a constraint on the number
of waves in a certain region. Speciﬁcally, the number
of waves has to be an integer in the deﬁned region.
The effect of the boundary becomes stronger when
there are only a small number of structures in the
deﬁned region.
Periodic boundary condition
We deﬁne the values of p(1, m) and p(N
total
, m) to
be equal for all m, so that the left-most tissue and
right-most tissue are connected.
Zero-ﬂux boundary condition
We assume that the boundary is impermeable; that
is, no material is transferred across it. In this case, we
have to deﬁne the change of p(1, m) and p(N
total
, m)
as follows:
p(1, m + 1) = p(1, m) + (f(p(1, m), q(1, m))
+ d
p
(p(2, m) p(1, m))/dx
2
) × dt [13]
p(N
total
, m + 1) = p (1, m) + (f(p(N
total
, m), q(N
total
, m))
+ d
p
(p(N
total
1, m) p (N
total
, m))/dx
2
) × dt
[14]
This condition is also called the Newmann boundary
condition.
Fixed boundary condition
In this case, we have to set p(1, m) and p(N
total
, m)
to speciﬁc values for all time. Hence:
p(1, m + 1) = α
p(N
total
, m + 1) = β [15]
for each iteration of the simulation where α and β
are ﬁxed (non-negative) numbers. This is also called
a Dirichlet condition.
Transformation to continuous equation
We can obtain the continuous differential equations
corresponding to the discrete equations above by
making dt and dx inﬁnitely small. At ﬁrst, we deﬁne
the concentration of activator and inhibitor as u(x, t)
and v(x, t), which are now continuous functions (x
and t are now real values instead of integers). Then,
the above discrete governing equations become:
[16]
[17]
d
pn m pnm
dx
dy dt
p
×
−−
××
(( , ) (, ))
1
d
pn m pnm
dx
dt
p
×
−−
×
(( , ) (, ))
1
2
d
pn m pn m pnm
dx
dt
p
×
++−−×
×
(( , ) ( , ) ( , ))
112
2
d
qn m qn m qn m
dx
dt
q
×
++ ×
×
(( , ) ( , ) ( , ))
112
2
uxt dt uxt
dt
fuxt vxt
d
dx
p
ux dxt uxt
dx
uxt ux dxt
dx
(, ) (, )
((, ), (, ))
(, )(, ) (, ) ( , )
+−
=+
+−
vxt dt vxt
dt
gux t vx t
d
dx
q
v
xdxt
v
xt
dx
v
xt
vx
dx t
dx
(, ) (, )
((, ), (, ))
(, )(, ) (, ) ( , )
+−
=+
+
−−

116 T. Miura and P.K. Maini
If dt and dx tend to zero, the left-hand side becomes a
ﬁrst-order time differentiation and the right-hand side
becomes a second-order spatial differentiation as
follows:
u(x, t )/
t = f(u(x, t ), v (x, t )) + d
p
2
u(x, t)/
x
2
v(x, t )/
t = g(u(x, t ), v (x, t )) + d
q
2
v(x, t)/
x
2
[18]
These are the governing equations of a reaction–
diffusion system.
Some readers may be puzzled to see the strange
character ‘
’. This is a ‘partial derivative’ and is
deﬁned as a derivative of a function of several
variables when all but one variable (the variable of
interest) are held ﬁxed during the differentiation. For
example, the left-hand side of the equation deals
with a rate of change with time of concentration at
a ﬁxed location in space, which means differentiating
with t while keeping x ﬁxed.
Moreover, if you remember the following conventions,
you can decipher equations 1: (i) ﬁrst-order time dif-
ferentiation
u/
t is sometimes written as u; (ii) ﬁrst-
order space differentiation
u/
x is sometimes written
as u (‘’ is called ‘del’ or ‘Nabla’); (iii) second-order
space differentiation
2
u/
x
2
is sometimes written as
u (‘’ is called Laplacian).
Equations 18 are based on the assumption that
the domain is one-dimensional. The equations can
be generalized to fully three-dimensional space,
where they take the following form:
[19]
The symbol
2
u is short-hand for ((
2
/
x
2
)u) +
((
2
/
y
2
)u) + ((
2
/
z
2
)u), sometimes also written as
. In this case, u = u(x, y, z, t) and v = v (x, y, z, t),
where (x, y, z) is position in three-dimensional space.
We describe the continuous equations only for
deciphering equation 1, but actually an analytical
treatment of the system is possible for the continuous
case (see Murray, 2003). This is beyond the scope
of the present paper so, from now on, we concen-
trate on the numerical calculation of the discrete
version of the model.
Numerical calculation of the reaction–
diffusion system
Parameters and equations used in the
numerical calculations
The reactiondiffusion system we now focus on is
as follows:
[20]
on spatial domain size [0, 1] with zero-ﬂux boundary
conditions.
Numerical calculation in Microsoft Excel
To construct the numerical calculation program, we
use Microsoft Excel, which is one of the most widely
available programs for biologists. As we will see
below, it is too time consuming to undertake all the
simulations with Excel, but it is worth trying at least
once by yourself to get a better understanding of the
numerical calculation and an appreciation of the pit-
falls involved.
Initial conditions
Suppose we are using the system in equations 20
with domain size 1 and we discretize it with dx = 0.05
and dt = 0.1. In that case, we have 20 pieces of tis-
sue (1/0.05 = 20) that contain both activator and
inhibitor, so we deﬁne the initial distribution of activ-
ator and inhibitor by assigning 20 random numbers
to each tissue piece. To describe this in an Excel
spreadsheet, we assume each column represents
the concentration of activator or inhibitor in a speciﬁc
tissue piece. So, the initial concentration of the activ-
ator and inhibitor can be expressed as a 20 × 2
matrix of numbers in the spreadsheet, as shown in
Fig. 5. In the case below, we use cell B1 U2 to
deﬁne the initial distribution of activator and inhibitor.
This set of random numbers represents the small
noise that should exist in the actual system and, as
we will see later, any set of random numbers as the
initial condition will eventually lead to a similar peri-
odic pattern in this simulation.
u
t
fuv d u
v
t
guv d v
p
q
(, )
(, )
=+
=+
2
2
u
t
uvu
u
x
v
t
uv
v
x
. .
. .
=−+
=−+
06 00002
15 2 001
3
2
2
2
2
Figure 5. Initial distribution of p and q in the Excel spreadsheet.

Citations
More filters
Journal ArticleDOI
24 Sep 2010-Science
TL;DR: The Turing or reaction-diffusion (RD) model is one of the best-known theoretical models used to explain self-regulated pattern formation in the developing animal embryo as mentioned in this paper.
Abstract: The Turing, or reaction-diffusion (RD), model is one of the best-known theoretical models used to explain self-regulated pattern formation in the developing animal embryo. Although its real-world relevance was long debated, a number of compelling examples have gradually alleviated much of the skepticism surrounding the model. The RD model can generate a wide variety of spatial patterns, and mathematical studies have revealed the kinds of interactions required for each, giving this model the potential for application as an experimental working hypothesis in a wide variety of morphological phenomena. In this review, we describe the essence of this theory for experimental biologists unfamiliar with the model, using examples from experimental studies in which the RD model is effectively incorporated.

1,309 citations

01 Jan 2010
TL;DR: The essence of this theory for experimental biologists unfamiliar with the response-diffusion model is described, using examples from experimental studies in which the RD model is effectively incorporated.
Abstract: Turing Model Explained The reaction-diffusion (Turing) model is a theoretical model used to explain self-regulated pattern formation in biology. Although many biologists have heard of this model, a better understanding of the concept would aid its application to many research projects and developmental principles. Kondo and Miura (p. 1616) now review the reaction-diffusion model. Despite the associated mathematics, the basic idea of the Turing model is relatively easy to understand and relates to morphogen gradients. In addition, user-friendly software makes it easy to understand how a whole variety of patterns can be produced by this simple mechanism. The Turing, or reaction-diffusion (RD), model is one of the best-known theoretical models used to explain self-regulated pattern formation in the developing animal embryo. Although its real-world relevance was long debated, a number of compelling examples have gradually alleviated much of the skepticism surrounding the model. The RD model can generate a wide variety of spatial patterns, and mathematical studies have revealed the kinds of interactions required for each, giving this model the potential for application as an experimental working hypothesis in a wide variety of morphological phenomena. In this review, we describe the essence of this theory for experimental biologists unfamiliar with the model, using examples from experimental studies in which the RD model is effectively incorporated.

204 citations

Journal ArticleDOI
TL;DR: It is shown that the expected morphologies that would arise during this relatively unconstrained "physical" stage of evolution correspond to the hollow, multilayered and segmented morphotypes seen in the gastrulation stage embryos of modern-day metazoa as well as in Ediacaran fossil deposits of approximately 600 Ma.
Abstract: By examining the formative role of physical processes in modern-day developmental systems, we infer that although such determinants are subject to constraints and rarely act in a "pure" fashion, they are identical to processes generic to all viscoelastic, chemically excitable media, non-living as well as living. The processes considered are free diffusion, immiscible liquid behavior, oscillation and multistability of chemical state, reaction-diffusion coupling and mechanochemical responsivity. We suggest that such processes had freer reign at early stages in the history of multicellular life, when less evolution had occurred of genetic mechanisms for stabilization and entrenchment of functionally successful morphologies. From this we devise a hypothetical scenario for pattern formation and morphogenesis in the earliest metazoa. We show that the expected morphologies that would arise during this relatively unconstrained "physical" stage of evolution correspond to the hollow, multilayered and segmented morphotypes seen in the gastrulation stage embryos of modern-day metazoa as well as in Ediacaran fossil deposits of ~600 Ma. We suggest several ways in which organisms that were originally formed by predominantly physical mechanisms could have evolved genetic mechanisms to perpetuate their morphologies.

183 citations

Journal ArticleDOI
TL;DR: Evidence is reported for the importance of the abiotic terrain feature ‘gradient’ in guiding the movements of African savannah elephants and global Positioning System tracking data overlaid onto digital elevation and surface gradient models show that elephants tend to avoid steep slopes.

173 citations

Journal ArticleDOI

TL;DR: The foundation of a unified, object-oriented, three-dimensional biomodelling environment, which allows us to integrate multiple submodels at scales from subcellular to those of tissues and organs, is presented.
Abstract: In this paper we present the foundation of a unified, object-oriented, three-dimensional biomodelling environment, which allows us to integrate multiple submodels at scales from subcellular to those of tissues and organs. Our current implementation combines a modified discrete model from statistical mechanics, the Cellular Potts Model, with a continuum reaction–diffusion model and a state automaton with well-defined conditions for cell differentiation transitions to model genetic regulation. This environment allows us to rapidly and compactly create computational models of a class of complex-developmental phenomena. To illustrate model development, we simulate a simplified version of the formation of the skeletal pattern in a growing embryonic vertebrate limb.

137 citations

Cites background from "Periodic pattern formation in react..."

• ...In the isotropic case, the boundary conditions in equation (4.7a) simplify to: vui vn̂ Z 0; (4.7b) The initial conditions are: uðx; 0ÞZ u initðxÞ: (4.8) For biological applications of RD see, among others, Meinhardt (1982), Murray (1993) and Miura & Maini (2004a)....

[...]

• ...Turing (1952) introduced the idea that interactions of reacting and diffusing chemicals (usually of two species) could form self-organizing instabilities that provide the basis for biological spatial patterning (e.g. animal coat patterning; see Murray 1993; Miura & Maini 2004a for reviews)....

[...]

• ...Several models attempt to account for RD of morphogenetic signalling during chondrogenesis in the limb (Newman & Frisch 1979; Hentschel et al. 2004) and its isolated mesenchymal tissue (Miura & Shiota 2000a,b; Miura et al. 2000; Miura & Maini 2004b)....

[...]

References
More filters
Journal ArticleDOI
TL;DR: In this article, it is suggested that a system of chemical substances, called morphogens, reacting together and diffusing through a tissue, is adequate to account for the main phenomena of morphogenesis.
Abstract: It is suggested that a system of chemical substances, called morphogens, reacting together and diffusing through a tissue, is adequate to account for the main phenomena of morphogenesis. Such a system, although it may originally be quite homogeneous, may later develop a pattern or structure due to an instability of the homogeneous equilibrium, which is triggered off by random disturbances. Such reaction-diffusion systems are considered in some detail in the case of an isolated ring of cells, a mathematically convenient, though biologically unusual system. The investigation is chiefly concerned with the onset of instability. It is found that there are six essentially different forms which this may take. In the most interesting form stationary waves appear on the ring. It is suggested that this might account, for instance, for the tentacle patterns on Hydra and for whorled leaves. A system of reactions and diffusion on a sphere is also considered. Such a system appears to account for gastrulation. Another reaction system in two dimensions gives rise to patterns reminiscent of dappling. It is also suggested that stationary waves in two dimensions could account for the phenomena of phyllotaxis. The purpose of this paper is to discuss a possible mechanism by which the genes of a zygote may determine the anatomical structure of the resulting organism. The theory does not make any new hypotheses; it merely suggests that certain well-known physical laws are sufficient to account for many of the facts. The full understanding of the paper requires a good knowledge of mathematics, some biology, and some elementary chemistry. Since readers cannot be expected to be experts in all of these subjects, a number of elementary facts are explained, which can be found in text-books, but whose omission would make the paper difficult reading.

9,015 citations

"Periodic pattern formation in react..." refers methods in this paper

• ...This model was originally proposed by the British mathematician Alan Turing (Turing, 1952) and a significant amount of work has been done using this idea in the field of mathematical biology (for reviews, see Bard, 1990; Meinhardt, 1995; Murray, 2003)....

[...]

Book
01 Jan 1997
TL;DR: In this paper, the Drosophila body plan was described and the Xenopus and zebrafish were shown to complete the body plan in the early stages of the development of the human embryo.
Abstract: 1. History and basic concepts 2. Development of the Drosophila body plan 3. Vertebrate development I: Life cycles and experimental techniques 4. Vertebrate development II: Xenopus and zebrafish 5. Vertebrate development III: Chick and mouse-completing the body plan 6. Development of nematodes and sea urchins 7. Plant development 8. Cell differentiation and stem cells 9. Morphogenesis: change in form in the early embryo 10. Germ cells, fertilization, and sex 11. Organogenesis 12. Development of the nervous system 13. Growth, post-embryonic development and regeneration 14. Evolution and development

1,138 citations

Journal ArticleDOI
31 Aug 1995-Nature
TL;DR: The striking similarity between the actual and simulated pattern rearrangement strongly suggests that a reactiondiffusion wave is a viable mechanism for the stripe pattern of Pomacanthus.
Abstract: IN 1952, Turing proposed a hypothetical molecular mechanism, called the reaction–diffusion system1, which can develop periodic patterns from an initially homogeneous state. Many theoretical models based on reaction–diffusion have been proposed to account for patterning phenomena in morphogenesis2–4, but, as yet, there is no conclusive experimental evidence for the existence of such a system in the field of biology5–8. The marine angelfish, Pomacanthus has stripe patterns which are not fixed in their skin. Unlike mammal skin patterns, which simply enlarge proportionally during their body growth, the stripes of Pomacanthus maintain the spaces between the lines by the continuous rearrangement of the patterns. Although the pattern alteration varies depending on the conformation of the stripes, a simulation program based on a Turing system can correctly predict future patterns. The striking similarity between the actual and simulated pattern rearrangement strongly suggests that a reaction–diffusion wave is a viable mechanism for the stripe pattern of Pomacanthus.

781 citations

"Periodic pattern formation in react..." refers background in this paper

• ...Examples include the skin pigment pattern in zebra (Bard, 1981; Murray, 2003), angelfish (Kondo & Asai, 1995; Shoji et al., 2003), zebrafish (Asai et al., 1999) and sea shells (Meinhardt, 1995), feather follicle formation (Jung et al., 1998), tooth development (SalazarCiudad & Jernvall, 2002) and…...

[...]

Book
01 Jan 1995
TL;DR: The patterns on the shells of tropical sea snails are not only compellingly beautiful but also tell a tale of biological development, which follows laws like those of dune formation or the spread of a flu epidemic.
Abstract: The patterns on the shells of tropical sea snails are not only compellingly beautiful but also tell a tale of biological development. The decorative patterns are records of their own genesis, which follows laws like those of dune formation or the spread of a flu epidemic. Hans Meinhardt has analyzed the dynamical processes that form these patterns and retraced them faithfully in computer simulations. His book is exciting not only for the astonishing scientific knowledge it reveals but also for its fascinating pictures. An accompanying CD-ROM with the corresponding algorithms offers wide scope to those who wish to try their hand at simulating and varying the patterns.

480 citations

"Periodic pattern formation in react..." refers background or methods in this paper

• ...This model was originally proposed by the British mathematician Alan Turing (Turing, 1952) and a significant amount of work has been done using this idea in the field of mathematical biology (for reviews, see Bard, 1990; Meinhardt, 1995; Murray, 2003)....

[...]

• ...…pattern in zebra (Bard, 1981; Murray, 2003), angelfish (Kondo & Asai, 1995; Shoji et al., 2003), zebrafish (Asai et al., 1999) and sea shells (Meinhardt, 1995), feather follicle formation (Jung et al., 1998), tooth development (SalazarCiudad & Jernvall, 2002) and digit formation during limb…...

[...]

• ...This results in a periodic pattern of activator and inhibitor peaks (Meinhardt, 1995; Kondo, 2002)....

[...]

Journal ArticleDOI
TL;DR: Large morphological effects frequently can be achieved by small changes, according to this model, and similar morphologies can be produced by different changes, consistent with why predicting the morphological outcomes of molecular experiments is challenging.
Abstract: Generation of morphological diversity remains a challenge for evolutionary biologists because it is unclear how an ultimately finite number of genes involved in initial pattern formation integrates with morphogenesis. Ideally, models used to search for the simplest developmental principles on how genes produce form should account for both developmental process and evolutionary change. Here we present a model reproducing the morphology of mammalian teeth by integrating experimental data on gene interactions and growth into a morphodynamic mechanism in which developing morphology has a causal role in patterning. The model predicts the course of tooth-shape development in different mammalian species and also reproduces key transitions in evolution. Furthermore, we reproduce the known expression patterns of several genes involved in tooth development and their dynamics over developmental time. Large morphological effects frequently can be achieved by small changes, according to this model, and similar morphologies can be produced by different changes. This finding may be consistent with why predicting the morphological outcomes of molecular experiments is challenging. Nevertheless, models incorporating morphology and gene activity show promise for linking genotypes to phenotypes.

384 citations