Exactly simulating stochastic chemical reaction networks in sub-constant time per reaction

2025-08-06

AI Summary

This work addresses the fundamental problem of sublinear-time exact simulation for stochastic chemical reaction networks: achieving, for the first time, strictly exact sampling—provably preserving the exact Gillespie algorithm’s stochastic distribution—while attaining subconstant amortized time per reaction. We introduce a novel algorithm grounded in the population protocol model and establish a theoretical time complexity of $O(ell / sqrt{n})$ for simulating $ell$ consecutive reactions, which asymptotically improves upon the classical $Omega(ell)$ linear lower bound when $ell geq n^{5/4}$. The algorithm is implemented via a high-performance Rust core with a user-friendly Python interface, enabling efficient large-scale network simulation. We release an open-source Python package that delivers provably exact, reproducible sublinear-time stochastic reaction sampling and event scheduling. This constitutes the first solution that simultaneously achieves theoretical optimality and practical engineering viability for computational systems biology and large-scale biochemical modeling.

Technology Category

• Search and Optimization: Sampling/Simulation-based Search
• Reasoning under Uncertainty: Stochastic Optimization
• Planning, Routing, and Scheduling: Optimization of Spatio-temporal Systems

Application Category

• Graph Algorithms and Modeling for the Web: Foundation models and LLMs for Web-related graphs
• Systems and Infrastructure for Web, Mobile and WoT: Experiences and lessons learnt from Web-based algorithms and system deployments
• Web Mining and Content Analysis: Web data generation and simulation

Abstract

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 o 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 $ell$ consecutive reactions, it requires total running time $Ω(ell)$. We give the first chemical reaction network stochastic simulation algorithm that can simulate $ell$ reactions, provably preserving the exact stochastic dynamics (sampling from precisely the same distribution as the Gillespie algorithm), yet using time provably sublinear in $ell$. Under reasonable assumptions, our algorithm can simulate $ell$ reactions among $n$ total molecules in time $O(ell/sqrt n)$ when $ell ge n^{5/4}$, and in time $O(ell/n^{2/5})$ when $n le ell le n^{5/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.

Problem

• Simulating stochastic chemical reactions in sublinear time
• Improving Gillespie algorithm's linear time complexity
• Enabling faster exact simulation for large reaction networks

Innovation

• Sublinear time stochastic simulation algorithm
• Exact dynamics preservation like Gillespie
• Adapted population protocol algorithm extension