Search

The complexity of mathematical models in biology has rendered model reduction an essential tool in the quantitative biologist’s toolkit. For stochastic reaction networks described using the Chemical Master Equation, commonly used methods include time-scale separation, Linear Mapping Approximation, and state-space lumping. Despite the success of these techniques, they appear to be rather disparate, and at present, no general-purpose approach to model reduction for stochastic reaction networks is known. In this paper, we show that most common model reduction approaches for the Chemical Master Equation can be seen as minimizing a well-known information-theoretic quantity between the full model and its reduction, the Kullback–Leibler divergence defined on the space of trajectories. This allows us to recast the task of model reduction as a variational problem that can be tackled using standard numerical optimization approaches. In addition, we derive general expressions for propensities of a reduced system that generalize those found using classical methods. We show that the Kullback–Leibler divergence is a useful metric to assess model discrepancy and to compare different model reduction techniques using three examples from the literature: an autoregulatory feedback loop, the Michaelis–Menten enzyme system, and a genetic oscillator.

Stochastic variation in gene products (“noise”) is an inescapable by-product of gene expression. Noise must be minimized to allow for the reliable execution of cellular functions. However, noise cannot be suppressed beyond an intrinsic lower limit. For constitutively expressed genes, this limit is believed to be Poissonian, meaning that the variance in mRNA numbers cannot be lower than their mean. Here, we show that several cell division genes in fission yeast have mRNA variances significantly below this limit, which cannot be explained by the classical gene expression model for low-noise genes. Our analysis reveals that multiple steps in both transcription and mRNA degradation are essential to explain this sub-Poissonian variance. The sub-Poissonian regime differs qualitatively from previously characterized noise regimes, a hallmark being that cytoplasmic noise is reduced when the mRNA export rate increases. Our study re-defines the lower limit of eukaryotic gene expression noise and identifies molecular requirements for ultra-low noise which are expected to support essential cell functions.

The simulation of spatial stochastic models is highly computationally expensive, an issue that has severely limited our understanding of the spatial nature of gene expression. Here we devise a graph neural network based method to learn, from stochastic trajectories in a small region of space, an effective master equation for the time-dependent marginal probability distributions of mRNA and protein numbers at sub-cellular resolution for every cell in a tissue. Numerical solution of this equation leads to accurate results in a small fraction of the computation time of standard simulation methods. Moreover its predictions can be extrapolated to a spatial organisation (a cell network topology) and regions of parameter space unseen in its neural network training. The scalability and accuracy of the method suggest it is a promising approach for whole cell modelling and for detailed comparisons of stochastic models with spatial genomics data.

High-throughput sequencing technology is central to our current understanding of the human methylome. The vast majority of studies use chemical conversion to analyse bulk-level patterns of DNA methylation across the genome from a population of cells. While this technology has been used to probe single-molecule methylation patterns, such analyses are limited to short reads of a few hundred basepairs. DNA methylation can also be directly detected using Nanopore sequencing which can generate reads measuring megabases in length. However, thus far these analyses have largely focused on bulk-level assessment of DNA methylation. Here, we analyse DNA methylation in single Nanopore reads with a mean length of 24.6kb, to show that bulk-level metrics underestimate large-scale heterogeneity in the methylome. We use the correlation in methylation state between neighbouring sites to quantify single-molecule heterogeneity and find that heterogeneity varies significantly across the human genome, with some regions having heterogeneous methylation patterns at the single-molecule level and others possessing more homogeneous methylation patterns. By comparing the genomic distribution of the correlation to epigenomic annotations, we find that the greatest heterogeneity in single-molecule patterns is observed within heterochromatic partially methylated domains (PMDs). In contrast, reads originating from euchromatic regions and gene bodies have more ordered DNA methylation patterns. By analysing the patterns of single molecules in more detail, we show the existence of a 185bp periodicity in DNA methylation that accounts for some of the heterogeneity we uncover in long single-molecule DNA methylation patterns. We find that this periodic structure is partially masked in bulk data in a manner that is consistent with imperfect phasing of nucleosomes between molecules. Our findings demonstrate the power of single-molecule analysis of long-read data to understand the structure of the human methylome.

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. 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) two-layer size control strategies such as sizer-timer or adder-timer are ideal because they maintain constant mean concentrations whilst minimising concentration noise; (iv) accurate concentration homeostasis requires a fine tuning of dosage compensation, replication timing, and size-dependent gene expression; (v) deviations from perfect concentration homeostasis show up as deviations of the concentration distribution from a gamma distribution. Some of these predictions are confirmed using data for E. coli, fission yeast, and budding yeast.

