Simulation and inference methods for non-Markovian stochastic reaction networks

Stochastic models of reaction networks are widely used to capture intrinsic noise in complex systems in the life sciences. Typical formulations of these models are based on Markov processes for which there is extensive research on efficient simulation and inference. However, there are complex processes in biology, such as gene transcription and translation, that introduce history dependent dynamics requiring non-Markovian processes to accurately capture the stochastic dynamics of the system. This greater realism comes with additional computational challenges for simulation and parameter inference. We develop efficient stochastic simulation algorithms for well-mixed non-Markovian stochastic reaction networks with stochastic delays that depend on system state and time. Our methods generalize the next reaction method and τ-leaping method to support arbitrary inter-event time distributions while preserving computational scalability. We also introduce a coupling scheme to generate exact non-Markovian sample paths that are positively correlated to an approximate non-Markovian τ-leaping sample path. This enables substantial computational gains for simulation and Bayesian inference through multilevel Monte Carlo and multifidelity schemes. We demonstrate the effectiveness of our approach using several non-Markovian examples, showing substantial gains in both simulation accuracy and inference efficiency. These results extend the practical applicability of non-Markovian models in systems biology and beyond.Author summary: To understand and predict how cells within living organisms behave we need mathematical models of biochemical processes that govern cell function. Accurate models of these biochemical processes often require inherent internal randomness to be explicitly described. Unfortunately this leads to substantial computational challenges when simulating these models and analysing them statistically. These challenges are especially true when the process is non-Markovian: that is, the history of the process is important to determine the future. We develop new and efficient methods for simulation and statistical inference of such non-Markovian models and demonstrate that our approach achieves a substantial reduction in computation cost without sacrificing accuracy.