| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Department of Bioengineering and Institute for Mathematical Sciences, Imperial College London, London, United Kingdom
Correspondence: Address reprint requests to M. Barahona, Tel.: 44-207-594-5189; E-mail: m.barahona{at}imperial.ac.uk.
 |
ABSTRACT |
|---|
|
INTRODUCTION |
|---|
–10). The inherent stochasticity of these biochemical processes is due to the low number of molecules involved in the reactions (11), amplified in some cases by the proximity to a critical point of the dynamics of the system (12). Low copy numbers lead to intrinsic fluctuations and to the breakdown of models based on differential equations (13). A number of stochastic models of gene regulation (11,14–21) and signaling (22–24) have been formulated and analyzed numerically with the aim of identifying the sources of randomness, and how noise is controlled and harnessed in the cellular environment.
The starting point for such stochastic descriptions is the master equation (ME), a conservation equation that gives the time evolution of the probability distribution of the state of the system. For a discrete state space, it can be written as (13)
 | (1) |
). This work has been extended by a number of authors (26,27). Equation 1 is usually referred to as the chemical master equation (CME) when it describes stochastic versions of law-of-mass-action systems. In this article, we extend the notation and apply the term CME to compound reaction rates (e.g., Michaelis-Menten) as well.
Although theoretically rigorous, there are very few systems for which the ME has been solved analytically. Direct numerical integration is often very difficult, due to the fact that the state space grows very rapidly, and also to the stiffness of the equations. One way to overcome these issues is to employ some kind of approximation scheme. In certain limits, one can consider continuum approximations leading to partial differential equations, such as van Kampen's linear-noise approximation and Fokker-Planck equations (13). However, these approximations disregard the discreteness of the state of the system, and can therefore give rise to significant deviations when the number of molecules is very small, as can be the case for gene regulatory networks (GRNs).
A different approach for dealing with the CME numerically is provided by the widely used stochastic simulation algorithm (SSA) by Gillespie (28). The SSA is a kinetic Monte Carlo algorithm, rigorously derived from the same assumptions as the CME, which gives realizations of the trajectory of a given CME. The SSA follows an idea that can be traced back to Doob (29) and provides an exact procedure for the kinetic sampling of the CME (30) in the sense that we obtain an unbiased and convergent estimate of the solution of the CME. Due to its biological applications, there has been considerable interest in the SSA in recent years and several extensions have been proposed (e.g., the explicit introduction of space (31,32) and time-delays (14)), as well as algorithmic speed-up improvements (33,34).
In many experiments of interest (3,7,9,29), and in related theoretical analyses (11,16,18,21–23,35,36), the system at hand is assumed to have reached the stationary distribution, i.e., the left-hand side of Eq. 1 is zero. Consider the ME (Eq. 1) in operator form,
 | (2) |
