algorithm.md

Gillespie's Direct Method

1. Mathematical model

The implementation targets time-homogeneous continuous-time Markov chains (CTMCs) represented as reaction networks. The state is a vector of nonnegative integer counts,

$$ \mathbf{x}(t) \in \mathbb{N}_0^D, $$

and the model contains $M$ reaction channels. Reaction $j$ is defined by

• a stoichiometric increment $\boldsymbol{\nu}_j \in \mathbb{Z}^D$;
• a propensity $a_j(\mathbf{x}) \ge 0$.

Conditioned on the current state $\mathbf{x}$, the infinitesimal transition law is

$$ \Pr\left\lbrace \mathbf{X}(t+dt)=\mathbf{x}+\boldsymbol{\nu}_j \mid \mathbf{X}(t)=\mathbf{x} \right\rbrace =a_j(\mathbf{x})\thinspace dt+o(dt). $$

The total propensity is

$$ a_0(\mathbf{x})=\sum_{j=1}^{M}a_j(\mathbf{x}). $$

The public propensity interface is propensity(state). This restriction makes the time-homogeneity assumption explicit: propensities may depend on the current state and fixed model parameters, but not continuously on absolute time between jumps. Nonhomogeneous processes require an integrated-hazard or thinning method and are outside the scope of this implementation.

2. Distribution of the next event

While the state remains fixed, $a_0$ is constant. The survival probability for the waiting time $T$ is therefore

$$ \Pr(T>\tau\mid\mathbf{x})=\exp[-a_0(\mathbf{x})\tau], $$

so that

$$ T\mid\mathbf{x}\sim\mathrm{Exp}(a_0), \qquad \mathbb{E}[T\mid\mathbf{x}]=\frac{1}{a_0}. $$

Conditioned on an event occurring, its channel is categorical:

$$ \Pr(J=j\mid\mathbf{x})=\frac{a_j(\mathbf{x})}{a_0(\mathbf{x})}. $$

Equivalently, the joint density of waiting time and channel is

$$ p(\tau,j\mid\mathbf{x}) =a_j(\mathbf{x})\exp[-a_0(\mathbf{x})\tau]. $$

The implementation samples $T$ from NumPy's exponential distribution with scale $1/a_0$. Given $u\in[0,1)$, it selects the first channel $\mu$ satisfying

$$ \sum_{j=1}^{\mu}a_j(\mathbf{x})>u\thinspace a_0(\mathbf{x}). $$

The strict inequality matches NumPy's half-open uniform interval and ensures that a zero-propensity channel is not selected when $u=0$.

3. Simulation algorithm

For initial state $\mathbf{x}(t_0)$, initial time $t_0$, and horizon $t_{\mathrm{end}}$, the Direct Method is implemented as follows:

t <- t_start
x <- copy(initial_state)
record(t, x)

while t < t_end:
    if accepted_events == max_events:
        terminate with max_events

    a <- evaluate_propensities(x)
    require finite(a) and a >= 0
    a0 <- sum(a)

    if a0 == 0:
        terminate with absorbing_state

    tau <- Exponential(rate=a0)
    t_next <- t + tau

    if t_next > t_end:
        terminate with time_horizon

    mu <- categorical_channel(weights=a)
    x_next <- x + nu[mu]
    validate(x_next)

    t <- t_next
    x <- x_next
    record(t, x, mu)

An event with t_next == t_end is accepted. An event with t_next > t_end is not selected or applied.

Termination semantics

Reason	Condition	Final boundary row
time_horizon	The next event lies beyond t_end, or an accepted event reaches it	Present at t_end unless the last event is already there
absorbing_state	$a_0(\mathbf{x})=0$	The absorbing state is extended to t_end
max_events	The accepted-event limit is reached before t_end	Not added, because the uncomputed continuation is unknown

reaction_indices contains only accepted reactions. A synthetic row at t_end records the state at the boundary but does not increment SimulationResult.event_count.

4. Validation and numerical invariants

The engine enforces the following conditions:

• the state is a one-dimensional int64 vector of nonnegative counts;
• each stoichiometric increment has the same dimension as the state;
• species names and reaction names are nonempty and unique;
• every propensity is scalar, finite, and nonnegative;
• the total propensity is finite;
• accepted state updates remain nonnegative and do not overflow int64;
• event times are strictly increasing;
• optional model-specific invariants hold initially and after every event.

A zero total propensity is a valid absorbing state, not an exception. A negative or nonfinite propensity is a model error. If $1/a_0$ overflows, the engine terminates at the requested horizon; if a positive waiting time is too small to advance the floating-point clock, it raises FloatingPointError rather than entering a non-progressing loop.

For a trajectory with $K$ accepted events, the current implementation stores the complete state after every event. Its storage cost is $O(KD)$. With $M$ reaction channels, a Direct Method step costs $O(M)$ propensity evaluations and a linear channel selection; this design favors a compact, auditable kernel over dependency-graph optimizations.

5. Trajectory representation and resampling

SimulationResult.times and SimulationResult.states represent a right-continuous path with left limits. For recorded event times

$$ t_0<t_1<\dots<t_K, $$

the state on a sampling grid is defined by previous-value interpolation:

$$ k(q)=\max\lbrace i:t_i\le q\rbrace, \qquad \mathbf{x}(q)=\mathbf{x}_{k(q)}. $$

At an event time, this convention returns the post-reaction state. Between events, it returns the most recently accepted state. Linear interpolation is not appropriate for integer-valued jump processes, and nearest-neighbor sampling can introduce look-ahead bias by exposing a state before its event.

Ensemble statistics are computed only after all trajectories have been mapped to a common grid with this convention. The grid changes the representation of the paths, not the simulated event process.

6. Closed SIR model

The example model uses the state

$$ \mathbf{x}=(S,I,R)^\mathsf{T}, \qquad N=S+I+R, $$

with two reaction channels.

Infection

$$ S+I\longrightarrow 2I, \qquad \boldsymbol{\nu}_{\mathrm{inf}}=(-1,1,0)^\mathsf{T}, $$

$$ a_{\mathrm{inf}}(\mathbf{x})=\beta\frac{SI}{N}. $$

Recovery

$$ I\longrightarrow R, \qquad \boldsymbol{\nu}_{\mathrm{rec}}=(0,-1,1)^\mathsf{T}, $$

$$ a_{\mathrm{rec}}(\mathbf{x})=\gamma I. $$

Both reactions conserve $N$. The resulting paths satisfy

$$ S(t+dt)\le S(t), \qquad R(t+dt)\ge R(t), \qquad S(t)+I(t)+R(t)=N. $$

When $I=0$, both propensities vanish and the state is absorbing. The model is a closed, frequency-dependent SIR process: births, deaths, waning immunity, and external forcing are not included.

7. Ensembles and reproducibility

simulate constructs a run-local NumPy generator from the supplied seed; it does not use global numpy.random state. run_ensemble expands a master seed with numpy.random.SeedSequence and assigns one child seed to each run. The ordering of the returned trajectories is deterministic for fixed inputs and a fixed master seed.

Reproducible analysis requires retaining

• model parameters and initial state;
• t_start, t_end, and max_events;
• the run seed or ensemble master seed;
• the package, Python, and NumPy versions.

The result arrays do not encode arbitrary parameters captured by propensity callbacks. Consequently, summarize_ensemble verifies species and reaction catalogs but cannot detect different parameterizations with identical names. Only trajectories generated from the same configured model should be combined in one summary.

8. Scope and performance

This implementation provides the exact Direct Method for finite, time-homogeneous reaction networks. It does not implement tau-leaping, time-dependent hazards, delayed reactions, spatial processes, or deterministic integration. Those methods require different mathematical assumptions and separate numerical kernels.

The benchmark in benchmarks/benchmark_core.py reports median wall time and event throughput without imposing a hardware-dependent pass/fail threshold. It is intended for regression measurements under a controlled software and hardware configuration.

References

• D. T. Gillespie, “A general method for numerically simulating the stochastic time evolution of coupled chemical reactions”, Journal of Computational Physics, 1976.
• D. T. Gillespie, “Exact stochastic simulation of coupled chemical reactions”, The Journal of Physical Chemistry, 1977.