An Algorithm for the Stochastic Simulation of Gene Expression and Heterogeneous Population Dynamics

We present an algorithm for the stochastic simulation of gene expression and heterogeneous population dynamics. The algorithm combines an exact method to simulate molecular-level fluctuations in single cells and a constant-number Monte Carlo method to simulate time-dependent statistical characteristics of growing cell populations. To benchmark performance, we compare simulation results with steadystate and time-dependent analytical solutions for several scenarios, including steadystate and time-dependent gene expression, and the effects on population heterogeneity of cell growth, division, and DNA replication. This comparison demonstrates that the algorithm provides an efficient and accurate approach to simulate how complex biological features influence gene expression. We also use the algorithm to model gene expression dynamics within “bet-hedging” cell populations during their adaption to environmental stress. These simulations indicate that the algorithm provides a framework suitable for simulating and analyzing realistic models of heterogeneous population dynamics combining molecular-level stochastic reaction kinetics, relevant physiological details and phenotypic variability.

💡 Research Summary

The paper introduces a hybrid computational framework that merges the exact Gillespie stochastic simulation algorithm (SSA) for intracellular reaction dynamics with a constant‑number Monte Carlo (CN‑MC) scheme for tracking population‑level statistics in growing cell cultures. The SSA component treats each biochemical reaction as a discrete stochastic event, computing propensities based on current molecular counts and adjusting higher‑order reaction rates by the instantaneous cell volume V(t). Cell volume follows an exponential growth law V₀·exp(ln 2·t_div/τ₀), and division is triggered when the volume reaches a predefined threshold (typically twice the birth volume). Upon division, molecular content is partitioned probabilistically, allowing both symmetric and asymmetric inheritance. Additional discrete events such as DNA replication, cell‑cycle checkpoints, and time‑delayed reactions are managed through a per‑cell event stack, enabling flexible modeling of gene‑copy number changes and transcriptional bursts.

The CN‑MC component solves the population balance problem without explicitly simulating every descendant. A fixed number N of “sample cells” represents the whole culture. When a cell divides, the newly created daughter is stored in a separate array; at regular intervals (t_restore) the oldest daughters replace randomly chosen mothers. This replacement strategy preserves the age (or volume) distribution while keeping the simulated population size constant, thus avoiding the exponential blow‑up of computational resources.

Implementation details include a shared‑memory parallelization where each cell’s SSA is executed independently on multiple cores. The main simulation loop proceeds as follows: (1) parallel SSA updates for all cells until the next sampling time, (2) volume updates and execution of any scheduled events, (3) division handling and daughter bookkeeping, (4) after the sampling interval, invoke the CN‑MC replacement routine, and (5) compute population statistics (means, variances, histograms). Benchmarks on a dual‑quad‑core IBM machine show that a 10⁵‑second simulation of 8 000 cells completes in roughly 81 seconds, demonstrating orders‑of‑magnitude speed‑up over naïve full‑population simulations.

To validate accuracy, the authors compare simulation outputs against analytical steady‑state and time‑dependent solutions for progressively more complex models: (i) a simple birth‑death gene expression system, (ii) a model incorporating cell growth, division, and DNA replication, and (iii) a “bet‑hedging” scenario where environmental stress induces a rapid shift in transcription rates. In all cases, simulated means, variances, and full distributions match theoretical predictions, confirming that the hybrid algorithm faithfully captures both intrinsic molecular noise and extrinsic population heterogeneity.

The study concludes that the SSA + CN‑MC combination provides an efficient, scalable, and biologically realistic tool for exploring stochastic gene expression in heterogeneous populations. Its modular design permits inclusion of additional physiological details (e.g., nutrient limitation, spatial effects) and can be extended to larger systems such as microbial consortia, tissue development, or drug‑response studies, where both single‑cell variability and population dynamics are critical.