in Markov theory) could in principle be obtained from the first left eigenvector of the Q-matrix (37). However, this approach is impractical due to the curse of dimensionality: the run-time and memory requirements scale exponentially with the number of types of molecules. For small state spaces, one can still consider a (truncated) finite subspace to obtain an approximate convergent solution. Our own accurate version of this approximate eigenvector method is used in this article to check the accuracy of our DCFTP-SSA algorithm when an analytical expression for the stationary distribution is not known.
To sample from the stationary solution of Eq. 1, one would have to run the SSA for an infinite time. Of course, the most common practical solution is to run the algorithm repeatedly for a very long time and hope that the system has reached the stationary distribution when the run is stopped (15). The lack of a termination certificate can lead to computational inefficiency if the runs are longer than necessary or, more importantly, to missampling from the wrong distribution if the runs are not long enough.
In this article, we present an algorithm (DCFTP-SSA) that guarantees perfect sampling of the stationary solution of the CME for the general class of biochemical networks formed by the interconnection of generalized Hill- and Monod-type reactions, the canonical models for GRNs and enzymatic networks. The DCFTP-SSA builds on the standard SSA and considers it in the light of Markov chain Monte Carlo algorithms, specifically within the dominated-coupling-from-the-past (DCFTP) framework introduced by Wilson (38) and Kendall (39). Henceforth, we introduce the algorithm and its properties in detail through the analysis of the building blocks of GRNs. We then apply the DCFTP-SSA to the study of two systems of experimental interest in synthetic biology: the toggle switch (6) and the repressilator (8).
|
THE DOMINATED COUPLING FROM THE PAST STOCHASTIC SIMULATION ALGORITHM (DCFTP-SSA) |
|---|
). The major issue in these schemes is when to stop the Markov chain. This problem was addressed by Propp and Wilson with their coupling from the past (CFTP) algorithm (38), a celebrated example of what are commonly referred to as perfect sampling algorithms (40). The version used here is the dominated coupling from the past (DCFTP) algorithm, the extension introduced by Kendall (39) to study continuous-time Markov processes on an unbounded state space. We now give a brief introduction to the algorithm. In principle, a Markov chain would have to be run for an infinite time to reach its stationary distribution. The CFTP algorithm, however, recasts the problem as a procedure in which the running time becomes a random variable but a certificate is issued if the stationary distribution has been reached. This is signaled by the coalescence of relevant coupled Markov chains. Markov chains are coupled if they have the same transition rules and use the same sequence of random numbers for their realization. Coupled Markov chains are said to coalesce when they meet. Clearly, coupled chains with the same transition rule (but started from different initial states) will have the same values for all future times after the coalescence time Tc. Propp and Wilson proposed to use Markov chains coupled from the past to ensure that the whole state space of initial conditions maps to the same state at present. By extending sufficiently far back into the past and using a fixed sequence of random numbers, we will eventually reach a time from which all paths map to the same state at t = 0. This condition is equivalent to the stationary distribution, since the starting condition is irrelevant for the current state.
In many systems, the state space is very large, making it infeasible to monitor all paths and their coalescence. However, the situation is much simpler if the state space is partially ordered and the partial ordering is maintained by the transition rules. This is the case if we have a monotone transition rule 
, such that (x, R) 
(y, R) if x 
y, where denotes the partial ordering, x and y are states and R is a random number used to determine the transition. The algorithm also works for anti-monotone transition rules: (x, R) 
(y, R) if x y (41). For a state space to be partially ordered, it must be reflexive, antisymmetric, and transitive. If the monotonicity and partial ordering hold, we only need to monitor two paths: an upper path U and a lower path L, since all other paths lie in between and will coalesce when U and L have done so. This sandwiching property is illustrated in Fig. 1. The guaranteed sampling of the CFTP scheme comes with a tradeoff: the run time is unbounded, since Tc is itself a random variable, but it is finite almost surely. If the run is prematurely terminated by an impatient user, the sample will be biased.
|
). The DCFTP relies on finding a reversible dominating process D, which bounds the original process from above and for which the stationary distribution is known. To create a dominating process, we exploit the reversibility and stationarity of D: start the chain 
at t = 0 from the stationary distribution and evolve it until t = T; then use 
as the dominating process on the interval [–T, 0]. It then follows that all chains of the original dominated process started at t = –
will be less than or equal to D at time –T and, consequently, chains coupled to D started from a state U–T D–T can be interpreted as random realizations of the original (dominated) process started from t = –.
Gillespie's SSA implements a continuous-time Markov process. As such, it can be reformulated in the CFTP framework if the state space is partially ordered and the propensity functions are (anti-)monotone. Indeed, one can show that a partial ordering based on the number of molecules in each species is maintained for systems of unimolecular reactions with (anti-) monotone propensity functions, like those generated by Hill or Monod birth rates (M. Hemberg and M. Barahona, unpublished). If such a system has a fixed (or bounded) number of molecules, the upper path is trivial and the standard CFTP can be directly applied to the SSA. However, in most situations, the number of molecules is unbounded and we must use the DCFTP algorithm (39).
Based on the discussion above, a perfect sampling DCFTP-SSA can be developed if two prerequisites are fulfilled. Firstly, the Markov process defined by the dominating process must be reversible. This means that, for every reaction in the system, another reaction must exist with products and reactants exchanged. Note that this is different to requiring that the underlying reactions be reversible and it typically means that, for every reaction creating a molecule, a degradation reaction must exist. Secondly, we must find a dominating process for the particular CME with a known stationary distribution. We will show below that this requirement can be fulfilled for genetic and enzymatic networks formed by the interconnection of Hill-type, Monod-type, and linear enzymatic unimolecular reactions. This is based on the fact that the CME of networks of Hill-type or Monod-type elements can be bounded from above by processes based on networks of first-order reactions, for which the stationary distribution is known to be the multivariate Poisson (42).
A brief outline of the DCFTP-SSA is as follows:
with known stationary distribution forward in time from t = 0 until t = T. to obtain the dominating process D such that 
from t = –T to t = 0. The main feature of the DCFTP-SSA is that it provides a certificate for correct sampling from the stationary distribution. Therefore, we can use the DCFTP-SSA to sample, numerically, the stationary distribution of the system with guaranteed Monte Carlo accuracy. This can be of importance in high-dimensional systems with complicated landscapes where the SSA simulations can converge slowly. In the first two sections below, we test the application of the DCFTP-SSA to the sampling of stationary distributions. In those cases, we only use the final value of each run as a sample, thus ensuring that the samples are independent.
In addition, we can use the algorithm to discard the burn-in period in stochastic simulations by running the DCFTP-SSA until a guaranteed sample from the stationary distribution is obtained and then using it to initiate a run of the ordinary SSA from time t = 0 onwards. This is important for the numerical study of stationary properties of systems with high variability, e.g., with underlying oscillatory or excitatory behavior. In the final section, we present examples of these applications to the toggle switch and the repressilator.
Detailed application to the first-order reaction
To illustrate the DCFTP-SSA, we study in detail the simple one-dimensional first-order reaction, for which a full time-dependent analytical solution of the CME is available. This allows us to study the convergence of the algorithm, as compared to the standard SSA, and to explore the difference between two sources of error, which often get entangled: finite sampling error and the missampling error due to the fact that we may not have reached the true stationary distribution.
Consider the simple first-order reaction, which has been used as a very simplified description of transcription and translation (2,21),
 | (3) |
) and degraded to a sink (). The corresponding CME is
 | (4) |
