Introduction: From Single Species to Complex Ecosystems¶

In the previous lectures (particularly Lectures 35–37), a powerful theoretical framework for describing stochastic processes of single species has been systematically established. Whether through the probability-based Kramers-Moyal expansion or through the operator-algebra-based coherent-state path integral, the master equation describing discrete particle birth-death processes has been precisely transformed into a path integral paradigm that can utilize field-theoretic tools.

However, the vast majority of interesting phenomena in nature—whether species competition in ecosystems, disease transmission in epidemiology, or reaction networks in biochemistry—involve multiple interacting species. Therefore, while single-species theory provides a solid foundation, it must be generalized to describe the complexity of the real world.

The core objective of this lecture is to generalize the previously developed path integral formalism from single-species systems to more complex and realistic multi-species systems. Through this generalization, a unified stochastic dynamics theory will be constructed for systems containing arbitrary numbers of species and arbitrarily complex reaction networks.

To concretely demonstrate the power of this generalization, a classic and non-trivial model will be used as a case study: the Rock-Paper-Scissors model, also known in ecology as the May-Leonard model. This model is a paradigm for studying cyclic dominance and nonequilibrium dynamics. Despite its simple rules, it can exhibit complex spatiotemporal dynamical behaviors such as sustained population oscillations and spiral wave formation—rich phenomena that simple equilibrium statistical mechanics cannot describe.

This lecture also marks the transition of the course from constructing abstract mathematical tools to their direct application in complex nonequilibrium systems. The theoretical machinery carefully constructed in previous lectures (Lectures 31–37)—including the Janssen-De Dominicis response functional, Onsager-Machlup functional, and coherent-state path integral—will demonstrate their necessity and power when dealing with nonequilibrium problems with complex nonlinear interactions like the Rock-Paper-Scissors model.

1. Review: Coherent-State Path Integral Framework for Single-Species Systems¶

To generalize the theory from single-species to multi-species systems, we must first review the core concepts and mathematical structures established for single-species stochastic processes in previous lectures (particularly Lectures 36 and 37). This coherent-state path integral framework serves as a solid bridge from the master equation to field-theoretic analysis.

1.1 From Master Equation to Generating Function¶

The theoretical starting point for stochastic birth-death processes is the forward master equation. It describes the time evolution of the probability $p(n, t \mid n_0, t_0)$ that the system is in a state with $n$ particles at time $t$:

\[ \partial_t p(n, t \mid n_0, t_0) = \sum_m Q(n, m) p(m, t \mid n_0, t_0) \]

$p(n, t \mid n_0, t_0)$ is the conditional probability, i.e., the probability of finding $n$ particles at time $t$ given that there were $n_0$ particles at the initial time $t_0$.
$Q(n, m)$ is an element of the transition rate matrix, representing the probability per unit time for the system to transition from a state with $m$ particles to a state with $n$ particles. This equation is essentially a detailed accounting of probability flow between different particle number states.

Directly solving this infinite-dimensional system of coupled ordinary differential equations is usually very difficult. To simplify the problem, Lecture 36 introduced the coherent-state representation, whose core is defining a marginalized distribution function. It is essentially a generating function with respect to the initial particle number $n_0$:

\[ p(n, t \mid x, \tau) := \sum_{n_0=0}^{\infty} p(n, t \mid n_0, \tau) \frac{x^{n_0}}{n_0!} e^{-x} \]

This transformation converts a family of functions depending on discrete initial values $n_0$ into a single function depending on a continuous variable $x$ through the Poisson kernel function $\frac{x^{n_0}}{n_0!} e^{-x}$.
The initial time is denoted as $\tau$ here. It is important to emphasize that the variable $x$ is not a physical quantity (such as particle density or position) but a mathematical auxiliary variable used for generating moments.

The power of this transformation lies in converting the master equation from an infinite-dimensional system of ordinary differential equations into a single, more manageable partial differential equation (PDE):

\[ \partial_{-\tau} p(\cdot \mid x, \tau) = \mathcal{Q} p(\cdot \mid x, \tau) \]

where $\mathcal{Q}$ is the evolution operator in the coherent-state representation, commonly known as the Liouville operator.

The Liouville operator is named after the French mathematician Joseph Liouville, whose conceptual origins can be traced back to his research on conservation laws in classical mechanics phase space in the 1830s. However, the systematic formalization of this operator in modern theoretical physics was completed by John von Neumann and others in the quantum statistical mechanics established in the 1930s. Its physical essence is to describe the generator of time evolution of probability distributions or density matrices in phase space (classical) or Hilbert space (quantum): in classical mechanics, the Liouville operator L is defined through the Hamiltonian via Poisson brackets as $L \cdot = i \{H, \cdot\}$, generating the time evolution of phase space flow; in quantum mechanics, it is defined through the Hamiltonian via commutation relations as $L \hat{\rho} = [H, \hat{\rho}] / i\hbar$, governing the dynamics of the density matrix $\hat{\rho}$ (i.e., the von Neumann equation). Its core applications span various frontiers of statistical physics, including serving as the foundation of linear response theory (such as the Kubo formula), constructing the master equation framework for nonequilibrium statistical physics, analyzing eigenvalues of chaotic system dynamics (Pollicott-Ruelle resonances), and describing decoherence processes in quantum information and open system dynamics.

