Reaction diffusion is a mathematical model typically used to simulate the behavior of two chemicals (A&B) reacting to one another. Through this reaction they transform from one to the other based on a feed and kill rate. The two chemicals are diffused or spread through space at a rate defined by a diffusion rate for both chemicals A and B. Based on the changes of these rates a reaction diffusion model will create a wide range of emergent behaviors and patterns that are also found in nature. There are a number of great resources online. Some of the basic math structure is shown below. The model used in the processing sketch above is the Gray Scott model. The code for the processing sketch abive can be found at the bottom of this page.

There is a great explanation of the Gray Scott model in Karl Sims Reaction Diffusion Tutorial. The two dimensional space is constantly updated using the equation above as explained on Sims’ tutorial.

The images above from Sims’ tutorial show the various patterns that emergent patterns when the feed, kill, and diffusion rates change throughout the two dimensional space.

Not all feed and kill rates produce behaviors that are stable. The graphs above show the range of feed and kill rates that produce stable patterns. The image on the left is a graph from Robert Munafo’s detailed explanation of Gray-Scott reaction-diffusion. The youtube video to the right is a great graph showing the actual dynamic behaviors in the stable zones.

The processing sketch below generates a lower resolution spatial model that is scaled up using the variable “space.” The pixels in between are interpolated. This helps speed up the simulation.

``````import controlP5.*;
ControlP5 cp5;

int W = 280;
int H = 160;

//System parameters
float diffU;
float diffV;
float paramF;
float paramK;

boolean rndInitCondition = false;

float[][] U = new float[W][H];
float[][] V = new float[W][H];
float[][] cH = new float[W][H];
float[][] cB = new float[W][H];

float[][] dU = new float[W][H];
float[][] dV = new float[W][H];

int[][] offset = new int[W*H][4];

int space = 3;

void setup() {
size(1040,680);
frameRate(25);
smooth();
colorMode(HSB,1.0);

cp5 = new ControlP5(this);

.setPosition(100,20)
.setRange(0,0.1)
.setSize(240,20)
.setValue(0.035)
.setColorActive(color(0.0,0.0,0.6))
.setColorBackground(color(0.0,0.0,0.2))
.setColorForeground(color(0.0,0.0,0.4))
;

.setPosition(100,50)
.setRange(0,0.1)
.setSize(240,20)
.setValue(0.06)
.setColorActive(color(0.0,0.0,0.6))
.setColorBackground(color(0.0,0.0,0.2))
.setColorForeground(color(0.0,0.0,0.4))
;

.setPosition(400,20)
.setRange(0,0.2)
.setSize(240,20)
.setValue(0.16)
.setColorActive(color(0.0,0.0,0.6))
.setColorBackground(color(0.0,0.0,0.2))
.setColorForeground(color(0.0,0.0,0.4))
;

.setPosition(400,50)
.setRange(0,0.2)
.setSize(240,20)
.setValue(0.08)
.setColorActive(color(0.0,0.0,0.6))
.setColorBackground(color(0.0,0.0,0.2))
.setColorForeground(color(0.0,0.0,0.4))
;

.setPosition(765,50)
.setSize(50,20)
.setColorActive(color(0.0,0.0,0.6))
.setColorBackground(color(0.0,0.0,0.2))
.setColorForeground(color(0.0,0.0,0.4))
;
.setPosition(825,20)
.setSize(50,20)
.setColorActive(color(0.0,0.0,0.6))
.setColorBackground(color(0.0,0.0,0.2))
.setColorForeground(color(0.0,0.0,0.4))
;

.setPosition(825,50)
.setSize(50,20)
.setColorActive(color(0.0,0.0,0.6))
.setColorBackground(color(0.0,0.0,0.2))
.setColorForeground(color(0.0,0.0,0.4))
;

.setPosition(885,20)
.setSize(50,20)
.setColorActive(color(0.0,0.0,0.6))
.setColorBackground(color(0.0,0.0,0.2))
.setColorForeground(color(0.0,0.0,0.4))
;

.setPosition(885,50)
.setSize(50,20)
.setColorActive(color(0.0,0.0,0.6))
.setColorBackground(color(0.0,0.0,0.2))
.setColorForeground(color(0.0,0.0,0.4))
;

.setPosition(100,600)
.setSize(50,20)
.setColorActive(color(0.0,0.0,0.6))
.setColorBackground(color(0.0,0.0,0.2))
.setColorForeground(color(0.0,0.0,0.4))
;

.setPosition(765,20)
.setSize(50,20)
.setColorActive(color(0.0,0.0,0.6))
.setColorBackground(color(0.0,0.0,0.2))
.setColorForeground(color(0.0,0.0,0.4))
;

//Set default parameters;
diffU = 0.16;
diffV = 0.08;
paramF = 0.035;
paramK = 0.06;

rndInitCondition = true;

//Populate U and V with initial data
generateInitialState();

//Set up offsets
int c = 0;
for (int i = 0; i < W; i++) {
for (int j = 0; j < H; j++) {
if(i > 0){
offset[c][0] = i - 1;
}else{
offset[c][0] = W-1;
}

if(i < W-1){
offset[c][1] = i + 1;
}else{
offset[c][1] = 0;
}
if(j > 0){
offset[c][2] = j - 1;
}else{
offset[c][2] = H-1;
}
if(j < H-1){
offset[c][3] = j + 1;
}else{
offset[c][3] = 0;
}
c++;
}
}
}