and –1 are the step operators defined by van Kampen (31): f(j) = f(j + 1) and –1f(j) = f(j – 1) for a function f(j).
Equation 4 is one of a few CMEs for which the analytical expression of the stationary solution is known (13,42): it is the Poisson distribution with parameter 
= k. More importantly for our purposes, one can obtain the full time-dependent solution of the CME (see Appendix 1). For the usual initial condition with 0 molecules, the time-dependent distribution P(t) is given by (see Fig. 2, inset)
 | (5) |
.
|
. As N is increased, with fixed tSSA, we sample with increasing Monte Carlo accuracy from P(tSSA) given by Eq. 5 but not from the true stationary distribution P*. We therefore reach an error floor that cannot be broken. Using the analytical expression (Eq. 5), the asymptotic error floor 
E*(tSSA) is shown to be
 | (6) |
and I0(x) is the modified Bessel function of the first kind. Fig. 2 shows that the leveling-off of the Euclidean error for PSSA(N, tSSA) is explained by the error floors calculated from Eq. 6.
This source of error is eliminated in a DCFTP-SSA formulation of this process. The key step is to find a reversible dominating process with known stationary distribution. In this example, it is clear that the best choice is to use the process (Eq. 4) itself, i.e., the U and D chains in Fig. 1 will be identical. Fig. 2 shows the Euclidean error for the DCFTP-SSA sample distribution, which shows no flooring and the expected N–1/2 scaling with the number of Monte Carlo runs (43). Equivalent results are obtained with other error measures for the sampled distribution, e.g., the Kullback-Leibler and Kolmogorov distances, and the 
2-goodness of fit.
|
THE BUILDING BLOCKS: NETWORKS AND NONLINEAR ELEMENTS |
|---|
Application to networks: simple gene model of two coupled first-order reactions
To showcase how to apply the DCFTP-SSA to networks of several species, we use a widely studied, simplified model of gene expression (17,19,21). In its simplest form, the model consists of two types of molecules: mRNA (M) and proteins (R), which are produced and depleted at constant rates. In addition, the mRNA catalyzes protein production through a linear enzymatic reaction:
 |
