Reaction–diffusion

Reaction and diffusion

Two chemicals spreading and eating each other. Four numbers decide whether you get spots, stripes, a maze, or blobs that split like cells. Turing wrote the equations in 1952 to ask how an animal decides where to put its markings.

EXP-021

Move the cursor over the canvas — it is part of the simulation.

Drag your finger across the canvas — it is part of the simulation.

What to look for

Draw on it with the cursor: you are injecting V, and whatever you draw grows into the pattern the parameters allow. Then move k slowly. Around F=0.037, k=0.06 you get spots that divide like bacteria. Push k up and they freeze into a static maze. Push it down and everything floods.

The maths

∂U/∂t = Du∇²U − UV² + F(1−U)

U diffuses, gets eaten by the reaction UV², and is topped back up toward 1 at rate F. That last term is the feed: without it the system runs out and dies flat.

∂V/∂t = Dv∇²V + UV² − (F+k)V

V is produced by the same reaction that consumes U, which makes it autocatalytic: V needs V to make more V. It is removed at rate F+k. The entire zoo of patterns lives in the balance between that self-feeding and that removal.

Du / Dv = 2

The one condition that makes patterns possible: the inhibitor must spread faster than the activator. If they diffused at the same rate the sheet would stay uniform forever. Turing called it diffusion-driven instability, and it is deeply counterintuitive, because diffusion is supposed to smooth things out.

∇²U ≈ U↑ + U↓ + U← + U→ − 4U

The Laplacian on a grid, five points. It measures how much a cell differs from the average of its neighbours: positive if it sits in a dip, negative on a bump. That single number is all the geometry the simulation ever knows.

How it is built

A 150×150 grid of two Float32Arrays, explicit Euler in time, five-point Laplacian in space. Two extra buffers hold the next state so a cell never reads a half-updated neighbour, and the pairs get swapped instead of copied. It runs several steps per frame because the timestep has to stay small for the integration to be stable.

The code, step by step

Cut straight from lab.js. These are the real lines that run above, not a simplified version.

01

Seed the disturbance

