2025/08/06 by Petrack, Joshua, Doty, David
#Data Structures and Algorithms (cs.DS) #FOS: Computer and information sciences
paper · doi:10.48550/arxiv.2508.04079
The model of chemical reaction networks is among the oldest and most widely studied and used in natural science. The model describes reactions among abstract chemical species, for instance A + B → C, which indicates that if a molecule of type A interacts with a molecule of type B (the reactants), they may stick together to form a molecule of type C (the product). The standard algorithm for simulating (discrete, stochastic) chemical reaction networks is the Gillespie algorithm [JPC 1977], which stochastically simulates one reaction at a time, so to simulate ℓ consecutive reactions, it requires total running time Ω(ℓ). We give the first chemical reaction network stochastic simulation algorithm that can simulate ℓ reactions, provably preserving the exact stochastic dynamics (sampling from precisely the same distribution as the Gillespie algorithm), yet using time provably sublinear in ℓ. Under reasonable assumptions, our algorithm can simulate ℓ reactions among n total molecules in time O(ℓ/√ n) when ℓ ≥ n5/4, and in time O(ℓ/n2/5) when n ≤ ℓ ≤ n5/4. Our work adapts an algorithm of Berenbrink, Hammer, Kaaser, Meyer, Penschuck, and Tran [ESA 2020] for simulating the distributed computing model known as population protocols, extending it (in a very nontrivial way) to the more general chemical reaction network setting. We provide an implementation of our algorithm as a Python package, with the core logic implemented in Rust, with remarkably fast performance in practice.