Skip Navigation


Bioinformatics Advance Access originally published online on February 4, 2005
Bioinformatics 2005 21(9):2136-2137; doi:10.1093/bioinformatics/bti308
This Article
Right arrow Abstract Freely available
Right arrow FREE Full Text (Print PDF) Freely available
Right arrow All Versions of this Article:
21/9/2136    most recent
bti308v1
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Add to My Personal Archive
Right arrow Download to citation manager
Right arrow Search for citing articles in:
ISI Web of Science (11)
Right arrowRequest Permissions
Google Scholar
Right arrow Articles by Chatterjee, A.
Right arrow Articles by Vlachos, D. G.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Chatterjee, A.
Right arrow Articles by Vlachos, D. G.
Social Bookmarking
Add to CiteULike   Add to Connotea   Add to Del.icio.us  
What's this?

© The Author 2005. Published by Oxford University Press. All rights reserved. For Permissions, please email: journals.permissions{at}oupjournals.org

Time accelerated Monte Carlo simulations of biological networks using the binomial {tau}-leap method

Abhijit Chatterjee , Kapil Mayawala , Jeremy S. Edwards and Dionisios G. Vlachos *

Department of Chemical Engineering, University of Delaware Newark, DE 19716, USA

*To whom correspondence should be addressed.


   Abstract
TOP
 Abstract
INTRODUCTION
 THE BINOMIAL {tau}-LEAP METHOD
 RESULTS
 REFERENCES
 

Summary: Developing a quantitative understanding of intracellular networks requires simulations and computational analyses. However, traditional differential equation modeling tools are often inadequate due to the stochasticity of intracellular reaction networks that can potentially influence the phenotypic characteristics. Unfortunately, stochastic simulations are computationally too intense for most biological systems. Herein, we have utilized the recently developed binomial {tau}-leap method to carry out stochastic simulations of the epidermal growth factor receptor induced mitogen activated protein kinase cascade. Results indicate that the binomial {tau}-leap method is computationally 100–1000 times more efficient than the exact stochastic simulation algorithm of Gillespie. Furthermore, the binomial {tau}-leap method avoids negative populations and accurately captures the species populations along with their fluctuations despite the large difference in their size.

Availability: http://www.dion.che.udel.edu/multiscale/Introduction.html. Fortran 90 code available for academic use by email.

Contact: vlachos{at}che.udel.edu

Supplementary information: Details about the binomial {tau}-leap algorithm, software and a manual are available at the above website.


   INTRODUCTION
TOP
 Abstract
 INTRODUCTION
THE BINOMIAL {tau}-LEAP METHOD
 RESULTS
 REFERENCES
 
The importance of stochasticity in biological systems is well established by several theoretical studies (Arkin et al., 1998; Morton-Firth and Bray, 1998; Resat et al., 2003; Meng et al., 2004) and experimental work (Elowitz et al., 2002; Blake et al., 2003). Stochasticity in chemical systems has been extensively studied using the exact stochastic simulation algorithm (SSA) of Gillespie (1976). However, SSA is limited to non-stiff systems on short times. Biological networks exhibit multiple time scales, ranging from microseconds to days (Edwards and Palsson, 2000), and demand efficient stochastic simulation algorithms (Lok, 2004). Such recently developed algorithms are reviewed in Turner et al. (2004) and Vlachos (2005). The family of {tau} -leap methods is the only coarse-grained stochastic simulation method that can provide accurate simulation results.

Here the binomial {tau}-leap method (Chatterjee et al., 2005), a variant of the original Poisson {tau}-leap method of Gillespie (2001) is employed to study the signaling pathway of EGF receptor (EGFR) activated mitogen activated protein (MAP) kinase cascade. Activated EGF receptors trigger a rich network of signaling pathways and regulate cell functions such as proliferation, differentiation and migration (Yarden and Gur, 2004). In this work, we have used a mathematical model of this signaling pathway to study the activation of the MAP kinase cascade from surface and internalized EGF receptors (Schoeberl et al., 2002). The published model is based on ordinary differential equations (ODEs) of 94 signaling species with 296 reactions (referred to as events in this paper) in a well-mixed environment. The complexity of the model renders the MAP kinase cascade an excellent benchmark problem for testing the applicability of the binomial {tau}-leap method for biological systems.


   THE BINOMIAL {tau}-LEAP METHOD
TOP
 Abstract
 INTRODUCTION
 THE BINOMIAL {tau}-LEAP METHOD
