Exact distributions for stochastic gene expression models with bursting and feedback

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.

Exact burst-size distributions for gene-expression models with complex promoter structure

In prokaryotic and eukaryotic cells, most genes are transcribed in a bursty fashion on one hand and complex gene regulations may lead to complex promoter structure on the other hand. This raises an unsolved issue: how does promoter structure shape transcriptional bursting kinetics characterized by burst size and frequency? Here we analyze stochastic models of gene transcription, which consider complex regulatory mechanisms. Notably, we develop an efficient method to derive exact burst-size distributions. The analytical results show that if the promoter of a gene contains only one active state, the burst size indeed follows a geometric distribution, in agreement with the previous result derived under certain limiting conditions. However, if it contains a multitude of active states, the burst size in general obeys a non-geometric distribution, which is a linearly weighted sum of geometric distributions. This superposition principle reveals the essential feature of bursting kinetics in complex cases of transcriptional regulation although it seems that there has been no direct experimental confirmation. The derived burst-size distributions not only highlight the importance of promoter structure in regulating bursting kinetics, but can be also used in the exact inference of this kinetics based on experimental data.

Solving stochastic gene-expression models using queueing theory: A tutorial review

Stochastic models of gene expression are typically formulated using the chemical master equation, which can be solved exactly or approximately using a repertoire of analytical methods. Here, we provide a tutorial review of an alternative approach based on queueing theory that has rarely been used in the literature of gene expression. We discuss the interpretation of six types of infinite-server queues from the angle of stochastic single-cell biology and provide analytical expressions for the stationary and nonstationary distributions and/or moments of mRNA/protein numbers and bounds on the Fano factor. This approach may enable the solution of complex models that have hitherto evaded analytical solution.

What can we learn when fitting a simple telegraph model to a complex gene expression model?

In experiments, the distributions of mRNA or protein numbers in single cells are often fitted to the random telegraph model which includes synthesis and decay of mRNA or protein, and switching of the gene between active and inactive states. While commonly used, this model does not describe how fluctuations are influenced by crucial biological mechanisms such as feedback regulation, non-exponential gene inactivation durations, and multiple gene activation pathways. Here we investigate the dynamical properties of four relatively complex gene expression models by fitting their steady-state mRNA or protein number distributions to the simple telegraph model. We show that despite the underlying complex biological mechanisms, the telegraph model with three effective parameters can accurately capture the steady-state gene product distributions, as well as the conditional distributions in the active gene state, of the complex models. Some effective parameters are reliable and can reflect realistic dynamic behaviors of the complex models, while others may deviate significantly from their real values in the complex models. The effective parameters can also be applied to characterize the capability for a complex model to exhibit multimodality. Using additional information such as single-cell data at multiple time points, we provide an effective method of distinguishing the complex models from the telegraph model. Furthermore, using measurements under varying experimental conditions, we show that fitting the mRNA or protein number distributions to the telegraph model may even reveal the underlying gene regulation mechanisms of the complex models. The effectiveness of these methods is confirmed by analysis of single-cell data for E. coli and mammalian cells. All these results are robust with respect to cooperative transcriptional regulation and extrinsic noise. In particular, we find that faster relaxation speed to the steady state results in more precise parameter inference under large extrinsic noise.

Solving the time-dependent protein distributions for autoregulated bursty gene expression using spectral decomposition

In this study, we obtain an exact time-dependent solution of the chemical master equation (CME) of an extension of the two-state telegraph model describing bursty or non-bursty protein expression in the presence of positive or negative autoregulation. Using the method of spectral decomposition, we show that the eigenfunctions of the generating function solution of the CME are Heun functions, while the eigenvalues can be determined by solving a continued fraction equation. Our solution generalizes and corrects a previous time-dependent solution for the CME of a gene circuit describing non-bursty protein expression in the presence of negative autoregulation [Ramos et al., Phys. Rev. E 83, 062902 (2011)]. In particular, we clarify that the eigenvalues are generally not real as previously claimed. We also investigate the relationship between different types of dynamic behavior and the type of feedback, the protein burst size, and the gene switching rate.

Power-law behavior of transcriptional bursting regulated by enhancer-promoter communication

Revealing how transcriptional bursting kinetics are genomically encoded is challenging because genome structures are stochastic at the organization level and are suggestively linked to gene transcription. To address this challenge, we develop a generic theoretical framework that integrates chromatin dynamics, enhancer-promoter (E-P) communication, and gene-state switching to study transcriptional bursting. The theory predicts that power law can be a general rule to quantitatively describe bursting modulations by E-P spatial communication. Specifically, burst frequency and burst size are up-regulated by E-P communication strength, following power laws with positive exponents. Analysis of the scaling exponents further reveals that burst frequency is preferentially regulated. Bursting kinetics are down-regulated by E-P genomic distance with negative power-law exponents, and this negative modulation desensitizes at large distances. The mutual information between burst frequency (or burst size) and E-P spatial distance further reveals essential characteristics of the information transfer from E-P communication to transcriptional bursting kinetics. These findings, which are in agreement with experimental observations, not only reveal fundamental principles of E-P communication in transcriptional bursting but also are essential for understanding cellular decision-making.

