Reaction Diffusion

3 min read•1 page

In 1952, Alan Turing (the father of modern computing) wrote a paper explaining how leopards get their spots and zebras get their stripes. He proposed that two invisible chemicals spread through the skin and react with each other to create the patterns.

Chemical A and Chemical B:

1. Imagine a grid filled with Chemical A (represented as white). 2. We drop a tiny seed of Chemical B (represented as black) into the center. 3. The Gray-Scott Equations define how they interact. 4. Chemical B "eats" Chemical A to replicate itself. But if Chemical B doesn't eat fast enough, it "dies" and vanishes.

python

1# Reaction-Diffusion Chemistry
2def gray_scott_equation(A, B, feed, kill):
3    """How the two chemicals interact."""
4    
5    # 1. Chemical A is added to the system at a constant 'feed' rate.
6    # 2. Chemical B is removed from the system at a constant 'kill' rate.
7    
8    # 3. The Reaction:
9    # 2 parts of B react with 1 part of A to create 3 parts of B!
10    # (A + 2B -> 3B)
11    
12    reaction_rate = A * B * B
13    
14    # 4. Change in A = Feed - Reaction
15    delta_A = feed * (1 - A) - reaction_rate
16    
17    # 5. Change in B = Reaction - Kill
18    delta_B = reaction_rate - (kill + feed) * B
19    
20    return delta_A, delta_B

Chemical State (Only A / Seed B / Reaction)

0.00

2 min read•1 page

Cellular Automata (like the Game of Life) are digital: a cell is either 1 or 0. Reaction-Diffusion is analog: it simulates continuous fluid dynamics using floating-point numbers.

Decimals vs Integers:

1. We use a 2D Grid, just like before. 2. But a cell doesn't just hold a state of Alive or Dead. 3. Every single cell holds two exact decimal numbers between 0.0 and 1.0. 4. The first number represents the concentration of Chemical A. The second represents Chemical B. 5. Because these are decimals, we can create incredibly smooth, organic gradients instead of blocky pixels.

python

1# Floats, not Booleans
2def create_concentration_grid(width, height):
3    """Setting up the petri dish."""
4    
5    grid = []
6    
7    for y in range(height):
8        row = []
9        for x in range(width):
10            
11            # Instead of a boolean (Dead/Alive)
12            # A cell stores two decimal numbers!
13            cell_state = {
14                "A": 1.0, # 100% full of A
15                "B": 0.0  # 0% full of B
16            }
17            
18            row.append(cell_state)
19            
20        grid.append(row)
21        
22    return grid

Concentration of B (0.0 -> 1.0)

0.00

3 min read•1 page

If we drop ink into a glass of water, it slowly bleeds outward until the entire glass is tinted. This is called Diffusion. In our grid, chemicals bleed into their neighbors using a Laplacian Convolution Matrix.

Weighted Neighborhoods:

1. Just like Conway, a cell looks at its 8 Moore Neighbors. 2. But instead of just counting them, it performs a weighted sum. 3. The center cell subtracts its own value (-1.0). 4. The 4 direct neighbors (Top/Bottom/Left/Right) contribute 0.2 each. 5. The 4 diagonal neighbors are slightly further away, so they only contribute 0.05 each. 6. This mathematical matrix (a 3x3 Convolution Kernel) perfectly simulates liquid diffusion!

python

1# 2D Spatial Diffusion
2def calculate_laplacian(grid, x, y, chemical):
3    """How the chemical bleeds into neighbors."""
4    
5    # 1. Start with the center cell (multiplied by -1)
6    sum_val = grid[y][x][chemical] * -1.0
7    
8    # 2. Add the direct neighbors (Top, Bottom, Left, Right)
9    # They get a weight of 0.2
10    sum_val += grid[y-1][x][chemical] * 0.2
11    sum_val += grid[y+1][x][chemical] * 0.2
12    sum_val += grid[y][x-1][chemical] * 0.2
13    sum_val += grid[y][x+1][chemical] * 0.2
14    
15    # 3. Add the diagonal neighbors
16    # They get a smaller weight of 0.05 because they are further away
17    sum_val += grid[y-1][x-1][chemical] * 0.05
18    sum_val += grid[y-1][x+1][chemical] * 0.05
19    sum_val += grid[y+1][x-1][chemical] * 0.05
20    sum_val += grid[y+1][x+1][chemical] * 0.05
21    
22    return sum_val

Weights (Center / Direct / Diagonal)

0.00

3 min read•1 page

We now have the two halves of our physics simulation: The Reaction (how the chemicals eat each other) and the Diffusion (how the chemicals bleed across the grid).

The Full Equation:

1. Just like Conway, we use two separate grids (a Read Buffer and a Write Buffer). 2. We loop through every single cell in the Read Buffer. 3. We calculate the Diffusion (Laplacian Matrix) for that cell. 4. We calculate the Reaction (Gray-Scott Equation) for that cell. 5. We add those two numbers together, multiply by a Time Step, and add it to the cell's current value. 6. We save the resulting decimal number to the Write Buffer!

