Cellular Automata II: Gray-Scott Reaction-Diffusion Model

03-12-2023

1 Introduction
2 Background
3 Interactive walkthrough

Introduction

During the last term of my M.S. in Computer Science program at Johns Hopkins, I took an Independent Study course with Dr. Wiley focused on Cellular Automata (CA). We investigated a total of five CA models, and this blog post serves as the second in a series of several documentations of those projects. The first post was on Thomas Schelling's model of segregation, and you can find it here. In this post, I will present a CA implementation of a particularly interesting reaction-diffusion system known as the Gray-Scott model.

Background

Reaction-Diffusion systems

The Wikipedia page on reaction-diffusion systems has a good high-level explanation that I will reformulate in my own words. Reaction-diffusion systems are mathematical models (usually involving partial differential equations) that are used to model physical phenomena in chemistry, biology, geology, physics, and ecology. Perhaps the most famous example of these systems is the evolution of the concentration of one or more chemical substances across space and time. These changes involve chemical reactions between the substances and their diffusion over a surface in space.

One motivation for studying reaction-diffusion systems is the formation of patterns that can arise from certain models. Just like the segregation of populations that arose from the simple rules in Schelling's model of segregation, so too can complex and beautiful patterns arise from reaction-diffusion systems like Turing's patterns, the Belousov-Zhabotinsky (BZ) reaction, and the Gray-Scott model. The patterns that arise from the simulation of these models often resemble things like animal coats and skin pigmentation in nature, leading to the argument that reaction-diffusion processes might be the foundation for certain morphogenic processes in biology. In addition, we are motivated to study certain reaction-diffusion systems due to the trippy images that get generated.

The Gray-Scott model

The Gray-Scott model is based on a chemical reaction model developed by Peter Gray and Steven Scott in the 1980s [1]. The chemical reaction takes one molecule of a substance $u$ and two molecules of a substance $v$ and turns it into three molecules of $v$. The chemical equation is

$u + 2v \to 3v$

$u$ is continuously replenished from an external source at a feed rate $F$, while $v$ is removed from the system at a rate $k$ slightly faster than $F$. Then, the dynamics of the Gray-Scott model are given by the system of equations

$\frac{\partial u}{\partial t} = F(1 - u) - uv^2 + D_u \nabla^2 u$

$\frac{\partial v}{\partial t} = -(F+k)v + uv^2 + D_v \nabla^2 v$

The second term in the two equations above describe the aforementioned chemical reaction in that $u$ is consumed to produce $v$. The third term in the two equations describe the diffusion of the chemical substances, such that $D_u$ and $D_v$ describe the diffusion rates of $u$ and $v$.

Cellular Automata implementation

To simulate the Gray-Scott model using CA, we will use a two-dimensional grid to represent the concentrations of $u$ and $v$ at each timestep. Each chemical will have their own grid, but the grids interact with each other according to the differential equations. Each cell represents some arbitrary point in space, and the value at the cell, ranging from 0 to 1, determines the concentration of the substance at that location. The main dynamics we want to observe is the diffusion of the substances across the grid over time. We will use von Neumann neighborhoods and periodic boundary conditions for the CA model.

The state of each cell is updated by solving the differential equations at each timestep using numerical techniques. A simple strategy we adopt is the Euler method, where we calculate the value of a function at the next timestep by adding the value of the function at the current timestep to the derivative of the function times a step size. In other words, to calculate the concentration of a substance at the next timestep, we add the current concentration to the differential equation that describes the rate of change of the substance times a step size. For our CA model, the update rule can be conceptually split into two parts. In the first part, we determine the diffusion at each cell by considering the cell's four closest neighbors (von Neumann neighborhood), namely

$D[u(x, y, t+1)] = u(x + 1, y, t) + u(x, y + 1, t) + u(x - 1, y, t) + u(x, y - 1, t) - 4u(x, y, t)$

In the equation above, $x$ and $y$ are the coordinates of a substance $u$, and $t$ describes the current timestep. In the second part, the chemical reaction terms $F(1 - u) - uv^2$ and $-(F+k)v + uv^2$ from the differential equations are computed as

$R[u(x, y, t+1)] = -u(x, y, t)\cdot v(x,y,t)^2 + F\cdot (1 - u(x,y,t))$

$R[v(x, y, t+1)] = u(x, y, t)\cdot v(x,y,t)^2 - (F + k)\cdot v(x,y,t)$

Then, the update rules can be summarized as

$u(x,y,t+1)=u(x,y,t)+D_u \cdot D[u(x, y, t+1)] + R[u(x, y, t+1)]$

$v(x,y,t+1)=v(x,y,t)+D_v \cdot D[v(x, y, t+1)] + R[v(x, y, t+1)]$

One final note we make is that these simulations can often take a very long time to run due to the large number of cells and the calculations involved. As a result, we borrow from @gouarin's efficient implementation of the Gray-Scott model using numpy (we use numjs instead). In addition, readers interested in the Gray-Scott model can find further information on Robert Munafo's blog post.

Sample patterns

Here are some of my favorite patterns and the parameters that generated them. Use the arrow buttons below to cycle through 36 different patterns.

Pattern 1. F: 0.045, k: 0.06, Du: 0.16, Dv: 0.008

Interactive walkthrough

To run the simulation, first use the input boxes on the left to tweak the parameters. The first four input boxes correspond to the feed rate, kill rate, and diffusion rates of the system. Since it can sometimes be hard to find a combination of parameters that produce an interesting pattern, users can also insert preset parameters using the Sample parameters dropdown menu. The Frames and Steps/frame parameters can be used to tweak the smoothness of the animation. Since diffusion often happens slowly, users are encouraged to use a larger number of timesteps per frame. The Delay control determines the delay between each frame when playing the animation. Finally, the Board size slider controls the number of locations in the configuration.

The Chemical type dropdown can be used to display the chemical type (V is usually the more interesting one). The Neighborhood type dropdown can be used to specify the type of CA neighborhoods. Users can also use the Initial conditions dropdown to choose between several different initial conditions. Users can use the Low color and High color options to specify the colors that correspond to the concentration of the chemicals at each location. The Autoplay parameter specifies whether the animation automatically plays after being generated.

After selecting the parameters, click on the Generate button to generate the animation. Since the simulation can take a long time to complete, this demo works by first generating all the frames for the animation and then playing them afterwards for a smoother experience. After the animation is done generating, press the Play button to play the animation, or press the Prev or Next buttons to step through the frames. The Pause button can be used to pause the animation, and the Cancel button can be used to stop the generation of an animation. There are several randomization options on the third row of buttons that allow for rapid experimentation.

DESKTOP USERS: Drag your mouse across the board to "draw" your initial conditions.

Idle - Press 'Generate' to generate the animation frames

$F$:

$k$:

$D_u$:

$D_v$:

Frames:

Steps/frame:

Delay (0ms)

Board size (100x100)