The Validity of Quasi-Steady-State Approximations in Discrete Stochastic Simulations

In biochemical networks, reactions often occur on disparate timescales and can be characterized as either fast or slow. The quasi-steady-state approximation (QSSA) utilizes timescale separation to project models of biochemical networks onto lower-dimensional slow manifolds. As a result, fast elementary reactions are not modeled explicitly, and their effect is captured by nonelementary reaction-rate functions (e.g., Hill functions). The accuracy of the QSSA applied to deterministic systems depends on how well timescales are separated. Recently, it has been proposed to use the nonelementary rate functions obtained via the deterministic QSSA to define propensity functions in stochastic simulations of biochemical networks. In this approach, termed the stochastic QSSA, fast reactions that are part of nonelementary reactions are not simulated, greatly reducing computation time. However, it is unclear when the stochastic QSSA provides an accurate approximation of the original stochastic simulation. We show that, unlike the deterministic QSSA, the validity of the stochastic QSSA does not follow from timescale separation alone, but also depends on the sensitivity of the nonelementary reaction rate functions to changes in the slow species. The stochastic QSSA becomes more accurate when this sensitivity is small. Different types of QSSAs result in nonelementary functions with different sensitivities, and the total QSSA results in less sensitive functions than the standard or the prefactor QSSA. We prove that, as a result, the stochastic QSSA becomes more accurate when nonelementary reaction functions are obtained using the total QSSA. Our work provides an apparently novel condition for the validity of the QSSA in stochastic simulations of biochemical reaction networks with disparate timescales.     

The Validity of Quasi-Steady-State Approximations in Discrete Stochastic Simulations

In biochemical networks, reactions often occur on disparate timescales and can be characterized as either fast or slow. The quasi-steady-state approximation (QSSA) utilizes timescale separation to project models of biochemical networks onto lower-dimensional slow manifolds. As a result, fast elementary reactions are not modeled explicitly, and their effect is captured by nonelementary reaction-rate functions (e.g., Hill functions). The accuracy of the QSSA applied to deterministic systems depends on how well timescales are separated. Recently, it has been proposed to use the nonelementary rate functions obtained via the deterministic QSSA to define propensity functions in stochastic simulations of biochemical networks. In this approach, termed the stochastic QSSA, fast reactions that are part of nonelementary reactions are not simulated, greatly reducing computation time. However, it is unclear when the stochastic QSSA provides an accurate approximation of the original stochastic simulation. We show that, unlike the deterministic QSSA, the validity of the stochastic QSSA does not follow from timescale separation alone, but also depends on the sensitivity of the nonelementary reaction rate functions to changes in the slow species. The stochastic QSSA becomes more accurate when this sensitivity is small. Different types of QSSAs result in nonelementary functions with different sensitivities, and the total QSSA results in less sensitive functions than the standard or the prefactor QSSA. We prove that, as a result, the stochastic QSSA becomes more accurate when nonelementary reaction functions are obtained using the total QSSA. Our work provides an apparently novel condition for the validity of the QSSA in stochastic simulations of biochemical reaction networks with disparate timescales.     

Stochastic simulation of enzyme-catalyzed reactions with disparate timescales

Many physiological characteristics of living cells are regulated by protein interaction networks. Because the total numbers of these protein species can be small, molecular noise can have significant effects on the dynamical properties of a regulatory network. Computing these stochastic effects is made difficult by the large timescale separations typical of protein interactions (e.g., complex formation may occur in fractions of a second, whereas catalytic conversions may take minutes). Exact stochastic simulation may be very inefficient under these circumstances, and methods for speeding up the simulation without sacrificing accuracy have been widely studied. We show that the “total quasi-steady-state approximation” for enzyme-catalyzed reactions provides a useful framework for efficient and accurate stochastic simulations. The method is applied to three examples: a simple enzyme-catalyzed reaction where enzyme and substrate have comparable abundances, a Goldbeter-Koshland switch, where a kinase and phosphatase regulate the phosphorylation state of a common substrate, and coupled Goldbeter-Koshland switches that exhibit bistability. Simulations based on the total quasi-steady-state approximation accurately capture the steady-state probability distributions of all components of these reaction networks. In many respects, the approximation also faithfully reproduces time-dependent aspects of the fluctuations. The method is accurate even under conditions of poor timescale separation.     

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.     

