Python Cellular Automata: From Wolfram Rules to the Game of Life

Cellular automata are computational models where a grid of cells evolves over time according to a fixed set of rules. Despite their simplicity, they produce extraordinarily complex behavior -- from fractal-like patterns in one dimension to self-replicating structures in two. Python and NumPy make it straightforward to build, visualize, and experiment with these systems from scratch.

This article walks through the core concepts behind cellular automata, then builds two complete implementations in Python: a one-dimensional elementary automaton using Wolfram's rule numbering system, and Conway's Game of Life in two dimensions. All code uses only Python standard libraries and NumPy, so there are no extra dependencies to install.

What Are Cellular Automata

A cellular automaton (CA) is a discrete computational system that operates on a grid of cells. Each cell holds a state -- typically a simple value like 0 or 1. At each time step, every cell updates its state simultaneously based on a rule that examines the cell's current state and the states of its neighbors.

The concept was pioneered by John von Neumann and Stanislaw Ulam in the 1940s. Von Neumann was investigating self-replicating systems and designed a theoretical automaton that could reproduce itself on a 2D grid. Decades later, Stephen Wolfram systematically classified one-dimensional cellular automata, and John Conway created the Game of Life, which demonstrated that a handful of rules on a grid could produce patterns capable of universal computation.

The defining characteristics of a cellular automaton are:

• A grid of cells: This can be one-dimensional (a row), two-dimensional (a flat grid), or higher. Each cell holds a discrete state.
• A neighborhood: A defined set of nearby cells that influence a given cell's next state. In 1D, this is typically the cell and its immediate left and right neighbors. In 2D, common neighborhoods include the Moore neighborhood (8 surrounding cells) and the von Neumann neighborhood (4 cardinal neighbors).
• A transition rule: A function that maps the current configuration of a cell and its neighborhood to the cell's next state.
• Synchronous updates: All cells update at the same time, based on the state from the previous step.

Elementary Cellular Automata and Wolfram Rules

An elementary cellular automaton is the simplest possible CA: a one-dimensional row of binary cells (each cell is either 0 or 1), where the next state of a cell depends on its own current state and the states of its two immediate neighbors. That gives three inputs, each binary, producing 2^3 = 8 possible neighborhood configurations.

For each of those 8 configurations, the rule specifies whether the resulting cell should be 0 or 1. Since there are 8 outputs and each is binary, there are 2^8 = 256 possible rules. Wolfram assigned each rule a number from 0 to 255 based on the binary representation of its output pattern.

Here is how Rule 30 works. The eight possible neighborhood patterns (left, center, right) are listed from 111 down to 000. Rule 30 in binary is 00011110. Reading left to right, the outputs for patterns 111, 110, 101, 100, 011, 010, 001, 000 are 0, 0, 0, 1, 1, 1, 1, 0 respectively.

Rule 30 is particularly interesting because it produces chaotic, aperiodic output from a single active cell. It was famously used as a pseudorandom number generator in Wolfram's Mathematica software. Rule 90 produces a Sierpinski triangle pattern. Rule 110, as mentioned, is Turing complete.

Building an Elementary CA in Python

The implementation below creates a function that accepts any Wolfram rule number (0-255), generates the corresponding rule lookup table, and evolves a 1D grid over a specified number of time steps. The result is a 2D NumPy array where each row represents the state of all cells at one time step.

import numpy as np


def wolfram_rule_table(rule_number):
    """Convert a Wolfram rule number (0-255) to a lookup dictionary.

    Each key is a 3-tuple representing (left, center, right) neighbor states.
    Each value is the resulting cell state (0 or 1).
    """
    rule_binary = format(rule_number, '08b')
    table = {}
    for i in range(8):
        # Convert index to a 3-bit neighborhood pattern
        pattern = tuple(int(b) for b in format(7 - i, '03b'))
        table[pattern] = int(rule_binary[i])
    return table


