Discrete-time stochastic modeling and simulation of biochemical networks

Since inherent randomness in chemically reacting systems is evident, stochastic modeling and simulation are exceedingly important for investigating complex biological networks. Within the most common stochastic approach a network is modeled by a continuous-time Markov chain governed by the chemical master equation. We show how the continuous-time Markov chain can be converted to a stochastically identical discrete-time Markov chain and obtain a discrete-time version of the chemical master equation. Simulating the discrete-time Markov chain is equivalent to the Gillespie algorithm but requires less effort in that it eliminates the generation of exponential random variables. Thus, exactness as possessed by the Gillespie algorithm is preserved while the simulation can be performed more efficiently.     

Hybrid stochastic simulations of intracellular reaction–diffusion systems

With the observation that stochasticity is important in biological systems, chemical kinetics have begun to receive wider interest. While the use of Monte Carlo discrete event simulations most accurately capture the variability of molecular species, they become computationally costly for complex reaction–diffusion systems with large populations of molecules. On the other hand, continuous time models are computationally efficient but they fail to capture any variability in the molecular species. In this study a hybrid stochastic approach is introduced for simulating reaction–diffusion systems. We developed an adaptive partitioning strategy in which processes with high frequency are simulated with deterministic rate-based equations, and those with low frequency using the exact stochastic algorithm of Gillespie. Therefore the stochastic behavior of cellular pathways is preserved while being able to apply it to large populations of molecules. We describe our method and demonstrate its accuracy and efficiency compared with the Gillespie algorithm for two different systems. First, a model of intracellular viral kinetics with two steady states and second, a compartmental model of the postsynaptic spine head for studying the dynamics of Ca+2 and NMDA receptors.     

Two classes of quasi-steady-state model reductions for stochastic kinetics

The quasi-steady-state approximation (QSSA) is a model reduction technique used to remove highly reactive species from deterministic models of reaction mechanisms. In many reaction networks the highly reactive intermediates (QSSA species) have populations small enough to require a stochastic representation. In this work we apply singular perturbation analysis to remove the QSSA species from the chemical master equation for two classes of problems. The first class occurs in reaction networks where all the species have small populations and the QSSA species sample zero the majority of the time. The perturbation analysis provides a reduced master equation in which the highly reactive species can sample only zero, and are effectively removed from the model. The reduced master equation can be sampled with the Gillespie algorithm. This first stochastic QSSA reduction is applied to several example reaction mechanisms (including Michaelis-Menten kinetics) [Biochem. Z. 49, 333 (1913)]. A general framework for applying the first QSSA reduction technique to new reaction mechanisms is derived. The second class of QSSA model reductions is derived for reaction networks where non-QSSA species have large populations and QSSA species numbers are small and stochastic. We derive this second QSSA reduction from a combination of singular perturbation analysis and the Omega expansion. In some cases the reduced mechanisms and reaction rates from these two stochastic QSSA models and the classical deterministic QSSA reduction are equivalent; however, this is not usually the case.     

Lattice microbes: High‐performance stochastic simulation method for the reaction‐diffusion master equation

Spatial stochastic simulation is a valuable technique for studying reactions in biological systems. With the availability of high‐performance computing (HPC), the method is poised to allow integration of data from structural, single‐molecule and biochemical studies into coherent computational models of cells. Here, we introduce the Lattice Microbes software package for simulating such cell models on HPC systems. The software performs either well‐stirred or spatially resolved stochastic simulations with approximated cytoplasmic crowding in a fast and efficient manner. Our new algorithm efficiently samples the reaction‐diffusion master equation using NVIDIA graphics processing units and is shown to be two orders of magnitude faster than exact sampling for large systems while maintaining an accuracy of ～ 0.1%. Display of cell models and animation of reaction trajectories involving millions of molecules is facilitated using a plug‐in to the popular VMD visualization platform. The Lattice Microbes software is open source and available for download at http://www.scs.illinois.edu/schulten/lm © 2012 Wiley Periodicals, Inc.     

