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

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.

Effects of mRNA conformational switching on translational noise in gene circuits

Intragenic translational heterogeneity describes the variation in translation at the level of transcripts for an individual gene. A factor that contributes to this source of variation is the mRNA structure. Both the composition of the thermodynamic ensemble, i.e., the stationary distribution of mRNA structures, and the switching dynamics between those play a role. The effect of the switching dynamics on intragenic translational heterogeneity remains poorly understood. We present a stochastic translation model that accounts for mRNA structure switching and is derived from a Markov model via approximate stochastic filtering. We assess the approximation on various timescales and provide a method to quantify how mRNA structure dynamics contributes to translational heterogeneity. With our approach, we allow quantitative information on mRNA switching from biophysical experiments or coarse-grain molecular dynamics simulations of mRNA structures to be included in gene regulatory chemical reaction network models without an increase in the number of species. Thereby, our model bridges a gap between mRNA structure kinetics and gene expression models, which we hope will further improve our understanding of gene regulatory networks and facilitate genetic circuit design.

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.

Density physics-informed neural networks reveal sources of cell heterogeneity in signal transduction

The transduction time between signal initiation and final response provides valuable information on the underlying signaling pathway, including its speed and precision. Furthermore, multi-modality in a transduction-time distribution indicates that the response is regulated by multiple pathways with different transduction speeds. Here, we developed a method called density physics-informed neural networks (Density-PINNs) to infer the transduction-time distribution from measurable final stress response time traces. We applied Density-PINNs to single-cell gene expression data from sixteen promoters regulated by unknown pathways in response to antibiotic stresses. We found that promoters with slower signaling initiation and transduction exhibit larger cell-to-cell heterogeneity in response intensity. However, this heterogeneity was greatly reduced when the response was regulated by slow and fast pathways together. This suggests a strategy for identifying effective signaling pathways for consistent cellular responses to disease treatments. Density-PINNs can also be applied to understand other time delay systems, including infectious diseases.

Epigenetic OCT4 regulatory network: stochastic analysis of cellular reprogramming

Experimental studies have shown that chromatin modifiers have a critical effect on cellular reprogramming, i.e., the conversion of differentiated cells to pluripotent stem cells. Here, we develop a model of the OCT4 gene regulatory network that includes genes expressing chromatin modifiers TET1 and JMJD2, and the chromatin modification circuit on which these modifiers act. We employ this model to compare three reprogramming approaches that have been considered in the literature with respect to reprogramming efficiency and latency variability. These approaches are overexpression of OCT4 alone, overexpression of OCT4 with TET1, and overexpression of OCT4 with JMJD2. Our results show more efficient and less variable reprogramming when also JMJD2 and TET1 are overexpressed, consistent with previous experimental data. Nevertheless, TET1 overexpression can lead to more efficient reprogramming compared to JMJD2 overexpression. This is the case when the recruitment of DNA methylation by H3K9me3 is weak and the methyl-CpG-binding domain (MBD) proteins are sufficiently scarce such that they do not hamper TET1 binding to methylated DNA. The model that we developed provides a mechanistic understanding of existing experimental results and is also a tool for designing optimized reprogramming approaches that combine overexpression of cell-fate specific transcription factors (TFs) with targeted recruitment of epigenetic modifiers.

Exact solution of a three-stage model of stochastic gene expression including cell-cycle dynamics

The classical three-stage model of stochastic gene expression predicts the statistics of single cell mRNA and protein number fluctuations as a function of the rates of promoter switching, transcription, translation, degradation and dilution. While this model is easily simulated, its analytical solution remains an unsolved problem. Here we modify this model to explicitly include cell-cycle dynamics and then derive an exact solution for the time-dependent joint distribution of mRNA and protein numbers. We show large differences between this model and the classical model which captures cell-cycle effects implicitly via effective first-order dilution reactions. In particular we find that the Fano factor of protein numbers calculated from a population snapshot measurement are underestimated by the classical model whereas the correlation between mRNA and protein can be either over- or underestimated, depending on the timescales of mRNA degradation and promoter switching relative to the mean cell-cycle duration time.

Inferring delays in partially observed gene regulation processes

MOTIVATION