Exact results for gene-expression models with general waiting-time distributions

Complex molecular details of transcriptional regulation can be coarse-grained by assuming that reaction waiting times for promoter-state transitions, the mRNA synthesis, and the mRNA degradation follow general distributions. However, how such a generalized two-state model is analytically solved is a long-standing issue. Here we first present analytical formulas of burst-size distributions for this model. Then, we derive an iterative equation for the mRNA moment-generating function, by which mRNA raw and binomial moments of any order can be conveniently calculated. The analytical results obtained in the special cases of phase-type waiting-time distributions not only provide insights into the mechanisms of complex transcriptional regulations but also bring conveniences for experimental data-based statistical inferences.

4D nucleome equation predicts gene expression controlled by long-range enhancer-promoter interaction

Recent experimental evidence strongly supports that three-dimensional (3D) long-range enhancer-promoter (E-P) interactions have important influences on gene-expression dynamics, but it is unclear how the interaction information is translated into gene expression over time (4D). To address this question, we developed a general theoretical framework (named as a 4D nucleome equation), which integrates E-P interactions on chromatin and biochemical reactions of gene transcription. With this equation, we first present the distribution of mRNA counts as a function of the E-P genomic distance and then reveal a power-law scaling of the expression level in this distance. Interestingly, we find that long-range E-P interactions can induce bimodal and trimodal mRNA distributions. The 4D nucleome equation also allows for model selection and parameter inference. When this equation is applied to the mouse embryonic stem cell smRNA-FISH data and the E-P genomic-distance data, the predicted E-P contact probability and mRNA distribution are in good agreement with experimental results. Further statistical inference indicates that the E-P interactions prefer to modulate the mRNA level by controlling promoter activation and transcription initiation rates. Our model and results provide quantitative insights into both spatiotemporal gene-expression determinants (i.e., long-range E-P interactions) and cellular fates during development.

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.

Uncovering the effect of RNA polymerase steric interactions on gene expression noise: Analytical distributions of nascent and mature RNA numbers

The telegraph model is the standard model of stochastic gene expression, which can be solved exactly to obtain the distribution of mature RNA numbers per cell. A modification of this model also leads to an analytical distribution of nascent RNA numbers. These solutions are routinely used for the analysis of single-cell data, including the inference of transcriptional parameters. However, these models neglect important mechanistic features of transcription elongation, such as the stochastic movement of RNA polymerases and their steric (excluded-volume) interactions. Here we construct a model of gene expression describing promoter switching between inactive and active states, binding of RNA polymerases in the active state, their stochastic movement including steric interactions along the gene, and their unbinding leading to a mature transcript that subsequently decays. We derive the steady-state distributions of the nascent and mature RNA numbers in two important limiting cases: constitutive expression and slow promoter switching. We show that RNA fluctuations are suppressed by steric interactions between RNA polymerases, and that this suppression can in some instances even lead to sub-Poissonian fluctuations; these effects are most pronounced for nascent RNA and less prominent for mature RNA, since the latter is not a direct sensor of transcription. We find a relationship between the parameters of our microscopic mechanistic model and those of the standard models that ensures excellent consistency in their prediction of the first and second RNA number moments over vast regions of parameter space, encompassing slow, intermediate, and rapid promoter switching, provided the RNA number distributions are Poissonian or super-Poissonian. Furthermore, we identify the limitations of inference from mature RNA data, specifically showing that it cannot differentiate between highly distinct RNA polymerase traffic patterns on a gene.

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.

Theoretical investigation of functional responses of bio-molecular assembly networks

Cooperative protein-protein and protein-DNA interactions form programmable complex assemblies, often performing non-linear gene regulatory operations involved in signal transductions and cell fate determination. The apparent structure of those complex assemblies is very similar, but their functional response strongly depends on the topology of the protein-DNA interaction networks. Here, we demonstrate how the coordinated self-assembly creates gene regulatory network motifs that corroborate the existence of a precise functional response at the molecular level using thermodynamic and dynamic analyses. Our theoretical and Monte Carlo simulations show that a complex network of interactions can form a decision-making loop, such as feedback and feed-forward circuits, only by a few molecular mechanisms. We characterize each possible network of interactions by systematic variations of free energy parameters associated with the binding among biomolecules and DNA looping. We also find that the higher-order networks exhibit alternative steady states from the stochastic dynamics of each network. We capture this signature by calculating stochastic potentials and attributing their multi-stability features. We validate our findings against the Gal promoter system in yeast cells. Overall, we show that the network topology is vital in phenotype diversity in regulatory circuits.

A Markovian arrival stream approach to stochastic gene expression in cells

We analyse a generalisation of the stochastic gene expression model studied recently in Fromion et al. (SIAM J Appl Math 73:195-211, 2013) and Robert (Probab Surv 16:277-332, 2019) that keeps track of the production of both mRNA and protein molecules, using techniques from the theory of point processes, as well as ideas from the theory of matrix-analytic methods. Here, both the activity of a gene and the creation of mRNA are modelled with an arbitrary Markovian Arrival Process governed by finitely many phases, and each mRNA molecule during its lifetime gives rise to protein molecules in accordance with a Poisson process. This modification is important, as Markovian Arrival Processes can be used to approximate many types of point processes on the nonnegative real line, meaning this framework allows us to further relax our assumptions on the overall process of transcription.