An equation-free probabilistic steady-state approximation: dynamic application to the stochastic simulation of biochemical reaction networks

Stochastic chemical kinetics more accurately describes the dynamics of "small" chemical systems, such as biological cells. Many real systems contain dynamical stiffness, which causes the exact stochastic simulation algorithm or other kinetic Monte Carlo methods to spend the majority of their time executing frequently occurring reaction events. Previous methods have successfully applied a type of probabilistic steady-state approximation by deriving an evolution equation, such as the chemical master equation, for the relaxed fast dynamics and using the solution of that equation to determine the slow dynamics. However, because the solution of the chemical master equation is limited to small, carefully selected, or linear reaction networks, an alternate equation-free method would be highly useful. We present a probabilistic steady-state approximation that separates the time scales of an arbitrary reaction network, detects the convergence of a marginal distribution to a quasi-steady-state, directly samples the underlying distribution, and uses those samples to accurately predict the state of the system, including the effects of the slow dynamics, at future times. The numerical method produces an accurate solution of both the fast and slow reaction dynamics while, for stiff systems, reducing the computational time by orders of magnitude. The developed theory makes no approximations on the shape or form of the underlying steady-state distribution and only assumes that it is ergodic. We demonstrate the accuracy and efficiency of the method using multiple interesting examples, including a highly nonlinear protein-protein interaction network. The developed theory may be applied to any type of kinetic Monte Carlo simulation to more efficiently simulate dynamically stiff systems, including existing exact, approximate, or hybrid stochastic simulation techniques.     

Developing Itô stochastic differential equation models for neuronal signal transduction pathways

Mathematical modeling and simulation of dynamic biochemical systems are receiving considerable attention due to the increasing availability of experimental knowledge of complex intracellular functions. In addition to deterministic approaches, several stochastic approaches have been developed for simulating the time-series behavior of biochemical systems. The problem with stochastic approaches, however, is the larger computational time compared to deterministic approaches. It is therefore necessary to study alternative ways to incorporate stochasticity and to seek approaches that reduce the computational time needed for simulations, yet preserve the characteristic behavior of the system in question. In this work, we develop a computational framework based on the It? stochastic differential equations for neuronal signal transduction networks. There are several different ways to incorporate stochasticity into deterministic differential equation models and to obtain It? stochastic differential equations. Two of the developed models are found most suitable for stochastic modeling of neuronal signal transduction. The best models give stable responses which means that the variances of the responses with time are not increasing and negative concentrations are avoided. We also make a comparative analysis of different kinds of stochastic approaches, that is the It? stochastic differential equations, the chemical Langevin equation, and the Gillespie stochastic simulation algorithm. Different kinds of stochastic approaches can be used to produce similar responses for the neuronal protein kinase C signal transduction pathway. The fine details of the responses vary slightly, depending on the approach and the parameter values. However, when simulating great numbers of chemical species, the Gillespie algorithm is computationally several orders of magnitude slower than the It? stochastic differential equations and the chemical Langevin equation. Furthermore, the chemical Langevin equation produces negative concentrations. The It? stochastic differential equations developed in this work are shown to overcome the problem of obtaining negative values.     

Efficient stochastic simulation of simultaneous reaction and diffusion in a gas-liquid interface

A stochastic simulation of simultaneous reaction and diffusion is proposed for the gas-liquid interface formed in the surface of a gas bubble within a liquid. The interface between a carbon dioxide bubble and an aqueous solution of calcium hydroxide was simulated as an application example, taken from the integrated production of calcium carbonate. First Gillespie’s stochastic simulation algorithm was applied in separate reaction and diffusion simulations. The results from these simulations were consistent with deterministic solutions based on differential equations. However it was observed that stochastic diffusion simulations are extremely slow. The sampling of diffusion events was accelerated applying a group molecule transfer scheme based on the binomial distribution function. Simulations of the reaction-diffusion in the gas-liquid interface based on the standard Gillespie’s stochastic algorithm were also slow. However the application of the binomial distribution function scheme allowed to compute the concentration profiles in the gas-liquid interface in a fraction of the time required with the standard Gillespie’s stochastic algorithm.     

