🧬 biology

Simulation and inference methods for non-Markovian stochastic biochemical reaction networks

This paper develops efficient stochastic simulation and inference methods for non-Markovian biochemical reaction networks by generalizing the next reaction and $\tau$ -leaping algorithms to handle arbitrary inter-event time distributions and introducing a coupling scheme that enables substantial computational gains through multilevel Monte Carlo and multifidelity approaches.

Original authors: Thomas P. Steele, David J. Warne

Published 2026-08-20

📖 6 min read🧠 Deep dive

CC BY 4.0

Original authors: Thomas P. Steele, David J. Warne

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ⚕️ This is an AI-generated explanation of a preprint that has not been peer-reviewed. It is not medical advice. Do not make health decisions based on this content. Read full disclaimer

Life inside a cell is a chaotic, noisy place. Molecules collide, react, and transform in a constant storm of randomness. For decades, scientists have used mathematical models to make sense of this chaos, treating these molecular interactions like a series of instantaneous snaps. In these traditional models, once a reaction starts, it finishes immediately, and the system moves forward based only on where it is right now, forgetting everything that happened before. This approach, known as a Markovian process, has been incredibly useful for understanding many biological systems. However, it fails to capture the reality of complex processes like gene transcription, where the journey from a DNA instruction to a finished protein takes a measurable amount of time. In these cases, the history of the reaction matters; the system remembers how long it has been working on a task, and that memory changes how likely the task is to finish.

When scientists try to model these "non-Markovian" systems, where time and history play a crucial role, the mathematics becomes incredibly difficult. The standard tools for simulating these networks become too slow to be useful, and figuring out the hidden rules of the system from experimental data becomes nearly impossible. This creates a bottleneck for researchers trying to understand everything from how cells differentiate to how diseases spread. A team of researchers at the Queensland University of Technology has now developed a new set of tools to break through this barrier. They created faster ways to simulate these complex, memory-dependent systems and a clever method to guess the hidden parameters of these models much more efficiently than before.

The core of the problem lies in how biological reactions are timed. In a simple model, a reaction happens the moment the conditions are right. In the real world, a reaction like the transcription of a gene into messenger RNA involves a series of steps that take time. The researchers realized that to model this accurately, they had to treat the reaction as a process with a start and a finish, where the time between them is not fixed but depends on the current state of the cell and how long the process has already been running. To handle this, they adapted a classic simulation technique called the "next reaction method." This method usually tracks when the next chemical event will happen in a system. The researchers modified it to keep track of every ongoing reaction, remembering when each one started and how long it has been running. This allows the computer to calculate exactly when a delayed reaction will complete, even if the rules for that completion change as the cell changes.

While this exact method is accurate, it is still very slow because it has to calculate every single step of every reaction individually. To speed things up, the team developed a second, approximate method. This approach, an extension of a technique known as "tau-leaping," looks at the system over small chunks of time rather than one instant at a time. Instead of tracking every single molecule collision, it estimates how many reactions will happen in that short window. This approximation is much faster but introduces a small amount of error. The breakthrough came when the researchers figured out how to link these two methods together. They created a "coupling scheme," a way to run the slow, exact simulation and the fast, approximate simulation side-by-side using the same random events. Because they are running on the same random path, the errors in the fast simulation are highly correlated with the exact one. This correlation is the key to efficiency.

By using this connection, the researchers could apply a powerful statistical trick called "multilevel Monte Carlo." Instead of running millions of expensive, exact simulations to get a precise answer, they run a few exact simulations and many more cheap, approximate ones. The fast simulations provide the bulk of the data, while the slow ones provide a correction that removes the bias. This allows them to get the accuracy of the slow method with a fraction of the computational cost. They tested this approach on two very different biological scenarios. The first was a gene regulation network where a protein inhibits its own production, a process that naturally creates oscillations or rhythms. The second was a model of disease spread where the time it takes for an infected person to recover is not fixed but follows a specific distribution. In both cases, their new methods produced results that matched the exact, slow simulations but did so much faster.

