Learn Creative Coding (#52) - Reaction-Diffusion: Turing Patterns

Last two episodes we built flocking simulations -- free-moving agents with separation, alignment, and cohesion producing emergent group behavior. Predators, obstacles, formations. Agent-based emergence at its finest. But boids are discrete entities moving through continuous space. Today we go back to the grid, except this time the cells hold continuous floating-point values and the update rules are based on real chemistry.

Reaction-diffusion systems. This is how zebra stripes form. How leopard spots appear. How coral grows its fractal branches. How the patterns on seashells get their ridges and spirals. Alan Turing figured this out in 1952 -- yes, THAT Turing, the codebreaking one -- in a paper called "The Chemical Basis of Morphogenesis." He asked: how does a ball of identical cells become a striped tiger or a spotted cow? His answer: two chemicals diffusing at different rates can spontaneously create stable spatial patterns from nothing.

The math is elegant and the visuals are genuinely stunning. And the implementation is basically a 2D grid update, same concept as Game of Life from episode 48. Different equations, same loop structure. If you followed ep048 and ep049 (continuous automata), you're ready for this.

The idea: two chemicals fighting

Imagine a flat surface covered in two chemicals. Call them A and B. Chemical A is a substrate -- it gets fed into the system from outside at a constant rate. Chemical B is a catalyst -- it converts A into more B via an autocatalytic reaction: when A meets two molecules of B, it becomes a third B molecule. So B grows by eating A. But B also decays -- it's removed from the system at a constant rate.

Here's the key: both chemicals diffuse (spread out over time, like ink in water), but A diffuses FASTER than B. Chemical A spreads quickly, smoothing itself out. Chemical B spreads slowly, forming local concentrations.

This asymmetry is everything. B can't spread fast enough to distribute evenly, so it forms clumps. Those clumps consume the local A supply. A diffuses in from surrounding areas to replace what was consumed. This creates a feedback loop: B concentrations become stable because they're constantly fed by incoming A, and they can't spread because they diffuse too slowly. The result: spots, stripes, rings, spirals -- depending on the exact parameters.

The Gray-Scott model

There are several reaction-diffusion models. We'll use the Gray-Scott model because it's well-studied and produces the most visually diverse patterns. Two parameters (feed rate and kill rate) control the entire pattern space, and different parameter values give you everything from static spots to pulsating blobs to self-replicating cells.

The equations:

A_new = A + (Da * laplacian(A) - A*B*B + f*(1-A)) * dt
B_new = B + (Db * laplacian(B) + A*B*B - (k+f)*B) * dt

Breaking it down piece by piece:

• Da * laplacian(A) -- A diffuses. The laplacian measures how different a cell is from its neighbors. If the cell has more A than its neighbors, the laplacian is negative and A flows out. If it has less, A flows in. Da controls the diffusion speed.
• Db * laplacian(B) -- same for B, but with a smaller diffusion constant (Db < Da).
• -A*B*B -- the reaction. When A meets B^2, A gets consumed. This term removes A from the system.
• +A*B*B -- the same reaction adds B. For every unit of A consumed, a unit of B is created.
• f*(1-A) -- the feed. A is replenished from outside. f is the feed rate. When A is low, more is added. When A is near 1, very little is added. This keeps A concentrations bounded.
• -(k+f)*B -- the kill. B decays at rate k+f. The f part is because B also gets "diluted" by the incoming feed. k is the kill rate.
• dt -- timestep for stability (usually 1.0 for this system).

Standard parameters: Da = 1.0, Db = 0.5 (B diffuses at half the rate). The interesting stuff happens when you vary f and k.

The laplacian: how different am I from my neighbors?

The laplacian of a 2D field at a point is basically "the average of my neighbors minus me." It tells you whether a value is a local peak (negative laplacian, I'm higher than my neighbors) or a local valley (positive laplacian, I'm lower). Diffusion moves stuff from peaks to valleys, smoothing everything out over time.

We approximate it with a 3x3 convolution kernel:

// laplacian kernel weights
// center = -1, direct neighbors = 0.2, diagonals = 0.05
// this is a weighted average that gives diagonals less influence

function laplacian(grid, x, y, w, h) {
  const c = grid[y * w + x];

  const t = grid[((y - 1 + h) % h) * w + x];
  const b = grid[((y + 1) % h) * w + x];
  const l = grid[y * w + ((x - 1 + w) % w)];
  const r = grid[y * w + ((x + 1) % w)];

  const tl = grid[((y - 1 + h) % h) * w + ((x - 1 + w) % w)];
  const tr = grid[((y - 1 + h) % h) * w + ((x + 1) % w)];
  const bl = grid[((y + 1) % h) * w + ((x - 1 + w) % w)];
  const br = grid[((y + 1) % h) * w + ((x + 1) % w)];

  return t * 0.2 + b * 0.2 + l * 0.2 + r * 0.2
       + tl * 0.05 + tr * 0.05 + bl * 0.05 + br * 0.05
       - c;
}

The weights sum to: 4 * 0.2 + 4 * 0.05 = 1.0. Then we subtract the center (multiply by -1). So the total kernel is: [[0.05, 0.2, 0.05], [0.2, -1, 0.2], [0.05, 0.2, 0.05]]. Direct neighbors get more weight because they're closer. Diagonals get less because they're further (root-2 times further).

The wrapping (% h, % w) gives us toroidal boundary conditions -- same as Game of Life. Patterns wrap around the edges smoothly.

Basic implementation: CPU version

Let's build it. Two flat arrays for A and B, a second pair for the "next" state (double buffering, same as ep048), and the update loop:

const canvas = document.getElementById('canvas');
const ctx = canvas.getContext('2d');
const W = 256;
const H = 256;
canvas.width = W;
canvas.height = H;

const Da = 1.0;
const Db = 0.5;
let f = 0.035;   // feed rate
let k = 0.065;   // kill rate
const dt = 1.0;

// two grids for current state
let gridA = new Float32Array(W * H);
let gridB = new Float32Array(W * H);
// two grids for next state
let nextA = new Float32Array(W * H);
let nextB = new Float32Array(W * H);

// initial conditions: A = 1 everywhere, B = 0 everywhere
gridA.fill(1.0);
gridB.fill(0.0);

// seed: drop some B in the center
for (let dy = -10; dy <= 10; dy++) {
  for (let dx = -10; dx <= 10; dx++) {
    const x = Math.floor(W / 2) + dx;
    const y = Math.floor(H / 2) + dy;
    if (x >= 0 && x < W && y >= 0 && y < H) {
      gridB[y * W + x] = 1.0;
      gridA[y * W + x] = 0.0;
    }
  }
}

The initial state is crucial. We fill A with 1.0 (substrate everywhere) and B with 0.0 (no catalyst). Then we "seed" a small square of B in the center. This seed is where the pattern will grow from. Without it, nothing happens -- the system is in equilibrium when B = 0 everywhere. The seed breaks the symmetry and the pattern radiates outward.

The update step

Each frame, we compute the laplacian for both A and B at every cell, then apply the Gray-Scott equations:

function step() {
  for (let y = 0; y < H; y++) {
    for (let x = 0; x < W; x++) {
      const idx = y * W + x;
      const a = gridA[idx];
      const b = gridB[idx];

      const lapA = laplacian(gridA, x, y, W, H);
      const lapB = laplacian(gridB, x, y, W, H);

      const reaction = a * b * b;

      nextA[idx] = a + (Da * lapA - reaction + f * (1.0 - a)) * dt;
      nextB[idx] = b + (Db * lapB + reaction - (k + f) * b) * dt;

      // clamp to [0, 1]
      if (nextA[idx] < 0) nextA[idx] = 0;
      if (nextA[idx] > 1) nextA[idx] = 1;
      if (nextB[idx] < 0) nextB[idx] = 0;
      if (nextB[idx] > 1) nextB[idx] = 1;
    }
  }

  // swap buffers
  const tmpA = gridA;
  const tmpB = gridB;
  gridA = nextA;
  gridB = nextB;
  nextA = tmpA;
  nextB = tmpB;
}

This is exactly the same double-buffer pattern from Game of Life. Read from current grid, write to next grid, swap. The clamping keeps values in [0, 1] -- without it, numerical instabilities can push values negative or beyond 1, which causes the simulation to explode.

The a * b * b reaction term is the heart of the system. It's nonlinear (B squared), which is what creates the interesting dynamics. If B is low, the reaction barely fires. If B is high, the reaction fires hard, consuming A quickly and creating even more B. This positive feedback is what makes B concentrations self-reinforcing -- a spot of B tends to grow rather than shrink, until the local A supply is depleted.

Rendering: mapping chemistry to color

The simplest rendering: use the B concentration as brightness. Where B is high, draw light pixels. Where B is low, draw dark pixels.

const imgData = ctx.createImageData(W, H);

function render() {
  for (let i = 0; i < W * H; i++) {
    const b = gridB[i];
    const val = Math.floor(b * 255);
    imgData.data[i * 4 + 0] = val;
    imgData.data[i * 4 + 1] = val;
    imgData.data[i * 4 + 2] = val;
    imgData.data[i * 4 + 3] = 255;
  }
  ctx.putImageData(imgData, 0, 0);
}

This gives you white patterns on a black background. The spots and stripes that form are white -- regions where B accumulated. The black areas are "substrate territory" -- A dominates there and B was killed off.

For something more interesting, map B to a color ramp:

function render() {
  for (let i = 0; i < W * H; i++) {
    const b = gridB[i];
    const a = gridA[i];

    // use the ratio for a smoother gradient
    const ratio = b / (a + b + 0.001);

    // cool palette: dark blue -> teal -> white
    let r, g, bl;
    if (ratio < 0.3) {
      const t = ratio / 0.3;
      r = Math.floor(10 + t * 20);
      g = Math.floor(15 + t * 60);
      bl = Math.floor(40 + t * 80);
    } else if (ratio < 0.6) {
      const t = (ratio - 0.3) / 0.3;
      r = Math.floor(30 + t * 80);
      g = Math.floor(75 + t * 130);
      bl = Math.floor(120 + t * 60);
    } else {
      const t = (ratio - 0.6) / 0.4;
      r = Math.floor(110 + t * 145);
      g = Math.floor(205 + t * 50);
      bl = Math.floor(180 + t * 75);
    }

    imgData.data[i * 4 + 0] = r;
    imgData.data[i * 4 + 1] = g;
    imgData.data[i * 4 + 2] = bl;
    imgData.data[i * 4 + 3] = 255;
  }
  ctx.putImageData(imgData, 0, 0);
}

Using the ratio B/(A+B) instead of raw B gives smoother gradients. The edges of patterns become visible as color transitions rather than sharp binary boundaries. Different palettes produce wildly different aesthetics from the same simulation -- coral tones for organic, neon for digital, earth tones for biological.

The parameter space: spots vs stripes vs chaos

The feed rate f and kill rate k determine what kind of pattern forms. Karl Sims created a famous map of the Gray-Scott parameter space -- a 2D grid where each pixel represents one (f, k) combination and shows the pattern that emerges. It's one of the most beautiful images in computational art.

Here are some parameter combos to try:

// spots (mitosis-like self-replication)
f = 0.035; k = 0.065;

// stripes and worms
f = 0.025; k = 0.060;

// spirals (rare, unstable)
f = 0.014; k = 0.054;

// coral growth
f = 0.055; k = 0.062;

// moving spots (solitons)
f = 0.030; k = 0.062;

// maze-like patterns
f = 0.029; k = 0.057;

// holes in a solid field
f = 0.039; k = 0.058;

The spots at f=0.035, k=0.065 are the classic starting point. You'll see the initial seed grow into a cluster of round spots that split (mitosis!) and spread across the grid. Each spot is a stable B concentration maintained by the balance between the reaction (creating B), the kill term (destroying B), and diffusion (spreading B outward).

Stripes form when spots merge along one axis. The B concentrations elongate into worms rather than staying round. This happens at slightly lower feed rates -- less A is available, so the spots can't sustain themselves individually and merge into stripes to share resources. It's an efficiency thing. Biologically, this is how animal stripe patterns work.

Multiple seeds and interactive seeding

One centered seed is fine for testing, but the creative potential explodes when you seed multiple points or let the user seed with the mouse:

canvas.addEventListener('mousemove', function(e) {
  if (!mouseDown) return;
  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));

  // paint a small circle of B at mouse position
  const radius = 3;
  for (let dy = -radius; dy <= radius; dy++) {
    for (let dx = -radius; dx <= radius; dx++) {
      if (dx * dx + dy * dy > radius * radius) continue;
      const px = (mx + dx + W) % W;
      const py = (my + dy + H) % H;
      const idx = py * W + px;
      gridB[idx] = 1.0;
      gridA[idx] = 0.0;
    }
  }
});

