Stochastic fluctuations of molecule numbers are ubiquitous in biological systems. Important examples include gene expression and enzymatic processes in living cells. Such systems are typically modelled as chemical reaction networks whose dynamics are governed by the chemical master equation. Despite its simple structure, no analytic solutions to the chemical master equation are known for most systems. Moreover, stochastic simulations are computationally expensive, making systematic analysis and statistical inference a challenging task. Consequently, significant effort has been spent in recent decades on the development of efficient approximation and inference methods. This article gives an introduction to basic modelling concepts as well as an overview of state of the art methods. First, we motivate and introduce deterministic and stochastic methods for modelling chemical networks, and give an overview of simulation and exact solution methods. Next, we discuss several approximation methods, including the chemical Langevin equation, the system size expansion, moment closure approximations, time-scale separation approximations and hybrid methods. We discuss their various properties and review recent advances and remaining challenges for these methods. We present a comparison of several of these methods by means of a numerical case study and highlight some of their respective advantages and disadvantages. Finally, we discuss the problem of inference from experimental data in the Bayesian framework and review recent methods developed the literature. In summary, this review gives a self-contained introduction to modelling, approximations and inference methods for stochastic chemical kinetics.

Abstract The stochasticity of gene expression is manifested in the fluctuations of mRNA and protein copy numbers within a cell lineage over time. While data of this type can be obtained for many generations, most mathematical models are unsuitable to interpret such data since they assume non-growing cells. Here we develop a theoretical approach that quantitatively links the frequency content of lineage data to subcellular dynamics. We elucidate how the position, height, and width of the peaks in the power spectrum provide a distinctive fingerprint that encodes a wealth of information about mechanisms controlling transcription, translation, replication, degradation, bursting, promoter switching, cell cycle duration, cell division, and gene dosage compensation. Predictions are confirmed by analysis of single-cell Escherichia coli data obtained using fluorescence microscopy. Furthermore, by matching the experimental and theoretical power spectra, we infer the temperature-dependent gene expression parameters, without the need of measurements relating fluorescence intensities to molecule numbers.

Abstract Intracellular reaction rates depend on concentrations and hence their levels are often regulated. However classical models of stochastic gene expression lack a cell size description and cannot be used to predict noise in concentrations. Here, we construct a model of gene product dynamics that includes a description of cell growth, cell division, size-dependent gene expression, gene dosage compensation, and size control mechanisms that can vary with the cell cycle phase. We obtain expressions for the approximate distributions and power spectra of concentration fluctuations which lead to insight into the emergence of concentration homeostasis. Furthermore, we find that (i) the conditions necessary to suppress cell division-induced concentration oscillations are difficult to achieve; (ii) mRNA concentration and number distributions can have different number of modes; (iii) certain size control strategies are ideal because they maintain constant mean concentrations whilst minimising concentration noise. Predictions are confirmed using lineage data for E. coli , fission yeast and budding yeast.

Abstract In the noisy cellular environment, RNAs and proteins are subject to considerable stochastic fluctuations in copy numbers over time. As a consequence, single cells within the same isoclonal population can differ in their expression profile and reside in different phenotypic states. The dynamic nature of this intercellular variation, where individual cells can transition between different states over time makes it a particularly hard phenomenon to characterize. Here we propose a novel fluctuation-test approach to infer the kinetics of transitions between cell states. More specifically, single cells are randomly drawn from the population and grown into cell colonies. After growth for a fixed number of generations, the number of cells residing in different states is assayed for each colony. In a simple system with reversible switching between two cell states, our analysis shows that the extent of colony-to-colony fluctuations in the fraction of cells in a given state is monotonically related to the switching kinetics. Several closed-form formulas for inferring the switching rates from experimentally quantified fluctuations are presented. We further extend this approach to multiple cell states where harnessing fluctuation signatures can reveal both the topology and the rates of cell-state switching. In summary, our analysis provides a powerful approach for dissecting cell-state transitions based on a single time point measurement. This is especially important for scenarios where a measurement involves killing the cell (for example, performing single-cell RNA-seq or assaying whether a microbial/cancer cell is in a drug-sensitive or drug-tolerant state), and hence the state of the same cell cannot be measured at different time points.

Abstract Transcriptional bursting is a major source of noise in gene expression. The telegraph model of gene expression, whereby transcription switches between “on” and “off” states, is the dominant model for bursting. Recently it was shown that the telegraph model cannot explain a number of experimental observations from perturbation data. Here we study an alternative model that is consistent with the data and which explicitly describes RNA polymerase recruitment and polymerase pause release, two steps necessary for mRNA production. We derive the exact steady-state distribution of mRNA numbers and an approximate steady-state distribution of protein numbers which are given by generalized hypergeometric functions. The theory is used to calculate the relative sensitivity of the coefficient of variation of mRNA fluctuations for thousands of genes in mouse fibroblasts. This indicates that the size of fluctuations is mostly sensitive to the rate of burst initiation and the mRNA degradation rate. Furthermore we show that (i) the time-dependent distribution of mRNA numbers is accurately approximated by a modified telegraph model with a Michaelis-Menten like dependence of the effective transcription rate on RNA polymerase abundance. (ii) the model predicts that if the polymerase recruitment rate is comparable or less than the pause release rate, then upon gene replication the mean number of RNA per cell remains approximately constant. This gene dosage compensation property has been experimentally observed and cannot be explained by the telegraph model with constant rates. Statement of Significance The random nature of gene expression is well established experimentally. Mathematical modelling provides a means of understanding the factors leading to the observed stochasticity. There is evidence that the classical two-state model of stochastic mRNA dynamics (the telegraph model) cannot describe perturbation experiments and a new model that includes polymerase dynamics has been proposed. In this paper, we present the first detailed study of this model, deriving an exact solution for the mRNA distribution in steady-state conditions, an approximate time-dependent solution and showing the model can explain gene dosage compensation. As well, we use the theory together with transcriptomic data, to deduce which parameters when perturbed lead to a maximal change in the size of mRNA fluctuations.