Cell function is regulated by gene regulatory networks (GRNs) defined by protein-mediated interaction between constituent genes. Despite advances in experimental techniques, we can still measure only a fraction of the processes that govern GRN dynamics. To infer the properties of GRNs using partial observation, unobserved sequential processes can be replaced with distributed time delays, yielding non-Markovian models. Inference methods based on the resulting model suffer from the curse of dimensionality.

RESULTS

We develop a simulation-based Bayesian MCMC method employing an approximate likelihood for the efficient and accurate inference of GRN parameters when only some of their products are observed. We illustrate our approach using a two-step activation model: an activation signal leads to the accumulation of an unobserved regulatory protein, which triggers the expression of observed fluorescent proteins. With prior information about observed fluorescent protein synthesis, our method successfully infers the dynamics of the unobserved regulatory protein. We can estimate the delay and kinetic parameters characterizing target regulation including transcription, translation, and target searching of an unobserved protein from experimental measurements of the products of its target gene. Our method is scalable and can be used to analyze non-Markovian models with hidden components.

AVAILABILITY AND IMPLEMENTATION

Our code is implemented in R and is freely available with a simple example data at https://github.com/Mathbiomed/SimMCMC.

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.

Fluctuations in tissue growth portray homeostasis as a critical state and long-time non-Markovian cell proliferation as Markovian

Tissue growth is an emerging phenomenon that results from the cell-level interplay between proliferation and apoptosis, which is crucial during embryonic development, tissue regeneration, as well as in pathological conditions such as cancer. In this theoretical article, we address the problem of stochasticity in tissue growth by first considering a minimal Markovian model of tissue size, quantified as the number of cells in a simulated tissue, subjected to both proliferation and apoptosis. We find two dynamic phases, growth and decay, separated by a critical state representing a homeostatic tissue. Since the main limitation of the Markovian model is its neglect of the cell cycle, we incorporated a refractory period that temporarily prevents proliferation immediately following cell division, as a minimal proxy for the cell cycle, and studied the model in the growth phase. Importantly, we obtained from this last model an effective Markovian rate, which accurately describes general trends of tissue size. This study shows that the dynamics of tissue growth can be theoretically conceptualized as a Markovian process where homeostasis is a critical state flanked by decay and growth phases. Notably, in the growing non-Markovian model, a Markovian-like growth process emerges at large time scales.

Bye bye, linearity, bye: quantification of the mean for linear CRNs in a random environment

Molecular reactions within a cell are inherently stochastic, and cells often differ in morphological properties or interact with a heterogeneous environment. Consequently, cell populations exhibit heterogeneity both due to these intrinsic and extrinsic causes. Although state-of-the-art studies that focus on dissecting this heterogeneity use single-cell measurements, the bulk data that shows only the mean expression levels is still in routine use. The fingerprint of the heterogeneity is present also in bulk data, despite being hidden from direct measurement. In particular, this heterogeneity can affect the mean expression levels via bimolecular interactions with low-abundant environment species. We make this statement rigorous for the class of linear reaction systems that are embedded in a discrete state Markov environment. The analytic expression that we provide for the stationary mean depends on the reaction rate constants of the linear subsystem, as well as the generator and stationary distribution of the Markov environment. We demonstrate the effect of the environment on the stationary mean. Namely, we show how the heterogeneous case deviates from the quasi-steady state (Q.SS) case when the embedded system is fast compared to the environment.

The σB alternative sigma factor circuit modulates noise to generate different types of pulsing dynamics

Single-cell approaches are revealing a high degree of heterogeneity, or noise, in gene expression in isogenic bacteria. How gene circuits modulate this noise in gene expression to generate robust output dynamics is unclear. Here we use the Bacillus subtilis alternative sigma factor σB as a model system for understanding the role of noise in generating circuit output dynamics. σB controls the general stress response in B. subtilis and is activated by a range of energy and environmental stresses. Recent single-cell studies have revealed that the circuit can generate two distinct outputs, stochastic pulsing and a single pulse response, but the conditions under which each response is generated are under debate. We implement a stochastic mathematical model of the σB circuit to investigate this and find that the system's core circuit can generate both response types. This is despite one response (stochastic pulsing) being stochastic in nature, and the other (single response pulse) being deterministic. We demonstrate that the main determinant for whichever response is generated is the degree with which the input pathway activates the core circuit, although the noise properties of the input pathway also biases the system towards one or the other type of output. Thus, our work shows how stochastic modelling can reveal the mechanisms behind non-intuitive gene circuit output dynamics.

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.

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.

