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

Single-cell data reveal the presence of biological stochasticity between cells of identical genome and environment, in particular highlighting the transcriptional bursting phenomenon. To account for this property, gene expression may be modeled as a continuous-time Markov chain where biochemical species are described in a discrete way, leading to Gillespie's stochastic simulation algorithm (SSA) which turns out to be computationally expensive for realistic mRNA and protein copy numbers. Alternatively, hybrid models based on piecewise-deterministic Markov processes (PDMPs) offer an effective compromise for capturing cell-to-cell variability, but their simulation remains limited to specialized mathematical communities. With a view to making them more accessible, we present here a simple simulation method that is reminiscent of SSA, while allowing for much lower computational cost. We detail the algorithm for a bursty PDMP describing an arbitrary number of interacting genes, and prove that it simulates exact trajectories of the model. As an illustration, we use the algorithm to simulate a two-gene toggle switch: this example highlights the fact that bimodal distributions as observed in real data are not explained by transcriptional bursting per se, but rather by distinct burst frequencies that may emerge from interactions between genes.