May 21, 2018, Monday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links

MBW:Pre-Pattern Formation Mechanisms for Animal Coat Markings

From MathBio

Jump to: navigation, search


In the following discussion, coat patterns of various animals are studied. By treating the system as a reaction diffusion model, realistic coat patterns can be created. Using both analytic and numerical solutions to the diffusion equation and imposing various boundary conditions, simulations of the various possible coat patterns are produced. This process includes nondimensionalizing and linearizing a set of differential equations and examining the stability of the analytic solutions, as well as the implementation of a finite difference method in order to obtain a numerical solution.

Categorization of Project

Mathematics Used: Finite differences, Numerical Stability, Stability theory with Linear Systems, Partial Differential Equations, Neuman Boundary conditions

Type of Model: Reaction Diffusion model with two chemical species

Biological System Studied: Animal coat pattern formation, Cows, Cheetahs, Zebras


Various animal coat markings show great variation and diversity from both species to species and individual to individual. It has been shown that by using a reaction diffusion model it is possible to qualitatively reproduce the coat markings found in nature on a wide variety of animals [1]. By setting up a reaction diffusion system for two chemical species such as an inhibitor and activator for melanin production it is possible to gain valuable insight about how animals produce their coat patterns. In this article we discuss both methods of obtaining analytic and numerical solutions for coat patterns. In addition, it may be possible to gain insight into other possible coat patterns that have never been observed before based on given chemical properties and initial distributions.


It is known that in mammals coloration is most often attributed to the presence of melanin in the epidermal layer of the skin or in the hair and are contained in pigments cells known as melanocytes. These melanocytes do not occur in a homogeneous concentration, and it is believed that this inhomogeneity is what leads to color variations in mammals[2]. The first major paper published on the topic of reaction diffusion in a biological application was done by A. Turing in 1952 [3] where the idea of using reaction diffusion as a mechanism for transforming the homogeneous distribution of morphogens to a pattern due to some form of an instability or perturbation from an unstable homogeneity was presented and supported. This publication did not directly discuss coat markings though it formed the framework on which more research could be done in the following decades. During the years 1976-1981 there was a vast amount of experimentation and mathematical inquiry performed to study and understand reaction diffusions systems and their applications[2]. Based on the wealth of knowledge from this spurt of progress in the field of reaction diffusion systems in 1981 J. Murray was able to publish “A Pre-pattern Formation Mechanism for Animal Coat Markings” in J. theor. Biol. [2]. In this paper Murray adapts the reaction diffusion equation to address how an activator and inhibitor system could produce the multitude of coat patterns found in mammal coats. In his paper Murray proposes that the system he describes is capable of producing qualitatively similar patterns to those seen in nature. Murray also proposes that this is likely the true mechanism, as from an evolutionary standpoint a simpler mechanism like reaction diffusion is the most likely candidate to be the true mechanism driving the formation of these complex coat patterns.

Two Species Reaction Diffusion

To form many of the coat patterns seen on various mammals we can study a system consisting of two species that are diffusing and reacting with one another on a two dimensional plane. For simplicity we will now consider this system existing on a rectangular domain. The general form of the equations is:

{\frac  {\partial S}{\partial t}}=U(S,A)+D_{{S}}\nabla ^{{2}}S
{\frac  {\partial A}{\partial t}}=V(S,A)+D_{{A}}\nabla ^{{2}}A

Where D_{{A}} is and D_{{S}} are the diffusion coefficients for chemical inhibitor and activator respectively, U(S,A) and V(S,A) are the reaction rates of the different chemicals, and \nabla ^{{2}}S and \nabla ^{{2}}A are the diffusion terms for the respective chemicals. By selecting specific solutions for this system we are able to construct any two tone pattern, due to the solutions being a sum of terms that make up a complete basis set.

Mathematical Analysis

Using the framework built in the previous section we are able to begin to study how reaction diffusion can lead to realistic coat pattern formation. In this example we will study the pattern formation on a rectangular region with Neumann Boundary Conditions. The nondimensionalized system is:

{\frac  {\partial s}{\partial t}}=g(s,a)+\nabla ^{{2}}s,
{\frac  {\partial a}{\partial t}}=f(s,a)+\beta \nabla ^{{2}}a,

