The slow-scale linear noise approximation: an accurate, reduced stochastic description of biochemical networks under timescale separation conditions

Sequestration of gene products by decoys enhances precision in the timing of intracellular events

Expressed gene products often interact ubiquitously with binding sites at nucleic acids and macromolecular complexes, known as decoys. The binding of transcription factors (TFs) to decoys can be crucial in controlling the stochastic dynamics of gene expression. Here, we explore the impact of decoys on the timing of intracellular events, as captured by the time taken for the levels of a given TF to reach a critical threshold level, known as the first passage time (FPT). Although nonlinearity introduced by binding makes exact mathematical analysis challenging, employing suitable approximations and reformulating FPT in terms of an alternative variable, we analytically assess the impact of decoys. The stability of the decoy-bound TFs against degradation impacts FPT statistics crucially. Decoys reduce noise in FPT, and stable decoy-bound TFs offer greater timing precision with less expression cost than their unstable counterparts. Interestingly, when both bound and free TFs decay at the same rate, decoy binding does not directly alter FPT noise. We verify these results by performing exact stochastic simulations. These results have important implications for the precise temporal scheduling of events involved in the functioning of biomolecular clocks, development processes, cell-cycle control, and cell-size homeostasis.

Analysis of a detailed multi-stage model of stochastic gene expression using queueing theory and model reduction

We introduce a biologically detailed, stochastic model of gene expression describing the multiple rate-limiting steps of transcription, nuclear pre-mRNA processing, nuclear mRNA export, cytoplasmic mRNA degradation and translation of mRNA into protein. The processes in sub-cellular compartments are described by an arbitrary number of processing stages, thus accounting for a significantly finer molecular description of gene expression than conventional models such as the telegraph, two-stage and three-stage models of gene expression. We use two distinct tools, queueing theory and model reduction using the slow-scale linear-noise approximation, to derive exact or approximate analytic expressions for the moments or distributions of nuclear mRNA, cytoplasmic mRNA and protein fluctuations, as well as lower bounds for their Fano factors in steady-state conditions. We use these to study the phase diagram of the stochastic model; in particular we derive parametric conditions determining three types of transitions in the properties of mRNA fluctuations: from sub-Poissonian to super-Poissonian noise, from high noise in the nucleus to high noise in the cytoplasm, and from a monotonic increase to a monotonic decrease of the Fano factor with the number of processing stages. In contrast, protein fluctuations are always super-Poissonian and show weak dependence on the number of mRNA processing stages. Our results delineate the region of parameter space where conventional models give qualitatively incorrect results and provide insight into how the number of processing stages, e.g. the number of rate-limiting steps in initiation, splicing and mRNA degradation, shape stochastic gene expression by modulation of molecular memory.

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.

First-passage time analysis of diffusion-controlled reactions in single-molecule detection

Single-molecule detection (SMD) aims to achieve the ultimate limit-of-detection (LOD) in biosensing. This method detects a countable number of targeted analyte molecules in solution, where the dynamics of molecule diffusion, capturing, identification and delivery greatly impact the SMD's efficiency and accuracy. In this study, we adopt the first-passage time method to investigate the diffusion-controlled reaction process in SMD. We analyze the influence of detection conditions on incubation time and the expected coefficient of variation (CV) under three SMD molecule capturing strategies, including solid-phase capturing (one-dimensional solid-liquid interface fixation), liquid-phase magnetic bead (MB) capturing, and liquid-phase direct fluorescence pair labeling. We find that inside a finite-sized reaction chamber, a finite average reaction time exists in all three capturing strategies, while the liquid-phase strategies are in general more efficient than the solid-phase approaches. CV can be estimated by averaging first-passage time solely in all three strategies, and the CV reduction is achievable given an extended reaction time. To further enable zeptomolar detection, extra treatments, such as adopting liquid-phase fluorescence pairs with high diffusion rates to label the molecule, or designing specific sensing devices with large effective sensing areas would be required. This framework provides solid theoretical support to guide the design of SMD sensing strategies and sensor structures to achieve desired measurement time and CV.

Model reduction for the Chemical Master Equation: An information-theoretic approach

The complexity of mathematical models in biology has rendered model reduction an essential tool in the quantitative biologist's toolkit. For stochastic reaction networks described using the Chemical Master Equation, commonly used methods include time-scale separation, Linear Mapping Approximation, and state-space lumping. Despite the success of these techniques, they appear to be rather disparate, and at present, no general-purpose approach to model reduction for stochastic reaction networks is known. In this paper, we show that most common model reduction approaches for the Chemical Master Equation can be seen as minimizing a well-known information-theoretic quantity between the full model and its reduction, the Kullback-Leibler divergence defined on the space of trajectories. This allows us to recast the task of model reduction as a variational problem that can be tackled using standard numerical optimization approaches. In addition, we derive general expressions for propensities of a reduced system that generalize those found using classical methods. We show that the Kullback-Leibler divergence is a useful metric to assess model discrepancy and to compare different model reduction techniques using three examples from the literature: an autoregulatory feedback loop, the Michaelis-Menten enzyme system, and a genetic oscillator.

Stochastic enzyme kinetics and the quasi-steady-state reductions: Application of the slow scale linear noise approximation à la Fenichel

The linear noise approximation models the random fluctuations from the mean-field model of a chemical reaction that unfolds near the thermodynamic limit. Specifically, the fluctuations obey a linear Langevin equation up to order [Formula: see text], where [Formula: see text] is the size of the chemical system (usually the volume). In the presence of disparate timescales, the linear noise approximation admits a quasi-steady-state reduction referred to as the slow scale linear noise approximation (ssLNA). Curiously, the ssLNAs reported in the literature are slightly different. The differences in the reported ssLNAs lie at the mathematical heart of the derivation. In this work, we derive the ssLNA directly from geometric singular perturbation theory and explain the origin of the different ssLNAs in the literature. Moreover, we discuss the loss of normal hyperbolicity and we extend the ssLNA derived from geometric singular perturbation theory to a non-classical singularly perturbed problem. In so doing, we disprove a commonly-accepted qualifier for the validity of stochastic quasi-steady-state approximation of the Michaelis -Menten reaction mechanism.

