Introduction
========================================================================
Complex patterns in nature can sometimes emerge from surprisingly simple processes.
The modeling of such patterns was one of the last things Alan Turing studied before his much too early death.
His last paper, published in 1952, describes how chemical systems can give rise to various stable patters found in nature [#turing1952].
If we also consider chemicals that spread out over time we get a reaction-diffusion system. Many patterns found on the skin of various animals, such as the puffer fish in Figure [header], can be modeled in this manner.
![Figure [header]: Turing pattern observed on the back of a pufferfish at the Denmark National Aquarium outside of Copenhagen.](files/figures/copenhagen-aquarium.png width=100%)
A reaction-diffusion model capable of generating such patterns is the Gray-Scott model.
This model is capable of creating patterns like the one in Figure [header] and simulating
the growth of certain slime molds and corals.
A typical such pattern resembling the puffer fish skin in Figure [header] can
be seen in Figure [fishpattern].
The model consists of two chemicals that move through
diffusion and react in such a way that one _eats_ the other.
![Figure [fishpattern]: Typical Turing pattern simulated using the Gray-Scott model.](files/images/fish-turing-pattern-flat-000099.png width=100%)
In this blog post the Gray-Scott reaction-diffusion model is introduced and I show how it can be implemented in the python programming language.
The two parameters of the model are mapped out and the two main patterns that can be created using the model are shown.
These two patterns include the coral pattern seen in Figure [fishpattern] and a mitosis pattern.
!!! Info
_**Additional resources on the Gray-Scott reaction-diffusion model**_
I initially discovered this model through some posts by [Phillip Compeau](https://twitter.com/PhillipCompeau) on twitter advertising a [course that was later put into a book](https://biologicalmodeling.org/) that also covers some other interesting topics in biological modeling [#compeau2022].
Another great resource is the [reaction-diffusion website](https://www.karlsims.com/rd.html) by Karl Sims from which I recreated many of the figures in this blog post.
This site also has an absolutely wonderful [interactive reaction-diffusion simulation tool](https://www.karlsims.com/rdtool.html) that runs
as a GPU shader in the web browser. I highly recomend checking it out [#sims].
The work by Jason Webb at the [reaction-diffusion playground](https://jasonwebb.github.io/reaction-diffusion-playground/) looks very comprehensive
and showcases some nice examples created using an interactive tool.
Finally, I also found the web page of [Robert Munafo](https://mrob.com/pub/comp/xmorphia/) that contains a lot of implementation details for different older CPU architectures and a list of
publications for similar systems that generate different patterns.
I am sure that there are many more resources out there that cover this topic.
The Gray-Scott reaction-diffusion model
========================================================================
The Gray-Scott model was introduced in a series of papers in the 1980s [#gray1983] [#gray1984] [#gray1985] and has become popular, among others, in the computer graphics community as tool for pattern generation.
It consists of two chemical compounds, $\A$ and $\B$, that diffuse and react with each other.
Apart from diffusion three things happen in the system:
1. Chemical $\A$ is added to the system at a given feed rate, $f$. The rate decreses the higher the local amount of $\A$ is.
2. Two $\B$ can react with one $\A$ and produce an additional $\B$.
3. Chemical $\B$ is removed at a given kill rate, $k$.
These three processes are illustrated in Figure [ill0].
![Figure [ill0]: The Grey-Scott reaction-diffusion model contains a set of three update rules. Chemical $\A$ is added at a feed rate. Two $\B$ particles can *eat* an $\A$ particle and produce a $\B$, and lastly $\B$ particles are removed at a kill rate.](files/illustrations/reaction-illustration-stix.svg width=100%)
The chemicals disperse through diffusion where the rate, the diffusivity, of $\A$ is twice that of $\B$ which is illustrated in Figure [ill1].
![Figure [ill1]: Both the $\A$ and $\B$ particles diffuse so they spread across the domain.](files/illustrations/diffusion-rate-illustration-stix.svg width=100%)
The model is not simulated by tracking the individual particles but instead, as we can see in Figure [ill2], a grid is used and each grid cell records the local amount of $\A$ and $\B$.
![Figure [ill2]: The simulation is performed on a grid where each cell contains an amount of $\A$ and one of $\B$.](files/illustrations/diffusion-grid-illustration-stix.svg width=90%)
The model can be written as a system of two partial differential equations
\begin{equation}
\left\{
\begin{split}
& \frac{\partial \A}{\partial_{t}} = g_{A} \nabla^{2} \A + f(1-\A) - r \A \B^{2}\\
& \frac{\partial \B}{\partial_{t}} = g_{B} \nabla^{2} \B - k \B + r \A \B^{2}
\end{split}
\right.
\label{eq:diffeqdef}
\end{equation}
where the different variables represent the following:
\begin{equation*}
\left\{
\begin{aligned}
\A: & ~\text{Chemical}~\A. \\
\B: & ~\text{Chemical}~\B. \\
g_{\A}: & ~\text{Diffusivity of chemical}~\A. \\
g_{\B}: & ~\text{Diffusivity of chemical}~\B~\text{. Set to } g_{\B} = 0.5 \cdot g_{\A}. \\
f: & ~\text{Feed rate for chemical}~\A. \\
k: & ~\text{Kill rate for chemical}~\B. \\
r: & ~\text{Reaction rate. Set to 1.} \\
\end{aligned}
\right.
\end{equation*}
A dynamic simulation can then be implemented using a numerical scheme such as finite differences
both to discretize the temporal and spatial derivatives.
Simulating the reaction-diffusion system
========================================================================
The most straight forward solution to simulating this model, as described by the partial differential equations in Equation [eq:diffeqdef] and discretized on a regular grid as illustrated in Figure [ill2], is to use finite differences.
We discretize the temporal derivative using forward Euler with the time step $\tau$ as
\begin{equation}
\label{eq:tder}
\frac{\partial \A}{\partial t} \approx \frac{\An - \Ap}{\tau}
\end{equation}
with the same process being performed for $\B$.
The laplacian operators are discretized using a spatial kernel and both equations are solved for $\An$ and $\Bn$ respectively.
\begin{equation}
\left\{
\begin{split}
& \An = \Ap + \tau \cdot \Bigg( g_{A} \cdot \mathbf{D} \ast \Ap + f(1-\Ap) - \Ap {\Bp}^{2} \Bigg)\\
& \Bn = \Bp + \tau \cdot \Bigg( g_{B} \cdot \mathbf{D} \ast \Bp - k \Bp + \Ap {\Bp}^{2} \Bigg)
\end{split}
\right.
\label{eq:discrete}
\end{equation}
Where $\Dk$ is a two-dimensional discrete Laplacian kernel:
\begin{equation}
\Dk =
\begin{pmatrix}
0 & 1 & 0 \\
1 & -4 & 1 \\
0 & 1 & 0
\end{pmatrix}
\label{eq:laplaciankernel}
\end{equation}
and $\ast$ is the convolution operator.
I cover these simulation and discretization aspects somewhat more thouroughly in my blog posts
on [nonlinear](../image_smoothing_by_nonlinear_diffusion/) and [anisotropic](../image_smearing_by_anisotropic_diffusion/) diffusion modeling for image smoothing.
Now the grid can be set to some initial value and then simulated forward in time by updating the grid values of $\A$ and $\B$ according to Equation [eq:discrete].
!!! Tip
**Solution as cellular automata**
In the book [Biological Modeling](https://biologicalmodeling.org/) by
Phillip Compeau both the concept of diffusion and the Gray-Scott model
are introduced without using complicated mathematical concepts that might
deter readers with some other background.
Instead the model is explained through the random motions of particles
that cause them to diffuse and the reactions they have with each other.
The simulations are then implemented as a cellular automaton on a regular grid
that effectively performs the same thing as a finite difference solution
of the partial differential equations but without the mathematical jargon!
I prefer the differential equation way of presenting the model since by using it
we get a lot of tools for free and it is more generalizable.
However, what is actually important is the connection between these two ways of thinking!
We should not forget that diffusion represented by a differential equation is in of
itself a _macroscopic_ model of a _microscopic_ phenomenon.
Parameter selection using a style map
========================================================================
The range of values of feed and kill rates $f$ and $k$ that actually produces any interesting patterns is rather small and outside of this range no patterns will form or they will be highly unstable.
By creating a simulation that varies the two paramater values spatially across the domain as in Figure [parameter_map] we can observe what style of pattern various pairs of $f$ and $k$ will create.
![Figure [parameter_map]: A simulation where the feed ($f$) and kill ($k$) rates varies across the domain. It can be seen that only a narrow band of values creates any interesting patterns. Suitable parameter values for the coral pattern can be observed at **a)** and for the mitosis pattern at **b)**.](files/video/varying-feed-kill-rates-pattern-dr-04-1-bw-flat-matplotlib-grid_.mp4 width=75%)
In Figure [parameter_map] we can see how two stable patterns form:
1. At the parameter values $f \in [0.04, 0.06]$ and $k \in [0.060, 0.065]$ indicated by **a)** we see a _coral_ pattern form.
2. With values of $f \in [0.02, 0.04]$ and $k \in [0.060, 0.067]$ indicated by **b)** we instead see a _mitosis_ pattern. Named for its visual similarity to a cell division.
The style map can also be seen as an image [here](files/images/varying-feed-kill-rates-pattern-dr-04-1-bw-flat-matplotlib-grid-000399.png).
Python implementation using the GPU with pytorch
========================================================================
Since the simulation update is [embarrassingly parallel](https://en.wikipedia.org/wiki/Embarrassingly_parallel) and my initial implementation using numpy was a bit slow I tried
putting my GPU to good use with pytorch.
Theoretically the speedup should depend primarily on the number of cells, _i.e._ pixels,
in the grid, assuming the GPU can execute across almost all grid cells simultaneously while the CPU needs to loop over all the elements.
For smaller grid sizes such as $100 \times 100$ doing the computations on the CPU was faster
since the memory communication starts to dominate over the computation time.
The speedup I got using an Nvidia RTX 2070 GPU was approximately 2 for a grid of $500 \times 500$ cells, 25 for a grid of $1000 \times 1000$ and 40 for a grid of $2000 \times 2000$ cells. This compared to using an Intel CPU with six physical cores and 12 threads.
This was tested using pytorch version 1.13.1 and cuda 11.7 in Ubuntu 22.04.
Getting a reduction in compute time by a factor 40 essentially for free is something I consider pretty good!
**Gray-Scott reaction-diffusion model**
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Python
import numpy as np
import scipy
import torch
import torch.nn
def printCudaInfo() -> None:
print(f"PyTorch version: {torch.__version__}")
print(f"CUDA available: {torch.cuda.is_available()}")
print(f"CUDA version: {torch.version.cuda}")
print(f"CUDA device: {torch.cuda.get_device_name(0)}")
return
class ReactionDiffusionPdeTorch:
"""
Gray-Scott reaction diffusion model.
"""
def __init__(self):
self.simulation_time_s: float = 0.0
self.f: float|torch.Tensor = 0.050 # feed rate
self.k: float|torch.Tensor = 0.062 # kill rate
self.dA: float|torch.Tensor = 1.0 # prey diffusion constant
self.dB: float|torch.Tensor = 0.5 * self.dA # predator diffusion constant
self.dt: float = 0.125 # time constant
## 3 x 3 Laplacian matrix to calculate diffusion.
laplacian = np.zeros((1,1,3,3), dtype=np.float32)
## Karl Sims laplacian
laplacian[0,0,:,:] = np.array([ [0.05, 0.20 ,0.05],
[0.20,-1.00, 0.20],
[0.05, 0.20, 0.05]], dtype=np.float32)
## Ordinary FD laplacian (works fine)
laplacian[0,0,:,:] = np.array([ [0.00, 1.00 ,0.00],
[1.00,-4.00, 1.00],
[0.00, 1.00, 0.00]], dtype=np.float32)
self.laplacian: torch.Tensor = torch.tensor(laplacian)
self.A: torch.Tensor = torch.zeros(1) # Prey matrix
self.B: torch.Tensor = torch.zeros(1) # Predator matrix
self.G: Optional[torch.Tensor] = None # Diffusivity matrix (not used)
## Cuda
self.__device_force_cpu: bool = False
self.__device = torch.device('cpu')
self.setDevice(self.__device_force_cpu)
return
def setDevice(
self,
device_force_cpu: bool = False,
) -> None:
"""
"""
self.__device_force_cpu = device_force_cpu
self.__device = torch.device('cpu')
if not self.__device_force_cpu and torch.cuda.is_available():
self.__device = torch.device('cuda')
return
def getPopulationFraction(self) -> np.ndarray:
"""
Return a scalar image representing the pixel-wise relative
populations of predators by total population in that pixel.
This is used to visualize the population patterns.
"""
frac: torch.Tensor = self.B / (self.A + self.B)
frac = frac.cpu()
frac = frac[0,0,:,:]
return frac.numpy()
def initializeDefault(
self,
grid_rows: int=251,
grid_cols: int=251,
seed_size: int=11, # Needs to be an odd number
) -> None:
"""
Initialize with a square blob of B predators at center of domain.
"""
## World is full of prey and empty of predators
self.A = torch.ones ((1, 1, grid_rows, grid_cols), dtype=torch.float32, device=self.__device)
self.B = torch.zeros((1, 1, grid_rows, grid_cols), dtype=torch.float32, device=self.__device)
## Seed the predators in square at center
self.B[0, 0, int(grid_rows/2)-int(seed_size/2):int(grid_rows/2)+int(seed_size/2)+1, \
int(grid_cols/2)-int(seed_size/2):int(grid_cols/2)+int(seed_size/2)+1] = \
torch.ones((seed_size, seed_size), dtype=torch.float32, device=self.__device)
return
def initializeRandom(
self,
grid_rows: int=251,
grid_cols: int=251,
number_of_seeds: int=10,
smoothing_sigma: float = 0.0,
seed: int = -1
) -> None:
"""
Initialize with random smeared blobs.
"""
## World is full of prey and empty of predators
self.A = torch.ones ((1, 1, grid_rows, grid_cols), dtype=torch.float32, device=self.__device)
self.B = torch.zeros((1, 1, grid_rows, grid_cols), dtype=torch.float32, device=self.__device)
## Seed the predators in square at center
if seed>=0:
np.random.seed(seed)
indices_rows = np.random.choice(self.B.shape[2], number_of_seeds, replace=False)
indices_cols = np.random.choice(self.B.shape[3], number_of_seeds, replace=False)
self.B[0, 0, indices_rows, indices_cols] = 1.0
if smoothing_sigma > 0.0:
B = self.B.cpu().numpy()
B = scipy.ndimage.gaussian_filter(B, sigma=smoothing_sigma)
max_v = np.amax(B)
scale = 1.0 / max_v
B *= scale
self.B = torch.tensor(B, device=self.__device)
assert( (self.A.shape == self.B.shape) )
return
def gpu(self) -> None:
"""
Send all simulation data arrays to GPU.
"""
if not self.__device_force_cpu and torch.cuda.is_available():
self.A = self.A.cuda(self.__device)
self.B = self.B.cuda(self.__device)
self.laplacian = self.laplacian.cuda(self.__device)
if isinstance(self.f, torch.Tensor):
self.f = self.f.cuda(self.__device)
if isinstance(self.k, torch.Tensor):
self.k = self.k.cuda(self.__device)
if isinstance(self.dA, torch.Tensor):
self.dA = self.dA.cuda(self.__device)
if isinstance(self.dB, torch.Tensor):
self.dB = self.dB.cuda(self.__device)
else:
self.cpu()
return
def cpu(self) -> None:
"""
Retrieve all simulation data arrays to CPU.
"""
self.A = self.A.cpu()
self.B = self.B.cpu()
self.laplacian = self.laplacian.cpu()
if isinstance(self.f, torch.Tensor):
self.f = self.f.cpu()
if isinstance(self.k, torch.Tensor):
self.k = self.k.cpu()
if isinstance(self.dA, torch.Tensor):
self.dA = self.dA.cpu()
if isinstance(self.dB, torch.Tensor):
self.dB = self.dB.cpu()
return
def update(self) -> None:
"""
Update the simulation by one time step according to the
system om partial differential equations.
"""
## Compute the update of the Gray-Scott model for both the A and B populations
laplacian_A: torch.Tensor = torch.nn.functional.conv2d(input=self.A, weight=self.laplacian, stride=1, padding='same')
laplacian_B: torch.Tensor = torch.nn.functional.conv2d(input=self.B, weight=self.laplacian, stride=1, padding='same')
f: float = self.f
k: float = self.k
A_new = self.A + self.dt * (self.dA * laplacian_A - (self.A * self.B * self.B) + (f * (1-self.A)))
B_new = self.B + self.dt * (self.dB * laplacian_B + (self.A * self.B * self.B) - ((k+f) * self.B))
self.A = torch.clone(A_new)
self.B = torch.clone(B_new)
return
def simulate(
self,
iterations: int,
save_interval: int = 1,
verbose: bool = False,
)-> List[np.ndarray]:
"""
Simulate for a number of iterations and store the resulting
image used for visualization.
"""
images = []
t_total = 0.0
t0 = time.time()
self.gpu()
t1 = time.time()
t_total += t1 - t0
for i in range(iterations):
if verbose and i%(iterations//100)==0:
print(f"Iteration {i+1: 7d}/{iterations}")
if (i%save_interval) == 0:
images.append(self.getPopulationFraction())
t0 = time.time()
self.update()
t1 = time.time()
t_total += t1 - t0
self.simulation_time_s = t_total
if verbose:
print(f"Total simulation time: {t_total:20.4f} [s]")
return images
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Example: Animation
========================================================================
I performed several simulations to test both the reaction-diffusion model and my implementation with various parameter settings.
The Phong lighting model described in my [previous blog post](../heightfield_shading/) was used to give the patterns an embossed 3D look.
A typical coral pattern can be seen in Figure [animation_wide_coral_pattern] and a coral pattern
where a variable diffusivity modulates the scale of the pattern can be seen in Figure [animation_varying_diffusivity].
A pattern that includes both coral and mitosis parts which was produced by varying the feed and kill rates can be seen in Figure [animation_varying_f_k].
Finally a coral and a mitosis pattern can be seen separately in the two subfigures of Figure 10.
![Figure [animation_wide_coral_pattern]: Coral pattern simulation of the Gray-Scott model ($f=0.055, k=0.062$).](files/video/random-seeded-coral-pattern-shaded.mp4 width=75%)
![Figure [animation_varying_diffusivity]: Coral pattern simulation of the Gray-Scott model. The diffusivity varies from lower at the left to higher on the right. The diffusivity modulates the scale of the pattern. ($f=0.055, k=0.062$).](files/video/wide-varying-diffusivity-rates-pattern-shaded.mp4 width=75%)
![Figure [animation_varying_f_k]: Varying the feed and kill rates across the domain creates a coral pattern on the left and a mitosis pattern in the bottom right corner.](files/video/wide-varying-feed-kill-rates-pattern-shaded.mp4 width=75%)
a) Coral pattern.
($f=0.055, k=0.062$).
b) Mitosis pattern.
($f=0.034, k=0.064$).
**Figure 10:** _Gray-Scott model simulations that create a coral **a)** and a mitosis **b)** pattern._
Conclusion
========================================================================
The Gray-Scott reaction-diffusion model can be used to generate quite interesting looking patterns.
Picking the parameters $f$ and $k$ can be a bit tricky but is facilitated by mapping the patterns like in Figure [parameter_map].
Using a GPU with pytorch to solve the finite-difference systems in this case was
very straight forward and gave a substantial improvement in computation time.
I've always been very curios as to how these kinds of patterns could be created
so reading about and implementing this method was a very interesting and rewarding project.
!!! Info
_**Use for geoemtric design and other generative models**_
The Gray-Scott model has been used in the [design of various products](https://n-e-r-v-o-u-s.com/projects/albums/reaction-products/) such as lamp shades.
This is done by solving the system of differential equations on the surface of a 3D shape instead of on a regular grid.
A collection of links and resources that describe other pattern generating algorithms, cellular automata, fractals and more
can be found at [That Creative Code Page](https://available-anaconda-10d.notion.site/That-Creative-Code-Page-c5550ef2f7574126bdc77b09ed76651b)
and [Morphogenesis Resources](https://github.com/jasonwebb/morphogenesis-resources).
I also found the [GollyGang](https://github.com/GollyGang/ready/) tool that primarily uses cellular-automata
to solve various reaction-diffusion systems on different geometries other than regular grids.
+++++
**Bibliography**:
[#turing1952]: Alan Turing. 1952.
**_The Chemical Basis of Morphogenesis_**
Philosophical Transactions of the Royal Society of London B. 237 (641): 37–72
[DOI: 10.1098/rstb.1952.0012](https://doi.org/10.1098/rstb.1952.0012)
[PDF](https://hopf.chem.brandeis.edu/members_content/yanglingfa/pattern/Turing/The%20Chemical%20Basis%20of%20Morphogenesis.pdf)
[#gray1983]: P. Gray and SK. Scott. 1983.
**_Autocatalytic reactions in the isothermal, continuous stirred tank reactor: Isolas and other forms of multistability_**
Chemical Engineering Science. 38 (1): 29-43
[DOI: 10.1016/0009-2509(83)80132-8]("https://doi.org/10.1016/0009-2509(83)80132-8")
[Publisher link](https://www.sciencedirect.com/science/article/pii/0009250983801328)
[#gray1984]: P. Gray and SK. Scott. 1984.
**_Autocatalytic reactions in the isothermal, continuous stirred tank reactor: Oscillations and instabilities in the system A + 2B → 3B; B → C_**
Chemical Engineering Science. 39 (6): 1087-1097
[DOI: 10.1016/0009-2509(84)87017-7]("https://doi.org/10.1016/0009-2509(84)87017-7")
[Publisher link](https://www.sciencedirect.com/science/article/pii/0009250984870177)
[#gray1985]: P. Gray and SK. Scott. 1985.
**_Sustained oscillations and other exotic patterns of behavior in isothermal reactions_**
The Journal of Physical Chemistry. 89 (1): 22-32
[#compeau2022]: Phillip Compeau. 2022.
**_Biological modeling_**
https://biologicalmodeling.org/prologue/
+++++
**Links**:
**Karl Sims**
https://www.karlsims.com/
**Karl Sims - Reaction-Diffusion Tutorial**
https://www.karlsims.com/rd.html
**Karl Sims - Reaction-Diffusion interactive tool**
https://www.karlsims.com/rdtool.html
**Jason Webb**
https://github.com/jasonwebb
**Jason Webb - Reaction-Diffusion Playground**
https://github.com/jasonwebb/reaction-diffusion-playground
**Jason Webb - Morphogenesis resources**
https://github.com/jasonwebb/morphogenesis-resources
**Robert Munafo - Reaction-Diffusion**
http://mrob.com/pub/comp/xmorphia/
**R. de Courville and T. Muhonen - That Creative Code Page**
https://available-anaconda-10d.notion.site/That-Creative-Code-Page-c5550ef2f7574126bdc77b09ed76651b
**Embarrassingly parallel**
https://en.wikipedia.org/wiki/Embarrassingly_parallel
**GollyGang cellular automata**
https://github.com/GollyGang/ready/
**Reaction-Diffusion for geometric design**
https://n-e-r-v-o-u-s.com/projects/albums/reaction-products/
Code for all figures and examples:
https://bitbucket.org/nils_olovsson/reaction_diffusion_and_turing_patterns