The effect of microRNA on protein variability and gene expression fidelity

Small regulatory RNA molecules such as microRNA modulate gene expression through inhibiting the translation of messenger RNA (mRNA). Such posttranscriptional regulation has been recently hypothesized to reduce the stochastic variability of gene expression around average levels. Here, we quantify noise in stochastic gene expression models with and without such regulation. Our results suggest that silencing mRNA posttranscriptionally will always increase, rather than decrease, gene expression noise when the silencing of mRNA also increases its degradation, as is expected for microRNA interactions with mRNA. In that regime, we also find that silencing mRNA generally reduces the fidelity of signal transmission from deterministically varying upstream factors to protein levels. These findings suggest that microRNA binding to mRNA does not generically confer precision to protein expression.

Concentration fluctuations in growing and dividing cells: Insights into the emergence of concentration homeostasis

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.

Buffering variability in cell regulation motifs close to criticality

Bistable biological regulatory systems need to cope with stochastic noise to fine tune their function close to bifurcation points. Here, we study stability properties of this regime in generic systems to demonstrate that cooperative interactions buffer system variability, hampering noise-induced regime shifts. Our analysis also shows that, in the considered cooperativity range, impending regime shifts can be generically detected by statistical early warning signals from distributional data. Our generic framework, based on minimal models, can be used to extract robustness and variability properties of more complex models and empirical data close to criticality.

Stochastic gene transcription with non-competitive transcription regulatory architecture

The transcription factors, such as activators and repressors, can interact with the promoter of gene either in a competitive or non-competitive way. In this paper, we construct a stochastic model with non-competitive transcriptional regulatory architecture and develop an analytical theory that re-establishes the experimental results with an improved data fitting. The analytical expressions in the theory allow us to study the nature of the system corresponding to any of its parameters and hence, enable us to find out the factors that govern the regulation of gene expression for that architecture. We notice that, along with transcriptional reinitiation and repressors, there are other parameters that can control the noisiness of this network. We also observe that, the Fano factor (at mRNA level) varies from sub-Poissonian regime to super-Poissonian regime. In addition to the aforementioned properties, we observe some anomalous characteristics of the Fano factor (at mRNA level) and that of the variance of protein at lower activator concentrations in the presence of repressor molecules. This model is useful to understand the architecture of interactions which may buffer the stochasticity inherent to gene transcription.

Coloured noise induces phenotypic diversity with energy dissipation

Genes integrate many different sources of noise to adapt their survival strategy with energy costs, but how this noise impacts gene phenotype switching is not fully understood. Here, we refine a mechanistic model with multiplicative and additive coloured noise and analyse the influence of noise strength (NS) and autocorrelation time (AT) on gene phenotypic diversity. Different from white noise, we found that in the autocorrelation time-scale plane, increasing the multiplicative noise will broaden the bimodal region of the gene product, and additive noise will induce bimodal region drift from the lower level to the higher level, while the AT will promote this transition. Specifically, the effect of AT on gene expression is similar to a feedback loop; that is, the AT of multiplicative noise will elongate the mean first passage time (MFPT) from the low stable state to the high stable state, but it will reduce the MFPT from the high stable state to the low stable state, and the opposite is true for additive noise. Moreover, these transitions will violate the detailed equilibrium and then consume energy. By effective topology network reconstruction, we found that when the NS is small, the more obvious the bimodality is, the lower the energy dissipation; however, when the NS is large, it will consume more energy with a tendency for bimodality. The overall analysis implies that living organisms will utilize noise strength and its autocorrelation time for better survival in complex and fluctuating environments.

Exploring the role of eRNA in regulating gene expression

eRNAs as the products of enhancers can regulate gene expression via various possible ways, but which regulation way is more reasonable is debatable in biology, and in particular, how eRNAs impact gene expression remains unclear. Here we introduce a mechanistic model of gene expression to address these issues. This model considers three possible regulation ways of eRNA: Type-I by which eRNA regulates transcriptional activity by facilitating the formation of enhancer-promoter (E-P) loop, Type-II by which eRNA directly promotes the mRNA production rate, and mixed regulation (i.e., the combination of Type-I and Type-II). We show that with the increase of the E-P loop length, mRNA distribution can transition from unimodality to bimodality or vice versa in all the three regulation cases. However, in contrast to the other two regulations, Type-II regulation can lead to the highest mean mRNA level and the lowest mRNA noise, independent of the E-P loop length. These results would not only reveal the essential mechanism of how eRNA regulates gene expression, but also imply a new mechanism for phenotypic switching, namely the E-P loop can induce phenotypic switching.

Exact distributions for stochastic gene expression models with arbitrary promoter architecture and translational bursting

