Gray–Scott Reaction–Diffusion Model

• The Gray–Scott model is a reaction–diffusion system defined by cubic autocatalysis, continuous feed, and removal rates that produce spots, stripes, and chaotic patterns.
• It employs both deterministic and stochastic frameworks to analyze rich bifurcation structures including Turing, Hopf, and saddle–node instabilities.
• The model extends to fractional, nonlocal, and discrete formulations, offering insights for control, data assimilation, and ultradiscrete cellular automata studies.

The Gray–Scott model is a two-species reaction–diffusion system associated with the autocatalytic chemistry $U+2V\to 3V$, decay/removal of $U$ and $V$, and continuous feed of $U$. In a common nondimensional form, the concentrations $u(x,t)$ and $v(x,t)$ satisfy

$$\frac{\partial u}{\partial t}=D_u\,\Delta u-u\,v^2+F\,(1-u),\qquad \frac{\partial v}{\partial t}=D_v\,\Delta v+u\,v^2-(F+k)\,v,$$

and the system exhibits spots, stripes, traveling waves, pulses, self-replication, labyrinths, and spatio-temporal chaos. At the deterministic level it is a paradigmatic two-species model; at the mesoscopic level it also admits a derivation from the chemical master equation that introduces composite fields and multiplicative internal noise (Kyrychko et al., 2012, Cooper et al., 2013).

1. Canonical reaction scheme and field equations

The underlying chemistry consists of four elementary steps: $$U+2V \xrightarrow{\lambda} 3V,\qquad V\xrightarrow{\mu}P,\qquad U\xrightarrow{\nu}Q,\qquad \varnothing\xrightarrow{f}U.$$ In the Doi–Peliti formulation one introduces stochastic density fields $\phi_v(x,t)$ and $\phi_u(x,t)$, with $U$0 and $U$1 (Cooper et al., 2013).

A frequently used dimensionless PDE form is

$U$2

where $U$3 is the constant feed rate of $U$4, $U$5 is the additional removal rate of the autocatalyst $U$6, and $U$7 are diffusion coefficients. In the alternative notation used above, $U$8 denotes the feed rate and $U$9 the kill rate of $V$0. The spatially homogeneous reduction is obtained by dropping the diffusion terms, yielding a two-dimensional autonomous ODE for $V$1 (Kyrychko et al., 2012, Delgado et al., 2017).

Boundary conditions vary with the analytical or numerical setting. Homogeneous Neumann conditions are standard in several PDE and data-assimilation formulations, periodic conditions are used in pattern computations, zero Dirichlet conditions occur in some fractional formulations, and nonlocal Dirichlet or Neumann-type constraints arise when Laplacian diffusion is replaced by integral operators (Randrianasolo, 4 Oct 2025, Wang et al., 2018, Cappanera et al., 2022).

2. Homogeneous equilibria and bifurcation structure

For the local kinetics

$V$2

there are up to three homogeneous steady states: the trivial state

$V$3

and, when

$V$4

two nontrivial states

$V$5

The nontrivial equilibria collide in a saddle–node bifurcation when $V$6. A Hopf bifurcation of $V$7 occurs along a curve $V$8 starting at the Takens–Bogdanov point $V$9, while a Turing instability of $U$0 occurs for sufficiently large diffusion ratio $U$1 when

$U$2

For $U$3, the Turing-unstable region lies mostly inside the Hopf-unstable region, and in the narrow parameter band bounded above by the Turing boundary and below by the saddle–node curve the system exhibits spatio-temporal chaos (Kyrychko et al., 2012).

In the $U$4-parameterization of the spatially homogeneous ODE,

$U$5

the trivial equilibrium $U$6 is always a stable node. The nontrivial branches are organized by explicit saddle–node and Hopf curves, a Takens–Bogdanov point at

$U$7

and a Bautin point at

$U$8