where a and s are the chemical concentrations of the nondimensionalized inhibitor and activator chemicals respectively, \beta is the ratio of the diffusion constants for a and s, and g and f are the nondimensionalized reaction rates for a and s. The specific forms of the reaction terms g and f are:

g(s,a)=\gamma [s_{{0}}-s-\rho F(s,a)],
f(s,a)=\gamma [\alpha (a_{{0}}-a)-\rho F(s,a)],
F(s,a)={\frac  {sa}{1+s+Ks^{{2}}}},

where \gamma , \rho , \alpha , and K are all positive, dimensionless constants. The F(s,a) term is the rate by which the given chemicals are consumed by the reaction and is closely related to Michaelis-Menton Kinetics. The \gamma term is acting as a spacial scaling parameter on the domain for our simulation. For more information on the derivation of this method see Appendix A in Murray (1981) [2] or the original derivation by Thomas and Kernevez [4].


It is proposed that the spatial patterns are a result of diffusion driven instability of the uniform steady state solutions {\tilde  {s}} and {\tilde  {a}}. The steady state solutions are found by letting g(s,a) and f(s,a) be zero and solving for {\tilde  {s}} and {\tilde  {a}}. Murray discusses that only a certain range of values of the inhibition parameter K will result in diffusive instability due to small spatially inhomogeneous perturbations [2]. It is also important that the parameter \beta >1 so that chemical a diffuses quicker than chemical s. For the purpose of the discussion of the paper assume all values of K and \beta are selected to ensure diffusive instability. After introducing a perturbation of the steady state the instability in the system allows for the evolution into a heterogeneous steady state that gives the spatial structure to coat patterns.

To study the stability lets analyze a rectangular grid with 0\leq x\leq a and 0\leq x\leq b and Neumann boundary conditions. The solutions to a small perturbation u(x,y,t)=s(x,y,t)-{\tilde  {s}} are given by

u(x,y,t)=\sum _{{m}}\sum _{{n}}U_{{m,n}}e^{{\lambda _{{m,n}}t}}cos({\frac  {n\pi }{a}}x)cos({\frac  {m\pi }{b}}y)

Where U_{{m,n}} are the constants that describe the magnitude of each eigenfunction and \lambda _{{m,n}} are the eigenvalues of the solution. If at least one \lambda _{{m,n}}>0 then the solution is system is considered unstable to the perturbation and will result in a pattern formation. A similar analysis can be made with respect to a(x,y,t) and it's solutions.


In appendix B of his paper Murray derives conditions on m, n, a, b for a given set of parameters that ensures at least one positive eigenvalue, hence leading to instability in the system. We will briefly discuss his methods here. To find the unstable modes for the reaction diffusion equations above we will linearize the equations about the stable points {\tilde  {s}} and {\tilde  {a}} and to do so we will use the matrices below:

{\mathbf  {w}}=\left({\begin{array}{c}s-{\tilde  {s}}\\a-{\tilde  {a}}\end{array}}\right),
D=\left({\begin{array}{cc}1&0\\0&\beta \end{array}}\right),
M=\left({\begin{array}{cc}{\frac  {\partial g}{\partial s}}&{\frac  {\partial g}{\partial a}}\\{\frac  {\partial f}{\partial s}}&{\frac  {\partial f}{\partial a}}\end{array}}\right)_{{{\tilde  {s}},{\tilde  {a}}}}.

Using this vectorized form of the system we are able to linearize the system to take the form:

{\frac  {\partial {\mathbf  {w}}}{\partial t}}=M{\mathbf  {w}}+D\nabla ^{{2}}{\mathbf  {w}}.

If we let {\mathbf  {w}}({\mathbf  {r}},t)=e^{{\lambda t}}{\mathbf  {W}}(r) the values of \lambda can be determined using:

\left|M-k^{{2}}D-\lambda I\right|=0

We then take the eigenvalues of {\mathbf  {W}} to be represented by k^{{2}} and are such that they satisfy the Helmholtz equation below:

\nabla ^{{2}}{\mathbf  {W}}+k^{{2}}{\mathbf  {W}}=0,

with the Neumann Boundary Conditions.

In the two dimensional case we are interested in we specifically are looking for solutions to:

{\frac  {\partial ^{{2}}{\mathbf  {W}}}{\partial x^{{2}}}}+{\frac  {\partial ^{{2}}{\mathbf  {W}}}{\partial y^{{2}}}}+k^{{2}}{\mathbf  {W}}=0

