Stochastic models of gene transcription with upstream drives: exact solution and sample path characterization

Moment-based parameter inference with error guarantees for stochastic reaction networks

Inferring parameters of biochemical kinetic models from single-cell data remains challenging because of the uncertainty arising from the intractability of the likelihood function of stochastic reaction networks. Such uncertainty falls beyond current error quantification measures, which focus on the effects of finite sample size and identifiability but lack theoretical guarantees when likelihood approximations are needed. Here, we propose a method for the inference of parameters of stochastic reaction networks that works for both steady-state and time-resolved data and is applicable to networks with non-linear and rational propensities. Our approach provides bounds on the parameters via convex optimization over sets constrained by moment equations and moment matrices by taking observations to form moment intervals, which are then used to constrain parameters through convex sets. The bounds on the parameters contain the true parameters under the condition that the moment intervals contain the true moments, thus providing uncertainty quantification and error guarantees. Our approach does not need to predict moments and distributions for given parameters (i.e., it avoids solving or simulating the forward problem) and hence circumvents intractable likelihood computations or computationally expensive simulations. We demonstrate its use for uncertainty quantification, data integration, and prediction of latent species statistics through synthetic data from common non-linear biochemical models including the Schlögl model and the toggle switch, a model of post-transcriptional regulation at steady state, and a birth-death model with time-dependent data.

Incorporating spatial diffusion into models of bursty stochastic transcription

The dynamics of gene expression are stochastic and spatial at the molecular scale, with messenger RNA (mRNA) transcribed at specific nuclear locations and then transported to the nuclear boundary for export. Consequently, the spatial distributions of these molecules encode their underlying dynamics. While mechanistic models for molecular counts have revealed numerous insights into gene expression, they have largely neglected now-available subcellular spatial resolution down to individual molecules. Owing to the technical challenges inherent in spatial stochastic processes, tools for studying these subcellular spatial patterns are still limited. Here, we introduce a spatial stochastic model of nuclear mRNA with two-state (telegraph) transcriptional dynamics. Observations of the model can be concisely described as following a spatial Cox process driven by a stochastically switching partial differential equation. We derive analytical solutions for spatial and demographic moments and validate them with simulations. We show that the distribution of mRNA counts can be accurately approximated by a Poisson-beta distribution with tractable parameters, even with complex spatial dynamics. This observation allows for efficient parameter inference demonstrated on synthetic data. Altogether, our work adds progress towards a new frontier of subcellular spatial resolution in inferring the dynamics of gene expression from static snapshot data.

Quantifying and correcting bias in transcriptional parameter inference from single-cell data

The snapshot distribution of mRNA counts per cell can be measured using single-molecule fluorescence in situ hybridization or single-cell RNA sequencing. These distributions are often fit to the steady-state distribution of the two-state telegraph model to estimate the three transcriptional parameters for a gene of interest: mRNA synthesis rate, the switching on rate (the on state being the active transcriptional state), and the switching off rate. This model assumes no extrinsic noise, i.e., parameters do not vary between cells, and thus estimated parameters are to be understood as approximating the average values in a population. The accuracy of this approximation is currently unclear. Here, we develop a theory that explains the size and sign of estimation bias when inferring parameters from single-cell data using the standard telegraph model. We find specific bias signatures depending on the source of extrinsic noise (which parameter is most variable across cells) and the mode of transcriptional activity. If gene expression is not bursty then the population averages of all three parameters are overestimated if extrinsic noise is in the synthesis rate; underestimation occurs if extrinsic noise is in the switching on rate; both underestimation and overestimation can occur if extrinsic noise is in the switching off rate. We find that some estimated parameters tend to infinity as the size of extrinsic noise approaches a critical threshold. In contrast when gene expression is bursty, we find that in all cases the mean burst size (ratio of the synthesis rate to the switching off rate) is overestimated while the mean burst frequency (the switching on rate) is underestimated. We estimate the size of extrinsic noise from the covariance matrix of sequencing data and use this together with our theory to correct published estimates of transcriptional parameters for mammalian genes.

Poisson representation: a bridge between discrete and continuous models of stochastic gene regulatory networks

Stochastic gene expression dynamics can be modelled either discretely or continuously. Previous studies have shown that the mRNA or protein number distributions of some simple discrete and continuous gene expression models are related by Gardiner's Poisson representation. Here, we systematically investigate the Poisson representation in complex stochastic gene regulatory networks. We show that when the gene of interest is unregulated, the discrete and continuous descriptions of stochastic gene expression are always related by the Poisson representation, no matter how complex the model is. This generalizes the results obtained in Dattani & Barahona (Dattani & Barahona 2017 J. R. Soc. Interface 14, 20160833 (doi:10.1098/rsif.2016.0833)). In addition, using a simple counter-example, we find that the Poisson representation in general fails to link the two descriptions when the gene is regulated. However, for a general stochastic gene regulatory network, we demonstrate that the discrete and continuous models are approximately related by the Poisson representation in the limit of large protein numbers. These theoretical results are further applied to analytically solve many complex gene expression models whose exact distributions are previously unknown.

Studying stochastic systems biology of the cell with single-cell genomics data

