Reaction Diffusion: How to Create Organic Turing Patterns With Code

The spots on a leopard, the stripes on a zebrafish, the branching pattern of coral, the whorls on your fingertips. None of these are painted on. They emerge from chemistry—two substances diffusing across a surface and reacting with each other. Alan Turing figured this out in 1952, in one of the last papers he published before his death. He called them “chemical morphogens.” We call the result reaction-diffusion.

The math is surprisingly simple. Two chemicals, A and B. A feeds itself and produces B. B kills A. Both spread outward at different rates. That asymmetry—the difference in diffusion speed—is what creates pattern. When B diffuses faster than A, stable spots and stripes appear from nothing. From uniform gray soup to leopard print in a few thousand iterations.

This guide builds eight reaction-diffusion systems from scratch. Every example runs in a single HTML file, no libraries, no dependencies. Each one produces patterns you could stare at for hours.

1. Basic Gray-Scott model — the foundation

The Gray-Scott model is the most widely used reaction-diffusion system in computer graphics. Two chemicals, U and V, interact on a 2D grid. U is the “food” chemical that gets replenished from outside. V is the “catalyst” that consumes U and reproduces. The equations:

• U′ = D_u ∇²U − UV² + f(1 − U)
• V′ = D_v ∇²V + UV² − (f + k)V

D_u and D_v are diffusion rates, f is the feed rate (how fast U is replenished), and k is the kill rate (how fast V decays). Different f and k values produce wildly different patterns—spots, stripes, waves, or chaos.

const canvas = document.createElement('canvas');
const W = 256, H = 256;
canvas.width = W; canvas.height = H;
canvas.style.width = '512px'; canvas.style.height = '512px';
canvas.style.imageRendering = 'pixelated';
document.body.appendChild(canvas);
document.body.style.margin = '0';
document.body.style.background = '#000';
document.body.style.display = 'flex';
document.body.style.justifyContent = 'center';
document.body.style.alignItems = 'center';
document.body.style.height = '100vh';
const ctx = canvas.getContext('2d');
const img = ctx.createImageData(W, H);

// Two grids: current and next
let U = new Float32Array(W * H).fill(1);
let V = new Float32Array(W * H).fill(0);
let nextU = new Float32Array(W * H);
let nextV = new Float32Array(W * H);

// Seed: small square of V in the center
for (let y = H / 2 - 10; y < H / 2 + 10; y++) {
  for (let x = W / 2 - 10; x < W / 2 + 10; x++) {
    const i = y * W + x;
    V[i] = 1;
    U[i] = 0.5;
  }
}
// Add random noise to break symmetry
for (let i = 0; i < W * H; i++) {
  if (Math.random() < 0.01) { V[i] = Math.random(); U[i] = Math.random(); }
}

const Du = 0.2097, Dv = 0.105;
const f = 0.037, k = 0.06;

function laplacian(arr, x, y) {
  const i = y * W + x;
  const l = x > 0 ? i - 1 : i + W - 1;
  const r = x < W - 1 ? i + 1 : i - W + 1;
  const u = y > 0 ? i - W : i + (H - 1) * W;
  const d = y < H - 1 ? i + W : i - (H - 1) * W;
  return arr[l] + arr[r] + arr[u] + arr[d] - 4 * arr[i];
}

function step() {
  for (let y = 0; y < H; y++) {
    for (let x = 0; x < W; x++) {
      const i = y * W + x;
      const u = U[i], v = V[i];
      const uvv = u * v * v;
      nextU[i] = u + Du * laplacian(U, x, y) - uvv + f * (1 - u);
      nextV[i] = v + Dv * laplacian(V, x, y) + uvv - (f + k) * v;
      nextU[i] = Math.max(0, Math.min(1, nextU[i]));
      nextV[i] = Math.max(0, Math.min(1, nextV[i]));
    }
  }
  [U, nextU] = [nextU, U];
  [V, nextV] = [nextV, V];
}

function draw() {
  for (let i = 0; i < W * H; i++) {
    const val = Math.floor((1 - V[i]) * 255);
    const p = i * 4;
    img.data[p] = val;
    img.data[p + 1] = val;
    img.data[p + 2] = Math.min(255, val + 30);
    img.data[p + 3] = 255;
  }
  ctx.putImageData(img, 0, 0);
}

function animate() {
  for (let s = 0; s < 8; s++) step();
  draw();
  requestAnimationFrame(animate);
}
animate();

After a few hundred frames, spots and worm-like stripes emerge from the seeded center. The pattern never repeats and never reaches a truly stable state—it keeps slowly shifting. Those f and k values (0.037 and 0.06) sit in the “mitosis” region of the Gray-Scott parameter space. Change them and you get entirely different morphologies.

2. Turing spots and stripes — parameter space explorer

The Gray-Scott model has a rich parameter space. Small changes in f and k produce dramatically different patterns: spots that split and multiply, labyrinthine stripes, pulsing solitons, or uniform extinction. This example divides the canvas into a 3×3 grid, each cell running its own f/k pair simultaneously.