The deterministic description is given by
 | (7) |
 | (8) |
In fact, Eq. 8 has been solved: the stationary state of the CME of any general network of linear reactions is a multivariate Poisson distribution (20,39). This means that, similarly to the one-dimensional first-order reaction, the DCFTP-SSA is directly applicable to Eq. 8, since the reaction network is reversible and the multivariate Poisson stationary solution is itself a dominating process for the system. The existence of an analytical solution allows us to check how the accuracy of the DCFTP-SSA depends on the dimensionality of the system. Fig. 3 (inset) shows that the marginal distribution (for M) converges like the one-dimensional reaction in Fig. 2, with similar error floors for the standard SSA. Note, however, that the joint distribution (Fig. 3) exhibits slower convergence of the standard SSA due to the higher dimensionality of the system, indicating the need for longer SSA runs to avoid missampling from the true distribution. The DCFTP-SSA shows N–1/2 Monte Carlo scaling with no error floors regardless of the dimensionality of the system.
|
,14,24). This model incorporates negative feedback, whereby the rate of creation of new molecules decreases as their concentration increases:
 | (9) |
Here, f(A;k, 
) = k/(
+ A) is the state-dependent reaction rate, k is the renormalized reaction constant, and is the cooperativity factor (44). Equation 9 can be derived from elementary law-of-mass-action kinetics through elimination of variables (18). The corresponding CME is given by
 | (10) |
It is usual to fix = 1 and we do so in the following. For the particular case = 1, Eq. 9 reduces to the familiar Michaelis-Menten equation, for which we can obtain the following analytical expression for the stationary distribution (see Appendix B): P*j = ckj/(j!)2, where 
is a normalization constant.
The Monod reaction is used to model gene upregulation, as referred in particular to the autocatalysis common in many eukaryotes (3,45,46). The reaction can be written as
 | (11) |
) = kA/( + A) is the state-dependent reaction rate that encapsulates the positive feedback. The corresponding CME is given by
 | (12) |
Again, we fix = 1 in the following. The stationary distribution of the Monod CME (Eq. 12) with = 1 can be obtained analytically (see Appendix B): P*j = ckj/j!, with normalization constant c = 1/(ek + k – 1).
For general , no explicit solution for the stationary distribution of the CMEs Eq. 10 and Eq. 12 is known, and numerical methods are therefore necessary. Fortunately, the DCFTP-SSA is directly applicable to the reversible Hill and Monod CMEs: the first-order reaction, studied in detail in the previous section, can be used to produce a dominating process for both processes. To see this, note that Eq. 10 and Eq. 12 have birth rates that are bounded from above by k at all times. It is thus straightforward to make sure that the upper and lower paths have rates that fulfill the requirements of the algorithm. As an illustration, Fig. 1 shows one DCFTP-SSA run for the Hill CME (Eq. 10), where D corresponds to the dominating first-order linear process (Eq. 4) and U and L are the coupled Hill processes (39,40). Our detailed simulations of the Hill and Monod CMEs (data not shown) show N–1/2 convergence of the DCFTP-SSA unaffected by the missampling errors that can appear when using the standard SSA with short stopping times.
|
APPLICATION TO NONLINEAR MODELS OF GENE REGULATORY NETWORKS |
|---|
), with a bimodal stationary distribution; and the repressilator (6), a simple synthetic oscillator.
Toggle switch: two coupled Hill equations
An important characteristic of biochemical networks is the possibility of multistability (1,8,15,16,18,21,47,48), as exemplified by the toggle switch designed by Gardner et al. by combining two mutually repressing genes (8,11,15,16,18). This leads to a bistable system with a sharp switching threshold:
 |
A deterministic model can be derived using two coupled Hill equations:
 | (13) |