Recent experimental developments in genome-wide RNA quantification hold considerable promise for systems biology. However, rigorously probing the biology of living cells requires a unified mathematical framework that accounts for single-molecule biological stochasticity in the context of technical variation associated with genomics assays. We review models for a variety of RNA transcription processes, as well as the encapsulation and library construction steps of microfluidics-based single-cell RNA sequencing, and present a framework to integrate these phenomena by the manipulation of generating functions. Finally, we use simulated scenarios and biological data to illustrate the implications and applications of the approach.

"Bet hedging" against climate change in developing and adult animals: roles for stochastic gene expression, phenotypic plasticity, epigenetic inheritance and adaptation

Animals from embryos to adults experiencing stress from climate change have numerous mechanisms available for enhancing their long-term survival. In this review we consider these options, and how viable they are in a world increasingly experiencing extreme weather associated with climate change. A deeply understood mechanism involves natural selection, leading to evolution of new adaptations that help cope with extreme and stochastic weather events associated with climate change. While potentially effective at staving off environmental challenges, such adaptations typically occur very slowly and incrementally over evolutionary time. Consequently, adaptation through natural selection is in most instances regarded as too slow to aid survival in rapidly changing environments, especially when considering the stochastic nature of extreme weather events associated with climate change. Alternative mechanisms operating in a much shorter time frame than adaptation involve the rapid creation of alternate phenotypes within a life cycle or a few generations. Stochastic gene expression creates multiple phenotypes from the same genotype even in the absence of environmental cues. In contrast, other mechanisms for phenotype change that are externally driven by environmental clues include well-understood developmental phenotypic plasticity (variation, flexibility), which can enable rapid, within-generation changes. Increasingly appreciated are epigenetic influences during development leading to rapid phenotypic changes that can also immediately be very widespread throughout a population, rather than confined to a few individuals as in the case of favorable gene mutations. Such epigenetically-induced phenotypic plasticity can arise rapidly in response to stressors within a generation or across a few generations and just as rapidly be "sunsetted" when the stressor dissipates, providing some capability to withstand environmental stressors emerging from climate change. Importantly, survival mechanisms resulting from adaptations and developmental phenotypic plasticity are not necessarily mutually exclusive, allowing for classic "bet hedging". Thus, the appearance of multiple phenotypes within a single population provides for a phenotype potentially optimal for some future environment. This enhances survival during stochastic extreme weather events associated with climate change. Finally, we end with recommendations for future physiological experiments, recommending in particular that experiments investigating phenotypic flexibility adopt more realistic protocols that reflect the stochastic nature of weather.

Studying stochastic systems biology of the cell with single-cell genomics data

Recent experimental developments in genome-wide RNA quantification hold considerable promise for systems biology. However, rigorously probing the biology of living cells requires a unified mathematical framework that accounts for single-molecule biological stochasticity in the context of technical variation associated with genomics assays. We review models for a variety of RNA transcription processes, as well as the encapsulation and library construction steps of microfluidics-based single-cell RNA sequencing, and present a framework to integrate these phenomena by the manipulation of generating functions. Finally, we use simulated scenarios and biological data to illustrate the implications and applications of the approach.

Gene expression model inference from snapshot RNA data using Bayesian non-parametrics

Gene expression models, which are key towards understanding cellular regulatory response, underlie observations of single-cell transcriptional dynamics. Although RNA expression data encode information on gene expression models, existing computational frameworks do not perform simultaneous Bayesian inference of gene expression models and parameters from such data. Rather, gene expression models-composed of gene states, their connectivities and associated parameters-are currently deduced by pre-specifying gene state numbers and connectivity before learning associated rate parameters. Here we propose a method to learn full distributions over gene states, state connectivities and associated rate parameters, simultaneously and self-consistently from single-molecule RNA counts. We propagate noise from fluctuating RNA counts over models by treating models themselves as random variables. We achieve this within a Bayesian non-parametric paradigm. We demonstrate our method on the Escherichia coli lacZ pathway and the Saccharomyces cerevisiae STL1 pathway, and verify its robustness on synthetic data.

Interpretable and tractable models of transcriptional noise for the rational design of single-molecule quantification experiments

The question of how cell-to-cell differences in transcription rate affect RNA count distributions is fundamental for understanding biological processes underlying transcription. Answering this question requires quantitative models that are both interpretable (describing concrete biophysical phenomena) and tractable (amenable to mathematical analysis). This enables the identification of experiments which best discriminate between competing hypotheses. As a proof of principle, we introduce a simple but flexible class of models involving a continuous stochastic transcription rate driving a discrete RNA transcription and splicing process, and compare and contrast two biologically plausible hypotheses about transcription rate variation. One assumes variation is due to DNA experiencing mechanical strain, while the other assumes it is due to regulator number fluctuations. We introduce a framework for numerically and analytically studying such models, and apply Bayesian model selection to identify candidate genes that show signatures of each model in single-cell transcriptomic data from mouse glutamatergic neurons.

Mother–Fetus Immune Cross-Talk Coordinates “Extrinsic”/“Intrinsic” Embryo Gene Expression Noise and Growth Stability