let mouseDown = false;
canvas.addEventListener('mousedown', function() { mouseDown = true; });
canvas.addEventListener('mouseup', function() { mouseDown = false; });

Click and drag to paint B chemical onto the grid. Watch patterns grow from where you draw. You can seed a straight line and watch stripes form. Seed scattered dots and watch spots grow and interact. Draw a circle and watch it break up into spots or collapse into a ring depending on the parameters.

This is where it becomes a creative tool rather than a simulation. You're sculpting the initial conditions and the chemistry does the rest.

Speed: multiple steps per frame

Reaction-diffusion is slow. You need hundreds or thousands of iterations before the pattern stabilizes. Running one step per frame at 60fps means waiting minutes for anything interesting. The fix: run multiple update steps per frame render:

const STEPS_PER_FRAME = 10;

function loop() {
  for (let i = 0; i < STEPS_PER_FRAME; i++) {
    step();
  }
  render();
  requestAnimationFrame(loop);
}

loop();

At 256x256 with 10 steps per frame, the CPU can handle this fine in JavaScript. You'll see the pattern develop in real time -- growing, splitting, connecting. For higher resolution (512x512 or beyond), you'd want a GPU implementation, which is essentially the same algorithm but running as a fragment shader with ping-pong framebuffers. The shader version runs at full 1080p no problem because every pixel updates in parallel.