A dynamical stochastic model of yeast translation across the cell cycle

Translation is a central step in gene expression, however its quantitative and time-resolved regulation is poorly understood. We developed a discrete, stochastic model for protein translation in S. cerevisiae in a whole-transcriptome, single-cell context. A "base case" scenario representing an average cell highlights translation initiation rates as the main co-translational regulatory parameters. Codon usage bias emerges as a secondary regulatory mechanism through ribosome stalling. Demand for anticodons with low abundancy is shown to cause above-average ribosome dwelling times. Codon usage bias correlates strongly both with protein synthesis rates and elongation rates. Applying the model to a time-resolved transcriptome estimated by combining data from FISH and RNA-Seq experiments, it could be shown that increased total transcript abundance during the cell cycle decreases translation efficiency at single transcript level. Translation efficiency grouped by gene function shows highest values for ribosomal and glycolytic genes. Ribosomal proteins peak in S phase while glycolytic proteins rank highest in later cell cycle phases.

Genome-wide inference reveals that feedback regulations constrain promoter-dependent transcriptional burst kinetics

Gene expression in mammalian cells is highly variable and episodic, resulting in a series of discontinuous bursts of mRNAs. A challenge is to understand how static promoter architecture and dynamic feedback regulations dictate bursting on a genome-wide scale. Although single-cell RNA sequencing (scRNA-seq) provides an opportunity to address this challenge, effective analytical methods are scarce. We developed an interpretable and scalable inference framework, which combined experimental data with a mechanistic model to infer transcriptional burst kinetics (sizes and frequencies) and feedback regulations. Applying this framework to scRNA-seq data generated from embryonic mouse fibroblast cells, we found Simpson's paradoxes, i.e. genome-wide burst kinetics exhibit different characteristics in two cases without and with distinguishing feedback regulations. We also showed that feedbacks differently modulate burst frequencies and sizes and conceal the effects of transcription start site distributions on burst kinetics. Notably, only in the presence of positive feedback, TATA genes are expressed with high burst frequencies and enhancer-promoter interactions mainly modulate burst frequencies. The developed inference method provided a flexible and efficient way to investigate transcriptional burst kinetics and the obtained results would be helpful for understanding cell development and fate decision.

Interplay of Cellular mRNA, miRNA and Viral miRNA during Infection of a Cell

The understanding of the kinetics of gene expression in cells infected by viruses is currently limited. As a rule, the corresponding models do not take viral microRNAs (miRNAs) into account. Such RNAs are, however, operative during the replication of some viruses, including, e.g., herpesvirus. To clarify the kinetics of this category (with emphasis on the information available for herpesvirus), I introduce a generic model describing the transient interplay of cellular mRNA, protein, miRNA and viral miRNA. In the absence of viral miRNA, the cellular miRNA is considered to suppress the populations of mRNA and protein due to association with mRNA and subsequent degradation. During infection, the viral miRNA suppresses the population of cellular miRNA and via this pathway makes the mRNA and protein populations larger. This effect becomes appreciable with the progress of intracellular viral replication. Using biologically reasonable parameters, I investigate the corresponding mean-field kinetics and show the scale of the effect of viral miRNAs on cellular miRNA and mRNA. The scale of fluctuations of the populations of these species is illustrated as well by employing Monte Carlo simulations.

Coupling gene expression dynamics to cell size dynamics and cell cycle events: Exact and approximate solutions of the extended telegraph model

The standard model describing the fluctuations of mRNA numbers in single cells is the telegraph model which includes synthesis and degradation of mRNA, and switching of the gene between active and inactive states. While commonly used, this model does not describe how fluctuations are influenced by the cell cycle phase, cellular growth and division, and other crucial aspects of cellular biology. Here, we derive the analytical time-dependent solution of an extended telegraph model that explicitly considers the doubling of gene copy numbers upon DNA replication, dependence of the mRNA synthesis rate on cellular volume, gene dosage compensation, partitioning of molecules during cell division, cell-cycle duration variability, and cell-size control strategies. Based on the time-dependent solution, we obtain the analytical distributions of transcript numbers for lineage and population measurements in steady-state growth and also find a linear relation between the Fano factor of mRNA fluctuations and cell volume fluctuations. We show that generally the lineage and population distributions in steady-state growth cannot be accurately approximated by the steady-state solution of extrinsic noise models, i.e. a telegraph model with parameters drawn from probability distributions. This is because the mRNA lifetime is often not small enough compared to the cell cycle duration to erase the memory of division and replication. Accurate approximations are possible when this memory is weak, e.g. for genes with bursty expression and for which there is sufficient gene dosage compensation when replication occurs.

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.