const canvas = document.createElement('canvas');
const cellSize = 128, cols = 3, rows = 3;
const W = cellSize * cols, H = cellSize * rows;
canvas.width = W; canvas.height = H;
canvas.style.width = (W * 2) + 'px'; canvas.style.height = (H * 2) + 'px';
canvas.style.imageRendering = 'pixelated';
document.body.appendChild(canvas);
document.body.style.margin = '0';
document.body.style.background = '#111';
document.body.style.display = 'flex';
document.body.style.justifyContent = 'center';
document.body.style.alignItems = 'center';
document.body.style.height = '100vh';
const ctx = canvas.getContext('2d');
const img = ctx.createImageData(W, H);

// Nine parameter sets spanning the Gray-Scott zoo
const params = [
  { f: 0.026, k: 0.051, label: 'Spots' },
  { f: 0.030, k: 0.055, label: 'Worms' },
  { f: 0.037, k: 0.060, label: 'Mitosis' },
  { f: 0.039, k: 0.058, label: 'Coral' },
  { f: 0.050, k: 0.065, label: 'Chaos' },
  { f: 0.014, k: 0.054, label: 'Moving spots' },
  { f: 0.018, k: 0.051, label: 'Pulsing' },
  { f: 0.029, k: 0.057, label: 'Labyrinth' },
  { f: 0.042, k: 0.063, label: 'Holes' },
];

const Du = 0.2097, Dv = 0.105;
let U = new Float32Array(W * H).fill(1);
let V = new Float32Array(W * H).fill(0);
let nU = new Float32Array(W * H);
let nV = new Float32Array(W * H);

// Seed each cell with a central blob
for (let ci = 0; ci < 9; ci++) {
  const cx = (ci % cols) * cellSize + cellSize / 2;
  const cy = Math.floor(ci / cols) * cellSize + cellSize / 2;
  for (let dy = -6; dy <= 6; dy++) {
    for (let dx = -6; dx <= 6; dx++) {
      const i = (cy + dy) * W + (cx + dx);
      V[i] = 1; U[i] = 0.5;
    }
  }
}

function getParam(x, y) {
  const ci = Math.floor(x / cellSize) + Math.floor(y / cellSize) * cols;
  return params[Math.min(ci, 8)];
}

function lap(arr, x, y) {
  const i = y * W + x;
  const l = x > 0 ? arr[i - 1] : arr[i];
  const r = x < W - 1 ? arr[i + 1] : arr[i];
  const u = y > 0 ? arr[i - W] : arr[i];
  const d = y < H - 1 ? arr[i + W] : arr[i];
  return l + r + u + d - 4 * arr[i];
}

function step() {
  for (let y = 0; y < H; y++) {
    for (let x = 0; x < W; x++) {
      const i = y * W + x;
      const p = getParam(x, y);
      const u = U[i], v = V[i], uvv = u * v * v;
      nU[i] = Math.max(0, Math.min(1, u + Du * lap(U, x, y) - uvv + p.f * (1 - u)));
      nV[i] = Math.max(0, Math.min(1, v + Dv * lap(V, x, y) + uvv - (p.f + p.k) * v));
    }
  }
  [U, nU] = [nU, U]; [V, nV] = [nV, V];
}

function draw() {
  for (let i = 0; i < W * H; i++) {
    const v = V[i];
    const r = Math.floor(v * 80);
    const g = Math.floor(v * 180);
    const b = Math.floor((1 - v) * 220);
    const p = i * 4;
    img.data[p] = r; img.data[p + 1] = g; img.data[p + 2] = b; img.data[p + 3] = 255;
  }
  ctx.putImageData(img, 0, 0);
  // Draw labels
  ctx.font = '11px monospace';
  ctx.fillStyle = '#fff';
  ctx.textAlign = 'center';
  for (let ci = 0; ci < 9; ci++) {
    const cx = (ci % cols) * cellSize + cellSize / 2;
    const cy = Math.floor(ci / cols) * cellSize + 14;
    ctx.fillText(params[ci].label, cx, cy);
  }
}

function animate() {
  for (let s = 0; s < 6; s++) step();
  draw();
  requestAnimationFrame(animate);
}
animate();

Run this and watch nine different universes evolve simultaneously. “Mitosis” produces cells that split like bacteria. “Labyrinth” grows maze-like corridors. “Coral” branches outward like a reef. All from the same two equations with slightly different constants. Karl Sims mapped this parameter space extensively in the 1990s—it remains one of the richest explorations of emergent pattern from minimal rules.

3. Coral growth — branching reaction-diffusion

Real coral grows through a reaction-diffusion process combined with nutrient flow. This simulation emphasizes the branching behavior by using an asymmetric seed and parameters tuned for finger-like protrusions. The result looks remarkably like actual Acropora coral.