Parameter animation: morphing patterns

Here's where it gets really creative. Slowly change f and k over time and watch the pattern morph between different regimes:

let time = 0;

function animateParams() {
  time += 0.001;

  // oscillate between spots and stripes
  f = 0.030 + Math.sin(time) * 0.008;
  k = 0.061 + Math.cos(time * 0.7) * 0.004;
}

function loop() {
  animateParams();
  for (let i = 0; i < STEPS_PER_FRAME; i++) {
    step();
  }
  render();
  requestAnimationFrame(loop);
}

Spots slowly stretch into stripes. Stripes break apart into spots. New spots grow from nothing. Old spots die and dissolve. The pattern is never static -- it's always in transition, always finding a new equilibrium that shifts before it fully stabilizes. It looks alive in a way that static parameters don't. Like watching bacteria multiply under a microscope.

You can also map mouse position to f and k. Move left for stripes, move right for spots. Move up for sparse patterns, move down for dense ones. The pattern responds in real time, morphing under your cursor. It's basically the Karl Sims parameter map but interactive and continuous.

Spatial parameter variation

Even more powerful: make f and k vary across the grid itself. Different parts of the surface have different chemistries, producing different patterns that meet at boundaries:

function getParams(x, y) {
  // left half: spots, right half: stripes
  const t = x / W;
  const localF = 0.025 + t * 0.015;  // 0.025 -> 0.040
  const localK = 0.058 + t * 0.008;  // 0.058 -> 0.066

  return { f: localF, k: localK };
}