const seed = () => {
  U = new Float32Array(N * N).fill(1); V = new Float32Array(N * N);
  U2 = new Float32Array(N * N); V2 = new Float32Array(N * N);
  for (let q = 0; q < 22; q++) {
    const cx = 10 + ((Math.random() * (N - 20)) | 0), cy = 10 + ((Math.random() * (N - 20)) | 0);
    for (let dy = -4; dy <= 4; dy++) for (let dx = -4; dx <= 4; dx++) {
      const k = (cy + dy) * N + (cx + dx);
      U[k] = 0.5; V[k] = 0.25;

U starts at 1 everywhere and V at 0, which is a perfectly stable state: leave it alone and nothing happens, forever. The patterns need a defect to grow from, so a handful of square patches get knocked to a different mix. The pattern is a response to damage.

02

One Euler step of both equations

const Du = 0.16, Dv = 0.08, F = P.F, kk = P.k;
for (let s = 0; s < (P.it | 0); s++) {
  for (let y = 1; y < N - 1; y++) {
    for (let x = 1; x < N - 1; x++) {
      const i = y * N + x;
      // Laplaciano de 5 puntos
      const lu = U[i - 1] + U[i + 1] + U[i - N] + U[i + N] - 4 * U[i];
      const lv = V[i - 1] + V[i + 1] + V[i - N] + V[i + N] - 4 * V[i];
      const uvv = U[i] * V[i] * V[i];
      U2[i] = U[i] + Du * lu - uvv + F * (1 - U[i]);
      V2[i] = V[i] + Dv * lv + uvv - (F + kk) * V[i];

The five-point Laplacian is the four neighbours minus four times the centre — the discrete version of ∇². Then the two equations, written literally. Du is twice Dv, and that ratio is the entire condition for patterns to be possible: the inhibitor has to outrun the activator.

03

Swap, do not copy

const tu = U; U = U2; U2 = tu;
const tv = V; V = V2; V2 = tv;

The new state is written into a second pair of buffers so no cell ever reads a neighbour that has already been updated this step, which would silently turn the explicit scheme into something else. Then the references are swapped, which is free, instead of copying 22.500 floats twice per iteration.

Why it matters

Turing published this the year before he died, and it is the least famous of his big ideas. It gives a mechanism for how an undifferentiated sheet of cells can end up with stripes in the right places without any cell being told the plan. Zebrafish stripe mutants have since been shown to behave the way the equations predict.

The whole thing

Everything above, in one piece. This is the entire experiment as it lives in lab.js.

FULL

lab.js · grayscott

{ id: 'grayscott',
      name: L('Reaction and diffusion', 'Reacción y difusión'),
      note: L('Two chemicals, four numbers. Spots, stripes, mazes, things that divide.',
              'Dos químicos, cuatro números. Manchas, rayas, laberintos, cosas que se dividen.'),
      params: [
        { key: 'F', label: L('Feed F', 'Alimentación F'), min: 0.01, max: 0.08, step: 0.001, def: 0.037 },
        { key: 'k', label: L('Kill k', 'Remoción k'), min: 0.045, max: 0.07, step: 0.0005, def: 0.06 },
        { key: 'it', label: L('Steps per frame', 'Pasos por frame'), min: 1, max: 14, step: 1, def: 6 }
      ],
      make() {
        const N = 150;
        let U = null, V = null, U2 = null, V2 = null;
        const seed = () => {
          U = new Float32Array(N * N).fill(1); V = new Float32Array(N * N);
          U2 = new Float32Array(N * N); V2 = new Float32Array(N * N);
          for (let q = 0; q < 22; q++) {
            const cx = 10 + ((Math.random() * (N - 20)) | 0), cy = 10 + ((Math.random() * (N - 20)) | 0);
            for (let dy = -4; dy <= 4; dy++) for (let dx = -4; dx <= 4; dx++) {
              const k = (cy + dy) * N + (cx + dx);
              U[k] = 0.5; V[k] = 0.25;
            }
          }
        };
        seed();
        return {
          reset() { seed(); },
          step(ctx, w, h, t, acc, P, M) {
            if (M.in) {                                  // sembrar con el cursor
              const S0 = Math.min(w, h) * 0.8, ox0 = (w - S0) / 2, oy0 = (h - S0) / 2;
              const cx = ((M.x - ox0) / S0 * N) | 0, cy = ((M.y - oy0) / S0 * N) | 0;
              for (let dy = -3; dy <= 3; dy++) for (let dx = -3; dx <= 3; dx++) {
                const x = cx + dx, y = cy + dy;
                if (x > 0 && x < N - 1 && y > 0 && y < N - 1) { U[y * N + x] = 0.5; V[y * N + x] = 0.25; }
              }
            }
            const Du = 0.16, Dv = 0.08, F = P.F, kk = P.k;
            for (let s = 0; s < (P.it | 0); s++) {
              for (let y = 1; y < N - 1; y++) {
                for (let x = 1; x < N - 1; x++) {
                  const i = y * N + x;
                  // Laplaciano de 5 puntos
                  const lu = U[i - 1] + U[i + 1] + U[i - N] + U[i + N] - 4 * U[i];
                  const lv = V[i - 1] + V[i + 1] + V[i - N] + V[i + N] - 4 * V[i];
                  const uvv = U[i] * V[i] * V[i];
                  U2[i] = U[i] + Du * lu - uvv + F * (1 - U[i]);
                  V2[i] = V[i] + Dv * lv + uvv - (F + kk) * V[i];
                }
              }
              const tu = U; U = U2; U2 = tu;
              const tv = V; V = V2; V2 = tv;
            }
            const img = ctx.createImageData(N, N), px = img.data;
            const rgb = acc.split(',').map(Number);
            for (let i = 0; i < N * N; i++) {
              const v = Math.min(1, Math.max(0, V[i] * 3.4));
              const o = i * 4;
              px[o] = rgb[0] * v; px[o+1] = rgb[1] * v; px[o+2] = rgb[2] * v; px[o+3] = 255;
            }
            const off = document.createElement('canvas');
            off.width = off.height = N; off.getContext('2d').putImageData(img, 0, 0);
            const S = Math.min(w, h) * 0.8;
            ctx.imageSmoothingEnabled = true;
            ctx.drawImage(off, (w - S) / 2, (h - S) / 2, S, S);

            ctx.font = '11px monospace';
            ctx.fillStyle = `rgba(${acc},0.55)`;
            ctx.fillText('∂U/∂t = Du∇²U − UV² + F(1−U)      ∂V/∂t = Dv∇²V + UV² − (F+k)V', 12, h - 28);
            ctx.fillStyle = `rgba(${acc},0.95)`;
            ctx.fillText('F = ' + F.toFixed(3) + '   k = ' + kk.toFixed(4), 12, h - 12);
          }
        };
      }
    },