const canvas = document.createElement('canvas');
const W = 300, H = 300;
canvas.width = W; canvas.height = H;
canvas.style.width = '600px'; canvas.style.height = '600px';
canvas.style.imageRendering = 'pixelated';
document.body.appendChild(canvas);
document.body.style.margin = '0';
document.body.style.background = '#001828';
document.body.style.display = 'flex';
document.body.style.justifyContent = 'center';
document.body.style.alignItems = 'center';
document.body.style.height = '100vh';
const ctx = canvas.getContext('2d');
const img = ctx.createImageData(W, H);

let U = new Float32Array(W * H).fill(1);
let V = new Float32Array(W * H).fill(0);
let nU = new Float32Array(W * H);
let nV = new Float32Array(W * H);

// Seed: scattered small blobs along the bottom
for (let s = 0; s < 5; s++) {
  const sx = W * 0.2 + Math.random() * W * 0.6;
  const sy = H * 0.7 + Math.random() * H * 0.15;
  for (let dy = -4; dy <= 4; dy++) {
    for (let dx = -4; dx <= 4; dx++) {
      const x = Math.floor(sx + dx), y = Math.floor(sy + dy);
      if (x >= 0 && x < W && y >= 0 && y < H) {
        const i = y * W + x;
        V[i] = 1; U[i] = 0.25;
      }
    }
  }
}

const Du = 0.16, Dv = 0.08;
const f = 0.035, k = 0.065;

function lap(arr, x, y) {
  const i = y * W + x;
  // Weighted Laplacian (includes diagonals)
  let sum = -arr[i] * 4;
  if (x > 0) sum += arr[i - 1];
  else sum += arr[i];
  if (x < W - 1) sum += arr[i + 1];
  else sum += arr[i];
  if (y > 0) sum += arr[i - W];
  else sum += arr[i];
  if (y < H - 1) sum += arr[i + W];
  else sum += arr[i];
  return sum;
}

function step() {
  for (let y = 0; y < H; y++) {
    for (let x = 0; x < W; x++) {
      const i = y * W + x;
      const u = U[i], v = V[i], uvv = u * v * v;
      // Slight upward bias to simulate nutrient flow
      const bias = y > 0 ? V[(y - 1) * W + x] * 0.001 : 0;
      nU[i] = Math.max(0, Math.min(1, u + Du * lap(U, x, y) - uvv + f * (1 - u)));
      nV[i] = Math.max(0, Math.min(1, v + Dv * lap(V, x, y) + uvv - (f + k) * v + bias));
    }
  }
  [U, nU] = [nU, U]; [V, nV] = [nV, V];
}

function draw() {
  for (let i = 0; i < W * H; i++) {
    const v = V[i], u = U[i];
    // Ocean-coral palette
    const r = Math.floor(v * 255 * 0.9 + (1 - u) * 40);
    const g = Math.floor(v * 120 + u * 60);
    const b = Math.floor((1 - v) * 180 + 40);
    const p = i * 4;
    img.data[p] = Math.min(255, r);
    img.data[p + 1] = Math.min(255, g);
    img.data[p + 2] = Math.min(255, b);
    img.data[p + 3] = 255;
  }
  ctx.putImageData(img, 0, 0);
}

function animate() {
  for (let s = 0; s < 10; s++) step();
  draw();
  requestAnimationFrame(animate);
}
animate();

The upward bias in the V equation simulates how real coral grows toward nutrient-rich water. The scattered seed points along the bottom compete for space, producing the characteristic branching pattern where fingers grow away from each other toward open water. Biologists have used these exact equations to model actual coral species—the resemblance is not coincidental.

4. Fingerprint pattern — anisotropic diffusion

Fingerprints form through a buckling process in fetal skin, but the resulting ridge pattern closely matches reaction-diffusion with directional bias. This example adds anisotropic diffusion—the chemicals spread faster along one direction than another, and that preferred direction slowly rotates across the surface. The result: concentric loops, whorls, and arches that look startlingly like real fingerprints.

const canvas = document.createElement('canvas');
const W = 200, H = 280;
canvas.width = W; canvas.height = H;
canvas.style.width = (W * 2.5) + 'px'; canvas.style.height = (H * 2.5) + 'px';
canvas.style.imageRendering = 'pixelated';
canvas.style.borderRadius = '40% 40% 35% 35%';
document.body.appendChild(canvas);
document.body.style.margin = '0';
document.body.style.background = '#1a1a2e';
document.body.style.display = 'flex';
document.body.style.justifyContent = 'center';
document.body.style.alignItems = 'center';
document.body.style.height = '100vh';
const ctx = canvas.getContext('2d');
const img = ctx.createImageData(W, H);

let U = new Float32Array(W * H).fill(1);
let V = new Float32Array(W * H).fill(0);
let nU = new Float32Array(W * H);
let nV = new Float32Array(W * H);