Gene expression in individual cells is inherently variable and sporadic, leading to cell-to-cell variability in mRNA and protein levels. Recent single-cell and single-molecule experiments indicate that promoter architecture and translational bursting play significant roles in controlling gene expression noise and generating the phenotypic diversity that life exhibits. To quantitatively understand the impact of these factors, it is essential to construct an accurate mathematical description of stochastic gene expression and find the exact analytical results, which is a formidable task. Here, we develop a stochastic model of bursty gene expression, which considers the complex promoter architecture governing the variability in mRNA expression and a general distribution characterizing translational burst. We derive the analytical expression for the corresponding protein steady-state distribution and all moment statistics of protein counts. We show that the total protein noise can be decomposed into three parts: the low-copy noise of protein due to probabilistic individual birth and death events, the noise due to stochastic switching between promoter states, and the noise resulting from translational busting. The theoretical results derived provide quantitative insights into the biochemical mechanisms of stochastic gene expression.

Volumetric compression develops noise-driven single-cell heterogeneity

Recent studies have revealed that extensive heterogeneity of biological systems arises through various routes ranging from intracellular chromosome segregation to spatiotemporally varying biochemical stimulations. However, the contribution of physical microenvironments to single-cell heterogeneity remains largely unexplored. Here, we show that a homogeneous population of non-small-cell lung carcinoma develops into heterogeneous subpopulations upon application of a homogeneous physical compression, as shown by single-cell transcriptome profiling. The generated subpopulations stochastically gain the signature genes associated with epithelial-mesenchymal transition (EMT; VIM, CDH1, EPCAM, ZEB1, and ZEB2) and cancer stem cells (MKI67, BIRC5, and KLF4), respectively. Trajectory analysis revealed two bifurcated paths as cells evolving upon the physical compression, along each path the corresponding signature genes (epithelial or mesenchymal) gradually increase. Furthermore, we show that compression increases gene expression noise, which interplays with regulatory network architecture and thus generates differential cell-fate outcomes. The experimental observations of both single-cell sequencing and single-molecule fluorescent in situ hybridization agrees well with our computational modeling of regulatory network in the EMT process. These results demonstrate a paradigm of how mechanical stimulations impact cell-fate determination by altering transcription dynamics; moreover, we show a distinct path that the ecology and evolution of cancer interplay with their physical microenvironments from the view of mechanobiology and systems biology, with insight into the origin of single-cell heterogeneity.

Detailed balance, local detailed balance, and global potential for stochastic chemical reaction networks

Detailed balance of a chemical reaction network can be defined in several different ways. Here we investigate the relationship among four types of detailed balance conditions: deterministic, stochastic, local, and zero-order local detailed balance. We show that the four types of detailed balance are equivalent when different reactions lead to different species changes and are not equivalent when some different reactions lead to the same species change. Under the condition of local detailed balance, we further show that the system has a global potential defined over the whole space, which plays a central role in the large deviation theory and the Freidlin–Wentzell-type metastability theory of chemical reaction networks. Finally, we provide a new sufficient condition for stochastic detailed balance, which is applied to construct a class of high-dimensional chemical reaction networks that both satisfies stochastic detailed balance and displays multistability.

Understanding the molecular mechanisms of transcriptional bursting

In recent years, it has been experimentally established that transcription, a fundamental biological process that involves the synthesis of messenger RNA molecules from DNA templates, does not proceed continuously as was expected. Rather, it exhibits a distinct dynamic behavior of alternating between productive phases when RNA molecules are actively synthesized and inactive phases when there is no RNA production at all. The bimodal transcriptional dynamics is now confirmed to be present in most living systems. This phenomenon is known as transcriptional bursting and it attracts significant amounts of attention from researchers in different fields. However, despite multiple experimental and theoretical investigations, the microscopic origin and biological functions of the transcriptional bursting remain unclear. Here we discuss the recent developments in uncovering the underlying molecular mechanisms of transcriptional bursting and our current understanding of them. Our analysis presents a physicochemical view of the processes that govern transcriptional bursting in living cells.

Steady-state joint distribution for first-order stochastic reaction kinetics

While the analytical solution for the marginal distribution of a stochastic chemical reaction network has been extensively studied, its joint distribution, i.e., the solution of a high-dimensional chemical master equation, has received much less attention. Here we develop an alternative method of computing the exact joint distributions of a wide class of first-order stochastic reaction systems in steady-state conditions. The effectiveness of our method is validated by applying it to four gene expression models of biological significance, including models with 2A peptides, nascent mRNA, gene regulation, translational bursting, and alternative splicing.

Noise suppression in stochastic genetic circuits using PID controllers

Inside individual cells, protein population counts are subject to molecular noise due to low copy numbers and the inherent probabilistic nature of biochemical processes. We investigate the effectiveness of proportional, integral and derivative (PID) based feedback controllers to suppress protein count fluctuations originating from two noise sources: bursty expression of the protein, and external disturbance in protein synthesis. Designs of biochemical reactions that function as PID controllers are discussed, with particular focus on individual controllers separately, and the corresponding closed-loop system is analyzed for stochastic controller realizations. Our results show that proportional controllers are effective in buffering protein copy number fluctuations from both noise sources, but this noise suppression comes at the cost of reduced static sensitivity of the output to the input signal. In contrast, integral feedback has no effect on the protein noise level from stochastic expression, but significantly minimizes the impact of external disturbances, particularly when the disturbance comes at low frequencies. Counter-intuitively, integral feedback is found to amplify external disturbances at intermediate frequencies. Next, we discuss the design of a coupled feedforward-feedback biochemical circuit that approximately functions as a derivate controller. Analysis using both analytical methods and Monte Carlo simulations reveals that this derivative controller effectively buffers output fluctuations from bursty stochastic expression, while maintaining the static input-output sensitivity of the open-loop system. In summary, this study provides a systematic stochastic analysis of biochemical controllers, and paves the way for their synthetic design and implementation to minimize deleterious fluctuations in gene product levels.

