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

Optimisation of gene expression noise for cellular persistence against lethal events

Bacterial cell persistence, crucial for survival under adverse conditions like antibiotic exposure, is intrinsically linked to stochastic fluctuations in gene expression. Certain genes, while inhibiting growth under normal circumstances, confer tolerance to antibiotics at elevated expression levels. The occurrence of antibiotic events lead to instantaneous cellular responses with varied survival probabilities correlated with gene expression levels. Notably, cells with lower protein concentrations face higher mortality rates. This study aims to elucidate an optimal strategy for protein expression conducive to cellular survival. Through comprehensive mathematical analysis, we determine the optimal burst size and frequency that maximise cell proliferation. Furthermore, we explore how the optimal expression distribution changes as the cost of protein expression to growth escalates. Our model reveals a hysteresis phenomenon, characterised by discontinuous transitions between deterministic and stochastic optima. Intriguingly, stochastic optima possess a noise floor, representing the minimal level of fluctuations essential for optimal cellular resilience.

Transcriptional bursting dynamics in gene expression

Gene transcription is a stochastic process that occurs in all organisms. Transcriptional bursting, a critical molecular dynamics mechanism, creates significant heterogeneity in mRNA and protein levels. This heterogeneity drives cellular phenotypic diversity. Currently, the lack of a comprehensive quantitative model limits the research on transcriptional bursting. This review examines various gene expression models and compares their strengths and weaknesses to guide researchers in selecting the most suitable model for their research context. We also provide a detailed summary of the key metrics related to transcriptional bursting. We compared the temporal dynamics of transcriptional bursting across species and the molecular mechanisms influencing these bursts, and highlighted the spatiotemporal patterns of gene expression differences by utilizing metrics such as burst size and burst frequency. We summarized the strategies for modeling gene expression from both biostatistical and biochemical reaction network perspectives. Single-cell sequencing data and integrated multiomics approaches drive our exploration of cutting-edge trends in transcriptional bursting mechanisms. Moreover, we examined classical methods for parameter estimation that help capture dynamic parameters in gene expression data, assessing their merits and limitations to facilitate optimal parameter estimation. Our comprehensive summary and review of the current transcriptional burst dynamics theories provide deeper insights for promoting research on the nature of cell processes, cell fate determination, and cancer diagnosis.

Identifying Key Regulatory Genes in Drug Resistance Acquisition: Modeling Pseudotime Trajectories of Breast Cancer Single-Cell Transcriptome

Single-cell RNA-sequencing (scRNA-seq) technology has provided significant insights into cancer drug resistance at the single-cell level. However, understanding dynamic cell transitions at the molecular systems level remains limited, requiring a systems biology approach. We present an approach that combines mathematical modeling with a pseudotime analysis using time-series scRNA-seq data obtained from the breast cancer cell line MCF-7 treated with tamoxifen. Our single-cell analysis identified five distinct subpopulations, including tamoxifen-sensitive and -resistant groups. Using a single-gene mathematical model, we discovered approximately 560-680 genes out of 6000 exhibiting multistable expression states in each subpopulation, including key estrogen-receptor-positive breast cancer cell survival genes, such as RPS6KB1. A bifurcation analysis elucidated their regulatory mechanisms, and we mapped these genes into a molecular network associated with cell survival and metastasis-related pathways. Our modeling approach comprehensively identifies key regulatory genes for drug resistance acquisition, enhancing our understanding of potential drug targets in breast cancer.

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.

Effects of bursty synthesis in organelle biogenesis

A fundamental question of cell biology is how cells control the number of organelles. The processes of organelle biogenesis, namely de novo synthesis, fission, fusion, and decay, are inherently stochastic, producing cell-to-cell variability in organelle abundance. In addition, experiments suggest that the synthesis of some organelles can be bursty. We thus ask how bursty synthesis impacts intracellular organelle number distribution. We develop an organelle biogenesis model with bursty de novo synthesis by considering geometrically distributed burst sizes. We analytically solve the model in biologically relevant limits and provide exact expressions for the steady-state organelle number distributions and their means and variances. We also present approximate solutions for the whole model, complementing with exact stochastic simulations. We show that bursts generally increase the noise in organelle numbers, producing distinct signatures in noise profiles depending on different mechanisms of organelle biogenesis. We also find different shapes of organelle number distributions, including bimodal distributions in some parameter regimes. Notably, bursty synthesis broadens the parameter regime of observing bimodality compared to the 'non-bursty' case. Together, our framework utilizes number fluctuations to elucidate the role of bursty synthesis in producing organelle number heterogeneity in cells.

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.

