Adaptive stochastic-deterministic chemical kinetic simulations

MOTIVATION: Biochemical signaling pathways and genetic circuits often involve very small numbers of key signaling molecules. Computationally expensive stochastic methods are necessary to simulate such chemical situations. Single-molecule chemical events often co-exist with much larger numbers of signaling molecules where mass-action kinetics is a reasonable approximation. Here, we describe an adaptive stochastic method that dynamically chooses between deterministic and stochastic calculations depending on molecular count and propensity of forward reactions. The method is fixed timestep and has first order accuracy. We compare the efficiency of this method with exact stochastic methods. RESULTS: We have implemented an adaptive stochastic-deterministic approximate simulation method for chemical kinetics. With an error margin of 5%, the method solves typical biologically constrained reaction schemes more rapidly than exact stochastic methods for reaction volumes >1-10 micro m(3). We have developed a test suite of reaction cases to test the accuracy of mixed simulation methods. AVAILABILITY: Simulation software used in the paper is freely available from http://www.ncbs.res.in/kinetikit/download.html     

Hybrid simulation of cellular behavior

MOTIVATION: To be valuable to biological or biomedical research, in silico methods must be scaled to complex pathways and large numbers of interacting molecular species. The correct method for performing such simulations, discrete event simulation by Monte Carlo generation, is computationally costly for large complex systems. Approximation of molecular behavior by continuous models fails to capture stochastic behavior that is essential to many biological phenomena. RESULTS: We present a novel approach to building hybrid simulations in which some processes are simulated discretely, while other processes are handled in a continuous simulation by differential equations. This approach preserves the stochastic behavior of cellular pathways, yet enables scaling to large populations of molecules. We present an algorithm for synchronizing data in a hybrid simulation and discuss the trade-offs in such simulation. We have implemented the hybrid simulation algorithm and have validated it by simulating the statistical behavior of the well-known lambda phage switch. Hybrid simulation provides a new method for exploring the sources and nature of stochastic behavior in cells.     

Parameter estimation in stochastic biochemical reactions

Gene regulatory, signal transduction and metabolic networks are major areas of interest in the newly emerging field of systems biology. In living cells, stochastic dynamics play an important role; however, the kinetic parameters of biochemical reactions necessary for modelling these processes are often not accessible directly through experiments. The problem of estimating stochastic reaction constants from molecule count data measured, with error, at discrete time points is considered. For modelling the system, a hidden Markov process is used, where the hidden states are the true molecule counts, and the transitions between those states correspond to reaction events following collisions of molecules. Two different algorithms are proposed for estimating the unknown model parameters. The first is an approximate maximum likelihood method that gives good estimates of the reaction parameters in systems with few possible reactions in each sampling interval. The second algorithm, treating the data as exact measurements, approximates the number of reactions in each sampling interval by solving a simple linear equation. Maximising the likelihood based on these approximations can provide good results, even in complex reaction systems.     

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 modelling of biochemical systems of multi-step reactions using a simplified two-variable model

Background

A fundamental issue in systems biology is how to design simplified mathematical models for describing the dynamics of complex biochemical reaction systems. Among them, a key question is how to use simplified reactions to describe the chemical events of multi-step reactions that are ubiquitous in biochemistry and biophysics. To address this issue, a widely used approach in literature is to use one-step reaction to represent the multi-step chemical events. In recent years, a number of modelling methods have been designed to improve the accuracy of the one-step reaction method, including the use of reactions with time delay. However, our recent research results suggested that there are still deviations between the dynamics of delayed reactions and that of the multi-step reactions. Therefore, more sophisticated modelling methods are needed to accurately describe the complex biological systems in an efficient way.

Results

This work designs a two-variable model to simplify chemical events of multi-step reactions. In addition to the total molecule number of a species, we first introduce a new concept regarding the location of molecules in the multi-step reactions, which is the second variable to represent the system dynamics. Then we propose a simulation algorithm to compute the probability for the firing of the last step reaction in the multi-step events. This probability function is evaluated using a deterministic model of ordinary differential equations and a stochastic model in the framework of the stochastic simulation algorithm. The efficiency of the proposed two-variable model is demonstrated by the realization of mRNA degradation process based on the experimentally measured data.

Conclusions

