Bioinformatics Advance Access originally published online on August 19, 2004
Bioinformatics 2005 21(2):208-217; doi:10.1093/bioinformatics/bth479
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Bioinformatics vol. 21 issue 2 © Oxford University Press 2005; all rights reserved.
A method for estimating stochastic noise in large genetic regulatory networks
Institute for Systems Biology 1441 North 34th Street, Seattle, WA 98103, USA
*To whom correspondence should be addressed.
 |
Abstract |
|---|
 TOP Abstract 1 INTRODUCTION2 ALGORITHM3 IMPLEMENTATION4 DISCUSSIONREFERENCES |
|---|
Motivation: Genetic regulatory networks are often affected by stochastic noise, due to the low number of molecules taking part in certain reactions. The networks can be simulated using stochastic techniques that model each reaction as a stochastic event. As models become increasingly large and sophisticated, however, the solution time can become excessive; particularly if one wishes to determine the effect on noise of changes to a series of parameters, or the model structure. Methods are therefore required to rapidly estimate stochastic noise.
Results: This paper presents an algorithm, based on error growth techniques from non-linear dynamics, to rapidly estimate the noise characteristics of genetic networks of arbitrary size. The method can also be used to determine analytical solutions for simple sub-systems. It is demonstrated on a number of cases, including a prototype model of the galactose regulatory pathway in yeast.
Availability: A software tool which incorporates the algorithm is available for use as part of the stochastic simulation package Dizzy. It is available for download at http://labs.systemsbiology.net/bolouri/software/Dizzy/
Contact: dorrell{at}systemsbiology.org
Supplementary information: A conceptual model of the regulatory part of the galactose utilization pathway in yeast, used as an example in the paper, is available at http://labs.systemsbiology.net/bolouri/models/galconcept.dizzy
|
1 INTRODUCTION |
|---|
|
TOPAbstract1 INTRODUCTION2 ALGORITHM3 IMPLEMENTATION4 DISCUSSIONREFERENCES |
|---|
Genetic regulatory networks can be modeled using a number of methods, including ordinary differential equations (Bower and Bolouri, 2001; Neves and Iyengar, 2002), and Monte Carlo methods such as the Gillespie algorithm (Gillespie, 1976). An advantage of the former is that they are fast and easy to compute. The latter techniques explicitly account for stochastic effects due to the finite number of molecules per cell (McAdams and Arkin, 1999; Thattai and van Oudenaarden, 2001), but computation times can be slow. This is particularly the case for models of large genetic networks, which attempt to simulate many interactions among various species.
Two types of questions that typically arise with stochastic simulations of large networks are how to derive simple estimates for the effect on noise of a particular reaction, in terms of the reaction parameters; and how to estimate the noise in the entire system. Suppose for example that, as part of a larger genetic network, some protein X regulates the production of a protein Y. We might wish to explore, in isolation, the effect on stochastic noise of dimerization of X; or the effect of multiple binding sites for X on the promoter region for Y. An analytical approach is often preferable to simply performing many stochastic simulations, since it can provide insight into generic properties and usually requires less computation.
While viewing such sub-systems in isolation certainly makes calculations more tractable, changes in one component of a complex system often have a non-linear and sometimes counter-intuitive effect on other parts of the network. One therefore requires the capability to test the system as a whole, perhaps while making changes to a series of parameters or model structure. If the model is highly detailed, such an analysis can be prohibitively slow, even with the use of refined algorithms that improve the speed of computation (Gibson and Bruck, 2000; Gillespie, 2001; Gillespie and Petzold, 2003).
In this paper, we combine the ODE and stochastic methods, and derive a technique for estimating the noise levels in genetic networks of arbitrary size. A requirement of the method is that the ODE have a steady state, so it cannot be used for example to calculate the noise in a transient solution, or in an oscillator system with a limit cycle attractor. However, many models of biological systems, including ones for which experimental noise data exists (Paulsson, 2004), do satisfy this requirement, and there is a clear need for a practical means to quickly assess the noise characteristics of such networks. The algorithm can be used to both derive analytical estimates for simple cases, and numerical estimates for complete systems.
The method is compared with other techniques for simple cases, and used to analyze the effect of multiple binding sites in the regulation of transcription, and dimerization of proteins. It is then applied to a more complicated prototype model based on the galactose regulatory pathway in yeast, which includes features such as dimerization, double repression and negative auto-regulation. The effect of multiple binding sites in downstream genes is also discussed. A software tool that incorporates the method is available for use as part of the stochastic simulation program Dizzy. We begin in the next section by defining the method, and computing it analytically for some simple cases.
|
2 ALGORITHM |
|---|
|
TOPAbstract1 INTRODUCTION2 ALGORITHM3 IMPLEMENTATION4 DISCUSSIONREFERENCES |
|---|
2.1 Definitions and methods
Suppose that a biochemical network consists of N molecular species with concentrations x i molecules/cell, and M reactions. We first make the following definitions:
- r j is the N-dimensional unit column vector that describes the direction in ‘species space’ of the j-th reaction. For example, if the reaction is of the form
, then species x a and x b lose one molecule, and species x c gains one molecule. Therefore r j is the column vector with the value 
in rows a and b, and 
in row c.
- s j is the reaction rate. For the above example, the rate is s j = k j x a x b .
- J is the Jacobian of the ODE system equations. If the system is at a stable equilibrium, then the eigenvalues of J will have negative real part (Alligood et al., 1997).
- P is the (possibly complex) matrix formed from the column eigenvectors of J.
- u j are the associated eigenvalues.
- T is the diagonal matrix with entries
.
- Q is the matrix defined as Q = PTP –1.
- w j is the vector with entries which are the square of those in Qr j .
- s j is the reaction rate. For the above example, the rate is s j = k j x a x b .
The column vector giving the variance of each species is then given by
 |
