In the last few decades, developmental biologists have extensively used the reaction-diffusion model to explain pattern formation in living organisms. The original model was proposed by Alan Turing in 1952 . The model is based on the idea that pattern formation results from two fundamental mechanisms: (1) coupled catalytic and autocatalytic reactions in a space element between two chemical species, an activator and an inhibitor, and (2) transfer of the interacting species to and from the neighboring space elements through a diffusional transport mechanism. Under appropriate reaction and diffusion conditions, a periodic pattern is formed from an initially homogeneous spatial distribution of activator and inhibitor [2, 3]. Examples of pattern formation can be found in biology, chemistry (the famous Belousov–Zhabotinskii reaction), physics, and mathematics [4, 5].
To illustrate the mechanism of pattern formation, consider the hypothetical activator-inhibitor reaction sequence:
The species , , and are supposed to be sufficiently abundant in the reaction mixture that their respective concentrations , , , and can be considered constant. The activator reacts with species to produce more by the autocatalytic reaction R1 (so that is reactant, catalyst, and product), but it also promotes the production of the inhibitor by the catalytic reaction R2, in which is a catalyst. Both activator and inhibitor decay with time (reactions R4 and R5).
Because of the equilibrium reaction R3, the concentrations of and are such that . Denoting by the reaction constant of reaction , the reaction rates of and are:
From the expression of the reaction rate of , , the species inhibits the activator production, hence its name: the larger its concentration, the lower the production rate of .
Introducing the variables and , the governing equations of the reaction-diffusion system in 2D can be written in the following nondimensional form:
with Set . For the formation of spatial patterns, the diffusion rates of activator and inhibitor should be very different: set the diffusion coefficients to be and .
For the numerical solution of the system of ODEs, assume (1) periodic boundary conditions as well as (2) the initial conditions:
where and are numbers taken randomly in the interval .
The Chebyshev orthogonal collocation method applied with nodes in both spatial directions transforms the system of two coupled nonlinear PDEs into a system of 392 nonlinear coupled first-order ordinary differential equations. This system of ODEs is solved using the built-in Mathematica command NDSolve. The snapshots show the results for the 2D inhibitor concentration distribution at times and at . It is interesting to note, in particular, the emergence of Turing patterns similar to leopard spots in the concentration distribution at .