# Learn Creative Coding (#52) - Reaction-Diffusion: Turing Patterns ![cc-banner](https://images.hive.blog/DQmX7rSRusqv16SWFmWy4UNwvjM6hq635UEA9766RR6qHLr/cc-banner-blue.png) 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: ```javascript // 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: ```javascript 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: ```javascript 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. ```javascript 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: ```javascript 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: ```javascript // 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: ```javascript 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: ```javascript 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: ```javascript 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: ```javascript 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: ```javascript 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 `

Time to stretch, kitty kitty.

Time to stretch, kitty kitty.