Numerical results suggest that the proposed new two-variable model produces predictions that match the multi-step chemical reactions very well. The successful realization of the mRNA degradation dynamics indicates that the proposed method is a promising approach to reduce the complexity of biological systems.

     

Cellular growth and division in the Gillespie algorithm

Recent experimental studies elucidating the importance of noise in gene regulation have ignited widespread interest in Gillespie's stochastic simulation technique for biochemical networks. We formulate modifications to the Gillespie algorithm which are necessary to correctly simulate chemical reactions with time-dependent reaction rates. We concentrate on time dependence of kinetic rates arising from the periodic process of growth and division of the cellular volume, and demonstrate that a careful re-derivation of the Gillespie algorithm is important when all stochastically simulated reactions have rates slower or comparable to the cellular growth rate. For an unregulated single-gene system, we illustrate our findings using recently proposed hybrid simulation techniques, and systematically compare our algorithm with analytic results obtained from the chemical master equation.     

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.     

Comparing Stochastic Differential Equations and Agent-Based Modelling and Simulation for Early-Stage Cancer

There is great potential to be explored regarding the use of agent-based modelling and simulation as an alternative paradigm to investigate early-stage cancer interactions with the immune system. It does not suffer from some limitations of ordinary differential equation models, such as the lack of stochasticity, representation of individual behaviours rather than aggregates and individual memory. In this paper we investigate the potential contribution of agent-based modelling and simulation when contrasted with stochastic versions of ODE models using early-stage cancer examples. We seek answers to the following questions: (1) Does this new stochastic formulation produce similar results to the agent-based version? (2) Can these methods be used interchangeably? (3) Do agent-based models outcomes reveal any benefit when compared to the Gillespie results? To answer these research questions we investigate three well-established mathematical models describing interactions between tumour cells and immune elements. These case studies were re-conceptualised under an agent-based perspective and also converted to the Gillespie algorithm formulation. Our interest in this work, therefore, is to establish a methodological discussion regarding the usability of different simulation approaches, rather than provide further biological insights into the investigated case studies. Our results show that it is possible to obtain equivalent models that implement the same mechanisms; however, the incapacity of the Gillespie algorithm to retain individual memory of past events affects the similarity of some results. Furthermore, the emergent behaviour of ABMS produces extra patters of behaviour in the system, which was not obtained by the Gillespie algorithm.     

Exact Hybrid Particle/Population Simulation of Rule-Based Models of Biochemical Systems

Detailed modeling and simulation of biochemical systems is complicated by the problem of combinatorial complexity, an explosion in the number of species and reactions due to myriad protein-protein interactions and post-translational modifications. Rule-based modeling overcomes this problem by representing molecules as structured objects and encoding their interactions as pattern-based rules. This greatly simplifies the process of model specification, avoiding the tedious and error prone task of manually enumerating all species and reactions that can potentially exist in a system. From a simulation perspective, rule-based models can be expanded algorithmically into fully-enumerated reaction networks and simulated using a variety of network-based simulation methods, such as ordinary differential equations or Gillespie''s algorithm, provided that the network is not exceedingly large. Alternatively, rule-based models can be simulated directly using particle-based kinetic Monte Carlo methods. This “network-free” approach produces exact stochastic trajectories with a computational cost that is independent of network size. However, memory and run time costs increase with the number of particles, limiting the size of system that can be feasibly simulated. Here, we present a hybrid particle/population simulation method that combines the best attributes of both the network-based and network-free approaches. The method takes as input a rule-based model and a user-specified subset of species to treat as population variables rather than as particles. The model is then transformed by a process of “partial network expansion” into a dynamically equivalent form that can be simulated using a population-adapted network-free simulator. The transformation method has been implemented within the open-source rule-based modeling platform BioNetGen, and resulting hybrid models can be simulated using the particle-based simulator NFsim. Performance tests show that significant memory savings can be achieved using the new approach and a monetary cost analysis provides a practical measure of its utility.     

Smoldyn on graphics processing units: massively parallel Brownian dynamics simulations