Quantifying biochemical reaction rates from static population variability within incompletely observed complex networks

Quantifying biochemical reaction rates within complex cellular processes remains a key challenge of systems biology even as high-throughput single-cell data have become available to characterize snapshots of population variability. That is because complex systems with stochastic and non-linear interactions are difficult to analyze when not all components can be observed simultaneously and systems cannot be followed over time. Instead of using descriptive statistical models, we show that incompletely specified mechanistic models can be used to translate qualitative knowledge of interactions into reaction rate functions from covariability data between pairs of components. This promises to turn a globally intractable problem into a sequence of solvable inference problems to quantify complex interaction networks from incomplete snapshots of their stochastic fluctuations.

The nonequilibrium mechanism of noise-enhanced drug synergy in HIV latency reactivation

Noise-modulating chemicals can synergize with transcriptional activators in reactivating latent HIV to eliminate latent HIV reservoirs. To understand the underlying biomolecular mechanism, we investigate a previous two-gene-state model and identify two necessary conditions for the synergy: an assumption of the inhibition effect of transcription activators on noise enhancers; and frequent transitions to the gene non-transcription-permissive state. We then develop a loop-four-gene-state model with Tat transcription/translation and find that drug synergy is mainly determined by the magnitude and direction of energy input into the genetic regulatory kinetics of the HIV promoter. The inhibition effect of transcription activators is actually a phenomenon of energy dissipation in the nonequilibrium gene transition system. Overall, the loop-four-state model demonstrates that energy dissipation plays a crucial role in HIV latency reactivation, which might be useful for improving drug effects and identifying other synergies on lentivirus latency reactivation.

•

The inhibition of Activator on Noise enhancer is necessary for their synergy in reactivating HIV

•

The drug synergy is a nonequilibrium phenomenon in the gene regulatory system

•

The magnitude and direction of energy input determine the drug synergy

•

This nonequilibrium mechanism is general without regarding molecular details

Systematic inference identifies a major source of heterogeneity in cell signaling dynamics: The rate-limiting step number

Identifying the sources of cell-to-cell variability in signaling dynamics is essential to understand drug response variability and develop effective therapeutics. However, it is challenging because not all signaling intermediate reactions can be experimentally measured simultaneously. This can be overcome by replacing them with a single random time delay, but the resulting process is non-Markovian, making it difficult to infer cell-to-cell heterogeneity in reaction rates and time delays. To address this, we developed an efficient and scalable moment-based Bayesian inference method (MBI) with a user-friendly computational package that infers cell-to-cell heterogeneity in the non-Markovian signaling process. We applied MBI to single-cell expression profiles from promoters responding to antibiotics and discovered a major source of cell-to-cell variability in antibiotic stress response: the number of rate-limiting steps in signaling cascades. This knowledge can help identify effective therapies that destroy all pathogenic or cancer cells, and the approach can be applied to precision medicine.

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.

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.

Pathway dynamics can delineate the sources of transcriptional noise in gene expression

Single-cell expression profiling opens up new vistas on cellular processes. Extensive cell-to-cell variability at the transcriptomic and proteomic level has been one of the stand-out observations. Because most experimental analyses are destructive we only have access to snapshot data of cellular states. This loss of temporal information presents significant challenges for inferring dynamics, as well as causes of cell-to-cell variability. In particular, we typically cannot separate dynamic variability from within cells ('intrinsic noise') from variability across the population ('extrinsic noise'). Here, we make this non-identifiability mathematically precise, allowing us to identify new experimental set-ups that can assist in resolving this non-identifiability. We show that multiple generic reporters from the same biochemical pathways (e.g. mRNA and protein) can infer magnitudes of intrinsic and extrinsic transcriptional noise, identifying sources of heterogeneity. Stochastic simulations support our theory, and demonstrate that 'pathway-reporters' compare favourably to the well-known, but often difficult to implement, dual-reporter method.