Numerical continuation of the local bifurcation curves adds a homoclinic branch and a limit-point-of-cycles branch, producing a global bifurcation map with regions that contain, depending on location, a stable equilibrium, an unstable equilibrium enclosed by one or two limit cycles, or no nontrivial cycles (Delgado et al., 2017).

This bifurcation geometry explains why the Gray–Scott equations support both diffusion-driven stationary structures and oscillatory or excitable dynamics. It also clarifies that “the” Gray–Scott phase diagram is not a single instability threshold, but a superposition of saddle–node, Hopf, Turing, homoclinic, and fold-of-cycles phenomena (Kyrychko et al., 2012, Delgado et al., 2017).

3. Pattern-forming regimes and localized structures

Away from instabilities, solutions relax to homogeneous states. Crossing the Turing boundary with Hopf stability gives stationary spatial patterns, while outside the saddle–node curve one finds coherent traveling fronts and pulses that may self-replicate. Strong interaction of Hopf and Turing modes yields mixed states and spatio-temporal chaos. Representative deterministic regimes include stationary Turing patterns, traveling waves and pulses, spatio-temporal chaos, and self-replicating pulses (Kyrychko et al., 2012).

A narrow parameter region with $U$9 supports stable localized two-dimensional structures of three classes: stationary, moving, and rotating. For

$u(x,t)$0

and for each $u(x,t)$1 a tiny interval of $u(x,t)$2 roughly $u(x,t)$3 wide, stable localized patterns occur just below the Turing-instability line and below the saddle–node bifurcation. At

$u(x,t)$4

one finds stationary multi-spot molecules, two distinct moving patterns, and rotating aggregates. Reported examples include a “U”-shaped front with speed

$u(x,t)$5

a three-spot cluster with speed

$u(x,t)$6

and a four-spot rotor with one full revolution in $u(x,t)$7. The proposed stability mechanism is constructive reinforcement of concentric standing-wave “halos” surrounding each localized feature, with preferred binding distances set by positive halo maxima (Munafo, 2014).

In the singularly perturbed regime of small diffusivity $u(x,t)$8 for one component, multi-spot quasi-equilibria can be analyzed by matched asymptotic expansions. A differential-algebraic ODE system for collective coordinates $u(x,t)$9 and $v(x,t)$0 describes slow spot drift, while fast instabilities are classified into spot self-replication, spot oscillation, and spot annihilation. For the non-radial $v(x,t)$1 mode, numerical solution of the local eigenvalue problem yields the self-replication threshold

$v(x,t)$2

so that source strengths above this value trigger peanut-splitting and eventual spot duplication (Chen et al., 2010).

Constructive existence results now complement numerical evidence. A Newton–Kantorovich framework on a specifically constructed Hilbert space has been used to prove the existence of smooth localized stationary patterns in the two-dimensional stationary Gray–Scott system, including spike, ring, and four-lobed “leaf” solutions, all satisfying $v(x,t)$3-symmetry (Cadiot et al., 2024).

Pattern diversity has also been quantified directly from simulated concentration fields. A parameter sweep over

$v(x,t)$4

using Shannon entropy, Simpson diversity, a PNG-based approximation of Lempel–Ziv complexity, and expressivity identified high-complexity bands near

$v(x,t)$5

and near

$v(x,t)$6

The highest-complexity regimes are composed of wave-fragments and traveling localizations (Adamatzky, 2016).

4. Stochastic derivation, composite fields, and broken symmetry

A first-principles derivation starts from a lattice chemical master equation for the reactions and diffusion processes. After the standard Doi mapping to bosonic creation and annihilation operators, the probability state $v(x,t)$7 evolves by

$v(x,t)$8

with a many-body Hamiltonian whose feed and decay terms break the $v(x,t)$9 symmetry associated with global particle-number conservation. Probability conservation requires

$$\frac{\partial u}{\partial t}=D_u\,\Delta u-u\,v^2+F\,(1-u),\qquad \frac{\partial v}{\partial t}=D_v\,\Delta v+u\,v^2-(F+k)\,v,$$0

