How the Leopard Got Its Spots at 60 FPS: Gray–Scott Reaction–Diffusion in WebGL

After building a multi-pass incompressible Navier–Stokes solver on ping-pong WebGL Framebuffer Objects (FBOs), a natural question arises: what happens when you strip away the velocity field and let two diffusing chemicals react directly on the grid?

Common sense says diffusion is a great equalizer – drop ink into still water and the gradients smooth out into a uniform blur. In 1952, Alan Turing proved that common sense is wrong. In his seminal paper "The Chemical Basis of Morphogenesis," Turing showed that if an autocatalytic activator diffuses more slowly than its substrate or inhibitor ($D_v \lt D_u$), diffusion itself destabilizes the uniform state and spontaneously carves stripes, spots, labyrinths, and dividing cells out of nothing.

This post builds a real-time WebGL solver for the Gray–Scott reaction–diffusion system (Pearson, 1993), executing 14 PDE substeps per animation frame (~840 chemical steps per second) on an isotropic 9-point Laplacian stencil and shading the concentration field as a 3D Blinn–Phong normal-mapped coral reef.

0. Anatomy of the Gray–Scott System

Every texel in the simulation above stores two dimensionless chemical concentrations in its red and green channels, $(u, v) \in [0, 1]^2$, governed by two irreversible reactions inside an open Continuously Stirred Tank Reactor (CSTR):

Autocatalysis ($U + 2V \to 3V$)

Two molecules of catalyst $V$ convert one molecule of substrate $U$ into a third molecule of $V$ at nonlinear rate $R(u, v) = u v^2$.

Feed Rate ($F$)

Fresh substrate $U$ is continuously replenished toward $u = 1$ at rate $F(1 - u)$, while both species are washed out at rate $F$.

Kill Rate ($k$)

Catalyst $V$ decays into an inert product $P$ ($V \to P$) at rate $k v$, giving a total removal rate $(F + k)v$.

Diffusion Ratio ($D_u / D_v \approx 2$)

With $D_u = 0.2097$ and $D_v = 0.1050$, substrate $U$ diffuses twice as fast as catalyst $V$, creating local lateral depletion zones around each growing $V$ peak.

Ping-Pong Substepping

Because each step is a single fullscreen pass with no Poisson solve, the GPU runs 14 ping-pong substeps per frame at 60 fps.

1. The Gray–Scott PDEs and Diffusion-Driven Instability

Combining Fickian diffusion ($D \nabla^2$) with the cubic autocatalytic reaction $u v^2$, substrate replenishment $F(1 - u)$, and catalyst decay $(F + k)v$ yields the coupled Gray–Scott partial differential equations:

Why does $D_u \gt D_v$ create intricate self-organizing structures instead of a uniform soup? Consider a tiny island where $v$ is locally high:

1. Local self-amplification: Inside the island, the cubic reaction $+u v^2$ rapidly converts $u$ into more $v$, trying to make the peak taller.
2. Fast lateral refueling ($D_u \gg D_v$): As $u$ is depleted at the center of the island, fast-diffusing substrate $U$ rushes in from the surrounding sea ($u \approx 1$) to feed the outer rim of the island, while slow-diffusing $V$ remains trapped in a sharp ridge.
3. Core starvation & bifurcation: Once the island grows wider than the diffusion length $\ell_u \sim \sqrt{D_u / F}$, fresh substrate $U$ can no longer reach the center before being consumed by the rim. The center starves and collapses back to $v \approx 0$, splitting a single spot into a dividing dumbbell (Mitosis) or a hollow branching ring (Coral Labyrinth)!

2. The Isotropic 9-Point Discrete Laplacian

If you discretize $\nabla^2 u$ with the standard 5-point cross stencil ($N + S + E + W - 4C$), its leading truncation error $\frac{h^2}{12}(\partial_x^4 + \partial_y^4)$ is strongly anisotropic: waves travel at different speeds along the grid axes ($0^\circ, 90^\circ$) than along the diagonals ($45^\circ$), turning circular Turing spots into diamond-shaped squares.

Blending the 5-point cardinal cross with the 4 diagonal corners ($NE, NW, SE, SW$) in a $4:1$ ratio eliminates the $\mathcal{O}(h^2)$ directional anisotropy, making the discrete Laplacian rotationally invariant to fourth order:

void main() {
  vec2 C  = texture2D(uState, vUv).rg;
  vec2 N  = texture2D(uState, vUv + vec2( 0.0,       uTexel.y)).rg;
  vec2 S  = texture2D(uState, vUv + vec2( 0.0,      -uTexel.y)).rg;
  vec2 E  = texture2D(uState, vUv + vec2( uTexel.x,  0.0     )).rg;
  vec2 W  = texture2D(uState, vUv + vec2(-uTexel.x,  0.0     )).rg;
  vec2 NE = texture2D(uState, vUv + vec2( uTexel.x,  uTexel.y)).rg;
  vec2 NW = texture2D(uState, vUv + vec2(-uTexel.x,  uTexel.y)).rg;
  vec2 SE = texture2D(uState, vUv + vec2( uTexel.x, -uTexel.y)).rg;
  vec2 SW = texture2D(uState, vUv + vec2(-uTexel.x, -uTexel.y)).rg;

  // Isotropic 9-point Laplacian (0.20 cardinal, 0.05 diagonal, -1.0 center)
  vec2 lap = 0.20 * (N + S + E + W) + 0.05 * (NE + NW + SE + SW) - C;

  float u   = C.r;
  float v   = C.g;
  float uvv = u * v * v;

  // Optional spatial (F, k) parameter map across the canvas
  float F = (uSpatialMap > 0.5) ? mix(0.018, 0.058, vUv.y) : uF;
  float k = (uSpatialMap > 0.5) ? mix(0.048, 0.065, vUv.x) : uK;

  float du = uDu * lap.r - uvv + F * (1.0 - u);
  float dv = uDv * lap.g + uvv - (F + k) * v;

  u = clamp(u + du * uDt, 0.0, 1.0);
  v = clamp(v + dv * uDt, 0.0, 1.0);

  gl_FragColor = vec4(u, v, 0.0, 1.0);
}

3. Turning a 2D Chemical Field into 3D Coral Relief

A raw grayscale plot of $v(\vec x)$ looks flat. Instead, the display shader treats the catalyst concentration $v(x, y)$ as a 3D heightfield $z = \alpha\,v(x, y)$ and reconstructs its unit surface normal $\hat n$ per pixel using central finite differences:

4. The Multi-Substep Ping-Pong Pipeline

Because the Gray–Scott step requires only a single 9-tap texture read, each animation frame loops the reaction–diffusion shader 14 times across a single RGBA16F double FBO before computing normals and lighting in the display pass:

5. A Tour of Pearson's $(F, k)$ Morphogenesis Zoo

All four preset buttons in the simulation run the exact same 30-line fragment shader – only the two scalar uniforms $F$ and $k$ change:

Click Mitosis Cells in the controls above and paint a single dot on an empty patch of the canvas – watching one droplet of $V$ swell, pinch its waist, and divide exponentially into a colony of hundreds of cells is the purest demonstration of how two lines of partial differential equations can encode biological form.