In the noisy cellular environment, biochemical species such as genes, RNAs and proteins that often occur at low molecular counts, are subject to considerable stochastic fluctuations in copy numbers over time. How cellular biochemical processes function reliably in the face of such randomness is an intriguing fundamental problem. Increasing evidence suggests that random fluctuations (noise) in protein copy numbers play important functional roles, such as driving genetically identical cells to different cell fates. Moreover, many disease states have been attributed to elevated noise levels in specific proteins. Here we systematically investigate design of biochemical systems that function as proportional, integral and derivative-based feedback controllers to suppress protein count fluctuations arising from bursty expression of the protein and external disturbance in protein synthesis. Our results show that different controllers are effective in buffering different noise components, and identify ranges of feedback gain for minimizing deleterious fluctuations in protein levels.

Optimal transition paths of phenotypic switching in a non-Markovian self-regulation gene expression

Gene expression is a complex biochemical process involving multiple reaction steps, creating molecular memory because the probability of waiting time between consecutive reaction steps no longer follows exponential distributions. What effect the molecular memory has on metastable states in gene expression remains not fully understood. Here, we study transition paths of switching between bistable states for a non-Markovian model of gene expression equipped with a self-regulation. Employing the large deviation theory for this model, we analyze the optimal transition paths of switching between bistable states in gene expression, interestingly finding that dynamic behaviors in gene expression along the optimal transition paths significantly depend on the molecular memory. Moreover, we discover that the molecular memory can prolong the time of switching between bistable states in gene expression along the optimal transition paths. Our results imply that the molecular memory may be an unneglectable factor to affect switching between metastable states in gene expression.

Analytic solutions for stochastic hybrid models of gene regulatory networks

Discrete-state stochastic models are a popular approach to describe the inherent stochasticity of gene expression in single cells. The analysis of such models is hindered by the fact that the underlying discrete state space is extremely large. Therefore hybrid models, in which protein counts are replaced by average protein concentrations, have become a popular alternative. The evolution of the corresponding probability density functions is given by a coupled system of hyperbolic PDEs. This system has Markovian nature but its hyperbolic structure makes it difficult to apply standard functional analytical methods. We are able to prove convergence towards the stationary solution and determine such equilibrium explicitly by combining abstract methods from the theory of positive operators and elementary ideas from potential analysis.

Queuing Models of Gene Expression: Analytical Distributions and Beyond

Activation of a gene is a multistep biochemical process, involving recruitments of transcription factors and histone kinases as well as modification of histones. Many of these intermediate reaction steps would have been unspecified by experiments. Therefore, classical two-state models of gene expression established based on the memoryless (or Markovian) assumption would not well describe the reality in gene expression. Recent experimental data have indicated that the inactive phases of gene promoters are differently distributed, showing strong memory. Here, we use a nonexponential waiting-time distribution to model the complex activation process of a gene, and then analyze a queuing model of stochastic transcription. We successfully derive the analytical expression of the stationary mRNA distribution, which provides insight into the effect of molecular memory created by complex activating events on the mRNA expression. We find that the reduction in the waiting-time noise may result in the increase in the mRNA noise, contrary to the previous conclusion. Based on the derived distribution, we also develop a method to infer the waiting-time distribution from a known mRNA distribution. Data analysis on a realistic example verifies the validity of this method.

Urn models for stochastic gene expression yield intuitive insights into the probability distributions of single-cell mRNA and protein counts

Fitting the probability mass functions from analytical solutions of stochastic models of gene expression to the single-cell count distributions of mRNA and protein molecules can yield valuable insights into mechanisms underlying gene expression. Solutions of chemical master equations are available for various kinetic schemes but, even for the basic ON-OFF genetic switch, they take complex forms with generating functions given as hypergeometric functions. Interpretation of gene expression dynamics in terms of bursts is not consistent with the complete range of parameters for these functions. Physical insights into the probability mass functions are essential to ensure proper interpretations but are lacking for models considering genetic switches. To fill this gap, we develop urn models for stochastic gene expression. We sample RNA polymerases or ribosomes from a master urn, which represents the cytosol, and assign them to recipient urns of two or more colors, which represent time intervals in which no switching occurs. Colors of the recipient urns represent sub-systems of the promoter states, and the assignments to urns of a specific color represent gene expression. We use elementary principles of discrete probability theory to solve a range of kinetic models without feedback, including the Peccoud-Ycart model, the Shahrezaei-Swain model, and models with an arbitrary number of promoter states. In the last case, we obtain a novel result for the protein distribution. For activated genes, we show that transcriptional lapses, which are events of gene inactivation for short time intervals separated by long active intervals, quantify the transcriptional dynamics better than bursts. We show that the intuition gained from our urn models may also be useful in understanding existing solutions for models with feedback. We contrast our models with urn models for related distributions, discuss a generalization of the Delaporte distribution for single-cell data analysis, and highlight the limitations of our models.

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.