This equation has a known Fourier series solution and the k^{{2}} are given by the following:

k_{{mn}}^{{2}}={\frac  {m^{{2}}\pi ^{{2}}}{b^{{2}}}}+{\frac  {n^{{2}}\pi ^{{2}}}{a^{{2}}}}

By requiring at least one value of \lambda such that its real component is greater than 0 under the presence of diffusion we can show only certain modes of the analytic solution will cause instability. The formula is quite involved and is explained in greater detail in Appendix B in Murray[2]. In our simulation section will we be concerned with a set of parameter values that ensure the following restriction holds:

0.036\gamma <{\frac  {m^{{2}}}{b^{{2}}}}+{\frac  {n^{{2}}}{a^{{2}}}}<0.115\gamma .

Examples Using Simulations

By using the analysis in the stability section we can check to see if our simulation is correct as well as demonstrate how size and time can affect the formation of pattern. To do this we will create a rectangular grid with the x dimension larger than the y dimension, but the spacing between grid points equal in all directions. Here we will stick with the notation of the original paper labeling the x-direction as a (not to be confused with the substrate a) and the y-direction as b. We use the same value for parameters as given in the paper: \alpha =1.5, \rho =13, s_{{0}}=102, a_{{0}}=77, K=0.1, \beta =5 and \gamma =1. Using these parameters the condition for there to be at least one unstable mode is:

0.036\gamma <{\frac  {m^{{2}}}{b^{{2}}}}+{\frac  {n^{{2}}}{a^{{2}}}}<0.115\gamma .

Allowable Modes

If we let the a = 3 and b = 1, hence creating a rectangular region with an aspect ratio of 3, then using the above equation we can see that the only instability inducing, valid modes of the analytic solution are m=0, n=1. Running our simulation long enough for a pattern to form, we obtain the following:


Any value that was above the steady state solution is represented as white. This is in agreement with the analytical solution where one would expect the perturbed solution to be

u(x,y,t)=e^{{\lambda _{{01}}t}}cos({\frac  {1*\pi }{3}}x).

If we change the boundaries to a = 10, b= 1 we can expect that the allowable modes will become m=0, n=1,2 or 3. The resulting pattern after a simulation is:


Remember that the color scheme is chosen so that values above or below the steady state represent white and black respectively. The solution is a sum of cosine terms, the sine terms drop out as they cannot satisfy the Neumann boundary conditions, but the resulting pattern is only concerned with the sign of the variations from the steady state, regardless it's amplitude.

Time Effect

We noticed that the random perturbation we applied quickly vanished and a pattern began to form. The rate of change of the analytic solution with respect to time is one of exponential decay, so it follows that changes to the pattern are made relatively early in the simulation with only refinements to the pattern made later in the simulation. In Murray's paper [2] it is suggested that the pattern formation of the fetus takes place during some time frame in the womb. We can simulate different time periods where reaction-diffusion would take place by simply limiting the number of time steps we allow our system to evolve. To show this we ran our a=3, b=0 simulation varying the total amount of time steps.

T200.jpg T1000.jpg

t=200, t=1000

T3000.jpg A3b1pattern.jpg

t=3000, t=5000

It is clear that different patterns are formed for different lengths of time. The physical analogy would be if two comparable sized animals had similar shapes but different gestation lengths then you might expect them to have qualitatively different coat patterns.

Size Effect

Murray described that the size of the fetus will affect the pattern formation. From the stability analysis we noticed that if a and b are increased more modes are permitted that can cause diffusive instabilities. To demonstrate this we ran four simulations with equal time steps (t=1000), varying a and b while keeping the aspect ratio of the region the same. The results of these various simulations are shown here:

T1000.jpg A6b2.jpg

a=3, b=1: a=6,b=2

A12b4.jpg A18b6.jpg

a=12, b=3: a=18,b=6

Clearly different patterns are formed with different sized domains. In nature we notice that different sized animals have different qualitative patterns. In fact if the domain is too small, no modes would cause instability. To test this we ran the simulation with a=1, b=1 which resulted in no pattern formation. For obvious regions we do not show an image of this. We could compare this to mice in which we do not generally notice spotted coat patterns.

Simple Examples Using Analytical Solution