Developmental instability (DI) is thought to be inversely related to a capacity of an organism to buffer its development against random genetic and environmental perturbations. DI is represented by a trait’s inter- and intra-individual variabilities. The inter-individual variability (inversely referred to as canalization) indicates the capability of organisms to reproduce a trait from individual to individual. The intra-individual variability reflects an organism’s capability to stabilize a trait internally under the same conditions, and, for symmetric traits, it is expressed as fluctuating asymmetry (FA). When representing a trait as a random variable conditioned on environmental fluctuations, it is clear that, in statistical terms, the DI partitions into “extrinsic” (canalization) and “intrinsic” (FA) components of a trait’s variance/noise. We established a simple statistical framework to dissect both parts of a symmetric trait variance/noise using a PCA (principal component analysis) projection of the left/right measurements on eigenvectors followed by GAMLSS (generalized additive models for location scale and shape) modeling of eigenvalues. The first eigenvalue represents “extrinsic” and the second—“intrinsic” DI components. We applied this framework to investigate the impact of mother–fetus major histocompatibility complex (MHC)-mediated immune cross-talk on gene expression noise and developmental stability. We showed that “intrinsic” gene noise for the entire transcriptional landscape could be estimated from a small subset of randomly selected genes. Using a diagnostic set of genes, we found that allogeneic MHC combinations tended to decrease “extrinsic” and “intrinsic” gene noise in C57BL/6J embryos developing in the surrogate NOD-SCID and BALB/c mothers. The “intrinsic” gene noise was negatively correlated with growth (embryonic mass) and the levels of placental growth factor (PLGF), but not vascular endothelial growth factor (VEGF). However, it was positively associated with phenotypic growth instability and noise in PLGF. In mammals, the mother–fetus MHC interaction plays a significant role in development, contributing to the fitness of the offspring. Our results demonstrate that a positive impact of distant MHC combinations on embryonic growth could be mediated by the reduction of “intrinsic” gene noise followed by the developmental stabilization of growth.

Quantifying how post-transcriptional noise and gene copy number variation bias transcriptional parameter inference from mRNA distributions

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. While 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.

RNA velocity unraveled

We perform a thorough analysis of RNA velocity methods, with a view towards understanding the suitability of the various assumptions underlying popular implementations. In addition to providing a self-contained exposition of the underlying mathematics, we undertake simulations and perform controlled experiments on biological datasets to assess workflow sensitivity to parameter choices and underlying biology. Finally, we argue for a more rigorous approach to RNA velocity, and present a framework for Markovian analysis that points to directions for improvement and mitigation of current problems.

Epigenetic fluctuations underlie gene expression timescales and variability

Isogenic populations of mammalian cells exhibit significant gene expression variability. This variability can be separated into two sources, cis, or allele-specific sources, and trans and global processes. Furthermore, each source of variability has its own timescale. Fast timescales will result in rapid fluctuation of gene expression whereas slow timescales will result in longer persistence of gene expression levels over time. Here we investigated sources of gene expression that are intrinsic, i.e. coming from cis-regulatory factors and follow slow timescales. To do so, we developed a reporter system that isolates allele-specific variability and measures its persistence in imaging and long-term fluctuation analysis experiments. Our results identify a new source of gene expression variability that is allele-specific but that is fluctuating on timescales of days. We hypothesized that allele-specific fluctuations of epigenetic regulatory factors are responsible for the newly discovered allele-specific and slow source of gene expression variability. Using mathematical modeling we showed that the addition of this effect to the two-state model is sufficient to account for all empirical observation. Furthermore, using direct assays of chromatin markers we find fluctuation in H3K4me3 levels that match the observed changes in gene expression levels providing direct experimental support of our model. Collectively, our work shows that slow fluctuations of regulatory chromatin modifications contribute to the variability in gene expression.

Modeling bursty transcription and splicing with the chemical master equation

Splicing cascades that alter gene products posttranscriptionally also affect expression dynamics. We study a class of processes and associated distributions that emerge from models of bursty promoters coupled to directed acyclic graphs of splicing. These solutions provide full time-dependent joint distributions for an arbitrary number of species with general noise behaviors and transient phenomena, offering qualitative and quantitative insights about how splicing can regulate expression dynamics. Finally, we derive a set of quantitative constraints on the minimum complexity necessary to reproduce gene coexpression patterns using synchronized burst models. We validate these findings by analyzing long-read sequencing data, where we find evidence of expression patterns largely consistent with these constraints.

Toll-Like Receptor 4 in Pain: Bridging Molecules-to-Cells-to-Systems

Pain impacts the lives of billions of people around the world - both directly and indirectly. It is complex and transcends beyond an unpleasant sensory experience to encompass emotional experiences. To date, there are no successful treatments for sufferers of chronic pain. Although opioids do not provide any benefit to chronic pain sufferers, they are still prescribed, often resulting in more complications such as hyperalgesia and dependence. In order to develop effective and safe medications to manage, and perhaps even treat pain, it is important to evaluate novel contributors to pain pathologies. As such, in this chapter we review the role of Toll-like receptor 4, a receptor of the innate immune system, that continues to gain substantial attention in the field of pain research. Positioned in the nexus of the neuro and immune systems, TLR4 may provide one of the missing pieces in understanding the complexities of pain. Here we consider how TLR4 enables a mechanistical understanding of pain as a multidimensional biopsychosocial state from molecules to cells to systems and back again.

Mean-field theory accurately captures the variation of copy number distributions across the mRNA life cycle

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 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 the 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.