The true power of this work shines when trying to infer the hidden rules of a system from real data. In biology, scientists often have measurements of how many molecules are present at certain times, but they do not know the exact rates at which reactions occur. To find these rates, they usually have to run thousands of simulations to see which set of rules produces data that looks like the real world. This is computationally exhausting. The researchers applied their new coupled simulation method to this problem, using a technique called "multifidelity approximate Bayesian computation." They compared their method against the standard, slow approach. The results were striking. Their new method achieved the same level of accuracy in estimating the hidden parameters as the standard method, but it did so with about five times less computing time. In some comparisons, the speed-up was even more dramatic, reaching an order of magnitude.

The researchers demonstrated that their approach works not just for simple models but for complex systems where the reaction times depend on the state of the system and the passage of time. They showed that by carefully managing the relationship between the fast and slow simulations, they could drastically reduce the time needed to understand these systems. This opens the door for applying these detailed, realistic models to larger and more complex problems, such as modeling entire cells or the spread of epidemics with more realistic recovery times. The work does not claim to have solved every problem in biological modeling, but it provides a robust and efficient foundation. It proves that by combining exact and approximate methods in a smart way, scientists can now simulate and understand the history-dependent, non-Markovian processes that are fundamental to life, without being held back by the limits of computing power.

Technical Summary: Simulation and Inference Methods for Non-Markovian Stochastic Reaction Networks

Problem Statement
Stochastic models of reaction networks are essential for capturing intrinsic noise in complex biological systems, such as gene transcription and translation. While standard formulations rely on Markov processes (where future dynamics depend only on the current state), many biological processes exhibit history-dependent dynamics that necessitate non-Markovian modeling. These processes often involve stochastic delays where the duration between reaction initiation and completion depends on both system state and time.

The primary challenge lies in the computational intractability of simulating and performing parameter inference on these non-Markovian systems.

Simulation: Exact simulation methods for non-Markovian systems (e.g., the non-Markovian Gillespie Algorithm) are computationally expensive, particularly when delays depend on system state. Furthermore, existing literature lacks efficient approximate simulation schemes (analogous to $\tau$ -leaping) that can be coupled with exact methods.
Inference: Bayesian inference for these systems is hindered by the intractability of the likelihood function, which depends on the Chemical Master Equation (CME). Consequently, researchers rely on likelihood-free, simulation-based inference methods (e.g., Approximate Bayesian Computation, ABC). However, these methods are computationally intensive, requiring vast numbers of simulations.
Multifidelity Gap: Acceleration strategies like Multilevel Monte Carlo (MLMC) and multifidelity schemes, which combine cheap low-fidelity approximations with expensive high-fidelity simulations, have been successfully applied to Markovian systems. However, their application to non-Markovian systems has been limited due to the lack of accurate approximate simulators and, crucially, the absence of exact coupling schemes required to generate positively correlated sample paths between exact and approximate models.

Methodology
The authors develop a comprehensive framework for simulating and inferring parameters in well-mixed non-Markovian stochastic reaction networks with state- and time-dependent delays.

Mathematical Formulation:
The authors extend the Kurtz random time-change representation to accommodate delay reactions. They define a system state vector that includes not only chemical species counts but also counts of ongoing delayed reactions ( $X^*_j(t)$ ). This allows the system to track the internal time of each initiated delay, enabling the modeling of non-exponential waiting time distributions.
Exact Stochastic Simulation (nM-NRM):
The authors generalize the Next Reaction Method (NRM) to the non-Markovian setting (Algorithm 3).
- Instead of sampling a global waiting time, the method maintains independent unit-rate Poisson processes for each reaction channel (both instantaneous and delay initiations) and for each active delay completion.
- It utilizes cumulative hazard functions that depend on the current state and the internal time of ongoing delays.
- The algorithm dynamically updates "remaining intervals" for all active channels following every reaction event, allowing for exact simulation of arbitrary inter-event time distributions dependent on state and time.
Approximate Stochastic Simulation (nM-TLM):
To enable acceleration, the authors derive a non-Markovian $\tau$ -leaping method (Algorithm 4).
- This method assumes the system state remains constant over a time step $\tau$ .
- It approximates the cumulative hazard functions by integrating propensities over the fixed interval $\tau$ .
- Crucially, it handles delay initiations in "groups" (bulk arrivals) at the end of the time step, treating multiple initiations within $\tau$ as having the same initiation time. This discretization allows for the generation of approximate sample paths that are computationally cheaper than the exact method.
Coupling Scheme:
The authors introduce a novel coupling mechanism (detailed in Appendix B) to generate positively correlated pairs of exact (nM-NRM) and approximate (nM-TLM) sample paths.
- Both simulations are driven by the same underlying inhomogeneous Poisson processes (common random numbers).
- This coupling ensures that the difference between the exact and approximate paths is minimized, which is a prerequisite for the variance reduction required by MLMC and multifidelity inference.
Inference Framework:
The authors apply these simulation tools to Multifidelity Approximate Bayesian Computation (MF-ABC).
- They use the nM-TLM as a low-fidelity simulator and the nM-NRM as a high-fidelity simulator.
- Using the coupling scheme, they implement a weighting strategy (Equation 11) that corrects the bias introduced by the low-fidelity approximation, yielding an unbiased estimator of the posterior distribution.