Non-Markovian 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-Markovian models by the solutions of much simpler time-inhomogeneous Markovian 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-Markovian model. We show using a variety of models, where the delays stem from transcriptional processes and feedback control, that the Markovian models learnt by the neural network accurately reflect the stochastic dynamics across parameter space. Cells are complex systems that make decisions biologists struggle to understand. Here, the authors use neural networks to approximate the solution of mathematical models that capture the history and randomness of biochemical processes in order to understand the principles of transcription control.

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.

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.

Temporal variation of environmental stimuli leads to changes in gene expression. Since the latter is noisy and since many reaction events occur between the birth and death of an mRNA molecule, it is of interest to understand how a stimulus affects the transcript numbers measured at various sub-cellular locations. Here, we construct a stochastic model describing the dynamics of signal-dependent gene expression and its propagation downstream of transcription. For any time-dependent stimulus and assuming bursty gene expression, we devise a procedure which allows us to obtain time-dependent distributions of mRNA numbers at various stages of its life-cycle, e.g. in its nascent form at the transcription site, post-splicing in the nucleus, and after it is exported to the cytoplasm. We also derive an expression for the error in the approximation whose accuracy is verified via stochastic simulation. We find that, depending on the frequency of oscillation and the time of measurement, a stimulus can lead to cytoplasmic amplification or attenuation of transcriptional noise.

Auto-regulatory feedback loops are one of the most common network motifs. A wide variety of stochastic models have been constructed to understand how the fluctuations in protein numbers in these loops are influenced by the kinetic parameters of the main biochemical steps. These models differ according to (i) which sub-cellular processes are explicitly modelled; (ii) the modelling methodology employed (discrete, continuous or hybrid); (iii) whether they can be analytically solved for the steady-state distribution of protein numbers. We discuss the assumptions and properties of the main models in the literature, summarize our current understanding of the relationship between them and highlight some of the insights gained through modelling., 12 pages, 3 figures. Submitted to Biophysical Journal

Significance The random nature of gene expression is well established experimentally. Mathematical modeling provides a means of understanding the factors leading to the observed stochasticity. In this article, we extend the classical two-state model of stochastic mRNA dynamics to include a considerable number of salient features of single-cell biology, such as cell division, replication, mRNA maturation, dosage compensation, and growth-dependent transcription. By means of biologically relevant approximations, we obtain expressions for the time-dependent distributions of mRNA and protein numbers. These provide insight into how fluctuations are modified and controlled by complex intracellular processes., The stochasticity of gene expression presents significant challenges to the modeling of genetic networks. A two-state model describing promoter switching, transcription, and messenger RNA (mRNA) decay is the standard model of stochastic mRNA dynamics in eukaryotic cells. Here, we extend this model to include mRNA maturation, cell division, gene replication, dosage compensation, and growth-dependent transcription. We derive expressions for the time-dependent distributions of nascent mRNA and mature mRNA numbers, provided two assumptions hold: 1) nascent mRNA dynamics are much faster than those of mature mRNA; and 2) gene-inactivation events occur far more frequently than gene-activation events. We confirm that thousands of eukaryotic genes satisfy these assumptions by using data from yeast, mouse, and human cells. We use the expressions to perform a sensitivity analysis of the coefficient of variation of mRNA fluctuations averaged over the cell cycle, for a large number of genes in mouse embryonic stem cells, identifying degradation and gene-activation rates as the most sensitive parameters. Furthermore, it is shown that, despite the model’s complexity, the time-dependent distributions predicted by our model are generally well approximated by the negative binomial distribution. Finally, we extend our model to include translation, protein decay, and auto-regulatory feedback, and derive expressions for the approximate time-dependent protein-number distributions, assuming slow protein decay. Our expressions enable us to study how complex biological processes contribute to the fluctuations of gene products in eukaryotic cells, as well as allowing a detailed quantitative comparison with experimental data via maximum-likelihood methods.

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.