equivalently $$\frac{\partial u}{\partial t}=D_u\,\Delta u-u\,v^2+F\,(1-u),\qquad \frac{\partial v}{\partial t}=D_v\,\Delta v+u\,v^2-(F+k)\,v,$$1 for all $$\frac{\partial u}{\partial t}=D_u\,\Delta u-u\,v^2+F\,(1-u),\qquad \frac{\partial v}{\partial t}=D_v\,\Delta v+u\,v^2-(F+k)\,v,$$2. After a Doi shift $$\frac{\partial u}{\partial t}=D_u\,\Delta u-u\,v^2+F\,(1-u),\qquad \frac{\partial v}{\partial t}=D_v\,\Delta v+u\,v^2-(F+k)\,v,$$3, $$\frac{\partial u}{\partial t}=D_u\,\Delta u-u\,v^2+F\,(1-u),\qquad \frac{\partial v}{\partial t}=D_v\,\Delta v+u\,v^2-(F+k)\,v,$$4, the saddle-point equations of the shifted path integral recover the usual deterministic Gray–Scott PDEs (Cooper et al., 2013).

The stochastic formulation introduces composite fields that expose the internal structure of the autocatalytic vertex. In the broken-symmetry path integral one defines

$$\frac{\partial u}{\partial t}=D_u\,\Delta u-u\,v^2+F\,(1-u),\qquad \frac{\partial v}{\partial t}=D_v\,\Delta v+u\,v^2-(F+k)\,v,$$5

Here $$\frac{\partial u}{\partial t}=D_u\,\Delta u-u\,v^2+F\,(1-u),\qquad \frac{\partial v}{\partial t}=D_v\,\Delta v+u\,v^2-(F+k)\,v,$$6 represents a bound state of two $$\frac{\partial u}{\partial t}=D_u\,\Delta u-u\,v^2+F\,(1-u),\qquad \frac{\partial v}{\partial t}=D_v\,\Delta v+u\,v^2-(F+k)\,v,$$7 molecules, $$\frac{\partial u}{\partial t}=D_u\,\Delta u-u\,v^2+F\,(1-u),\qquad \frac{\partial v}{\partial t}=D_v\,\Delta v+u\,v^2-(F+k)\,v,$$8 is a three-body composite that mediates autocatalysis, and $$\frac{\partial u}{\partial t}=D_u\,\Delta u-u\,v^2+F\,(1-u),\qquad \frac{\partial v}{\partial t}=D_v\,\Delta v+u\,v^2-(F+k)\,v,$$9 represents a $$U+2V \xrightarrow{\lambda} 3V,\qquad V\xrightarrow{\mu}P,\qquad U\xrightarrow{\nu}Q,\qquad \varnothing\xrightarrow{f}U.$$0–$$U+2V \xrightarrow{\lambda} 3V,\qquad V\xrightarrow{\mu}P,\qquad U\xrightarrow{\nu}Q,\qquad \varnothing\xrightarrow{f}U.$$1 bound state that enters fluctuation-induced elastic re-scattering. Converting the action so that it is at most quadratic in the conjugate fields leads to Langevin equations with multiplicative noise for $$U+2V \xrightarrow{\lambda} 3V,\qquad V\xrightarrow{\mu}P,\qquad U\xrightarrow{\nu}Q,\qquad \varnothing\xrightarrow{f}U.$$2, $$U+2V \xrightarrow{\lambda} 3V,\qquad V\xrightarrow{\mu}P,\qquad U\xrightarrow{\nu}Q,\qquad \varnothing\xrightarrow{f}U.$$3, and $$U+2V \xrightarrow{\lambda} 3V,\qquad V\xrightarrow{\mu}P,\qquad U\xrightarrow{\nu}Q,\qquad \varnothing\xrightarrow{f}U.$$4, with nontrivial cross-correlations such as