def evolve_1d(rule_number, width=101, steps=50):
    """Evolve a 1D elementary cellular automaton.

    Args:
        rule_number: Wolfram rule number between 0 and 255.
        width: Number of cells in the row.
        steps: Number of time steps to simulate.

    Returns:
        A 2D NumPy array of shape (steps + 1, width). Each row is one
        generation. Row 0 is the initial state with a single active cell
        in the center.
    """
    table = wolfram_rule_table(rule_number)
    grid = np.zeros((steps + 1, width), dtype=int)
    grid[0, width // 2] = 1  # seed: single cell in the center

    for t in range(steps):
        for i in range(width):
            left = grid[t, (i - 1) % width]
            center = grid[t, i]
            right = grid[t, (i + 1) % width]
            grid[t + 1, i] = table[(left, center, right)]

    return grid

The function wolfram_rule_table converts a rule number into its 8-bit binary representation and maps each 3-cell neighborhood pattern to an output. The evolve_1d function initializes a grid with a single active cell in the center, then iterates through each time step, applying the rule to every cell. Periodic (wrapping) boundary conditions are used, so the leftmost and rightmost cells treat each other as neighbors.

Visualizing the Result

Use Matplotlib to render the evolution as a black-and-white image, where each pixel corresponds to a cell at a given time step:

import matplotlib.pyplot as plt


def plot_1d_ca(grid, rule_number):
    """Display a 1D CA evolution as a binary image."""
    plt.figure(figsize=(10, 5))
    plt.imshow(grid, cmap='binary', interpolation='nearest')
    plt.title(f'Rule {rule_number}')
    plt.xlabel('Cell Index')
    plt.ylabel('Time Step')
    plt.tight_layout()
    plt.show()


# Generate and display several well-known rules
for rule in [30, 90, 110, 184]:
    grid = evolve_1d(rule, width=201, steps=100)
    plot_1d_ca(grid, rule)

Vectorized Version with NumPy

The nested loop in evolve_1d works well for learning, but it becomes slow at larger scales. Here is a vectorized version that processes all cells in a row simultaneously using NumPy's roll function:

def evolve_1d_fast(rule_number, width=201, steps=100):
    """Vectorized 1D elementary CA evolution using NumPy.

    Instead of looping over each cell, this function computes the
    neighborhood index for all cells at once using bit shifting,
    then looks up the result in a precomputed array.
    """
    # Build a lookup array: index is the 3-bit neighborhood value (0-7),
    # value is the output state
    rule_bits = np.array(
        [int(b) for b in format(rule_number, '08b')][::-1],
        dtype=int
    )

    grid = np.zeros((steps + 1, width), dtype=int)
    grid[0, width // 2] = 1

    for t in range(steps):
        left = np.roll(grid[t], 1)
        center = grid[t]
        right = np.roll(grid[t], -1)
        # Combine into a 3-bit index: left*4 + center*2 + right
        neighborhood_index = left * 4 + center * 2 + right
        grid[t + 1] = rule_bits[neighborhood_index]

    return grid

This version encodes each neighborhood as a 3-bit integer (left * 4 + center * 2 + right) and uses it to index directly into the rule_bits array. The np.roll function handles the periodic boundary conditions by shifting the array. On a grid of 1,000 cells over 500 steps, the vectorized version runs roughly 50 to 100 times faster than the cell-by-cell loop.

Conway's Game of Life

Conway's Game of Life is a two-dimensional cellular automaton invented by mathematician John Conway in 1970. It uses a grid of binary cells with the Moore neighborhood (all 8 surrounding cells). The rules are:

1. Underpopulation: A live cell with fewer than 2 live neighbors dies.
2. Survival: A live cell with 2 or 3 live neighbors stays alive.
3. Overpopulation: A live cell with more than 3 live neighbors dies.
4. Reproduction: A dead cell with exactly 3 live neighbors becomes alive.

These four rules produce an astonishing variety of emergent behavior. Stable patterns ("still lifes") like the block and beehive remain unchanged. Oscillators like the blinker and pulsar repeat a cycle of states. Spaceships like the glider move across the grid. And complex structures called "guns" can produce an infinite stream of gliders. The Game of Life has been proven to be Turing complete, meaning it can simulate any computation.

Building the Game of Life in Python

The implementation below uses a 2D NumPy array and scipy.signal.convolve2d to count neighbors efficiently. This avoids the need for nested loops over every cell and every neighbor:

import numpy as np
from scipy.signal import convolve2d


def game_of_life_step(grid):
    """Compute one generation of Conway's Game of Life.

    Uses 2D convolution with a 3x3 kernel of ones (center zeroed)
    to count live neighbors for every cell simultaneously.

    Args:
        grid: 2D NumPy array of 0s and 1s.

    Returns:
        A new 2D NumPy array representing the next generation.
    """
    # Kernel counts all 8 neighbors (Moore neighborhood)
    kernel = np.array([[1, 1, 1],
                       [1, 0, 1],
                       [1, 1, 1]])

    neighbor_count = convolve2d(grid, kernel, mode='same', boundary='wrap')

    # Apply the four rules using boolean logic
    new_grid = np.zeros_like(grid)
    # Survival: live cell with 2 or 3 neighbors stays alive
    new_grid[(grid == 1) & ((neighbor_count == 2) | (neighbor_count == 3))] = 1
    # Reproduction: dead cell with exactly 3 neighbors becomes alive
    new_grid[(grid == 0) & (neighbor_count == 3)] = 1

    return new_grid


def run_game_of_life(height=100, width=100, steps=200, density=0.3):
    """Run the Game of Life simulation.

    Args:
        height: Number of rows in the grid.
        width: Number of columns in the grid.
        steps: Number of generations to simulate.
        density: Probability that each cell starts alive.

    Returns:
        A list of 2D NumPy arrays, one per generation.
    """
    # Random initial state
    grid = (np.random.random((height, width)) < density).astype(int)
    history = [grid.copy()]

    for _ in range(steps):
        grid = game_of_life_step(grid)
        history.append(grid.copy())

    return history

The key insight is the convolution approach. The kernel is a 3x3 matrix of ones with a zero in the center. When convolved with the grid, it produces a new array where each cell contains the count of its live neighbors. This replaces what would otherwise be a double-nested loop checking all 8 neighbors of every cell.

Visualizing with Matplotlib Animation

To watch the Game of Life evolve in real time, use Matplotlib's animation module:

import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation


def animate_life(history, interval=50):
    """Animate the Game of Life history as a Matplotlib animation.

    Args:
        history: List of 2D arrays from run_game_of_life.
        interval: Milliseconds between frames.
    """
    fig, ax = plt.subplots(figsize=(8, 8))
    ax.set_axis_off()
    img = ax.imshow(history[0], cmap='binary', interpolation='nearest')

    def update(frame):
        img.set_data(history[frame])
        ax.set_title(f'Generation {frame}')
        return [img]

    anim = FuncAnimation(
        fig,
        update,
        frames=len(history),
        interval=interval,
        blit=True
    )
    plt.tight_layout()
    plt.show()
    return anim


# Run and animate
history = run_game_of_life(height=80, width=80, steps=300, density=0.3)
anim = animate_life(history)

Placing Known Patterns

Instead of starting from a random grid, you can place well-known patterns at specific positions. Here is a helper function and a few classic patterns:

def place_pattern(grid, pattern, row, col):
    """Place a 2D pattern onto the grid at the given position.

    Args:
        grid: The target 2D NumPy array.
        pattern: A 2D list or array of 0s and 1s.
        row: Top-left row position for the pattern.
        col: Top-left column position for the pattern.
    """
    pattern = np.array(pattern)
    h, w = pattern.shape
    grid[row:row + h, col:col + w] = pattern


# Classic patterns
GLIDER = [
    [0, 1, 0],
    [0, 0, 1],
    [1, 1, 1]
]

BLINKER = [
    [1, 1, 1]
]

LIGHTWEIGHT_SPACESHIP = [
    [0, 1, 0, 0, 1],
    [1, 0, 0, 0, 0],
    [1, 0, 0, 0, 1],
    [1, 1, 1, 1, 0]
]

GLIDER_GUN = [
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1],
    [0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1],
    [1,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [1,1,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,1,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
]


# Example: place a Gosper glider gun on a grid
grid = np.zeros((60, 60), dtype=int)
place_pattern(grid, GLIDER_GUN, 10, 5)
history = [grid.copy()]
for _ in range(300):
    grid = game_of_life_step(grid)
    history.append(grid.copy())
anim = animate_life(history)

Extending Your Automata

Once you have the basic framework running, there are several directions to explore:

Multiple states. Instead of binary cells, use k possible states with totalistic rules, where the next state depends on the sum of neighbor values rather than their exact arrangement. Wolfram's totalistic rules generalize elementary CA to any number of colors.

Larger neighborhoods. Increase the radius r so each cell considers neighbors further away. A radius-2 elementary CA examines 5 cells (2 on each side plus center), producing 2^(2^5) = over 4 billion possible rules.

Continuous values. Replace discrete states with floating-point values between 0.0 and 1.0, and use averaging or threshold functions as transition rules. This leads to reaction-diffusion systems and Lenia, a continuous generalization of the Game of Life that produces lifelike, organic-looking structures.

Alternative rule sets. Try HighLife (adds a reproduction rule for 6 neighbors), Brian's Brain (three states with a "dying" phase), or WireWorld (used to simulate electronic circuits). Each uses the same grid-and-neighbor framework with different transition logic.

Performance at scale. For very large grids, consider HashLife, an algorithm that uses memoization to skip ahead exponentially many generations. The golly application implements HashLife and can simulate trillions of generations in seconds for structured patterns.

Key Takeaways

1. Simple rules, complex behavior: Cellular automata demonstrate that a small set of deterministic rules applied to a grid can produce unpredictable, chaotic, and even computationally universal behavior.
2. Wolfram's numbering system: All 256 elementary 1D cellular automata can be encoded as a single integer from 0 to 255, making it easy to enumerate, compare, and experiment with every possible rule.
3. NumPy vectorization matters: Replacing cell-by-cell loops with array operations like np.roll and convolve2d can yield speedups of 50x or more, making it practical to simulate large grids over hundreds of generations.
4. The Game of Life is more than a toy: Conway's creation is Turing complete and has an active community that continues to discover new patterns, including self-replicating structures and programmable logic gates built entirely from gliders.
5. Endless room to explore: Multi-state rules, continuous-valued Lenia systems, genetic algorithm-driven rule evolution, and alternative neighborhoods all build on the same core framework covered in this article.

Cellular automata sit at the intersection of computer science, mathematics, and art. The code in this article gives you a working foundation -- a set of functions you can import into any project, extend with new rules, and use to explore one of the richest areas of computational modeling.