Abstract Non-Markov models of stochastic biochemical kinetics often incorporate explicit time delays to effectively model large numbers of intermediate biochemical processes. Analysis and simulation of these models, as well as the inference of their parameters from data, are fraught with difficulties because the dynamics depends on the system’s history. Here we use an artificial neural network to approximate the time-dependent distributions of non-Markov models by the solutions of much simpler time-inhomogeneous Markov models; the approximation does not increase the dimensionality of the model and simultaneously leads to inference of the kinetic parameters. The training of the neural network uses a relatively small set of noisy measurements generated by experimental data or stochastic simulations of the non-Markov model. We show using a variety of models, where the delays stem from transcriptional processes and feedback control, that the Markov models learnt by the neural network accurately reflect the stochastic dynamics across parameter space.

A bstract Estimating uncertainty in model predictions is a central task in quantitative biology. Biological models at the single-cell level are intrinsically stochastic and nonlinear, creating formidable challenges for their statistical estimation which inevitably has to rely on approximations that trade accuracy for tractability. Despite intensive interest, a sweet spot in this trade off has not been found yet. We propose a flexible procedure for uncertainty quantification in a wide class of reaction networks describing stochastic gene expression including those with feedback. The method is based on creating a tractable coarse-graining of the model that is learned from simulations, a synthetic model , to approximate the likelihood function. We demonstrate that synthetic models can substantially outperform state-of-the-art approaches on a number of nontrivial systems and datasets, yielding an accurate and computationally viable solution to uncertainty quantification in stochastic models of gene expression.

Abstract Recent advances in fluorescence microscopy have made it possible to measure the fluctuations of nascent (actively transcribed) RNA. These closely reflect transcription kinetics, as opposed to conventional measurements of mature (cellular) RNA, whose kinetics is affected by additional processes downstream of transcription. Here, we formulate a stochastic model which describes promoter switching, initiation, elongation, premature detachment, pausing, and termination while being analytically tractable. By computational binning of the gene into smaller segments, we derive exact closed-form expressions for the mean and variance of nascent RNA fluctuations in each of these segments, as well as for the total nascent RNA on a gene. We also derive exact expressions for the first two moments of mature RNA fluctuations, and approximate distributions for total numbers of nascent and mature RNA. Our results, which are verified by stochastic simulation, uncover the explicit dependence of the statistics of both types of RNA on transcriptional parameters and potentially provide a means to estimate parameter values from experimental data.

Fluctuations in the number of nascent RNA accurately reflect transcriptional activity. However, mathematical models predicting their distributions are difficult to solve analytically due to their non-Markovian nature stemming from transcriptional elongation. Here we circumvent this problem by deriving an exact relationship between the steady-state distribution of nascent RNA and the distribution of initiation times, which can be computed for any general initiation mechanism described by a set of first-order reactions. We test our theory using simulations and live cell imaging data.

Abstract Transcriptional rates are often estimated by fitting the distribution of mature mRNA numbers measured using smFISH (single molecule fluorescence in situ hybridization) with the distribution predicted by the telegraph model of gene expression, which defines two promoter states of activity and inactivity. However, fluctuations in mature mRNA numbers are strongly affected by processes downstream of transcription. In addition, the telegraph model assumes one gene copy, but in experiments cells may have two gene copies as cells replicate their genome during the cell cycle. Whilst it is often presumed that post-transcriptional noise and gene copy number variation affect transcriptional parameter estimation, the size of the error introduced remains unclear. To address this issue, here we measure both mature and nascent mRNA distributions of GAL10 in yeast cells using smFISH and classify each cell according to its cell cycle phase. We infer transcriptional parameters from mature and nascent mRNA distributions, with and without accounting for cell cycle phase and compare the results to live-cell transcription measurements of the same gene. We find that: (i) correcting for cell cycle dynamics decreases the promoter switching rates and the initiation rate, and increases the fraction of time spent in the active state, as well as the burst size; (ii) additional correction for post-transcriptional noise leads to further increases in the burst size and to a large reduction in the errors in parameter estimation. Furthermore, we outline how to correctly adjust for measurement noise in smFISH due to uncertainty in transcription site localisation when introns cannot be labelled. Simulations with parameters estimated from nascent smFISH data, which is corrected for cell cycle phases and measurement noise, leads to autocorrelation functions that agree with those obtained from live-cell imaging.