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 $X_1\ldots X_N$, in two dimensions $(x,y)$,

• A diffusion process is characterized by the Fick equations

where $D_{X_i}$ (for $1\leq i \leq N$) are the diffusions constants with units of $[distance]^2/[time]$.

The Laplace operator $\nabla^2$ (also referred to as $\Delta$) for a function $X(x,y,t)$ that depends on $x$ and $y$ and possibly other variables (like time $t$ in this case) stands for

• A reaction process is characterized by the deterministic equations

where the functions $F_i(X_1,\ldots, X_N; k)$ depend on the species, and a number of rate constants that we designate generically with $k$.

Combining the two sources of variability together for species $X_i$ 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 $A$ and $B$ interact such that 2 molecules of $B$ interact with one molecule of $A$ in order to transform $A$ into a $B$ molecule (see Figure 1),

Figure 1. The Gray-Scott diffusion-reaction model. Adapted from K. Sims (http://www.karlsims.com/rd.html).

The $A$ molecules are produced and destroyed with the same rate constant ($f$), that is,

While $B$ can only be produced by the above reactions involving 2 $B$s and one $A$, $B$ decays on its own with a rate ($f+g$) larger than the one controlling the decay of $A$.

The $A$ and $B$ molecules are also subject to diffusion (in both directions $x$ and $y$) with diffusion constants $D_A$ and $D_B$ respectively. The Gray-Scott model assumes that molecules of species $A$ diffuse faster that those of $B$ ($D_A>D_B$).

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 $\Delta t$ and distance $\Delta x$ and $\Delta y$ intervals.

In discrete form, using the Euler algorithm, we can write

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 $\Delta x = \Delta y = 1$

Which can be re-written, introducing a $3\times 3$ 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 $3\times 3$ 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 $t$, we know the concentrations of both species for all points of the square,

then, we can calculate the concentrations for all points at time $t+\Delta t$,

using the equations

The parameters I have used in the simulation presented in Figure 2 are

Parameter Value Description
$k$ $1.0$ rate of producing species $B$
$f$ $0.033$ birth (and death) rate for species $A$
$g$ $0.062$ $g+f$ is the death rate for species $B$
$D_A$ $1.0$ diffusion coefficient for species $A$
$D_B$ $0.2$ diffusion coefficient for species $B$
$T$ $3000$ total time sampled
$\Delta t$ $1$ time increment
$L$ $100$ total distance sampled in each direction
$\Delta x$ $1$ space increment (both in x and y directions)

We assume that we have a lattice with $100\times 100$ points. At each point $(x,y)$ and at each time $t$, there is one concentration for each species $A(x,y,t)$ and $B(x,y,t)$.

The initial conditions (Figure 2 at $t=0$) correspond to having $A=1$ everywhere in the square (in blue), and $B=0$ everywhere in the square expect for four small (magenta) areas in which $B=1$. 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:

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

Consider a general reaction with two components $A$ and $B$, as in the case of the Gray-Scott model. Let’s consider which are the conditions under which a stable set of reactions between $A$ and $B$, 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 $A_0(x,y,t)$ and $B_0(x,y,t)$ of the reaction-only system

that is, the steady-state solution is characterized by

this solution $(A^\ast, B^\ast)$ is stable if the Jacobian matrix $J$ defined by

is such that the real part of the eigenvalues is negative.

For a $2\times 2$ 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 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 $A$ and $B$, subject to a reaction-diffusion process, their concentrations depending on time and space $A(x,y,t)$ and $B(x,y,t)$ satisfy some equations

We assume that we have found some stable solution of the non-diffusion process $A^\ast$ and $B^\ast$ 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 $\mathbf{u}(x,y,t)$ has a constant term $\mathbf{u_0}$, and such that the dependencies in time and space are given by

for some positive constants $\alpha$ and $\beta$.

Notice that we are assuming that both reactants $A$ and $B$ have the same dependencies in time and space. This is an important simplification for the results that we will find.

Under these assumptions for $\mathbf{u}(x,y,t)$, the linear differential equation becomes

or

And since this results has to be true for all values of $t$ and $x,y$, we have

This has now become a 2D eigenvalue problem for the 2D matrix

$\mathbf{u_0}$ has two solutions, which are the two eigenvectors of $M$ corresponding to its two eigenvalues

where the trace of $M$ is

and the determinant of $M$ is given by

The eigen values $\alpha_{1,2}$ are the two solutions to the quadratic equation in $\alpha$ resulting from imposing the condition

For a solution to be stable we need two conditions

• (a) The solution is real, which requires $det(M) > 0$

• (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 $M_0 = M(\beta=0)$, 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 $Tr(M)$ is negative as long as $T(M_0)$ is negative.

While $det(M_0) > 0$, diffusion has to be capable of achieving the condition $% $ in order to produce instability.

Notice that

The $det(M)$ is a parabola as a function of $\beta^2$, and it is positive for $\beta = 0$.

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 $det(M)$ at the minimum is

The $det(M)$ being negative at the bottom of the parabola then requires

• $\beta^2_{min} > 0$ or
$$D_B F^\ast_{AA} + D_A F^\ast_{BB} > 0.$$
• $% $ or

which can be re-written 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.