$$U+2V \xrightarrow{\lambda} 3V,\qquad V\xrightarrow{\mu}P,\qquad U\xrightarrow{\nu}Q,\qquad \varnothing\xrightarrow{f}U.$$5

$$U+2V \xrightarrow{\lambda} 3V,\qquad V\xrightarrow{\mu}P,\qquad U\xrightarrow{\nu}Q,\qquad \varnothing\xrightarrow{f}U.$$6

where $$U+2V \xrightarrow{\lambda} 3V,\qquad V\xrightarrow{\mu}P,\qquad U\xrightarrow{\nu}Q,\qquad \varnothing\xrightarrow{f}U.$$7. Because the amplitude of $$U+2V \xrightarrow{\lambda} 3V,\qquad V\xrightarrow{\mu}P,\qquad U\xrightarrow{\nu}Q,\qquad \varnothing\xrightarrow{f}U.$$8 is itself stochastic, perturbative expansion generates higher-order $$U+2V \xrightarrow{\lambda} 3V,\qquad V\xrightarrow{\mu}P,\qquad U\xrightarrow{\nu}Q,\qquad \varnothing\xrightarrow{f}U.$$9 elastic scattering processes for $\phi_v(x,t)$0 (Cooper et al., 2013).

At the classical level, the composite interpretation can be made explicit by introducing an auxiliary field $\phi_v(x,t)$1 with

$\phi_v(x,t)$2

Together with

$\phi_v(x,t)$3

this yields a three-field reaction–diffusion system in which the heavy composite $\phi_v(x,t)$4 relaxes toward $\phi_v(x,t)$5. Under adiabatic elimination, $\phi_v(x,t)$6, and the standard Gray–Scott equations are recovered exactly. The late-time dynamics of the three-field model reproduces the same pattern formation as the two-field model for suitable composite diffusion parameters (Dawson et al., 2019).

A common simplification is therefore to regard the Gray–Scott model purely as a cubic two-field PDE. The master-equation derivation shows that this is the classical limit of a mesoscopic theory in which probability conservation, broken $\phi_v(x,t)$7 symmetry, composite intermediates, and multiplicative internal noise are structurally significant (Cooper et al., 2013).

5. Fractional, nonlocal, and reversible extensions

Anomalous diffusion is introduced by replacing the Laplacian with a fractional Laplacian of order $\phi_v(x,t)$8: $\phi_v(x,t)$9 For this fractional Gray–Scott system, continuous solutions are unique, a Crank–Nicolson time discretization combined with weighted shifted Grünwald spatial differences is unconditionally stable, and benchmark computations confirm second-order convergence: $\phi_u(x,t)$0 Simulations show that decreasing $\phi_u(x,t)$1 from $\phi_u(x,t)$2 to $\phi_u(x,t)$3 produces smaller, more densely packed spots, earlier replication, and modified wave propagation. Steady-state spot arrangements can be quantified by the radial distribution function $\phi_u(x,t)$4, whose first peak satisfies the empirical scaling law

$\phi_u(x,t)$5

with fitted $\phi_u(x,t)$6 depending on $\phi_u(x,t)$7 (Wang et al., 2018).

A mixed fractional–classical formulation replaces only the $\phi_u(x,t)$8-equation’s Laplacian by a fractional operator of order $\phi_u(x,t)$9, while retaining the classical Laplacian in the $U$00-equation. Using semigroup methods and duality estimates, one obtains unique global-in-time, componentwise nonnegative solutions with uniform-in-time bounds such as

$U$01

Numerically, decreasing $U$02 produces coarser, less regular spot-like structures, larger characteristic wavelengths, and suppression of small-scale instabilities; as $U$03, classical Fickian patterns are recovered (Alam, 28 Sep 2025).

Nonlocal diffusion replaces $U$04 by an integral operator

$U$05

or, more generally,

$U$06

