Emergent Lévy behavior in single-cell stochastic gene expression

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.

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.

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.

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.

Spatially and Temporally Distributed Complexity-A Refreshed Framework for the Study of GRN Evolution

Irrespective of the heuristic value of interpretations of developmental processes in terms of gene regulatory networks (GRNs), larger-angle views often suffer from: (i) an inadequate understanding of the relationship between genotype and phenotype; (ii) a predominantly zoocentric vision; and (iii) overconfidence in a putatively hierarchical organization of animal body plans. Here, we constructively criticize these assumptions. First, developmental biology is pervaded by adultocentrism, but development is not necessarily egg to adult. Second, during development, many unicells undergo transcriptomic profile transitions that are comparable to those recorded in pluricellular organisms; thus, their study should not be neglected from the GRN perspective. Third, the putatively hierarchical nature of the animal body is mirrored in the GRN logic, but in relating genotype to phenotype, independent assessments of the dynamics of the regulatory machinery and the animal’s architecture are required, better served by a combinatorial than by a hierarchical approach. The trade-offs between spatial and temporal aspects of regulation, as well as their evolutionary consequences, are also discussed. Multicellularity may derive from a unicell’s sequential phenotypes turned into different but coexisting, spatially arranged cell types. In turn, polyphenism may have been a crucial mechanism involved in the origin of complex life cycles.

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.

A data-driven method to learn a jump diffusion process from aggregate biological gene expression data

Dynamic models of gene expression are urgently required. In this paper, we describe the time evolution of gene expression by learning a jump diffusion process to model the biological process directly. Our algorithm needs aggregate gene expression data as input and outputs the parameters of the jump diffusion process. The learned jump diffusion process can predict population distributions of gene expression at any developmental stage, obtain long-time trajectories for individual cells, and offer a novel approach to computing RNA velocity. Moreover, it studies biological systems from a stochastic dynamic perspective. Gene expression data at a time point, which is a snapshot of a cellular process, is treated as an empirical marginal distribution of a stochastic process. The Wasserstein distance between the empirical distribution and predicted distribution by the jump diffusion process is minimized to learn the dynamics. For the learned jump diffusion process, its trajectories correspond to the development process of cells, the stochasticity determines the heterogeneity of cells, its instantaneous rate of state change can be taken as "RNA velocity", and the changes in scales and orientations of clusters can be noticed too. We demonstrate that our method can recover the underlying nonlinear dynamics better compared to previous parametric models and the diffusion processes driven by Brownian motion for both synthetic and real world datasets. Our method is also robust to perturbations of data because the computation involves only population expectations.

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.

Dissecting transition cells from single-cell transcriptome data through multiscale stochastic dynamics

Advances in single-cell technologies allow scrutinizing of heterogeneous cell states, however, detecting cell-state transitions from snap-shot single-cell transcriptome data remains challenging. To investigate cells with transient properties or mixed identities, we present MuTrans, a method based on multiscale reduction technique to identify the underlying stochastic dynamics that prescribes cell-fate transitions. By iteratively unifying transition dynamics across multiple scales, MuTrans constructs the cell-fate dynamical manifold that depicts progression of cell-state transitions, and distinguishes stable and transition cells. In addition, MuTrans quantifies the likelihood of all possible transition trajectories between cell states using coarse-grained transition path theory. Downstream analysis identifies distinct genes that mark the transient states or drive the transitions. The method is consistent with the well-established Langevin equation and transition rate theory. Applying MuTrans to datasets collected from five different single-cell experimental platforms, we show its capability and scalability to robustly unravel complex cell fate dynamics induced by transition cells in systems such as tumor EMT, iPSC differentiation and blood cell differentiation. Overall, our method bridges data-driven and model-based approaches on cell-fate transitions at single-cell resolution.

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.

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.

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.

Efficient analysis of stochastic gene dynamics in the non-adiabatic regime using piecewise deterministic Markov processes

Single-cell experiments show that gene expression is stochastic and bursty, a feature that can emerge from slow switching between promoter states with different activities. In addition to slow chromatin and/or DNA looping dynamics, one source of long-lived promoter states is the slow binding and unbinding kinetics of transcription factors to promoters, i.e. the non-adiabatic binding regime. Here, we introduce a simple analytical framework, known as a piecewise deterministic Markov process (PDMP), that accurately describes the stochastic dynamics of gene expression in the non-adiabatic regime. We illustrate the utility of the PDMP on a non-trivial dynamical system by analysing the properties of a titration-based oscillator in the non-adiabatic limit. We first show how to transform the underlying chemical master equation into a PDMP where the slow transitions between promoter states are stochastic, but whose rates depend upon the faster deterministic dynamics of the transcription factors regulated by these promoters. We show that the PDMP accurately describes the observed periods of stochastic cycles in activator and repressor-based titration oscillators. We then generalize our PDMP analysis to more complicated versions of titration-based oscillators to explain how multiple binding sites lengthen the period and improve coherence. Last, we show how noise-induced oscillation previously observed in a titration-based oscillator arises from non-adiabatic and discrete binding events at the promoter site.

Model reduction methods for population dynamics with fast-switching environments: Reduced master equations, stochastic differential equations, and applications

We study stochastic population dynamics coupled to fast external environments and combine expansions in the inverse switching rate of the environment and a Kramers-Moyal expansion in the inverse size of the population. This leads to a series of approximation schemes, capturing both intrinsic and environmental noise. These methods provide a means of efficient simulation and we show how they can be used to obtain analytical results for the fluctuations of population dynamics in switching environments. We place the approximations in relation to existing work on piecewise-deterministic and piecewise-diffusive Markov processes. Finally, we demonstrate the accuracy and efficiency of these model-reduction methods in different research fields, including systems in biology and a model of crack propagation.

Stochastic fluctuations can reveal the feedback signs of gene regulatory networks at the single-molecule level

Understanding the relationship between spontaneous stochastic fluctuations and the topology of the underlying gene regulatory network is of fundamental importance for the study of single-cell stochastic gene expression. Here by solving the analytical steady-state distribution of the protein copy number in a general kinetic model of stochastic gene expression with nonlinear feedback regulation, we reveal the relationship between stochastic fluctuations and feedback topology at the single-molecule level, which provides novel insights into how and to what extent a feedback loop can enhance or suppress molecular fluctuations. Based on such relationship, we also develop an effective method to extract the topological information of a gene regulatory network from single-cell gene expression data. The theory is demonstrated by numerical simulations and, more importantly, validated quantitatively by single-cell data analysis of a synthetic gene circuit integrated in human kidney cells.