On the origins of approximations for stochastic chemical kinetics

This paper considers the derivation of approximations for stochastic chemical kinetics governed by the discrete master equation. Here, the concepts of (1) partitioning on the basis of fast and slow reactions as opposed to fast and slow species and (2) conditional probability densities are used to derive approximate, partitioned master equations, which are Markovian in nature, from the original master equation. Under different conditions dictated by relaxation time arguments, such approximations give rise to both the equilibrium and hybrid (deterministic or Langevin equations coupled with discrete stochastic simulation) approximations previously reported. In addition, the derivation points out several weaknesses in previous justifications of both the hybrid and equilibrium systems and demonstrates the connection between the original and approximate master equations. Two simple examples illustrate situations in which these two approximate methods are applicable and demonstrate the two methods' efficiencies.     

A reduction method for multiple time scale stochastic reaction networks

In this paper we develop a reduction method for multiple time scale stochastic reaction networks. When the transition-rate matrix between different states of the species is available, we obtain systems of reduced equations, whose solutions can successively approximate, to any degree of accuracy, the exact probability that the reaction system be in any particular state. For the case when the transition-rate matrix is not available, one needs to rely on the chemical master equation. For this case, we obtain a corresponding reduced master equation with first-order accuracy. We illustrate the accuracy and efficiency of both approaches by simulating several motivating examples and comparing the results of our simulations with the results obtained by the exact method. Our examples include both linear and nonlinear reaction networks as well as a three time scale stochastic reaction-diffusion model arising from gene expression.     

Stochastic hybrid modeling of intracellular calcium dynamics

Deterministic models of biochemical processes at the subcellular level might become inadequate when a cascade of chemical reactions is induced by a few molecules. Inherent randomness of such phenomena calls for the use of stochastic simulations. However, being computationally intensive, such simulations become infeasible for large and complex reaction networks. To improve their computational efficiency in handling these networks, we present a hybrid approach, in which slow reactions and fluxes are handled through exact stochastic simulation and their fast counterparts are treated partially deterministically through chemical Langevin equation. The classification of reactions as fast or slow is accompanied by the assumption that in the time-scale of fast reactions, slow reactions do not occur and hence do not affect the probability of the state. Our new approach also handles reactions with complex rate expressions such as Michaelis-Menten kinetics. Fluxes which cannot be modeled explicitly through reactions, such as flux of Ca(2+) from endoplasmic reticulum to the cytosol through inositol 1,4,5-trisphosphate receptor channels, are handled deterministically. The proposed hybrid algorithm is used to model the regulation of the dynamics of cytosolic calcium ions in mouse macrophage RAW 264.7 cells. At relatively large number of molecules, the response characteristics obtained with the stochastic and deterministic simulations coincide, which validates our approach in the limit of large numbers. At low doses, the response characteristics of some key chemical species, such as levels of cytosolic calcium, predicted with stochastic simulations, differ quantitatively from their deterministic counterparts. These observations are ubiquitous throughout dose response, sensitivity, and gene-knockdown response analyses. While the relative differences between the peak-heights of the cytosolic [Ca(2+)] time-courses obtained from stochastic (mean of 16 realizations) and deterministic simulations are merely 1%-4% for most perturbations, it is specially sensitive to levels of G(βγ) (relative difference as large as 90% at very low G(βγ)).     

受高斯白噪声影响的有限化学反应体系的随机热力学——外噪声对化学反应体系影响的有效热力学量度

通过分析噪声对跃迁概率的不同影响,借助Novikov定理及MSR理论,建立了受多源白噪声影响的有限化学反应体系的有效主方程及有效熵平衡方程,导出非平衡定态时这类噪声体系熵产生的一般表达式,揭示涨落熵产生的统计内涵及噪声贡献,并针对非宏观量级的外噪声,借助扰动按分布参数分离法及有效主方程的Kramers-Moyal展开,进一步对简单加合性噪声建立了非平衡定态宏观稳定性判据的随机模拟,论证了噪声对化学反应体系定态稳定性的弱化作用.     

