Gray-Scott Reaction-Diffusion: Equations and Parameters
Two chemicals, two diffusion rates, and two hand-tunable numbers — from this deceptively small recipe, the Gray-Scott model (Peter Gray and Stephen Scott, 1983) produces an entire zoo of self-organising patterns: pulsating spots, branching coral, labyrinthine mazes, and travelling worms. It has become the standard entry point for anyone exploring Turing-style reaction-diffusion systems, precisely because its parameter space is small enough to explore by hand yet rich enough to never run out of surprises.
1. The Gray-Scott Equations
The model describes a continuously fed, well-mixed reactor containing two chemical species U and V, coupled by a single irreversible autocatalytic reaction:
Written as coupled reaction-diffusion partial differential equations over a 2D domain:
Each term has a direct physical reading. In the U equation: Dᵤ·∇²u spreads U by diffusion; −u·v² consumes U in the autocatalytic reaction; F·(1 − u) is a constant replenishment (feed) that pulls u back toward 1 wherever it has been depleted. In the V equation: D_v·∇²v spreads V; +u·v² is the same reaction term, now producing V; −(F+k)·v removes V both by the feed's outflow and by an additional kill rate k.
2. The Four Parameters
| Parameter | Typical range | Role |
|---|---|---|
| F (feed rate) | 0.01–0.09 | How fast U is replenished from outside; controls overall pattern density |
| k (kill rate) | 0.03–0.07 | Extra removal rate of V; together with F selects the pattern regime |
| Dᵤ | ≈ 0.16 | Diffusion coefficient of U — usually about twice D_v |
| D_v | ≈ 0.08 | Diffusion coefficient of V — must be smaller than Dᵤ for patterns to form |
The crucial asymmetry is D_v < Dᵤ: V (the "activator", which catalyses its own production) diffuses more slowly than U (the "substrate/inhibitor" supply). This combination — local self-activation with longer-range depletion of the substrate — is exactly the recipe Alan Turing identified in 1952 as capable of breaking spatial symmetry.
3. The F-k Parameter Map
Sweeping F and k while holding Dᵤ and D_v fixed reveals distinct named regions, each corresponding to a different long-term attractor of the dynamics:
| Pattern | Approx. (F, k) | Behaviour |
|---|---|---|
| Spots (pearls) | 0.035, 0.065 | Isolated, roughly circular V-rich spots on a U background |
| Mitosis | 0.028, 0.062 | Spots elongate and split in two as they grow — resembling cell division |
| Worms | 0.078, 0.061 | Elongated, mobile stripe segments that wander and occasionally split |
| Stripes / labyrinth | 0.029, 0.057 | Stationary maze-like, space-filling stripe patterns |
| Coral / chaos | 0.082, 0.059 | Branching, tree-like growth with chaotic tip-splitting |
| Waves | 0.014, 0.054 | Spiral and target waves propagating without settling |
There is no closed-form boundary between these regions — the F-k plane is explored empirically and by linear stability analysis near the boundaries. Small changes in F or k (as little as 0.002–0.005) can push the system from one qualitative behaviour to another, which is part of what makes the parameter map so rewarding to explore interactively.
Reaction-Diffusion Simulation
Sweep F and k live and watch spots become worms, mazes, and coral in real time on a GPU-accelerated grid.
4. Discretisation and Stability
On a grid with spacing h, the Laplacian is approximated with a 5-point (or the more isotropic 9-point) finite-difference stencil, and the PDE is advanced with explicit forward Euler in time:
u[i,j] ← u[i,j] + Δt·( Dᵤ·∇²u − u·v² + F·(1−u) ) v[i,j] ← v[i,j] + Δt·( D_v·∇²v + u·v² − (F+k)·v )
Explicit Euler diffusion is stable only when the timestep satisfies the Courant-Friedrichs-Lewy (CFL) condition for the diffusion term:
Exceeding this bound produces numerical oscillations that grow without limit — the grid visibly "explodes" into noise within a few dozen steps. Because of this constraint, most implementations run several sub-steps per rendered frame rather than a single large step, and GPU implementations (via WebGL ping-pong framebuffers) parallelise the stencil evaluation across every pixel simultaneously to sustain real-time frame rates even at high resolution.
5. Why Patterns Form: Turing Instability
Linearising the equations around the uniform state and applying a perturbation of wavenumber k (spatial frequency, not to be confused with the kill-rate parameter) yields a dispersion relation σ(k) describing how quickly a perturbation of that wavelength grows or decays. Diffusion-driven (Turing) instability requires the growth rate σ(k) to be positive for a band of non-zero wavenumbers, even though the k = 0 (uniform) mode is stable. The necessary condition is that the inhibiting species (here, effectively U's depletion) diffuses fast enough relative to the self-activating species V — matching the D_v < Dᵤ requirement identified in Section 2.
The wavenumber at peak growth rate sets the characteristic spacing of the resulting pattern — larger Dᵤ/D_v ratios and lower F generally produce coarser, more widely spaced features, while smaller ratios produce finer textures.
6. Real-World Analogues
Animal Coat Patterns
Reaction-diffusion systems of this general activator-inhibitor type are the leading explanation for periodic patterns on animal skins — leopard spots, zebra stripes, and pufferfish patterns. The Gray-Scott model itself is a simplified caricature; more detailed morphogen models (e.g. Meinhardt's activator-inhibitor systems) are used for quantitative comparison to biology, but Gray-Scott captures the same qualitative pattern zoo with far fewer parameters.
Chemical Oscillators
The Gray-Scott chemistry (though hypothetical as literally written) belongs to the same broad family as the Belousov-Zhabotinsky reaction: an autocatalytic step coupled to a depleting substrate. Both systems illustrate how simple, local chemical rules generate global spatial order far from equilibrium.
Numerical Methods Benchmark
Because Gray-Scott is stiff (fast reaction, slower diffusion), it is a common benchmark for testing implicit and semi-implicit PDE solvers, adaptive timestepping schemes, and GPU parallelisation strategies before applying them to more complex reaction-diffusion or fluid systems.
Frequently Asked Questions
Why does D_v need to be smaller than Dᵤ for patterns to form?
If V diffused as fast as or faster than U, any local excess of V would spread out before the autocatalytic reaction could reinforce it, and the system would simply relax back to the uniform state. With D_v smaller, a local patch of V stays concentrated long enough to consume nearby U and grow, while the depleted U region spreads outward faster than V, creating a self-reinforcing spot surrounded by a "moat" of low U. This is the short-range activation, long-range inhibition mechanism at the heart of Turing pattern formation.
Why do so many named regions (spots, worms, mazes) exist in the F-k plane?
The nonlinear reaction term u·v² means the system's local dynamics are bistable or excitable depending on exactly how F and k balance feed and removal. Near certain boundaries, spots become unstable to elongation (producing worms); near others, growing spots split rather than simply expand (mitosis); near yet others, the uniform state itself becomes fully unstable, producing space-filling stripe mazes rather than isolated features. There is no single mechanism — each named region corresponds to a distinct bifurcation of the local reaction kinetics combined with the pattern-forming instability.
Is the Gray-Scott model a realistic chemical reaction?
Not directly — it is a minimal, idealised model designed to be the simplest possible system exhibiting a rich pattern zoo, rather than a literal description of a specific chemical mixture. Its value is pedagogical and computational: because it has only two species and four parameters, it is tractable to explore exhaustively, while still producing the same qualitative phenomena (spots, stripes, travelling waves) seen in genuine reaction-diffusion chemistry and biology.
How is Gray-Scott related to the Belousov-Zhabotinsky reaction?
Both are reaction-diffusion systems built around an autocatalytic step, but BZ is a real, experimentally realised chemical reaction with a well-characterised mechanism (the Oregonator model), producing temporal oscillations and travelling spiral waves in an excitable medium. Gray-Scott is a purely mathematical model, without a specific real chemistry attached, tuned to explore the broadest possible range of stationary and quasi-stationary Turing-type patterns rather than oscillatory behaviour.
What happens if you increase the grid resolution?
Increasing resolution (more grid points per unit length) while keeping the physical domain size fixed changes the numerical grid spacing h. Because pattern wavelength is set by the physical diffusion coefficients and reaction rates, not by the grid, higher resolution reveals finer detail within the same pattern rather than changing the pattern itself — provided the timestep is rescaled according to the CFL condition (Δt ∝ h²) to maintain numerical stability.
Can Gray-Scott patterns be computed in 3D?
Yes — the same PDE extends directly to three dimensions with a 3D Laplacian stencil (7-point or higher-order). 3D Gray-Scott produces analogous structures: spots become blobs, stripes become sheets or tubes, and mitosis-like splitting occurs in three dimensions. 3D simulation is significantly more computationally expensive (an N³ grid versus N² in 2D) and is usually visualised via isosurface extraction or volumetric rendering rather than a direct colour map.
Who invented the Gray-Scott model and why?
Peter Gray and Stephen Scott, chemical engineers at the University of Leeds, introduced the model in a 1983-84 series of papers analysing a hypothetical continuously stirred tank reactor (CSTR) undergoing cubic autocatalysis. Their original interest was in chemical reactor stability and oscillation, not pattern formation per se; the model's remarkable spatial pattern-forming behaviour under diffusion was explored more thoroughly by later researchers, notably John Pearson's 1993 Science paper mapping the F-k parameter space and naming many of the pattern classes still used today.
How does Gray-Scott differ from the Fitzhugh-Nagumo model?
Fitzhugh-Nagumo is a simplified excitable-medium model (derived from the Hodgkin-Huxley neuron equations) tuned to produce travelling action-potential-like pulses and spiral waves, with a fast activator and slow recovery variable. Gray-Scott is instead tuned toward stationary and quasi-stationary Turing patterns from a chemical reaction perspective. Both are two-variable reaction-diffusion PDEs and share the same underlying Turing-instability mathematics, but they emphasise different regions of the qualitative pattern space (excitable waves versus static patterning).