void draw(){
background(color(0.0,0.0,0.8));
for (int k = 0; k < 20; k++) {
timestep(paramF, paramK, diffU, diffV);
}

// Draw points
noStroke();
for (int i = 0; i < W; i++) {
for (int j = 0; j < H; j++) {
cH[i][j] = (float)(1.3-U[i][j]*1.3);
cB[i][j] = (float)(1-U[i][j]);

color c = color((float)(1.3-U[i][j]*1.3), 0.87, (float)(1-U[i][j]));
set(i*space+100,j*space+100,c);
}
}

for (int i = 0; i < W*space-space; i++) {
for (int j = 0; j < H*space-space; j++) {
if((i%space != 0) || (j%space !=0)){
float x = i;
float x1 = int(i/space)*space;
float x2 = x1 + space;

float y = j;
float y1 = int(j/space)*space;
float y2 = y1 + space;

float cH1 = (x2 - x)/(x2 - x1)*cH[int(i/space)][int(j/space)] + ((x - x1)/(x2 - x1))*cH[int(i/space)+1][int(j/space)] ;
float cH2 = ((x2 - x)/(x2 - x1))*cH[int(i/space)][int(j/space)+1] + ((x - x1)/(x2 - x1))*cH[int(i/space)+1][int(j/space)+1];

float cHH = ((y2 - y)/(y2 - y1))*cH1 + ((y - y1)/(y2 - y1))*cH2;

float cB1 = (x2 - x)/(x2 - x1)*cB[int(i/space)][int(j/space)] + ((x - x1)/(x2 - x1))*cB[int(i/space)+1][int(j/space)] ;
float cB2 = ((x2 - x)/(x2 - x1))*cB[int(i/space)][int(j/space)+1] + ((x - x1)/(x2 - x1))*cB[int(i/space)+1][int(j/space)+1];

float cBB = ((y2 - y)/(y2 - y1))*cB1 + ((y - y1)/(y2 - y1))*cB2;
color c = color(cHH,0.87,cBB);
set(i+100,j+100,c);
}
}
}
}

void timestep(float F, float K, float diffU, float diffV) {
int c = 0;
for (int i = 0; i < W; i++) {
for (int j = 0; j < H; j++) {

float u = U[i][j];
float v = V[i][j];

int left = offset[c][0];
int right = offset[c][1];
int up = offset[c][2];
int down = offset[c][3];

float uvv = u*v*v;

float lapU = (U[left][j] + U[right][j] + U[i][up] + U[i][down] - 4*u);
float lapV = (V[left][j] + V[right][j] + V[i][up] + V[i][down] - 4*v);

dU[i][j] = diffU*lapU  - uvv + F*(1 - u);
dV[i][j] = diffV*lapV + uvv - (K+F)*v;
c++;
}
}

for (int i= 0; i < W; i++) {
for (int j = 0; j < H; j++){
U[i][j] += dU[i][j];
V[i][j] += dV[i][j];
}
}
}

void generateInitialState() {
for (int i = 0; i < W; i++) {
for (int j = 0; j < H; j++) {
U[i][j] = 1.0;
V[i][j] = 0.0;
cH[i][j] = (float)(1.3-U[i][j]*1.3);
cB[i][j] = (float)(1-U[i][j]);
}
}
}

void state_1(){
cp5.getController("diffU").setValue(0.16);
cp5.getController("diffV").setValue(0.08);
cp5.getController("paramF").setValue(0.035);
cp5.getController("paramK").setValue(0.06);
}

void state_2(){
cp5.getController("diffU").setValue(0.16);
cp5.getController("diffV").setValue(0.08);
cp5.getController("paramF").setValue(0.042);
cp5.getController("paramK").setValue(0.065);
}

void state_3(){
cp5.getController("diffU").setValue(0.18);
cp5.getController("diffV").setValue(0.13);
cp5.getController("paramF").setValue(0.025);
cp5.getController("paramK").setValue(0.056);
}

void state_4(){
cp5.getController("diffU").setValue(0.18);
cp5.getController("diffV").setValue(0.09);
cp5.getController("paramF").setValue(0.02);
cp5.getController("paramK").setValue(0.056);
}

void state_5(){
cp5.getController("diffU").setValue(0.14);
cp5.getController("diffV").setValue(0.06);
cp5.getController("paramF").setValue(0.035);
cp5.getController("paramK").setValue(0.065);
}

void state_6(){
cp5.getController("diffU").setValue(0.16);
cp5.getController("diffV").setValue(0.08);
cp5.getController("paramF").setValue(0.05);
cp5.getController("paramK").setValue(0.065);
}

void reset(){
generateInitialState();
}

void mouseDragged() {
for (int i = 0; i < W; i++) {
for (int j = 0; j < H; j++) {
if(dist(i,j,(mouseX-100)/space,(mouseY-100)/space) < 4){
U[i][j] = 0.5*(1 + random(-1, 1));
V[i][j] = 0.25*(1 + random(-1, 1));
}
}
}
}

void mousePressed() {
for (int i = 0; i < W; i++) {
for (int j = 0; j < H; j++) {
if(dist(i,j,(mouseX-100)/space,(mouseY-100)/space) < 4){
U[i][j] = 0.5*(1 + random(-1, 1));
V[i][j] = 0.25*(1 + random(-1, 1));
}
}
}
}``````