The effect of memory in the stochastic master equation analyzed using the stochastic Liouville equation of motion. Electronic energy migration transfer between reorienting donor-donor, donor-acceptor chromophores

This paper discusses the process of energy migration transfer within reorientating chromophores using the stochastic master equation (SME) and the stochastic Liouville equation (SLE) of motion. We have found that the SME over-estimates the rate of the energy migration compared to the SLE solution for a case of weakly interacting chromophores. This discrepancy between SME and SLE is caused by a memory effect occurring when fluctuations in the dipole-dipole Hamiltonian (H(t)) are on the same timescale as the intrinsic fast transverse relaxation rate characterized by (1/T(2)). Thus the timescale critical for energy-transfer experiments is T(2) approximately 10(-13) s. An extended SME is constructed, accounting for the memory effect of the dipole-dipole Hamiltonian dynamics. The influence of memory on the interpretation of experiments is discussed.     

Spatiotemporal laser inactivation of inositol 1,4,5-trisphosphate receptors using synthetic small-molecule probes

A malachite green-conjugated inositol 1,4,5-trisphosphate (MGIP(3)) induces specific inactivation of IP(3) receptor (IP(3)R) in tissue samples upon laser irradiation. To verify potential usefulness of the method for studies of cellular Ca(2+) signaling, we conducted laser inactivation at the single-cell level and show that IP(3)R was inactivated with extremely high spatiotemporal resolution. In the presence of MGIP(3), the Ca(2+) release function of IP(3)R in single B lymphoma cells decayed exponentially with increasing duration of laser irradiation with a time constant of 3.4 s. Moreover, by confining laser irradiation to a spatially distinct region of differentiated PC12 cells, subcellular inactivation of IP(3)R was attained, as revealed by a loss of local Ca(2+) signal. Such real-time inactivation of IP(3)R only within a subcellular region may provide a powerful method for investigating spatiotemporal dynamics of Ca(2+) signaling.     

Optical control of calcium affinity in a spiroamido-rhodamine based calcium chelator

An optically controlled Ca(2+)-chelator 1 was developed to mimic natural calcium oscillations. Compound 1, a spiroamido-rhodamine derivative of 1,2-bis(o-aminophenoxy)ethane-N,N,N',N'-tetraacetic acid (BAPTA), underwent cycles of reversible transitions between a colorless closed state and a fluorescent open form. The closed-state exhibited a high affinity for Ca(2+) (K(d): 509 nM) with excellent selectivity over Mg(2+) (K(d): 19 mM). The open isomer had a 350-fold lower Ca(2+) affinity (K(d): 181 μM), while the Mg(2+) affinity was not significantly affected (K(d): 14 mM).     

Nested stochastic simulation algorithm for chemical kinetic systems with disparate rates

An efficient simulation algorithm for chemical kinetic systems with disparate rates is proposed. This new algorithm is quite general, and it amounts to a simple and seamless modification of the classical stochastic simulation algorithm (SSA), also known as the Gillespie [J. Comput. Phys. 22, 403 (1976); J. Phys. Chem. 81, 2340 (1977)] algorithm. The basic idea is to use an outer SSA to simulate the slow processes with rates computed from an inner SSA which simulates the fast reactions. Averaging theorems for Markov processes can be used to identify the fast and slow variables in the system as well as the effective dynamics over the slow time scale, even though the algorithm itself does not rely on such information. This nested SSA can be easily generalized to systems with more than two separated time scales. Convergence and efficiency of the algorithm are discussed using the established error estimates and illustrated through examples.     

"All possible steps" approach to the accelerated use of Gillespie's algorithm

