2D Reaction-Diffusion (Gray-Scott) - Feed/Kill Map and Turing Patterns

Presets

Live readouts

Graphs

Feed/kill map: click or drag to choose f, k (green = Turing unstable, amber = oscillatory, blue = stable, dark = no nontrivial steady state)

Linear growth rate sigma(q) of a perturbation, and the measured wavenumber

Mean and contrast (std) of chemical B against time

Radial power spectrum of the pattern (peak = pattern wavenumber)

Numerical checks

running...

How to use

Pick a preset or click in the feed/kill map. Paint chemical B with the mouse (left button) or erase it (right button or the Erase mode) directly on the dish. Keys: Space pause, R restart.

Presets: spots, mitosis (one seed that keeps dividing), stripes, coral, holes, moving spots, U-skate and a linearly unstable Turing case.
Layers: pick chemical B or A, a colour map, lighting (shades the concentration as a height field) and contour lines at v = 0.15.
Boundary: the dish is periodic (a torus), so patterns wrap around the edges.
The spectrum graph is the radial power spectrum of the pattern; its peak is compared with the linear-theory wavenumber.

The science

Two chemicals A (concentration u) and B (concentration v) diffuse and react. The Gray-Scott model (1983, used by Pearson 1993) is

du/dt = Du lap(u) - u v^2 + f (1 - u) dv/dt = Dv lap(v) + u v^2 - (f + k) v

The term u v^2 is the autocatalytic step A + 2B -> 3B. A is fed in at rate f (towards u = 1) and B is removed at rate f + k. Defaults: Du = 0.2097, Dv = Du/2, lengths in pattern units (one lattice cell of the 2D page is 1).

Steady states. The trivial state is u = 1, v = 0. A nontrivial state satisfies u v = f + k and (f+k) v^2 - f v + f (f+k) = 0, so it exists only if f >= 4 (f+k)^2 (the dashed white curve on the map). With the Jacobian J of the reaction terms at that state:

J = [[-v^2 - f, -2uv], [v^2, 2uv - (f+k)]]

Linear stability (Turing). A mode cos(q x) grows at the rate sigma(q) = largest real part of the eigenvalues of J - q^2 diag(Du, Dv). If sigma(0) < 0 (stable without diffusion) but sigma(q) > 0 for some q, diffusion destabilises the state: a Turing instability. It needs different diffusivities (here the activator B diffuses more slowly than A). The fastest mode q* predicts the wavelength 2 pi / q*; the graph shows sigma(q) and the measured wavenumber.

Surfaces. On a curved surface the Laplacian becomes the Laplace-Beltrami operator. The mesh uses the cotangent formula with lumped vertex areas:

(lap f)_i = (1 / A_i) sum_j (cot a_ij + cot b_ij) / 2 * (f_j - f_i)

where a, b are the angles opposite edge ij and A_i is one third of the area of the triangles around vertex i. Time stepping is explicit Euler with dt chosen from the Gershgorin bound (stable for every mesh here). The slider h sets how many pattern units one mesh edge represents, i.e. how big the surface is compared with the pattern wavelength.

Numerical checks (run in the page). (1) a 1-D ring of the full nonlinear PDE is perturbed by 1e-9 along the unstable eigenvector and the growth rate of the Fourier mode is measured and compared with sigma(q); (2) the cotangent operator applied to the vertex coordinates must equal the exact mean-curvature vector, lap X = -(k1 + k2) n (sphere: -2x/R^2; torus: exact principal curvatures); (3) the Laplacian of a constant vanishes. The measured pattern wavelength is the Rayleigh-quotient estimate q^2 = <dv, -lap dv> / <dv, dv>.

Honest simplifications. Most of the famous patterns (spots, mitosis, coral, holes) live where no nontrivial steady state exists, or where it is stable: they appear by finite-amplitude invasion of the trivial state, so the linear wavelength is only a guide there (the measured value is typically within a few tens of percent and drifts as spots divide or merge). Only the green region of the map is a genuine Turing instability. Mesh edge lengths are not perfectly uniform (especially on the torus), explicit Euler is first order in time, and chemical concentrations are dimensionless. The two surfaces also differ topologically: stripes on a sphere cannot all line up (hairy-ball theorem), a torus can carry parallel stripes.

References: P. Gray and S. K. Scott, Chem. Eng. Sci. 38 (1983) 29; J. E. Pearson, Science 261 (1993) 189; A. M. Turing, Phil. Trans. R. Soc. B 237 (1952) 37; R. Munafo, mrob.com/pub/comp/xmorphia.