🧬 biology

Efficient stochastic simulation of gene regulatory networks using hybrid models of transcriptional bursting

This paper introduces a computationally efficient, SSA-like simulation algorithm for piecewise-deterministic Markov processes that accurately models transcriptional bursting in gene regulatory networks, demonstrating through a toggle switch example that bimodal distributions arise from distinct burst frequencies driven by gene interactions rather than bursting alone.

Original authors: Mathilde Gaillard, Ulysse Herbach

Published 2026-05-19

📖 4 min read☕ Coffee break read

CC BY 4.0

Original authors: Mathilde Gaillard, Ulysse Herbach

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

Imagine a cell as a busy factory. Inside this factory, there are machines (genes) that produce blueprints (mRNA) and then build products (proteins). For a long time, scientists thought these machines worked like a smooth, steady assembly line: a little bit of blueprint comes out, then a little bit more, constantly.

But new, high-tech cameras (single-cell data) revealed something surprising: the machines don't work smoothly. Instead, they work in bursts. Imagine a machine that sits silent for a long time, then suddenly goes into overdrive, churning out a huge pile of blueprints in a few seconds, before shutting down again. This is called "transcriptional bursting."

The Problem: The "Slow Motion" Camera

To understand how these factories work, scientists use computer simulations. The traditional way to do this is like watching a movie frame-by-frame, counting every single blueprint and protein molecule as it appears or disappears.

The Issue: In a real cell, there are thousands of these molecules. Counting every single one, one by one, is like trying to count every grain of sand on a beach while the tide is coming in. It takes a computer an incredibly long time to run these simulations, making it too slow for studying complex networks of many genes.

The Solution: The "Hybrid" Shortcut

The authors of this paper, Mathilde Gaillard and Ulysse Herbach, propose a smarter way to simulate these factories. They call it a hybrid model.

Think of it this way:

The Old Way (SSA): You count every single grain of sand (molecule) individually.
The New Way (PDMP): You treat the sand as a flowing river (continuous) but acknowledge that the water is being dumped in by buckets (discrete bursts).

In their model:

The Flow: The levels of mRNA and protein that are already inside the cell are treated like a smooth, flowing river. They rise and fall predictably, like water draining from a bathtub or filling up a cup.
The Buckets: The "bursts" are treated as sudden, random events where a bucket of new material is dumped into the river.

The Algorithm: The "Ghost Jump" Trick

The paper introduces a specific computer recipe (an algorithm) to simulate this hybrid system. Here is how it works, using a simple analogy:

Imagine you are walking down a hallway, and you want to know where you will be in 10 minutes.

The Clock: You have a timer that rings at random intervals (like a Poisson process).
The Walk: Between rings, you walk in a straight line (the smooth flow of the river).
The Decision: When the timer rings, you have to decide: "Did I actually get pushed by a bucket, or was it just a false alarm?"
- The computer picks a random number. If it's a "false alarm" (a phantom jump), you keep walking straight. Nothing changes.
- If it's a "real event" (a true jump), you stop, and a bucket of mRNA is dumped into your system, changing your path.

Why is this cool?

Speed: Because the computer doesn't have to count every single molecule, it runs about 200 times faster than the old method.
Accuracy: The authors mathematically proved that even though they are using this "shortcut" with phantom jumps, the final result is exactly the same as if they had counted every single molecule. It's like taking a shortcut through a park that leads to the exact same destination as the long road, without getting lost.

The Big Discovery: Why Cells Look Different

The authors used their fast, new simulation to study a "toggle switch" (a system with two genes that turn each other off). They found something interesting about why cells look so different from one another (variability).

Old Idea: People thought the "on/off" switching of the gene itself caused the big differences.
New Insight: The simulation showed that the big differences (bimodal distributions) actually come from the frequency of the bursts. Sometimes a gene bursts very rarely; other times, it bursts often. The proteins in the cell act like a buffer, smoothing out these bursts, but the different "bursting rhythms" create two distinct groups of cells.

Summary

This paper gives scientists a fast, accurate, and easy-to-use tool to simulate how genes behave in single cells. Instead of getting bogged down counting every single molecule, they can now watch the "river flow" and the "bucket drops." This allows them to study complex genetic networks much faster while still capturing the real, messy, bursty nature of life inside a cell.

The authors have even made the code for this tool available for free, so other scientists can use it immediately.

Technical Summary: Efficient Stochastic Simulation of Gene Regulatory Networks Using Hybrid Models of Transcriptional Bursting

