MCB111: Mathematics in Biology (Spring 2018)


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 “Reaction-diffusion” tutorial by Karl Sims. Another reference is the report “The Reaction-Diffusion Equations”.

There is even a Reaction-Diffusion media wall at the Boston Museum of Science!

A 2D Diffusion-Reaction process

Consider a system with species , in two dimensions ,

Combining the two sources of variability together for species result in the diffusion-reaction equation

The Gray-Scott model

The Gray-Scott model is a 2D model of diffusion-reaction 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 Gray-Scott diffusion-reaction 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 Gray-Scott model assumes that molecules of species diffuse faster that those of ().

The Gray-Scott continuous differential equations

The 2D reaction equation for the Gray-Scott 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 re-written, 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 diffusion-reaction equations of the Gray-Scott 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 Gray-Scott diffusion-reaction 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:

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 two-component system

Consider a general reaction with two components and , as in the case of the Gray-Scott 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 diffusion-reaction system has the general form

Now consider a steady-state solution and of the reaction-only system

that is, the steady-state 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

Remember that we want to find steady solutions of the reaction that become unstable due to the diffusion process. Instability of the diffusion-reaction process for a 2-component 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 reacting-diffusing 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 diffusion-reaction start producing patterns?

For two reactants and , subject to a reaction-diffusion process, their concentrations depending on time and space and satisfy some equations

We assume that we have found some stable solution of the non-diffusion 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

Remember that the Turing instability occurs when

That is

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.