// Direction field: swirling pattern centered on the finger
const angle = new Float32Array(W * H);
const coreX = W * 0.5, coreY = H * 0.38;
for (let y = 0; y < H; y++) {
  for (let x = 0; x < W; x++) {
    const dx = x - coreX, dy = y - coreY;
    angle[y * W + x] = Math.atan2(dy, dx) + Math.PI / 2;
  }
}

// Seed: small random patches
for (let s = 0; s < 20; s++) {
  const sx = Math.floor(Math.random() * W);
  const sy = Math.floor(Math.random() * H);
  for (let dy = -3; dy <= 3; dy++) {
    for (let dx = -3; dx <= 3; dx++) {
      const nx = sx + dx, ny = sy + dy;
      if (nx >= 0 && nx < W && ny >= 0 && ny < H) {
        const i = ny * W + nx;
        V[i] = 1; U[i] = 0.5;
      }
    }
  }
}

const f = 0.037, k = 0.06;

function step() {
  for (let y = 1; y < H - 1; y++) {
    for (let x = 1; x < W - 1; x++) {
      const i = y * W + x;
      const a = angle[i];
      const cosA = Math.cos(a), sinA = Math.sin(a);
      // Anisotropic Laplacian: stronger diffusion along ridge direction
      const along = 0.28, across = 0.08;
      const Du = along, Dv = across * 0.5;
      const l = U[i - 1], r = U[i + 1], u = U[i - W], d = U[i + W];
      const lapU = (l + r) * (0.5 + cosA * cosA * 0.3) +
                   (u + d) * (0.5 + sinA * sinA * 0.3) - 4 * U[i] * 0.65;
      const lV = V[i - 1], rV = V[i + 1], uV = V[i - W], dV = V[i + W];
      const lapV = (lV + rV) * (0.5 + cosA * cosA * 0.3) +
                   (uV + dV) * (0.5 + sinA * sinA * 0.3) - 4 * V[i] * 0.65;
      const uv = U[i], vv = V[i], uvv = uv * vv * vv;
      nU[i] = Math.max(0, Math.min(1, uv + Du * lapU - uvv + f * (1 - uv)));
      nV[i] = Math.max(0, Math.min(1, vv + Dv * lapV + uvv - (f + k) * vv));
    }
  }
  [U, nU] = [nU, U]; [V, nV] = [nV, V];
}

function draw() {
  for (let i = 0; i < W * H; i++) {
    const v = V[i];
    // Fingerprint ink palette: dark ridges on light skin
    const ridge = v > 0.25 ? 1 : 0;
    const val = ridge ? 40 : 210;
    const p = i * 4;
    img.data[p] = val + 15; img.data[p + 1] = val + 5; img.data[p + 2] = val;
    img.data[p + 3] = 255;
  }
  ctx.putImageData(img, 0, 0);
}

let frame = 0;
function animate() {
  for (let s = 0; s < 8; s++) step();
  frame++;
  if (frame % 3 === 0) draw();
  if (frame < 3000) requestAnimationFrame(animate);
  else draw();
}
animate();

The direction field creates the whorl: all ridge lines curve around a central core, just like a real fingerprint’s “type lines” converge on the delta. Forensic scientists classify prints as loops, whorls, and arches—you can produce all three by changing the direction field. Move the core off-center for a loop. Add two cores for a double-loop whorl. Use parallel lines for an arch.

5. Mitosis animation — cells that split

In the Gray-Scott parameter space, there is a narrow region where spots grow, elongate, and split in two. Then each daughter spot does the same. It looks exactly like cell division. This example starts with a single spot and watches it multiply into a colony.

const canvas = document.createElement('canvas');
const W = 256, H = 256;
canvas.width = W; canvas.height = H;
canvas.style.width = '512px'; canvas.style.height = '512px';
canvas.style.imageRendering = 'pixelated';
document.body.appendChild(canvas);
document.body.style.margin = '0';
document.body.style.background = '#0a0a1a';
document.body.style.display = 'flex';
document.body.style.justifyContent = 'center';
document.body.style.alignItems = 'center';
document.body.style.height = '100vh';
const ctx = canvas.getContext('2d');
const img = ctx.createImageData(W, H);

let U = new Float32Array(W * H).fill(1);
let V = new Float32Array(W * H).fill(0);
let nU = new Float32Array(W * H);
let nV = new Float32Array(W * H);

// Single seed spot in center
for (let dy = -5; dy <= 5; dy++) {
  for (let dx = -5; dx <= 5; dx++) {
    if (dx * dx + dy * dy <= 25) {
      const i = (H / 2 + dy) * W + (W / 2 + dx);
      V[i] = 1; U[i] = 0;
    }
  }
}

// Mitosis parameters
const Du = 0.2097, Dv = 0.105;
const f = 0.0367, k = 0.0649;

function lap(arr, x, y) {
  const i = y * W + x;
  const l = x > 0 ? arr[i - 1] : arr[i + W - 1];
  const r = x < W - 1 ? arr[i + 1] : arr[i - W + 1];
  const u = y > 0 ? arr[i - W] : arr[i + (H - 1) * W];
  const d = y < H - 1 ? arr[i + W] : arr[i - (H - 1) * W];
  return l + r + u + d - 4 * arr[i];
}