RESULTS
 REFERENCES
 
The binomial {tau}-leap method follows the stochastic evolution of N species among M reactions as an approximate Markov process. The algorithm requires the initial population size of all species (i.e. the number of molecules) and the reaction kinetics as inputs. This information is used to compute transition probabilities per unit time of all events at time t. A ‘bundle’ of events, sampled from a binomial distribution, are allowed to trigger in a time interval of size {tau}, and the time is incremented to t+{tau}. This process is repeated as the algorithm leaps along the time axis. The maximum {tau} is chosen to prevent substantial changes in populations (Gillespie, 2001). This is automatically taken care of using a simple, computationally inexpensive criterion (Chatterjee et al., 2005) according to which the user can control the speed-up by inputting the coarse-graining factor r of time stepping in comparison to the SSA. This criterion is chosen here because r is constrained by mass conservation, 0 < r≤ 1, but future work is needed to develop adaptive time step methods. Previous experience indicates that r up to 0.2 provides accurate simulations that are much faster than the SSA. A similar binomial {tau}-leap method has been developed independently by Tian and Burrage (2004).

The algorithm conserves mass at both the single reaction and the entire network levels by constraining the number of reaction firings, thereby eliminating for the first time negative populations in {tau}-leaping (Chatterjee et al., 2005). Numerical simulations using simple reaction networks (Chatterjee et al., 2005) demonstrated that the method is more accurate than the original Poisson {tau}-leap method of Gillespie (2001).


   RESULTS
TOP
 Abstract
 INTRODUCTION
 THE BINOMIAL {tau}-LEAP METHOD
 RESULTS
REFERENCES
 
