Reaction-diffusion: how chemistry paints leopard spots

1. Reaction, diffusion and self-organisation

У 1952 році Алан Тюрінг опублікував статтю «Хімічна основа морфогенезу», де показав, що дифузна нестабільність може виникати навіть у системі, яка однорідна та стабільна без дифузії. Це здається парадоксом: дифузія зазвичай вирівнює концентрації — але якщо два реагенти дифундують з різними швидкостями, може виникнути просторовий паттерн.

Загальний механізм: активатор (A) стимулює своє власне виробництво та виробництво інгібітора (B); інгібітор пригнічує активатор; при цьому інгібітор дифундує швидше. Результат — локальне збудження активатора, оточене зоною придушення: плями, смуги, спіралі.

2. General reaction-diffusion equations

Загальна форма для двох компонент — це система двох PDE (рівнянь у частинних похідних):

Різні вибори f і g дають різні моделі: Бріссель (Prigogine), Гіерера-Майнгардта, Фіцгью-Нагумо (збудливі нейрони), Ленгтон (цифровий антропоморфізм) — і нашу найпростішу, найзручнішу для симуляції: модель Грея-Скотта.

3. The Gray-Scott model

Грей та Скотт (1984, 1985) описали реакцію у відкритому реакторі, де речовина A подається ззовні і реагує з B за схемою: A + 2B → 3B, B → P (P — неактивний продукт).

Розберемо доданки:

• Dₐ·∇²A — дифузія A у просторі.
• −A·B² — A витрачається на реакцію (автокаталіз: B каталізує само-виробництво).
• +F·(1−A) — зовнішнє поповнення A до рівноважної концентрації 1.
• +A·B² — B виробляється тою самою реакцією.
• −(F+k)·B — B виводиться (розбавляється + перетворюється на P).

4. Parameters F and k — pattern zoo

Невеликі зміни у F та k кардинально змінюють морфологію паттернів. Нижче — класифікація основних режимів:

5. FTCS discretisation

FTCS (Forward Time, Central Space) — найпростіша явна схема різниць. Лапласіан апроксимується п'ятиточковим шаблоном (4-зв'язність):

Умова стійкості Куранта: dt · D / h² ≤ 0.5. При Dₐ=0.21, h=1, dt=1 — 0.21·1/1 = 0.21 < 0.5. Але для більш стійкого результату часто роблять кілька кроків FTCS за один "відображений" кадр.

6. JavaScript Canvas implementation

const W = 256, H = 256;
const Da = 0.21, Db = 0.105;
const F  = 0.035, k  = 0.065;  // ← spots
const dt = 1.0;

// Grid initialisation
let A = new Float32Array(W * H).fill(1);
let B = new Float32Array(W * H).fill(0);
let nA = new Float32Array(W * H);
let nB = new Float32Array(W * H);

// Initial excitation
for (let seed = 0; seed < 5; seed++) {
  const cx = Math.floor(Math.random()*W);
  const cy = Math.floor(Math.random()*H);
  for (let dy = -5; dy <= 5; dy++) {
    for (let dx = -5; dx <= 5; dx++) {
      const x = (cx+dx+W) % W, y = (cy+dy+H) % H;
      A[y*W+x] = 0.5;
      B[y*W+x] = 0.25 + Math.random()*0.1;
    }
  }
}

function step() {
  for (let y = 0; y < H; y++) {
    for (let x = 0; x < W; x++) {
      const i = y*W + x;

      // Laplacian (wrap-around)
      const la = A[((y-1+H)%H)*W+x] + A[((y+1)%H)*W+x]
               + A[y*W+(x-1+W)%W] + A[y*W+(x+1)%W] - 4*A[i];
      const lb = B[((y-1+H)%H)*W+x] + B[((y+1)%H)*W+x]
               + B[y*W+(x-1+W)%W] + B[y*W+(x+1)%W] - 4*B[i];

      const ab2 = A[i] * B[i] * B[i];
      nA[i] = A[i] + dt*(Da*la - ab2 + F*(1-A[i]));
      nB[i] = B[i] + dt*(Db*lb + ab2 - (F+k)*B[i]);

      // Clamp [0,1]
      if (nA[i] < 0) nA[i] = 0;
      if (nB[i] < 0) nB[i] = 0;
    }
  }
  [A, nA] = [nA, A]; [B, nB] = [nB, B]; // swap buffers
}