For integrable kernels satisfying a mass-balance condition, the nonlocal Gray–Scott system is globally well posed in $U$07, the semiflow preserves nonnegativity, and solutions converge to the classical Gray–Scott PDE in a rigorous diffusive limit. In one-dimensional pulse computations, the boundary treatment matters strongly when kernels are wide: exponential kernels can produce mesa-shaped pulses, algebraic kernels can produce cat-ear pulses, Dirichlet constraints sharpen edges, and periodic conditions can induce spurious oscillations (Laurençot et al., 2023, Cappanera et al., 17 Apr 2025).

Thermodynamically consistent closed-system variants enlarge the chemistry to include reverse reactions and auxiliary species. A reversible four-species Gray–Scott type model derived by the energetic variational approach introduces species $U$08 and $U$09, admits a natural entropy structure, and converges to an irreversible Gray–Scott type system as backward reaction coefficients tend to zero (Liang et al., 2021). A related variational model with species $U$10 has the form

$U$11

and the classical two-species subsystem is recovered formally as $U$12. In one spatial dimension, stationary classical patterns persist as long-lived transients for small $U$13, while oscillatory and traveling-wave-like patterns have persistence time of order $U$14 (Hao et al., 2024).

6. Control, discretization, and data-driven reconstruction

The model is sufficiently structured to support explicit feedback-control design. A local time-delayed feedback control term

$U$15

has been studied in single-species, diagonal, and mixed configurations. Linearization about the nontrivial homogeneous state $U$16 yields stability boundaries in the $U$17-plane. Activator control requires $U$18 to stabilize $U$19, inhibitor control requires $U$20, diagonal control is essentially ineffective for spatio-temporal chaos or traveling waves, and mixed control unlocks stable mixed Turing–Hopf states, coarsening, and bistability between a trivial steady state and traveling waves (Kyrychko et al., 2012).

Data assimilation has been formulated through nudging with coarse cell averages. For the Gray–Scott system on $U$21 with homogeneous Neumann conditions, the reconstructed fields $U$22 satisfy

$U$23

$U$24

where $U$25 is a finite-volume interpolant on a coarse observation grid. Under conditions linking observation resolution, nudging gains, diffusion, and time step, the $U$26-error decays exponentially for both the continuous problem and a fully discrete semi-implicit finite-volume scheme. Numerical experiments in a labyrinthine regime show that observing only the $U$27-field can suffice to recover the full $U$28 dynamics (Randrianasolo, 4 Oct 2025).

Discrete and ultradiscrete realizations show that Gray–Scott phenomenology survives severe reduction of the state space. A nonstandard discrete-time, discrete-space system reproduces traveling pulses, self-replication, and steady states. After ultradiscretization, the resulting max–plus cellular automata include a regime equivalent to Elementary Cellular Automaton Rule 90,

$U$29

which generates a Sierpinski gasket from a single-site seed. In two dimensions, the ultradiscrete system also produces ring patterns, self-replication, and chaotic patchy states (Matsuya et al., 2013).

Numerically, one-dimensional occupation by self-replicating pulses has been computed by two distinct high-order approaches: sinc differential quadrature combined with third-fourth order implicit Rosenbrock time stepping, and exponential B-spline collocation with Crank–Nicolson. These methods capture regimes of propagating replication, standing pulses, and domain covering by multiple initial pulses (Korkmaz et al., 2016).

Taken together, these developments place the Gray–Scott model at the intersection of nonlinear pattern theory, stochastic reaction–diffusion field theory, anomalous and nonlocal transport, thermodynamically consistent extensions, and modern control and inference. The same cubic autocatalytic core supports deterministic bifurcation analysis, constructive existence theorems, master-equation derivations, coarse-observation reconstruction, and ultradiscrete cellular-automaton limits (Delgado et al., 2017, Cadiot et al., 2024, Cooper et al., 2013, Randrianasolo, 4 Oct 2025, Matsuya et al., 2013).