Stochastic simulations on a model of circadian rhythm generation

Biological phenomena are often modeled by differential equations, where states of a model system are described by continuous real values. When we consider concentrations of molecules as dynamical variables for a set of biochemical reactions, we implicitly assume that numbers of the molecules are large enough so that their changes can be regarded as continuous and they are described deterministically. However, for a system with small numbers of molecules, changes in their numbers are apparently discrete and molecular noises become significant. In such cases, models with deterministic differential equations may be inappropriate, and the reactions must be described by stochastic equations. In this study, we focus a clock gene expression for a circadian rhythm generation, which is known as a system involving small numbers of molecules. Thus it is appropriate for the system to be modeled by stochastic equations and analyzed by methodologies of stochastic simulations. The interlocked feedback model proposed by Ueda et al. as a set of deterministic ordinary differential equations provides a basis of our analyses. We apply two stochastic simulation methods, namely Gillespie's direct method and the stochastic differential equation method also by Gillespie, to the interlocked feedback model. To this end, we first reformulated the original differential equations back to elementary chemical reactions. With those reactions, we simulate and analyze the dynamics of the model using two methods in order to compare them with the dynamics obtained from the original deterministic model and to characterize dynamics how they depend on the simulation methodologies.     

Timescale analysis of rule-based biochemical reaction networks

The flow of information within a cell is governed by a series of protein–protein interactions that can be described as a reaction network. Mathematical models of biochemical reaction networks can be constructed by repetitively applying specific rules that define how reactants interact and what new species are formed on reaction. To aid in understanding the underlying biochemistry, timescale analysis is one method developed to prune the size of the reaction network. In this work, we extend the methods associated with timescale analysis to reaction rules instead of the species contained within the network. To illustrate this approach, we applied timescale analysis to a simple receptor–ligand binding model and a rule‐based model of interleukin‐12 (IL‐12) signaling in naïve CD4+ T cells. The IL‐12 signaling pathway includes multiple protein–protein interactions that collectively transmit information; however, the level of mechanistic detail sufficient to capture the observed dynamics has not been justified based on the available data. The analysis correctly predicted that reactions associated with Janus Kinase 2 and Tyrosine Kinase 2 binding to their corresponding receptor exist at a pseudo‐equilibrium. By contrast, reactions associated with ligand binding and receptor turnover regulate cellular response to IL‐12. An empirical Bayesian approach was used to estimate the uncertainty in the timescales. This approach complements existing rank‐ and flux‐based methods that can be used to interrogate complex reaction networks. Ultimately, timescale analysis of rule‐based models is a computational tool that can be used to reveal the biochemical steps that regulate signaling dynamics. © 2011 American Institute of Chemical Engineers Biotechnol. Prog., 2012     

Robust permanence and impermanence for stochastic replicator dynamics

Garay and Hofbauer (SIAM J. Math. Anal. 34 (2003)) proposed sufficient conditions for robust permanence and impermanence of the deterministic replicator dynamics. We reconsider these conditions in the context of the stochastic replicator dynamics, which is obtained from its deterministic analogue by introducing Brownian perturbations of payoffs. When the deterministic replicator dynamics is permanent and the noise level small, the stochastic dynamics admits a unique ergodic distribution whose mass is concentrated near the maximal interior attractor of the unperturbed system; thus, permanence is robust against small unbounded stochastic perturbations. When the deterministic dynamics is impermanent and the noise level small or large, the stochastic dynamics converges to the boundary of the state space at an exponential rate.     

Are Quasi-Steady-State Approximated Models Suitable for Quantifying Intrinsic Noise Accurately?