Maintenance and regeneration of adult tissues rely on the self-renewal of stem cells. Regeneration without over-proliferation requires precise regulation of the stem cell proliferation and differentiation rates. The nature of such regulatory mechanisms in different tissues, and how to incorporate them in models of stem cell population dynamics, is incompletely understood. The critical birth-death (CBD) process is widely used to model stem cell populations, capturing key phenomena, such as scaling laws in clone size distributions. However, the CBD process neglects regulatory mechanisms. Here, we propose the birth-death process with volume exclusion (vBD), a variation of the birth-death process that considers crowding effects, such as may arise due to limited space in a stem cell niche. While the deterministic rate equations predict a single non-trivial attracting steady state, the master equation predicts extinction and transient distributions of stem cell numbers with three possible behaviours: long-lived quasi-steady state, and short-lived bimodal or unimodal distributions. In all cases, we approximate solutions to the vBD master equation using a renormalised system-size expansion, quasi-steady state approximation and the WKB method. Our study suggests that the size distribution of a stem cell population bears signatures that are useful to detect negative feedback mediated via volume exclusion., Comment: 23 pages, 5 figures

The accurate establishment and maintenance of DNA methylation patterns is vital for mammalian development and disruption to these processes causes human disease. Our understanding of DNA methylation mechanisms has been facilitated by mathematical modelling, particularly stochastic simulations. Megabase-scale variation in DNA methylation patterns is observed in development, cancer and ageing and the mechanisms generating these patterns are little understood. However, the computational cost of stochastic simulations prevents them from modelling such large genomic regions. Here, we test the utility of three different mean-field models to predict summary statistics associated with large-scale DNA methylation patterns. By comparison to stochastic simulations, we show that a cluster mean-field model accurately predicts the statistical properties of steady-state DNA methylation patterns, including the mean and variance of methylation levels calculated across a system of CpG sites, as well as the covariance and correlation of methylation levels between neighbouring sites. We also demonstrate that a cluster mean-field model can be used within an approximate Bayesian computation framework to accurately infer model parameters from data. As mean-field models can be solved numerically in a few seconds, our work demonstrates their utility for understanding the processes underpinning large-scale DNA methylation patterns.

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.

Funding Information: S.B. and R.G. were supported by a Leverhulme Trust grant no. (RPG-2018-423). J.H. was supported by a BBSRC EASTBIO PhD studentship. Two-state models (telegraph-like models) have a successful history of predicting distributions of cellular and nascent mRNA numbers that can well fit experimental data. These models exclude key rate limiting steps, and hence it is unclear why they are able to accurately predict the number distributions. To answer this question, here we compare these models to a novel stochastic mechanistic model of transcription in mammalian cells that presents a unified description of transcriptional factor, polymerase and mature mRNA dynamics. We show that there is a large region of parameter space where the first, second and third moments of the distributions of the waiting times between two consecutively produced transcripts (nascent or mature) of two-state and mechanistic models exactly match. In this region: (i) one can uniquely express the two-state model parameters in terms of those of the mechanistic model, (ii) the models are practically indistinguishable by comparison of their transcript numbers distributions, and (iii) they are distinguishable from the shape of their waiting time distributions. Our results clarify the relationship between different gene expression models and identify a means to select between them from experimental data. Publisher PDF

We consider a stochastic model where a gene switches between two states, an mRNA transcript is released in the active state and subsequently it undergoes an arbitrary number of sequential unimolecular steps before being degraded. The reactions effectively describe various stages of the mRNA life cycle such as initiation, elongation, termination, splicing, export and degradation. We construct a novel mean-field approach that leads to closed-form steady-state distributions for the number of transcript molecules at each stage of the mRNA life cycle. By comparison with stochastic simulations, we show that the approximation is highly accurate over all of parameter space, independent of the type of expression (constitutive or bursty) and of the shape of the distribution (unimodal, bimodal and nearly bimodal). The theory predicts that in a population of identical cells, any bimodality is gradually washed away as the mRNA progresses through its life cycle.