Nonmodular oscillator and switch based on RNA decay drive regeneration of multimodal gene expression

Periodic gene expression dynamics are key to cell and organism physiology. Studies of oscillatory expression have focused on networks with intuitive regulatory negative feedback loops, leaving unknown whether other common biochemical reactions can produce oscillations. Oscillation and noise have been proposed to support mammalian progenitor cells’ capacity to restore heterogenous, multimodal expression from extreme subpopulations, but underlying networks and specific roles of noise remained elusive. We use mass-action-based models to show that regulated RNA degradation involving as few as two RNA species—applicable to nearly half of human protein-coding genes—can generate sustained oscillations without explicit feedback. Diverging oscillation periods synergize with noise to robustly restore cell populations’ bimodal expression on timescales of days. The global bifurcation organizing this divergence relies on an oscillator and bistable switch which cannot be decomposed into two structural modules. Our work reveals surprisingly rich dynamics of post-transcriptional reactions and a potentially widespread mechanism underlying development, tissue regeneration, and cancer cell heterogeneity.

Pathway-resolved decomposition demonstrates correlation and noise dependencies of redundant information processing in recurrent feed-forward topologies

In a biochemical assay that converts fan-in networks into feed-forward loops (FFLs), we show that the inter-regulator redundant information about the output gene product can be decomposed into finer components, mediated by the constituent pathways. Variance-based information within the linear noise regime facilitates quantifying these submodular redundancies. Contrary to the conventional wisdom on information decomposition, we report that information redundancy depends nontrivially on inter-regulator correlation. For the type-1 coherent (C1) and incoherent (I1) FFLs, the direct regulatory path-mediated redundancy is certainly correlation independent. However, components induced by the indirect regulatory path and interpathway interference are correlation dependent in (non)linear fashion. The trade-off between information redundancy and similarly decomposable extrinsic noise from input to output node has been demonstrated for the pathways and full motifs. Our analyses suggest that the interpathway cross redundancy positively and negatively influences the superposition of elementary redundancies in the C1- and I1-FFLs, respectively. Their corresponding total extrinsic noise is produced by the weighted sum and difference of the pathway-specific components. We find that the I1-FFL is able to manufacture more varied redundancy and extrinsic noise responses compared to the C1-FFL. Underlying the differing characteristics of the composite metrics across FFL variants, there exist uniformly behaving pathway-dependent elements. The decomposition framework has been meticulously explored in biologically rational parametric realizations through analytical estimates and stochastic simulations.

Periodic synchronization of isolated network elements facilitates simulating and inferring gene regulatory networks including stochastic molecular kinetics

Background

The temporal progression of many fundamental processes in cells and organisms, including homeostasis, differentiation and development, are governed by gene regulatory networks (GRNs). GRNs balance fluctuations in the output of their genes, which trace back to the stochasticity of molecular interactions. Although highly desirable to understand life processes, predicting the temporal progression of gene products within a GRN is challenging when considering stochastic events such as transcription factor–DNA interactions or protein production and degradation.

Results

We report a method to simulate and infer GRNs including genes and biochemical reactions at molecular detail. In our approach, we consider each network element to be isolated from other elements during small time intervals, after which we synchronize molecule numbers across all network elements. Thereby, the temporal behaviour of network elements is decoupled and can be treated by local stochastic or deterministic solutions. We demonstrate the working principle of this modular approach with a repressive gene cascade comprising four genes. By considering a deterministic time evolution within each time interval for all elements, our method approaches the solution of the system of deterministic differential equations associated with the GRN. By allowing genes to stochastically switch between on and off states or by considering stochastic production of gene outputs, we are able to include increasing levels of stochastic detail and approximate the solution of a Gillespie simulation. Thereby, CaiNet is able to reproduce noise-induced bi-stability and oscillations in dynamically complex GRNs. Notably, our modular approach further allows for a simple consideration of deterministic delays. We further infer relevant regulatory connections and steady-state parameters of a GRN of up to ten genes from steady-state measurements by identifying each gene of the network with a single perceptron in an artificial neuronal network and using a gradient decent method originally designed to train recurrent neural networks. To facilitate setting up GRNs and using our simulation and inference method, we provide a fast computer-aided interactive network simulation environment, CaiNet.

Conclusion

We developed a method to simulate GRNs at molecular detail and to infer the topology and steady-state parameters of GRNs. Our method and associated user-friendly framework CaiNet should prove helpful to analyze or predict the temporal progression of reaction networks or GRNs in cellular and organismic biology. CaiNet is freely available at https://gitlab.com/GebhardtLab/CaiNet.

Supplementary Information

The online version contains supplementary material available at 10.1186/s12859-021-04541-6.

On the Validity of the Stochastic Quasi-Steady-State Approximation in Open Enzyme Catalyzed Reactions: Timescale Separation or Singular Perturbation?

The quasi-steady-state approximation is widely used to develop simplified deterministic or stochastic models of enzyme catalyzed reactions. In deterministic models, the quasi-steady-state approximation can be mathematically justified from singular perturbation theory. For several closed enzymatic reactions, the homologous extension of the quasi-steady-state approximation to the stochastic regime, known as the stochastic quasi-steady-state approximation, has been shown to be accurate under the analogous conditions that permit the quasi-steady-state reduction in the deterministic counterpart. However, it was recently demonstrated that the extension of the stochastic quasi-steady-state approximation to an open Michaelis-Menten reaction mechanism is only valid under a condition that is far more restrictive than the qualifier that ensures the validity of its corresponding deterministic quasi-steady-state approximation. In this paper, we suggest a possible explanation for this discrepancy from the lens of geometric singular perturbation theory. In so doing, we illustrate a misconception in the application of the quasi-steady-state approximation: timescale separation does not imply singular perturbation.