Approximating gene transcription dynamics using steady-state formulas

Understanding how genes in a single cell respond to dynamically changing signals has been a central question in stochastic gene transcription research. Recent studies have generated massive steady-state or snapshot mRNA distribution data of individual cells, and inferred a large spectrum of kinetic transcription parameters under varying conditions. However, there have been few algorithms to convert these static data into the temporal variation of kinetic rates. Real-time imaging has been developed to monitor stochastic transcription processes at the single-cell level, but the immense technicality has prevented its application to most endogenous loci in mammalian cells. In this article, we introduced a stochastic gene transcription model with variable kinetic rates induced by unstable cellular conditions. We approximated the transcription dynamics using easily obtained steady-state formulas in the model. We tested the approximation against experimental data in both prokaryotic and eukaryotic cells and further solidified the conditions that guarantee the robustness of the method. The method can be easily implemented to provide convenient tools for quantifying dynamic kinetics and mechanisms underlying the widespread static transcription data, and may shed a light on circumventing the limitation of current bursting data on transcriptional real-time imaging.

A Hill type equation can predict target gene expression driven by p53 pulsing

Many factors determine target gene expression dynamics under p53 pulsing. In this study, I sought to determine the mechanism by which duration, frequency, binding affinity and maximal transcription rate affect the expression dynamics of target genes. Using an analytical method to solve a simple model, I found that the fold change of target gene expression increases relative to the number of p53 pulses, and the optimal frequency, 0.18 h-1 , from two real p53 pulses drives the maximal fold change with a decay rate of 0.18 h-1 . Moreover, p53 pulses may also lead to a higher fold change than sustained p53. Finally, I discovered that a Hill-type equation, including these effect factors, can characterise target gene expression. The average error between the theoretical predictions and experiments was 23%. Collectively, this equation advances the understanding of transcription factor dynamics, where duration and frequency play a significant role in the fine regulation of target gene expression with higher binding affinity.

Mutational pressure by host APOBEC3s more strongly affects genes expressed early in the lytic phase of herpes simplex virus-1 (HSV-1) and human polyomavirus (HPyV) infection

Herpes-Simplex Virus 1 (HSV-1) infects most humans when they are young, sometimes with fatal consequences. Gene expression occurs in a temporal order upon lytic HSV-1 infection: immediate early (IE) genes are expressed, then early (E) genes, followed by late (L) genes. During this infection cycle, the HSV-1 genome has the potential for exposure to APOBEC3 (A3) proteins, a family of cytidine deaminases that cause C>U mutations on single-stranded DNA (ssDNA), often resulting in a C>T transition. We developed a computational model for the mutational pressure of A3 on the lytic cycle of HSV-1 to determine which viral kinetic gene class is most vulnerable to A3 mutations. Using in silico stochastic methods, we simulated the infectious cycle under varying intensities of A3 mutational pressure. We found that the IE and E genes are more vulnerable to A3 than L genes. We validated this model by analyzing the A3 evolutionary footprints in 25 HSV-1 isolates. We find that IE and E genes have evolved to underrepresent A3 hotspot motifs more so than L genes, consistent with greater selection pressure on IE and E genes. We extend this model to two-step infections, such as those of polyomavirus, and find that the same pattern holds for over 25 human Polyomavirus (HPyVs) genomes. Genes expressed earlier during infection are more vulnerable to mutations than those expressed later.

Transcriptional bursts explain autosomal random monoallelic expression and affect allelic imbalance

Transcriptional bursts render substantial biological noise in cellular transcriptomes. Here, we investigated the theoretical extent of allelic expression resulting from transcriptional bursting and how it compared to the amount biallelic, monoallelic and allele-biased expression observed in single-cell RNA-sequencing (scRNA-seq) data. We found that transcriptional bursting can explain the allelic expression patterns observed in single cells, including the frequent observations of autosomal monoallelic gene expression. Importantly, we identified that the burst frequency largely determined the fraction of cells with monoallelic expression, whereas the burst size had little effect on monoallelic observations. The high consistency between the bursting model predictions and scRNA-seq observations made it possible to assess the heterogeneity of a group of cells as their deviation in allelic observations from the expected. Finally, both burst frequency and size contributed to allelic imbalance observations and reinforced that studies of allelic imbalance can be confounded from the inherent noise in transcriptional bursting. Altogether, we demonstrate that allele-level transcriptional bursting renders widespread, although predictable, amounts of monoallelic and biallelic expression in single cells and cell populations.

Genes are transcribed into RNA and further translated into proteins. The maternal and paternal copy of each gene are typically transcribed independently, and transcription itself occur in discrete stochastic bursts (transcriptional bursts). Pioneering single-cell analysis of RNA across cells revealed abundant fluctuations in the amounts of maternal and paternal RNA in cells, with frequent observations of RNA from only the maternal or paternal gene copy (monoallelic expression). In this study, we investigated to which extent the observed monoallelic expression across single cells can be explained by transcriptional bursting. We demonstrate that the process of transcriptional bursting is sufficient to explain the amount of monoallelic expression, and we further demonstrate that the frequency of bursts mainly determines the frequency of monoallelic observations. Furthermore, we show that transcriptional bursts may lead to false positive observations of monoallelic expression across cell populations. Therefore, stochastic transcription renders large fluctuations in allelic origin of RNA in cells over time, including frequent monoallelic observations when profiling single cells.