1.2 Single-Species Path Integral Action¶

Starting from this partial differential equation that is formally isomorphic to the (imaginary-time) Schrödinger equation, through standard path integral construction methods (such as time slicing), one can obtain the path integral representation of the system's evolution propagator. Its core is the action $\mathcal{S}$, which assigns a complex weight to every possible path of the system in the state space composed of auxiliary variables $x(\tau)$:

\[ \mathcal{S}[x(\tau), iq(\tau)] = \int_{t_{0}}^{t} d\tau \left[ iq(\tau) \partial_{\tau} x(\tau) - \mathcal{Q}(x, iq + 1) \right] \]

The structure of this action consists of two parts, which has a similar logical structure to the J-D action derived for continuous Langevin systems in Lecture 32:

Kinetic Term $iq \partial_\tau x$: This term is often called the "Berry phase" term in field theory. It does not depend on the specific reactions of the system but arises from the mathematical overlap between coherent-state basis vectors at different times. It ensures that the path integral correctly describes continuous time evolution and dynamically connects the field $x(\tau)$ with its conjugate response field $q(\tau)$.
Hamiltonian/Liouvillian $\mathcal{Q}(x, iq + 1)$: This term is the "classical" counterpart of the master equation evolution operator $\hat{\mathcal{Q}}$ in the path integral (mathematically called the "symbol" of the operator). It contains all the dynamical information of the system, i.e., all birth-death reactions and their rates. Its specific form is determined by the system's reaction network. For example, for a generic reaction $\sum kA \to \sum lA$, the corresponding Liouville operator part is:

$$ \mathcal{Q} = \tilde{\gamma} \left( (a^\dagger)^k (a^l - a^k) \right) $$ where $a^\dagger$ and $a$ are the creation and annihilation operators in the coherent-state representation, which are replaced by $x$ and $(iq+1)$ respectively in the path integral.

This derivation flow from "master equation → generating function → partial differential equation → path integral action" establishes a complete theoretical framework for single-species systems. The core task ahead is to show how this method can be systematically generalized to systems containing arbitrary numbers of species.

2. Generalization to Multi-Species Systems¶

After reviewing single-species path integrals, we can now begin the core generalization work. The key insight is that the mathematical framework of coherent-state path integrals has powerful universality. The transition from single-species to multi-species is essentially elevating scalars (such as particle numbers $n$, auxiliary fields $x$) in the theory to vectors, while keeping the basic structure of the action unchanged.

2.1 Multi-Species State Space and Operators¶

First, we need to extend the mathematical space describing the system state.

In multi-species systems, the microscopic state of the system is no longer described by a single integer $n$, but by a particle number vector $\boldsymbol{n} = (n_1, n_2, \ldots, n_S)$, where $n_\sigma$ is the number of particles of species $\sigma$.

Physical meaning: The system's state space thus expands from a one-dimensional integer axis to an $S$-dimensional discrete lattice space $\mathbb{N}_0^S$. Each point in this high-dimensional space represents a complete "census" of an ecosystem or chemical reaction network—a specific population composition.

Correspondingly, the coherent-state representation also needs to be generalized. Multi-species coherent states can be constructed as tensor products of single-species coherent states. This means each species has its own independent coherent-state basis, and the entire system's basis is their combination:

\[ |\boldsymbol{n}\rangle_\boldsymbol{x} = |n_1, n_2, \ldots, n_S\rangle_{\{x_1, x_2, \ldots, x_S\}} = |n_1\rangle_{x_1} \otimes \cdots \otimes |n_S\rangle_{x_s} \]

$\boldsymbol{x} = (x_1, x_2, \ldots, x_S)$ is a vector containing $S$ formal variables.
The tensor product $\otimes$ is the natural mathematical language for combining independent spaces, meaning that when interactions are not considered, each species' state is independent.

Creation and annihilation operators now also need to carry species labels $\sigma$. Their action on multi-species Poisson bases is defined as:

\[ a^\dagger_\sigma = x_\sigma, \quad a_\sigma = \partial_{x_\sigma} + 1 \]

Physical meaning: $a^\dagger_\sigma$ acts to add a particle of species $\sigma$ to the system (equivalent to multiplying by $x_\sigma$ in the Poisson representation), while $a_\sigma$ acts to remove a particle of species $\sigma$ (equivalent to acting with a differential operator in the Poisson representation).

These operators only act on their corresponding species' subspace. Operators belonging to different species commute, e.g., $[a_\sigma, a_\tau^\dagger] = \delta_{\sigma\tau}$, which again reflects the independence of each species' degrees of freedom.

2.2 Multi-Species Action¶

With vectorized state spaces and operators, following the blueprint of the single-species case, we can directly write down the path integral action for multi-species systems. Its form is remarkably similar to the single-species case, with all scalars replaced by vectors:

\[ S[\boldsymbol{x}(\tau), i\boldsymbol{q}(\tau)] = \int_{t_0}^{t} d\tau \left[ i\boldsymbol{q}(\tau)^\top \partial_\tau \boldsymbol{x}(\tau) - \mathcal{Q}(\boldsymbol{x}, i\boldsymbol{q} + 1) \right] \]

$\boldsymbol{x}(\tau) = (x_1(\tau), \ldots, x_S(\tau))$ and $\boldsymbol{q}(\tau) = (q_1(\tau), \ldots, q_S(\tau))$ are now $S$-dimensional vector fields evolving in time.

$i\boldsymbol{q}^\top \partial_\tau \boldsymbol{x}$ is the inner product of vectors, i.e., $\sum_{\sigma=1}^S i q_\sigma(\tau) \partial_\tau x_\sigma(\tau)$, representing the sum of kinetic terms for all species.

$\mathcal{Q}(\boldsymbol{x}, i\boldsymbol{q} + 1)$ is the multi-species Liouville quantity, which is now a function of $S$ pairs of variables $(\boldsymbol{x}, \boldsymbol{q})$ and contains all information about inter-species interactions.

The transition from scalars to vectors is not merely a notational change; it embodies profound physical and geometric imagery. The state space of a single-species system is one-dimensional. For a three-species system (like Rock-Paper-Scissors), the state space is a three-dimensional population lattice space $\mathbb{N}_0^3$. The essence of the path integral is to perform weighted summation over all possible random trajectories $\boldsymbol{n}(t)$ traversing this high-dimensional "population space." The fields $\boldsymbol{x}(\tau)$ and $\boldsymbol{q}(\tau)$ in the coherent-state path integral are precisely the continuous field-theoretic counterparts of these discrete trajectories. Therefore, the action $S[\boldsymbol{x}, i\boldsymbol{q}]$ describes the "cost" or probability weight of a specific evolutionary path in this high-dimensional population space. This geometric intuition provides powerful support for understanding the abstract field-theoretic formalism.

This generalization reveals a core principle: All the complexity of multi-species interactions is elegantly encapsulated in the specific form of the multi-species Liouville operator $\mathcal{Q}$, while the basic framework of the path integral itself remains unchanged. This means that we don't need to reinvent a theory for each new multi-species system; we only need to learn how to correctly "construct" the corresponding Liouville operator for a given reaction network. This is precisely the core task of the next section.

3. Constructing Liouville Operators for Chemical Reaction Networks¶

In the previous section, the theoretical framework has been successfully generalized from single-species scalar form to multi-species vector form. This generalization reveals a core principle: all the complexity of multi-species interactions is elegantly encapsulated in the specific form of the multi-species Liouville operator $\mathcal{Q}$.

The question now is: how to systematically construct this crucial Liouville operator $\mathcal{Q}$ based on a given reaction network? This lecture will provide a universal "recipe" or "algorithm" that can precisely translate any physical, chemical, or biological process (manifested as a series of reaction rules) into a mathematical object (i.e., the Liouville operator) that generates its complete stochastic dynamics.

Consider a general chemical reaction where $k_\sigma$ molecules of species $A_\sigma$ as reactants are converted to $l_\sigma$ molecules of species $A_\sigma$ as products at rate $\gamma$:

\[ \sum_{\sigma=1}^{S} k_\sigma A_\sigma \xrightarrow{\gamma} \sum_{\sigma=1}^{S} l_\sigma A_\sigma \]

3.1 From Reaction Rules to Transition Matrices¶

Before constructing the operator, we first review how the transition rate matrix $Q(\boldsymbol{n}, \boldsymbol{m})$ in the master equation is generated from basic reaction rules.

Suppose the system's current state is $\boldsymbol{m} = (m_1, \ldots, m_S)$. For a reaction to occur, we need to select $k_1$ from $m_1$ $A_1$ particles, $k_2$ from $m_2$ $A_2$ particles, and so on. The total reaction rate is proportional to the product of these combinatorial ways. After the reaction occurs, the system state jumps from $\boldsymbol{m}$ to $\boldsymbol{n} = \boldsymbol{m} - \boldsymbol{k} + \boldsymbol{l}$. Combined, the transition matrix element has the form:

\[ Q(\boldsymbol{n}, \boldsymbol{m}) = \tilde{\gamma} (m_1)_{k_1} \cdots (m_S)_{k_S} \left[ \delta_{\boldsymbol{n}, \boldsymbol{m}-\boldsymbol{k}+\boldsymbol{l}} - \delta_{\boldsymbol{n},\boldsymbol{m}} \right] \]

$(m_\sigma)_{k_\sigma} = m_\sigma(m_\sigma - 1)\cdots(m_\sigma - k_\sigma + 1)$ is the falling factorial, representing combinatorial selection numbers.

$\tilde{\gamma} = \gamma / (\prod_\sigma k_\sigma!)$ is the rate constant including combinatorial symmetry factors.