Stochastic pH Oscillations in a Model of the Urea-Urease Reaction Confined to Lipid Vesicles

The urea-urease clock reaction is a pH switch from acid to basic that can turn into a pH oscillator if it occurs inside a suitable open reactor. We numerically study the confinement of the reaction to lipid vesicles, which permit the exchange with an external reservoir by differential transport, enabling the recovery of the pH level and yielding a constant supply of urea molecules. For microscopically small vesicles, the discreteness of the number of molecules requires a stochastic treatment of the reaction dynamics. Our analysis shows that intrinsic noise induces a significant statistical variation of the oscillation period, which increases as the vesicles become smaller. The mean period, however, is found to be remarkably robust for vesicle sizes down to approximately 200 nm, but the periodicity of the rhythm is gradually destroyed for smaller vesicles. The observed oscillations are explained as a canard-like limit cycle that differs from the wide class of conventional feedback oscillators.

Information transmission in a two-step cascade: interplay of activation and repression

We present an information-theoretic formalism to study signal transduction in four architectural variants of a model two-step cascade with increasing input population. Our results categorize these four types into two classes depending upon the effect of activation and repression on mutual information, net synergy, and signal-to-noise ratio. Using the Gaussian framework and linear noise approximation, we derive the analytic expressions for these metrics to establish their underlying relationships in terms of the biochemical parameters. We also verify our approximations through stochastic simulations.

Universally valid reduction of multiscale stochastic biochemical systems using simple non-elementary propensities

Biochemical systems consist of numerous elementary reactions governed by the law of mass action. However, experimentally characterizing all the elementary reactions is nearly impossible. Thus, over a century, their deterministic models that typically contain rapid reversible bindings have been simplified with non-elementary reaction functions (e.g., Michaelis-Menten and Morrison equations). Although the non-elementary reaction functions are derived by applying the quasi-steady-state approximation (QSSA) to deterministic systems, they have also been widely used to derive propensities for stochastic simulations due to computational efficiency and simplicity. However, the validity condition for this heuristic approach has not been identified even for the reversible binding between molecules, such as protein-DNA, enzyme-substrate, and receptor-ligand, which is the basis for living cells. Here, we find that the non-elementary propensities based on the deterministic total QSSA can accurately capture the stochastic dynamics of the reversible binding in general. However, serious errors occur when reactant molecules with similar levels tightly bind, unlike deterministic systems. In that case, the non-elementary propensities distort the stochastic dynamics of a bistable switch in the cell cycle and an oscillator in the circadian clock. Accordingly, we derive alternative non-elementary propensities with the stochastic low-state QSSA, developed in this study. This provides a universally valid framework for simplifying multiscale stochastic biochemical systems with rapid reversible bindings, critical for efficient stochastic simulations of cell signaling and gene regulation. To facilitate the framework, we provide a user-friendly open-source computational package, ASSISTER, that automatically performs the present framework.

Variance Reduction with Array-RQMC for Tau-Leaping Simulation of Stochastic Biological and Chemical Reaction Networks

We explore the use of Array-RQMC, a randomized quasi-Monte Carlo method designed for the simulation of Markov chains, to reduce the variance when simulating stochastic biological or chemical reaction networks with [Formula: see text]-leaping. The task is to estimate the expectation of a function of molecule copy numbers at a given future time T by the sample average over n sample paths, and the goal is to reduce the variance of this sample-average estimator. We find that when the method is properly applied, variance reductions by factors in the thousands can be obtained. These factors are much larger than those observed previously by other authors who tried RQMC methods for the same examples. Array-RQMC simulates an array of realizations of the Markov chain and requires a sorting function to reorder these chains according to their states, after each step. The choice of sorting function is a key ingredient for the efficiency of the method, although in our experiments, Array-RQMC was never worse than ordinary Monte Carlo, regardless of the sorting method. The expected number of reactions of each type per step also has an impact on the efficiency gain.

Stationary distributions via decomposition of stochastic reaction networks

We examine reaction networks (CRNs) through their associated continuous-time Markov processes. Studying the dynamics of such networks is in general hard, both analytically and by simulation. In particular, stationary distributions of stochastic reaction networks are only known in some cases. We analyze class properties of the underlying continuous-time Markov chain of CRNs under the operation of join and examine conditions such that the form of the stationary distributions of a CRN is derived from the parts of the decomposed CRNs. The conditions can be easily checked in examples and allow recursive application. The theory developed enables sequential decomposition of the Markov processes and calculations of stationary distributions. Since the class of processes expressible through such networks is big and only few assumptions are made, the principle also applies to other stochastic models. We give examples of interest from CRN theory to highlight the decomposition.

Extrinsic noise of the target gene governs abundance pattern of feed-forward loop motifs

Feed-forward loop (FFL) is found to be a recurrent structure in bacterial and yeast gene transcription regulatory networks. In a generic FFL, transcription factor (TF) S regulates production of another TF X while both of these TFs regulate production of final gene-product Y. Depending upon the regulatory programs (activation or repression), FFLs are grouped into two broad classes: coherent (C) and incoherent (I), each class containing four distinct types (C1-C4 and I1-I4). These FFL types are experimentally observed to occur with varied frequencies, C1 and I1 being the abundant ones. Here we present a stochastic framework singling out the absolute value of the normalized covariance of X and Y to be the determining factor behind the abundance of FFLs while considering differential promoter activities of X and Y. Our theoretical construct employs two possible signal integration mechanisms (additive and multiplicative) to synthesize Y while steady-state population level of S remains fixed or becomes tunable reflecting two possible environmental signaling scenarios. Our model categorically points out that abundant FFLs exhibit higher amount of the designated metric which has a biophysical connotation of extrinsic noise for the target gene Y. Our predictions emanating from an overarching analytical expression utilizing biologically plausible parametric conditions are substantiated by stochastic simulation.