Many physical and biological processes are stochastic in nature. Computational models and simulations of such processes are a mathematical and computational challenge. The basic stochastic simulation algorithm was published by Gillespie about three decades ago [J. Phys. Chem. 81, 2340 (1977)]. Since then, intensive work has been done to make the algorithm more efficient in terms of running time. All accelerated versions of the algorithm are aimed at minimizing the running time required to produce a stochastic trajectory in state space. In these simulations, a necessary condition for reliable statistics is averaging over a large number of simulations. In this study the author presents a new accelerating approach which does not alter the stochastic algorithm, but reduces the number of required runs. By analysis of collected data the author demonstrates high precision levels with fewer simulations. Moreover, the suggested approach provides a good estimation of statistical error, which may serve as a tool for determining the number of required runs.     

非恒稳含时电极过程中涨落-耗散的电化学效应——恒电势滴汞电极体系的随机热力学模型

以滴汞电极体系为模型,对非恒稳恒电势动态不可逆电极过程中的耗散-涨落效应进行了系统的研究.基于滴汞电极体系的电化学特征,提出了一个简化的含时随机热力学模型,从而可能对这类重要的含时物理化学过程进行涨落和耗散的定量分析.借助该简化的模型,成功地建立了恒电势滴汞电极过程的基本随机热力学公式,由此推出耗散-涨落效应的理论极谱曲线.在滴汞电极生长缓慢及扩散步骤严重滞后情况下,含时的滴汞电极过程将趋于在有效扩散层厚度演化的慢流型上的准定态过程.在这种准定态近似下,具体分析了涨落对极谱曲线的影响.结果表明,在涨落影响可以忽略的近平衡区,从耗散-涨落导出的极谱方程与从平衡态Nernst公式导出的极谱方程完全吻合.还计算了一个涨落诱导的极谱曲线偏离的典型范例.     

Transient, sparsely populated compact states of apo and calcium-loaded calmodulin probed by paramagnetic relaxation enhancement: interplay of conformational selection and induced fit

Calmodulin (CaM) is the universal calcium sensor in eukaryotes, regulating the function of numerous proteins. Crystallography and NMR show that free CaM-4Ca(2+) exists in an extended conformation with significant interdomain separation, but clamps down upon target peptides to form a highly compact structure. NMR has revealed substantial interdomain motions in CaM-4Ca(2+), enabled by a flexible linker. In one instance, CaM-4Ca(2+) has been crystallized in a compact configuration; however, no direct evidence for transient interdomain contacts has been observed in solution, and little is known about how large-scale interdomain motions contribute to biological function. Here, we use paramagnetic relaxation enhancement (PRE) to characterize transient compact states of free CaM that are too sparsely populated to observe by traditional NMR methods. We show that unbound CaM samples a range of compact structures, populated at 5-10%, and that Ca(2+) dramatically alters the distribution of these configurations in favor of states resembling the peptide-bound structure. In the absence of Ca(2+), the target peptide binds only to the C-terminal domain, and the distribution of compact states is similar with and without peptide. These data suggest an alternative pathway of CaM action in which CaM remains associated with its kinase targets even in the resting state. Only CaM-4Ca(2+), however, shows an innate propensity to form the physiologically active compact structures, suggesting that Ca(2+) activates CaM not only through local structural changes within each domain but also through more global remodeling of interdomain interactions. Thus, these findings illustrate the subtle interplay between conformational selection and induced fit.     

非平衡定态化学反应体系对细致平衡偏离的不可逆性之随机热力学量度

在随机描述的层次上,成功地构筑了非平衡定态化学反应体系对细致平衡偏离的不可逆性之随机热力学判据.基于Fokker-Planck方程建立了连续Markov过程系统的随机熵平衡方程,发现在随机态空间中随机熵的时间变率亦可分为源项和流项两部分贡献.对于化学反应体系,作为态空间源项的随机熵产生可作为偏离细致平衡不可逆性的合适度量,此泛函量的零值标志着细致平衡.进一步借助按系统广度量的Ω展开法,通过对随机力及共轭的随机流的分析揭示态空间中的随机熵产生仅源于状态对细致平衡的偏离,并主要来自非Poisson涨落的贡献,因而可作为随机热力学量去判别和量度化学反应体系对细致平衡的偏离.