A data-integrated method for analyzing stochastic biochemical networks

Variability and fluctuations among genetically identical cells under uniform experimental conditions stem from the stochastic nature of biochemical reactions. Understanding network function for endogenous biological systems or designing robust synthetic genetic circuits requires accounting for and analyzing this variability. Stochasticity in biological networks is usually represented using a continuous-time discrete-state Markov formalism, where the chemical master equation (CME) and its kinetic Monte Carlo equivalent, the stochastic simulation algorithm (SSA), are used. These two representations are computationally intractable for many realistic biological problems. Fitting parameters in the context of these stochastic models is particularly challenging and has not been accomplished for any but very simple systems. In this work, we propose that moment equations derived from the CME, when treated appropriately in terms of higher order moment contributions, represent a computationally efficient framework for estimating the kinetic rate constants of stochastic network models and subsequent analysis of their dynamics. To do so, we present a practical data-derived moment closure method for these equations. In contrast to previous work, this method does not rely on any assumptions about the shape of the stochastic distributions or a functional relationship among their moments. We use this method to analyze a stochastic model of a biological oscillator and demonstrate its accuracy through excellent agreement with CME/SSA calculations. By coupling this moment-closure method with a parameter search procedure, we further demonstrate how a model's kinetic parameters can be iteratively determined in order to fit measured distribution data.     

State-dependent doubly weighted stochastic simulation algorithm for automatic characterization of stochastic biochemical rare events

In recent years there has been substantial growth in the development of algorithms for characterizing rare events in stochastic biochemical systems. Two such algorithms, the state-dependent weighted stochastic simulation algorithm (swSSA) and the doubly weighted SSA (dwSSA) are extensions of the weighted SSA (wSSA) by H. Kuwahara and I. Mura [J. Chem. Phys. 129, 165101 (2008)]. The swSSA substantially reduces estimator variance by implementing system state-dependent importance sampling (IS) parameters, but lacks an automatic parameter identification strategy. In contrast, the dwSSA provides for the automatic determination of state-independent IS parameters, thus it is inefficient for systems whose states vary widely in time. We present a novel modification of the dwSSA--the state-dependent doubly weighted SSA (sdwSSA)--that combines the strengths of the swSSA and the dwSSA without inheriting their weaknesses. The sdwSSA automatically computes state-dependent IS parameters via the multilevel cross-entropy method. We apply the method to three examples: a reversible isomerization process, a yeast polarization model, and a lac operon model. Our results demonstrate that the sdwSSA offers substantial improvements over previous methods in terms of both accuracy and efficiency.     

Stochastic simulation of chemically reacting systems using multi-core processors

In recent years, computer simulations have become increasingly useful when trying to understand the complex dynamics of biochemical networks, particularly in stochastic systems. In such situations stochastic simulation is vital in gaining an understanding of the inherent stochasticity present, as these models are rarely analytically tractable. However, a stochastic approach can be computationally prohibitive for many models. A number of approximations have been proposed that aim to speed up stochastic simulations. However, the majority of these approaches are fundamentally serial in terms of central processing unit (CPU) usage. In this paper, we propose a novel simulation algorithm that utilises the potential of multi-core machines. This algorithm partitions the model into smaller sub-models. These sub-models are then simulated, in parallel, on separate CPUs. We demonstrate that this method is accurate and can speed-up the simulation by a factor proportional to the number of processors available.     

A constant-time kinetic Monte Carlo algorithm for simulation of large biochemical reaction networks

The time evolution of species concentrations in biochemical reaction networks is often modeled using the stochastic simulation algorithm (SSA) [Gillespie, J. Phys. Chem. 81, 2340 (1977)]. The computational cost of the original SSA scaled linearly with the number of reactions in the network. Gibson and Bruck developed a logarithmic scaling version of the SSA which uses a priority queue or binary tree for more efficient reaction selection [Gibson and Bruck, J. Phys. Chem. A 104, 1876 (2000)]. More generally, this problem is one of dynamic discrete random variate generation which finds many uses in kinetic Monte Carlo and discrete event simulation. We present here a constant-time algorithm, whose cost is independent of the number of reactions, enabled by a slightly more complex underlying data structure. While applicable to kinetic Monte Carlo simulations in general, we describe the algorithm in the context of biochemical simulations and demonstrate its competitive performance on small- and medium-size networks, as well as its superior constant-time performance on very large networks, which are becoming necessary to represent the increasing complexity of biochemical data for pathways that mediate cell function.     