The first term $\delta_{\boldsymbol{n}, \boldsymbol{m}-\boldsymbol{k}+\boldsymbol{l}}$ in $[...]$ describes the "gain" process where the state jumps into the new state $\boldsymbol{n}$ from $\boldsymbol{m}$.

The second term $-\delta_{\boldsymbol{n},\boldsymbol{m}}$ describes the "loss" process where probability flows out from the original state $\boldsymbol{m}$.

3.2 General Form of the Liouville Operator¶

The crucial step is to translate the above discrete transition matrix into an operator form that acts on continuous functions in the coherent-state representation. By utilizing the algebraic properties of creation operators $a^\dagger_\sigma$ and annihilation operators $a_\sigma$, we can obtain an extremely elegant general formula:

\[ \hat{\mathcal{Q}} = \tilde{\gamma} \left( (a^\dagger_1)^{k_1} \cdots (a^\dagger_S)^{k_S} \right) \left( a_1^{l_1} \cdots a_S^{l_S} - a_1^{k_1} \cdots a_S^{k_S} \right) \]

The structure of this formula can be intuitively interpreted as: Rate Factor × (Jump Operation).

Rate Factor: $\tilde{\gamma} (a^\dagger_1)^{k_1} \cdots (a^\dagger_S)^{k_S}$

This part corresponds to the macroscopic rate of the reaction. In the coherent-state representation, creation operators $a^\dagger_\sigma$ ultimately correspond to auxiliary fields $x_\sigma$ in the path integral. Therefore, this part becomes $\tilde{\gamma} x_1^{k_1} \cdots x_S^{k_S}$ in the action, which precisely reproduces the law of mass action in deterministic chemical kinetics. It shows that the total propensity for a reaction to occur is proportional to the power product of reactant concentrations.

Jump Operation: $(a_1^{l_1} \cdots a_S^{l_S} - a_1^{k_1} \cdots a_S^{k_S})$

This part is the core of executing the reaction, precisely encoding the change in the particle number vector. It first annihilates the state composed of reactants through the term $-a_1^{k_1} \cdots a_S^{k_S}$, then creates the state composed of products through the term $a_1^{l_1} \cdots a_S^{l_S}$. This "subtract-add" structure corresponds to the "loss-gain" form in the master equation.

This formula turns the construction of Liouville operators into a "universal compiler." Any reaction list based on empirical observations (the "source code" of physics) can be mechanically and unambiguously translated into a unique mathematical operator (the "executable file" of dynamics) through this formula, and this operator generates all the stochastic dynamics of the system.

3.3 Case Analysis: Translating Rules into Operators¶

Through several simple examples given by the professor in class, we can better understand this construction process.

Death: $A \xrightarrow{\mu} \emptyset$

Identification: Species A: $k_A = 1$, $l_A = 0$. Rate $\tilde{\gamma} = \mu$.

Apply formula: $\hat{\mathcal{Q}}_{\text{death}} = \mu (a^\dagger)^1 (a^0 - a^1) = \mu a^\dagger (1 - a)$.

Birth: $B \xrightarrow{\sigma} B + B$

Identification: Species B: $k_B = 1$, $l_B = 2$. Rate $\tilde{\gamma} = \sigma$.

Apply formula: $\hat{\mathcal{Q}}_{\text{birth}} = \sigma (b^\dagger)^1 (b^2 - b^1) = \sigma b^\dagger (b^2 - b)$.

Predation: $A + B \xrightarrow{\lambda} A + A$

Identification: Species A: $k_A = 1, l_A = 2$. Species B: $k_B = 1, l_B = 0$. Rate $\tilde{\gamma} = \lambda$.

Apply formula: $\hat{\mathcal{Q}}_{\text{pred}} = \lambda (a^\dagger)^1 (b^\dagger)^1 (a^2 b^0 - a^1 b^1) = \lambda a^\dagger b^\dagger (a^2 - ab)$.

For a complex system containing multiple reactions, the total Liouville operator is simply the sum of all individual reaction operators: $$\hat{\mathcal{Q}}_{\text{total}} = \sum_{\text{all reactions}} \hat{\mathcal{Q}}_{\text{reaction}}$$ This modular property makes the framework extremely powerful, capable of systematically handling arbitrarily complex reaction networks, laying the foundation for analyzing the "Rock-Paper-Scissors" model in the next section.

4. Rock-Paper-Scissors Dynamics of the May-Leonard Model¶

In the previous section, a powerful "compiler" was established that can translate any reaction network into its corresponding Liouville operator. Now it is time to apply this formalism to a concrete, non-trivial ecological model: the May-Leonard model.