, ß > 1 and reaction constants kA and kB, Eq. 13 has two stable fixed points (8).
The corresponding CME is given by
 | (14) |
Based on our discussion above, it is easy to see that linear process with constant birth rates kA and kB will be a dominating process for Eq. 14. Since the Hill function is anti-monotone, we must use a crossover scheme for the updates to make sure that all initial conditions are sandwiched between the upper and lower paths, as described in Häggström and Nelander (41). We have applied the DCFTP-SSA to Eq. 14 to obtain the stationary distribution of the system, which is clearly bimodal (Fig. 4). Because there is no analytical solution in this case, we have checked the accuracy and convergence of the DCFTP-SSA against the approximate eigenvector method with excellent results (not shown).
|
In their experiments, Gardner et al. (21) were able to induce faster A# B# switching by reducing the downregulation capability of A in the cellular environment. This can be mimicked by a modified CME in which the Aß-term in Eq. 14 is removed. If we run the SSA of this modified CME starting from the stationary distribution of the original system (Eq. 14), the escape time 
becomes approximately two orders-of-magnitude smaller, in broad agreement with the experiment. In fact, it can be shown that the stationary distribution of the modified CME becomes unimodal, with no observable peak for species A.
Repressilator: six-dimensional system of linear enzymatic and Hill equations
The repressilator, a synthetic system of transcriptional regulators, is perhaps the simplest biochemical oscillator that has been implemented experimentally (6). It consists of three genes in a loop, in which the expression of one gene is inhibited by the product of another gene in succession (Fig. 5 a):
 | (15) |
|
The corresponding CME is given by
 | (16) |
and i mod 3. We have simulated this six-dimensional system using the crossover scheme described above: we use the DCFTP-SSA to discard the burn-in period, i.e., the time it takes for the Markov process to reach the stationary distribution, such that the initial conditions for the SSA runs are sampled from the stationary distribution. As a measure of the stationarity of our DCFTP-SSA initialized sampling, we have checked its ergodicity by establishing through an F-test that the average period and amplitude of several short simulations are statistically indistinguishable from those obtained from one long simulation.
As can be seen in the sample trajectory in Fig. 5 b, the stochastic system is extremely noisy but it retains the overall oscillation of the genes in succession. The oscillations do not die out even for very long simulations and we have collected statistics of such long runs. Fig. 5 c shows a robust feature of the oscillator: the distribution of the period is well fitted by a generalized Gamma distribution, which is characteristic of excitatory systems with a refractory period. The oscillatory behavior is present for a wide range of the parameters kM, kR, and dM (44), which could be modified experimentally. Fig. 5 d summarizes our investigation of the robustness of the repressilator to parametric variation. Note the good agreement of the calculated mean period with the corresponding deterministic values. The simulations show that changes in parameters kM and kR produce a similar, small effect on the period. The system is, however, more sensitive to parametric changes in dM. The evaluation of how these parametric variations could be used for the design of tunable and reliable oscillators will be the subject of further study (49).
|
DISCUSSION |
|---|
The guaranteed sampling provided by the DCFTP-SSA comes at an extra computational cost. In general, the coalescence time will depend on the topology of the network and the form of the reactions. However, the CPU overhead is in no way prohibitive and the DCFTP-SSA can be run on an ordinary desktop computer for the systems presented in this article. An important theoretical feature of the DCFTP-SSA is the fact that its runtime is not bounded. As our simulations show, this feature does not seem to impose practical restrictions on the algorithm. If necessary, however, this issue can be resolved by using FMMR (50), another perfect sampling algorithm, which dispenses with variable stopping times by running the Markov chains for a fixed time and using rejection sampling to discard those simulations that have not reached the stationary distribution.
In its current form, the DCFTP-SSA can be applied to reversible systems of linear and nonlinear (anti-)monotonic unimolecular reactions for which there exists a dominating process with known stationary distribution. We note that the algorithm is applicable not only to networks of linear enzymatic, Hill and Monod birth reactions of standard form, but also to reactions whose rate equation is a rational function that can be dominated by a linear process. This condition is equivalent to requiring that the polynomial in the denominator of the propensity function have a higher degree than the polynomial in the numerator. These types of rational functions are frequently obtained as compound unimolecular rate laws to represent more complex enzymatic (45) or gene regulatory reactions (51). An example of such functions is the model for the -phage lysis-lysogeny switch (11,16), which we have also simulated with the DCFTP-SSA (results not shown). Bimolecular and two-substrate enzymatic reactions are not included in this group and extending the DCFTP or CFTP to provide perfect sampling for the generic SSA of second-order reactions will entail further research.
The DCFTP-SSA introduced here can be viewed as a reformulation and extension of the Gillespie algorithm that provides guaranteed sampling from stationarity for a generic class of systems relevant in GRNs, thereby removing a source of uncontrolled uncertainty from the simulations. By eliminating these extraneous sources of error, comparisons between the predictions of different stochastic models at stationarity can be performed more meaningfully. In the case of the repressilator, we provided an accurate numerical characterization of the stationary distribution for the period of this stochastic oscillator. Because stationarity is guaranteed, the characteristics of this distribution can be correlated with changes in the parameter values (as shown in Fig. 5 d) or compared with those of other genetic oscillators. This could be a useful tool for the analysis and design of synthetic transcriptional oscillators and will be the subject of future research.
|
APPENDIX A: TIME-DEPENDENT SOLUTION OF THE MASTER EQUATION FOR THE ONE-DIMENSIONAL FIRST-ORDER REACTION |
|---|
-queue with Poisson arrival and exponential service time (37). Associated with Eq. 4, there exists a family of orthogonal polynomials 
with the three-term recurrence relation (52):
 | (17) |