function draw(ctx, imageData) {
  const d = imageData.data;
  for (let i = 0, p = 0; i < W*H; i++, p+=4) {
    const v = Math.floor(B[i] * 255);
    d[p]=v; d[p+1]=v>>1; d[p+2]=255-v; d[p+3]=255;
  }
  ctx.putImageData(imageData, 0, 0);
}

// Loop: 8 steps per frame for visible progress
function loop() {
  for (let s = 0; s < 8; s++) step();
  draw(ctx, imageData);
  requestAnimationFrame(loop);
}

При 256×256 і 8 кроків за кадр отримуємо ~350 000 оновлень клітин за кадр. На сучасному браузері це виконується за <8 мс (менше ніж бюджет 16 мс для 60 FPS).

7. WebGL2: GPU simulation

Для роздільної здатності 512×512+ CPU-підхід повільний. WebGL2 дозволяє виконати кроки FTCS у фрагментному шейдері паралельно для всіх пікселів:

// Fragment shader (GLSL ES 3.0)
#version 300 es
precision highp float;
uniform sampler2D u_state;   // rgba: r=A, g=B
uniform vec2      u_res;
out vec4 fragColor;

void main() {
  vec2 uv = gl_FragCoord.xy / u_res;
  vec2 texel = 1.0 / u_res;

  vec2 cur = texture(u_state, uv).rg;
  vec2 top = texture(u_state, uv + vec2(0, texel.y)).rg;
  vec2 bot = texture(u_state, uv - vec2(0, texel.y)).rg;
  vec2 lft = texture(u_state, uv - vec2(texel.x, 0)).rg;
  vec2 rgt = texture(u_state, uv + vec2(texel.x, 0)).rg;

  vec2 lap = top + bot + lft + rgt - 4.0 * cur;

  float F = 0.035, k = 0.065;
  float Da = 0.21, Db = 0.105;
  float ab2 = cur.r * cur.g * cur.g;

  float newA = cur.r + Da*lap.r - ab2 + F*(1.0 - cur.r);
  float newB = cur.g + Db*lap.g + ab2 - (F+k)*cur.g;

  fragColor = vec4(clamp(newA, 0.0, 1.0),
                   clamp(newB, 0.0, 1.0),
                   0.0, 1.0);
}

Техніка ping-pong буферів: два framebuffer-и (FBO) чергуються — один є вхідним texture, інший — вихідним render target. Це дозволяє робити 30–50 кроків за кадр при 512×512 без будь-якого гальмування.

8. Connection to Turing and biology

Тюрінг (1952) довів: якщо інгібітор дифундує швидше за активатора, однорідний стан стає нестійким і виникає просторовий паттерн. Gray-Scott реалізує саме такий механізм: Dₐ > D_b і A грає роль «активатора» (само-каталітично утворює B).

Реальні приклади «тюрінгових паттернів»:

• Шкіра риби фугу — чорно-білі смуги, передбачені Gray-Scott за F≈0.06, k≈0.06.
• Плями гепарда — ізольовані острівці активатора пігменту, F≈0.03, k≈0.065.
• Пальці руки — позиція зачатків пальців у ебріоні поза shumway залежить від реакційно-дифузійних паттернів Hox-генів.
• BZ реакція (Белоусов–Жаботинський) — реальна хімічна система, що генерує спіральні хвилі подібні до моделі Фіцгью-Нагумо.

9. Extensions

• 3D симуляція: той самий алгоритм у кубічній сітці породжує 3D паттерни (губчасті структури, трубчасті лабіринти). Вимагає 3D текстур у WebGL або WebGPU.
• Нерівномірні параметри: F та k можуть бути полями — зображення-текстура задає «карту паттернів», де різні ділянки ростуть по-різному.
• Multi-species: додаємо третю компоненту C, яка пригнічує B — складніша екосистемна динаміка (Lotka-Volterra на просторовій сітці).
• Додаток до клітинного автомата: замість неперервних концентрацій — дискретні стани клітин. SmoothLife — узагальнення Life Конвея з розмитими сусідами, що близьке до реакційно-дифузійної динаміки.