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 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.     

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.     

Computational methods for diffusion-influenced biochemical reactions

MOTIVATION: We compare stochastic computational methods accounting for space and discrete nature of reactants in biochemical systems. Implementations based on Brownian dynamics (BD) and the reaction-diffusion master equation are applied to a simplified gene expression model and to a signal transduction pathway in Escherichia coli. RESULTS: In the regime where the number of molecules is small and reactions are diffusion-limited predicted fluctuations in the product number vary between the methods, while the average is the same. Computational approaches at the level of the reaction-diffusion master equation compute the same fluctuations as the reference result obtained from the particle-based method if the size of the sub-volumes is comparable to the diameter of reactants. Using numerical simulations of reversible binding of a pair of molecules we argue that the disagreement in predicted fluctuations is due to different modeling of inter-arrival times between reaction events. Simulations for a more complex biological study show that the different approaches lead to different results due to modeling issues. Finally, we present the physical assumptions behind the mesoscopic models for the reaction-diffusion systems. AVAILABILITY: Input files for the simulations and the source code of GMP can be found under the following address: http://www.cwi.nl/projects/sic/bioinformatics2007/     

An approximate stochastic optimal control framework to simulate nonlinear neuro-musculoskeletal models in the presence of noise

Optimal control simulations have shown that both musculoskeletal dynamics and physiological noise are important determinants of movement. However, due to the limited efficiency of available computational tools, deterministic simulations of movement focus on accurately modelling the musculoskeletal system while neglecting physiological noise, and stochastic simulations account for noise while simplifying the dynamics. We took advantage of recent approaches where stochastic optimal control problems are approximated using deterministic optimal control problems, which can be solved efficiently using direct collocation. We were thus able to extend predictions of stochastic optimal control as a theory of motor coordination to include muscle coordination and movement patterns emerging from non-linear musculoskeletal dynamics. In stochastic optimal control simulations of human standing balance, we demonstrated that the inclusion of muscle dynamics can predict muscle co-contraction as minimal effort strategy that complements sensorimotor feedback control in the presence of sensory noise. In simulations of reaching, we demonstrated that nonlinear multi-segment musculoskeletal dynamics enables complex perturbed and unperturbed reach trajectories under a variety of task conditions to be predicted. In both behaviors, we demonstrated how interactions between task constraint, sensory noise, and the intrinsic properties of muscle influence optimal muscle coordination patterns, including muscle co-contraction, and the resulting movement trajectories. Our approach enables a true minimum effort solution to be identified as task constraints, such as movement accuracy, can be explicitly imposed, rather than being approximated using penalty terms in the cost function. Our approximate stochastic optimal control framework predicts complex features, not captured by previous simulation approaches, providing a generalizable and valuable tool to study how musculoskeletal dynamics and physiological noise may alter neural control of movement in both healthy and pathological movements.     

On the validity of the quasi-steady state approximation of bimolecular reactions in solution

Two-step binding kinetics are extensively used to study the relative importance of diffusion in biochemical reactions. Classical analysis of this problem assumes ad hoc that the encounter complex is at quasi-steady state (QSS). Using scaling arguments we derive a criterion for the validity of this assumption in the limit of irreversible product formation. We find that the QSS approximation (QSSA) of two-step binding is only valid if the total ligand and receptor concentrations are much smaller than (k2+k-1)/k1, where k1 and k-1 are, respectively, the forward and reverse diffusion encounter rate constants and k2 is the chemical association rate constant. This criterion can be shown to imply that the average time between encounters is much longer than the half-life of the encounter complex and also guarantees that the concentration of the encounter complex is negligible compared to the reactant and product concentrations. Numerical examples of irreversible and reversible cases corroborate our analysis and illustrate that the QSS may be invalid even if k-2相似文献   

Robustness Analysis of Stochastic Biochemical Systems

We propose a new framework for rigorous robustness analysis of stochastic biochemical systems that is based on probabilistic model checking techniques. We adapt the general definition of robustness introduced by Kitano to the class of stochastic systems modelled as continuous time Markov Chains in order to extensively analyse and compare robustness of biological models with uncertain parameters. The framework utilises novel computational methods that enable to effectively evaluate the robustness of models with respect to quantitative temporal properties and parameters such as reaction rate constants and initial conditions. We have applied the framework to gene regulation as an example of a central biological mechanism where intrinsic and extrinsic stochasticity plays crucial role due to low numbers of DNA and RNA molecules. Using our methods we have obtained a comprehensive and precise analysis of stochastic dynamics under parameter uncertainty. Furthermore, we apply our framework to compare several variants of two-component signalling networks from the perspective of robustness with respect to intrinsic noise caused by low populations of signalling components. We have successfully extended previous studies performed on deterministic models (ODE) and showed that stochasticity may significantly affect obtained predictions. Our case studies demonstrate that the framework can provide deeper insight into the role of key parameters in maintaining the system functionality and thus it significantly contributes to formal methods in computational systems biology.     