Hierarchical Bayesian models of transcriptional and translational regulation processes with delays

MOTIVATION

Simultaneous recordings of gene network dynamics across large populations have revealed that cell characteristics vary considerably even in clonal lines. Inferring the variability of parameters that determine gene dynamics is key to understanding cellular behavior. However, this is complicated by the fact that the outcomes and effects of many reactions are not observable directly. Unobserved reactions can be replaced with time delays to reduce model dimensionality and simplify inference. However, the resulting models are non-Markovian, and require the development of new inference techniques.

RESULTS

We propose a non-Markovian, hierarchical Bayesian inference framework for quantifying the variability of cellular processes within and across cells in a population. We illustrate our approach using a delayed birth-death process. In general, a distributed delay model, rather than a popular fixed delay model, is needed for inference, even if only mean reaction delays are of interest. Using in silico and experimental data we show that the proposed hierarchical framework is robust and leads to improved estimates compared to its non-hierarchical counterpart. We apply our method to data obtained using time-lapse microscopy and infer the parameters that describe the dynamics of protein production at the single cell and population level. The mean delays in protein production are larger than previously reported, have a coefficient of variation of around 0.2 across the population, and are not strongly correlated with protein production or growth rates.

AVAILABILITY AND IMPLEMENTATION

Accompanying code in Python is available at https://github.com/mvcortez/Bayesian-Inference.

CONTACT

kresimir.josic@gmail.com or jaekkim@kaist.ac.kr or cbskust@korea.ac.kr.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Correlation between external regulators governs the mean-noise relationship in stochastic gene expression

Gene transcription in single cells is inherently a probabilistic process. The relationship between variance ($ \sigma^{2} $) and mean expression ($ \mu $) is of paramount importance for investigations into the evolutionary origins and consequences of noise in gene expression. It is often formulated as $ \log \left({{{\sigma}^{2}}}/{{{\mu}^{2}}}\; \right) = \beta\log\mu+\log\alpha $, where $ \beta $ is a key parameter since its sign determines the qualitative dependence of noise on mean. We reveal that the sign of $ \beta $ is controlled completely by external regulation, but independent of promoter structure. Specifically, it is negative if regulators as stochastic variables are independent but positive if they are correlated. The essential mechanism revealed here can well interpret diverse experimental phenomena underlying expression noise. Our results imply that external regulation rather than promoter sequence governs the mean-noise relationship.

Coordination of gene expression noise with cell size: analytical results for agent-based models of growing cell populations

The chemical master equation and the Gillespie algorithm are widely used to model the reaction kinetics inside living cells. It is thereby assumed that cell growth and division can be modelled through effective dilution reactions and extrinsic noise sources. We here re-examine these paradigms through developing an analytical agent-based framework of growing and dividing cells accompanied by an exact simulation algorithm, which allows us to quantify the dynamics of virtually any intracellular reaction network affected by stochastic cell size control and division noise. We find that the solution of the chemical master equation-including static extrinsic noise-exactly agrees with the agent-based formulation when the network under study exhibits stochastic concentration homeostasis, a novel condition that generalizes concentration homeostasis in deterministic systems to higher order moments and distributions. We illustrate stochastic concentration homeostasis for a range of common gene expression networks. When this condition is not met, we demonstrate by extending the linear noise approximation to agent-based models that the dependence of gene expression noise on cell size can qualitatively deviate from the chemical master equation. Surprisingly, the total noise of the agent-based approach can still be well approximated by extrinsic noise models.

Cell size distribution of lineage data: analytic results and parameter inference

Recent advances in single-cell technologies have enabled time-resolved measurements of the cell size over several cell cycles. These data encode information on how cells correct size aberrations so that they do not grow abnormally large or small. Here, we formulate a piecewise deterministic Markov model describing the evolution of the cell size over many generations, for all three cell size homeostasis strategies (timer, sizer, and adder). The model is solved to obtain an analytical expression for the non-Gaussian cell size distribution in a cell lineage; the theory is used to understand how the shape of the distribution is influenced by the parameters controlling the dynamics of the cell cycle and by the choice of cell tracking protocol. The theoretical cell size distribution is found to provide an excellent match to the experimental cell size distribution of E. coli lineage data collected under various growth conditions.

