Evolving Digital Universes: When Genetic Algorithms Breed the Laws of Physics

What if you could breed the laws of physics? Not tweak a gravitational constant or adjust the speed of light—but start with thousands of random, nonsensical rule sets governing a tiny universe, let the ones that produce interesting behavior survive, and iterate for hundreds of generations until something extraordinary emerges from nothing. That's what happens when you combine two of computer science's deepest ideas: cellular automata and genetic algorithms. One builds universes. The other evolves them. Together, they generate emergent phenomena rarely before encountered.

This article assumes comfort with Python, familiarity with both cellular automata and evolutionary computation, and a willingness to think about computation as something that can be grown rather than engineered. Every concept discussed here is grounded in verifiable research—from Christopher Langton's lambda parameter to the density classification problem that stumped human designers for decades until genetic algorithms solved it.

We're not building a toy. We're building an evolutionary engine that discovers rules capable of producing computation, self-organization, and emergent behavior—all from a randomly initialized search.

The Architecture of Tiny Universes

A cellular automaton is deceptively simple: a grid of cells, each in one of a finite number of states, updated simultaneously according to a local rule that considers only a cell's immediate neighborhood. No global coordination. No central controller. Just cells blindly following a lookup table. And yet, from this simplicity, structures emerge that move, replicate, compute, and self-organize.

The canonical starting point is the one-dimensional elementary cellular automaton, studied systematically by Stephen Wolfram beginning in 1983–1984. Each cell has two possible states (0 or 1), and its next state is determined by its own current state plus the states of its two immediate neighbors. Three binary inputs produce eight possible input configurations, and a binary output for each means the entire rule is a single 8-bit number. Exactly 256 possible universes.

Wolfram cataloged all 256 and found that they fall into four behavioral classes: Class I (static, dead), Class II (periodic, repetitive), Class III (chaotic, pseudorandom), and—the rare and fascinating Class IV—complex. Class IV rules produce localized structures that interact, persist, and compute. Rule 110, a Class IV automaton, is the only elementary CA for which Turing completeness has been directly proven, meaning any algorithm that can be computed at all can, in principle, be computed within Rule 110. Matthew Cook developed the proof while working as Wolfram's research assistant, identifying the core construction in the mid-1990s by showing Rule 110 could emulate a cyclic tag system—a known universal computational model. Details were refined and the proof completed by 1998. Cook presented the completed proof publicly at the Santa Fe Institute conference CA98 in 1998—an act that led Wolfram Research to obtain a temporary restraining order blocking the conference proceedings from publication, on the grounds that Cook had violated his NDA; the publisher ultimately complied with Wolfram Research's demands even after the order expired. The proof's existence was now known, but formally unpublishable. Cook's paper finally appeared in Complex Systems in 2004 after years of legal limbo. A universe in 8 bits, its universality suppressed for six years.

Here's the engine that runs any elementary CA:

import numpy as np