Size-Independent Differences between the Mean of Discrete Stochastic Systems and the Corresponding Continuous Deterministic Systems

In this paper, it is shown that for a class of reaction networks, the discrete stochastic nature of the reacting species and reactions results in qualitative and quantitative differences between the mean of exact stochastic simulations and the prediction of the corresponding deterministic system. The differences are independent of the number of molecules of each species in the system under consideration. These reaction networks are open systems of chemical reactions with no zero-order reaction rates. They are characterized by at least two stationary points, one of which is a nonzero stable point, and one unstable trivial solution (stability based on a linear stability analysis of the deterministic system). Starting from a nonzero initial condition, the deterministic system never reaches the zero stationary point due to its unstable nature. In contrast, the result presented here proves that this zero-state is a stable stationary state for the discrete stochastic system, and other finite states have zero probability of existence at large times. This result generalizes previous theoretical studies and simulations of specific systems and provides a theoretical basis for analyzing a class of systems that exhibit such inconsistent behavior. This result has implications in the simulation of infection, apoptosis, and population kinetics, as it can be shown that for certain models the stochastic simulations will always yield different predictions for the mean behavior than the deterministic simulations.     

Stochastic Sensitivity Analysis and Kernel Inference via Distributional Data

Cellular processes are noisy due to the stochastic nature of biochemical reactions. As such, it is impossible to predict the exact quantity of a molecule or other attributes at the single-cell level. However, the distribution of a molecule over a population is often deterministic and is governed by the underlying regulatory networks relevant to the cellular functionality of interest. Recent studies have started to exploit this property to infer network states. To facilitate the analysis of distributional data in a general experimental setting, we introduce a computational framework to efficiently characterize the sensitivity of distributional output to changes in external stimuli. Further, we establish a probability-divergence-based kernel regression model to accurately infer signal level based on distribution measurements. Our methodology is applicable to any biological system subject to stochastic dynamics and can be used to elucidate how population-based information processing may contribute to organism-level functionality. It also lays the foundation for engineering synthetic biological systems that exploit population decoding to more robustly perform various biocomputation tasks, such as disease diagnostics and environmental-pollutant sensing.     

Stochastic Sensitivity Analysis and Kernel Inference via Distributional Data

Cellular processes are noisy due to the stochastic nature of biochemical reactions. As such, it is impossible to predict the exact quantity of a molecule or other attributes at the single-cell level. However, the distribution of a molecule over a population is often deterministic and is governed by the underlying regulatory networks relevant to the cellular functionality of interest. Recent studies have started to exploit this property to infer network states. To facilitate the analysis of distributional data in a general experimental setting, we introduce a computational framework to efficiently characterize the sensitivity of distributional output to changes in external stimuli. Further, we establish a probability-divergence-based kernel regression model to accurately infer signal level based on distribution measurements. Our methodology is applicable to any biological system subject to stochastic dynamics and can be used to elucidate how population-based information processing may contribute to organism-level functionality. It also lays the foundation for engineering synthetic biological systems that exploit population decoding to more robustly perform various biocomputation tasks, such as disease diagnostics and environmental-pollutant sensing.     

Order Reduction of the Chemical Master Equation via Balanced Realisation

We consider a Markov process in continuous time with a finite number of discrete states. The time-dependent probabilities of being in any state of the Markov chain are governed by a set of ordinary differential equations, whose dimension might be large even for trivial systems. Here, we derive a reduced ODE set that accurately approximates the probabilities of subspaces of interest with a known error bound. Our methodology is based on model reduction by balanced truncation and can be considerably more computationally efficient than solving the chemical master equation directly. We show the applicability of our method by analysing stochastic chemical reactions. First, we obtain a reduced order model for the infinitesimal generator of a Markov chain that models a reversible, monomolecular reaction. Later, we obtain a reduced order model for a catalytic conversion of substrate to a product (a so-called Michaelis-Menten mechanism), and compare its dynamics with a rapid equilibrium approximation method. For this example, we highlight the savings on the computational load obtained by means of the reduced-order model. Furthermore, we revisit the substrate catalytic conversion by obtaining a lower-order model that approximates the probability of having predefined ranges of product molecules. In such an example, we obtain an approximation of the output of a model with 5151 states by a reduced model with 16 states. Finally, we obtain a reduced-order model of the Brusselator.     