Asymmetry of nanoparticle inheritance upon cell division: Effect on the coefficient of variation

Several previous studies have shown that when a cell that has taken up nanoparticles divides, the nanoparticles are inherited by the two daughter cells in an asymmetrical fashion, with one daughter cell receiving more nanoparticles than the other. This interesting observation is typically demonstrated either indirectly using mathematical modelling of high-throughput experimental data or more directly by imaging individual cells as they divide. Here we suggest that measurements of the coefficient of variation (standard deviation over mean) of the number of nanoparticles per cell over the cell population is another means of assessing the degree of asymmetry. Using simulations of an evolving cell population, we show that the coefficient of variation is sensitive to the degree of asymmetry and note its characteristic evolution in time. As the coefficient of variation is readily measurable using high-throughput techniques, this should allow a more rapid experimental assessment of the degree of asymmetry.

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.

Stochastic reaction networks in dynamic compartment populations

Many biochemical processes in living systems take place in compartmentalized environments, where individual compartments can interact with each other and undergo dynamic remodeling. Studying such processes through mathematical models poses formidable challenges because the underlying dynamics involve a large number of states, which evolve stochastically with time. Here we propose a mathematical framework to study stochastic biochemical networks in compartmentalized environments. We develop a generic population model, which tracks individual compartments and their molecular composition. We then show how the time evolution of this system can be studied effectively through a small number of differential equations, which track the statistics of the population. Our approach is versatile and renders an important class of biological systems computationally accessible.

Compartmentalization of biochemical processes underlies all biological systems, from the organelle to the tissue scale. Theoretical models to study the interplay between noisy reaction dynamics and compartmentalization are sparse, and typically very challenging to analyze computationally. Recent studies have made progress toward addressing this problem in the context of specific biological systems, but a general and sufficiently effective approach remains lacking. In this work, we propose a mathematical framework based on counting processes that allows us to study dynamic compartment populations with arbitrary interactions and internal biochemistry. We derive an efficient description of the dynamics in terms of differential equations which capture the statistics of the population. We demonstrate the relevance of our approach by analyzing models inspired by different biological processes, including subcellular compartmentalization and tissue homeostasis.

Fluctuation relations and fitness landscapes of growing cell populations

We construct a pathwise formulation of a growing population of cells, based on two different samplings of lineages within the population, namely the forward and backward samplings. We show that a general symmetry relation, called fluctuation relation relates these two samplings, independently of the model used to generate divisions and growth in the cell population. These relations lead to estimators of the population growth rate, which can be very efficient as we demonstrate by an analysis of a set of mother machine data. These fluctuation relations lead to general and important inequalities between the mean number of divisions and the doubling time of the population. We also study the fitness landscape, a concept based on the two samplings mentioned above, which quantifies the correlations between a phenotypic trait of interest and the number of divisions. We obtain explicit results when the trait is the age or the size, for age and size-controlled models.

Effects of cell cycle variability on lineage and population measurements of messenger RNA abundance

Many models of gene expression do not explicitly incorporate a cell cycle description. Here, we derive a theory describing how messenger RNA (mRNA) fluctuations for constitutive and bursty gene expression are influenced by stochasticity in the duration of the cell cycle and the timing of DNA replication. Analytical expressions for the moments show that omitting cell cycle duration introduces an error in the predicted mean number of mRNAs that is a monotonically decreasing function of η, which is proportional to the ratio of the mean cell cycle duration and the mRNA lifetime. By contrast, the error in the variance of the mRNA distribution is highest for intermediate values of η consistent with genome-wide measurements in many organisms. Using eukaryotic cell data, we estimate the errors in the mean and variance to be at most 3% and 25%, respectively. Furthermore, we derive an accurate negative binomial mixture approximation to the mRNA distribution. This indicates that stochasticity in the cell cycle can introduce fluctuations in mRNA numbers that are similar to the effect of bursty transcription. Finally, we show that for real experimental data, disregarding cell cycle stochasticity can introduce errors in the inference of transcription rates larger than 10%.

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.