function stepWithSpatialParams() {
  for (let y = 0; y < H; y++) {
    for (let x = 0; x < W; x++) {
      const idx = y * W + x;
      const a = gridA[idx];
      const b = gridB[idx];
      const params = getParams(x, y);

      const lapA = laplacian(gridA, x, y, W, H);
      const lapB = laplacian(gridB, x, y, W, H);
      const reaction = a * b * b;

      nextA[idx] = a + (Da * lapA - reaction + params.f * (1.0 - a)) * dt;
      nextB[idx] = b + (Db * lapB + reaction - (params.k + params.f) * b) * dt;

      if (nextA[idx] < 0) nextA[idx] = 0;
      if (nextA[idx] > 1) nextA[idx] = 1;
      if (nextB[idx] < 0) nextB[idx] = 0;
      if (nextB[idx] > 1) nextB[idx] = 1;
    }
  }

  const tmpA = gridA; const tmpB = gridB;
  gridA = nextA; gridB = nextB;
  nextA = tmpA; nextB = tmpB;
}

The left side produces stripes. The right side produces spots. In between, you get a transiton zone where stripes break into spots or spots merge into stripes. The boundary region is where the most interesting things happen -- neither pattern dominates, and you get novel hybrid structures that neither parameter set produces alone.

You could map the parameters to a Perlin noise field (from episode 12) instead of a linear gradient. Then you'd get islands of different pattern types scattered across the grid, with organic boundaries between them. The patterns would look like a satellite photo of different biomes meeting at their edges.