Simulations spanning 15 min of real time were performed for the MAP kinase cascade using the binomial {tau}-leap method, with different levels of coarse-graining r=0.1–0.5, and the SSA. Reaction dynamics and initial populations are specified on the Website (http://www.dion.che.udel.edu/multiscale/Introduction.html). Intracellular populations, namely free dimerized EGF bound EGFR ((EGF-EGFR)2), doubly phosphorylated ERK (ERK-PP), doubly phosphorylated MEK (MEK-PP) and free and activated RaF (RaF*), are plotted against time in Figure 1 for selected species whose populations span from small to large ones. The binomial {tau}-leap method accurately captures the network kinetics for all times despite the large disparity in species population size in this complex network. The Poisson {tau}-leap method on the other hand encounters negative populations since small populations are present (Gillespie, 2001).



View larger version (29K):
[in this window]
[in a new window]
  		Fig. 1 Time-dependent populations, denoted as Xi (number of molecules) for species i, using the stochastic simulation algorithm (SSA; thick solid line) and the binomial {tau}-leap method (thin solid line, r = 0.2 and squares, r = 0.5).

 
Figure 2 compares fluctuations of the two methods using the probability density functions obtained at time t=1 min (with ensemble averages over 500 trajectories for the SSA and 1000 for the binomial {tau}-leap method). The mean and standard deviation for representative species, namely internalized dimerized EGF bound EGFR ((EGF-EGFRi)2), internalized EGFR (EGFRi) and activated (EGF-EGFR)2 (EGF-EGFR*2), that span many orders of magnitude in number, are tabulated in the inset of Figure 2. The consistency between the methods for the mean, standard deviation and higher order moments is indicative of the accuracy that can be obtained using the binomial {tau}-leap method.



View larger version (37K):
[in this window]
[in a new window]
  		Fig. 2 Probability density functions (pdf) using the stochastic simulation algorithm (SSA, circles) and the binomial {tau}-leap method at time t=1 min (squares, r=0.2; diamonds, r=0.5). Populations for species i are denoted as Xi (number of molecules). The solid lines are curve fits for visual aid. Average populations and standard deviation are tabulated in the inset.

 
The SSA required 3 days on a DEC alpha 833 MHz processor for the simulation depicted in Figure 1 as compared to 6–30 min taken using the {tau}-leap method. The binomial {tau}-leap is 100–1000 times more efficient than the SSA over the range of r values explored. The low computational requirements of the binomial {tau}-leap method enable accurate long time simulations that are beyond the reach of SSA (see example on the Website). Due to its accuracy, robustness, simplicity and substantial speed-up, the binomial {tau}-leap method appears to be a particularly promising, approximate stochastic simulation method for computational biology.


   Acknowledgments
 
This work was partially supported by the NSF through CTS-0312117.

Received on January 9, 2005; revised on February 2, 2005; accepted on February 3, 2005

   REFERENCES
TOP
 Abstract
 INTRODUCTION
 THE BINOMIAL {tau}-LEAP METHOD
 RESULTS
 REFERENCES
 

    Arkin, A.P., et al. (1998) Stochastic kinetic analysis of developmental pathway bifurcation in phage l-infected Escherichia coli cells. Genetics, 149, 1633–1648[Abstract/Free Full Text].

    Blake, W.J., et al. (2003) Noise in eukaryotic gene expression. Nature, 422, 633–637[CrossRef][Medline].

    Chatterjee, A., et al. (2005) Binomial distribution based {tau}-leap accelerated stochastic simulation. J. Chem. Phys., 122, 024112[Medline].

    Edwards, J.S. and Palsson, B.O. (2000) Multiple steady states in kinetic models of red cell metabolism. J. Theoret. Biol., 207, 125–127[CrossRef][Medline].

    Elowitz, M.B., et al. (2002) Stochastic gene expression in a single cell. Science, 297, 1183–1186[Abstract/Free Full Text].

    Gillespie, D.T. (1976) A general method for numerically simulating the stochastic evolution of coupled chemical reactions. J. Comput. Phys., 22, 403–434[CrossRef].

    Gillespie, D.T. (2001) Approximate accelerated stochastic simulation of chemically reacting systems. J. Chem. Phys., 115, 1716–1733[CrossRef].

    Lok, L. (2004) The need for speed in stochastic simulation. Nat. Biotechnol., 22, 964–965[Medline].

    Meng, T.C., et al. (2004) Modeling and simulation of biological systems with stochasticity. In Silico Biol., 4, 0024.

    Morton-Firth, C.J. and Bray, D. (1998) Predicting temporal fluctuations in an intracellular signalling pathway. J. Theoret. Biol., 192, 117–128[CrossRef][ISI][Medline].

    Resat, H., et al. (2003) An integrated model of epidermal growth factor receptor trafficking and signal transduction. Biophys. J., 85, 730–743[Abstract/Free Full Text].

    Schoeberl, B., et al. (2002) Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized receptors. Nat. Biotechnol., 20, 370–375[CrossRef][ISI][Medline].

    Tian, T. and Burrage, K. (2004) Binomial leap methods for simulating stochastic chemical kinetics. J. Chem. Phys., 121, 10356–10364[CrossRef][Medline].

    Turner, T.E., et al. (2004) Stochastic approaches for modelling in vivo reactions. Comput. Biol. Chem., 28, 165–178[CrossRef][ISI][Medline].

    Vlachos, D.G. (2005) The emerging field of multiscale analysis: a review with examples from systems biology, materials engineering, and fluid-surface interacting systems. Adv. Chem. Eng, in press.

    Yarden, Y. and Gur, G. (2004) Enlightened receptor dynamics. Nat. Biotechnol., 22, 169–170[Medline].


Add to CiteULike CiteULike   Add to Connotea Connotea   Add to Del.icio.us Del.icio.us    What's this?


This article has been cited by other articles:


Home page
 Biophys. JHome page
A.-T. Dinh, C. Pangarkar, T. Theofanous, and S. Mitragotri
Understanding Intracellular Transport Processes Pertinent to Synthetic Gene Delivery via Stochastic Simulations and Sensitivity Analyses
Biophys. J., February 1, 2007; 92(3): 831 - 846.
[Abstract] [Full Text] [PDF]

Home page
 BioinformaticsHome page
T. Tian, S. Xu, J. Gao, and K. Burrage
Simulated maximum likelihood method for estimating kinetic rates in gene expression
Bioinformatics, January 1, 2007; 23(1): 84 - 91.
[Abstract] [Full Text] [PDF]

This Article
Right arrow Abstract Freely available
Right arrow FREE Full Text (Print PDF) Freely available
Right arrow All Versions of this Article:
21/9/2136    most recent
bti308v1
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Add to My Personal Archive
Right arrow Download to citation manager
Right arrow Search for citing articles in:
ISI Web of Science (11)
Right arrowRequest Permissions
Google Scholar
Right arrow Articles by Chatterjee, A.
Right arrow Articles by Vlachos, D. G.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Chatterjee, A.
Right arrow Articles by Vlachos, D. G.
Social Bookmarking
Add to CiteULike   Add to Connotea   Add to Del.icio.us  
What's this?