Hybrid modeling and simulation of stochastic effects on progression through the eukaryotic cell cycle

The eukaryotic cell cycle is regulated by a complicated chemical reaction network. Although many deterministic models have been proposed, stochastic models are desired to capture noise in the cell resulting from low numbers of critical species. However, converting a deterministic model into one that accurately captures stochastic effects can result in a complex model that is hard to build and expensive to simulate. In this paper, we first apply a hybrid (mixed deterministic and stochastic) simulation method to such a stochastic model. With proper partitioning of reactions between deterministic and stochastic simulation methods, the hybrid method generates the same primary characteristics and the same level of noise as Gillespie's stochastic simulation algorithm, but with better efficiency. By studying the results generated by various partitionings of reactions, we developed a new strategy for hybrid stochastic modeling of the cell cycle. The new approach is not limited to using mass-action rate laws. Numerical experiments demonstrate that our approach is consistent with characteristics of noisy cell cycle progression, and yields cell cycle statistics in accord with experimental observations.     

Automatic identification of model reductions for discrete stochastic simulation

Multiple time scales in cellular chemical reaction systems present a challenge for the efficiency of stochastic simulation. Numerous model reductions have been proposed to accelerate the simulation of chemically reacting systems by exploiting time scale separation. However, these are often identified and deployed manually, requiring expert knowledge. This is time-consuming, prone to error, and opportunities for model reduction may be missed, particularly for large models. We propose an automatic model analysis algorithm using an adaptively weighted Petri net to dynamically identify opportunities for model reductions for both the stochastic simulation algorithm and tau-leaping simulation, with no requirement of expert knowledge input. Results are presented to demonstrate the utility and effectiveness of this approach.     

Binomial tau-leap spatial stochastic simulation algorithm for applications in chemical kinetics

In cell biology, cell signaling pathway problems are often tackled with deterministic temporal models, well mixed stochastic simulators, and/or hybrid methods. But, in fact, three dimensional stochastic spatial modeling of reactions happening inside the cell is needed in order to fully understand these cell signaling pathways. This is because noise effects, low molecular concentrations, and spatial heterogeneity can all affect the cellular dynamics. However, there are ways in which important effects can be accounted without going to the extent of using highly resolved spatial simulators (such as single-particle software), hence reducing the overall computation time significantly. We present a new coarse grained modified version of the next subvolume method that allows the user to consider both diffusion and reaction events in relatively long simulation time spans as compared with the original method and other commonly used fully stochastic computational methods. Benchmarking of the simulation algorithm was performed through comparison with the next subvolume method and well mixed models (MATLAB), as well as stochastic particle reaction and transport simulations (CHEMCELL, Sandia National Laboratories). Additionally, we construct a model based on a set of chemical reactions in the epidermal growth factor receptor pathway. For this particular application and a bistable chemical system example, we analyze and outline the advantages of our presented binomial tau-leap spatial stochastic simulation algorithm, in terms of efficiency and accuracy, in scenarios of both molecular homogeneity and heterogeneity.     

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.     

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

The sorting direct method for stochastic simulation of biochemical systems with varying reaction execution behavior

A key to advancing the understanding of molecular biology in the post-genomic age is the development of accurate predictive models for genetic regulation, protein interaction, metabolism, and other biochemical processes. To facilitate model development, simulation algorithms must provide an accurate representation of the system, while performing the simulation in a reasonable amount of time. Gillespie's stochastic simulation algorithm (SSA) accurately depicts spatially homogeneous models 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, we examine the performance of different versions of the SSA when applied to several biochemical models. Through our analysis, we discover that transient changes in reaction execution frequencies, which are typical of biochemical models with gene induction and repression, can dramatically affect simulator performance. To account for these shifts, we propose a new algorithm called the sorting direct method that maintains a loosely sorted order of the reactions as the simulation executes. Our measurements show that the sorting direct method performs favorably when compared to other well-known exact stochastic simulation algorithms.     

Modeling signal transduction networks: a comparison of two stochastic kinetic simulation algorithms

Computational efficiency of stochastic kinetic algorithms depend on factors such as the overall species population, the total number of reactions, and the average number of nodal interactions or connectivity in a network. These size measures of the network model can have a significant impact on computational efficiency. In this study, two scalable biological networks are used to compare the size scaling efficiencies of two popular and conceptually distinct stochastic kinetic simulation algorithms--the random substrate method of Firth and Bray (FB), and the Gillespie algorithm as implemented using the Gibson-Bruck method (GGB). The arithmetic computational efficiencies of these two algorithms, respectively, scale with the square of the total species population and the logarithm of the total number of active reactions. The two scalable models considered are the size scalable model (SSM), a four compartment reaction model for a signal transduction network involving receptors with single phosphorylation binding sites, and the variable connectivity model (VCM), a single compartment model where receptors possess multiple phosphorylation binding sites. The SSM has fixed species connectivity while the connectivity between species in VCM increases with the number of phosphorylation sites. For SSM, we find that, as the total species population is increased over four orders of magnitude, the GGB algorithm performs significantly better than FB for all three SSM compartment models considered. In contrast, for VCM, we find that as the overall species population decreases while the number of phosphorylation sites increases (implying an increase in network linkage) there exists a crossover point where the computational demands of the GGB method exceed that of the FB.     

Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions

The dynamical solution of a well-mixed, nonlinear stochastic chemical kinetic system, described by the Master equation, may be exactly computed using the stochastic simulation algorithm. However, because the computational cost scales with the number of reaction occurrences, systems with one or more "fast" reactions become costly to simulate. This paper describes a hybrid stochastic method that partitions the system into subsets of fast and slow reactions, approximates the fast reactions as a continuous Markov process, using a chemical Langevin equation, and accurately describes the slow dynamics using the integral form of the "Next Reaction" variant of the stochastic simulation algorithm. The key innovation of this method is its mechanism of efficiently monitoring the occurrences of slow, discrete events while simultaneously simulating the dynamics of a continuous, stochastic or deterministic process. In addition, by introducing an approximation in which multiple slow reactions may occur within a time step of the numerical integration of the chemical Langevin equation, the hybrid stochastic method performs much faster with only a marginal decrease in accuracy. Multiple examples, including a biological pulse generator and a large-scale system benchmark, are simulated using the exact and proposed hybrid methods as well as, for comparison, a previous hybrid stochastic method. Probability distributions of the solutions are compared and the weak errors of the first two moments are computed. In general, these hybrid methods may be applied to the simulation of the dynamics of a system described by stochastic differential, ordinary differential, and Master equations.     

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.     

Exact on-lattice stochastic reaction-diffusion simulations using partial-propensity methods

Stochastic reaction-diffusion systems frequently exhibit behavior that is not predicted by deterministic simulation models. Stochastic simulation methods, however, are computationally expensive. We present a more efficient stochastic reaction-diffusion simulation algorithm that samples realizations from the exact solution of the reaction-diffusion master equation. The present algorithm, called partial-propensity stochastic reaction-diffusion (PSRD) method, uses an on-lattice discretization of the reaction-diffusion system and relies on partial-propensity methods for computational efficiency. We describe the algorithm in detail, provide a theoretical analysis of its computational cost, and demonstrate its computational performance in benchmarks. We then illustrate the application of PSRD to two- and three-dimensional pattern-forming Gray-Scott systems, highlighting the role of intrinsic noise in these systems.     

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.     

Moment estimation for chemically reacting systems by extended Kalman filtering

In stochastic models of chemically reacting systems that contain bimolecular reactions, the dynamics of the moments of order up to n of the species populations do not form a closed system, in the sense that their time-derivatives depend on moments of order n + 1. To close the dynamics, the moments of order n + 1 are generally approximated by nonlinear functions of the lower order moments. If the molecule counts of some of the species have a high probability of becoming zero, such approximations may lead to imprecise results and stochastic simulation is the only viable alternative for system analysis. Stochastic simulation can produce exact realizations of chemically reacting systems, but tends to become computationally expensive, especially for stiff systems that involve reactions at different time scales. Further, in some systems, important stochastic events can be very rare and many simulations are necessary to obtain accurate estimates. The computational cost of stochastic simulation can then be prohibitively large. In this paper, we propose a novel method for estimating the moments of chemically reacting systems. The method is based on closing the moment dynamics by replacing the moments of order n + 1 by estimates calculated from a small number of stochastic simulation runs. The resulting stochastic system is then used in an extended Kalman filter, where estimates of the moments of order up to n, obtained from the same simulation, serve as outputs of the system. While the initial motivation for the method was improving over the performance of stochastic simulation and moment closure methods, we also demonstrate that it can be used in an experimental setting to estimate moments of species that cannot be measured directly from time course measurements of the moments of other species.     

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.     

Look before you leap: a confidence-based method for selecting species criticality while avoiding negative populations in τ-leaping

The stochastic simulation algorithm was introduced by Gillespie and in a different form by Kurtz. There have been many attempts at accelerating the algorithm without deviating from the behavior of the simulated system. The crux of the explicit τ-leaping procedure is the use of Poisson random variables to approximate the number of occurrences of each type of reaction event during a carefully selected time period, τ. This method is acceptable providing the leap condition, that no propensity function changes "significantly" during any time-step, is met. Using this method there is a possibility that species numbers can, artificially, become negative. Several recent papers have demonstrated methods that avoid this situation. One such method classifies, as critical, those reactions in danger of sending species populations negative. At most, one of these critical reactions is allowed to occur in the next time-step. We argue that the criticality of a reactant species and its dependent reaction channels should be related to the probability of the species number becoming negative. This way only reactions that, if fired, produce a high probability of driving a reactant population negative are labeled critical. The number of firings of more reaction channels can be approximated using Poisson random variables thus speeding up the simulation while maintaining the accuracy. In implementing this revised method of criticality selection we make use of the probability distribution from which the random variable describing the change in species number is drawn. We give several numerical examples to demonstrate the effectiveness of our new method.     

Dual stacked partial least squares for analysis of near-infrared spectra

A new ensemble learning algorithm is presented for quantitative analysis of near-infrared spectra. The algorithm contains two steps of stacked regression and Partial Least Squares (PLS), termed Dual Stacked Partial Least Squares (DSPLS) algorithm. First, several sub-models were generated from the whole calibration set. The inner-stack step was implemented on sub-intervals of the spectrum. Then the outer-stack step was used to combine these sub-models. Several combination rules of the outer-stack step were analyzed for the proposed DSPLS algorithm. In addition, a novel selective weighting rule was also involved to select a subset of all available sub-models. Experiments on two public near-infrared datasets demonstrate that the proposed DSPLS with selective weighting rule provided superior prediction performance and outperformed the conventional PLS algorithm. Compared with the single model, the new ensemble model can provide more robust prediction result and can be considered an alternative choice for quantitative analytical applications.     

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.     

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.     

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.