0(x) = 1 and j(x) = 0 for j < 0. Karlin and McGregor showed that the recurrence relation (Eq. 17) leads to a spectral representation for Pji(t), the transition probability of going from state i to state j in time t,
 | (18) |
(x) is a positive measure on x and Pj* is the stationary distribution.
The recurrence relation (Eq. 17) is satisfied by the Charlier polynomials with (x) equal to the Poisson distribution (52), i.e., the Charlier polynomials are a family of polynomials which are orthogonal with respect to the Poisson measure on a discrete lattice on the nonnegative integers (53). The Charlier polynomial of order i, parameter a, and argument x is denoted Ci(x;a).
Using the following bilinear generating form of the Charlier polynomials (54)
 | (19) |
Lee (55) showed that Eq. 18 can be simplified to give
 | (20) |
The important special case with initial condition given by a 
-distribution at the origin can be simplified even further. If Pj(0) = (0), then C0(j;a) = 1 and Eq. 20 becomes
 | (21) |
 | (22) |
|
APPENDIX B: STATIONARY SOLUTION OF THE HILL AND MONOD CMES WITH = 1
|
|---|
= 1, which is equivalent to the Michaelis-Menten reaction, can be obtained analytically. This reaction can be viewed as a nonlinear birth-death process, which allows us to write the stationary distribution Pj* as (37)
 | (23) |
i and µi are the birth and death rates of state i and c is a normalization constant. Inserting the values from the CME (Eq. 10), we obtain
 | (24) |
 |
(x) are Bessel functions. For the special case of = 1, we use a recurrence relation for Bessel functions to simplify the above result to 
Consider now the Monod CME (12) with = 1:
 |
Note that the birth rate of the state j = 0 is modified to avoid it becoming absorbing. Using the same strategy as for the Michaelis-Menten CME above, one can find the stationary distribution
 | (25) |
