Spectral solutions to stochastic models of gene expression with bursts and regulation

Signal-processing molecules inside cells are often present at low copy number, which necessitates probabilistic models to account for intrinsic noise. Probability distributions have traditionally been found using simulation-based approaches which then require estimating the distributions from many samples. Here we present in detail an alternative method for directly calculating a probability distribution by expanding in the natural eigenfunctions of the governing equation, which is linear. We apply the resulting spectral method to three general models of stochastic gene expression: a single gene with multiple expression states (often used as a model of bursting in the limit of two states), a gene regulatory cascade, and a combined model of bursting and regulation. In all cases we find either analytic results or numerical prescriptions that greatly outperform simulations in efficiency and accuracy. In the last case, we show that bimodal response in the limit of slow switching is not only possible but optimal in terms of information transmission.

💡 Research Summary

The paper addresses the challenge of quantifying intrinsic noise in gene expression when molecular copy numbers are low. Traditional approaches rely on stochastic simulations (e.g., Gillespie algorithm) followed by histogram estimation, which become computationally expensive and statistically noisy for large parameter sweeps or rare events. The authors propose a fundamentally different strategy: they solve the master equation analytically by expanding the probability distribution in the natural eigenfunctions of the linear operator governing the system. Because the master equation is linear, its eigenfunctions form a complete basis—typically a mixture of Poisson and geometric forms—each associated with an eigenvalue that reflects a characteristic relaxation time scale. By projecting the initial condition onto this basis, the time evolution reduces to simple exponential decay or growth of the expansion coefficients, yielding an exact or highly accurate representation of the full probability distribution without any Monte‑Carlo sampling.

The method is applied to three increasingly complex models.

1. Bursting gene model (multiple expression states).
   The gene switches among several promoter states, each with its own transcription rate. In the two‑state limit the model reproduces transcriptional bursts. Using the spectral expansion the authors obtain a closed‑form probability mass function for mRNA copy number. Computationally, the spectral solution is orders of magnitude faster than stochastic simulation while preserving exact normalization and allowing straightforward sensitivity analysis with respect to switching rates and transcription rates.

2. Gene‑regulatory cascade.
   Here a protein product of an upstream gene modulates the transcription rate of a downstream gene, creating a chain of stochastic reactions (transcription, translation, degradation). The authors treat each layer as a linear master equation, solve it in its own eigenbasis, and then combine the solutions via tensor products to obtain the joint distribution of all molecular species. This approach remains numerically stable even when feedback loops are introduced, and it provides direct access to higher‑order moments and cross‑correlations that would be costly to estimate from simulation data.

3. Combined bursting and regulation model.
   This most general case couples slow promoter switching (bursting) with a nonlinear feedback where the downstream protein concentration feeds back on the upstream transcription rate. The spectral expansion reveals two dominant modes corresponding to the “on” and “off” promoter states; when switching is sufficiently slow the distribution becomes bimodal. The authors go further by evaluating information‑theoretic metrics—Fisher information and Shannon channel capacity—showing that the bimodal regime maximizes the transmission of information about an external signal (e.g., a stimulus that modulates the switching propensity). They identify a specific region in the parameter space (slow but not infinitely slow switching, strong nonlinear regulation) where the channel capacity peaks, demonstrating that cellular systems can exploit stochastic bursting to achieve optimal information flow.

Beyond the specific models, the paper discusses the broader applicability of the spectral method. While strictly limited to linear master equations, many biologically relevant nonlinearities can be linearized around steady states or treated perturbatively, allowing the eigenfunction expansion to serve as a powerful approximation tool. Moreover, the authors suggest integrating Bayesian inference with the spectral coefficients to fit experimental single‑cell data directly, turning the analytical framework into a practical parameter‑estimation pipeline.

In summary, the authors provide a rigorous, efficient, and versatile analytical technique for solving stochastic gene‑expression models. The spectral method outperforms simulation in speed and accuracy, yields exact probability distributions, and offers deep insight into how bursting and regulatory feedback shape cellular noise and information transmission. This work bridges theoretical biophysics and quantitative systems biology, offering a valuable toolbox for researchers designing synthetic circuits, interpreting single‑cell transcriptomics, or exploring the fundamental limits of cellular communication.