The May-Leonard model was jointly proposed by ecologists Robert May and Joseph Leonard in 1975, with its historical background stemming from reflection on the "competitive exclusion principle"—that species with identical ecological niches cannot stably coexist. The physical essence of this model is to construct a three-dimensional (or higher-dimensional) nonlinear dynamical system (system of ordinary differential equations) describing the cyclic competitive relationships between three species (i.e., A beats B, B beats C, C beats A, forming "Rock-Paper-Scissors"-style intransitive competition), whose core feature is that under the deterministic differential equation framework, due to the absence of stable fixed points, the system exhibits nonequilibrium behaviors such as sustained oscillations and spontaneous random alternation of species abundances. The specific applications of this model have far exceeded theoretical ecology, becoming a paradigmatic tool for explaining biodiversity and spatial self-organization patterns, widely applied in microbial experimental systems (such as cyclic competition in E. coli), coral reef ecosystem evolution analysis, and mechanism research for maintaining cooperative behavior in evolutionary game theory, providing a key mathematical model for nonlinear dynamics and nonequilibrium statistical physics to demonstrate how neutral stability, limit cycles, and stochasticity jointly maintain biodiversity.

This simple game rule—rock beats scissors, scissors beat paper, paper beats rock—has widespread analogies in nature. For example, in certain E. coli strains, some strains produce toxins (rock) that kill sensitive strains (scissors); sensitive strains grow faster and can outcompete toxin-producing but slow-growing strains (paper); while a third strain is immune to toxins but grows slower than sensitive strains (paper beats rock).

The May-Leonard model is renowned for producing rich, far-from-equilibrium spatiotemporal structures and is an ideal "touchstone" for testing stochastic population dynamics theory.

The model contains three species A, B, C (analogous to rock, paper, scissors), whose interactions follow cyclic dominance rules. According to the PPT shown by the professor in class, these interactions can be divided into three categories:

Competition (Contest Competition / Predation): Predator meets prey, prey is removed.
- $A + B \xrightarrow{\sigma} A + \emptyset$ (A preys on B)
- $B + C \xrightarrow{\sigma} B + \emptyset$ (B preys on C)
- $C + A \xrightarrow{\sigma} C + \emptyset$ (C preys on A)
Reproduction (Scramble Competition / Reproduction): Individuals produce offspring on available space (empty sites $\emptyset$).
- $A + \emptyset \xrightarrow{\mu} A + A$
- $B + \emptyset \xrightarrow{\mu} B + B$
- $C + \emptyset \xrightarrow{\mu} C + C$
Mobility (Mobility / Diffusion): Individuals move in space, exchanging positions.
- $A_i \leftrightarrow A_j$ (where $i, j$ represent spatial lattice sites)

4.1 Mean-Field Dynamics: The Deterministic "Death Spiral"¶

Before delving into the complete stochastic dynamics, we first analyze a simplified version: a mean-field or well-mixed system. This approximation assumes that particles in the system can interact instantaneously with any other particle, thus ignoring spatial structure and random fluctuations.

According to the law of mass action, we can write the deterministic rate equations for population densities $a, b, c$:

\[ \partial_t a = a[\mu(1 - \rho) - \sigma c] \\ \partial_t b = b[\mu(1 - \rho) - \sigma a] \\ \partial_t c = c[\mu(1 - \rho) - \sigma b] \]

$\rho = a+b+c$ is the total density.
$a\mu(1 - \rho)$ represents the logistic growth of species A in available space $(1 - \rho)$ (reproduction term).
$-\sigma a c$ represents species A being consumed by its predator C (competition term).

Performing stability analysis on this dynamical system, we can find several fixed points:

Complete extinction fixed point: $(0, 0, 0)$.
Single-species absorbing fixed points: $(1, 0, 0), (0, 1, 0), (0, 0, 1)$. These are stable fixed points.
Coexistence fixed point: There exists a three-species coexistence fixed point at the center of phase space, but it is unstable.

Physical meaning: Heteroclinic cycle¶

Since the central coexistence point is unstable, any small perturbation will cause the system state to leave this point. The system's phase space trajectories exhibit outwardly diverging spiral lines, a structure called a heteroclinic cycle. Physically, this means that population numbers will experience oscillations with increasing amplitude: an increase in A's numbers leads to a decrease in its prey B; B's decrease reduces C's natural enemies, so C begins to increase; C's increase leads to A's decrease, completing a cycle.

In the deterministic mean-field description, this oscillation will continue indefinitely with increasing amplitude, and the trajectory will asymptotically approach the boundary composed of three absorbing states, leading to the extinction of two species and leaving only one winner. In finite but large populations, the time scale for such random extinction can be estimated as $T \sim \ln N$.

4.2 Spatial Stochastic Dynamics: Emergence of Spiral Waves and Species Coexistence¶

Mean-field theory (Section 4.1) depicts the "skeleton" of the system—a deterministic "death spiral" leading to species extinction. However, this theory has two fatal simplifications: it ignores demographic noise (random birth-death events in finite populations) and spatial structure (locality of interactions).

To explore how these two factors fundamentally change the system's fate, we need to upgrade from a 0-dimensional "well-mixed" model to a 2-dimensional spatial stochastic model. We will no longer track the total species numbers of the entire system, but simulate on a two-dimensional grid where each lattice site is occupied by a certain species and how they react with their neighbors. We use a cellular automaton approach that captures the dynamics of such local interactions.