function step() {
  for (let y = 0; y < H; y++) {
    for (let x = 0; x < W; x++) {
      const i = y * W + x;
      const u = U[i], v = V[i], uvv = u * v * v;
      nU[i] = Math.max(0, Math.min(1, u + Du * lap(U, x, y) - uvv + f * (1 - u)));
      nV[i] = Math.max(0, Math.min(1, v + Dv * lap(V, x, y) + uvv - (f + k) * v));
    }
  }
  [U, nU] = [nU, U]; [V, nV] = [nV, V];
}

let frame = 0;
function draw() {
  for (let i = 0; i < W * H; i++) {
    const v = V[i];
    // Cell-like palette: green organelles on dark background
    const r = Math.floor(v * 40);
    const g = Math.floor(v * 255);
    const b = Math.floor(v * 100);
    const p = i * 4;
    img.data[p] = r; img.data[p + 1] = g; img.data[p + 2] = b; img.data[p + 3] = 255;
  }
  ctx.putImageData(img, 0, 0);
  ctx.font = '12px monospace';
  ctx.fillStyle = '#4f8';
  ctx.fillText('Frame: ' + frame, 8, 18);
}

function animate() {
  for (let s = 0; s < 10; s++) { step(); frame++; }
  draw();
  requestAnimationFrame(animate);
}
animate();

Watch the single green dot stretch into an oval, pinch in the middle, and split into two. Each half grows, stretches, splits again. After a few minutes you have dozens of spots bouncing off each other in a packed colony. The parameters (f=0.0367, k=0.0649) sit at the exact edge of the mitosis region. Nudge k up by 0.001 and the spots stop splitting. Nudge it down and they dissolve into stripes. Biology works the same way—cell division is triggered by crossing a threshold.

6. Chemical waves — target patterns and spirals

The Fitzhugh-Nagumo model produces excitable waves rather than static patterns. Drop a stimulus anywhere and a ring of activation expands outward, followed by a refractory zone where the medium cannot be re-excited for a while. The result: expanding target patterns, collision annihilation, and spiral waves.

const canvas = document.createElement('canvas');
const W = 300, H = 300;
canvas.width = W; canvas.height = H;
canvas.style.width = '600px'; canvas.style.height = '600px';
canvas.style.imageRendering = 'pixelated';
document.body.appendChild(canvas);
document.body.style.margin = '0';
document.body.style.background = '#000';
document.body.style.display = 'flex';
document.body.style.justifyContent = 'center';
document.body.style.alignItems = 'center';
document.body.style.height = '100vh';
const ctx = canvas.getContext('2d');
const img = ctx.createImageData(W, H);

// u = activator, v = inhibitor (recovery variable)
let u = new Float32Array(W * H).fill(0);
let v = new Float32Array(W * H).fill(0);
let nu = new Float32Array(W * H);
let nv = new Float32Array(W * H);

const Du = 1.0, a = 0.75, b = 0.01, epsilon = 0.02, dt = 0.4;

function lap(arr, x, y) {
  const i = y * W + x;
  let s = -4 * arr[i];
  s += x > 0 ? arr[i - 1] : arr[i];
  s += x < W - 1 ? arr[i + 1] : arr[i];
  s += y > 0 ? arr[i - W] : arr[i];
  s += y < H - 1 ? arr[i + W] : arr[i];
  return s;
}

function step() {
  for (let y = 0; y < H; y++) {
    for (let x = 0; x < W; x++) {
      const i = y * W + x;
      const ui = u[i], vi = v[i];
      // Cubic reaction term (bistable)
      const reaction = ui * (1 - ui) * (ui - a) - vi;
      nu[i] = ui + dt * (Du * lap(u, x, y) + reaction);
      nv[i] = vi + dt * epsilon * (ui - b * vi);
      nu[i] = Math.max(-0.5, Math.min(1.5, nu[i]));
      nv[i] = Math.max(-0.5, Math.min(1.5, nv[i]));
    }
  }
  [u, nu] = [nu, u]; [v, nv] = [nv, v];
}

// Click to drop a stimulus
canvas.addEventListener('click', (e) => {
  const rect = canvas.getBoundingClientRect();
  const mx = Math.floor((e.clientX - rect.left) / rect.width * W);
  const my = Math.floor((e.clientY - rect.top) / rect.height * H);
  for (let dy = -8; dy <= 8; dy++) {
    for (let dx = -8; dx <= 8; dx++) {
      if (dx * dx + dy * dy > 64) continue;
      const x = mx + dx, y = my + dy;
      if (x >= 0 && x < W && y >= 0 && y < H) {
        u[y * W + x] = 1;
      }
    }
  }
});

// Initial stimulus
for (let dy = -10; dy <= 10; dy++) {
  for (let dx = -10; dx <= 10; dx++) {
    if (dx * dx + dy * dy <= 100) {
      u[(H / 2 + dy) * W + (W / 2 + dx)] = 1;
    }
  }
}

