Mathematics · Chemistry · Pattern Formation

Gray-Scott Reaction-Diffusion

Watch Turing patterns self-organise from two coupled PDEs — coral reefs, stripes, mitosis blobs and labyrinthine waves emerge from microscopic chemical rules.

Preset  
F      
k      
Frame   0

Gray-Scott Reaction-Diffusion

In 1952 Alan Turing published The Chemical Basis of Morphogenesis, proposing that a pair of reacting, diffusing chemicals could spontaneously break spatial symmetry and produce the patterns we see in nature — from leopard spots to zebra stripes. The Gray-Scott model (1983) is the cleanest two-species realisation of this idea:

dU/dt = Du·∇²U − U·V² + F·(1−U)
dV/dt = Dv·∇²V + U·V² − (F+k)·V

Species U (activator, shown in low concentrations) is continuously fed into the domain at rate F. Species V (inhibitor) is autocatalytically produced from U (the U·V² term is the irreversible reaction U + 2V → 3V) and killed at rate F+k. The ratio Du/Dv = 2 ensures V diffuses twice as fast — the classic Turing condition for creating structure.

The feed-kill parameter pair (F, k) acts as a 2-D phase space. By sweeping it you can find coral-like growths, self-replicating spots, oscillating worm-like domains, travelling solitons and labyrinths.

Turing Instability

A uniform steady state becomes unstable to small spatial perturbations when Du/Dv is large enough, triggering spontaneous pattern formation at a characteristic wavenumber.

Self-Organisation

No external template guides the pattern. Only initial random noise, local chemistry, and the diffusion asymmetry determine whether spots, stripes, or waves appear.

Phase Diagram

Pearson (1993) mapped the (F, k) plane into 12 qualitatively different regimes. The simulator lets you traverse this landscape in real-time.

Biological Relevance

Turing mechanisms are now confirmed in fish skin pigmentation (Ocelots, Zebrafish), finger digit spacing, palatal rugae formation, and hair-follicle patterning.

Preset Chemical Regimes

Each preset corresponds to a distinct region of the (F, k) phase space. The Du and Dv values are kept at their standard Gray-Scott ratio of 2:1.

PresetF valuek valuePattern TypeBiological Analogue
Coral Reef0.05450.0620Dendritic coral-like networksCoral polyp growth morphology
Spots0.03500.0650Stable polka-dot domainsJaguar / cheetah spots
Stripes0.04000.0600Quasi-periodic stripe bandsZebra stripe pigmentation
Mitosis0.03670.0649Self-replicating spot divisionCell division topology
Solitons0.03000.0620Travelling solitary V pulsesNerve impulse propagation
Labyrinthine0.06300.0609Branching maze-like tunnelsMammalian cortical folding

Numerical Method

For a 256×256 grid the simulator uses explicit forward-Euler time integration with periodic boundary conditions:

Δt = 1.0  (non-dimensionalised)
∇²U[i,j] ≈ U[i+1,j] + U[i−1,j] + U[i,j+1] + U[i,j−1] − 4·U[i,j]

Two Float32Array buffers are ping-ponged to avoid in-place updating artefacts. Each animation frame runs 12 substeps to achieve the visual impression of continuous evolution. Values are clamped to [0, 1] to prevent blow-up near the Turing instability boundary.

A V-concentration lookup table (256 entries × RGBA) maps the scalar field to colour without per-pixel branching, allowing smooth interactive framerates on a 65 536-cell grid.

Grid Size

256 × 256 = 65 536 cells. Each cell stores U and V as 32-bit floats, totalling ≈ 0.5 MB per double buffer.

Stability

The explicit Euler scheme is stable for Δt = 1 with standard Gray-Scott Du = 0.2097, Dv = 0.105. Larger Δt values cause oscillation blow-up.

Substeps

12 PDE steps per rAF frame give roughly 720 simulation steps/second at 60 fps — enough to see pattern convergence within seconds.

Colour LUT

Four pre-built 256-entry RGBA LUTs (Purple, Fire, Greyscale, Ocean) map V ∈ [0, 0.5] to visible colour with O(1) lookup.