The circadian oscillator analysed at the single-transcript level

The circadian clock is an endogenous and self-sustained oscillator that anticipates daily environmental cycles. While rhythmic gene expression of circadian genes is well-described in populations of cells, the single-cell mRNA dynamics of multiple core clock genes remain largely unknown. Here we use single-molecule fluorescence in situ hybridisation (smFISH) at multiple time points to measure pairs of core clock transcripts, Rev-erbα (Nr1d1), Cry1 and Bmal1, in mouse fibroblasts. The mean mRNA level oscillates over 24 h for all three genes, but mRNA numbers show considerable spread between cells. We develop a probabilistic model for multivariate mRNA counts using mixtures of negative binomials, which accounts for transcriptional bursting, circadian time and cell-to-cell heterogeneity, notably in cell size. Decomposing the mRNA variability into distinct noise sources shows that clock time contributes a small fraction of the total variability in mRNA number between cells. Thus, our results highlight the intrinsic biological challenges in estimating circadian phase from single-cell mRNA counts and suggest that circadian phase in single cells is encoded post-transcriptionally.

Computation of Single-Cell Metabolite Distributions Using Mixture Models

Metabolic heterogeneity is widely recognized as the next challenge in our understanding of non-genetic variation. A growing body of evidence suggests that metabolic heterogeneity may result from the inherent stochasticity of intracellular events. However, metabolism has been traditionally viewed as a purely deterministic process, on the basis that highly abundant metabolites tend to filter out stochastic phenomena. Here we bridge this gap with a general method for prediction of metabolite distributions across single cells. By exploiting the separation of time scales between enzyme expression and enzyme kinetics, our method produces estimates for metabolite distributions without the lengthy stochastic simulations that would be typically required for large metabolic models. The metabolite distributions take the form of Gaussian mixture models that are directly computable from single-cell expression data and standard deterministic models for metabolic pathways. The proposed mixture models provide a systematic method to predict the impact of biochemical parameters on metabolite distributions. Our method lays the groundwork for identifying the molecular processes that shape metabolic heterogeneity and its functional implications in disease.

Application of the Goodwin model to autoregulatory feedback for stochastic gene expression

In this paper we analyse stochastic expression of a single gene with its dynamics given by the classical Goodwin model with mRNA and protein contribution. We compare the effect of the presence of positive and negative feedback on the transcription regulation. In such cases we observe two qualitatively different types of asymptotic behaviour. In the case of a negative feedback loop, under sufficient conditions, one can find a stationary density for mRNA and protein molecules. In the case of a positive feedback loop we observe extinction of both types of molecules with time.

Order Through Disorder: The Characteristic Variability of Systems

Randomness characterizes many processes in nature, and therefore its importance cannot be overstated. In the present study, we investigate examples of randomness found in various fields, to underlie its fundamental processes. The fields we address include physics, chemistry, biology (biological systems from genes to whole organs), medicine, and environmental science. Through the chosen examples, we explore the seemingly paradoxical nature of life and demonstrate that randomness is preferred under specific conditions. Furthermore, under certain conditions, promoting or making use of variability-associated parameters may be necessary for improving the function of processes and systems.

Extrinsic Noise and Heavy-Tailed Laws in Gene Expression

Noise in gene expression is one of the hallmarks of life at the molecular scale. Here we derive analytical solutions to a set of models describing the molecular mechanisms underlying transcription of DNA into RNA. Our ansatz allows us to incorporate the effects of extrinsic noise-encompassing factors external to the transcription of the individual gene-and discuss the ramifications for heterogeneity in gene product abundance that has been widely observed in single cell data. Crucially, we are able to show that heavy-tailed distributions of RNA copy numbers cannot result from the intrinsic stochasticity in gene expression alone, but must instead reflect extrinsic sources of variability.

Exact solution of stochastic gene expression models with bursting, cell cycle and replication dynamics

The bulk of stochastic gene expression models in the literature do not have an explicit description of the age of a cell within a generation and hence they cannot capture events such as cell division and DNA replication. Instead, many models incorporate the cell cycle implicitly by assuming that dilution due to cell division can be described by an effective decay reaction with first-order kinetics. If it is further assumed that protein production occurs in bursts, then the stationary protein distribution is a negative binomial. Here we seek to understand how accurate these implicit models are when compared with more detailed models of stochastic gene expression. We derive the exact stationary solution of the chemical master equation describing bursty protein dynamics, binomial partitioning at mitosis, age-dependent transcription dynamics including replication, and random interdivision times sampled from Erlang or more general distributions; the solution is different for single lineage and population snapshot settings. We show that protein distributions are well approximated by the solution of implicit models (a negative binomial) when the mean number of mRNAs produced per cycle is low and the cell cycle length variability is large. When these conditions are not met, the distributions are either almost bimodal or else display very flat regions near the mode and cannot be described by implicit models. We also show that for genes with low transcription rates, the size of protein noise has a strong dependence on the replication time, it is almost independent of cell cycle variability for lineage measurements, and increases with cell cycle variability for population snapshot measurements. In contrast for large transcription rates, the size of protein noise is independent of replication time and increases with cell cycle variability for both lineage and population measurements.