function draw() {
  for (let i = 0; i < W * H; i++) {
    const ui = Math.max(0, Math.min(1, u[i]));
    const vi = Math.max(0, Math.min(1, v[i] * 4));
    const p = i * 4;
    img.data[p] = Math.floor(ui * 255);
    img.data[p + 1] = Math.floor(ui * 100 + vi * 80);
    img.data[p + 2] = Math.floor(vi * 255);
    img.data[p + 3] = 255;
  }
  ctx.putImageData(img, 0, 0);
  ctx.font = '12px monospace';
  ctx.fillStyle = '#fff';
  ctx.fillText('Click to add waves', 8, 18);
}

function animate() {
  for (let s = 0; s < 4; s++) step();
  draw();
  requestAnimationFrame(animate);
}
animate();

Click anywhere to drop a new wave. Watch what happens when two rings collide—they annihilate each other. This is a fundamental property of excitable media: the refractory zone behind each wave prevents it from being re-triggered. The same physics governs cardiac tissue (where spiral waves cause arrhythmias), forest fire spread, and slime mold aggregation. Belousov first observed these waves in a beaker of cerium and citric acid in the 1950s. Nobody believed him for a decade.

7. Belousov-Zhabotinsky spirals — rotating chemical pinwheels

The BZ reaction is the most photogenic demonstration of chemical self-organization. In a real BZ reaction, spiral waves rotate indefinitely, annihilate when they collide, and fill any container with a mesmerizing pattern of spinning pinwheels. This simulation uses the Oregonator model, a simplified version of the BZ reaction kinetics.

const canvas = document.createElement('canvas');
const W = 250, H = 250;
canvas.width = W; canvas.height = H;
canvas.style.width = '500px'; canvas.style.height = '500px';
canvas.style.imageRendering = 'pixelated';
document.body.appendChild(canvas);
document.body.style.margin = '0';
document.body.style.background = '#000';
document.body.style.display = 'flex';
document.body.style.justifyContent = 'center';
document.body.style.alignItems = 'center';
document.body.style.height = '100vh';
const ctx = canvas.getContext('2d');
const img = ctx.createImageData(W, H);

let a = new Float32Array(W * H).fill(0);
let b = new Float32Array(W * H).fill(0);
let c = new Float32Array(W * H).fill(0);
let na = new Float32Array(W * H);
let nb = new Float32Array(W * H);
let nc = new Float32Array(W * H);

// Seed spiral by creating a broken wave front
// Half-plane of activator + perpendicular half of inhibitor
for (let y = 0; y < H; y++) {
  for (let x = 0; x < W; x++) {
    const i = y * W + x;
    if (x > W / 2 - 20 && x < W / 2 + 20 && y < H / 2) {
      a[i] = 1;
    }
    if (y > H / 2 - 20 && y < H / 2 + 20 && x < W / 2) {
      b[i] = 1;
    }
    // Random noise to break symmetry
    a[i] += Math.random() * 0.05;
    b[i] += Math.random() * 0.05;
  }
}

const Da = 1.0, Db = 0.5;
const alpha = 1.0, beta = 1.0, gamma = 1.0;
const dt = 0.15;

function lap(arr, x, y) {
  const i = y * W + x;
  let s = -4 * arr[i];
  s += x > 0 ? arr[i - 1] : arr[i + W - 1];
  s += x < W - 1 ? arr[i + 1] : arr[i - W + 1];
  s += y > 0 ? arr[i - W] : arr[i + (H - 1) * W];
  s += y < H - 1 ? arr[i + W] : arr[i - (H - 1) * W];
  return s;
}

function step() {
  for (let y = 0; y < H; y++) {
    for (let x = 0; x < W; x++) {
      const i = y * W + x;
      const ai = a[i], bi = b[i], ci = c[i];
      // Three-variable Oregonator-like model
      const ra = alpha * (ai * (1 - ai) - bi * (ai - 0.1) / (ai + 0.1));
      const rb = beta * ai - bi;
      na[i] = Math.max(0, Math.min(1, ai + dt * (Da * lap(a, x, y) + ra)));
      nb[i] = Math.max(0, Math.min(1, bi + dt * (Db * lap(b, x, y) + rb * 0.01)));
    }
  }
  [a, na] = [na, a]; [b, nb] = [nb, b];
}

function draw() {
  for (let i = 0; i < W * H; i++) {
    const ai = Math.max(0, Math.min(1, a[i]));
    const bi = Math.max(0, Math.min(1, b[i]));
    const p = i * 4;
    // BZ reaction colors: red-orange-blue like the real thing
    img.data[p] = Math.floor(ai * 255);
    img.data[p + 1] = Math.floor(ai * 80 + bi * 60);
    img.data[p + 2] = Math.floor(bi * 255);
    img.data[p + 3] = 255;
  }
  ctx.putImageData(img, 0, 0);
}