Spatial Simulation Method and Code¶

Unlike the Gillespie algorithm that tracks the next reaction of the entire system, cellular automata update each lattice site's state in discrete time steps based on the states of its surrounding neighbors. The core rules are as follows:

Initialization: Create a two-dimensional grid and randomly place a species (A, B, C) or leave empty ($\emptyset$) at each lattice site.
Evolution: At each time step, calculate state changes for all lattice sites on the grid simultaneously:
- Predation: The more a lattice site is surrounded by its "natural enemy" neighbors, the higher the probability it becomes empty in this step.
- Reproduction: The more a lattice site is surrounded by neighbors of a certain species, the higher the probability it is occupied by that species in this step.
- Mobility: Individuals on lattice sites have a certain probability of exchanging positions with neighboring empty sites.
Update: Apply all calculated changes to a new grid, completing one time step of evolution.

The following Python code implements this spatial model and has been optimized to stably reproduce the spiral wave animations shown in class.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.colors import ListedColormap
from scipy.signal import convolve2d

def simulate_spatial_rps_spirals(grid_size, rates, steps):
    """
    Simulates the spatial Rock-Paper-Scissors model.

    Args:
        grid_size (int): The width and height of the grid.
        rates (dict): A dictionary of probabilities {'repro', 'pred', 'mobil'}.
        steps (int): The number of animation frames to generate.

    Returns:
        list: A history of the grid states for animation.
    """
    # Grid states: 0=Empty, 1=A (Rock/Red), 2=B (Paper/Blue), 3=C (Scissors/Yellow)
    grid = np.random.randint(0, 4, size=(grid_size, grid_size))
    history = [grid.copy()]

    # Kernel for counting neighbors
    kernel = np.array([[1, 1, 1],
                       [1, 0, 1],
                       [1, 1, 1]])

    for step in range(steps):
        # Create a copy to store the changes for this step
        new_grid = grid.copy()

        # --- 1. Reaction Phase ---
        # Count neighbors for each species
        neighbors_A = convolve2d((grid == 1), kernel, mode='same', boundary='wrap')
        neighbors_B = convolve2d((grid == 2), kernel, mode='same', boundary='wrap')
        neighbors_C = convolve2d((grid == 3), kernel, mode='same', boundary='wrap')

        # Identify predators for each species
        # Predators: C(3) for A(1), A(1) for B(2), B(2) for C(3)
        predators_for_A = neighbors_C
        predators_for_B = neighbors_A
        predators_for_C = neighbors_B

        # Calculate probabilities of events for each cell
        # Predation happens to non-empty cells
        prob_predation_A = rates['pred'] * predators_for_A
        prob_predation_B = rates['pred'] * predators_for_B
        prob_predation_C = rates['pred'] * predators_for_C

        # Reproduction happens in empty cells
        prob_repro_A = rates['repro'] * neighbors_A
        prob_repro_B = rates['repro'] * neighbors_B
        prob_repro_C = rates['repro'] * neighbors_C

        # Generate random numbers for all cells at once
        rand_field = np.random.rand(grid_size, grid_size)

        # Apply predation rules
        new_grid[(grid == 1) & (rand_field < prob_predation_A)] = 0
        new_grid[(grid == 2) & (rand_field < prob_predation_B)] = 0
        new_grid[(grid == 3) & (rand_field < prob_predation_C)] = 0

        # Apply reproduction rules (only to cells that are still empty)
        # Find empty cells and rank reproduction probabilities
        empty_mask = (new_grid == 0)
        repro_chances = np.stack([prob_repro_A, prob_repro_B, prob_repro_C], axis=-1)

        # The species with the most neighbors has the best chance to reproduce
        colonizer = np.argmax(repro_chances, axis=-1) + 1

        # Check if the max probability is > 0 and beats a random threshold
        max_prob = np.max(repro_chances, axis=-1)
        new_grid[empty_mask & (rand_field < max_prob)] = colonizer[empty_mask & (rand_field < max_prob)]

        # --- 2. Mobility Phase ---
        # A fraction of individuals try to move into empty spaces
        for _ in range(int(grid_size * grid_size * rates['mobil'])):
            r, c = np.random.randint(0, grid_size, 2)
            # Only non-empty sites can move
            if new_grid[r, c] != 0:
                # Pick a random neighbor
                dr, dc = [(0,1), (0,-1), (1,0), (-1,0)][np.random.randint(4)]
                nr, nc = (r + dr) % grid_size, (c + dc) % grid_size
                # If neighbor is empty, swap
                if new_grid[nr, nc] == 0:
                    new_grid[r, c], new_grid[nr, nc] = new_grid[nr, nc], new_grid[r, c]

        grid = new_grid
        history.append(grid.copy())

        if (step + 1) % 50 == 0:
            print(f"Step {step + 1}/{steps} completed.")

    return history