Efficient formulation of the stochastic simulation algorithm for chemically reacting systems

In this paper we examine the different formulations of Gillespie's stochastic simulation algorithm (SSA) [D. Gillespie, J. Phys. Chem. 81, 2340 (1977)] with respect to computational efficiency, and propose an optimization to improve the efficiency of the direct method. Based on careful timing studies and an analysis of the time-consuming operations, we conclude that for most practical problems the optimized direct method is the most efficient formulation of SSA. This is in contrast to the widely held belief that Gibson and Bruck's next reaction method [M. Gibson and J. Bruck, J. Phys. Chem. A 104, 1876 (2000)] is the best way to implement the SSA for large systems. Our analysis explains the source of the discrepancy.     

Abstract Next Subvolume Method: a logical process-based approach for spatial stochastic simulation of chemical reactions

The spatial stochastic simulation of biochemical systems requires significant calculation efforts. Parallel discrete-event simulation is a promising approach to accelerate the execution of simulation runs. However, achievable speedup depends on the parallelism inherent in the model. One of our goals is to explore this degree of parallelism in the Next Subvolume Method type simulations. Therefore we introduce the Abstract Next Subvolume Method, in which we decouple the model representation from the sequential simulation algorithms, and prove that state trajectories generated by its executions statistically accord with those generated by the Next Subvolume Method. The experimental performance analysis shows that optimistic synchronization algorithms, together with careful controls over the speculative execution, are necessary to achieve considerable speedup and scalability in parallel spatial stochastic simulation of chemical reactions. Our proposed method facilitates a flexible incorporation of different synchronization algorithms, and can be used to select the proper synchronization algorithm to achieve the efficient parallel simulation of chemical reactions.     

The finite state projection algorithm for the solution of the chemical master equation

This article introduces the finite state projection (FSP) method for use in the stochastic analysis of chemically reacting systems. One can describe the chemical populations of such systems with probability density vectors that evolve according to a set of linear ordinary differential equations known as the chemical master equation (CME). Unlike Monte Carlo methods such as the stochastic simulation algorithm (SSA) or tau leaping, the FSP directly solves or approximates the solution of the CME. If the CME describes a system that has a finite number of distinct population vectors, the FSP method provides an exact analytical solution. When an infinite or extremely large number of population variations is possible, the state space can be truncated, and the FSP method provides a certificate of accuracy for how closely the truncated space approximation matches the true solution. The proposed FSP algorithm systematically increases the projection space in order to meet prespecified tolerance in the total probability density error. For any system in which a sufficiently accurate FSP exists, the FSP algorithm is shown to converge in a finite number of steps. The FSP is utilized to solve two examples taken from the field of systems biology, and comparisons are made between the FSP, the SSA, and tau leaping algorithms. In both examples, the FSP outperforms the SSA in terms of accuracy as well as computational efficiency. Furthermore, due to very small molecular counts in these particular examples, the FSP also performs far more effectively than tau leaping methods.     

A cutoff phenomenon in accelerated stochastic simulations of chemical kinetics via flow averaging (FLAVOR-SSA)

We present a simple algorithm for the simulation of stiff, discrete-space, continuous-time Markov processes. The algorithm is based on the concept of flow averaging for the integration of stiff ordinary and stochastic differential equations and ultimately leads to a straightforward variation of the the well-known stochastic simulation algorithm (SSA). The speedup that can be achieved by the present algorithm [flow averaging integrator SSA (FLAVOR-SSA)] over the classical SSA comes naturally at the expense of its accuracy. The error of the proposed method exhibits a cutoff phenomenon as a function of its speed-up, allowing for optimal tuning. Two numerical examples from chemical kinetics are provided to illustrate the efficiency of the method.     

K-leap method for accelerating stochastic simulation of coupled chemical reactions