MomentClosure.jl is a Julia package providing automated derivation of the time-evolution equations of the moments of molecule numbers for virtually any chemical reaction network using a wide range of moment closure approximations. It extends the capabilities of modelling stochastic biochemical systems in Julia and can be particularly useful when exact analytic solutions of the chemical master equation are unavailable and when Monte Carlo simulations are computationally expensive. MomentClosure.jl is freely accessible under the MIT license. Source code and documentation are available at https://github.com/augustinas1/MomentClosure.jl, Comment: 2 pages, 1 figure

Unlike many single-celled organisms, the growth of fission yeast cells within a cell cycle is not exponential. It is rather characterized by three distinct phases (elongation, septation and fission), each with a different growth rate. Experiments also show that the distribution of cell size in a lineage is often bimodal, unlike the unimodal distributions measured for the bacterium Escherichia coli. Here we construct a detailed stochastic model of cell size dynamics in fission yeast. The theory leads to analytic expressions for the cell size and the birth size distributions, and explains the origin of bimodality seen in experiments. In particular our theory shows that the left peak in the bimodal distribution is associated with cells in the elongation phase while the right peak is due to cells in the septation and fission phases. We show that the size control strategy, the variability in the added size during a cell cycle and the fraction of time spent in each of the three cell growth phases have a strong bearing on the shape of the cell size distribution. Furthermore we infer all the parameters of our model by matching the theoretical cell size and birth size distributions to those from experimental single cell time-course data for seven different growth conditions. Our method provides a much more accurate means of determining the cell size control strategy (timer, adder or sizer) than the standard method based on the slope of the best linear fit between the birth and division sizes. We also show that the variability in added size and the strength of cell size control of fission yeast depend weakly on the temperature but strongly on the culture medium.Author summaryAdvances in microscopy enable us to follow single cells over long timescales from which we can understand how their size varies with time and the nature of innate strategies developed to control cell size. This data shows that in many cell types growth is exponential and the distribution of cell sizes has one peak, namely there is a single characteristic cell size. However data for fission yeast shows remarkable differences: growth is non-exponential and the distribution of cell sizes has two peaks, meaning two characteristic cell sizes exist. Here we construct the first mathematical model of this organism; by solving the model analytically we show that it is able to predict the two peaked distributions of cell size seen in data and provides an explanation for each peak in terms of the various growth phases of the single-celled organism. Furthermore by fitting the model to the data, we infer values for the rates of all microscopic processes in our model. This method is shown to provide a much more reliable inference than current methods and sheds light on how the strategy used by fission yeast cells to control their size varies with external conditions.

Propensity functions of the Hill-type are commonly used to model transcriptional regulation in stochastic models of gene expression. This leads to an effective reduced master equation for the mRNA and protein dynamics only. Based on deterministic considerations, it is often stated or tacitly assumed that such models are valid in the limit of rapid promoter switching. Here, starting from the chemical master equation describing promoter-protein interactions, mRNA transcription, protein translation and decay, we prove that in the limit of fast promoter switching, the distribution of protein numbers is different than that given by standard stochastic models with Hill-type propensities. We show the differences are pronounced whenever the protein-DNA binding rate is much larger than the unbinding rate, a special case of fast promoter switching. Furthermore we show using both theory and simulations that use of the standard stochastic models leads to drastically incorrect predictions for the switching properties of positive feedback loops and that these differences decrease with increasing mean protein burst size. Our results confirm that commonly used stochastic models of gene regulatory networks are only accurate in a subset of the parameter space consistent with rapid promoter switching.Statement of SignificanceA large number of models of gene regulatory networks in the literature assume that since promoter switching is fast then transcriptional regulation can be effectively modeled using Hill functions. While this approach can be rigorously justified for deterministic models, it is presently unclear if it is also the case for stochastic models. In this article we prove that this is not the case, i.e. stochastic models of gene regulatory systems, namely those with feedback loops, describing transcriptional regulation using Hill functions are only valid in a subset of parameter conditions consistent with fast promoter switching. We identify parameter regimes where these models are correct and where their predictions cannot be trusted.