Reaction-Diffusion: Turing Patterns on the GPU

The Gray–Scott Equations

The Gray–Scott model tracks two chemical species, U and V, across a 2D grid. U is fed continuously into the system; V is consumed. Their interaction is the reaction U + 2V → 3V:

Two parameters, f and k, control everything. Tiny shifts in their values produce radically different pattern types: spots, stripes, mazes, holes, coral, fingerprints.

Why the CPU Is the Wrong Place for This

For a 512×512 grid, that's 262,144 cells. Each cell needs the values at its four neighbours (the Laplacian), then two floating-point multiplications, and two additions. Running this in JavaScript at 60 FPS requires ~31 million operations per second — achievable, but just barely, and with no headroom for a larger grid or any post-processing.

The GPU solution is obvious: the Gray–Scott update is a classic embarrassingly parallel operation — every cell is independent. A fragment shader that reads from one texture and writes to another runs the same update on all 262,144 cells simultaneously.

Ping-Pong Framebuffers

You cannot read and write the same texture in a single render pass (undefined behaviour in WebGL). The solution is ping-ponging: maintain two render targets and alternate which one you read from and write to each frame.

// Setup: two identical render targets
const size = 512;
const opts = {
  type: THREE.FloatType,        // float32 for precision
  format: THREE.RGBAFormat,     // R = U concentration, G = V concentration
  minFilter: THREE.LinearFilter,
  magFilter: THREE.LinearFilter
};
let rtA = new THREE.WebGLRenderTarget(size, size, opts);
let rtB = new THREE.WebGLRenderTarget(size, size, opts);

// Each frame: swap
function tick() {
  // Step 1: render Gray-Scott update (read A, write B)
  uniforms.tState.value = rtA.texture;
  renderer.setRenderTarget(rtB);
  renderer.render(quadScene, orthoCamera);

  // Step 2: display result on screen
  displayUniforms.tState.value = rtB.texture;
  renderer.setRenderTarget(null);
  renderer.render(displayScene, orthoCamera);

  // Swap
  [rtA, rtB] = [rtB, rtA];
  requestAnimationFrame(tick);
}

The Update Fragment Shader

The GLSL update shader samples the current state texture and applies the Gray–Scott equations. The discrete Laplacian uses a 5-point stencil — the classic finite-difference approximation.

uniform sampler2D tState;   // current RG texture (R=U, G=V)
uniform vec2 uTexelSize;    // 1/width, 1/height
uniform float uF;           // feed rate
uniform float uK;           // kill rate
uniform float uDu;          // diffusion U
uniform float uDv;          // diffusion V
uniform float uDt;          // time step per frame

void main() {
  vec2 uv = vUv;
  vec2 ts = uTexelSize;

  // Sample current cell and 4 neighbours
  vec2 curr  = texture2D(tState, uv).rg;
  vec2 top   = texture2D(tState, uv + vec2( 0, ts.y)).rg;
  vec2 bot   = texture2D(tState, uv + vec2( 0,-ts.y)).rg;
  vec2 left  = texture2D(tState, uv + vec2(-ts.x, 0)).rg;
  vec2 right = texture2D(tState, uv + vec2( ts.x, 0)).rg;

  // 5-point discrete Laplacian
  vec2 lap = top + bot + left + right - 4.0 * curr;

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

  // Gray-Scott reaction-diffusion equations
  float du = uDu * lap.r - uvv + uF * (1.0 - u);
  float dv = uDv * lap.g + uvv - (uF + uK) * v;

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

The Parameter Zoo

The most fascinating aspect of Gray–Scott is the phase diagram. Keeping D_u = 0.21 and D_v = 0.105 fixed, here are the patterns that emerge from varying just F and k:

Seeding: Where Patterns Come From

The initial condition matters enormously. Starting from a uniform state produces no pattern — it's a stable equilibrium. To kick the system into pattern-forming mode, I seed small "drops" of high-V concentration in a field of high-U:

// Seed initial texture (on CPU, uploaded once)
const data = new Float32Array(size * size * 4);
for (let i = 0; i < size * size; i++) {
  data[i * 4 + 0] = 1.0; // U = 1 everywhere
  data[i * 4 + 1] = 0.0; // V = 0 everywhere
}
// Drop 20 random seeds of V
for (let s = 0; s < 20; s++) {
  const cx = Math.floor(Math.random() * size);
  const cy = Math.floor(Math.random() * size);
  for (let dx = -4; dx <= 4; dx++) {
    for (let dy = -4; dy <= 4; dy++) {
      const idx = ((cy + dy) * size + (cx + dx)) * 4;
      data[idx + 0] = 0.5; // reduce U
      data[idx + 1] = 0.25; // inject V
    }
  }
}
const initTexture = new THREE.DataTexture(data, size, size, THREE.RGBAFormat, THREE.FloatType);

Interactive Painting

The most satisfying feature: click-and-drag to inject V directly onto the texture at runtime. A mouse event handler converts screen coordinates to UV space, then renders a small "brush" of high-V into the current framebuffer. The patterns you paint grow and evolve organically, merging with the existing structure.

This required care around coordinate transforms — the quad is drawn in normalised device coordinates but the mouse position is in CSS pixels. One getBoundingClientRect() call and a viewport-to-UV conversion solved it cleanly.

What Surprised Me Most

Two things. First, how slow the patterns are on biological timescales — the simulation runs thousands of Euler steps per second, yet real organism patterning happens over hours to days of embryonic development. The stunningly complex structure you see crystallise in seconds in the simulation takes an organism a full developmental cycle to produce.

Second, that the same equations describe phenomena across seven orders of magnitude of scale: from sub-millimetre cell signalling patches to dune fields kilometres wide. The Laplacian + reaction nonlinearity combination appears to be a universal pattern-forming machine that nature discovered early and never stopped using.