Large gene regulatory networks (GRN) are often modeled with quasi-steady-state approximation (QSSA) to reduce the huge computational time required for intrinsic noise quantification using Gillespie stochastic simulation algorithm (SSA). However, the question still remains whether the stochastic QSSA model measures the intrinsic noise as accurately as the SSA performed for a detailed mechanistic model or not? To address this issue, we have constructed mechanistic and QSSA models for few frequently observed GRNs exhibiting switching behavior and performed stochastic simulations with them. Our results strongly suggest that the performance of a stochastic QSSA model in comparison to SSA performed for a mechanistic model critically relies on the absolute values of the mRNA and protein half-lives involved in the corresponding GRN. The extent of accuracy level achieved by the stochastic QSSA model calculations will depend on the level of bursting frequency generated due to the absolute value of the half-life of either mRNA or protein or for both the species. For the GRNs considered, the stochastic QSSA quantifies the intrinsic noise at the protein level with greater accuracy and for larger combinations of half-life values of mRNA and protein, whereas in case of mRNA the satisfactory accuracy level can only be reached for limited combinations of absolute values of half-lives. Further, we have clearly demonstrated that the abundance levels of mRNA and protein hardly matter for such comparison between QSSA and mechanistic models. Based on our findings, we conclude that QSSA model can be a good choice for evaluating intrinsic noise for other GRNs as well, provided we make a rational choice based on experimental half-life values available in literature.     

Resonance of plankton communities with temperature fluctuations

The interplay between intrinsic population dynamics and environmental variation is still poorly understood. It is known, however, that even mild environmental noise may induce large fluctuations in population abundances. This is due to a resonance effect that occurs in communities on the edge of stability. Here, we use a simple predator-prey model to explore the sensitivity of plankton communities to stochastic environmental fluctuations. Our results show that the magnitude of resonance depends on the timescale of intrinsic population dynamics relative to the characteristic timescale of the environmental fluctuations. Predator-prey communities with an intrinsic tendency to oscillate at a period T are particularly responsive to red noise characterized by a timescale of τ = T/2π. We compare these theoretical predictions with the timescales of temperature fluctuations measured in lakes and oceans. This reveals that plankton communities will be highly sensitive to natural temperature fluctuations. More specifically, we demonstrate that the relatively fast temperature fluctuations in shallow lakes fall largely within the range to which rotifers and cladocerans are most sensitive, while marine copepods and krill will tend to resonate more strongly with the slower temperature variability of the open ocean.     

Macroscopic Equations Governing Noisy Spiking Neuronal Populations with Linear Synapses

Deriving tractable reduced equations of biological neural networks capturing the macroscopic dynamics of sub-populations of neurons has been a longstanding problem in computational neuroscience. In this paper, we propose a reduction of large-scale multi-population stochastic networks based on the mean-field theory. We derive, for a wide class of spiking neuron models, a system of differential equations of the type of the usual Wilson-Cowan systems describing the macroscopic activity of populations, under the assumption that synaptic integration is linear with random coefficients. Our reduction involves one unknown function, the effective non-linearity of the network of populations, which can be analytically determined in simple cases, and numerically computed in general. This function depends on the underlying properties of the cells, and in particular the noise level. Appropriate parameters and functions involved in the reduction are given for different models of neurons: McKean, Fitzhugh-Nagumo and Hodgkin-Huxley models. Simulations of the reduced model show a precise agreement with the macroscopic dynamics of the networks for the first two models.     

Noise-reduction through interaction in gene expression and biochemical reaction processes

We demonstrate that interaction in gene expression and biochemical reaction processes has a significant influence on reducing fluctuations. Especially, we have found that the interaction between synthesized proteins and background molecules can reduce the fluctuation level in gene expression, which is a counter example to the intuition that background factors disturb information processing in genetic networks by increasing the noise level. This fact also indicates that the macromolecular crowding observed in actual cells can contribute to reduce the noise level. In addition, the noise-reduction phenomenon is not limited to the interaction between the proteins and background molecules, but can be applied to other reactions such as a dimerization process and the coupling of reactions with large fluctuations by intrinsic noise. Finally, on the basis of these results, we propose a new and plausible method for reducing the fluctuations generated in synthesized genetic networks, and also discuss the applicability of this method to the stabilization of system dynamics by using a toggle switch model.     