Using steady-state formula to estimate time-dependent parameters of stochastic gene transcription models

When studying stochastic gene transcription, it is important to understand how system parameters are temporally modulated in response to varying environments. Experimentally, the dynamic distribution data of RNA copy numbers measured at multiple time points are often fitted to stochastic transcription models to estimate time-dependent parameters. However, current methods require determining which parameters are time-dependent, as well as their analytical formulas, before the optimal fit. In this study, we developed a method to estimate time-dependent parameters in a classical two-state model without prior assumptions regarding the system parameters. At each measured time point, the method fitted the dynamic distribution data using a steady-state distribution formula, in which the estimated constant parameters were approximated as time-dependent parameter values at the measured time point. The accuracy of this method can be guaranteed for RNA molecules with relatively high degradation rates and genes with relatively slow responses to induction. We quantify the accuracy of the method and implemented this method on two sets of dynamic distribution data from prokaryotic and eukaryotic cells, and revealed the temporal modulation of transcription burst size in response to environmental changes.

Effects of microRNA-mediated negative feedback on gene expression noise

MicroRNAs (miRNAs) are small noncoding RNAs that regulate gene expression post-transcriptionally in eukaryotes by binding with target mRNAs and preventing translation. miRNA-mediated feedback motifs are ubiquitous in various genetic networks that control cellular decision making. A key question is how such a feedback mechanism may affect gene expression noise. To answer this, we have developed a mathematical model to study the effects of a miRNA-dependent negative-feedback loop on mean expression and noise in target mRNAs. Combining analytics and simulations, we show the existence of an expression threshold demarcating repressed and expressed regimes in agreement with earlier studies. The steady-state mRNA distributions are bimodal near the threshold, where copy numbers of mRNAs and miRNAs exhibit enhanced anticorrelated fluctuations. Moreover, variation of negative-feedback strength shifts the threshold locations and modulates the noise profiles. Notably, the miRNA-mRNA binding affinity and feedback strength collectively shape the bimodality. We also compare our model with a direct auto-repression motif, where a gene produces its own repressor. Auto-repression fails to produce bimodal mRNA distributions as found in miRNA-based indirect repression, suggesting the crucial role of miRNAs in creating phenotypic diversity. Together, we demonstrate how miRNA-dependent negative feedback modifies the expression threshold and leads to a broader parameter regime of bimodality compared to the no-feedback case.

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.

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.

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.

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.

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.

Effect of high variation in transcript expression on identifying differentially expressed genes in RNA-seq analysis

Great efforts have been made on the algorithms that deal with RNA-seq data to enhance the accuracy and efficiency of differential expression (DE) analysis. However, no consensus has been reached on the proper threshold values of fold change and adjusted p-value for filtering differentially expressed genes (DEGs). It is generally believed that the more stringent the filtering threshold, the more reliable the result of a DE analysis. Nevertheless, by analyzing the impact of both adjusted p-value and fold change thresholds on DE analyses, with RNA-seq data obtained for three different cancer types from the Cancer Genome Atlas (TCGA) database, we found that, for a given sample size, the reproducibility of DE results became poorer when more stringent thresholds were applied. No matter which threshold level was applied, the overlap rates of DEGs were generally lower for small sample sizes than for large sample sizes. The raw read count analysis demonstrated that the transcript expression of the same gene in different samples, whether in tumor groups or in normal groups, showed high variations, which resulted in a drastic fluctuation in fold change values and adjustedp-values when different sets of samples were used. Overall, more stringent thresholds did not yield more reliable DEGs due to high variations in transcript expression; the reliability of DEGs obtained with small sample sizes was more susceptible to these variations. Therefore, less stringent thresholds are recommended for screening DEGs. Moreover, large sample sizes should be considered in RNA-seq experimental designs to reduce the interfering effect of variations in transcript expression on DEG identification.

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.