Key Contributions

Generalized Exact Simulation: A non-Markovian extension of the Next Reaction Method (nM-NRM) capable of handling delays that depend on both time and system state, generalizing previous work restricted to fixed or time-only dependent delays.
Non-Markovian $\tau$ -Leaping: The first approximate stochastic simulation algorithm (nM-TLM) for non-Markovian networks that supports arbitrary inter-event time distributions, providing a computationally efficient low-fidelity alternative.
Exact Coupling Scheme: A novel method to generate positively correlated exact and approximate sample paths for non-Markovian systems. This fills a critical gap, enabling the application of MLMC and multifidelity methods to non-Markovian settings for the first time.
Multifidelity Inference: Demonstration that these simulation tools can be integrated into MF-ABC to significantly accelerate parameter inference for non-Markovian models.

Results
The authors validate their methods using two case studies:

Delayed Auto-Inhibition Gene Regulation: A model of mRNA transcription and protein translation with state-dependent delays.
- Convergence: The nM-NRM correctly reduces to the standard Markovian NRM when the delay distribution becomes exponential ( $\alpha \to 1$ ).
- Approximation Error: The nM-TLM exhibits expected weak ( $O(\tau)$ ) and strong ( $O(\tau^{1/2})$ ) convergence rates. While the approximation introduces path-wise errors (e.g., overestimating peak protein levels with large $\tau$ ), the coupling scheme effectively reduces variance when comparing exact and approximate paths.
Non-Markovian SIS Epidemic Model: A model where infection and recovery processes follow Weibull waiting time distributions.
- MLMC Efficiency: Using the unbiased MLMC estimator, the authors achieved order-of-magnitude improvements in computational efficiency (MSE vs. cost) compared to direct Monte Carlo sampling. Smaller time steps $\tau$ were found to be more efficient in this setting due to stronger coupling effects.
Inference Performance:
- In the gene regulation model, the MF-ABC scheme (using nM-TLM for low-fidelity and nM-NRM for correction) achieved a 5-fold speed-up in computational cost compared to standard ABC rejection sampling for equivalent posterior mean accuracy.
- The MF-ABC method successfully recovered posterior marginal probability densities for unknown parameters ( $\beta_{m*}, \beta_p, \alpha$ ) with a 10-fold computational speed-up compared to a gold-standard ABC estimate.

Significance and Claims
The paper claims to provide a foundational framework that removes a critical barrier to the practical implementation of non-Markovian models in systems biology and beyond. By extending the Next Reaction Method and $\tau$ -leaping to non-Markovian settings and, most importantly, establishing a valid coupling scheme, the authors enable the use of powerful variance reduction techniques (MLMC and multifidelity methods) for these complex systems.

The authors emphasize that their contributions are not limited to the specific ABC implementation demonstrated. The simulation algorithms and coupling schemes are general and applicable to any simulation-based inference method (including Sequential Monte Carlo, neural posterior estimation, or other MLMC variants) that relies on coupling exact and approximate simulations. The work suggests that these methods could enable the simulation of whole-cell models and other computationally expensive quantitative biology problems that were previously intractable due to the cost of non-Markovian stochastic simulation.

Technical Summary: Simulation and Inference Methods for Non-Markovian Stochastic Reaction Networks

More like this