Equation (1) can be interpreted as a weighted sum of the reaction rates. A natural application of the equation, which will be discussed in a future work, is to determine those reactions which contribute the most to noise in a particular species.
2.2 Derivation
The above result is a generalization of the error growth method presented by Orrell and Bolouri (2004). The linear propagator (Strang, 1986) is given by
 |
From the eigenvalue decomposition of the Jacobian, we can write
 |
where 
is the diagonal matrix consisting of the eigenvalues.
The propagated drift (Orrell, 2002) used to estimate the expected difference between the ODE and an ensemble of stochastic simulations, due to variance in reaction j with direction r j and rate s j , is given by
 |
where y j = P –1 r j s j and 
j is a Gaussian random variable with SD = 1 (Orrell and Bolouri, 2004).
The n-th entry of the propagated drift vector is
 |
Taking the limit as time goes to infinity, the square of the integral asymptotes to a variance
 |
Summing over the reaction index, and comparing with Equation (1), gives the desired result.
Estimates of the variance can also be obtained directly using techniques such as the 
-expansion (van Kampen, 1992) or a Langevin approach (Thattai and van Oudenaarden, 2002); however, in complicated models solving the equations for the variance is not straightforward. Another advantage of the error growth approach is that it incorporates the results into a general theory of error growth in non-linear systems. For simple cases, one would expect the results from our algorithm to be the same or similar to those obtained from other methods. In the next section, we benchmark the algorithm against other results for the simple case where one species regulates production of another, and investigate some implications for network behavior.
|
3 IMPLEMENTATION |
|---|
|
TOPAbstract1 INTRODUCTION2 ALGORITHM3 IMPLEMENTATION4 DISCUSSIONREFERENCES |
|---|
3.1 Noise due to system components
One use of the estimation technique is to analyze stochastic effects due to particular components of a larger genetic regulatory network. As an example, consider the following system consisting of two species x 1, x 2, and four reactions:
 |
where 
denotes the empty set.
The ODE equations are
 |
 |
At equilibrium, x 1 = v 1/d 1 and x 2 = [v 2 f(x 1)]/d 2. Letting the prime denote differentiation with respect to x 1, and using the equilibrium conditions, we then calculate:
 |
 |
 |
 |
Suppose now that f(x 1) = x 1 and d 1 >> d 2, as in the case where x 1 represents RNA and x 2 is a translated protein. The variance in x 2 is given by the second element v(2) of the vector v:
 |
where ß2 = v 2/d 1 gives the average number of proteins produced by one molecule of RNA before it decays. This agrees with the result for regulation of a single gene by Thattai and van Oudenaarden (2001).
A shorthand method to represent the process of transcription and translation is therefore to write it as two reactions in a single species:
 |
