Reaction-Diffusion

Reaction-diffusion is a group of models that can describe how chemicals react with one another in a grid. Among other things, they’ve been used in biology to describe how Turing patterns like zebra stripes and those on seashells form. I am intrigued with creating natural, organic-feeling shapes and patterns algorithmically, and reaction-diffusion can produce beautiful results in that regard. When I saw Zach Lieberman tweet several experiments with this algorithm, I became inspired to write my own implementation.

This wasn’t the first time I looked into this algorithm, but previously I was intimidated and discouraged by the math involved. But with a slightly different explanation of the equations involved, it turned out to actually be pretty straightforward. Because of this experience, I’d like to share a more code-oriented explanation here.

Writing our own

The first and best resource I found on reaction-diffusion is this tutorial, and I highly recommend reading it for a more conceptual explanation. The tutorial describes a system to simulate two chemicals on a 2D grid. The basic structure of the simulation is as follows:

• We keep a 2D grid of cells. For each cell, we store two values \(A\) and \(B\) to represent the chemicals.
• We initialize the values by setting \(A\) to 1 in every cell, and \(B\) to 0 in most cells.
• For each update step, we loop through every cell and calculate new values for \(A\) and \(B\) based on the values in the current cell and the cells around it.

In javascript, this could look something like this:

let dim = 256; // We will simulate on a square grid with dim*dim cells
let As = new Array(dim*dim).fill(1);
let Bs = new Array(dim*dim).fill(0);

function update() {
  let nextAs = new Array(dim*dim);
  let nextBs = new Array(dim*dim);

  for (let y = 0; y < dim; y++) {
    for (let x = 0; x < dim; x++) {
      let cursor = x + y*dim;

      nextAs[cursor] = /* some equation here */;
      nextBs[cursor] = /* some equation here */;
    }
  }

  As = nextAs;
  Bs = nextBs;
}

We will want to initialize some cells with a different value for \(B\) later, but let’s look at the update equations first. Karls tutorial describes them as:

\[ \begin{aligned} A' = A + (D_{A}\nabla^{2}A - AB^{2} + f (1 - A))\Delta t \\ B' = B + (D_{B}\nabla^{2}B + AB^{2} - (k + f)B)\Delta t \end{aligned} \]

I recommend reading the article to understand what each of the parameters and variables mean. For the purpose of implementing this simulation, we can think of the terms in these equations as follows:

• \(A\) and \(B\) are the current state of a cell, while \(A'\) and \(B'\) are the updated values.
• \(D_{A}\), \(D_{B}\), \(f\), \(k\) and \(\Delta t\) are global parameters. The values we choose changes the behavior of the simulation. We will use DA, DB, f and k as constants in our code, and disregard the \(\Delta t\) altogether.
• The \(AB^{2}\) represents the reaction in reaction-diffusion, while \(\nabla^{2}A\) and \(\nabla^{2}B\) describes the diffusion.

If you’re like me, a more conceptual explanation of the last term is helpful. If \(A\) in a cell is much higher than in the cells around it, it should spread (diffuse) into those cells. If \(A\) is higher in the surrounding cells, they should spread into the current cell. This is what \(\nabla^2A\) does. The tutoral suggests to take the difference between \(A\) of the current cell and the (weighted) average of the 8 cells around it around it. We will simplify it even more and only consider the four cells that are directly adjacent.

let diffuseA = (As[cursor - dim]       // above
              + As[cursor +   1]       // to the right
              + As[cursor + dim]       // below
              + As[cursor -   1]) / 4  // to the left
              - As[cursor];            // current cell

After doing the same for \(B\), we can now implement the equations in our update() function:

let A = As[cursor], B = Bs[cursor];
let reaction = A*B*B;

nextAs[cursor] = A + DA*diffusionA - reaction + (1 - A)*f;
nextBs[cursor] = B + DB*diffusionB + reaction - (k + f)*B;

That’s almost all there is to it. There are just two small details to work out before we have a working simulation. First, you might have realized that in calculating the diffusion values, we run into trouble for the cells on the edge of the grid. One way to get around this is to wrap our coordinate system. A slightly easier way is to never update cells at the edge of the grid:

for(let y = 0; y < dim; y++) {
  for(let x = 0; x < dim; x++) {
    let cursor = x + y*dim;

    if (x == 0 || y == 0 || x == dim - 1 || y == dim -1) {
      nextAs[cursor] = 1;
      nextBs[cursor] = 0;
      continue;
    }


    /* ... */
  }
}

Finally, we need to improve our initialization. We want to set \(B\) in at least a few cells to something other than zero. An easy thing to implement is to set the cells near the center to \(B = 1\).

// Right after initializing Bs
for (let y = dim/2 - 5; y <= dim/2 + 5; y++) {
  for (let x = dim/2 - 5; x <= dim/2 + 5; x++) {
    Bs[x + y*dim] = 1;
  }
}

Now we have everything to run a reaction-diffusion simulation. The following example all uses code we wrote for the simulation, the Canvas API for drawing (using setPixelData()), and requestAnimationFrame() to keep updating the simulation.

Closing words

We barely scratched the surface with this tiny simulation. Zach Liebermans experiments shows there are a lot of interesting ways to draw the results of the simulation. This interactive simulation proves that our implementation can be improved a lot in terms of speed, if we translate it into a WebGL shader. The simulation could also be extended to more dimensions and substances.

There are also some really creative applications of reaction-diffusion out there. n-e-r-v-o-u-s.com uses them for generative jewelry and other beautiful physical products, like the lovely cup below. Another application that I’d like to experiment with is using a camera to interact with the simulation, like the example on the right:

It is safe to say that I am not bored yet with this algorithm. I wrote a cool WebGL implementation for kaesve.nl, and I have many more ideas I’d like to explore. If you have any questions or comments, or if you found this article interesting and would like to read more, please email or tweet at me: mail@kaesve.nl / http://twitter.com/kaesve.