# --- Simulation Parameters ---
GRID_SIZE = 150
# These rates are chosen to favor spiral formation
RATES = {'repro': 0.1, 'pred': 0.3, 'mobil': 0.05}
ANIMATION_STEPS = 500

# --- Run Simulation ---
print("Starting spatial simulation for spirals...")
simulation_history = simulate_spatial_rps_spirals(GRID_SIZE, RATES, ANIMATION_STEPS)
print("Simulation finished. Creating animation...")

# --- Create and Save Animation ---
fig, ax = plt.subplots(figsize=(8, 8))
cmap = ListedColormap(['black', 'red', 'blue', 'yellow'])
im = ax.imshow(simulation_history[0], cmap=cmap, interpolation='nearest', vmin=0, vmax=3)
ax.set_xticks([])
ax.set_yticks([])
title = ax.set_title("Spatial Rock-Paper-Scissors (Step 0)")

def update(frame):
    im.set_array(simulation_history[frame])
    title.set_text(f"Spatial Rock-Paper-Scissors (Step {frame})")
    return [im, title]

ani = animation.FuncAnimation(fig, update, frames=len(simulation_history), interval=30, blit=True)
ani.save('spatial_rps_spirals.gif', writer='pillow', fps=30)
print("Animation saved as 'spatial_rps_spirals.gif'")
plt.show()

From chaos to order: The simulation starts from a randomly distributed "salt-and-pepper noise" state. Within a short time, species begin to form small patches of aggregation.
Formation of chasing wave fronts: At the boundaries between different species aggregation regions, predation reactions begin to occur. One can see a red "wave front" beginning to erode blue regions, while blue wave fronts erode yellow regions, and yellow wave fronts in turn erode red regions, forming a "chase-escape" dynamic chain.
Emergence of spiral waves: Due to random fluctuations, these straight wave fronts become unstable and begin to curl. Eventually, they self-organize into magnificent, stably propagating spiral waves. The center of the spiral wave acts like a pacemaker, continuously emitting cyclic species waves outward.

Unlike mean-field and well-mixed models that ultimately lead to species extinction, spatial models achieve stable species coexistence through spiral wave formation. Why does space have such magic?

Locality creates "refuges": In spatial models, predators can only attack their direct neighbors. This means that species groups far from their natural enemies can safely exist in "refuges." A species does not face "encirclement" from all natural enemies in the entire system.
Spiral waves are dynamic refuges: Spiral wave structures are dynamic, self-sustaining refuges. Wave propagation ensures that no species can occupy the entire space, and there are always prey that can "escape" to regions not yet reached by predators, thus providing opportunities for species continuation.

Therefore, the spatial dimension is not a dispensable detail; it fundamentally changes the system's dynamical behavior and ultimate fate, transforming a system destined to lose diversity into a complex ecosystem that can maintain dynamic balance and species coexistence.

In the next section, we will formally introduce the spatial dimension into our theoretical framework, laying the foundation for understanding the formation mechanisms of beautiful spatiotemporal patterns like spiral waves observed in simulations.

5. Introducing Space into Theory: Diffusion and Locality¶

In the numerical simulations of Section 4.2, we saw that when spatial dimensions were introduced, the three species that were destined for extinction in "well-mixed" systems achieved stable coexistence through spiral wave formation. This indicates that space is not a dispensable detail but a key factor that fundamentally changes the system's macroscopic behavior and ultimate fate.

Therefore, to enable path integral theory to describe these rich spatiotemporal phenomena, we must formally incorporate spatial structure and local interactions into the model. The goal of this lecture is to extend the 0-dimensional theory that only considers total particle numbers into a complete theory defined on spatial lattices.

5.1 Lattice Discretization and New State Description¶

To incorporate spatial structure into the model, we define dynamics on a (e.g., two-dimensional square) lattice.

The system's state description also changes accordingly. The state is no longer a single particle number vector $\boldsymbol{n}$, but is jointly described by the number of particles $n_i^\sigma$ of each species $\sigma$ at each lattice site $i$.

Physical meaning:

The system's complete state is now a very high-dimensional vector containing information about particle numbers of all species at all positions—a snapshot of the entire "ecological map."
Previously defined "on-site" reactions (such as predation, reproduction) are now understood as local events occurring within the same lattice site $i$ or between adjacent lattice sites $i$ and $j$.

5.2 Formalizing Diffusion as "Chemical Reactions"¶

In addition to local reactions, we need to introduce the process of particles moving between different lattice sites, i.e., diffusion or mobility. In our theoretical framework, diffusion can be elegantly formalized as another type of "chemical reaction"—namely, a particle of species $\sigma$ "jumping" from lattice site $i$ to its neighboring lattice site $j$:

\[ A_i^\sigma \xrightarrow{D_{i \to j}^\sigma} A_j^\sigma \]

$A_i^\sigma$ represents a particle of species $\sigma$ located at lattice site $i$.
$D_{i \to j}^\sigma$ is the hopping rate for this jumping event. The specific form of this rate depends on our assumptions about the system's physical properties.

Different Choices for Hopping Rate $D_{i \to j}^\sigma$¶

