Creating a Reaction-Diffusion Simulation

Alexander Schoch,Fri Aug 16 2024•math python english

In this article, we'll look into how to create this animation (with a slightly lower resolution):

This animation shows a reaction-diffusion system, a system that contains multiple components, such as chemicals or bacteria, that can react with each other and diffuse over the spacial domain. For the simulation in this article, we'll look specifically at the Gray-Scott model.

The Gray-Scott model describes a system with the following rules:

The system knows two components, $A$ and $B$ , with their concentrations $a(x, y, t)$ and $b(x, y, t)$ at each pixel,
The component $A$ is constantly fed into the system,
The compoment $B$ is constantly removed from the system,
One particle $A$ and two particles $B$ react to from three particles $B$ ,
Both components diffuse with different rates over time.

This system can be expressed using the following set of partial differential equation (Source (opens in a new tab)):

\frac{\partial a}{\partial t} = D_\text{A} \nabla^2 a - ab^2 + F(1 - a)

\frac{\partial b}{\partial t} = D_\text{B} \nabla^2 b + ab^2 - (F + k) b

Let's break these equations down and look at their individual parts:

$\partial a / \partial t$ is the change of $A$ 's concentration over time.
$D_\text{A} \partial^2 a / \partial x^2$ is the diffusion of $A$ and corresponds to Fick's second law (opens in a new tab) (analoguous to the heat equation).
$ab^2$ is the reaction, given as a rate equation (opens in a new tab) (we just multiply the concentrations with the stoichiometric factor as an exponent). Here, $A$ is consumed (hence the $-$ ) and $B$ is produced (hence the $+$ ).
$F(1-a)$ is the feed of $A$ . The $(1-a)$ factor makes sure that the concentration of $A$ does not exceed 1 by making this term smaller for larger values of $a$ .
$(F+k)b$ is the kill of $B$ . The $(F+k)$ term makes sure that the kill of $B$ is always larger than the feed of $A$ , and the $b$ scaling factor increases the kill for higher values of $b$ .

This equation can be discretized for our purposes (Source (opens in a new tab)):

\Delta a = (D_\text{A} \nabla^2 a - ab^2 + F(1-a)) \Delta t

\Delta b = (D_\text{B} \nabla^2 b + ab^2 - (F+k)b) \Delta t

where $\Delta a$ is the change in concentration over one time step. $\nabla^2$ (Laplace operator, second derivative in space) can be calculated using a 3x3 Gaussian kernel (also called convolution or blur). We'll cover this later on.

Now, let's look at how to implement this in python.

The Code

First, we create our virtual environment and install all required packages:

python -m venv venv
source venv/bin/activate
 
pip install numpy matplotlib alive_progress

Next, we import these packages into our newly created python file:

import numpy as np
import matplotlib.pyplot as plt
import os
 
from alive_progress import alive_bar

Define the parameters that control the result of the simulation. I've taken the default ones from Karlsims (opens in a new tab) here. R_init describes the radius of the cirlce in pixel that we'll use as an initial seed. num_frames defines how many frames we want to render out. Here, 1200 frames with a frame rate of 60 FPS results in 20 seconds of animation.

N = 400           # Number of pixels in both dimensions
T = 5000          # Number of time steps
DA = 1            # Diffusivity of A
DB = 0.5          # Diffusivity of B
f = 0.055         # Feed rate of A
k = 0.062         # Kill rate of B
dt = 1            # Time Step
R_init = 5        # Radius of initial condition circle
num_frames = 1200 # Number of frames to be generated

Now, we calculate a few parameters we will use for image generation. frame_interval is an integer that stores the interval for image generation (i.e. "we want to output an image every 12th frame"). // is the python integer division that will make sure that frame_interval is an integer. num_digits stores the maximum number of digits the frame number can hold.

frame_interval = T // num_frames
frame_nr = 0
num_digits = int(np.log10(num_frames)) + 1 # For filename generation