where 1F1 is Gauss's hypergeometric function. |
ACKNOWLEDGEMENTS |
|---|
|
FOOTNOTES |
|---|
Submitted on October 11, 2006; accepted for publication March 20, 2007.
|
REFERENCES |
|---|
2. Austin, D., M. Allen, J. McCollum, R. Dar, J. Wilgus, G. Sayler, N. Samatova, C. Cox, and M. Simpson. 2006. Gene network shaping of inherent noise spectra. Nature. 439:608–611.[CrossRef][Medline]
3. Beckskei, A., B. Séraphin, and L. Serrano. 2001. Positive feedback in eukaryotic gene networks: cell differentiation by graded to binary response conversion. EMBO J. 20:2528–2535.[CrossRef][Medline]
4. Blake, W. J., M. Kaern, C. R. Cantor, and J. Collins. 2003. Noise in eukaryotic gene expression. Nature. 422:633–637.[CrossRef][Medline]
5. Cai, L., N. Friedman, and X. S. Xie. 2006. Stochastic protein expression in individual cells at the single molecule level. Nature. 440:358–362.[CrossRef][Medline]
6. Elowitz, M. B., and S. Leibler. 2000. A synthetic oscillatory network of transcriptional regulators. Nature. 403:335–339.[CrossRef][Medline]
7. Fung, E., W. W. Wong, J. K. Suen, T. Butler, S. G. Lee, and J. C. Liao. 2005. A synthetic gene-metabolic oscillator. Nature. 435:118–122.[CrossRef][Medline]
8. Gardner, T. S., C. R. Cantor, and J. J. Collins. 2000. Construction of a genetic toggle switch in Escherichia coli. Nature. 403:339–343.[CrossRef][Medline]
9. Pedraza, J. M., and A. van Oudenaarden. 2005. Noise propagation in gene networks. Science. 307:1965–1969.
10. Volfson, D., J. Marciniak, W. J. Blake, N. Ostroff, L. S. Tsimring, and J. Hasty. 2006. Origins of extrinsic variability in eukaryotic gene expression. Nature. 439:861–864.[CrossRef][Medline]
11. McAdams, H. H., and A. Arkin. 1997. Stochastic mechanisms in gene expression. Proc. Natl. Acad. Sci. USA. 94:814–819.
12. Elf, J., J. Paulsson, O. G. Berg, and M. Ehrenberg. 2003. Near-critical phenomena in intracellular metabolite pools. Biophys. J. 84:154–170.
13. van Kampen, N. G. 1992. Stochastic Processes in Physics and Chemistry, 2nd Ed. Elsevier, Dordrecht.
14. Bratsun, D., D. Volfson, L. S. Tsimring, and J. Hasty. 2005. Delay-induced stochastic oscillations in gene regulation. Proc. Natl. Acad. Sci. USA. 102:14593–14598.
15. Erban, R., I. G. Kevrekidis, D. Adalsteinsson, and T. C. Elston. 2006. Gene regulatory networks: a coarse-grained, equation-free approach to multiscale computation. J. Chem. Phys. 124:084106.[CrossRef][Medline]
16. Hasty, J., J. Pradines, M. Dolnik, and J. J. Collins. 2000. Noise based switches and amplifiers for gene expression. Proc. Natl. Acad. Sci. USA. 97:2075–2080.
17. Kaern, M., T. C. Elston, W. J. Blake, and J. J. Collins. 2005. Stochasticity in gene expression: from theories to phenotypes. Nature Rev. 6:451–464.[CrossRef]
18. Kepler, T. B., and T. C. Elston. 2001. Stochasticity in transcriptional regulation: origins, consequences and mathematical representations. Biophys. J. 81:3116–3136.
19. Paulsson, J. 2004. Summing up the noise in gene networks. Nature. 427:415–419.[CrossRef][Medline]
20. Swain, P. S., and M. B. Elowitz. 2002. Intrinsic and extrinsic contributions to stochasticity in gene expression. Proc. Natl. Acad. Sci. USA. 99:12795–12800.
21. Thattai, M., and A. van Oudenaarden. 2001. Intrinsic noise in gene regulatory networks. Proc. Natl. Acad. Sci. USA. 98:8614–8619.
22. Elf, J., and M. Ehrenberg. 2003. Fast evaluation of fluctuations in biochemical networks with the linear noise approximation. Genome Res. 13:2475–2484.
23. Lai, K., M. J. Robertson, and D. Schaffer. 2004. The sonic hedgehog signaling system as a bistable genetic switch. Biophys. J. 86:2748–2757.
24. Samoilov, M., S. Plyasunov, and A. P. Arkin. 2005. Stochastic amplification and signaling in enzymatic futile cycles through noise-induced bistability with oscillations. Proc. Natl. Acad. Sci. USA. 102:2310–2315.
25. Delbrück, M. 1940. Statistical fluctuations in autocatalytic reactions. J. Chem. Phys. 8:120–124.[CrossRef]
26. McQuarrie, D. 1967. Stochastic approach to chemical kinetics. J. Appl. Probab. 4:413–478.[CrossRef]
27. Paulsson, J. 2005. Models of stochastic gene expression. Phys. Life Rev. 2:157–175.[CrossRef]
28. Gillespie, D. T. 1976. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys. 22:403–434.[CrossRef]
29. Doob, J. L. 1945. Markoff chains—denumerable case. Trans. Am. Math. Soc. 58:455–473.[CrossRef]
30. Gillespie, D. T. 1992. A rigorous derivation of the chemical master equation. Physica A. 188:404–425.[CrossRef]
31. Elf, J., and M. Ehrenberg. 2004. Spontaneous separation of bi-stable biochemical systems into spatial domains of opposite phases. IEE Sys. Biol. 1:230–236.
32. Stundzia, A. B., and C. J. Lumsden. 1996. Stochastic simulation of coupled reaction-diffusion processes. J. Comput. Phys. 127:196–207.[CrossRef]
33. Cao, Y., D. T. Gillespie, and L. R. Petzold. 2006. Efficient step size selection for the 
-leaping simulation method. J. Chem. Phys. 124:044109.[CrossRef][Medline]
34. Gibson, M. A., and J. Bruck. 2000. Efficient exact stochastic simulation of chemical systems with many species and many channels. J. Phys. Chem. A. 104:1876–1899.[CrossRef]
35. Kummer, U., B. Krajnc, J. Pahle, A. K. Green, C. J. Dixon, and M. Marhl. 2005. Transition from stochastic to deterministic behavior in calcium oscillations. Biophys. J. 89:1603–1611.
36. Tomioka, R., H. Kimura, T. J. Kobayashi, and K. Aihara. 2004. Multivariate analysis of noise in genetic regulatory networks. J. Theor. Biol. 229:501–521.[CrossRef][Medline]
37. Norris, J. R. 1999. Markov Chains. Cambridge University Press, Cambridge.
38. Propp, J. G., and D. B. Wilson. 1996. Exact sampling with coupled Markov chains and applications to statistical mechanics. Random Structures Algorithms. 9:223–252.[CrossRef]
39. Kendall, W. S. 1997. Perfect simulation for spatial point processes. In Bulletin ISI, 51st Session Proceedings, Vol. 3. Istanbul.
40. Thönnes, E. 2000. A primer on perfect simulation. In Springer Lecture Notes in Physics, Vol. 554. K. R. Mecke and D. Stoyan, editors. Springer-Verlag, Berlin.
41. Häggström, O., and K. Nelander. 1998. Exact sampling from anti-monotone systems. Stat. Neerl. 52:360–380.[CrossRef]
42. Gadgil, C., C.-H. Lee, and H. G. Othmer. 2005. A stochastic analysis of first-order reaction networks. Bull. Math. Biol. 67:901–946.[CrossRef][Medline]
43. Shreider, Y. A. 1966. The Monte Carlo Method. Pergamon Press, New York.
TOP
ABSTRACT
INTRODUCTION
[–1, 0] is the same as above. In this case, coalescence occurs at Tc = –0.17 and the value U*0 = L*0 is guaranteed to be a sample from the stationary solution of the Hill CME (
), 2(
), 3(
), and 15(
). Note that the error of the SSA runs converges to the asymptotic error floors for the different stopping times given by 
and 
of finding the system in each basin computed with the DCFTP-SSA compare very well with those obtained with the approximate eigenvector method (in parenthesis). Similarly, the expected time to reach the separators averaged over all the states in each basin (37
-distribution while the dotted line shows a Poisson distribution with the same mean. (d) Change of the period of the repressilator as a function of parameter variation. There is good agreement between the sensitivity of the stochastic system with respect to a ±5% change in kM (*), kP (