Analytical results for non-Markovian models of bursty gene expression

Modeling stochastic gene expression has long relied on Markovian hypothesis. In recent years, however, this hypothesis is challenged by the increasing availability of time-resolved data. Correspondingly, there is considerable interest in understanding how non-Markovian reaction kinetics of gene expression impact protein variations across a population of genetically identical cells. Here, we analyze a stochastic model of gene expression with arbitrary waiting-time distributions, which includes existing gene models as its special cases. We find that stationary probabilistic behavior of this non-Markovian system is exactly the same as that of an equivalent Markovian system with the same substrates. Based on this fact, we derive analytical results, which provide insight into the roles of feedback regulation and molecular memory in controlling the protein noise and properties of the steady states, which are inaccessible via existing methodology. Our results also provide quantitative insight into diverse cellular processes involving stochastic sources of gene expression and molecular memory.

Stochastic Modeling of Autoregulatory Genetic Feedback Loops: A Review and Comparative Study

Autoregulatory 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 1) which subcellular processes are explicitly modeled, 2) the modeling methodology employed (discrete, continuous, or hybrid), and 3) 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 modeling.

Ab initio scaling laws between noise and mean of gene expression

Gene expression is inherently noisy due to fluctuations occurring at the molecular level. From a top-down perspective, noise has been traditionally decomposed into an intrinsic component that scales inversely with the mean expression level and an extrinsic component that is constant in absence of regulatory changes. Here, we adopt a bottom-up approach to reveal that extrinsic noise, by itself, can follow the aforementioned decomposition, which entails that one component of it can be confounded with intrinsic noise. Analytical expressions of the noise-mean relationship were derived for different scenarios, which were in part supported by numerical simulations.

Different effects of fast and slow input fluctuations on output in gene regulation

An important task in the post-gene era is to understand the role of stochasticity in gene regulation. Here, we analyze a cascade model of stochastic gene expression, where the upstream gene stochastically generates proteins that regulate, as transcription factors, stochastic synthesis of the downstream output. We find that in contrast to fast input fluctuations that do not change the behavior of the downstream system qualitatively, slow input fluctuations can induce different modes of the distribution of downstream output and even stochastic focusing or defocusing of the downstream output level, although the regulatory protein follows the same distribution in both cases. This finding is counterintuitive but can have broad biological implications, e.g., slow input rather than fast fluctuations may both increase the survival probability of cells and enhance the sensitivity of intracellular regulation. In addition, we find that input fluctuations can minimize the output noise.

Constraining the complexity of promoter dynamics using fluctuations in gene expression

Gene expression is an inherently stochastic process with transcription of mRNAs often occurring in bursts: short periods of activity followed by typically longer periods of inactivity. While a simple model involving switching between two promoter states has been widely used to analyze transcription dynamics, recent experimental observations have provided evidence for more complex kinetic schemes underlying bursting. Specifically, experiments provide evidence for complexity in promoter dynamics during the switch from the transcriptionally inactive to the transcriptionally active state. An open question in the field is: what is the minimal complexity needed to model promoter dynamics and how can we determine this? Here, we show that measurements of mRNA fluctuations can be used to set fundamental bounds on the complexity of promoter dynamics. We study models wherein the switching time distribution from transcriptionally inactive to active states is described by a general waiting-time distribution. Using approaches from renewal theory and queueing theory, we derive analytical expressions which connect the Fano factor of mRNA distributions to the waiting-time distribution for promoter switching between inactive and active states. The results derived lead to bounds on the minimal number of promoter states and thus allow us to derive bounds on the minimal complexity of promoter dynamics based on single-cell measurements of mRNA levels.

Single-cell stochastic gene expression kinetics with coupled positive-plus-negative feedback

Here we investigate single-cell stochastic gene expression kinetics in a minimal coupled gene circuit with positive-plus-negative feedback. A triphasic stochastic bifurcation is observed upon increasing the ratio of the positive and negative feedback strengths, which reveals a strong synergistic interaction between positive and negative feedback loops. We discover that coupled positive-plus-negative feedback amplifies gene expression mean but reduces gene expression noise over a wide range of feedback strengths when promoter switching is relatively slow, stabilizing gene expression around a relatively high level. In addition, we study two types of macroscopic limits of the discrete chemical master equation model: the Kurtz limit applies to proteins with large burst frequencies and the Lévy limit applies to proteins with large burst sizes. We derive the analytic steady-state distributions of the protein abundance in a coupled gene circuit for both the discrete model and its two macroscopic limits, generalizing the results obtained by Liu et al. [Chaos 26, 043108 (2016)CHAOEH1054-150010.1063/1.4947202]. We also obtain the analytic time-dependent protein distribution for the classical Friedman-Cai-Xie random bursting model [Friedman, Cai, and Xie, Phys. Rev. Lett. 97, 168302 (2006)PRLTAO0031-900710.1103/PhysRevLett.97.168302]. Our analytic results are further applied to study the structure of gene expression noise in a coupled gene circuit, and a complete decomposition of noise in terms of five different biophysical origins is provided.