Problem Statement
Single-cell data has revealed significant biological stochasticity between cells with identical genomes and environments, most notably the phenomenon of transcriptional bursting. To model this, gene expression is often treated as a continuous-time Markov chain with discrete biochemical species, typically simulated using Gillespie's Stochastic Simulation Algorithm (SSA). However, SSA becomes computationally prohibitive for realistic mRNA and protein copy numbers. While Piecewise-Deterministic Markov Processes (PDMPs) offer a hybrid framework that captures cell-to-cell variability by treating abundant species (mRNA, protein) as continuous variables and rare events (transcriptional bursts) as discrete jumps, their simulation has remained limited to specialized mathematical communities. Existing PDMP simulation methods often rely on numerical integration of survival functions, which can be computationally expensive and require time discretization.

Methodology
The authors propose a simple, exact simulation algorithm for bursty PDMPs describing interacting gene networks. The model assumes a "bursty regime" where promoter active periods are short compared to inactive periods, allowing the promoter state to be discarded. In this regime, mRNA is produced via instantaneous bursts with a frequency $kon,i$ that depends on protein levels, while mRNA and protein degradation and translation follow deterministic ordinary differential equations (ODEs).

The core of the methodology is an acceptance-rejection algorithm reminiscent of SSA but adapted for continuous state spaces:

Global Clock: A global upper bound $\lambda$ is established for the sum of all burst frequencies ( $\sum kon,i$ ).
Deterministic Flow: Between jumps, the system evolves deterministically according to the analytical solution of the ODEs governing mRNA and protein dynamics.
Event Simulation: The algorithm simulates a homogeneous Poisson process with intensity $\lambda$ to generate candidate jump times.
Acceptance-Rejection: At each candidate time, the algorithm decides whether a "true" jump occurs or if it is a "phantom" jump (rejection). This decision is based on the ratio of the actual burst frequency at the current state to the global bound $\lambda$ .
State Update: If a true jump occurs, a specific gene is selected based on the relative burst frequencies, and an exponentially distributed burst size is added to that gene's mRNA level. The process then restarts from this new state.

Key Contributions

Algorithmic Efficiency: The paper presents a simulation method that avoids the numerical integration required by standard PDMP simulation techniques. By using an acceptance-rejection step with a global bound, it bypasses the need to invert cumulative distribution functions or discretize time, significantly reducing computational cost.
Exactness Proof: The authors provide a rigorous mathematical proof (Section 5) demonstrating that the constructed process is an exact realization of the target PDMP. They show that the infinitesimal generator of the simulated process matches the theoretical generator of the bursty PDMP model, ensuring that the trajectories are statistically correct.
Modularity and Accessibility: The algorithm is designed to be easily implemented and extended. It accommodates arbitrary numbers of interacting genes and can be adapted to include additional variables (e.g., distinguishing mature mRNA from pre-mRNA) without altering the fundamental computational structure.
Biological Insight via Simulation: Using the algorithm to simulate a two-gene toggle switch, the authors demonstrate that bimodal distributions observed in real data are not solely a result of transcriptional bursting itself. Instead, bimodality emerges from distinct burst frequencies driven by gene-gene interactions. The study highlights that the bursty PDMP model correctly captures distribution tails observed in single-cell data, whereas standard two-state promoter models may fail to do so unless specific, often unrealistic, parameter regimes are met.

Results

Computational Performance: The authors report that for realistic parameter values, the PDMP simulation is approximately 200 times faster than the exact SSA, making it feasible to simulate large networks and collect extensive single-cell data trajectories.
Model Validation: The simulation of a toggle switch network successfully reproduced a bistable pattern. The results indicated that proteins, due to slower degradation rates, buffer mRNA bursts, creating robust bimodal distributions.
Theoretical Consistency: The proof confirms that the algorithm generates trajectories with the correct infinitesimal generator, validating the method as a rigorous alternative to discrete SSA for systems where species counts are high enough to justify continuous approximation but low enough to require stochastic treatment of bursts.

Significance and Claims
The paper argues that transcriptional bursting represents an ideal scenario for hybrid modeling, where abundant chemical species are described via ODEs while discrete events related to rare species drive the stochasticity. The authors claim that their work makes the simulation of bursty PDMPs accessible to a broader systems biology community, offering a "simple" yet "exact" alternative to both inefficient SSA and complex, specialized PDMP solvers.

The authors emphasize that their contribution is not a novel theoretical framework (as the algorithm is a specific case of general PDMP simulation methods) but rather a dedicated, practical instance for gene regulatory networks. They explicitly state that the primary motivation is to facilitate the efficient simulation of large networks and to provide a detailed proof of the algorithm's exactness, which they found lacking in the existing literature. The paper concludes that this approach is particularly well-suited for interpreting single-cell transcriptomics data, where capturing the correct distribution tails is crucial.