Leap methods are very promising for accelerating stochastic simulation of a well stirred chemically reacting system, while providing acceptable simulation accuracy. In Gillespie's tau-leap method [D. Gillespie, J. Phys. Chem. 115, 1716 (2001)], the number of firings of each reaction channel during a leap is a Poisson random variable, whose sample values are unbounded. This may cause large changes in the populations of certain molecular species during a leap, thereby violating the leap condition. In this paper, we develop an alternative leap method called the K-leap method, in which we constrain the total number of reactions occurring during a leap to be a number K calculated from the leap condition. As the number of firings of each reaction channel during a leap is upper bounded by a properly chosen number, our K-leap method can better satisfy the leap condition, thereby improving simulation accuracy. Since the exact stochastic simulation algorithm (SSA) is a special case of our K-leap method when K=1, our K-leap method can naturally change from the exact SSA to an approximate leap method during simulation, whenever the leap condition allows to do so.     

Efficient Stochastic Simulation Algorithm for Chemically Reacting Systems Based on Support Vector Regression

The stochastic simulation algorithm （SSA） accurately depicts spatially homogeneous wellstirred chemically reacting systems with small populations of chemical species and properly represents noise, but it is often abandoned when modeling larger systems because of its computational complexity. In this work, a twin support vector regression based stochastic simulations algorithm （TS^3A） is proposed by combining the twin support vector regression and SSA, the former is a well-known robust regression method in machine learning. Numerical results indicate that this proposed algorithm can be applied to a wide range of chemically reacting systems and obtain significant improvements on efficiency and accuracy with fewer simulating runs over the existing methods.     

An accelerated algorithm for discrete stochastic simulation of reaction-diffusion systems using gradient-based diffusion and tau-leaping

Stochastic simulation of reaction-diffusion systems enables the investigation of stochastic events arising from the small numbers and heterogeneous distribution of molecular species in biological cells. Stochastic variations in intracellular microdomains and in diffusional gradients play a significant part in the spatiotemporal activity and behavior of cells. Although an exact stochastic simulation that simulates every individual reaction and diffusion event gives a most accurate trajectory of the system's state over time, it can be too slow for many practical applications. We present an accelerated algorithm for discrete stochastic simulation of reaction-diffusion systems designed to improve the speed of simulation by reducing the number of time-steps required to complete a simulation run. This method is unique in that it employs two strategies that have not been incorporated in existing spatial stochastic simulation algorithms. First, diffusive transfers between neighboring subvolumes are based on concentration gradients. This treatment necessitates sampling of only the net or observed diffusion events from higher to lower concentration gradients rather than sampling all diffusion events regardless of local concentration gradients. Second, we extend the non-negative Poisson tau-leaping method that was originally developed for speeding up nonspatial or homogeneous stochastic simulation algorithms. This method calculates each leap time in a unified step for both reaction and diffusion processes while satisfying the leap condition that the propensities do not change appreciably during the leap and ensuring that leaping does not cause molecular populations to become negative. Numerical results are presented that illustrate the improvement in simulation speed achieved by incorporating these two new strategies.     

Accelerated stochastic simulation of the stiff enzyme-substrate reaction

The enzyme-catalyzed conversion of a substrate into a product is a common reaction motif in cellular chemical systems. In the three reactions that comprise this process, the intermediate enzyme-substrate complex is usually much more likely to decay into its original constituents than to produce a product molecule. This condition makes the reaction set mathematically "stiff." We show here how the simulation of this stiff reaction set can be dramatically speeded up relative to the standard stochastic simulation algorithm (SSA) by using a recently introduced procedure called the slow-scale SSA. The speedup occurs because the slow-scale SSA explicitly simulates only the relatively rare conversion reactions, skipping over occurrences of the other two less interesting but much more frequent reactions. We describe, explain, and illustrate this simulation procedure for the isolated enzyme-substrate reaction set, and then we show how the procedure extends to the more typical case in which the enzyme-substrate reactions occur together with other reactions and species. Finally, we explain the connection between this slow-scale SSA approach and the Michaelis-Menten [Biochem. Z. 49, 333 (1913)] formula, which has long been used in deterministic chemical kinetics to describe the enzyme-substrate reaction.     

Comparison of approaches for parameter estimation on stochastic models: Generic least squares versus specialized approaches