Moment fitting for parameter inference in repeatedly and partially observed stochastic biological models

The inference of reaction rate parameters in biochemical network models from time series concentration data is a central task in computational systems biology. Under the assumption of well mixed conditions the network dynamics are typically described by the chemical master equation, the Fokker Planck equation, the linear noise approximation or the macroscopic rate equation. The inverse problem of estimating the parameters of the underlying network model can be approached in deterministic and stochastic ways, and available methods often compare individual or mean concentration traces obtained from experiments with theoretical model predictions when maximizing likelihoods, minimizing regularized least squares functionals, approximating posterior distributions or sequentially processing the data. In this article we assume that the biological reaction network can be observed at least partially and repeatedly over time such that sample moments of species molecule numbers for various time points can be calculated from the data. Based on the chemical master equation we furthermore derive closed systems of parameter dependent nonlinear ordinary differential equations that predict the time evolution of the statistical moments. For inferring the reaction rate parameters we suggest to not only compare the sample mean with the theoretical mean prediction but also to take the residual of higher order moments explicitly into account. Cost functions that involve residuals of higher order moments may form landscapes in the parameter space that have more pronounced curvatures at the minimizer and hence may weaken or even overcome parameter sloppiness and uncertainty. As a consequence both deterministic and stochastic parameter inference algorithms may be improved with respect to accuracy and efficiency. We demonstrate the potential of moment fitting for parameter inference by means of illustrative stochastic biological models from the literature and address topics for future research.     

Synchronized dynamics and non-equilibrium steady states in a stochastic yeast cell-cycle network

Applying the mathematical circulation theory of Markov chains, we investigate the synchronized stochastic dynamics of a discrete network model of yeast cell-cycle regulation where stochasticity has been kept rather than being averaged out. By comparing the network dynamics of the stochastic model with its corresponding deterministic network counterpart, we show that the synchronized dynamics can be soundly characterized by a dominant circulation in the stochastic model, which is the natural generalization of the deterministic limit cycle in the deterministic system. Moreover, the period of the main peak in the power spectrum, which is in common use to characterize the synchronized dynamics, perfectly corresponds to the number of states in the main cycle with dominant circulation. Such a large separation in the magnitude of the circulations, between a dominant, main cycle and the rest, gives rise to the stochastic synchronization phenomenon.     

A method for estimating stochastic noise in large genetic regulatory networks

MOTIVATION: Genetic regulatory networks are often affected by stochastic noise, due to the low number of molecules taking part in certain reactions. The networks can be simulated using stochastic techniques that model each reaction as a stochastic event. As models become increasingly large and sophisticated, however, the solution time can become excessive; particularly if one wishes to determine the effect on noise of changes to a series of parameters, or the model structure. Methods are therefore required to rapidly estimate stochastic noise. RESULTS: This paper presents an algorithm, based on error growth techniques from non-linear dynamics, to rapidly estimate the noise characteristics of genetic networks of arbitrary size. The method can also be used to determine analytical solutions for simple sub-systems. It is demonstrated on a number of cases, including a prototype model of the galactose regulatory pathway in yeast. AVAILABILITY: A software tool which incorporates the algorithm is available for use as part of the stochastic simulation package Dizzy. It is available for download at http://labs.systemsbiology.net/bolouri/software/Dizzy/ CONTACT: dorrell@systemsbiology.org SUPPLEMENTARY INFORMATION: A conceptual model of the regulatory part of the galactose utilization pathway in yeast, used as an example in the paper, is available at http://labs.systemsbiology.net/bolouri/models/galconcept.dizzy     

Robust filtering circuit design for stochastic gene networks under intrinsic and extrinsic molecular noises

