Reaction Diffusion Equations and Animal Coat Patterns

Humans have long had a fascination with the coat patterns of animals. Alan Turing was the first to articulate an explanation of how the patterns of animals like leopards, jaguars and zebras are determined. Turing asserted that the patterns can arise as a result of instabilities in the diffusion of morphogenetic chemicals in the animals' skins during the embryonic stage of development.

Let C be the vector of morphogen concentrations. The vector equation giving the spatial and temporal dynamics of morphogen concentrations is:

∂C/∂t = F(C) + D∇²C

where F(C) is a nonlinear vector function and D is diagonal matrix.

The equations can be put into a form such that the diagonal matrix D is of the form
D = diag(1,d). Thus d represents the relative magnitude of the diffusion coefficient of one morphogen compared to the other.

Turing's Principle

Turing said in effect that if D=0 and C goes to the steady state solution of F(C*)=0, then under some circumstance the introduction of nonzero D can create spatial variations that produce patterns.

For the problem to be properly formulated there must be boundary and initial conditions specified. Let B be the domain of the problem and let ∂B be the boundary of the domain. Then the boundary and intial conditions are:

n·∇C = 0 on ∂B
C(r,0) = G(r).

Let C* be a solution to F(C)=0. C* is a solution which is uniform over space and thus ∇²C* = 0. Now consider a linearization of the equation about the point C*. Thus

∂C/∂t = F(C) + D∇²C
∂C*/∂t=F(C*)+D∇²C*
∂c/∂t = γAc+ D∇²c

where c = C - C* and γA = [∂F/∂C], where γ is a nondimensional scale parameter that is proportional to the square of the linear dimension of a one dimensional system and the area of a two dimensional system. The existence of γ and its relation to the scale of the system requires further analysis.¹

In the absence of spatial variation in c the behavior of the system for |c| small is the same as the solution to the system:

∂z/∂t = γAz

This system is stable if all of the roots of the determinant equation

det[λI - γA] = 0

have negative real parts. In general, the product of all the roots of this equation when A is n×n is equal to γⁿdet[A] and the sum of the roots is equal to γtr(A).

For a two-component system the equation for the eigenvalues is:

λ² - γtr(A)λ + γ²det(A)= 0

The condition that the real parts of both roots be negative is that the det(A) be positive and the trace of A be negative. These conditions are equivalent to both roots having the same sign and the sum of the roots being negative.

Now let us go back to the equation that includes spatial variation and consider a separation of variables; i.e., c(r,t)= W(r)T(t) where W(r) is a scalar-valued function and
T(t) is a vector function. Thus,

W(r)∂T/∂t =
W(r)γAT(t) + ∇²W(r)DT(t)

In order to proceed further it is necessary to limit W to functions which are solutions to the boundary value problem:

∇²W(r)= -k²W(r)
r·∇W = 0 for r on ∂B.

Now the system of equation is

W(r)[∂T/∂t - γAT + k²DT] = 0

For the solution for T(t) to be of the form exp(λt)T, we must have

W(r)[λI - γA + k²D]T = 0

and thus

det[λI - γA + k²D] = 0

The λ's are the eigenvalues of the matrix M = γA - k²D. For spatial instability one of the eigenvalues for some value of k² must have a positive real part. For k² = 0 none of the eigenvalues has a real part.

For a two-component system the equation satisfied by the eigenvalues is:

λ² - tr(M)λ + det(M) = 0

In order for one of the solutions of the quadratic equation to have a positive real part for some value of k it is necessary that either:

• the coefficient of λ is negative
• the constant term is negative

Since -tr(M) = -γtr(A) + k²tr(D) and tr(A) was assumed to be negative for steady state stability there is no way that the coefficient of λ, -tr(M), can be negative. Therefore the only way that at least one root can have a positive real value is that det(M) be negative for some value of k².

The graph shows the case that det(M) is negative for some value of k². For this case to hold there must be a value of d such that for some value of k² det(M)=0 and the derivative of det(M) with respect to k² is also zero. The value of d such that det(M) is tangent to the k² axis is called the critical value of d, d_c.

Natural Patterns and Simulations

For a periodic solution in an essentially linear domain, as for example in a snake skin, coloration of all regions above a threshold level would result in stripes. In more rectangular domains there could also be stripes, but there can also be spots as shown on the accompanying reproductions from J.D. Murray's book.

Murray's method of solution makes use of solutions to the Helmholtz equation,

∇²W(r)= -k²W(r).

He presents the patterns that occur on the special domains shaped as hexagons, rhombuses and triangles. The results shown in the following are evocative of types of patterns that occur on reptiles.

The most striking of Murray's simulations are those for a domain shaped like an animal skin with four legs, a head and a tail. He examines what the equations give for the junction of a horizontal stripe pattern on a leg with a vertical stripe pattern on the torso. The results are consistent with the pattern of zebra stripes at the junction of leg and torso.

Proto-types of Patterns of Animal Tails

Source: J.D. Murray, Mathematical Biology

Simulations of Animal Coat Patterns

Source: J.D. Murray, Mathematical Biology

References

1. James D. Murray, Mathematical Biology, 2nd Corrected Edition, Springer, 1993.
2. J.D. Murray and M. R. Myerscough, "Pigmentation Pattern Formation on Snakes," Journal of Theoretical Biology, vol. 149 (1991), pp. 339-360.
3. J.D. Murray, "Parameter Space for Turing Instability in Reaction Diffusion Mechanisms: A Comparison of Models," Journal of Theoretical Biology, vol. 98 (1982), pp. 143-163.