where x = x 2, d = d 2, v = v 1 v 2/d 1, and the stochasticity of the forward reaction is enhanced by a factor ß. Thus J = –d, 
, P = 1, 
, w 1 = w 2 = 1/d, and the variance is v = ßx + x = x(1 + ß) as before. One method to measure noise magnification is with the ratio of the variance to the mean [v(2)]/x 2, termed the Fano factor, which is unity for a purely Poissonian process. The term ß therefore represents the enhancement of internal stochasticity in the production process, above that which would exist in a Poissonian process.
Similarly, the case where x 1 represents one protein which regulates a second downstream protein x 2 becomes
 |
This case agrees with the corresponding result by Paulsson (2004), which was derived using van Kampen's -expansion (van Kampen, 1992). The -expansion technique involves finding the moments of a solution to a multivariate linear Fokker–Planck equation. In contrast, our algorithm is based on the error growth approach described by Orrell and Bolouri (2004). Although the two methods approach the problem very differently, they each require linearization of the dynamics, and should therefore agree at least for simple cases.
3.2 Multiple binding sites
A typical feature of genetic networks is the existence of multiple binding sites for regulatory proteins. It is interesting to ask whether the number of binding sites has any effect on the noise characteristics of the system. We here apply the above results to two cases: one where a protein x 1 represses the formation of a downstream protein x 2, and one where it promotes. For simplicity, we assume that the two proteins decay at the same rate, and are produced in a Poissonian fashion.
Suppose first that transcription occurs whenever any one of the binding sites is not bound with a protein, so x 1 represses the formation of x 2. The fractional saturation function, defined as the probability of the gene being in a promoted state, is then
 |
where K is the binding dissociation constant and N is the number of binding sites. Calculating the noise in the downstream protein from Equation (11) then gives
 |
For N = 1 binding site, the Fano factor achieves a maximum value of 2v 2/27dK at 
, for which the saturation function has f(x 1) = 
. Plotting the above function shows that the peak variance decreases with the number of binding sites.
Now suppose that transcription occurs whenever any one of the binding sites is bound with a protein, so x 1 promotes the formation of x 2. The saturation function is
 |
and
 |
The maximum Fano factor is v 2 N/2dK at x 1 = 0, with f(x 1) = 0. Therefore for positive regulation, the maximum increase in noise scales directly with the number of binding sites.
Figure 1 compares the positive and negative regulation scenarios for the case with N = 5 binding sites. Estimates using the above formulae agree well with detailed stochastic simulations from the program Dizzy. For the choices of parameters used to generate the figure, the case with positive regulation introduces somewhat more noise at low levels of induction (i.e. values of x 1 to the left of the square symbol in the right panels). This is intuitively reasonable, since the regulating protein is itself more affected by noise when present in low concentrations. The effect is reduced if there is a small degree of basal transcription.
|
3.3 Dimerization
Another important component of many genetic networks is dimerization of regulatory proteins. Suppose that a species can exist either as a monomer, with concentration in molecules per cell of x 1, or a homo-dimer, with concentration x 2. Assuming first-order creation and decay for the protein, the reactions and corresponding rates are:
 |
The ODE equations are
 |
 |
The predicted variance can then be estimated using the technique of Section 2. Assuming that dimerization occurs at a much faster rate than creation or decay of protein, so d << k 2, gives a variance around the ODE steady state of
 |
For any given choice of parameters, the variance in the dimer x 2 is less than the mean. To see this, let 
= k 1/k 2. Then the variance in x 2 is:
 |