Zebra Stripes

Referring to the analytical solution of u(x,y,t) we could form the fully evolved patterns by fixing the constants and allowed modes. A pattern that we are able to simply create with only one U_{{m,n}} term is that of the stripes on a zebra. Setting U_{{20,1}}=1and setting all other terms to zero we are able to generate the plot below which is a highly simplified. To visualize how this pattern goes on a zebra we consider the horizontal center line to be the spine of the zebra, then fold the two halves of the plane to cover the sides of the animal.


A pattern similar to that of a zebra coat but more complex and realistic is created by just setting U_{{9,4}}=U_{{8,8}}=1. If you think about taking this image and mapping it to the body of a zebra the same way as done in the previous example it can be seen how this coat pattern could go on a zebra.


More Complex Pattern Formation

These patterns are interesting novel examples, but to produce more complex patterns the use of numerical simulations is of great significance. In his paper[2] Murray has a number of realistic examples of these patterns based on the changing of the diffusion constant \beta as well as by changing the size parameter \gamma . A few examples of this are shown below are such that all constants are held constant except \gamma which is different for each example and shown in the plot (for other settings see [2]).


We experimented with creating numerical simulations to reproduce these results and our findings are outlined below.

Numerical Simulation and Results

Finite Differencing , Boundary Conditions, and Numerical Stability

The next step after being able to analytically predict patterns is to develop a numerical algorithm by which the reaction diffusion system can be tested. To do this we used a discretization technique on the system:

{\frac  {\partial s}{\partial t}}=g(s,a)+\nabla ^{{2}}s,
{\frac  {\partial a}{\partial t}}=f(s,a)+\beta \nabla ^{{2}}a.

where f and g are defined above in the generic form section. To create a finite difference method we used a centered difference for the second derivative in the spatial dimensions and a forward difference for the first derivative in the temporal dimension shown below in one dimension.

{\frac  {\partial ^{{2}}s}{\partial x^{{2}}}}={\frac  {s_{{i+1}}+s_{{i-1}}+s_{{i}}}{h^{{2}}}} (Centered Difference)
{\frac  {\partial s}{\partial t}}={\frac  {s_{{i+1}}-s_{{i}}}{{h}}} (Forward Difference)

We then set our parameters to the same values as given in Murray [2] which are also described in our analytic solution section. In addition, we applied Neumann boundary conditions to the boundaries of the region by setting the boundaries to have the value of their nearest neighbor in the direction of the interior of the region. This is justified by the fact that along the boundary we have the derivative equal zero, so in one dimensions it can be shown that:

{\frac  {\partial s}{\partial x}}={\frac  {s_{{i+1}}-s_{{i}}}{h}}=0\Rightarrow s_{{i+1}}=s_{{i}}.

The final consideration when setting up the numerical simulation is the stability based on spatial and temporal time steps. This stability condition is known as Courant–Friedrichs–Lewy (CFL) condition and states that for a finite difference simulation of a second order system that the following condition must be met,

{\frac  {\Delta t}{\Delta x^{{2}}}}<{\frac  {1}{2}}.

So when constructing our simulations we were careful to ensure that this condition was met to ensure stability as we move our simulation forward in time.

Simulation Results

The type of pattern resulting from our simulation are controlled by the value we set for \gamma as well as the aspect ratio between the of our simulation domain. \gamma acts as a scaling value for size, so as \gamma goes up the size of our redimensionalized domain goes up by a factor of \gamma ^{{2}}. So by experimenting with both altering \gamma and the aspect ratio of our domain we were able to arrive at various interesting results. Three of these results that correspond to well known animals are described here. It is important to remember that the spacial distribution of the chemicals is a smooth curve, but we use a binary color scheme to show contrast between locations where pigment can or cannot form.


To create the pattern found on Holstein Cattle as shown here we used an aspect ratio of one (i.e. a square region) and a intermediate value for \gamma of 250. The aspect ratio of one should allow for the same eigenfunctions to be stable in both the x and y directions. The resulting pattern consists of rather large sections of either black or white and no small spots or stripes form in the domain. CowSim.jpg