Mathematical Background

Linear Stability Analysis

The uniform steady state (u*, v*) of the Gray-Scott equations is found by setting ∂U/∂t = ∂V/∂t = 0. Perturbing this steady state by δu = εexp(σt + iq·x) yields a 2×2 Jacobian whose eigenvalues σ(q) determine stability. The Turing condition is satisfied when a purely real mode q_c > 0 first becomes unstable — giving a band of growing wavenumbers and thus a preferred spatial scale.

Pearson's Phase Diagram

John Pearson's 1993 Science paper catalogued 12 dynamical regimes by sweeping the (F, k) plane at fixed Du/Dv = 2, labelling them with letters A–L. The simulator's presets sample regimes from the "self-replicating" (μ and ζ) and "oscillatory" (β and γ) zones.

Relationship to Activator-Inhibitor Systems

Though written as a chemical feed-kill system, Gray-Scott is mathematically equivalent to Gierer-Meinhardt activator-inhibitor kinetics under a variable change. Both belong to the broader class of reaction-diffusion systems whose stationary spatial solutions are called Turing patterns.

The characteristic length of Turing patterns is given by λ ≈ 2π / q_c, where q_c is the most unstable wavenumber. It scales as √(D/k), so slower reactions produce larger patterns.

Algorithm Summary

StepOperationComplexity
1Compute Laplacians ∇²U, ∇²V via 5-point stencil with periodic BCO(N²)
2Evaluate nonlinear reaction term U·V²O(N²)
3Forward-Euler update → write to secondary bufferO(N²)
4Swap primary/secondary Float32Array referencesO(1)
5Map V field to RGBA via LUT, call putImageDataO(N²)
6Scale canvas to CSS display size (pixel art mode)O(1)

Applications & Real-World Pattern Formation

DomainTuring MechanismEvidence
Developmental BiologyFinger digit spacing, limb bud patterningSheth et al. 2012 (Science)
Skin PigmentationZebrafish stripe formation via Fgf signallingKondo & Miura 2010 (Science)
Hair Follicle SpacingWnt activator + Dkk inhibitor diffusionSick et al. 2006 (Science)
Palatal RugaeRidge spacing in mammalian palate formationEconomou et al. 2012 (Nature Genetics)
ElectrochemistryCatalyst surface CO oxidation patternsErtl Nobel Lecture 2007
Material SciencePeriodic nanoscale surface structures (dewetting)Karma & Rappel 1999

Curriculum Context

LevelPrerequisiteConcepts Illustrated
A-Level / Sixth FormBasic calculus, differential equationsCoupled ODEs, concept of steady state
Undergraduate Year 2Partial differential equationsReaction-diffusion PDEs, Laplacian, stability
Undergraduate Year 3Linear algebra, Fourier analysisTuring instability, dispersion relation, wavenumber
Postgraduate / ResearchBifurcation theoryPearson phase diagram, pattern selection, weakly-nonlinear analysis

Frequently Asked Questions

Why do I see different patterns every time I reset?
The initial seeds are placed at random positions with random noise. Because the Turing instability is sensitive to initial conditions, the final pattern topology — how many spots form, where stripes nucleate — differs each run even though the global statistics (average domain size, orientation correlations) are reproducible.
What happens if I set F very low or k very high?
If F → 0 the feed is cut off and U depletes; V dies away and the system approaches the trivial (U=0, V=0) state. If k is too large, the kill rate outpaces production and V cannot sustain itself. Near these extinction boundaries you may observe travelling pulses that fade rather than stabilise.
How does this relate to the Lorenz or Ising models?
All three are canonical nonlinear dynamical systems that exhibit symmetry breaking. The Ising model breaks spin symmetry in equilibrium; the Lorenz system breaks time-reversal symmetry in a chaotic ODE; Gray-Scott breaks spatial translational symmetry in a PDE, producing Turing patterns. Together they exemplify different universality classes of self-organisation.

Related Simulations & Articles