On rapid oscillations driving biological processes at disparate timescales

We consider a generic biological process described by a dynamical system, subject to an input signal with a high-frequency periodic component. The rapid oscillations of the input signal induce inherently multiscale dynamics, motivating order-reduction techniques. It is intuitive that the system behaviour is well approximated by its response to the averaged input signal. However, changes to the high-frequency component that preserve the average signal are beyond the reach of such intuitive reasoning. In this study, we explore system response under the influence of such an input signal by exploiting the timescale separation between high-frequency input variations and system response time. Employing the asymptotic method of multiple scales, we establish that, in some circumstances, the intuitive approach is simply the leading-order asymptotic contribution. We focus on higher-order corrections that capture the response to the details of the high-frequency component beyond its average. This approach achieves a reduction in system complexity while providing valuable insight into the structure of the response to the oscillations. We develop the general theory for nonlinear systems, while highlighting the important case of systems affine in the state and the input signal, presenting examples of both discrete and continuum state spaces. Importantly, this class of systems encompasses biochemical reaction networks described by the chemical master equation and its continuum approximations. Finally, we apply the framework to a nonlinear system describing mRNA translation and protein expression previously studied in the literature. The analysis shines new light on several aspects of the system quantification and both extends and simplifies results previously obtained.

Stochastic time-dependent enzyme kinetics: Closed-form solution and transient bimodality

We derive an approximate closed-form solution to the chemical master equation describing the Michaelis-Menten reaction mechanism of enzyme action. In particular, assuming that the probability of a complex dissociating into an enzyme and substrate is significantly larger than the probability of a product formation event, we obtain expressions for the time-dependent marginal probability distributions of the number of substrate and enzyme molecules. For delta function initial conditions, we show that the substrate distribution is either unimodal at all times or else becomes bimodal at intermediate times. This transient bimodality, which has no deterministic counterpart, manifests when the initial number of substrate molecules is much larger than the total number of enzyme molecules and if the frequency of enzyme-substrate binding events is large enough. Furthermore, we show that our closed-form solution is different from the solution of the chemical master equation reduced by means of the widely used discrete stochastic Michaelis-Menten approximation, where the propensity for substrate decay has a hyperbolic dependence on the number of substrate molecules. The differences arise because the latter does not take into account enzyme number fluctuations, while our approach includes them. We confirm by means of a stochastic simulation of all the elementary reaction steps in the Michaelis-Menten mechanism that our closed-form solution is accurate over a larger region of parameter space than that obtained using the discrete stochastic Michaelis-Menten approximation.

Misuse of the Michaelis-Menten rate law for protein interaction networks and its remedy

For over a century, the Michaelis-Menten (MM) rate law has been used to describe the rates of enzyme-catalyzed reactions and gene expression. Despite the ubiquity of the MM rate law, it accurately captures the dynamics of underlying biochemical reactions only so long as it is applied under the right condition, namely, that the substrate is in large excess over the enzyme-substrate complex. Unfortunately, in circumstances where its validity condition is not satisfied, especially so in protein interaction networks, the MM rate law has frequently been misused. In this review, we illustrate how inappropriate use of the MM rate law distorts the dynamics of the system, provides mistaken estimates of parameter values, and makes false predictions of dynamical features such as ultrasensitivity, bistability, and oscillations. We describe how these problems can be resolved with a slightly modified form of the MM rate law, based on the total quasi-steady state approximation (tQSSA). Furthermore, we show that the tQSSA can be used for accurate stochastic simulations at a lower computational cost than using the full set of mass-action rate laws. This review describes how to use quasi-steady state approximations in the right context, to prevent drawing erroneous conclusions from in silico simulations.

Basal leakage in oscillation: Coupled transcriptional and translational control using feed-forward loops

The circadian clock is a complex system that plays many important roles in most organisms. Previously, many mathematical models have been used to sharpen our understanding of the Arabidopsis clock, which brought to light the roles of each transcriptional and post-translational regulations. However, the presence of both regulations, instead of either transcription or post-translation, raised curiosity of whether the combination of these two regulations is important for the clock’s system. In this study, we built a series of simplified oscillators with different regulations to study the importance of post-translational regulation (specifically, 26S proteasome degradation) in the clock system. We found that a simple transcriptional-based oscillator can already generate sustained oscillation, but the oscillation can be easily destroyed in the presence of transcriptional leakage. Coupling post-translational control with transcriptional-based oscillator in a feed-forward loop will greatly improve the robustness of the oscillator in the presence of basal leakage. Using these general models, we were able to replicate the increased variability observed in the E3 ligase mutant for both plant and mammalian clocks. With this insight, we also predict a plausible regulator of several E3 ligase genes in the plant’s clock. Thus, our results provide insights into and the plausible importance in coupling transcription and post-translation controls in the clock system.

For circadian clocks, several current models had successfully captured the essential dynamic behavior of the clock system mainly with transcriptional regulation. Previous studies have shown that the 26S proteasome degradation controls are important in maintaining the stability of circadian rhythms. However, how the loss-of-function or over-expression mutant of this targeted degradations lead to unstable oscillation is still unclear. In this work, we investigate the importance of coupled transcriptional and post-translational feedback loop in the circadian oscillator. With general models our study indicate that the unstable behavior of degradation mutants could be caused by the increase in the basal level of the clock genes. We found that coupling a non-linear degradation control into this transcriptional based oscillator using feed-forward loop improves the robustness of the oscillator. Using this finding, we further predict some plausible regulators of Arabidopsis’s E3 ligase protein such as COP1 and SINAT5. Hence, our results provide insights on the importance of coupling transcription and post-translation controls in the clock system.

Enhancement of gene expression noise from transcription factor binding to genomic decoy sites