To generate the cheetah coat pattern we altered our \gamma value to be 2500, 100 times larger than the \gamma value used for our cow. Again we have an aspect ratio of 1 but have refined the grid size slightly to reduce the gain in the image. Similar patterns formed for our simulation with a more coarse grid size, but the higher resolution makes the pattern have smoother and more realistic boundaries. Our simulation created a large number of similarly sized spots in a uniform distribution across the coat, with a few sections where spots had not yet fully separated from one another resulting in short localized stripes.



The final pattern that we formed was that of zebra stripes. To do this we adjusted both \gamma and the aspect ratio, setting them to 10 and 50 respectively. This resulted in a domain that is long and thin. This aspect ratio forces our solution to prefer having stable horizontal eigenfunctions as opposed to vertical ones as there is more space in which they can form, hence resulting in the striping pattern seen below.


Future Work

The work presented here can be used as a framework on which to base more advanced work. One possible project would be to use the same numerical scheme described above but applied to an animal shaped domain as shown below. The boundary conditions are not Neumann along the entire boundary, but rather the head and body have Neumann boundary condition, the limbs and tail have periodic boundary conditions along their length, and the hoofs and the tail’s end have Neumann conditions again. The pattern shown on this animal shape is not numerically derived but is rather an example of a pattern which could possibly form on the domain.

AnimalDomain.jpgbe wrapped around to form a conic section which could represent a tail or tapered limb. An example from Murray of this type of domain is shown below. It is interesting to note how as the domain tapers down in size the same pattern that form spots begin to form stripes as the spot is in fact wrapping all the way around the tail to form stripes.



Based on the results in the previous sections we are able to see how basic coat patterns are formed from only considering a reaction diffusion system with two chemical species. These results can be extrapolated to form any pattern based on a two color scheme since our solution is a sum of basis vectors which form a complete basis. In addition, we were able to produce realistic results which agreed with those of Murray using a finite differencing scheme. Using what has been shown here it is also possible to map the spot or stripping patterns found on the tails of various animals by taking the same solution and having them transformed onto a conic surfaces to recreate tail or limb patterns. To do this numerically we would need to create a non-rectilinear grid on which our finite differencing method could be performed.

Further Reading

In Murray’s paper, mammalian coat patterns are analyzed as a reaction-diffusion system. Some, however, disagree that this is the best modeling method. In the paper “Integrating Shape and Pattern in Mammalian Models,"[5] a method called Clonal Mosaic (CM) is used instead. This model postulates that coat patterns are just an arrangement of cells that all originated from the same source; different colors simply arise from different cells.

Murray’s paper is cited within the context of coat pattern creation. The pattern is created in two different phases. In the first phase, both growth and pattern formation occur. This happens during the fetal stage of the animal. In the second phase, the pattern has been formed and can now only be changed due to the growth of the animal. This happens both before and after the animal is born. In his paper, Murray provided an estimate of how long the pattern formation phase lasted. This was based on a reaction-diffusion model, however, and may or may not apply to different models like clonal mosaic. Instead of using Murray’s estimate, this model assumes the first phase of pattern creation for giraffes (the model animal used) lasts about half of the gestation period.

Murray's paper is also cited in the paper "Generation of Diverse Biological Forms through Combinatorial Interactions between Tissue Polarity and Growth," by Richard Kennaway, Enrico Coen, Amelia Green, Andrew Bangham. Kennaways article covers how genes control the formation of a Snapdagon flower. The theory they are showing is a polarity based system. This system creates signals based on local polarities. Kennaway's paper cites Murray's paper in the introduction section stating that there have been many studies and papers discussing how the patterns are formed, but there paper is different because the model also depends on mechanics of the neighboring regions. [1]


  1. N. Britton, Essential Mathematical Biology. (2003).
  2. 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 J. Murray, "Pre-Pattern Formation Mechanisms for Animal Coat Markings". 'J. theor. Biol' 88 (1981) 161-199.
  3. A. Turing, “The Chemical Basis of Morphogenesis”. 'Phil. Trans. R. Soc. Lond. B' 237 (1952). 37-72.
  4. D. Thomas and J. kernvenez "Numerical Analysis and Control of some Biochemical Systems". 'Applied Mathematics & Optimization' 1 (1976) 222-285.
  5. Marcelo Walter , Alain Fournier , Daniel Menevaux, "Integrating shape and pattern in mammalian models", Proceedings of the 28th annual conference on Computer graphics and interactive techniques, p.317-326, August 2001.