For our initial condition, we define two matrices A and B, where A is filled with 1 and B with 0. We then create a mask that is True whenever that pixel is less than R_init away from the center of the image and False otherwise. We then set the concentration of $B$ to 1 in this circle.

A = np.ones((N, N))
B = np.zeros((N, N))
X, Y = np.meshgrid(range(N), range(N))
mask = np.sqrt((X - N//2)**2 + (Y - N//2)**2) < R_init
B[mask] = 1

A short demonstration on how this works, exactly:

In [1]: import numpy as np
 
In [2]: X, Y = np.meshgrid(range(5), range(5))
 
In [3]: X
Out[3]: 
array([[0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4]])
 
In [4]: Y
Out[4]: 
array([[0, 0, 0, 0, 0],
       [1, 1, 1, 1, 1],
       [2, 2, 2, 2, 2],
       [3, 3, 3, 3, 3],
       [4, 4, 4, 4, 4]])
 
In [5]: N = 5
 
In [6]: np.sqrt((X - N//2)**2 + (Y - N//2)**2) < 2
Out[6]: 
array([[False, False, False, False, False],
       [False,  True,  True,  True, False],
       [False,  True,  True,  True, False],
       [False,  True,  True,  True, False],
       [False, False, False, False, False]])

In order to make sure that this does what we want it to do, let's quickly plot it.

plt.imshow(B)
plt.colorbar()
plt.show()

Initial Seed

Now, we need to implement the mathematical expression described above. Note that it is important to first calculate dA and dB and to add them to A and B afterwards. Otherwise, we would already update A when we still need the old values for the calculation of B.

t % frame_interval == 0 checks if the current frame number is a multiple of frame_interval. If so, save an image of the current concentration, calculated as $B / (A + B)$ . This gives a number between 0 and 1 for each pixel.

convolution and save_image are functions we will write below.

with alive_bar(T) as bar: # progress bar
    for t in range(T):
        dA = DA * convolution(A) - A * B * B + f * (1 - A)
        dB = DB * convolution(B) + A * B * B - (k + f) * B
        A += dA * dt
        B += dB * dt
 
        if(t % frame_interval == 0):
            save_image(B / (A + B), frame_nr)
            frame_nr += 1
 
        bar()

Convolution

The 3x3 Gaussian kernel mentioned above basically describes the weighted average of all cells within a 3x3 square around any pixel. Here, the specific weights are

-1 for the pixel itself,
0.2 for the four direct neighbors, and
0.05 for the four diagonal pixels.

For example, the center pixel of the array

[[ 1, 2, 3],
 [ 4, 5, 6],
 [ 7, 8, 9]]

has a convolution of

\text{conv}(\text{middle}) = -5 + 0.2(2 + 4 + 6 + 8) + 0.05(1 + 3 + 7 + 9) = 0

Before implementing this, there is one more thing to think about: What do we do at the bounary of the domain? At the leftmost edge, there is no left neighbor, so we need to think about what to do in this case. Two simple options include

ignoring the boundary cells and just applying the algorithm to all non-border pixels
looping the domain around, such that the left neighbor of a left edge cell is a right border cell

In this article, we'll use the latter, as a very nice function exists for this exact purpose: np.roll (opens in a new tab). This function will take an array, a shift and an axis as inputs and will loop the values around by this shift in the axis dimension.

In [1]: X
Out[1]: 
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])
 
In [2]: np.roll(X, 1, axis=0) # Roll vertically
Out[2]: 
array([[7, 8, 9],
       [1, 2, 3],
       [4, 5, 6]])
 
In [3]: np.roll(X, 1, axis=1) # Roll horizontally
Out[3]: 
array([[3, 1, 2],
       [6, 4, 5],
       [9, 7, 8]])

Using this, we can now write our convolution function. Note that this function needs to be defined above the main loop.

def convolution(X):
    return -X
    + 0.2 * (
        np.roll(X, 1, axis=0) 
        + np.roll(X, 1, axis=1) 
        + np.roll(X, -1, axis=0) 
        + np.roll(X, -1, axis=1)) 
    + 0.05 * (
        np.roll(X, (1, 1), axis=(0, 1)) 
        + np.roll(X, (-1, 1), axis=(0, 1)) 
        + np.roll(X, (1, -1), axis=(0, 1)) 
        + np.roll(X, (-1, -1), axis=(0, 1)))

Video Generation

Lastly, the save_image function plots the concentration of $B$ using matplotlib.pyplot.imshow with a nice color map and saves it into e.g. frames/0001.png. The top four lines of the function remove the white margin around the plot and set a 1:1 aspect ratio. This function needs to be defined above the main loop as well. Also, make sure that the frames folder exists!

def save_image(X, i):
    fig = plt.figure(figsize=(1, 1))
    ax = plt.Axes(fig, [0., 0., 1., 1.])
    ax.set_axis_off()
    fig.add_axes(ax)
 
    ax.imshow(X, cmap='inferno')
    filename = "%s.png" % str(i).zfill(num_digits)
    plt.savefig(os.path.join('frames', filename), dpi=N)
    plt.close('all')

After running the script and generating all frames, we can combine them into a single video using ffmpeg.

ffmpeg -r 60 -pattern_type glob -i 'frames/????.png' -c:v libx264 video.mp4

That's it!

Varying the parameters $D_\text{A}$ , $D_\text{B}$ , $F$ and $k$ allows for a very wide range of patterns. It is also possible to change parameters over time or space to get really cool animations!

Please feel free to share your finished animation with me! You can send it via email (schochal[at]ethz.ch) or post it on Mastodon and tag me (@alexander_schoch@linuxrocks.online). Also, I'd love to hear your feedback about this post!

Full Code

import numpy as np
import matplotlib.pyplot as plt
import os
 
from alive_progress import alive_bar
 
N = 400
T = 20000
DA = 1
DB = 0.5
f = 0.055
k = 0.062
dt = 1
R_init = 5
num_frames = 1200
 
frame_interval = T // num_frames
frame_nr = 0
num_digits = int(np.log10(num_frames)) + 1 # For filename generation
 
A = np.ones((N, N))
B = np.zeros((N, N))
X, Y = np.meshgrid(range(N), range(N))
mask = np.sqrt((X - N//2)**2 + (Y - N//2)**2) < R_init
B[mask] = 1
 
def convolution(X):
    return -X
    + 0.2 * (
        np.roll(X, 1, axis=0) 
        + np.roll(X, 1, axis=1) 
        + np.roll(X, -1, axis=0) 
        + np.roll(X, -1, axis=1)) 
    + 0.05 * (
        np.roll(X, (1, 1), axis=(0, 1)) 
        + np.roll(X, (-1, 1), axis=(0, 1)) 
        + np.roll(X, (1, -1), axis=(0, 1)) 
        + np.roll(X, (-1, -1), axis=(0, 1))) 
 
def save_image(X, i):
    fig = plt.figure(figsize=(1, 1))
    ax = plt.Axes(fig, [0., 0., 1., 1.])
    ax.set_axis_off()
    fig.add_axes(ax)
 
    ax.imshow(X, cmap='inferno')
    filename = "%s.png" % str(i).zfill(num_digits)
    plt.savefig(os.path.join('frames', filename), dpi=N)
    plt.close('all')
 
 
with alive_bar(T) as bar:
    for t in range(T):
        dA = DA * convolution(A) - A * B * B + f * (1 - A)
        dB = DB * convolution(B) + A * B * B - (k + f) * B
        A += dA * dt
        B += dB * dt
 
        if(t % frame_interval == 0):
            save_image(B / (A + B), frame_nr)
            frame_nr += 1
 
        bar()

Comments

Subscribe via RSS.