Decoding the grammar of transcriptional regulation from RNA polymerase measurements: models and their applications

The genomic revolution has indubitably brought about a paradigm shift in the field of molecular biology, wherein we can sequence, write and re-write genomes. In spite of achieving such feats, we still lack a quantitative understanding of how cells integrate environmental and intra-cellular signals at the promoter and accordingly regulate the production of messenger RNAs. This current state of affairs is being redressed by recent experimental breakthroughs which enable the counting of RNA polymerase molecules (or the corresponding nascent RNAs) engaged in the process of transcribing a gene at the single-cell level. Theorists, in conjunction, have sought to unravel the grammar of transcriptional regulation by harnessing the various statistical properties of these measurements. In this review, we focus on the recent progress in developing falsifiable models of transcription that aim to connect the molecular mechanisms of transcription to single-cell polymerase measurements. We discuss studies where the application of such models to the experimental data have led to novel mechanistic insights into the process of transcriptional regulation. Such interplay between theory and experiments will likely contribute towards the exciting journey of unfurling the governing principles of transcriptional regulation ranging from bacteria to higher organisms.

Revisiting the Reduction of Stochastic Models of Genetic Feedback Loops with Fast Promoter Switching

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.

Exact and efficient hybrid Monte Carlo algorithm for accelerated Bayesian inference of gene expression models from snapshots of single-cell transcripts

Single cells exhibit a significant amount of variability in transcript levels, which arises from slow, stochastic transitions between gene expression states. Elucidating the nature of these states and understanding how transition rates are affected by different regulatory mechanisms require state-of-the-art methods to infer underlying models of gene expression from single cell data. A Bayesian approach to statistical inference is the most suitable method for model selection and uncertainty quantification of kinetic parameters using small data sets. However, this approach is impractical because current algorithms are too slow to handle typical models of gene expression. To solve this problem, we first show that time-dependent mRNA distributions of discrete-state models of gene expression are dynamic Poisson mixtures, whose mixing kernels are characterized by a piecewise deterministic Markov process. We combined this analytical result with a kinetic Monte Carlo algorithm to create a hybrid numerical method that accelerates the calculation of time-dependent mRNA distributions by 1000-fold compared to current methods. We then integrated the hybrid algorithm into an existing Monte Carlo sampler to estimate the Bayesian posterior distribution of many different, competing models in a reasonable amount of time. We demonstrate that kinetic parameters can be reasonably constrained for modestly sampled data sets if the model is known a priori. If there are many competing models, Bayesian evidence can rigorously quantify the likelihood of a model relative to other models from the data. We demonstrate that Bayesian evidence selects the true model and outperforms approximate metrics typically used for model selection.

Bounding the stationary distributions of the chemical master equation via mathematical programming

The stochastic dynamics of biochemical networks are usually modeled with the chemical master equation (CME). The stationary distributions of CMEs are seldom solvable analytically, and numerical methods typically produce estimates with uncontrolled errors. Here, we introduce mathematical programming approaches that yield approximations of these distributions with computable error bounds which enable the verification of their accuracy. First, we use semidefinite programming to compute increasingly tighter upper and lower bounds on the moments of the stationary distributions for networks with rational propensities. Second, we use these moment bounds to formulate linear programs that yield convergent upper and lower bounds on the stationary distributions themselves, their marginals, and stationary averages. The bounds obtained also provide a computational test for the uniqueness of the distribution. In the unique case, the bounds form an approximation of the stationary distribution with a computable bound on its error. In the nonunique case, our approach yields converging approximations of the ergodic distributions. We illustrate our methodology through several biochemical examples taken from the literature: Schlögl's model for a chemical bifurcation, a two-dimensional toggle switch, a model for bursty gene expression, and a dimerization model with multiple stationary distributions.

Stochastic Modeling of Gene Regulation by Noncoding Small RNAs in the Strong Interaction Limit

Noncoding small RNAs (sRNAs) are known to play a key role in regulating diverse cellular processes, and their dysregulation is linked to various diseases such as cancer. Such diseases are also marked by phenotypic heterogeneity, which is often driven by the intrinsic stochasticity of gene expression. Correspondingly, there is significant interest in developing quantitative models focusing on the interplay between stochastic gene expression and regulation by sRNAs. We consider the canonical model of regulation of stochastic gene expression by sRNAs, wherein interaction between constitutively expressed sRNAs and mRNAs leads to stoichiometric mutual degradation. The exact solution of this model is analytically intractable given the nonlinear interaction term between sRNAs and mRNAs, and theoretical approaches typically invoke the mean-field approximation. However, mean-field results are inaccurate in the limit of strong interactions and low abundances; thus, alternative theoretical approaches are needed. In this work, we obtain analytical results for the canonical model of regulation of stochastic gene expression by sRNAs in the strong interaction limit. We derive analytical results for the steady-state generating function of the joint distribution of mRNAs and sRNAs in the limit of strong interactions and use the results derived to obtain analytical expressions characterizing the corresponding protein steady-state distribution. The results obtained can serve as building blocks for the analysis of genetic circuits involving sRNAs and provide new insights into the role of sRNAs in regulating stochastic gene expression in the limit of strong interactions.