where we used 
, from Equation (18b) and the assumption d << k 2. This expression has a maximum of about 0.68x 2 at x 1 = 
, and asymptotes to a minimum of 0.5x 2 in the limit as x 1 goes to infinity. Therefore, dimerization reduces noise to below that which would occur in a purely Poissonian process, for which the variance is equal to the mean x 2. It follows that dimerization can play an important role in reducing noise in genetic networks (Bundschuh et al., 2003).
Figure 2 compares the estimates with numerical results from Dizzy. The reaction rates are v = 10, d = 0.1, k 1 = 0.1, with k 2 ranging from 1 to 100. Units of concentration are molecules per cell, time units are arbitrary. The method tends to slightly underestimate the true noise. This is probably due to the strong nonlinearity caused by the second-order dependence of the dimerization rate on x 1. However, the method captures the dependence of noise on the rate k 2. Therefore even in cases where the estimation technique is not completely accurate, it may still be possible to obtain valuable results about the effect of changing model parameters or boundary conditions.
|
3.4 Noise in complex systems
While the estimation technique can be used to derive analytical expressions for noise in simple systems or sub-systems, as soon as the model becomes even moderately complex it is necessary to calculate the noise numerically. This section analyzes the noise properties of a prototype model based on the regulatory portion of the galactose pathway in yeast (Johnston, 1987; Lohr et al., 1995; Peng and Hopper, 2000, 2002). A more detailed version of this model, including the basis for the choice of parameter values, is presented by de Atauri et al. (2004); our aim here is to exploit it as an example of a complex genetic network with a number of interesting, non-linear features. The biological function of the network is to sense the presence of external galactose, and turn on various structural genes that control the galactose metabolic pathway. The protein gal4p promotes transcription by binding with the specific upstream activation sequence of the GAL genes, while the GAL80 protein gal80p represses transcription by binding with gal4p. In the absence of galactose, the concentration of gal80p is sufficient for the transcription of structural genes to occur only at a low basal level. When the GAL3 protein gal3p is activated by internal galactose, it restores transcription by relieving the repressive action of gal80p. Thus the presence of internal galactose has the effect of switching on the system by repressing the repressor.
This double repression network is further complicated by dimerization of the proteins gal4p, gal80p and gal7p; negative auto-regulation of gal80p; cooperative binding of the gal80p dimer to the DNA complex; and multiple binding sites (two) for the downstream gene gal7p. (As discussed below, other downstream genes have up to five binding sites.) The model therefore reflects a degree of the complexity that is found in typical genetic networks.
The model was first simulated using Dizzy, running the Gillespie algorithm (Gillespie, 1976). The steady-state results are summarized in Tables 1 and 2 for two cases: no induction, where internal galactose is not present, and low induction, where the system is partially activated. Previous simulations showed that noise is not a concern at high levels of induction, since the large concentrations of structural proteins mean that stochastic effects become negligible. The statistics were obtained over a long simulation time of 10 000 min (for the purposes of this simulation, the changing volume of the yeast cell, as it grows and divides, was not considered).
|
|
The stochastic noise was then estimated using the new tool in Dizzy, based on the methods of this paper. The steady state was first obtained by running the program's ODE solver until it reached equilibrium. The noise estimate was then computed automatically, using the option ‘–computeFluctuations’ in batch mode. Figure 3 compares the dimensionless ratio of SD with the mean from the stochastic simulations (open circles) to the estimates (+symbols). The two are in reasonable agreement, though the method tends to slightly underestimate the true noise for the species that experience homo-dimerization reactions. The method does a good job of capturing the differences between the species and between the two different states of induction.
|
A full model would include other structural proteins, such as gal1p, gal2p and gal10p, each of which have multiple binding sites for gal4p. For example, gal2p is thought to have five binding sites. The effect of the number of binding sites can be estimated simply by using gal7p as a generic protein downstream of the main regulatory apparatus. For example, if the number of binding sites in the GAL7 gene is increased to N = 5, then the estimate of the ratio of SD to mean in the gal7p dimer with low induction decreases from 0.1783 to 0.1204. However, this improvement is primarily due to an increase in the mean from 329.85 to 723.44. Therefore, the functional importance of multiple binding sites is likely to lie more with the dynamic properties, such as basal transcription, than with the noise characteristics. This would be consistent with the results of Section 3, where it was seen that noise magnification is relatively low in systems operating through repression.
The application of our noise estimation method thus permits the identification of some properties that help to control noise. These include dimerization of some of the key proteins; the use of negative regulation; and also negative auto-regulation in the gal80p protein, which can be shown to reduce noise (Thattai and van Oudenaarden, 2001).
|
4 DISCUSSION |
|---|
|
TOPAbstract1 INTRODUCTION2 ALGORITHM3 IMPLEMENTATION4 DISCUSSIONREFERENCES |
|---|
The above examples show that the estimation technique can be used both analytically to determine the noise levels for simple systems, and also to estimate stochastic effects in large networks. It is capable of dealing with non-linear features such as dimerization, feedback loops and cooperative binding. Other effects can also be handled so long as they do not affect the basic model structure. One example is compartmentalization, where reactions are assumed to take place in distinct compartments such as the nucleus and cytoplasm. In Dizzy, each species is assigned a different name depending on which compartment it is located in, and diffusion is modeled by a first-order reaction. Therefore, the noise due to this process can be estimated in the same way as any other reaction. Similarly the effect of dimension reduction, when reactions occur on membranes, does not alter the basic model structure and can be handled using this method.
Of course, the technique can fail whenever the linearization of the dynamics is not a good approximation to the true dynamics. This can occur for example near bifurcation points (Guckenheimer and Holmes, 1983; Borisuk and Tyson, 1998) of the ODE. If the system ODE has multiple steady states, as in a bistable switch (Isaacs et al., 2003; Ozbudak et al., 2004), then the presence of noise may result in basin jumping (Alligood et al., 1997) from one state to another. Also, the technique can only estimate the magnitude of the variance, and does not produce time-series plots. For these reasons, final results should always be confirmed by a detailed stochastic simulation.
One situation where the method can give errors is in systems which feature zeroth order ultrasensitivity (Berg et al., 2000). Consider a system in which a target protein is switched between two states x and y, so the total concentration x + y = T is conserved. The concentration y can always be deduced from y = T – x, so the system can be described in terms of x alone. The reactions are:
 |