Bayesian inference on stochastic gene transcription from flow cytometry data

Motivation

Transcription in single cells is an inherently stochastic process as mRNA levels vary greatly between cells, even for genetically identical cells under the same experimental and environmental conditions. We present a stochastic two-state switch model for the population of mRNA molecules in single cells where genes stochastically alternate between a more active ON state and a less active OFF state. We prove that the stationary solution of such a model can be written as a mixture of a Poisson and a Poisson-beta probability distribution. This finding facilitates inference for single cell expression data, observed at a single time point, from flow cytometry experiments such as FACS or fluorescence in situ hybridization (FISH) as it allows one to sample directly from the equilibrium distribution of the mRNA population. We hence propose a Bayesian inferential methodology using a pseudo-marginal approach and a recent approximation to integrate over unobserved states associated with measurement error.

Results

We provide a general inferential framework which can be widely used to study transcription in single cells from the kind of data arising in flow cytometry experiments. The approach allows us to separate between the intrinsic stochasticity of the molecular dynamics and the measurement noise. The methodology is tested in simulation studies and results are obtained for experimental multiple single cell expression data from FISH flow cytometry experiments.

Availability and implementation

All analyses were implemented in R. Source code and the experimental data are available at https://github.com/SimoneTiberi/Bayesian-inference-on-stochastic-gene-transcription-from-flow-cytometry-data.

Supplementary information

Supplementary data are available at Bioinformatics online.

Gene expression variability: the other dimension in transcriptome analysis

Transcriptome sequencing is a powerful technique to study molecular changes that underlie the differences in physiological conditions and disease progression. A typical question that is posed in such studies is finding genes with significant changes between sample groups. In this respect expression variability is regarded as a nuisance factor that is primarily of technical origin and complicates the data analysis. However, it is becoming apparent that the biological variation in gene expression might be an important molecular phenotype that can affect physiological parameters. In this review we explore the recent literature on technical and biological variability in gene expression, sources of expression variability, (epi-)genetic hallmarks, and evolutionary constraints in genes with robust and variable gene expression. We provide an overview of recent findings on effects of external cues, such as diet and aging, on expression variability and on other biological phenomena that can be linked to it. We discuss metrics and tools that were developed for quantification of expression variability and highlight the importance of future studies in this direction. To assist the adoption of expression variability analysis, we also provide a detailed description and computer code, which can easily be utilized by other researchers. We also provide a reanalysis of recently published data to highlight the value of the analysis method.

A multiscale model of epigenetic heterogeneity-driven cell fate decision-making

The inherent capacity of somatic cells to switch their phenotypic status in response to damage stimuli in vivo might have a pivotal role in ageing and cancer. However, how the entry-exit mechanisms of phenotype reprogramming are established remains poorly understood. In an attempt to elucidate such mechanisms, we herein introduce a stochastic model of combined epigenetic regulation (ER)-gene regulatory network (GRN) to study the plastic phenotypic behaviours driven by ER heterogeneity. To deal with such complex system, we additionally formulate a multiscale asymptotic method for stochastic model reduction, from which we derive an efficient hybrid simulation scheme. Our analysis of the coupled system reveals a regime of tristability in which pluripotent stem-like and differentiated steady-states coexist with a third indecisive state, with ER driving transitions between these states. Crucially, ER heterogeneity of differentiation genes is for the most part responsible for conferring abnormal robustness to pluripotent stem-like states. We formulate epigenetic heterogeneity-based strategies capable of unlocking and facilitating the transit from differentiation-refractory (stem-like) to differentiation-primed epistates. The application of the hybrid numerical method validates the likelihood of such switching involving solely kinetic changes in epigenetic factors. Our results suggest that epigenetic heterogeneity regulates the mechanisms and kinetics of phenotypic robustness of cell fate reprogramming. The occurrence of tunable switches capable of modifying the nature of cell fate reprogramming might pave the way for new therapeutic strategies to regulate reparative reprogramming in ageing and cancer.

Certain modifications of the structure and functioning of the protein/DNA complex called chromatin can allow adult, fully differentiated, cells to adopt a stem cell-like pluripotent state in a purely epigenetic manner, not involving changes in the underlying DNA sequence. Such reprogramming-like phenomena may constitute an innate reparative route through which human tissues respond to injury and could also serve as a novel regenerative strategy in human pathological situations in which tissue or organ repair is impaired. However, it should be noted that in vivo reprogramming would be capable of maintaining tissue homeostasis provided the acquisition of pluripotency features is strictly transient and accompanied by an accurate replenishment of the specific cell types being lost. Crucially, an excessive reprogramming in the absence of controlled re-differentiation would impair the repair or the replacement of damaged cells, thereby promoting pathological alterations of cell fate. A mechanistic understanding of how the degree of chromatin plasticity dictates the reparative versus pathological behaviour of in vivo reprogramming to rejuvenate aged tissues while preventing tumorigenesis is urgently needed, including especially the intrinsic epigenetic heterogeneity of the tissue resident cells being reprogrammed. We here introduce a novel method that mathematically captures how epigenetic heterogeneity is actually the driving force that governs the routes and kinetics to entry into and exit from a pathological stem-like state. Moreover, our approach computationally validates the likelihood of unlocking chronic, unrestrained plastic states and drive their differentiation down the correct path by solely manipulating the intensity and direction of few epigenetic control switches. Our approach could inspire new therapeutic approaches based on in vivo cell reprogramming for efficient tissue regeneration and rejuvenation and cancer treatment.