python

1# Calculating the Next Frame
2def update_cell_concentration(grid, x, y, feed, kill, diff_A, diff_B, time_step):
3    
4    # 1. Look at the current amount of chemicals in this specific cell
5    A = grid[y][x]["A"]
6    B = grid[y][x]["B"]
7    
8    # 2. How much liquid is bleeding in/out from the neighbors?
9    laplace_A = calculate_laplacian(grid, x, y, "A")
10    laplace_B = calculate_laplacian(grid, x, y, "B")
11    
12    # 3. How much are the chemicals reacting with each other?
13    delta_A, delta_B = gray_scott_equation(A, B, feed, kill)
14    
15    # 4. Add it all together to get the new state!
16    # Formula: New = Old + (Diffusion + Reaction) * TimeStep
17    
18    new_A = A + (diff_A * laplace_A + delta_A) * time_step
19    new_B = B + (diff_B * laplace_B + delta_B) * time_step
20    
21    return new_A, new_B

Process (Reaction / Diffusion / Update)

0.00

2 min read•1 page

When calculating the Laplacian Matrix, a cell needs to look at its Top, Bottom, Left, and Right neighbors. But what happens if the cell is on the absolute edge of the grid? It doesn't have a Left neighbor!

Toroidal Wrapping:

1. If we just ignore the missing neighbors, the chemicals will hit the edge of the grid and behave strangely, creating ugly square borders. 2. Instead, we use Toroidal Boundary Conditions. 3. If a cell on the far Left edge looks for its left neighbor, we wrap around and give it the cell on the far Right edge. 4. Mathematically, this folds our flat 2D grid into a 3D Donut (Torus). The simulation becomes infinite and seamless!

python

1# Wrapping the Edges
2def wrap_coordinate(value, max_size):
3    """Making the grid infinite (Toroidal)."""
4    
5    # If we go off the right edge (value >= max_size)
6    # The modulo operator (%) wraps us back to 0!
7    
8    # If we go off the left edge (value < 0)
9    # Modulo wraps us back to max_size - 1!
10    
11    return value % max_size
12
13# Usage inside the Laplacian function:
14# ny = wrap_coordinate(y - 1, grid.height)
15# nx = wrap_coordinate(x + 1, grid.width)
16# neighbor_value = grid[ny][nx]

Topology (Flat Grid / Torus)

0.00

3 min read•1 page

Our Reaction-Diffusion simulation is just a giant grid of floating-point numbers. To make it useful for Computational Geometry, we must map those numbers onto a physical 3D mesh.

Displacement Mapping:

1. We generate a flat 3D Plane with thousands of vertices, perfectly matching the resolution of our grid. 2. We read the Chemical A and Chemical B values for every vertex. 3. We use the mathematical difference between the two chemicals to drive the Z-height (Vertical Displacement) of the vertex. 4. The invisible chemical reactions are instantly transformed into physical ridges, valleys, and bumps!

python

1# Data Visualization
2def render_grid_to_mesh(grid, mesh):
3    """Converting chemical data into physical shape."""
4    
5    # 1. We have a flat 3D Plane Mesh with lots of vertices
6    for vertex in mesh.vertices:
7        
8        # 2. Get the chemical concentration at this exact coordinate
9        A = grid[vertex.y][vertex.x]["A"]
10        B = grid[vertex.y][vertex.x]["B"]
11        
12        # 3. We calculate a "Difference" score
13        # If A is high and B is low, the number is positive.
14        # If B is high and A is low, the number is negative.
15        score = A - B
16        
17        # 4. We push the vertex UP or DOWN based on the score!
18        vertex.z = score * 5.0
19        
20    return mesh

Vertical Displacement Multiplier

0.00

9 min read•2 pages

By simulating Gray-Scott equations inside a hardware-accelerated WebGL GPGPU shader, we solve chemical diffusion in real-time. Varying the Feed and Kill rates yields self-organizing organic patterns resembling cell division, coral branching, or leopard spots.

Simulation Rules:

1. Mitosis (0.0367, 0.0649): Concentrated spots split outward as they grow, mirroring biological cell division. 2. Coral (0.0545, 0.0620): Generates branching, maze-like lines that fill the space. 3. Spots (0.0300, 0.0620): Produces stable, isolated circular spots of concentration B.

python

1def step_gray_scott(A, B, dA, dB, feed, kill):
2    nextA = copy(A)
3    nextB = copy(B)
4    for y in range(1, height-1):
5        for x in range(1, width-1):
6            # 1. Calculate laplacians
7            lapA = laplacian(A, x, y)
8            lapB = laplacian(B, x, y)
9            
10            # 2. Apply Gray-Scott update rules
11            a, b = A[y][x], B[y][x]
12            nextA[y][x] = a + dA * lapA - a*b*b + feed * (1.0 - a)
13            nextB[y][x] = b + dB * lapB + a*b*b - (kill + feed) * b
14    return nextA, nextB

Gray-Scott Pattern

Run Simulation