A deterministic map of Waddington's epigenetic landscape for cell fate specification
© Bhattacharya et al; licensee BioMed Central Ltd. 2011
Received: 17 February 2011
Accepted: 27 May 2011
Published: 27 May 2011
The image of the "epigenetic landscape", with a series of branching valleys and ridges depicting stable cellular states and the barriers between those states, has been a popular visual metaphor for cell lineage specification - especially in light of the recent discovery that terminally differentiated adult cells can be reprogrammed into pluripotent stem cells or into alternative cell lineages. However the question of whether the epigenetic landscape can be mapped out quantitatively to provide a predictive model of cellular differentiation remains largely unanswered.
Here we derive a simple deterministic path-integral quasi-potential, based on the kinetic parameters of a gene network regulating cell fate, and show that this quantity is minimized along a temporal trajectory in the state space of the gene network, thus providing a marker of directionality for cell differentiation processes. We then use the derived quasi-potential as a measure of "elevation" to quantitatively map the epigenetic landscape, on which trajectories flow "downhill" from any location. Stochastic simulations confirm that the elevation of this computed landscape correlates to the likelihood of occurrence of particular cell fates, with well-populated low-lying "valleys" representing stable cellular states and higher "ridges" acting as barriers to transitions between the stable states.
This quantitative map of the epigenetic landscape underlying cell fate choice provides mechanistic insights into the "forces" that direct cellular differentiation in the context of physiological development, as well as during artificially induced cell lineage reprogramming. Our generalized approach to mapping the landscape is applicable to non-gradient gene regulatory systems for which an analytical potential function cannot be derived, and also to high-dimensional gene networks. Rigorous quantification of the gene regulatory circuits that govern cell lineage choice and subsequent mapping of the epigenetic landscape can potentially help identify optimal routes of cell fate reprogramming.
In the quantitative view of a cell as a dynamical system governed by genetic interaction networks , an intuitive association can be made between the valleys ("creodes" in Waddington's terminology) on the epigenetic landscape and the trajectories leading to the attractors, or stable steady states, of the gene networks that regulate cell fate [7–9]. But can we quantitatively map the undulating surface of the landscape, thereby providing a predictive model of the "directionality" of cellular differentiation? Waddington himself cautioned that the epigenetic landscape, while useful as a "rough and ready picture" of development, "cannot be interpreted rigorously ". The mathematician René Thom, in his formulation of catastrophe theory inspired by Waddington's ideas, proposed that a generalized "potential surface" could be derived for any dynamical system [2, 10]. However Thom's later writings suggest that he did not believe it possible to quantify the epigenetic landscape . This view has been echoed by other authors, who have described the landscape as a "colorful metaphor " with "no grounding in physical reality ".
Huang, Wang and colleagues have recently proposed a probabilistic "pseudo-potential" to quantify the epigenetic landscape for a gene network regulating cell fate, where the elevation of the surface is inversely related to the likelihood of occurrence of a particular state in phase space [8, 12, 13]. In this formulation a stochastic potential energy landscape is characterized for a gene network, based on a Hartree mean-field approximation of the underlying master equation . Such stochastic formulations have also been used to derive probabilistic potential landscapes for the lysis-lysogeny switch in bacteriophage lambda [15–17], the mitogen-activated protein kinase (MAPK) signal transduction network , biochemical oscillations , and the predator-prey system .
Here we propose a simple numerical method to map the epigenetic landscape that is not based on a probabilistic or master-equation approach. Instead, a quasi-potential surface (Figure 1B) is derived directly from the deterministic rate equations governing the dynamic behavior of a gene regulatory circuit. We then use stochastic simulations to show that the elevation of this computed landscape correlates to the likelihood of occurrence of particular cell fates, with well-populated low-lying valleys representing stable cellular states and higher ridges acting as barriers to transitions between the stable states.
Finally, we discuss ways in which this quantitatively mapped landscape may help predict the efficiency of cellular de-differentiation or trans-differentiation, and identify optimal routes of cell fate reprogramming. Recent discoveries have challenged the dogma of cell fate determination as a unidirectional and irreversible process. Even terminally differentiated adult cells have now been shown to retain considerable phenotypic plasticity and the ability to be reprogrammed into pluripotent stem cell-like states [21–27] or into alternative differentiated lineages [28–34] by forced expression of a single gene or a small number of genes. These findings have led to a resurgence of interest in Waddington's ideas about cell lineage choice, with several authors invoking the image of the epigenetic landscape [7–9, 35–39]. However the theoretical basis of plasticity in cell fate is still not fully understood, and the efficiency of reprogramming in these studies is often quite low . A quantitative understanding of the "forces" that drive cell differentiation, and the "barriers" that separate stable cell states, is urgently needed. Such understanding may eventually enable us to predict the relative ease or difficulty of de-differentiation or trans-differentiation among multiple cellular states.
Results and Discussion
Derivation of the quasi-potential landscape
In general, condition (3) will not be valid for an arbitrary circuit of two genes x and y that regulate each other as per Eq. 1, making it impossible to derive a closed-form potential function.
where Δx and Δy are sufficiently small increments along the trajectory such that and can be assumed to remain unchanged over the interval [(x, x+Δx); (y, y+Δy)]. The quantities Δx and Δy are obtained as the products and , respectively, where Δt is the time increment. We use the term "quasi-potential" to describe V q , to emphasize its distinction from a closed-form potential function.
For positive increments in time Δt, Δ V q is thus always negative along an evolving trajectory, ensuring that trajectories flow "downhill" along a putative "quasi-potential surface". Stable steady states of the system (dx/dt = 0; dy/dt = 0) would correspond to local minima on this quasi-potential surface, given that at these states Δ V q = 0 (per Eq. 5). The overall change in the quasi-potential along a trajectory can then be calculated by numerically integrating the quantity Δ V q in Eq. 4 from a given initial configuration up to a stable steady state, thereby allowing us to map out a temporal trajectory along the putative quasi-potential surface (Figure 2B). The quasi-potential thus defined is a measurable quantity that is minimized along a trajectory from any initial condition to an attractor in the phase space of the two genes, and is in effect a Liapunov function of the dynamical system represented by the two-gene circuit .
The same procedure can be applied to systems with more than two stable steady states - for instance, a "tristable" system produced by a circuit of two genes that induce their own expression, in addition to mutual inhibition (Figure 3C, D). This system has three steady states - two of which represent alternative differentiated cell lineages, while the third state depicts the common progenitor cell of the two lineages [8, 13, 42].
Quantitative interpretation of the quasi-potential landscape
The "third dimension" (elevation) of the landscape represented by the quasi-potential, although directly derived from the dynamic rate equations without any additional information, thus yields an interpretation of cellular stability not immediately apparent from two-dimensional phase portrait analysis. The analysis above supports the contention that the length of the "least action trajectory" along the contours of the epigenetic landscape is more important in predicting transitions between alternative cellular states than the simple "aerial distance" in state space . It is also interesting to note that the contours of the quantitatively mapped epigenetic landscape act as a constraint on the extent of stochastic fluctuations in protein levels, with simulated cells "smeared out" on the surface of a shallower basin (Figure 4A, middle panel) compared to a tighter distribution of cells on a deeper valley (Figure 4D, middle panel).
These results suggest that calculating the relative heights of the ridges and valleys on the computed epigenetic landscape of a multi-gene system can help predict the probability of trans-differentiation from one cell lineage to another, or de-differentiation of a particular cell type to its progenitor state. Current efforts to reprogram cell fate with potential application in regenerative medicine suffer from a low rate of successful reprogramming  and a trial-and-error approach to choice of a reprogramming strategy . Computing the epigenetic landscape for the critical gene interactions regulating the transition between two cellular states may indicate particular genetic manipulations that would lower the barriers separating the two states, thereby increasing the efficiency of the reprogramming process. It can also help characterize the relative ease or difficulty of alternative routes of cell fate transition [7, 9, 36]. For instance, comparison of the elevation of the barriers separating two terminally-differentiated cell lineages on the epigenetic landscape might suggest that de-differentiation of cells of one lineage to the common progenitor cell of the two lineages followed by redirection to the second lineage would lead to more efficient reprogramming than direct trans-differentiation (Figure S2, Additional File 1).
A dynamic landscape
Interestingly, this flexibility of the quasi-potential surface under gene manipulation gives a quantitative interpretation of the revised image of the epigenetic landscape proposed by Waddington (Figure S3, Additional File 1), which showed an array of pegs representing genes, holding up a sheet of fabric (the landscape) through a network of guy ropes (gene interactions) - meant to convey the idea that "the modelling of the epigenetic landscape ... is controlled by the pull of these numerous guy-ropes which are ultimately anchored to the genes ". Similar changes in the shape of the epigenetic landscape may also be brought about by external signals - for example endogenous cytokines or environmental chemicals - which by transiently altering the landscape could have an instructive effect on cell fate choice.
In this work, we have defined a deterministic quasi-potential that is minimized along a temporal trajectory followed by a gene network, and used it to quantitatively derive the corresponding epigenetic landscape. A gene network not being a mechanical system, this quasi-potential should not be confused with a potential energy function. It is rather a Liapunov function of the dynamical system represented by the gene network, along which trajectories flow monotonically "downhill" towards the steady states of the network . Other investigators have used a term analogous to the quasi-potential difference Δ V q in Eq. 4 to calculate the "energy landscape" for concentrations of one component in a gene network [46, 47]. Here we have used the concept of alignment of multiple trajectories to interpolate the epigenetic landscape of a two-variable system.
This novel and simple process for deriving the surface of the landscape from a path-integral quasi-potential is not restricted to two-gene systems. While the landscape cannot be visually rendered for circuits with more than two genes, the rates of transition across the potential barriers between multiple steady states in the system can still be computed to predict optimum routes of cell fate reprogramming.
However, many binary branching points in development, particularly in blood cell lineage specification, are governed by mutual antagonism of only two transcription factors associated with alternative lineage choices . Mapping the epigenetic landscape of pairs of such cross-inhibitory "master regulators" should therefore be of particular interest in understanding both normal development and induced cell fate reprogramming, and can be greatly aided by detailed quantitative characterization of the interactions between these regulators.
Bistable network model
where variables x and y represent the concentrations of the two gene products, and parameters B X and B Y denote the basal (constitutive) expression rates of genes x and y, respectively. The parameters fold YX and fold XY represent the rate constants, and K DYX and K DXY the effective affinity constants, for the suppressive effects of gene y on gene x, and of gene x on gene y, respectively. The mutual suppression of the two genes is quantified by the Hill-coefficient n H (the interaction is ultrasensitive for values of n H > 1). Parameters deg X and deg Y represent the first-order degradation rate constants for the two gene products x and y, respectively. For this simplified model we used dimensionless parameters with the following values: fold YX = fold XY = 2; K DYX = 0.7; K DXY = 0.5; B X = B Y = 0.2; deg X = deg Y = 1; n H = 4. These values were tuned to ensure bistable switching behavior in the model.
Tristable network model
where the new parameters fold XX and fold YY represent the rate constants, and K DXX and K XYY the effective affinity constants, for the positive autoregulation of genes x and y, respectively. The default parameter values chosen to ensure three robust stable states in this model were as follows: fold XX = fold YY = fold YX = fold XY = 10; K DXX = K DYY = K DYX = K DXY = 4; B X = B Y = 0; deg X = deg Y = 1; n H = 4. This system has been modeled previously [42, 48, 49] in the context of mutual inhibition of the transcription factors PU.1 and GATA1 in common myeloid progenitor (CMP) cells, which gives rise to either bipotential granulocyte/macrophage progenitor (GMP) cells or megakaryocyte/erythroid progenitor (MEP) cells.
To evaluate the change in the quasi-potential along each trajectory in x-y phase space by numerical integration, the initial level of the quasi-potential at time t = 0 at the origin of the trajectory was arbitrarily set to zero (the same initial quasi-potential level was used for all trajectories so that the drop in the quasi-potential along each trajectory could be compared and used as a basis for alignment of multiple trajectories along a basin of attraction).
Thereafter, at each time step:
The rates and were updated to the current value of x and y according to Eqs. 8 and 9.
The above steps were repeated until the quasi-potential V q converged to a minimum (decided by a pre-set tolerance). Multiple trajectories thus obtained were aligned into basins of attraction according to the process described in the main text. The quasi-potential surface was then derived by linear interpolation among the aligned trajectories.
Software platforms used
The deterministic models were implemented and simulated on the MATLAB® (R2009a, The MathWorks, Inc., Natick, MA) platform, while the BioNetS program , based on the Gillespie algorithm [51, 52], was used for stochastic simulations. All graphics were rendered on MATLAB®.
Visualization of stochastic simulation results
The stochastically simulated "cells" (i.e. individual realizations of the stochastic network model) were overlaid on the quasi-potential surface at the x and y values predicted for each cell. Since stochastic simulations yield integral values, we added a small random "deviation" term [= (rand*0.5) where rand is a MATLAB® function that draws pseudorandom values from the standard uniform distribution on the open interval (0,1)] to each simulated x and y value to visualize multiple cells situated at the same point in x-y phase space. The appropriate "elevation" for each cell on the quasi-potential surface was calculated by linear interpolation between the two points on the deterministic trajectories "closest to" the location of the cell in x-y phase space. Source code for the model in MATLAB® format is appended in Additional File 2.
We thank J. D. Schroeter, C. Woods, R. B. Conolly, H. J. Clewell III, J.E. Trosko, M. Thattai, J. M. Haugh and B. Howell for critical discussions and reading of the manuscript. This work was supported by the Superfund Research Program of the U.S. National Institute of Environmental Health Sciences, and by the American Chemistry Council.
- Gilbert S: Epigenetic landscaping: Waddington's use of cell fate bifurcation diagrams. Biology and Philosophy. 1991, 6 (2): 135-154. 10.1007/BF02426835.View Article
- Slack JMW: Conrad Hal Waddington: the last renaissance biologist?. Nat Rev Genet. 2002, 3 (11): 889-895. 10.1038/nrg933View ArticlePubMed
- Waddington CH: An introduction to modern genetics. 1939, London: George Allen & Unwin
- Waddington CH: Organisers and Genes. 1940, Cambridge, UK: Cambridge University Press
- Waddington CH: The Strategy of the Genes. 1957, London: George Allen & Unwin
- Kauffman SA: Metabolic Stability and Epigenesis in Randomly Constructed Genetic Nets. J Theor Biol. 1969, 22 (3): 437-467. 10.1016/0022-5193(69)90015-0View ArticlePubMed
- Enver T, Pera M, Peterson C, Andrews PW: Stem Cell States, Fates, and the Rules of Attraction. Cell Stem Cell. 2009, 4 (5): 387-397. 10.1016/j.stem.2009.04.011View ArticlePubMed
- Huang S: Reprogramming cell fates: reconciling rarity with robustness. Bioessays. 2009, 31 (5): 546-560. 10.1002/bies.200800189View ArticlePubMed
- MacArthur BD, Ma'ayan A, Lemischka IR: Systems biology of stem cell fate and cellular reprogramming. Nat Rev Mol Cell Biol. 2009, 10 (10): 672-681.PubMed CentralPubMed
- Thom R: Structural Stability and Morphogenesis. 1976, Reading, MA: W. A. Benjamin, Inc
- Thom R: An Inventory of Waddingtonian Concepts. Theorical Biology: Epigenetic and Evolutionary Order from Complex Systems. Edited by: Goodwin B, Saunders P. 1989, 1-7. Edinburgh: Edinburgh University Press
- Wang J, Xu L, Wang E, Huang S: The Potential Landscape of Genetic Circuits Imposes the Arrow of Time in Stem Cell Differentiation. Biophys J. 2010, 99 (1): 29-39. 10.1016/j.bpj.2010.03.058PubMed CentralView ArticlePubMed
- Zhou JX, Huang S: Understanding gene circuits at cell-fate branch points for rational cell reprogramming. Trends Genet. 2011, 27 (2): 55-62. 10.1016/j.tig.2010.11.002View ArticlePubMed
- Kim KY, Wang J: Potential energy landscape and robustness of a gene regulatory network: Toggle switch. PLoS Comput Biol. 2007, 3 (3): 565-577.View Article
- Cao Y, Liang J: Optimal enumeration of state space of finitely buffered stochastic molecular networks and exact computation of steady state landscape probability. BMC Syst Biol. 2008, 2:
- Zhu XM, Yin L, Hood L, Ao P: Robustness, stability and efficiency of phage λ genetic switch: Dynamical structure analysis. J Bioinform Comput Biol. 2004, 2 (4): 785-817. 10.1142/S0219720004000946View ArticlePubMed
- Zhu XM, Yin L, Hood L, Ao P: Calculating biological behaviors of epigenetic states in the phage λ life cycle. Funct Integr Genomics. 2004, 4 (3): 188-195.View ArticlePubMed
- Lapidus S, Han B, Wang J: Intrinsic noise, dissipation cost, and robustness of cellular networks: The underlying energy landscape of MAPK signal transduction. Proc Natl Acad Sci USA. 2008, 105 (16): 6039-6044. 10.1073/pnas.0708708105PubMed CentralView ArticlePubMed
- Wang J, Xu L, Wang EK: Potential landscape and flux framework of nonequilibrium networks: Robustness, dissipation, and coherence of biochemical oscillations. Proc Natl Acad Sci USA. 2008, 105 (34): 12271-12276. 10.1073/pnas.0800579105PubMed CentralView ArticlePubMed
- Li C, Wang E, Wang J: Potential Landscape and Probabilistic Flux of a Predator Prey Network. PLoS ONE. 2011, 6 (3): e17888- 10.1371/journal.pone.0017888PubMed CentralView ArticlePubMed
- Takahashi K, Yamanaka S: Induction of pluripotent stem cells from mouse embryonic and adult fibroblast cultures by defined factors. Cell. 2006, 126 (4): 663-676. 10.1016/j.cell.2006.07.024View ArticlePubMed
- Takahashi K, Tanabe K, Ohnuki M, Narita M, Ichisaka T, Tomoda K, Yamanaka S: Induction of pluripotent stem cells from adult human fibroblasts by defined factors. Cell. 2007, 131 (5): 861-872. 10.1016/j.cell.2007.11.019View ArticlePubMed
- Wernig M, Meissner A, Foreman R, Brambrink T, Ku MC, Hochedlinger K, Bernstein BE, Jaenisch R: In vitro reprogramming of fibroblasts into a pluripotent ES-cell-like state. Nature. 2007, 448 (7151): 318-U312. 10.1038/nature05944View ArticlePubMed
- Maherali N, Sridharan R, Xie W, Utikal J, Eminli S, Arnold K, Stadtfeld M, Yachechko R, Tchieu J, Jaenisch R, Plath K, Hochedlinger K: Directly reprogrammed fibroblasts show global epigenetic remodeling and widespread tissue contribution. Cell Stem Cell. 2007, 1 (1): 55-70. 10.1016/j.stem.2007.05.014View ArticlePubMed
- Yu JY, Vodyanik MA, Smuga-Otto K, Antosiewicz-Bourget J, Frane JL, Tian S, Nie J, Jonsdottir GA, Ruotti V, Stewart R, Slukvin II, Thomson JA: Induced pluripotent stem cell lines derived from human somatic cells. Science. 2007, 318 (5858): 1917-1920. 10.1126/science.1151526View ArticlePubMed
- Dimos JT, Rodolfa KT, Niakan KK, Weisenthal LM, Mitsumoto H, Chung W, Croft GF, Saphier G, Leibel R, Goland R, Wichterle H, Henderson CE, Eggan K: Induced pluripotent stem cells generated from patients with ALS can be differentiated into motor neurons. Science. 2008, 321 (5893): 1218-1221. 10.1126/science.1158799View ArticlePubMed
- Hanna J, Markoulaki S, Schorderet P, Carey BW, Beard C, Wernig M, Creyghton MP, Steine EJ, Cassady JP, Foreman R, Lengner CJ, Dausman JA, Jaenisch R: Direct reprogramming of terminally differentiated mature B lymphocytes to pluripotency. Cell. 2008, 133 (2): 250-264. 10.1016/j.cell.2008.03.028PubMed CentralView ArticlePubMed
- Shen CN, Slack JMW, Tosh D: Molecular basis of transdifferentiation of pancreas to liver. Nat Cell Biol. 2000, 2 (12): 879-887. 10.1038/35046522View ArticlePubMed
- Xie HF, Ye M, Feng R, Graf T: Stepwise reprogramming of B cells into macrophages. Cell. 2004, 117 (5): 663-676. 10.1016/S0092-8674(04)00419-2View ArticlePubMed
- Cobaleda C, Jochum W, Busslinger M: Conversion of mature B cells into T cells by dedifferentiation to uncommitted progenitors. Nature. 2007, 449 (7161): 473-477. 10.1038/nature06159View ArticlePubMed
- Zhou Q, Brown J, Kanarek A, Rajagopal J, Melton DA: In vivo reprogramming of adult pancreatic exocrine cells to beta-cells. Nature. 2008, 455 (7213): 627-U630. 10.1038/nature07314View ArticlePubMed
- Vierbuchen T, Ostermeier A, Pang ZP, Kokubu Y, Sudhof TC, Wernig M: Direct conversion of fibroblasts to functional neurons by defined factors. Nature. 2010, 463 (7284): 1035-U1050. 10.1038/nature08797PubMed CentralView ArticlePubMed
- Ieda M, Fu JD, Delgado-Olguin P, Vedantham V, Hayashi Y, Bruneau BG, Srivastava D: Direct Reprogramming of Fibroblasts into Functional Cardiomyocytes by Defined Factors. Cell. 2010, 142 (3): 375-386. 10.1016/j.cell.2010.07.002PubMed CentralView ArticlePubMed
- Tursun B, Patel T, Kratsios P, Hobert O: Direct conversion of C. elegans germ cells into specific neuron types. Science. 2011, 331 (6015): 304-308. 10.1126/science.1199082PubMed CentralView ArticlePubMed
- Zhou Q, Melton DA: Extreme Makeover: Converting One Cell into Another. Cell Stem Cell. 2008, 3 (4): 382-388. 10.1016/j.stem.2008.09.015View ArticlePubMed
- Yamanaka S: Elite and stochastic models for induced pluripotent stem cell generation. Nature. 2009, 460 (7251): 49-52. 10.1038/nature08180View ArticlePubMed
- Graf T, Enver T: Forcing cells to change lineages. Nature. 2009, 462 (7273): 587-594. 10.1038/nature08533View ArticlePubMed
- Huang S: Non-genetic heterogeneity of cells in development: more than just noise. Development. 2009, 136 (23): 3853-3862. 10.1242/dev.035139PubMed CentralView ArticlePubMed
- Hochedlinger K, Plath K: Epigenetic reprogramming and induced pluripotency. Development. 2009, 136 (4): 509-523. 10.1242/dev.020867PubMed CentralView ArticlePubMed
- Gardner TS, Cantor CR, Collins JJ: Construction of a genetic toggle switch in Escherichia coli. Nature. 2000, 403 (6767): 339-342. 10.1038/35002131View ArticlePubMed
- Strogatz S: Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering (Studies in Nonlinearity). 2001, Cambridge, MA: Perseus Books Group
- Huang S, Guo YP, May G, Enver T: Bifurcation dynamics in lineage-commitment in bipotent progenitor cells. Dev Biol. 2007, 305 (2): 695-713. 10.1016/j.ydbio.2007.02.036View ArticlePubMed
- Wilkinson DJ: Stochastic modelling for quantitative description of heterogeneous biological systems. Nat Rev Genet. 2009, 10 (2): 122-133. 10.1038/nrg2509View ArticlePubMed
- Kaern M, Elston TC, Blake WJ, Collins JJ: Stochasticity in gene expression: From theories to phenotypes. Nat Rev Genet. 2005, 6 (6): 451-464. 10.1038/nrg1615View ArticlePubMed
- Roeder I, Radtke F: Stem cell biology meets systems biology. Development. 2009, 136 (21): 3525-3530. 10.1242/dev.040758View ArticlePubMed
- Hasty J, Pradines J, Dolnik M, Collins JJ: Noise-based switches and amplifiers for gene expression. Proc Natl Acad Sci USA. 2000, 97 (5): 2075-2080. 10.1073/pnas.040411297PubMed CentralView ArticlePubMed
- Acar M, Becskei A, van Oudenaarden A: Enhancement of cellular memory by reducing stochastic transitions. Nature. 2005, 435 (7039): 228-232. 10.1038/nature03524View ArticlePubMed
- Roeder I, Glauche I: Towards an understanding of lineage specification in hematopoietic stem cells: A mathematical model for the interaction of transcription factors GATA-1 and PU.1. J Theor Biol. 2006, 241 (4): 852-865.View ArticlePubMed
- Chickarmane V, Enver T, Peterson C: Computational Modeling of the Hematopoietic Erythroid-Myeloid Switch Reveals Insights into Cooperativity, Priming, and Irreversibility. PLoS Comput Biol. 2009, 5 (1): e1000268- 10.1371/journal.pcbi.1000268PubMed CentralView ArticlePubMed
- Adalsteinsson D, McMillen D, Elston TC: Biochemical Network Stochastic Simulator (BioNetS): software for stochastic modeling of biochemical networks. BMC Bioinformatics. 2004, 5: 24- 10.1186/1471-2105-5-24PubMed CentralView ArticlePubMed
- Gillespie DT: A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. Journal of Computational Physics. 1976, 22 (4): 403-434. 10.1016/0021-9991(76)90041-3.View Article
- Gillespie DT: Exact Stochastic Simulation of Coupled Chemical-Reactions. J Phys Chem. 1977, 81 (25): 2340-2361. 10.1021/j100540a008.View Article
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.