Parameter estimation for models with intrinsic stochasticity poses specific challenges that do not exist for deterministic models. Therefore, specialized numerical methods for parameter estimation in stochastic models have been developed. Here, we study whether dedicated algorithms for stochastic models are indeed superior to the naive approach of applying the readily available least squares algorithm designed for deterministic models.We compare the performance of the recently developed multiple shooting for stochastic systems (MSS) method designed for parameter estimation in stochastic models, a stochastic differential equations based Bayesian approach and a chemical master equation based techniques with the least squares approach for parameter estimation in models of ordinary differential equations (ODE). As test data, 1000 realizations of the stochastic models are simulated. For each realization an estimation is performed with each method, resulting in 1000 estimates for each approach. These are compared with respect to their deviation to the true parameter and, for the genetic toggle switch, also their ability to reproduce the symmetry of the switching behavior. Results are shown for different set of parameter values of a genetic toggle switch leading to symmetric and asymmetric switching behavior as well as an immigration-death and a susceptible-infected-recovered model. This comparison shows that it is important to choose a parameter estimation technique that can treat intrinsic stochasticity and that the specific choice of this algorithm shows only minor performance differences.     

Reaction Brownian dynamics and the effect of spatial fluctuations on the gain of a push-pull network

Brownian Dynamics algorithms have been widely used for simulating systems in soft-condensed matter physics. In recent times, their application has been extended to the simulation of coarse-grained models of biochemical networks. In these models, components move by diffusion and interact with one another upon contact. However, when reactions are incorporated into a Brownian dynamics algorithm, care must be taken to avoid violations of the detailed-balance rule, which would introduce systematic errors in the simulation. We present a Brownian dynamics algorithm for simulating reaction-diffusion systems that rigorously obeys detailed balance for equilibrium reactions. By comparing the simulation results to exact analytical results for a bimolecular reaction, we show that the algorithm correctly reproduces both equilibrium and dynamical quantities. We apply our scheme to a "push-pull" network in which two antagonistic enzymes covalently modify a substrate. Our results highlight that spatial fluctuations of the network components can strongly reduce the gain of the response of a biochemical network.     

Michaelis-Menten speeds up tau-leaping under a wide range of conditions

This paper examines the benefits of Michaelis-Menten model reduction techniques in stochastic tau-leaping simulations. Results show that although the conditions for the validity of the reductions for tau-leaping remain the same as those for the stochastic simulation algorithm (SSA), the reductions result in a substantial speed-up for tau-leaping under a different range of conditions than they do for SSA. The reason of this discrepancy is that the time steps for SSA and for tau-leaping are determined by different properties of system dynamics.     

Overcoming stiffness in stochastic simulation stemming from partial equilibrium: a multiscale Monte Carlo algorithm

In this paper the problem of stiffness in stochastic simulation of singularly perturbed systems is discussed. Such stiffness arises often from partial equilibrium or quasi-steady-state type of conditions. A multiscale Monte Carlo method is discussed that first assesses whether partial equilibrium is established using a simple criterion. The exact stochastic simulation algorithm (SSA) is next employed to sample among fast reactions over short time intervals (microscopic time steps) in order to compute numerically the proper probability distribution function for sampling the slow reactions. Subsequently, the SSA is used to sample among slow reactions and advance the time by large (macroscopic) time steps. Numerical examples indicate that not only long times can be simulated but also fluctuations are properly captured and substantial computational savings result.     

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.     

An adaptive stepsize method for the chemical Langevin equation

Mathematical and computational modeling are key tools in analyzing important biological processes in cells and living organisms. In particular, stochastic models are essential to accurately describe the cellular dynamics, when the assumption of the thermodynamic limit can no longer be applied. However, stochastic models are computationally much more challenging than the traditional deterministic models. Moreover, many biochemical systems arising in applications have multiple time-scales, which lead to mathematical stiffness. In this paper we investigate the numerical solution of a stochastic continuous model of well-stirred biochemical systems, the chemical Langevin equation. The chemical Langevin equation is a stochastic differential equation with multiplicative, non-commutative noise. We propose an adaptive stepsize algorithm for approximating the solution of models of biochemical systems in the Langevin regime, with small noise, based on estimates of the local error. The underlying numerical method is the Milstein scheme. The proposed adaptive method is tested on several examples arising in applications and it is shown to have improved efficiency and accuracy compared to the existing fixed stepsize schemes.