function animate() {
  for (let s = 0; s < 6; s++) step();
  draw();
  requestAnimationFrame(animate);
}
animate();

The broken wavefront curls into a spiral. Given enough time, the tip of the spiral traces a circle (or a more complex trajectory called a “meander”). Multiple spirals compete for territory. Where two spiral arms collide, they annihilate. Arthur Winfree won the MacArthur “genius grant” partly for connecting these chemical spirals to cardiac arrhythmias—the same rotating wave that makes a BZ reaction beautiful can kill you if it forms in heart muscle.

8. Generative reaction-diffusion art — combined composition

This final example combines multiple techniques: Gray-Scott reaction-diffusion with dynamic parameters that slowly drift over time, HSL color mapping based on both chemical concentrations and local gradient direction, and a click interaction that seeds new growth. The result is an ever-changing organic painting.

const canvas = document.createElement('canvas');
const W = 300, H = 300;
canvas.width = W; canvas.height = H;
canvas.style.width = '600px'; canvas.style.height = '600px';
canvas.style.imageRendering = 'pixelated';
document.body.appendChild(canvas);
document.body.style.margin = '0';
document.body.style.background = '#0a0a0f';
document.body.style.display = 'flex';
document.body.style.justifyContent = 'center';
document.body.style.alignItems = 'center';
document.body.style.height = '100vh';
const ctx = canvas.getContext('2d');
const img = ctx.createImageData(W, H);

let U = new Float32Array(W * H).fill(1);
let V = new Float32Array(W * H).fill(0);
let nU = new Float32Array(W * H);
let nV = new Float32Array(W * H);

// Multiple seed points
for (let s = 0; s < 8; s++) {
  const cx = 40 + Math.random() * (W - 80);
  const cy = 40 + Math.random() * (H - 80);
  for (let dy = -6; dy <= 6; dy++) {
    for (let dx = -6; dx <= 6; dx++) {
      if (dx * dx + dy * dy > 36) continue;
      const x = Math.floor(cx + dx), y = Math.floor(cy + dy);
      if (x >= 0 && x < W && y >= 0 && y < H) {
        const i = y * W + x;
        V[i] = 1; U[i] = 0.25;
      }
    }
  }
}

const Du = 0.2097, Dv = 0.105;
let time = 0;

function lap(arr, x, y) {
  const i = y * W + x;
  const l = x > 0 ? arr[i - 1] : arr[i + W - 1];
  const r = x < W - 1 ? arr[i + 1] : arr[i - W + 1];
  const u = y > 0 ? arr[i - W] : arr[i + (H - 1) * W];
  const d = y < H - 1 ? arr[i + W] : arr[i - (H - 1) * W];
  return l + r + u + d - 4 * arr[i];
}

function step() {
  // Parameters drift slowly over time
  const f = 0.034 + Math.sin(time * 0.0001) * 0.006;
  const k = 0.059 + Math.cos(time * 0.00013) * 0.005;
  for (let y = 0; y < H; y++) {
    for (let x = 0; x < W; x++) {
      const i = y * W + x;
      const u = U[i], v = V[i], uvv = u * v * v;
      nU[i] = Math.max(0, Math.min(1, u + Du * lap(U, x, y) - uvv + f * (1 - u)));
      nV[i] = Math.max(0, Math.min(1, v + Dv * lap(V, x, y) + uvv - (f + k) * v));
    }
  }
  [U, nU] = [nU, U]; [V, nV] = [nV, V];
  time++;
}

// Click to seed
canvas.addEventListener('click', (e) => {
  const rect = canvas.getBoundingClientRect();
  const mx = Math.floor((e.clientX - rect.left) / rect.width * W);
  const my = Math.floor((e.clientY - rect.top) / rect.height * H);
  for (let dy = -8; dy <= 8; dy++) {
    for (let dx = -8; dx <= 8; dx++) {
      if (dx * dx + dy * dy > 64) continue;
      const x = mx + dx, y = my + dy;
      if (x >= 0 && x < W && y >= 0 && y < H) {
        const i = y * W + x;
        V[i] = 1; U[i] = 0;
      }
    }
  }
});

function hsl2rgb(h, s, l) {
  const c = (1 - Math.abs(2 * l - 1)) * s;
  const x = c * (1 - Math.abs((h / 60) % 2 - 1));
  const m = l - c / 2;
  let r = 0, g = 0, b = 0;
  if (h < 60) { r = c; g = x; }
  else if (h < 120) { r = x; g = c; }
  else if (h < 180) { g = c; b = x; }
  else if (h < 240) { g = x; b = c; }
  else if (h < 300) { r = x; b = c; }
  else { r = c; b = x; }
  return [(r + m) * 255, (g + m) * 255, (b + m) * 255];
}