The choice of hopping rate profoundly reflects whether the system is in equilibrium or nonequilibrium.

Simplest nonequilibrium model: In many studies of nonequilibrium models (including the May-Leonard model), the simplest assumption is that the hopping rate is a constant $D$, independent of species, position, and direction. $$D_{i \to j}^{\sigma} = D$$ This assumption means particles perform unbiased, pure random walks.
Equilibrium model (satisfying detailed balance):

If the system has an energy landscape (e.g., each lattice site $i$ has a potential energy $U_i$), and we require the system to relax to thermodynamic equilibrium (i.e., Boltzmann distribution) after long-time evolution, then the hopping rate must satisfy the detailed balance condition.

Detailed balance requires that the forward probability flow from $j$ to $i$ equals the reverse probability flow from $i$ to $j$, i.e., $D_{ij} e^{-\beta U_j} = D_{ji} e^{-\beta U_i}$. There are various choices for hopping rates satisfying this condition, for example:

\[D_{ij} = D \exp\left[-\frac{1}{2}\beta(U_i - U_j)\right]\]

This form ensures that particles are more likely to jump from high-energy positions to low-energy positions, enabling the system to eventually reach the lowest-energy equilibrium state.

Important distinction: It must be emphasized that the May-Leonard model, due to its cyclic predation reaction rules, is an intrinsically nonequilibrium system that does not have a minimizable "energy" function, so its dynamics does not obey detailed balance. Therefore, when studying its nonequilibrium properties, adopting a simple constant hopping rate $D$ is a reasonable and effective starting point.

5.3 How Space Changes System Fate: Refuges and Pattern Formation¶

Introducing spatial dimensions and local interactions fundamentally changes the system's behavior.

In well-mixed models, the dominant species can indiscriminately attack all prey in the system, quickly leading to global population collapse and species extinction.
But in spatial models, an individual can only reproduce or prey with its neighboring individuals.

This locality allows prey to find "refuges" in regions far from predator populations. A species' patches can flourish in regions not yet reached by their natural enemies.

The combination of such local interactions with diffusion ultimately prevents the global extinction seen in well-mixed systems and gives rise to rich self-organized spatiotemporal patterns, such as the spiral waves we saw in the simulations of Section 4.2. In these spiral waves, the three species dynamically coexist on a global scale, although they are constantly replacing each other on a local scale.

Therefore, the transition from well-mixed models to spatial models is not merely a detail correction but a qualitative change in the system's dynamical behavior. After understanding how to formalize spatial dynamics, we have gathered all the necessary ingredients to construct a complete spatiotemporal path integral.

Conclusion¶

This lecture successfully generalized the path integral formalism from describing single-species stochastic processes to handling arbitrarily complex multi-species interaction systems. The core of this generalization lies in elevating scalar state variables and operators to vector form while keeping the basic structure of the action unchanged.

Through this process, a universal, algorithmic procedure was established for constructing the corresponding Liouville operator based on any given reaction network. This "compiler"-style tool directly translates reaction rules from physics, chemistry, or biology into mathematical generators of stochastic dynamics, greatly enhancing the practicality and universality of this theoretical framework.

Through in-depth analysis of the May-Leonard (Rock-Paper-Scissors) model, the powerful capabilities of this theoretical framework were fully demonstrated. Different levels of theory reveal dramatically different system fates:

Mean-field theory reveals the deterministic skeleton of the system—a heteroclinic cycle leading to species extinction.
Spatial stochastic simulation shows a completely opposite picture: local interactions and demographic noise jointly promote the emergence of spiral waves, enabling stable coexistence of species.

This stark contrast clearly demonstrates that for understanding the fate of such nonequilibrium systems, a complete stochastic description that includes spatial structure and random fluctuations is indispensable. In summary, this lecture not only provides a powerful theoretical tool—multi-species path integrals—but also deepens understanding of the roles of stochasticity, interactions, and spatial structure in complex systems through a classic nonequilibrium model.

In this lecture, all necessary physical ingredients—local reactions and diffusion between lattice sites—have been formally defined on discrete lattices. This naturally leads to the next question: How to construct a complete path integral on this discrete lattice system and derive macroscopic field-theoretic descriptions from it?

Next lecture (Lecture 39) "Path Integrals on Lattices: From Hopping to Continuous Field Theory" will directly answer this question. The multi-species action developed in this lecture will be combined with the concept of spatial lattices to construct a path integral describing spatiotemporal stochastic processes.

Furthermore, the course will explore the transition from this discrete, particle-based description to a continuous, field-based description. By performing a continuum limit—i.e., coarse-graining lattice spacing and particle numbers—we will ultimately derive stochastic partial differential equations (SPDEs) describing spatiotemporal fluctuations of population density fields and their corresponding field theory.

This will complete the entire theoretical construction journey from microscopic particle rules (master equation) to macroscopic continuous fields (stochastic field theory), opening the door for using powerful modern field-theoretic tools like renormalization groups to analyze critical behavior and phase transitions in nonequilibrium systems.