The genome contains several high-affinity non-functional binding sites for transcription factors (TFs) creating a hidden and unexplored layer of gene regulation. We investigate the role of such “decoy sites” in controlling noise (random fluctuations) in the level of a TF that is synthesized in stochastic bursts. Prior studies have assumed that decoy-bound TFs are protected from degradation, and in this case decoys function to buffer noise. Relaxing this assumption to consider arbitrary degradation rates for both bound/unbound TF states, we find rich noise behaviors. For low-affinity decoys, noise in the level of unbound TF always monotonically decreases to the Poisson limit with increasing decoy numbers. In contrast, for high-affinity decoys, noise levels first increase with increasing decoy numbers, before decreasing back to the Poisson limit. Interestingly, while protection of bound TFs from degradation slows the time-scale of fluctuations in the unbound TF levels, the decay of bound TFs leads to faster fluctuations and smaller noise propagation to downstream target proteins. In summary, our analysis reveals stochastic dynamics emerging from nonspecific binding of TFs and highlights the dual role of decoys as attenuators or amplifiers of gene expression noise depending on their binding affinity and stability of the bound TF.

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.

Deciphering noise amplification and reduction in open chemical reaction networks

The impact of fluctuations on the dynamical behaviour of complex biological systems is a longstanding issue, whose understanding would elucidate how evolutionary pressure tends to modulate intrinsic noise. Using the Itō stochastic differential equation formalism, we performed analytic and numerical analyses of model systems containing different molecular species in contact with the environment and interacting with each other through mass-action kinetics. For networks of zero deficiency, which admit a detailed- or complex-balanced steady state, all molecular species are uncorrelated and their Fano factors are Poissonian. Systems of higher deficiency have non-equilibrium steady states and non-zero reaction fluxes flowing between the complexes. When they model homo-oligomerization, the noise on each species is reduced when the flux flows from the oligomers of lowest to highest degree, and amplified otherwise. In the case of hetero-oligomerization systems, only the noise on the highest-degree species shows this behaviour.

Markovian approaches to modeling intracellular reaction processes with molecular memory

Many cellular processes are governed by stochastic reaction events. These events do not necessarily occur in single steps of individual molecules, and, conversely, each birth or death of a macromolecule (e.g., protein) could involve several small reaction steps, creating a memory between individual events and thus leading to nonmarkovian reaction kinetics. Characterizing this kinetics is challenging. Here, we develop a systematic approach for a general reaction network with arbitrary intrinsic waiting-time distributions, which includes the stationary generalized chemical-master equation (sgCME), the stationary generalized Fokker-Planck equation, and the generalized linear-noise approximation. The first formulation converts a nonmarkovian issue into a markovian one by introducing effective transition rates (that explicitly decode the effect of molecular memory) for the reactions in an equivalent reaction network with the same substrates but without molecular memory. Nonmarkovian features of the reaction kinetics can be revealed by solving the sgCME. The latter 2 formulations can be used in the fast evaluation of fluctuations. These formulations can have broad applications, and, in particular, they may help us discover new biological knowledge underlying memory effects. When they are applied to generalized stochastic models of gene-expression regulation, we find that molecular memory is in effect equivalent to a feedback and can induce bimodality, fine-tune the expression noise, and induce switch.

Investigation of chemical noise in multisite phosphorylation chain using linear noise approximation

Quantitative and qualitative nature of chemical noise propagation in biochemical reaction networks depend crucially on the topology of the networks. Multisite reversible phosphorylation-dephosphorylation of target proteins is one such recurrently found topology that regulates host of key functions in living cells. Here we analytically calculated the stochasticity in multistep reversible chemical reactions by determining variance of phosphorylated species at the steady state using linear noise approximation to investigate the effect of mass action and Michaelis-Menten kinetics on the noise of phosphorylated species. We probed the dependence of noise on the number of phosphorylation sites and the equilibrium constants of the reaction equilibria to investigate the chemical noise propagation in the multisite phosphorylation chain.

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.

Bayesian inference of distributed time delay in transcriptional and translational regulation

MOTIVATION

Advances in experimental and imaging techniques have allowed for unprecedented insights into the dynamical processes within individual cells. However, many facets of intracellular dynamics remain hidden, or can be measured only indirectly. This makes it challenging to reconstruct the regulatory networks that govern the biochemical processes underlying various cell functions. Current estimation techniques for inferring reaction rates frequently rely on marginalization over unobserved processes and states. Even in simple systems this approach can be computationally challenging, and can lead to large uncertainties and lack of robustness in parameter estimates. Therefore we will require alternative approaches to efficiently uncover the interactions in complex biochemical networks.

RESULTS

We propose a Bayesian inference framework based on replacing uninteresting or unobserved reactions with time delays. Although the resulting models are non-Markovian, recent results on stochastic systems with random delays allow us to rigorously obtain expressions for the likelihoods of model parameters. In turn, this allows us to extend MCMC methods to efficiently estimate reaction rates, and delay distribution parameters, from single-cell assays. We illustrate the advantages, and potential pitfalls, of the approach using a birth-death model with both synthetic and experimental data, and show that we can robustly infer model parameters using a relatively small number of measurements. We demonstrate how to do so even when only the relative molecule count within the cell is measured, as in the case of fluorescence microscopy.

AVAILABILITY AND IMPLEMENTATION

Accompanying code in R is available at https://github.com/cbskust/DDE_BD.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Comprehensive Review of Models and Methods for Inferences in Bio-Chemical Reaction Networks