The connection to nature

This isn't just abstract math. Reaction-diffusion is a real physical process. The patterns it produces match real biological patterns with remarkable accuracy:

• Leopard spots: f ~ 0.035, k ~ 0.065
• Zebra stripes: f ~ 0.025, k ~ 0.060
• Giraffe patches: similar to spots but with different boundary conditions
• Coral branches: f ~ 0.055, k ~ 0.062 (the branching growth regime)
• Fingerprint ridges: stripe-forming parameters on a curved surface
• Seashell patterns: 1D reaction-diffusion along the shell's growing edge, accumulated over time

Turing proposed this in 1952. It took until the 1990s for biologists to find the actual morphogen molecules in animal skin cells. He was right -- the mechanisms he described mathematically are the actual mechanisms used by nature to create pattern. Two chemicals, different diffusion rates, autocatalytic reaction. That's all it takes.

The fact that we can recreate these natural patterns in 50 lines of JavaScript running on a canvas element is... kind of wild to think about. The same math that generates the stripes on a zebrafish is running in your browser right now. Makes sense, right? Physics doesn't care what's computing it :-)

Creative exercise: the full interactive version

Put it all together. Interactive seeding, parameter control, color mapping:

const canvas = document.getElementById('canvas');
const ctx = canvas.getContext('2d');
const W = 200, H = 200;
canvas.width = W; canvas.height = H;
canvas.style.width = '600px'; canvas.style.height = '600px';
canvas.style.imageRendering = 'pixelated';

let gA = new Float32Array(W * H).fill(1);
let gB = new Float32Array(W * H).fill(0);
let nA = new Float32Array(W * H);
let nB = new Float32Array(W * H);
let f = 0.035, k = 0.065;
const img = ctx.createImageData(W, H);

function lap(g, x, y) {
  const i = y * W + x;
  const u = ((y-1+H)%H)*W+x, d = ((y+1)%H)*W+x;
  const l = y*W+((x-1+W)%W), r = y*W+((x+1)%W);
  return g[u]*0.2 + g[d]*0.2 + g[l]*0.2 + g[r]*0.2
       + g[((y-1+H)%H)*W+((x-1+W)%W)]*0.05
       + g[((y-1+H)%H)*W+((x+1)%W)]*0.05
       + g[((y+1)%H)*W+((x-1+W)%W)]*0.05
       + g[((y+1)%H)*W+((x+1)%W)]*0.05
       - g[i];
}

function step() {
  for (let y = 0; y < H; y++) for (let x = 0; x < W; x++) {
    const i = y * W + x;
    const a = gA[i], b = gB[i];
    const r = a * b * b;
    nA[i] = Math.max(0, Math.min(1, a + lap(gA,x,y) - r + f*(1-a)));
    nB[i] = Math.max(0, Math.min(1, b + lap(gB,x,y)*0.5 + r - (k+f)*b));
  }
  [gA, nA] = [nA, gA]; [gB, nB] = [nB, gB];
}