How to design a robust gene network to tolerate more intrinsic kinetic parameter variations and to attenuate more extrinsic environmental noises to achieve a desired filtering level will be an important topic for systems biology and synthetic biology. At present, there is no good systematic design method to achieve robust gene network design. In this study, a gene network suffering from intrinsic kinetic parameter fluctuations and extrinsic environmental noises is modeled as a Langevin equation with state-dependent stochastic noises. Based on the nonlinear stochastic filtering theory, a systematic gene circuit design method is proposed to make gene networks improve their robustness to tolerate more intrinsic noises and to attenuate extrinsic noises to a prescribed filtering level. The robust gene network design principles have not only yielded a comprehensive design theory of robust gene networks, but also gained valuable insights into the molecular noise filtering of gene networks from the systematic perspective.     

In-phase and anti-phase synchronization in noisy Hodgkin–Huxley neurons

We numerically investigate the influence of intrinsic channel noise on the dynamical response of delay-coupling in neuronal systems. The stochastic dynamics of the spiking is modeled within a stochastic modification of the standard Hodgkin–Huxley model wherein the delay-coupling accounts for the finite propagation time of an action potential along the neuronal axon. We quantify this delay-coupling of the Pyragas-type in terms of the difference between corresponding presynaptic and postsynaptic membrane potentials. For an elementary neuronal network consisting of two coupled neurons we detect characteristic stochastic synchronization patterns which exhibit multiple phase-flip bifurcations: The phase-flip bifurcations occur in form of alternate transitions from an in-phase spiking activity towards an anti-phase spiking activity. Interestingly, these phase-flips remain robust for strong channel noise and in turn cause a striking stabilization of the spiking frequency.     

The development of a noisy brain

Early in life, brain development carries with it a large number of structural changes that impact the functional interactions of distributed neuronal networks. Such changes enhance information processing capacity, moving the brain from a deterministic system to one that is more stochastic. The evidence from empirical studies with EEG and functional MRI suggests that this stochastic property is a result of an increased number of possible functional network configurations for a given situation. This is captured in the variability of endogenous and evoked responses or "brain noise ". In empirical data from infants and children, brain noise increases with maturation and correlates positively with stable behavior and accuracy. The noise increase is best explained through increased noise from network level interactions with a concomitant decrease of local noise. In old adults, brain noise continues to change, although the pattern of changes is not as global as in early development. The relation between high brain noise and stable behavior is maintained, but the relationships differ by region, suggesting changes in local dynamics that then impact potential network configurations. These data, when considered in concert with our extant modeling work, suggest that maturational changes in brain noise represent the enhancement offunctional network potential--the brain's dynamic repertoire.     

The Interplay of Intrinsic and Extrinsic Bounded Noises in Biomolecular Networks

After being considered as a nuisance to be filtered out, it became recently clear that biochemical noise plays a complex role, often fully functional, for a biomolecular network. The influence of intrinsic and extrinsic noises on biomolecular networks has intensively been investigated in last ten years, though contributions on the co-presence of both are sparse. Extrinsic noise is usually modeled as an unbounded white or colored gaussian stochastic process, even though realistic stochastic perturbations are clearly bounded. In this paper we consider Gillespie-like stochastic models of nonlinear networks, i.e. the intrinsic noise, where the model jump rates are affected by colored bounded extrinsic noises synthesized by a suitable biochemical state-dependent Langevin system. These systems are described by a master equation, and a simulation algorithm to analyze them is derived. This new modeling paradigm should enlarge the class of systems amenable at modeling. We investigated the influence of both amplitude and autocorrelation time of a extrinsic Sine-Wiener noise on: the Michaelis-Menten approximation of noisy enzymatic reactions, which we show to be applicable also in co-presence of both intrinsic and extrinsic noise, a model of enzymatic futile cycle and a genetic toggle switch. In and we show that the presence of a bounded extrinsic noise induces qualitative modifications in the probability densities of the involved chemicals, where new modes emerge, thus suggesting the possible functional role of bounded noises.