def elementary_ca(
    rule_number: int,
    width: int = 201,
    steps: int = 100,
    initial_state: np.ndarray | None = None,
) -> np.ndarray:
    """Simulate a 1D elementary cellular automaton.

    The rule_number encodes all 8 output bits of the transition function.
    Each bit position corresponds to a specific 3-cell neighborhood pattern.
    """
    # Decode the rule number into its 8-bit lookup table
    rule_binary = np.array(
        [int(b) for b in format(rule_number, '08b')[::-1]],
        dtype=np.uint8
    )

    grid = np.zeros((steps, width), dtype=np.uint8)

    if initial_state is not None:
        grid[0] = initial_state[:width]
    else:
        grid[0, width // 2] = 1  # single seed cell

    for t in range(1, steps):
        for i in range(1, width - 1):
            # Convert the 3-cell neighborhood into an index (0-7)
            neighborhood = (
                (grid[t-1, i-1] << 2)
                | (grid[t-1, i] << 1)
                | grid[t-1, i+1]
            )
            grid[t, i] = rule_binary[neighborhood]

    return grid

The key insight is that rule_binary is a genome. It's a compact encoding of every decision the universe will ever make. Change one bit and you get an entirely different reality. But 256 rules can be explored exhaustively. The real challenge—and the reason genetic algorithms earn their keep—arrives when we scale up.

A two-state, three-neighbor rule needs 8 bits. A three-state, five-neighbor rule needs 3⁵ = 243 entries, each with three possible outputs—that's 3²⁴³ possible rules, a number so astronomically large it exceeds the number of particles in the observable universe by hundreds of orders of magnitude. You can't search this space by hand. You need evolution.

Langton's Lambda: Mapping the Edge of Chaos

Before we evolve rules blindly, it helps to understand the landscape we're searching. In 1990, Christopher Langton introduced a single parameter—lambda (λ)—that characterizes a cellular automaton rule by measuring what fraction of its lookup table entries map to a non-quiescent (typically non-zero) output state. A rule where every entry maps to the quiescent state has λ = 0. A rule where no entries do has λ = 1. Langton's central empirical finding was that complex, Class IV behavior tends to cluster near a phase transition between the ordered (λ near 0) and chaotic (λ near 1) regimes—what he called the "edge of chaos."

This is a profound but contested idea. Langton demonstrated it empirically using K=4, N=5 CAs (four states, neighborhood of five). The critical λ value at which this transition occurs is not a universal constant: it shifts depending on the number of states, neighborhood structure, and the path taken through rule space. Mitchell, Hraber, and Crutchfield later showed in 1993 that high-performance rules for the density classification task are not reliably concentrated near critical λ values—evolution finds them, but lambda alone does not predict them. Lambda is a useful navigational prior, not a guarantee of complexity.

def langton_lambda(genome: list[int], num_states: int) -> float:
    """Calculate Langton's lambda parameter for a rule genome.

    Lambda measures the fraction of non-quiescent transitions.
    The quiescent state is conventionally state 0.

    Returns a value between 0.0 (all-dead) and 1.0 (no dead transitions).
    Langton's 1990 paper showed empirically that complex behavior tends to
    cluster near a phase-transition lambda value that is context-dependent
    (varies with state count, neighborhood size, and the path through rule
    space). There is no single universal critical lambda.
    """
    quiescent_count = sum(1 for g in genome if g == 0)
    return 1.0 - (quiescent_count / len(genome))


def lambda_critical(num_states: int) -> float:
    """Heuristic approximation of a useful lambda target for initialization.

    Langton's empirical work suggested complex behavior concentrates near
    a phase transition. The value 1 - 1/num_states is a common heuristic
    approximation for this region (e.g., ~0.5 for binary, ~0.667 for ternary).
    It is NOT a derived theorem: Mitchell et al. (1993) showed that
    high-performing evolved rules need not cluster at this lambda value.
    Use it as a biased starting point, not a guarantee.
    """
    return 1.0 - (1.0 / num_states)

Here's how to use lambda-biased initialization in practice:

def create_lambda_biased_genome(
    num_states: int,
    neighborhood_size: int,
    target_lambda: float | None = None,
    spread: float = 0.1,
) -> list[int]:
    """Generate a genome biased toward a target lambda value.

    If target_lambda is None, defaults to the critical value.
    The 'spread' parameter controls how tightly genomes cluster
    around the target (smaller = tighter).
    """
    if target_lambda is None:
        target_lambda = lambda_critical(num_states)

    num_entries = num_states ** neighborhood_size

    # Sample an actual lambda from a normal distribution around the target
    actual_lambda = np.clip(
        np.random.normal(target_lambda, spread), 0.05, 0.95
    )

    # Determine how many entries should be non-quiescent
    num_nonquiescent = int(round(actual_lambda * num_entries))

    genome = [0] * num_entries

    # Randomly select which positions get non-quiescent values
    nonquiescent_positions = np.random.choice(
        num_entries, size=num_nonquiescent, replace=False
    )

    for pos in nonquiescent_positions:
        # Assign a random non-zero state
        genome[pos] = np.random.randint(1, num_states)

    return genome

Encoding Laws of Physics as DNA

A genetic algorithm treats candidate solutions as organisms. Each organism carries a genome—a data structure encoding a complete solution—and is subjected to selection pressure based on how well it performs against a fitness criterion. The strong survive. The weak are culled. Survivors recombine and mutate to produce the next generation.

For cellular automata, the genome is the rule table itself: a flat array where each index represents a neighborhood configuration and each value represents the output state. Each RuleOrganism isn't a solution to a math problem or an optimized parameter set. It's a complete set of physical laws for a computational universe. Mutation doesn't tweak a variable—it rewrites a law. Crossover doesn't blend two answers—it splices two realities together and asks whether the chimera can survive.

from dataclasses import dataclass, field
import random

@dataclass
class RuleOrganism:
    """A candidate cellular automaton rule, treated as a living organism.

    The genome is a lookup table: genome[neighborhood_index] = output_state.
    Its length is num_states ** neighborhood_size.
    """

    num_states: int
    neighborhood_size: int
    genome: list[int] = field(default_factory=list)
    fitness: float = 0.0

    def __post_init__(self):
        if not self.genome:
            # Lambda-biased initialization by default
            self.genome = create_lambda_biased_genome(
                self.num_states, self.neighborhood_size
            )

    @property
    def lam(self) -> float:
        """This organism's Langton lambda value."""
        return langton_lambda(self.genome, self.num_states)

    def mutate(self, mutation_rate: float = 0.02):
        """Point mutations: randomly alter individual rule entries.

        Mutations are lambda-aware: they swap between states rather than
        always randomizing, which prevents large lambda shifts from
        single mutation events.
        """
        for i in range(len(self.genome)):
            if random.random() < mutation_rate:
                # Swap to a different state (never mutate to current value)
                current = self.genome[i]
                options = [s for s in range(self.num_states) if s != current]
                self.genome[i] = random.choice(options)

    def crossover(self, other: 'RuleOrganism') -> 'RuleOrganism':
        """Uniform crossover: each gene independently chosen from a parent.

        Uniform crossover preserves lambda better than single-point or
        two-point crossover because it doesn't create large contiguous
        blocks from one parent.
        """
        child_genome = [
            self.genome[i] if random.random() < 0.5 else other.genome[i]
            for i in range(len(self.genome))
        ]
        return RuleOrganism(
            num_states=self.num_states,
            neighborhood_size=self.neighborhood_size,
            genome=child_genome,
        )

The Hardest Question: What Does "Fit" Mean?

Here's where many implementations collapse—and where the real intellectual work lives. Fitness is everything in a genetic algorithm. It is the selection pressure. It is the environment. It determines what "good" means for a universe.

If you want static rules, measure fitness by how quickly the grid reaches a fixed point. If you want chaos, reward high entropy. But if you want complexity—the deep, structured, Class IV behavior that produces computation and emergent phenomena—you need a fitness function that can distinguish interesting from merely noisy.

Several approaches have proven effective in published research, and the strongest implementations combine multiple signals:

import zlib
from scipy.stats import entropy as scipy_entropy

def compression_ratio(grid: np.ndarray) -> float:
    """Kolmogorov complexity proxy via compression ratio.

    This exploits a deep connection: the compressibility of a string
    is related to its algorithmic information content. Random data
    is incompressible (ratio near 1.0). Repetitive data compresses
    to almost nothing (ratio near 0.0). Complex, structured data
    with internal regularities sits in a sweet spot around 0.3-0.5.

    Empirically consistent with known Class IV rules: Rule 110 and
    similar complex CAs produce grids that typically compress to
    roughly 40-60% of raw size, varying with grid dimensions and
    initial conditions. This range is a useful target, not a
    hard constant.
    """
    raw_bytes = grid.tobytes()
    compressed = zlib.compress(raw_bytes, level=9)
    return len(compressed) / len(raw_bytes)


def spatial_entropy(grid: np.ndarray) -> float:
    """Shannon entropy of the final row's state distribution."""
    final_row = grid[-1]
    _, counts = np.unique(final_row, return_counts=True)
    return float(scipy_entropy(counts, base=2))


def temporal_diversity(grid: np.ndarray) -> float:
    """Fraction of time steps that produce unique row configurations.

    Static universes: ~1/steps (near 0). Chaotic universes: ~1.0.
    Complex universes with gliders and interactions: 0.5-0.8.
    """
    unique_rows = len(set(map(tuple, grid)))
    return unique_rows / len(grid)


def input_entropy(grid: np.ndarray) -> float:
    """Measure sensitivity to initial conditions.

    Run the same rule from two slightly different initial conditions
    and measure how rapidly the outputs diverge. This is effectively
    a discrete Lyapunov exponent.

    High sensitivity = chaotic. Zero sensitivity = ordered.
    Moderate, distance-dependent sensitivity = complex.
    """
    rows, cols = grid.shape
    row_entropies = []
    for t in range(rows):
        _, counts = np.unique(grid[t], return_counts=True)
        row_entropies.append(float(scipy_entropy(counts, base=2)))

    if len(row_entropies) < 2:
        return 0.0

    # Measure the variance of entropy over time
    # Constant entropy = boring. Wildly varying = complex dynamics.
    return float(np.std(row_entropies))

The compression ratio deserves particular attention. It exploits the theoretical connection between Kolmogorov complexity—the length of the shortest program that produces a given string—and compressibility. While Kolmogorov complexity is uncomputable in general (proven by Gregory Chaitin), compression algorithms like zlib provide a practical upper bound. The key insight is that complex cellular automata produce output that contains real structure (making it more compressible than random noise) but not trivially repetitive structure (making it less compressible than periodic patterns). This sweet spot is detectable and measurable.

A composite fitness function targeting Class IV behavior combines these signals with empirically calibrated weights:

def complexity_fitness(grid: np.ndarray) -> float:
    """Composite fitness targeting structured complexity.

    The ideal organism produces behavior that is:
    - Neither fully compressible nor incompressible (structured)
    - Spatially varied but not uniformly random (patterned)
    - Temporally diverse but not maximally chaotic (dynamic)
    - Entropically variable over time (alive)
    """
    cr = compression_ratio(grid)
    se = spatial_entropy(grid)
    td = temporal_diversity(grid)
    ie = input_entropy(grid)

    # Target the complexity sweet spot for each metric
    cr_score = 1.0 - abs(cr - 0.45) * 3.0
    se_score = min(se / 1.0, 1.0)          # normalize for binary states
    td_score = 1.0 - abs(td - 0.65) * 2.5
    ie_score = min(ie * 8.0, 1.0)          # reward entropic variability

    return max(0.0, (
        cr_score * 0.30
        + se_score * 0.20
        + td_score * 0.25
        + ie_score * 0.25
    ))

The Density Classification Problem: A Concrete Benchmark

Abstract fitness functions are powerful, but there's a way to make the evolutionary approach concrete and verifiable: solve the density classification problem (DCP). This is one of the most studied benchmarks in cellular automata research. The GKL rule itself—named after Gacs, Kurdyumov, and Levin, who published it in a 1978 Soviet information-theory journal as part of research on fault-tolerant computation—was the reigning champion hand-designed solution. Melanie Mitchell, James Crutchfield, and Patrik Hraber then reframed the GKL rule's task as a formal benchmark in the early 1990s, launching the research program that genetic algorithms would eventually crack.

The problem is simple to state: given a one-dimensional binary CA initialized with a random configuration of 0s and 1s, the rule must eventually converge to all-1s if the initial density of 1s exceeds 50%, and to all-0s otherwise. The CA must perform a global computation—determining majority—using only local rules. Each cell can only see its immediate neighbors, yet the collective must reach a consensus about the entire initial state.

No finite-radius two-state CA can solve this perfectly. Land and Belew proved this in a 1995 paper in Physical Review Letters: for any binary CA with a fixed neighborhood radius, there exist adversarially chosen initial configurations near the 50% boundary that it will misclassify. This is a worst-case result, not an upper bound on performance over uniformly random inputs—a CA can achieve higher accuracy on random initial conditions while still failing these adversarial edge cases. But genetic algorithms have discovered rules with remarkable performance on random initial configurations. The GKL rule itself achieves approximately 76–78% accuracy on uniformly random initial conditions over grids of width 149. Mitchell, Hraber, and Crutchfield evolved radius-3 rules that reached comparable accuracy while using qualitatively different computational strategies—strategies that, as Mitchell and Crutchfield's computational mechanics framework later revealed, had never been conceived by human designers.

def density_classification_fitness(
    organism: 'RuleOrganism',
    width: int = 149,
    steps: int = 320,
    num_trials: int = 50,
) -> float:
    """Fitness for the density classification problem.

    Tests the organism's rule against multiple random initial
    configurations and measures classification accuracy.

    Uses odd width (149) by convention so that exact 50/50 splits
    are impossible, guaranteeing a definite majority.

    The step count (320, slightly above the M=2N=298 convention used by
    Mitchell et al.) gives the CA enough time to converge while providing
    a small buffer for near-threshold configurations. To replicate their
    exact baseline, use steps=298.
    """
    correct = 0

    for _ in range(num_trials):
        # Generate random initial configuration
        initial = np.random.randint(0, 2, size=width).astype(np.uint8)
        density = initial.mean()

        # Run the CA
        grid = simulate_ca(organism, width, steps, initial_state=initial)
        final_row = grid[-1]

        # Check: did it converge to the correct uniform state?
        if density > 0.5 and final_row.mean() > 0.95:
            correct += 1
        elif density < 0.5 and final_row.mean() < 0.05:
            correct += 1

    return correct / num_trials

What makes the DCP particularly revealing is what the evolved rules look like. In their 1994 Physica D paper "Evolving Cellular Automata to Perform Computations," Mitchell, Crutchfield, and Hraber analyzed the best evolved solutions and found that they work by generating particles—boundaries between regions of uniform state that propagate at different velocities. When particles from low-density and high-density regions collide, the collision dynamics determine which state wins. The CA performs a distributed vote through wave-like information transfer. Mitchell and Crutchfield developed a formal framework they called computational mechanics to analyze this, identifying the specific domain boundaries and particle interactions that perform the computation. It was a strategy that no human had conceived before evolution discovered it—and it remains one of the most compelling documented examples of evolutionary computation genuinely surprising its creators.

The Evolutionary Loop

With genomes defined and fitness measurable, the evolutionary engine itself is architecturally straightforward—but the details of selection strategy, diversity maintenance, and adaptive control determine whether evolution discovers anything interesting or collapses into mediocrity:

from typing import Callable

def evolve_universes(
    population_size: int = 200,
    generations: int = 500,
    num_states: int = 2,
    neighborhood_size: int = 7,     # radius-3 for DCP
    grid_width: int = 149,
    simulation_steps: int = 320,
    elite_fraction: float = 0.05,
    tournament_size: int = 5,
    mutation_rate: float = 0.02,
    fitness_fn: Callable = complexity_fitness,
) -> list[RuleOrganism]:
    """Evolve a population of cellular automaton rules.

    Uses tournament selection with elitism and adaptive mutation
    to balance exploration and exploitation.
    """
    population = [
        RuleOrganism(num_states, neighborhood_size)
        for _ in range(population_size)
    ]

    best_ever = None
    stagnation_counter = 0
    history = []  # track fitness over generations

    for gen in range(generations):

        # === EVALUATE ===
        for organism in population:
            grid = simulate_ca(organism, grid_width, simulation_steps)
            organism.fitness = fitness_fn(grid)

        population.sort(key=lambda o: o.fitness, reverse=True)

        # Track progress
        gen_best = population[0].fitness
        gen_mean = np.mean([o.fitness for o in population])
        gen_diversity = np.std([o.lam for o in population])
        history.append((gen, gen_best, gen_mean, gen_diversity))

        # Stagnation detection
        if best_ever is None or gen_best > best_ever.fitness:
            best_ever = RuleOrganism(
                num_states=num_states,
                neighborhood_size=neighborhood_size,
                genome=list(population[0].genome),
            )
            best_ever.fitness = gen_best
            stagnation_counter = 0
        else:
            stagnation_counter += 1

        # Adaptive mutation: escalate when evolution stalls
        current_mutation = min(
            mutation_rate * (1 + stagnation_counter * 0.15),
            0.25  # cap to prevent total randomization
        )

        # === SELECTION & REPRODUCTION ===
        num_elite = max(2, int(population_size * elite_fraction))
        next_generation = [
            RuleOrganism(
                num_states=o.num_states,
                neighborhood_size=o.neighborhood_size,
                genome=list(o.genome),
            )
            for o in population[:num_elite]
        ]

        while len(next_generation) < population_size:
            parent_a = tournament_select(population, tournament_size)
            parent_b = tournament_select(population, tournament_size)
            child = parent_a.crossover(parent_b)
            child.mutate(current_mutation)
            next_generation.append(child)

        population = next_generation

        # Progress reporting
        if gen % 50 == 0:
            print(
                f"Gen {gen:4d} | Best: {gen_best:.4f} | "
                f"Mean: {gen_mean:.4f} | "
                f"Lambda diversity: {gen_diversity:.4f} | "
                f"Mutation rate: {current_mutation:.4f}"
            )

    return sorted(population, key=lambda o: o.fitness, reverse=True)


def tournament_select(
    population: list[RuleOrganism],
    k: int,
) -> RuleOrganism:
    """Select a parent by tournament among k random individuals.

    Tournament selection provides controllable selection pressure:
    larger k = stronger pressure (favors the best more aggressively).
    k=2 is mild. k=7+ is aggressive. k=5 is a common middle ground.
    """
    competitors = random.sample(population, k)
    return max(competitors, key=lambda o: o.fitness)

The adaptive mutation rate is a critical survival mechanism. Evolutionary search suffers from premature convergence: the population collapses around a local optimum, and without sufficient variation, it can't escape. By tracking stagnation and escalating mutation pressure when progress stalls, the algorithm maintains its ability to explore. This approach has an analog in biology. Bacteria under environmental stress upregulate error-prone DNA polymerases via the SOS response, effectively increasing their mutation rate when their current strategy isn't working—a phenomenon called stress-induced mutagenesis, well documented in E. coli and other species. Whether this is a selected adaptation or a side effect of DNA-repair mechanisms remains actively debated among evolutionary biologists, which itself says something interesting: evolution may not "intend" to escalate variation under pressure, but does so anyway as a consequence of chemistry. Our algorithm does it on purpose.

The lambda diversity metric in the progress output is equally important. If lambda diversity drops below 0.05, the population is converging prematurely and you should consider injecting random immigrants (fresh random organisms) into the population to restore diversity.

The Simulation Engine

The simulate_ca function bridges the genome to the grid, translating an organism's rule table into actual spatiotemporal dynamics. The critical design decision here is the neighborhood-to-index encoding, which generalizes cleanly from binary to arbitrary state counts:

def simulate_ca(
    organism: RuleOrganism,
    width: int,
    steps: int,
    initial_state: np.ndarray | None = None,
) -> np.ndarray:
    """Run a cellular automaton using the organism's genome as its rule.

    Neighborhoods are converted to genome indices by treating each cell
    in the neighborhood as a digit in a base-N number, where N is the
    number of states. This generalizes beyond binary CAs to arbitrary
    state counts and neighborhood sizes.

    Examples of neighborhood-to-index conversion:
      Binary, 3-cell:   [1, 0, 1] -> 1*4 + 0*2 + 1*1 = 5
      Ternary, 5-cell:  [2, 0, 1, 1, 0] -> 2*81 + 0*27 + 1*9 + 1*3 + 0 = 174
    """
    ns = organism.num_states
    radius = organism.neighborhood_size // 2

    grid = np.zeros((steps, width), dtype=np.uint8)

    if initial_state is not None:
        grid[0, :len(initial_state)] = initial_state[:width]
    else:
        grid[0, width // 2] = 1

    for t in range(1, steps):
        for i in range(radius, width - radius):
            neighborhood_index = 0
            for offset in range(-radius, radius + 1):
                neighborhood_index = (
                    neighborhood_index * ns + grid[t-1, i + offset]
                )
            grid[t, i] = organism.genome[neighborhood_index]

    return grid

This encoding means the genome grows exponentially with the neighborhood size. A radius-1 binary CA has 2³ = 8 genome entries. A radius-3 binary CA has 2⁷ = 128. A radius-3 ternary CA has 3⁷ = 2,187. The search space expands astronomically at each step, which is precisely why exhaustive search fails and evolutionary methods become essential.

Beyond One Dimension

The real power of evolutionary rule discovery becomes apparent in two dimensions, where intuition about rules breaks down entirely. Conway's Game of Life, the famous 2D CA, was hand-tuned over weeks of experimentation. But the Game of Life uses an outer-totalistic rule—one that depends only on the center cell's current state and the sum of neighboring states, not their individual arrangement. This dramatically reduces the search space. A non-totalistic 2D rule with a Moore neighborhood (8 neighbors + center cell) has 2⁹ = 512 input configurations and therefore 2⁵¹² possible rules. To put that in perspective: 2⁵¹² is approximately 10¹⁵⁴, which exceeds the estimated number of atoms in the observable universe (roughly 10⁸⁰) by a factor of about 10⁷⁴. Random exploration is computationally meaningless. This is the regime where evolution stops being a convenience and becomes a necessity.

def make_2d_organism(
    num_states: int = 2,
    neighborhood: str = "moore",
) -> RuleOrganism:
    """Create a rule organism for a 2D cellular automaton.

    For a full (non-totalistic) Moore neighborhood with 2 states,
    the genome has 2^9 = 512 entries. Even this modest configuration
    has 2^512 possible rules, a number (~10^154) that exceeds the
    estimated atoms in the observable universe (~10^80) by ~10^74.

    For tractability, consider using outer-totalistic rules instead,
    where the output depends on the center cell's state and the SUM
    of its neighbors' states. This reduces the genome to
    num_states * (num_states * neighbor_count + 1) entries.
    """
    if neighborhood == "moore":
        neighbor_count = 9  # 8 neighbors + center
    else:  # von_neumann
        neighbor_count = 5  # 4 neighbors + center

    return RuleOrganism(
        num_states=num_states,
        neighborhood_size=neighbor_count,
    )


def make_2d_totalistic_organism(
    num_states: int = 2,
    num_neighbors: int = 8,
) -> dict:
    """Create a totalistic 2D CA organism (more tractable).

    In a totalistic rule, the output depends on:
    - The center cell's current state
    - The sum of all neighbor states

    For binary with Moore neighborhood:
    center has 2 states, sum ranges 0-8, giving 2 * 9 = 18 entries.
    This is a vastly more searchable space than the full 512-entry table.
    """
    max_sum = num_neighbors * (num_states - 1)
    genome_size = num_states * (max_sum + 1)

    genome = [
        random.randint(0, num_states - 1)
        for _ in range(genome_size)
    ]

    return {
        "num_states": num_states,
        "num_neighbors": num_neighbors,
        "genome": genome,
        "genome_size": genome_size,
    }

The totalistic reduction is not just a convenience—it's also the space where Conway found the Game of Life. A binary totalistic rule with a Moore neighborhood has only 18 genome entries (2 center states times 9 possible neighbor sums). That's 2¹⁸ = 262,144 possible rules—large enough to contain surprises, small enough for a genetic algorithm to explore thoroughly. The non-totalistic space is where the truly alien, structurally rich rules hide, but the totalistic space is where you can validate your evolutionary framework against known results.

Novelty Search: Evolution Without a Goal

Everything discussed so far uses objective-based fitness: we define what "good" means and evolve toward it. But there's a radical alternative that often produces more interesting results: abandon objectives entirely.

Novelty search, first introduced by Joel Lehman and Kenneth Stanley at the Artificial Life XI conference in 2008 and formally published in the journal Evolutionary Computation in 2011, replaces fitness with novelty. An organism is rewarded not for being good at anything, but for being different from everything the evolutionary process has seen before. The algorithm maintains an archive of previously encountered behaviors, and each new organism is scored by how far its behavior is from its nearest neighbors in that archive.

The core insight is counterintuitive: explicitly searching for a goal can prevent you from finding it. Objective-based evolution gets trapped by deceptive fitness landscapes where the path to the global optimum requires temporarily decreasing fitness. Novelty search sidesteps this trap entirely by not caring about fitness at all—it just rewards exploration.

def behavioral_fingerprint(grid: np.ndarray) -> np.ndarray:
    """Extract a compact behavioral descriptor from a CA run.

    This fingerprint captures the essential dynamics of the CA
    without being so detailed that every run is trivially unique.

    The descriptor includes:
    - Density over time (coarse-grained to 20 bins)
    - Final-row entropy
    - Compression ratio
    - Temporal autocorrelation at multiple lags
    """
    steps = grid.shape[0]
    bin_size = max(1, steps // 20)

    # Coarse-grained density trajectory
    density_trajectory = [
        grid[i:i + bin_size].mean()
        for i in range(0, steps - bin_size + 1, bin_size)
    ]

    # Temporal autocorrelation (how self-similar is the time series?)
    densities = grid.mean(axis=1)
    autocorr = []
    for lag in [1, 5, 10, 20]:
        if lag < len(densities):
            c = np.corrcoef(densities[:-lag], densities[lag:])[0, 1]
            autocorr.append(c if not np.isnan(c) else 0.0)
        else:
            autocorr.append(0.0)

    features = (
        density_trajectory
        + [spatial_entropy(grid)]
        + [compression_ratio(grid)]
        + autocorr
    )

    return np.array(features, dtype=np.float64)


class NoveltyArchive:
    """Maintains a growing archive of behavioral fingerprints.

    Novelty is measured as the mean distance to the k nearest
    neighbors in the archive. This rewards organisms that behave
    differently from anything previously observed.
    """

    def __init__(self, k_nearest: int = 15, archive_threshold: float = 0.1):
        self.archive: list[np.ndarray] = []
        self.k_nearest = k_nearest
        self.archive_threshold = archive_threshold

    def novelty_score(self, fingerprint: np.ndarray) -> float:
        """Score how novel this behavior is relative to the archive."""
        if len(self.archive) < self.k_nearest:
            return float('inf')  # everything is novel at first

        distances = [
            np.linalg.norm(fingerprint - archived)
            for archived in self.archive
        ]
        distances.sort()
        return float(np.mean(distances[:self.k_nearest]))

    def maybe_add(self, fingerprint: np.ndarray, novelty: float):
        """Add to archive if sufficiently novel."""
        if novelty > self.archive_threshold:
            self.archive.append(fingerprint.copy())

Novelty search applied to cellular automata evolution is particularly potent because the behavioral space of CAs is extraordinarily rich. A novelty-driven search doesn't ask "find me a complex rule"—it asks "find me rules I haven't seen before." In practice, this naturally discovers rules across all four of Wolfram's behavioral classes, including Class IV rules that an objective-based search might miss because they sit in narrow, hard-to-reach regions of the fitness landscape.

The practical approach is often to blend novelty and fitness in a weighted combination, getting the goal-directedness of fitness-based search with the exploration benefits of novelty:

def blended_score(
    fitness: float,
    novelty: float,
    novelty_weight: float = 0.4,
) -> float:
    """Combine fitness and novelty for selection.

    Higher novelty_weight favors exploration.
    Lower novelty_weight favors exploitation.

    A value of 0.3-0.5 typically works well for CA evolution.
    """
    # Normalize novelty to a comparable scale
    normalized_novelty = min(novelty / 2.0, 1.0)
    return (
        (1.0 - novelty_weight) * fitness
        + novelty_weight * normalized_novelty
    )

What Evolution Discovers—And What It Means

When you run this system, a characteristic pattern unfolds. The first few generations produce rules that are almost universally boring: all-dead, all-alive, simple period-2 blinkers. Fitness scores cluster near zero and lambda values are scattered randomly.

Around generation 30-50, the population's lambda distribution tightens around the critical value. Pockets of moderate fitness appear. Rules emerge that sustain activity without dying out—but the activity is typically chaotic, a seething noise with no discernible structure.

By generation 150-250, if your fitness function is well-calibrated, the population begins converging on rules that produce genuine structure: propagating signals, oscillating regions, and interaction dynamics. The compression ratio stabilizes in the target band. Temporal diversity levels off. You're watching evolution discover complexity.

The rules evolution finds are often ones that human designers would be unlikely to generate. They lack the elegant symmetry of Rule 110 or Conway's Game of Life. They're messy, asymmetric, and counterintuitive—but they work. In the density classification problem, the evolved solutions use strategies qualitatively unlike anything previously devised: they generate particle-like boundaries between uniform regions that propagate at different speeds, and the collision dynamics of these particles perform the computation. Mitchell, Crutchfield, and Hraber documented this in detail in their 1994 Physica D paper and developed computational mechanics as a formal framework to analyze it—a framework that remains one of the most rigorous available for understanding emergent computation in dynamical systems.

This mirrors a deeper truth about biological evolution. The solutions it discovers are often inelegant by engineering standards, full of redundancy and seemingly arbitrary decisions. But they function. They're robust. And they reach regions of the solution space that top-down design rarely visits—because the path there requires temporarily sacrificing performance, and no fitness gradient points that way.

The question is not whether cellular automata can compute. The question is whether evolution can discover how they compute—and the answer, decisively demonstrated by decades of research from Langton to Mitchell, Crutchfield, and Hraber to Lehman and Stanley, is yes.

Where to Push Next

If this framework has captured your attention, several directions extend it into genuinely open research territory.

Coevolution and adversarial dynamics. Instead of a fixed fitness function, evolve two populations of rules simultaneously, where each population's fitness depends on its interaction with the other. One population evolves "pattern generators" while the other evolves "pattern classifiers." This produces arms-race dynamics and eliminates the need to design fitness functions by hand—the populations define each other's selection pressure.

Open-ended evolution. Combine novelty search with a quality-diversity framework like MAP-Elites (Mouret and Clune, 2015, arXiv:1504.04909). MAP-Elites maintains an archive that covers the entire behavioral feature space, keeping the highest-performing solution in each region. Unlike a standard GA that converges to one good answer, MAP-Elites produces an illumination of the search space: a diverse library of high-quality solutions spanning every combination of behavioral features you choose to track. Applied to CA evolution, this produces a library of high-quality rules spanning the full spectrum from ordered to chaotic, with complex rules naturally appearing along the phase boundary between them—and all of it discovered in a single run.

Multi-scale and hierarchical rules. Evolve rules that operate at multiple spatial or temporal scales simultaneously, where coarse-grained behavior feeds back into fine-grained dynamics. This approaches the hierarchical organization seen in biological systems, where molecular rules produce cellular behavior that produces tissue-level dynamics.

Transfer and generalization. Test whether rules evolved for the density classification problem also exhibit interesting behavior on other tasks or initial conditions they were never trained on. This probes whether evolution discovers general computational principles or narrow, task-specific hacks—a question with deep implications for our understanding of both artificial and biological evolution.

The undecidability frontier. Here is something that rarely appears in tutorials on this topic but deserves to. Because Rule 110 is Turing complete, questions about its long-term behavior are formally undecidable—no algorithm can determine, for arbitrary configurations in general, whether a given Rule 110 configuration will ever reach a particular state — this is a consequence of the halting problem, which applies to any Turing-complete system. Not every question about Rule 110 is undecidable, but the existence of undecidable questions is guaranteed. This undecidability propagates to evolutionary search: if you evolve a population toward Class IV rules, you are evolving into a region of rule space where the rules themselves are, in a precise mathematical sense, beyond analysis. You cannot predict what they will do; you can only run them. Evolution may be one of the only processes capable of navigating this territory, because it does not predict—it tries. This transforms the exercise from optimization into something closer to exploration of a mathematically unreachable landscape.

Each of these directions pushes deeper into the central question of this article: not "what can we build?" but "what can we grow?" The code above gives you a working evolutionary engine. The research literature gives you benchmarks and baselines. Everything else is unexplored territory.

The frontier is open. Run the code. Evolve some universes. See what emerges.