Space is a very important aspect in the simulation of biochemical systems; recently, the need for simulation algorithms able to cope with space is becoming more and more compelling. Complex and detailed models of biochemical systems need to deal with the movement of single molecules and particles, taking into consideration localized fluctuations, transportation phenomena, and diffusion. A common drawback of spatial models lies in their complexity: models can become very large, and their simulation could be time consuming, especially if we want to capture the systems behavior in a reliable way using stochastic methods in conjunction with a high spatial resolution. In order to deliver the promise done by systems biology to be able to understand a system as whole, we need to scale up the size of models we are able to simulate, moving from sequential to parallel simulation algorithms. In this paper, we analyze Smoldyn, a widely diffused algorithm for stochastic simulation of chemical reactions with spatial resolution and single molecule detail, and we propose an alternative, innovative implementation that exploits the parallelism of Graphics Processing Units (GPUs). The implementation executes the most computational demanding steps (computation of diffusion, unimolecular, and bimolecular reaction, as well as the most common cases of molecule-surface interaction) on the GPU, computing them in parallel on each molecule of the system. The implementation offers good speed-ups and real time, high quality graphics output     

A multi-scaled approach for simulating chemical reaction systems

In this paper we give an overview of some very recent work, as well as presenting a new approach, on the stochastic simulation of multi-scaled systems involving chemical reactions. In many biological systems (such as genetic regulation and cellular dynamics) there is a mix between small numbers of key regulatory proteins, and medium and large numbers of molecules. In addition, it is important to be able to follow the trajectories of individual molecules by taking proper account of the randomness inherent in such a system. We describe different types of simulation techniques (including the stochastic simulation algorithm, Poisson Runge–Kutta methods and the balanced Euler method) for treating simulations in the three different reaction regimes: slow, medium and fast. We then review some recent techniques on the treatment of coupled slow and fast reactions for stochastic chemical kinetics and present a new approach which couples the three regimes mentioned above. We then apply this approach to a biologically inspired problem involving the expression and activity of LacZ and LacY proteins in E. coli, and conclude with a discussion on the significance of this work.     

Markov chain models of coupled calcium channels: Kronecker representations and iterative solution methods

Mathematical models of calcium release sites derived from Markov chain models of intracellular calcium channels exhibit collective gating reminiscent of the experimentally observed phenomenon of stochastic calcium excitability (i.e., calcium puffs and sparks). Calcium release site models are stochastic automata networks that involve many functional transitions, that is, the transition probabilities of each channel depend on the local calcium concentration and thus the state of the other channels. We present a Kronecker-structured representation for calcium release site models and perform benchmark stationary distribution calculations using both exact and approximate iterative numerical solution techniques that leverage this structure. When it is possible to obtain an exact solution, response measures such as the number of channels in a particular state converge more quickly using the iterative numerical methods than occupation measures calculated via Monte Carlo simulation. In particular, multi-level methods provide excellent convergence with modest additional memory requirements for the Kronecker representation of calcium release site models. When an exact solution is not feasible, iterative approximate methods based on the power method may be used, with performance similar to Monte Carlo estimates. This suggests approximate methods with multi-level iterative engines as a promising avenue of future research for large-scale calcium release site models.     

Complex discrete dynamics from simple continuous population models

Nonoverlapping generations have been classically modelled as difference equations in order to account for the discrete nature of reproductive events. However, other events such as resource consumption or mortality are continuous and take place in the within-generation time. We have realistically assumed a hybrid ODE bidimensional model of resources and consumers with discrete events for reproduction. Numerical and analytical approaches showed that the resulting dynamics resembles a Ricker map, including the doubling route to chaos. Stochastic simulations with a handling-time parameter for indirect competition of juveniles may affect the qualitative behaviour of the model.     

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.     

Mathematical models of eye movements in reading: a possible role for autonomous saccades

 An efficient method for the exact numerical simulation of semi-Markov processes is used to study minimal models of the control of eye movements in reading. When we read a text, typical sequences of fixations form a rather complicated trajectory – almost like a random walk. Mathematical models of eye movement control can account for this behavior using stochastic transition rules between few discrete internal states, which represent combinations of certain stages of lexical access and saccade programs. We show that experimentally observed fixation durations can be explained by residence-time-dependent transition probabilities. Stochastic processes with this property are known as semi-Markov processes. For our numerical simulations we use the minimal process method (Gillespie algorithm), which is an exact and efficient simulation algorithm for this class of stochastic processes. Within this mathematical framework, we study different forms of coupling between eye movements and shifts of covert attention in reading. Our model lends support to the existence of autonomous saccades, i.e., the hypothesis that initiations of saccades are not completely determined by lexical access processes. Received: 21 March 2000 / Accepted in revised form: 10 January 2001