MCB111: Mathematics in Biology (Spring 2018)
 A 2D DiffusionReaction process
 The GrayScott model
 Patterns
 Biological Validation
 Turing instability for a twocomponent system
 Epilogue
week 13:
Pattern formation and Turing instability
Alan Turing had the brilliant intuition (brilliant because it was against common thinking, yet extremely simple) that a reaction at a stable steady state solution, could be destabilized by an added diffusion process and ultimately produce periodic spatial patterns.
The original Turing work is “On the chemical basis of morphogenesis”. A good recent source is this “Reactiondiffusion” tutorial by Karl Sims. Another reference is the report “The ReactionDiffusion Equations”.
There is even a ReactionDiffusion media wall at the Boston Museum of Science!
A 2D DiffusionReaction process
Consider a system with species , in two dimensions ,

A diffusion process is characterized by the Fick equations
where (for ) are the diffusions constants with units of .
The Laplace operator (also referred to as ) for a function that depends on and and possibly other variables (like time in this case) stands for

A reaction process is characterized by the deterministic equations
where the functions depend on the species, and a number of rate constants that we designate generically with .
Combining the two sources of variability together for species result in the diffusionreaction equation
The GrayScott model
The GrayScott model is a 2D model of diffusionreaction that can produce interesting patterns. Two species and interact such that 2 molecules of interact with one molecule of in order to transform into a molecule (see Figure 1),
Figure 1. The GrayScott diffusionreaction model. Adapted from K. Sims (http://www.karlsims.com/rd.html).
The molecules are produced and destroyed with the same rate constant (), that is,
While can only be produced by the above reactions involving 2 s and one , decays on its own with a rate () larger than the one controlling the decay of .
The and molecules are also subject to diffusion (in both directions and ) with diffusion constants and respectively. The GrayScott model assumes that molecules of species diffuse faster that those of ().
The GrayScott continuous differential equations
The 2D reaction equation for the GrayScott model are
We cannot solve these equations analytically, but we can do a simulation using discrete time and distance and intervals.
In discrete form, using the Euler algorithm, we can write
What about the Laplacian term?
The Laplacian operator: discrete approximation
The laplacian operator involves the second derivative respect to the space coordinates,
In discrete space, we can write the derivative of a function as
thus and if we assume
Which can be rewritten, introducing a kernel matrix
as
This is called a convolution. In the image analysis lingo, this would be referred to as “approximating the laplacian function by a convolution”.
In image analysis, they usual use convolution kernels that smooth the derivative, which result in image reconstruction with softer edges. A typical kernel, which is the one I am using in this simulation is
You can find functions that perform this convolution for matlab and others, yet is it just a collection of sums, so it is easy to implement from scratch.
Patterns
Now, we are ready to simulate the diffusionreaction equations of the GrayScott model. At a given time , we know the concentrations of both species for all points of the square,
then, we can calculate the concentrations for all points at time ,
using the equations
The parameters I have used in the simulation presented in Figure 2 are
Parameter  Value  Description 

rate of producing species  
birth (and death) rate for species  
is the death rate for species  
diffusion coefficient for species  
diffusion coefficient for species  
total time sampled  
time increment  
total distance sampled in each direction  
space increment (both in x and y directions) 
We assume that we have a lattice with points. At each point and at each time , there is one concentration for each species and .
The initial conditions (Figure 2 at ) correspond to having everywhere in the square (in blue), and everywhere in the square expect for four small (magenta) areas in which . Using these parameters, I obtained the patterns in Figure 2.
Figure 2. Patterns generated with the GrayScott diffusionreaction model as a function of time, using the parameters provided in the text. We depict A molecules in cyan, and B molecules in magenta.
My images (Figure 2) are not as beautiful as Karl Sims’s are, but definitely patterns have emerged! Changing the values of the parameters, you can obtain many different patterns.
There are ways of making the pattern more complicated:
 Make diffusion occur with different constant in different directions.
 Change some the birth/death parameters as a function of the position in the plane.
 Speed up or down the reaction rate.
Biological Validation
Turing’s work was purely theoretical. You may want to check out this work by Fraden & Epstein at Brandais University, “Turing theory of morphogenesis validated”, where they observed experimentally all patterns that Turing predicted, plus seven new ones! here is the actual paper
Turing instability for a twocomponent system
Consider a general reaction with two components and , as in the case of the GrayScott model. Let’s consider which are the conditions under which a stable set of reactions between and , if subject to diffusion would become instable and form patterns. This phenomenon is referred to as a Turing bifurcation or Turing instability.
The diffusionreaction system has the general form
Now consider a steadystate solution and of the reactiononly system
that is, the steadystate solution is characterized by
this solution is stable if the Jacobian matrix defined by
is such that the real part of the eigenvalues is negative.
For a matrix, that is equivalent to the two conditions
 The trace in negative
 And the determinant is positive
Remember that we want to find steady solutions of the reaction that become unstable due to the diffusion process. Instability of the diffusionreaction process for a 2component system requires two additional conditions involving the diffusion coefficients
*
*
Derivation of these results can be found in this tutorial “Notes on turing instability and chemical instabilities” by Michael Cross (page 6), which I reproduce below.
Derivation of conditions for Turing instability for two reactingdiffusing chemicals
Turner’s discovery was that a spatially uniform and stable state (in the absence of diffusion) can become unstable due to non uniformities appearing as a result of diffusion, and patterns start to appear.
When will a diffusionreaction start producing patterns?
For two reactants and , subject to a reactiondiffusion process, their concentrations depending on time and space and satisfy some equations
We assume that we have found some stable solution of the nondiffusion process and which satisfy conditions
Let’s now linearize around those stable solutions
Introducing,
we can write,
Introducing
we have the linearized equations
Introducing vectorial notation
and
we have the linearized equations in vectorial form,
We can make an ansatz that has a constant term , and such that the dependencies in time and space are given by
for some positive constants and .
Notice that we are assuming that both reactants and have the same dependencies in time and space. This is an important simplification for the results that we will find.
Under these assumptions for , the linear differential equation becomes
or
And since this results has to be true for all values of and , we have
This has now become a 2D eigenvalue problem for the 2D matrix
has two solutions, which are the two eigenvectors of corresponding to its two eigenvalues
where the trace of is
and the determinant of is given by
The eigen values are the two solutions to the quadratic equation in resulting from imposing the condition
For a solution to be stable we need two conditions

(a) The solution is real, which requires

(b) The solution is negative, which requires .
Remember that the Turing instability occurs when

(1) There is a stable solution in the absence of diffusion.

(2) The solution becomes instable when diffusion is added.
That is

(1) In the absence of diffusion, we have , such that
and
And stability requires
* (a) $$F^\ast_{AA} + F^\ast_{BB} < 0$$ * (b) $$F^\ast_{AA}\ F^\ast_{BB}  F^\ast_{AB}\ F^\ast_{BA} >0$$ 
(2) When we add diffusion, we want the stability condition to break.
The modified trace is
so is negative as long as is negative.
While , diffusion has to be capable of achieving the condition in order to produce instability.
Notice that
The is a parabola as a function of , and it is positive for .
The instability condition will be achieved if the minimum of the parabola is negative.
The minimum of the parabola corresponds to
and the value of at the minimum is
The being negative at the bottom of the parabola then requires
 or
$$ D_B F^\ast_{AA} + D_A F^\ast_{BB} > 0. $$
or
which can be rewritten as
or
$$ 2 \sqrt{D_A D_B (F^\ast_{AA} F^\ast_{BB}  F^\ast_{AB} F^\ast_{BA})} < D_B F^\ast_{AA} + D_A F^\ast_{BB}. $$
Thus, we have justified the four conditions we stated without proof in the previous section.
Epilogue
Taking a quote from Turing’s paper on morphogenesis, and transferring it to our situation,
“It must be admitted that the biological examples which it has been possible to give in the present [course] are very limited. This can be ascribed quite simply to the fact that biological phenomena are usually very complicated. Taking this in combination with the relatively elementary mathematics used in this [course] one could hardly expect to find that many observed biological phenomena would be covered. [I hope], however, that the imaginary biological systems which have been treated, and the principles which have been discussed, should be of some help in interpreting real biological forms.”
This could directly apply to this course. With the difference that Turing’s “elementary mathematics” are much further advanced than mine.