function draw() {
  for (let i = 0; i < W*H; i++) {
    const v = Math.floor((1 - gA[i] + gB[i]) * 0.5 * 255);
    img.data[i*4] = v * 0.3;
    img.data[i*4+1] = v * 0.7;
    img.data[i*4+2] = v;
    img.data[i*4+3] = 255;
  }
  ctx.putImageData(img, 0, 0);
}

// seed center
for (let dy = -5; dy <= 5; dy++)
  for (let dx = -5; dx <= 5; dx++) {
    const i = (H/2+dy)*W + W/2+dx;
    gB[i] = 1; gA[i] = 0;
  }

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

This is a compact version you can paste into an HTML file with a <canvas id="canvas"> element. It'll show you the classic spot-forming reaction-diffusion pattern growing from a center seed. Tweak f and k to explore different pattern regimes. Add mouse seeding with the handler from earlier. Map the output to different color palettes.

The connection forward

Reaction-diffusion is a continuous-field emergent system. Unlike boids (discrete agents in continuous space) or Game of Life (discrete agents on a discrete grid), reaction-diffusion operates on continuous values over a continuous field. The grid is just our numerical approximation -- the real system is a PDE (partial differential equation) that describes concentration changes at every point in space simultaneously.

This pattern -- local interactions producing global structure -- keeps showing up. Episode 47 through 51 covered discrete agents (cells, boids). Now we're doing continuous fields. Coming up we'll look at L-systems, where simple grammar rules grow complex branching structures like trees and ferns. Different math, same principle: local rules, no master plan, complex emergent output.

Reaction-diffusion also connects forward to the GPU shader work we did in episodes 32-46. The update loop here is a perfect fit for a fragment shader -- each pixel computes its own update based on neighbor values, reads from one texture, writes to another (ping-pong). The GPU version runs at 1080p in real time because every pixel updates in parallel. If you want to explore high-res reaction-diffusion, port this to GLSL. The shader is almost identical to the JavaScript, just with texture reads instead of array indexing.

Six episodes into the emergent systems arc now. We went from 1D binary automata (47) to 2D binary automata (48) to continuous automata (49) to discrete agents (50-51) to continuous chemical fields (this one). Each step explores a different combination of discrete vs continuous, grid vs free, agent vs field. And they all produce that same magic: patterns that nobody programmed, growing from rules that fit in a tweet :-)

't Komt erop neer...

• Reaction-diffusion simulates two chemicals (A and B) diffusing and reacting on a 2D grid. A is fed in from outside, B is killed off, and A converts to B via an autocatalytic reaction (A + 2B -> 3B). The asymmetry in diffusion rates (A fast, B slow) causes B to form stable spatial concentrations
• The Gray-Scott model uses two parameters: feed rate (f) and kill rate (k). Different (f, k) values produce spots, stripes, spirals, coral growth, maze patterns, and moving solitons. Karl Sims' parameter map visualizes the entire pattern space
• The laplacian approximation uses a 3x3 weighted kernel (center = -1, edges = 0.2, corners = 0.05) to measure how different a cell is from its neighbors. Diffusion moves values from peaks toward valleys, smoothing the field
• Implementation uses double-buffered Float32Arrays, same pattern as Game of Life. Initial condition: A = 1 everywhere, B = 0 everywhere, with a small seed of B to break symmetry. Multiple steps per frame (8-10) for reasonable speed
• Interactive seeding (mouse painting B onto the grid) turns it into a creative tool. Parameter animation (slowly varying f and k) morphs patterns between regimes in real time. Spatial parameter variation creates zones with different pattern types meeting at boundaries
• Alan Turing proposed this mechanism in 1952 to explain biological pattern formation. The same math produces leopard spots, zebra stripes, coral branches, and fingerprint ridges. The morphogen molecules were confirmed experimentally in the 1990s
• The technique connects backward to GPU shaders (ping-pong textures for high-res real-time) and forward to other emergent systems like L-systems and agent-based growth simulations

Sallukes! Thanks for reading.

X

Hive account@femdev