Stochastic modelling reveals mechanisms of metabolic heterogeneity

Phenotypic variation is a hallmark of cellular physiology. Metabolic heterogeneity, in particular, underpins single-cell phenomena such as microbial drug tolerance and growth variability. Much research has focussed on transcriptomic and proteomic heterogeneity, yet it remains unclear if such variation permeates to the metabolic state of a cell. Here we propose a stochastic model to show that complex forms of metabolic heterogeneity emerge from fluctuations in enzyme expression and catalysis. The analysis predicts clonal populations to split into two or more metabolically distinct subpopulations. We reveal mechanisms not seen in deterministic models, in which enzymes with unimodal expression distributions lead to metabolites with a bimodal or multimodal distribution across the population. Based on published data, the results suggest that metabolite heterogeneity may be more pervasive than previously thought. Our work casts light on links between gene expression and metabolism, and provides a theory to probe the sources of metabolite heterogeneity.

First-passage time statistics of stochastic transcription process for time-dependent reaction rates

Transcription in gene expression is an intrinsically noisy process which involves production and degradation of mRNAs. An important quantity to describe this stochastic process is the first-passage time (FPT), i.e., the time taken by mRNAs to reach a particular threshold. The process of transcription can be modelled as a simple birth-death process, assuming that the promoter is always in an active state and to encode the stochastic environment we consider the transcription rate to be time dependent. This generalization is suitable to capture bursty mRNA dynamics usually modelled as an ON-Off model and simplifies the calculation of FPT statistics for a cell population. We study the role of periodic modulation of the transcription rate on different moments of FPT distribution of a population of cells. Our calculation shows that for sinusoidal modulation there exists an extremal value of mean FPT as a function of the time period and phase of the transcription signal. However, for the square wave modulation of transcription rates simulation results show that the extremal value of the MFPT behaves monotonically with the variation of the phase.

Intrinsic and extrinsic noise of gene expression in lineage trees

Cell-to-cell heterogeneity is driven by stochasticity in intracellular reactions and the population dynamics. While these sources are usually studied separately, we develop an agent-based framework that accounts for both factors while tracking every single cell of a growing population. Apart from the common intrinsic variability, the framework also predicts extrinsic noise without the need to introduce fluctuating rate constants. Instead, extrinsic fluctuations are explained by cell cycle fluctuations and differences in cell age. We provide explicit formulas to quantify mean molecule numbers, intrinsic and extrinsic noise statistics in two-colour experiments. We find that these statistics differ significantly depending on the experimental setup used to observe the cells. We illustrate this fact using (i) averages over an isolated cell lineage tracked over many generations as observed in the mother machine, (ii) population snapshots with known cell ages as recorded in time-lapse microscopy, and (iii) snapshots with unknown cell ages as measured from static images or flow cytometry. Applying the method to models of stochastic gene expression and feedback regulation elucidates that isolated lineages, as compared to snapshot data, can significantly overestimate the mean number of molecules, overestimate extrinsic noise but underestimate intrinsic noise and have qualitatively different sensitivities to cell cycle fluctuations.

Sensitivity of asymmetric rate-dependent critical systems to initial conditions: Insights into cellular decision making

The work reported here aims to address the effects of time-dependent parameters and stochasticity on decision making in biological systems. We achieve this by extending previous studies that resorted to simple bifurcation normal forms, although in the present case we focus primarily on the issue of the system's sensitivity to initial conditions in the presence of two different noise distributions, Gaussian and Lévy. In addition, we also assess the impact of two-way sweeping at different rates through the critical region of a canonical Pitchfork bifurcation with a constant external asymmetry. The parallel with decision making in biocircuits is performed on this simple system since it is equivalent in its available states and dynamics to more complex genetic circuits published previously. Overall we verify that rate-dependent effects, previously reported as being important features of bifurcating systems, are specific to particular initial conditions. Processing of each starting state, which for the normal form underlying this study is akin to a classification task, is affected by the balance between sweeping speed through critical regions and the type of fluctuations added. For the heavy-tailed noise, two-way dynamic bifurcations are more efficient in processing the external signals, here understood to be jointly represented by the critical parameter profile and the external asymmetry amplitude, when compared to the system relying on escape dynamics. This is particular to the case when the system starts at an attractor not favored by the asymmetry and, in conjunction, when the sweeping amplitude is large.

Stabilization of dynamics of oscillatory systems by nonautonomous perturbation

Synchronization and stability under periodic oscillatory driving are well understood, but little is known about the effects of aperiodic driving, despite its abundance in nature. Here, we consider oscillators subject to driving with slowly varying frequency, and investigate both short-term and long-term stability properties. For a phase oscillator, we find that, counterintuitively, such variation is guaranteed to enlarge the Arnold tongue in parameter space. Using analytical and numerical methods that provide information on time-variable dynamical properties, we find that the growth of the Arnold tongue is specifically due to the growth of a region of intermittent synchronization where trajectories alternate between short-term stability and short-term neutral stability, giving rise to stability on average. We also present examples of higher-dimensional nonlinear oscillators where a similar stabilization phenomenon is numerically observed. Our findings help support the case that in general, deterministic nonautonomous perturbation is a very good candidate for stabilizing complex dynamics.