Extending the linear-noise approximation to biochemical systems influenced by intrinsic noise and slow lognormally distributed extrinsic noise

It is well known that the kinetics of an intracellular biochemical network is stochastic. This is due to intrinsic noise arising from the random timing of biochemical reactions in the network as well as due to extrinsic noise stemming from the interaction of unknown molecular components with the network and from the cell's changing environment. While there are many methods to study the effect of intrinsic noise on the system dynamics, few exist to study the influence of both types of noise. Here we show how one can extend the conventional linear-noise approximation to allow for the rapid evaluation of the molecule numbers statistics of a biochemical network influenced by intrinsic noise and by slow lognormally distributed extrinsic noise. The theory is applied to simple models of gene regulatory networks and its validity confirmed by comparison with exact stochastic simulations. In particular, we consider three important biological examples. First, we investigate how extrinsic noise modifies the dependence of the variance of the molecule number fluctuations on the rate constants. Second, we show how the mutual information between input and output of a network motif is affected by extrinsic noise. And third, we study the robustness of the ubiquitously found feed-forward loop motifs when subjected to extrinsic noise.

A stochastic model for post-transcriptional regulation of rare events in gene expression

Post-transcriptional regulation of gene expression is often a critical component of cellular processes involved in cell-fate decisions. Correspondingly, considerable efforts have focused on modeling post-transcriptional regulation of stochastic gene expression and on quantifying its impact on the mean and variance of protein levels. However, the impact of post-transcriptional regulation on rare events corresponding to large deviations in gene expression is less well understood. Here, we study a simple model involving post-transcriptional control of gene expression and characterize the impact of regulation on large deviations in activity fluctuations. We derive analytical results for the large deviation function for protein production rate and for the corresponding driven process which characterizes system fluctuations conditional on the rare event. Our results suggest that fluctuations in burst size rather than frequency play a dominant role in controlling rare events. The results derived also provide insight into how post-transcriptional regulation can be used to fine-tune the probability of rare fluctuations and to thereby control phenotypic variation in a population of cells.

Accuracy of parameter estimation for auto-regulatory transcriptional feedback loops from noisy data

Bayesian and non-Bayesian moment-based inference methods are commonly used to estimate the parameters defining stochastic models of gene regulatory networks from noisy single cell or population snapshot data. However, a systematic investigation of the accuracy of the predictions of these methods remains missing. Here, we present the results of such a study using synthetic noisy data of a negative auto-regulatory transcriptional feedback loop, one of the most common building blocks of complex gene regulatory networks. We study the error in parameter estimation as a function of (i) number of cells in each sample; (ii) the number of time points; (iii) the highest-order moment of protein fluctuations used for inference; (iv) the moment-closure method used for likelihood approximation. We find that for sample sizes typical of flow cytometry experiments, parameter estimation by maximizing the likelihood is as accurate as using Bayesian methods but with a much reduced computational time. We also show that the choice of moment-closure method is the crucial factor determining the maximum achievable accuracy of moment-based inference methods. Common likelihood approximation methods based on the linear noise approximation or the zero cumulants closure perform poorly for feedback loops with large protein-DNA binding rates or large protein bursts; this is exacerbated for highly heterogeneous cell populations. By contrast, approximating the likelihood using the linear-mapping approximation or conditional derivative matching leads to highly accurate parameter estimates for a wide range of conditions.

Multi-modality in gene regulatory networks with slow promoter kinetics

Phenotypical variability in the absence of genetic variation often reflects complex energetic landscapes associated with underlying gene regulatory networks (GRNs). In this view, different phenotypes are associated with alternative states of complex nonlinear systems: stable attractors in deterministic models or modes of stationary distributions in stochastic descriptions. We provide theoretical and practical characterizations of these landscapes, specifically focusing on stochastic Slow Promoter Kinetics (SPK), a time scale relevant when transcription factor binding and unbinding are affected by epigenetic processes like DNA methylation and chromatin remodeling. In this case, largely unexplored except for numerical simulations, adiabatic approximations of promoter kinetics are not appropriate. In contrast to the existing literature, we provide rigorous analytic characterizations of multiple modes. A general formal approach gives insight into the influence of parameters and the prediction of how changes in GRN wiring, for example through mutations or artificial interventions, impact the possible number, location, and likelihood of alternative states. We adapt tools from the mathematical field of singular perturbation theory to represent stationary distributions of Chemical Master Equations for GRNs as mixtures of Poisson distributions and obtain explicit formulas for the locations and probabilities of metastable states as a function of the parameters describing the system. As illustrations, the theory is used to tease out the role of cooperative binding in stochastic models in comparison to deterministic models, and applications are given to various model systems, such as toggle switches in isolation or in communicating populations, a synthetic oscillator, and a trans-differentiation network.