Monte Carlo simulations of VEGF binding to cell surface receptors in vitro

The vascular endothelial growth factor (VEGF) family binds multiple endothelial cell surface receptors. Our goal is to build comprehensive models of these interactions for the purpose of simulating angiogenesis. In view of low concentrations of growth factors in vivo and in vitro, stochastic modeling of molecular interactions may be necessary. Here, we compare Monte Carlo simulations of the stochastic binding of VEGF and two of its major receptors on cells in vitro to equivalent deterministic simulations. In the range of typical VEGF concentrations, the stochastic and deterministic models are in agreement. However, we observe significant variability in receptor binding, which may be linked to biological stochastic events, e.g., blood vessel sprout initiation. We study patches of cell surface of varying sizes to investigate spatial integration of the signal by the cell, which impacts directly the variability of binding, and find significant variability up to the single-cell level. Dimerization of VEGF receptors does not significantly alter the variability in ligand binding. A 'sliding window' approach demonstrated no reduction in the variability of binding by temporal integration. The variability is expected to be more prominent in in vivo situations where the number of ligand molecules available for binding is less.     

Revising Berg-Purcell for finite receptor kinetics

From nutrient uptake to chemoreception to synaptic transmission, many systems in cell biology depend on molecules diffusing and binding to membrane receptors. Mathematical analysis of such systems often neglects the fact that receptors process molecules at finite kinetic rates. A key example is the celebrated formula of Berg and Purcell for the rate that cell surface receptors capture extracellular molecules. Indeed, this influential result is only valid if receptors transport molecules through the cell wall at a rate much faster than molecules arrive at receptors. From a mathematical perspective, ignoring receptor kinetics is convenient because it makes the diffusing molecules independent. In contrast, including receptor kinetics introduces correlations between the diffusing molecules because, for example, bound receptors may be temporarily blocked from binding additional molecules. In this work, we present a modeling framework for coupling bulk diffusion to surface receptors with finite kinetic rates. The framework uses boundary homogenization to couple the diffusion equation to nonlinear ordinary differential equations on the boundary. We use this framework to derive an explicit formula for the cellular uptake rate and show that the analysis of Berg and Purcell significantly overestimates uptake in some typical biophysical scenarios. We confirm our analysis by numerical simulations of a many-particle stochastic system.     

Individual-based vs deterministic models for macroparasites: host cycles and extinction

Our understanding of the qualitative dynamics of host-macroparasite systems is mainly based on deterministic models. We study here an individual-based stochastic model that incorporates the same assumptions as the classical deterministic model. Stochastic simulations, using parameter values based on some case studies, preserve many features of the deterministic model, like the average value of the variables and the approximate length of the cycles.An important difference is that, even when deterministic models yield damped oscillations, stochastic simulations yield apparently sustained oscillations. The amplitude of such oscillations may be so large to threaten parasites' persistence.With density-dependence in parasite demographic traits, persistence increases somewhat. Allowing instead for infections from an external parasite reservoir, we found that host extinction may easily occur. However, the extinction probability is almost independent of the level of external infection over a wide intermediate parameter region.     

Dynamic partitioning for hybrid simulation of the bistable HIV-1 transactivation network

MOTIVATION: The stochastic kinetics of a well-mixed chemical system, governed by the chemical Master equation, can be simulated using the exact methods of Gillespie. However, these methods do not scale well as systems become more complex and larger models are built to include reactions with widely varying rates, since the computational burden of simulation increases with the number of reaction events. Continuous models may provide an approximate solution and are computationally less costly, but they fail to capture the stochastic behavior of small populations of macromolecules. RESULTS: In this article we present a hybrid simulation algorithm that dynamically partitions the system into subsets of continuous and discrete reactions, approximates the continuous reactions deterministically as a system of ordinary differential equations (ODE) and uses a Monte Carlo method for generating discrete reaction events according to a time-dependent propensity. Our approach to partitioning is improved such that we dynamically partition the system of reactions, based on a threshold relative to the distribution of propensities in the discrete subset. We have implemented the hybrid algorithm in an extensible framework, utilizing two rigorous ODE solvers to approximate the continuous reactions, and use an example model to illustrate the accuracy and potential speedup of the algorithm when compared with exact stochastic simulation. AVAILABILITY: Software and benchmark models used for this publication can be made available upon request from the authors.