The key processes in biological and chemical systems are described by networks of chemical reactions. From molecular biology to biotechnology applications, computational models of reaction networks are used extensively to elucidate their non-linear dynamics. The model dynamics are crucially dependent on the parameter values which are often estimated from observations. Over the past decade, the interest in parameter and state estimation in models of (bio-) chemical reaction networks (BRNs) grew considerably. The related inference problems are also encountered in many other tasks including model calibration, discrimination, identifiability, and checking, and optimum experiment design, sensitivity analysis, and bifurcation analysis. The aim of this review paper is to examine the developments in literature to understand what BRN models are commonly used, and for what inference tasks and inference methods. The initial collection of about 700 documents concerning estimation problems in BRNs excluding books and textbooks in computational biology and chemistry were screened to select over 270 research papers and 20 graduate research theses. The paper selection was facilitated by text mining scripts to automate the search for relevant keywords and terms. The outcomes are presented in tables revealing the levels of interest in different inference tasks and methods for given models in the literature as well as the research trends are uncovered. Our findings indicate that many combinations of models, tasks and methods are still relatively unexplored, and there are many new research opportunities to explore combinations that have not been considered-perhaps for good reasons. The most common models of BRNs in literature involve differential equations, Markov processes, mass action kinetics, and state space representations whereas the most common tasks are the parameter inference and model identification. The most common methods in literature are Bayesian analysis, Monte Carlo sampling strategies, and model fitting to data using evolutionary algorithms. The new research problems which cannot be directly deduced from the text mining data are also discussed.

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.

Excitable dynamics through toxin-induced mRNA cleavage in bacteria

Toxin-antitoxin (TA) systems in bacteria and archaea are small genetic elements consisting of the genes coding for an intracellular toxin and an antitoxin that can neutralize this toxin. In various cases, the toxins cleave the mRNA. In this theoretical work we use deterministic and stochastic modeling to explain how toxin-induced cleavage of mRNA in TA systems can lead to excitability, allowing large transient spikes in toxin levels to be triggered. By using a simplified network where secondary complex formation and transcriptional regulation are not included, we show that a two-dimensional, deterministic model captures the origin of such toxin excitations. Moreover, it allows to increase our understanding by examining the dynamics in the phase plane. By systematically comparing the deterministic results with Gillespie simulations we demonstrate that even though the real TA system is intrinsically stochastic, toxin excitations can be accurately described deterministically. A bifurcation analysis of the system shows that the excitable behavior is due to a nearby Hopf bifurcation in the parameter space, where the system becomes oscillatory. The influence of stress is modeled by varying the degradation rate of the antitoxin and the translation rate of the toxin. We find that stress increases the frequency of toxin excitations. The inclusion of secondary complex formation and transcriptional regulation does not fundamentally change the mechanism of toxin excitations. Finally, we show that including growth rate suppression and translational inhibition can lead to longer excitations, and even cause excitations in cases when the system would otherwise be non-excitable. To conclude, the deterministic model used in this work provides a simple and intuitive explanation of toxin excitations in TA systems.

Comparison of Deterministic and Stochastic Regime in a Model for Cdc42 Oscillations in Fission Yeast