The first reaction corresponds to the modification of x to the state y, while the second reaction corresponds to demodification, so the sink can be identified with y. The corresponding ODE equation is
 |
This equation has the property, for small values of K, that the equilibrium concentration depends sensitively on whether k 1 is smaller or larger than k 2, with a sudden switch in the equilibrium value of x between nearly 0 and nearly T. As discussed by Berg et al. (2000), a stochastic simulation shows a more moderate transition. However, if we apply the estimation technique to the case k 1 = k 2, then the Jacobian at the equilibrium is
 |
For small K the Jacobian approaches zero, with a minimum magnitude at 
of about 8k 1 K/T 2, so the estimated variance in Equation (6) tends to blow up. For example, with T = 100 molecules, the estimated SD at k 1 = k 2 is 35.7, while the actual value from a stochastic simulation is 22.3. Errors are smaller at other operating points, and also depend on different factors such as the number of molecules T.
It would clearly be desirable to detect when, and by how much, the method is likely to fail. One possibility is to launch an ensemble of short trajectories from initial conditions that are perturbed according to the expected fluctuation magnitude. If the perturbed trajectories do not decay towards the steady state in a manner consistent with the estimated decay rates, or do not decay at all, then the linearization is not valid. An example would be a system with multiple steady states, for which perturbed trajectories might be attracted towards a different equilibrium. A topic of future work is to incorporate such a test as an automatic feature in Dizzy.
The main advantage of fast estimation techniques is that they allow rapid assessment of effects that would be extremely time consuming to perform with a full simulation. For example, a full stochastic simulation of the galactose pathway can take hours to compute on a desktop machine, which makes it difficult to explore the effect of altering structure or parameter values on noise properties. In contrast, estimation of the noise level takes a matter of seconds. It can therefore be used to quickly assess the likely effect of a change to model parameters or model structure; or to indicate when stochastic noise is expected to make ODE simulations inaccurate (Vasudeva and Bhalla, 2004).
Stochastic effects due to low numbers of molecules are known to play an important role in genetic networks. Organisms have presumably therefore evolved to both control and exploit the effects of noise. As models of such systems grow increasingly complex and sophisticated, and incorporate non-linear effects such as feedback loops, the consequences of noise become increasingly difficult to predict or comprehend. Estimation techniques such as that presented here can help both to understand the effects of noise in relatively simple systems, and give guidance to the likely effects in large networks.
Received on March 3, 2004; revised on July 16, 2004; accepted on August 11, 2004
|
REFERENCES |
|---|
|
TOPAbstract1 INTRODUCTION2 ALGORITHM3 IMPLEMENTATION4 DISCUSSIONREFERENCES |
|---|
Alligood, K.T., Sauer, T.D., Yorke, J.A. Chaos: An Introduction to Dynamical Systems, (1997) , New York Springer-Verlag.
Berg, O.G., Paulsson, J., Ehrenperg, M. (2000) Fluctuations and quality of control in biological cells: zero-order ultrasensitivity reinvestigated. Biophys. J., 79, , pp. 1228–1236.
Borisuk, M.T. and Tyson, J.J. (1998) Bifurcation analysis of a model of mitotic control in frog eggs. J. Theor. Biol., 195, 69–85[CrossRef][ISI][Medline].
(Eds.). Computational Modeling of Genetic and Biochemical Networks, (2001) , Cambridge, MA MIT Press.
Bundschuh, R., Hayot, F., Jayaprakash, C. (2003) The role of dimerization in noise reduction of simple genetic networks. J. Theor. Biol., 220, , pp. 261–269.
de Atauri, P. and Orrell, D., et al. (2004) Evolution of ‘design’ principles in biochemical networks. IEE Syst. Biol., 1, 28–40[CrossRef].
Gibson, M.A. and Bruck, J. (2000) Efficient exact stochastic simulation of chemical species with many species and many channels. J. Phys. Chem., 104, 1876–1889.
Gillespie, D.T. (1976) A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys., 22, 435–450[CrossRef].
Gillespie, D.T. (2001) Approximate accelerated stochastic simulation of chemically reacting systems. J. Chem. Phys., 115, 1716–1733[CrossRef].
Gillespie, D.T. and Petzold, L.R. (2003) Improved leap-size selection for accelerated stochastic simulation. J. Chem. Phys., 119, 8229–8234[CrossRef].
Guckenheimer, J. and Holmes, P. Nonlinear Oscillations, Dynamic Systems, and Bifurcations of Vector Fields, (1983) , New York Springer-Verlag.
Isaacs, F.J., Hasty, J., Cantor, C.R., Collins, J.J. (2003) Prediction and measurement of an autoregulatory genetic module. Proc. Natl Acad. Sci. USA, 100, , pp. 7714–7719
Johnston, M. (1987) A model fungal gene regulatory mechanism: the GAL genes of Saccharomyces cerevisiae . Microbiol. Rev., 51, 458–476
Lohr, D., Venkov, P., Zlatanova, J. (1995) Transcriptional regulation in the yeast GAL gene family: a complex genetic network. FASEB J., 9, 777–787[Abstract].
McAdams, H.H. and Arkin, A. (1999) It's a noisy business! Genetic regulation at the nanomolar scale. Trends Genet., 15, 65–69[CrossRef][ISI][Medline].
Neves, S.R. and Iyengar, R. (2002) Modeling of signaling networks. BioEssays, 24, 1110–1117[CrossRef][ISI][Medline].
Orrell, D. (2002) Role of the metric in forecast error growth: how chaotic is the weather?. Tellus, 54A, 350–362.
Orrell, D. and Bolouri, H. (2004) Control of internal and external noise in genetic regulatory networks. J. Theor. Biol., 230, 301–312[CrossRef][ISI][Medline].
Ozbudak, E.M., Thattai, M., Lim, H.N., Shraiman, B.I., Van Oudenaarden, A. (2004) Multistability in the lactose utilization network of Escherichia coli . Nature, 427, 737–740[CrossRef][Medline].
Paulsson, J. (2004) Summing up the noise in genetic networks. Nature, 427, 415–418[CrossRef][Medline].
Peng, G. and Hopper, J. (2000) Evidence for Gal30's cytoplasmic location and Gal80p's dual cytoplasmic-nuclear location implicates new mechanisms for controlling Gal4p activity in Saccharomyces cerevisiae . Mol. Cell. Biol., 20, 5140–5148
Peng, G. and Hopper, J. (2002) Gene activation by interaction of an inhibitor with a cytoplasmic signaling protein. Proc. Natl Acad. Sci. USA, 99, 8548–8553
Strang, G. Introduction to Applied Mathematics, (1986) Wellesley-Cambridge Press.
Thattai, M. and van Oudenaarden, A. (2001) Intrinsic noise in gene regulatory networks. Proc. Natl Acad. Sci. USA, 98, , pp. 8614–8619
Thattai, M. and van Oudenaarden, A. (2002) Attenuation of noise in ultrasensitive signaling cascades. Biophys. J., 82, 2943–2950
van Kampen, N.G. Stochastic Processes in Physics and Chemistry, (1992) , North-Holland, Amsterdam .
Vasudeva, K. and Bhalla, U.S. (2004) Adaptive stochastic-deterministic chemical kinetic simulations. Bioinformatics, 20, , pp. 78–84
CiteULike 
Connotea 
Del.icio.us What's this?
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||