Gillespie-based simulation and inference for non-Markovian stochastic reaction networks

Discrete stochastic processes are widespread across physics, chemistry, ecology, and beyond. In computational biology and epidemiology, however, most simulators still assume Markovian kinetics with memoryless dynamics, despite growing evidence for history-dependent effects in gene regulation, RNA transcription, cell differentiation, and infection. This reliance on Markovian models limits the routine use and comparison of more realistic, memory-aware descriptions. Here, we develop and benchmark a unified framework for simulating non-Markovian reaction networks using Gillespie-based algorithms. We implement multiple algorithmic classes, including exact, rejection-based, delay-based, and hybrid Markovian/non-Markovian schemes, and compare them across representative biological models. Across three case studies, we show that non-Markovian waiting times can qualitatively change population-level predictions, and that delay-based approximations can break down when intrinsic system timescales approach the imposed delays. We further show how population-level measurements can be used to infer otherwise inaccessible waiting-time distributions, and how non-Markovian structure can be leveraged for sensitivity analysis and statistical inference. To support broad use of these approaches, we release NoMaSS (Non-Markovian Stochastic Simulations), an open-source Python library that provides a unified interface to a wide range of non-Markovian Gillespie algorithms and enables hybrid Markovian/non-Markovian models (https://github.com/AI-SysBio/NoMaSS).