function draw() {
  const hueShift = time * 0.02;
  for (let i = 0; i < W * H; i++) {
    const v = V[i], u = U[i];
    // Gradient direction for hue variation
    const x = i % W, y = Math.floor(i / W);
    let gx = 0, gy = 0;
    if (x > 0 && x < W - 1) gx = V[i + 1] - V[i - 1];
    if (y > 0 && y < H - 1) gy = V[i + W] - V[i - W];
    const gradAngle = Math.atan2(gy, gx);
    const hue = ((gradAngle / Math.PI * 180 + 360 + hueShift) % 360);
    const sat = 0.6 + v * 0.4;
    const light = v * 0.6 + 0.05;
    const [r, g, b] = hsl2rgb(hue, sat, light);
    const p = i * 4;
    img.data[p] = Math.min(255, Math.floor(r));
    img.data[p + 1] = Math.min(255, Math.floor(g));
    img.data[p + 2] = Math.min(255, Math.floor(b));
    img.data[p + 3] = 255;
  }
  ctx.putImageData(img, 0, 0);
  ctx.font = '11px monospace';
  ctx.fillStyle = 'rgba(255,255,255,0.6)';
  ctx.fillText('Click to seed growth', 8, 16);
}

function animate() {
  for (let s = 0; s < 8; s++) step();
  draw();
  requestAnimationFrame(animate);
}
animate();

The drifting parameters push the system through different morphological regimes over time—spots might emerge, merge into stripes, fragment back into spots, or dissolve and regrow. The gradient-based coloring maps local flow direction to hue, producing rainbow-like color bands along pattern edges. Click anywhere to introduce new V chemical and watch it interact with existing structures.

The math behind the patterns

Turing’s 1952 paper (“The Chemical Basis of Morphogenesis”) proved something counterintuitive: diffusion, which normally smooths things out, can create pattern when two substances diffuse at different rates. The fast-diffusing inhibitor B can suppress activator A at a distance while A remains strong locally. This creates a local activation, long-range inhibition dynamic. The wavelength of the resulting pattern is roughly proportional to the square root of the ratio of diffusion coefficients: λ ∝ √(D_B/D_A).

The Gray-Scott model adds feeding and killing terms to Turing’s framework. The feed rate f controls how fast the substrate U is replenished from an external reservoir. The kill rate k controls how fast the catalyst V decays. Together they determine which of the roughly dozen distinct pattern classes the system settles into. Pearson (1993) cataloged these into Greek-letter labeled regions. The website Robert Munafo’s xmorphia page maps the complete parameter space with interactive exploration.

The Laplacian operator ∇² measures how much a cell differs from its neighbors—essentially a discrete blur. In 2D on a grid, the simplest version is: ∇²f(x,y) = f(x+1,y) + f(x−1,y) + f(x,y+1) + f(x,y−1) − 4f(x,y). More accurate stencils include diagonal neighbors with appropriate weights.

Reaction-diffusion in nature

For decades, Turing patterns were a mathematical curiosity with no direct biological evidence. That changed in 2012 when Shigeru Kondo’s lab at Osaka University demonstrated that zebrafish stripes form through a genuine Turing mechanism involving interactions between melanophores and xanthophores. Since then, reaction-diffusion has been confirmed in:

• Mouse hair follicle spacing — WNT and DKK proteins act as activator and inhibitor
• Limb digit formation — BMP and SOX9 create the periodic pattern of fingers
• Seashell pigmentation — Hans Meinhardt modeled cone shell patterns in the 1990s
• Coral branching — calcium carbonate deposition follows RD dynamics
• Brain cortex folding — differential growth of gray and white matter produces gyri and sulci

The connection to art runs deep too. Andy Lomas creates award-winning digital sculptures using reaction-diffusion growth (“Morphogenetic Creations”). Nervous System uses RD algorithms to design jewelry and puzzles. Sage Jenson generates otherworldly organisms from Gray-Scott simulations with 3D extensions.

Performance tips

Reaction-diffusion is computationally heavy. Each pixel needs to read its four neighbors, compute two multiplications and several additions, and write back two values. For a 256×256 grid running 8 substeps per frame, that is over half a million operations per animation frame. Some ways to keep it smooth:

• Use typed arrays. Float32Array is 3–5x faster than regular arrays for this kind of work. Avoid object allocations inside the inner loop.
• Ping-pong buffers. Swap two arrays instead of copying. The [U, nU] = [nU, U] trick costs nothing.
• Skip frames for rendering. Run 8–16 simulation steps per rendered frame. The simulation bottleneck is usually the pixel-by-pixel draw, not the math.
• Move to WebGL. Reaction-diffusion maps perfectly to fragment shaders. Each pixel computes independently. A 1024×1024 Gray-Scott simulation runs at 60fps on any modern GPU, where CPU struggles with 256×256.
• Use ImageData directly. Write RGBA bytes into a Uint8ClampedArray and call putImageData once. Avoid per-pixel fillRect or fillStyle changes.

Explore more generative art on Lumitree, where every branch is a unique micro-world built from code. For related techniques, see the cellular automata guide, the noise texture tutorial, or the fluid simulation deep-dive.