Oscillations occur in a wide variety of essential cellular processes, such as cell cycle progression, circadian clocks and calcium signaling in response to stimuli. It remains unclear how intrinsic stochasticity can influence these oscillatory systems. Here, we focus on oscillations of Cdc42 GTPase in fission yeast. We extend our previous deterministic model by Xu and Jilkine to construct a stochastic model, focusing on the fast diffusion case. We use SSA (Gillespie's algorithm) to numerically explore the low copy number regime in this model, and use analytical techniques to study the long-time behavior of the stochastic model and compare it to the equilibria of its deterministic counterpart. Numerical solutions suggest noisy limit cycles exist in the parameter regime in which the deterministic system converges to a stable limit cycle, and quasi-cycles exist in the parameter regime where the deterministic model has a damped oscillation. Near an infinite period bifurcation point, the deterministic model has a sustained oscillation, while stochastic trajectories start with an oscillatory mode and tend to approach deterministic steady states. In the low copy number regime, metastable transitions from oscillatory to steady behavior occur in the stochastic model. Our work contributes to the understanding of how stochastic chemical kinetics can affect a finite-dimensional dynamical system, and destabilize a deterministic steady state leading to oscillations.

Linear mapping approximation of gene regulatory networks with stochastic dynamics

The presence of protein-DNA binding reactions often leads to analytically intractable models of stochastic gene expression. Here we present the linear-mapping approximation that maps systems with protein-promoter interactions onto approximately equivalent systems with no binding reactions. This is achieved by the marriage of conditional mean-field approximation and the Magnus expansion, leading to analytic or semi-analytic expressions for the approximate time-dependent and steady-state protein number distributions. Stochastic simulations verify the method's accuracy in capturing the changes in the protein number distributions with time for a wide variety of networks displaying auto- and mutual-regulation of gene expression and independently of the ratios of the timescales governing the dynamics. The method is also used to study the first-passage time distribution of promoter switching, the sensitivity of the size of protein number fluctuations to parameter perturbation and the stochastic bifurcation diagram characterizing the onset of multimodality in protein number distributions.

Memory functions reveal structural properties of gene regulatory networks

Gene regulatory networks (GRNs) control cellular function and decision making during tissue development and homeostasis. Mathematical tools based on dynamical systems theory are often used to model these networks, but the size and complexity of these models mean that their behaviour is not always intuitive and the underlying mechanisms can be difficult to decipher. For this reason, methods that simplify and aid exploration of complex networks are necessary. To this end we develop a broadly applicable form of the Zwanzig-Mori projection. By first converting a thermodynamic state ensemble model of gene regulation into mass action reactions we derive a general method that produces a set of time evolution equations for a subset of components of a network. The influence of the rest of the network, the bulk, is captured by memory functions that describe how the subnetwork reacts to its own past state via components in the bulk. These memory functions provide probes of near-steady state dynamics, revealing information not easily accessible otherwise. We illustrate the method on a simple cross-repressive transcriptional motif to show that memory functions not only simplify the analysis of the subnetwork but also have a natural interpretation. We then apply the approach to a GRN from the vertebrate neural tube, a well characterised developmental transcriptional network composed of four interacting transcription factors. The memory functions reveal the function of specific links within the neural tube network and identify features of the regulatory structure that specifically increase the robustness of the network to initial conditions. Taken together, the study provides evidence that Zwanzig-Mori projections offer powerful and effective tools for simplifying and exploring the behaviour of GRNs.

Gene regulatory networks are essential for cell fate specification and function. But the recursive links that comprise these networks often make determining their properties and behaviour complicated. Computational models of these networks can also be difficult to decipher. To reduce the complexity of such models we employ a Zwanzig-Mori projection approach. This allows a system of ordinary differential equations, representing a network, to be reduced to an arbitrary subnetwork consisting of part of the initial network, with the rest of the network (bulk) captured by memory functions. These memory functions account for the bulk by describing signals that return to the subnetwork after some time, having passed through the bulk. We show how this approach can be used to simplify analysis and to probe the behaviour of a gene regulatory network. Applying the method to a transcriptional network in the vertebrate neural tube reveals previously unappreciated properties of the network. By taking advantage of the structure of the memory functions we identify interactions within the network that are unnecessary for sustaining correct patterning. Upon further investigation we find that these interactions are important for conferring robustness to variation in initial conditions. Taken together we demonstrate the validity and applicability of the Zwanzig-Mori projection approach to gene regulatory networks.

Hill functions for stochastic gene regulatory networks from master equations with split nodes and time-scale separation

The deterministic Hill function depends only on the average values of molecule numbers. To account for the fluctuations in the molecule numbers, the argument of the Hill function needs to contain the means, the standard deviations, and the correlations. Here we present a method that allows for stochastic Hill functions to be constructed from the dynamical evolution of stochastic biocircuits with specific topologies. These stochastic Hill functions are presented in a closed analytical form so that they can be easily incorporated in models for large genetic regulatory networks. Using a repressive biocircuit as an example, we show by Monte Carlo simulations that the traditional deterministic Hill function inaccurately predicts time of repression by an order of two magnitudes. However, the stochastic Hill function was able to capture the fluctuations and thus accurately predicted the time of repression.

Beyond the Michaelis-Menten equation: Accurate and efficient estimation of enzyme kinetic parameters

Examining enzyme kinetics is critical for understanding cellular systems and for using enzymes in industry. The Michaelis-Menten equation has been widely used for over a century to estimate the enzyme kinetic parameters from reaction progress curves of substrates, which is known as the progress curve assay. However, this canonical approach works in limited conditions, such as when there is a large excess of substrate over enzyme. Even when this condition is satisfied, the identifiability of parameters is not always guaranteed, and often not verifiable in practice. To overcome such limitations of the canonical approach for the progress curve assay, here we propose a Bayesian approach based on an equation derived with the total quasi-steady-state approximation. In contrast to the canonical approach, estimates obtained with this proposed approach exhibit little bias for any combination of enzyme and substrate concentrations. Importantly, unlike the canonical approach, an optimal experiment to identify parameters with certainty can be easily designed without any prior information. Indeed, with this proposed design, the kinetic parameters of diverse enzymes with disparate catalytic efficiencies, such as chymotrypsin, fumarase, and urease, can be accurately and precisely estimated from a minimal amount of timecourse data. A publicly accessible computational package performing such accurate and efficient Bayesian inference for enzyme kinetics is provided.

Late-Arriving Signals Contribute Less to Cell-Fate Decisions

Gene regulatory networks are largely responsible for cellular decision-making. These networks sense diverse external signals and respond by adjusting gene expression, enabling cells to reach environment-dependent decisions crucial for their survival or reproduction. However, information-carrying signals may arrive at variable times. Besides the intrinsic strength of these signals, their arrival time (timing) may also carry information about the environment and can influence cellular decision-making in ways that are poorly understood. For example, it is unclear how the timing of individual phage infections affects the lysis-lysogeny decision of bacteriophage λ despite variable infection times being likely in the wild and even in laboratory conditions. In this work, we combine mathematical modeling with experimentation to address this question. We develop an experimentally testable theory, which reveals that late-infecting phages contribute less to cellular decision-making. This implies that infection delays lower the probability of lysogeny compared to simultaneous infections. Furthermore, we show that infection delays reduce lysogenization by providing insufficient CII for threshold crossing during the critical decision-making period. We find evidence for a cutoff time after which subsequent infections cannot influence the cellular decision. We derive an intuitive formula that approximates the probability of lysogeny for variable infection times by a time-weighted average of probabilities for simultaneous infections. We validate these theoretical predictions experimentally. Similar concepts and simplifying modeling approaches may help elucidate the mechanisms underlying other cellular decisions.

Statistical physics approaches to subnetwork dynamics in biochemical systems

We apply a Gaussian variational approximation to model reduction in large biochemical networks of unary and binary reactions. We focus on a small subset of variables (subnetwork) of interest, e.g. because they are accessible experimentally, embedded in a larger network (bulk). The key goal is to write dynamical equations reduced to the subnetwork but still retaining the effects of the bulk. As a result, the subnetwork-reduced dynamics contains a memory term and an extrinsic noise term with non-trivial temporal correlations. We first derive expressions for this memory and noise in the linearized (Gaussian) dynamics and then use a perturbative power expansion to obtain first order nonlinear corrections. For the case of vanishing intrinsic noise, our description is explicitly shown to be equivalent to projection methods up to quadratic terms, but it is applicable also in the presence of stochastic fluctuations in the original dynamics. An example from the epidermal growth factor receptor signalling pathway is provided to probe the increased prediction accuracy and computational efficiency of our method.

Reduction of multiscale stochastic biochemical reaction networks using exact moment derivation

Biochemical reaction networks (BRNs) in a cell frequently consist of reactions with disparate timescales. The stochastic simulations of such multiscale BRNs are prohibitively slow due to high computational cost for the simulations of fast reactions. One way to resolve this problem uses the fact that fast species regulated by fast reactions quickly equilibrate to their stationary distribution while slow species are unlikely to be changed. Thus, on a slow timescale, fast species can be replaced by their quasi-steady state (QSS): their stationary conditional expectation values for given slow species. As the QSS are determined solely by the state of slow species, such replacement leads to a reduced model, where fast species are eliminated. However, it is challenging to derive the QSS in the presence of nonlinear reactions. While various approximation schemes for the QSS have been developed, they often lead to considerable errors. Here, we propose two classes of multiscale BRNs which can be reduced by deriving an exact QSS rather than approximations. Specifically, if fast species constitute either a feedforward network or a complex balanced network, the reduced model based on the exact QSS can be derived. Such BRNs are frequently observed in a cell as the feedforward network is one of fundamental motifs of gene or protein regulatory networks. Furthermore, complex balanced networks also include various types of fast reversible bindings such as bindings between transcriptional factors and gene regulatory sites. The reduced models based on exact QSS, which can be calculated by the computational packages provided in this work, accurately approximate the slow scale dynamics of the original full model with much lower computational cost.

An open-loop approach to calculate noise-induced transitions

Bistability permits the co-existence of two distinct cell fates in a population of genetically identical cells. Noise induced transitions between two fates of a bistable system are difficult to calculate due to the intricate interplay between nonlinear dynamics and noise in bistable positive feedback loops. Here we opened multivariable feedback loops at the slowest variable to obtain the open-loop function and the fluctuations in the open-loop output. By the subsequent reclosing of the loop, we calculated the mean first passage time (MFPT) using the Fokker-Planck equation in good agreement with the exact stochastic simulation. When an external component interacts with a feedback component, it amplifies the extrinsic noise in the loop. Consequently, the open-loop function is shifted and the transition rates between the two states in the closed loop are increased. Despite this shift, the open-loop output reflects the system faithfully to predict the MFPT in the feedback loop. Therefore, the open-loop approach can help theoretical analysis. Furthermore, the measurement of the mean value, variance, and the reaction time-scale of the open-loop output permits the prediction of MFPT simply from experimental data, which underscores the practical value of the stochastic open-loop approach.

Noise Decomposition Principle in a Coherent Feed-Forward Transcriptional Regulatory Loop

Coherent feed-forward loops exist extensively in realistic biological regulatory systems, and are common signaling motifs. Here, we study the characteristics and the propagation mechanism of the output noise in a coherent feed-forward transcriptional regulatory loop that can be divided into a main road and branch. Using the linear noise approximation, we derive analytical formulae for the total noise of the full loop, the noise of the branch, and the noise of the main road, which are verified by the Gillespie algorithm. Importantly, we find that (i) compared with the branch motif or the main road motif, the full motif can effectively attenuate the output noise level; (ii) there is a transition point of system state such that the noise of the main road is dominated when the underlying system is below this point, whereas the noise of the branch is dominated when the system is beyond the point. The entire analysis reveals the mechanism of how the noise is generated and propagated in a simple yet representative signaling module.

Product-Form Stationary Distributions for Deficiency Zero Networks with Non-mass Action Kinetics

In many applications, for example when computing statistics of fast subsystems in a multiscale setting, we wish to find the stationary distributions of systems of continuous-time Markov chains. Here we present a class of models that appears naturally in certain averaging approaches whose stationary distributions can be computed explicitly. In particular, we study continuous-time Markov chain models for biochemical interaction systems with non-mass action kinetics whose network satisfies a certain constraint. Analogous with previous related results, the distributions can be written in product form.

Stochastic modeling of biochemical systems with multistep reactions using state-dependent time delay

To deal with the growing scale of molecular systems, sophisticated modelling techniques have been designed in recent years to reduce the complexity of mathematical models. Among them, a widely used approach is delayed reaction for simplifying multistep reactions. However, recent research results suggest that a delayed reaction with constant time delay is unable to describe multistep reactions accurately. To address this issue, we propose a novel approach using state-dependent time delay to approximate multistep reactions. We first use stochastic simulations to calculate time delay arising from multistep reactions exactly. Then we design algorithms to calculate time delay based on system dynamics precisely. To demonstrate the power of proposed method, two processes of mRNA degradation are used to investigate the function of time delay in determining system dynamics. In addition, a multistep pathway of metabolic synthesis is used to explore the potential of the proposed method to simplify multistep reactions with nonlinear reaction rates. Simulation results suggest that the state-dependent time delay is a promising and accurate approach to reduce model complexity and decrease the number of unknown parameters in the models.

Stochastic Simulation of Biomolecular Networks in Dynamic Environments

Simulation of biomolecular networks is now indispensable for studying biological systems, from small reaction networks to large ensembles of cells. Here we present a novel approach for stochastic simulation of networks embedded in the dynamic environment of the cell and its surroundings. We thus sample trajectories of the stochastic process described by the chemical master equation with time-varying propensities. A comparative analysis shows that existing approaches can either fail dramatically, or else can impose impractical computational burdens due to numerical integration of reaction propensities, especially when cell ensembles are studied. Here we introduce the Extrande method which, given a simulated time course of dynamic network inputs, provides a conditionally exact and several orders-of-magnitude faster simulation solution. The new approach makes it feasible to demonstrate—using decision-making by a large population of quorum sensing bacteria—that robustness to fluctuations from upstream signaling places strong constraints on the design of networks determining cell fate. Our approach has the potential to significantly advance both understanding of molecular systems biology and design of synthetic circuits.

Simulation algorithms have become indispensable tools in modern quantitative biology, providing deep insight into many biochemical systems, including gene regulatory networks. However, current stochastic simulation approaches handle the effects of fluctuating extracellular signals and upstream processes poorly, either failing to give qualitatively reliable predictions or being very inefficient computationally. Here we introduce the Extrande method, a novel approach for simulation of biomolecular networks embedded in the dynamic environment of the cell and its surroundings. The method is accurate and computationally efficient, and hence fills an important gap in the field of stochastic simulation. In particular, we employ it to study a bacterial decision-making network and demonstrate that robustness to fluctuations from upstream signaling places strong constraints on the design of networks determining cell fate.