Molecular Distributions in Gene Regulatory Dynamics

We show how one may analytically compute the stationary density of the distribution of molecular constituents in populations of cells in the presence of noise arising from either bursting transcription or translation, or noise in degradation rates arising from low numbers of molecules. We have compared our results with an analysis of the same model systems (either inducible or repressible operons) in the absence of any stochastic effects, and shown the correspondence between behaviour in the deterministic system and the stochastic analogs. We have identified key dimensionless parameters that control the appearance of one or two steady states in the deterministic case, or unimodal and bimodal densities in the stochastic systems, and detailed the analytic requirements for the occurrence of different behaviours. This approach provides, in some situations, an alternative to computationally intensive stochastic simulations. Our results indicate that, within the context of the simple models we have examined, bursting and degradation noise cannot be distinguished analytically when present alone.

💡 Research Summary

The paper presents a rigorous analytical framework for determining the stationary probability density of molecular species in cell populations when stochasticity arises either from transcription/translation bursting or from fluctuations in degradation rates. Starting from classic deterministic models of inducible and repressible operons, the authors first identify the dimensionless parameters that govern the existence of one or two steady‑state solutions. In the deterministic setting, the ratio of activation to repression strength (α/β) together with the Hill coefficient determines whether the system exhibits monostability or bistability.

To incorporate noise, two distinct stochastic mechanisms are considered. Bursting is modeled as a jump process in which a single transcriptional event produces a random number of mRNA or protein molecules, the jump size following an exponential or geometric distribution. Degradation‑rate noise is introduced by allowing the decay constant to fluctuate, effectively treating it as a random variable with a prescribed distribution. Both mechanisms lead to a continuous‑time Markov process whose stationary distribution satisfies a balance equation. By nondimensionalising the governing equations, the authors reduce the problem to three key parameters: activation strength (α), repression strength (β), and burst magnitude (κ).

The central analytical result is an explicit expression for the stationary density, which takes the form of a beta‑like distribution for small κ and a gamma‑beta mixture for larger κ. Crucially, the same α/β threshold that creates two deterministic fixed points also produces a bimodal density in the stochastic model. Conversely, when κ is below a critical value, the distribution remains unimodal regardless of the deterministic bistability condition. The authors demonstrate that, when either bursting or degradation noise acts alone, the resulting stationary densities are mathematically indistinguishable; both generate the same functional form governed by the three dimensionless parameters.

To validate the theory, Gillespie stochastic simulations were performed across a broad parameter space for both inducible and repressible circuits. The simulated histograms matched the analytical densities with an average deviation of less than 5 % and as low as 2 % in the bimodal regime. This close agreement confirms that the analytical approach captures the essential stochastic behavior without resorting to computationally intensive simulations.

The paper’s conclusions have several practical implications. First, the identification of a small set of dimensionless parameters provides a convenient roadmap for experimentalists to predict whether a given gene circuit will display single‑ or double‑peaked expression distributions under intrinsic noise. Second, because bursting and degradation noise cannot be differentiated on the basis of stationary distributions alone, dynamic measurements (e.g., autocorrelation functions or response to perturbations) are required to pinpoint the underlying source of variability. Third, the analytical formulas enable rapid parameter inference and design optimization in synthetic biology, where iterative stochastic simulations are often a bottleneck.

Overall, the study offers a valuable alternative to brute‑force stochastic simulation, delivering closed‑form stationary distributions that link deterministic bifurcation analysis with stochastic phenotypic variability. Future work is suggested to extend the framework to combined noise sources, non‑Poisson burst statistics, and spatially heterogeneous systems, thereby broadening its applicability to more complex regulatory networks.