The Chemical Fluctuation Theorem governing gene expression

Gene expression is a complex stochastic process composed of numerous enzymatic reactions with rates coupled to hidden cell-state variables. Despite advances in single-cell technologies, the lack of a theory accurately describing the gene expression process has restricted a robust, quantitative understanding of gene expression variability among cells. Here we present the Chemical Fluctuation Theorem (CFT), providing an accurate relationship between the environment-coupled chemical dynamics of gene expression and gene expression variability. Combined with a general, accurate model of environment-coupled transcription processes, the CFT provides a unified explanation of mRNA variability for various experimental systems. From this analysis, we construct a quantitative model of transcription dynamics enabling analytic predictions for the dependence of mRNA noise on the mRNA lifetime distribution, confirmed against stochastic simulation. This work suggests promising new directions for quantitative investigation into cellular control over biological functions by making complex dynamics of intracellular reactions accessible to rigorous mathematical deductions.

A unified framework to understand gene expression noise is still lacking. Here the authors derive a universal theorem relating the biological noise with dynamics of birth and death processes and present a model of transcription dynamics, allowing analytical prediction of the dependence of mRNA noise on mRNA lifetime variability.

Inferring gene regulatory networks from single-cell data: a mechanistic approach

Background

The recent development of single-cell transcriptomics has enabled gene expression to be measured in individual cells instead of being population-averaged. Despite this considerable precision improvement, inferring regulatory networks remains challenging because stochasticity now proves to play a fundamental role in gene expression. In particular, mRNA synthesis is now acknowledged to occur in a highly bursty manner.

Results

We propose to view the inference problem as a fitting procedure for a mechanistic gene network model that is inherently stochastic and takes not only protein, but also mRNA levels into account. We first explain how to build and simulate this network model based upon the coupling of genes that are described as piecewise-deterministic Markov processes. Our model is modular and can be used to implement various biochemical hypotheses including causal interactions between genes. However, a naive fitting procedure would be intractable. By performing a relevant approximation of the stationary distribution, we derive a tractable procedure that corresponds to a statistical hidden Markov model with interpretable parameters. This approximation turns out to be extremely close to the theoretical distribution in the case of a simple toggle-switch, and we show that it can indeed fit real single-cell data. As a first step toward inference, our approach was applied to a number of simple two-gene networks simulated in silico from the mechanistic model and satisfactorily recovered the original networks.

Conclusions

Our results demonstrate that functional interactions between genes can be inferred from the distribution of a mechanistic, dynamical stochastic model that is able to describe gene expression in individual cells. This approach seems promising in relation to the current explosion of single-cell expression data.

Electronic supplementary material

The online version of this article (doi:10.1186/s12918-017-0487-0) contains supplementary material, which is available to authorized users.

Effect of transcription factor resource sharing on gene expression noise

Gene expression is intrinsically a stochastic (noisy) process with important implications for cellular functions. Deciphering the underlying mechanisms of gene expression noise remains one of the key challenges of regulatory biology. Theoretical models of transcription often incorporate the kinetics of how transcription factors (TFs) interact with a single promoter to impact gene expression noise. However, inside single cells multiple identical gene copies as well as additional binding sites can compete for a limiting pool of TFs. Here we develop a simple kinetic model of transcription, which explicitly incorporates this interplay between TF copy number and its binding sites. We show that TF sharing enhances noise in mRNA distribution across an isogenic population of cells. Moreover, when a single gene copy shares it’s TFs with multiple competitor sites, the mRNA variance as a function of the mean remains unaltered by their presence. Hence, all the data for variance as a function of mean expression collapse onto a single master curve independent of the strength and number of competitor sites. However, this result does not hold true when the competition stems from multiple copies of the same gene. Therefore, although previous studies showed that the mean expression follows a universal master curve, our findings suggest that different scenarios of competition bear distinct signatures at the level of variance. Intriguingly, the introduction of competitor sites can transform a unimodal mRNA distribution into a multimodal distribution. These results demonstrate the impact of limited availability of TF resource on the regulation of noise in gene expression.

Genetically identical cells, even when they are exposed to the same environmental conditions, display incredible diversity. Gene expression noise is attributed to be a key source of this phenotypic diversity. Transcriptional dynamics is a dominant source of expression noise. Although scores of theoretical and experimental studies have explored how noise is regulated at the level of transcription, most of them focus on the gene specific, cis regulatory elements, such as the number of transcription factor (TF) binding sites, their binding strength, etc. However, how the global properties of transcription, such as the limited availability of TFs impact noise in gene expression remains rather elusive. Here we build a theoretical model that incorporates the effect of limiting TF pool on gene expression noise. We find that competition between genes for TFs leads to enhanced variability in mRNA copy number across an isogenic population. Moreover, for gene copies sharing TFs with other competitor sites, mRNA variance as a function of the mean shows distinct imprints for one gene copy and multiple gene copies respectively. This stands in sharp contrast to the universal behavior found in mean expression irrespective of the different scenarios of competition. An interesting feature of competition is that introduction of competitor sites can transform a unimodal mRNA distribution into a multimodal distribution, which could lead to phenotypic variability.