Amplitude distribution of stochastic oscillations in biochemical networks due to intrinsic noise
 Moritz Lang^{1},
 Steffen Waldherr^{1}Email author and
 Frank Allgöwer^{1}
Received: 6 March 2009
Accepted: 17 November 2009
Published: 17 November 2009
Abstract
Intrinsic noise is a common phenomenon in biochemical reaction networks and may affect the occurence and amplitude of sustained oscillations in the states of the network. To evaluate properties of such oscillations in the time domain, it is usually required to conduct longterm stochastic simulations, using for example the Gillespie algorithm. In this paper, we present a new method to compute the amplitude distribution of the oscillations without the need for longterm stochastic simulations. By the derivation of the method, we also gain insight into the structural features underlying the stochastic oscillations. The method is applicable to a wide class of nonlinear stochastic differential equations that exhibit stochastic oscillations. The application is exemplified for the MAPK cascade, a fundamental element of several biochemical signalling pathways. This example shows that the proposed method can accurately predict the amplitude distribution for the stochastic oscillations even when using further computational approximations.
PACS Codes: 87.10.Mn, 87.18.Tt, 87.18.Vf
MSC Codes: 92B05, 60G10, 65C30
1. Introduction
Oscillations are a widely occurring phenomenon in the dynamics of biological systems. On the intracellular level, oscillations occur for example in the activity of various genes or signalling proteins. To gain more insight into the processes related to these oscillations, mathematical models for the underlying biochemical and genetic networks are commonly constructed. Such models have proven helpful in connecting the underlying biochemistry with the temporal characteristics of emerging oscillations and the associated biological function. Examples for this type of results include the intracellular circadian clock [1] or the developmental process of somitogenesis [2].
In dynamical models for biological systems, oscillations are due to limit cycles or more complex attractors in deterministic frameworks, or may emerge from various stochastic effects in stochastic frameworks. Concerning the biological function, one would not expect that it makes a significant difference whether oscillations are due to a deterministic limit cycle or arise from stochastic effects. For the biological function, rather the temporal characteristics of oscillations are relevant, such as their frequency or amplitude. For deterministic limit cycle oscillations, these characteristics are easily computed by a numerical simulation of the model. For stochastic models, computing the temporal characteristics of oscillations is much more involved. Two approaches are common for solving stochastic systems. In the first approach, the so called chemical master equation (CME) is used [3]. The solution to the CME (or approximations thereof) yields the probabilities for each of the possible configurations of the system to be attained at a given point in time. However, since the temporal correlations of these probabilities are usually not obtained from the CME, the temporal characteristics of oscillations are not directly accessible from such a solution. In the second approach, a large number of realisations of the stochastic process are computed. These realisations can then be used to compute various temporal characteristics of the system, in particular oscillation amplitude and frequency distributions. Yet, the computation of a sufficiently large number of realisations typically entails a large computational effort.
In this paper, we develop a new method to compute the amplitude distribution for systems exhibiting stochastic oscillations. We thereby focus on systems where the deterministic part has a weakly stable equilibrium point (EP), typically with damped oscillations, and the stochastic effects induce sustained oscillations around this EP. In such systems, we can distinguish two mechanisms by which stochastic oscillations occur. For the first mechanism, the system needs to have the property that a certain small perturbation away from the equilibrium point leads to a large excursion in the state space before the system returns to the proximity of the EP [4]. If the noise is strong enough to reliably generate such a perturbation, but does not disturb the system during the round trip too much, we recognize regular oscillations with a well defined amplitude and frequency. This effect is called Coherence Resonance [5] or Stochastic Resonance [6, p. 149 ff.]. For the second mechanism, the deterministic part of the system already gives rise to damped oscillations, and the noise just serves to sustain the oscillations [7, 8]. In this case, the density distribution of realizations in the state space will typically be similar to a Gaussianlike distribution with a maximum at the EP, while individual trajectories show pronounced oscillations with a frequency close to the frequency of the deterministic part. The amplitudes of this type of oscillations increase with the noise power, in contrast to the other type of systems, which mainly exhibit oscillations of a fixed amplitude independent of the noise power.
For the second mechanism, it is of interest to compute the amplitude distribution of the stochastic oscillations. Several methods to solve this problems were developed previously. In [9] the average change of a stochastic Lyapunov function was determined. By setting this average change to zero one may identify an orbit around which a stochastic realization fluctuates, and thereby estimate a typical amplitude. In the context of this article the results of Kuske et al. [7, 8] are very interesting. They apply a multiscale ansatz to separate the properties of the oscillations mainly determined by the deterministic part of the analyzed system (the frequency) from those determined to a high percentage by the stochastic part (the amplitude). Therewith it is possible to calculate both the amplitude distribution and the mean frequency of the oscillations. Nevertheless this method has the disadvantage of requiring the discussed systems to be almost linear and the frequency of the oscillations not to be noteworthy disturbed by the stochastic effects.
In this article we focus on biochemical systems that can be modeled as a set of interconnected, possibly nonlinear, stochastic differential equations (SDEs). With the help of the FokkerPlanck equation (see e.g. [3, p. 193 ff.] or [10, p. 120 ff.]), we calculate the stationary density distribution in the state space. We provide a theorem allowing to calculate the amplitude distribution with only the knowledge of the density distribution and the corresponding SDE. The theorem allows not only to consider linear, but also a wide class of nonlinear systems and therefore makes it possible to analyze not only systems with oscillations being of low amplitude compared to average concentrations but also of intermediate and high amplitudes. Although stochastic oscillations are a mainly two dimensional effect, we show that they also occur in higher order systems and give an example on how to analyze them then. As far as we are aware of, this paper introduces for the first time a method to analytically analyse higher order possibly nonlinear systems being affected with intermediate to high amounts of noise, and thereby showing stochastic oscillations.
The paper is structured as follows. To give an easily accessible introduction to the material, we first present our results on the calculation of the amplitude distribution of stochastic oscillations and afterwards state the underlying assumptions and theorems. Then we explain the application of our theory by calculating the amplitude distribution of a simple example, the damped harmonic oscillator in the presence of additive noise, to get a first insight into the reasons for stochastic oscillations. We want to remark that this example simplifies, due to its structure, the necessary calculations significantly. As a more realistic example, we discuss oscillations in the MAP kinase cascade with an incorporated negative feedback with limited amounts of entities of each molecular species. For this system, the stationary probability distribution can be estimated by a linear approximation, or it can be computed numerically. We compare the amplitude distributions predicted by our method, based on these two approaches, to an "experimental" amplitude distribution obtained from a longterm stochastic simulation.
2. Results and Discussion
2.1. Amplitude Distribution of the Stochastic Oscillations
In this article we develop a method to derive the amplitude distribution of stochastic oscillations from the knowledge of the stationary density distribution of a stochastic differential equation (SDE). Here we first give a short outline of our results for which we present the corresponding theorems and proofs in the next sections.
with the state vector x ∈ ℛ^{2}, the system's dynamics f = (f_{1}, f_{2})^{ T }, where f_{1} and f_{2} are smooth functions (∈ C^{∞}), Σ a 2by2 matrix of smooth functions and Γ = (Γ_{1}, Γ_{2})^{ T }, where Γ_{1} and Γ_{2} are uncorrelated, statistical independent Gaussian white noise with zero mean and variance of one.
Hereby A is defined as A = {x ∈ ℛ^{2}⟨ x ⟩ = χ } and is the unit vector in the direction where the amplitude is measured.
The product on the right hand side of (2) represents the steady state flux of reactions through a certain state . Intuitively speaking, is chosen so that most of the realizations going through will also reach their maximal value in a small neighborhood of , so that under certain assumptions (see following sections) it is justifiable to approximate the amplitude of such a realization by the value reached at .
The apparent complicated definition of is a result of the freedom of choice in which direction the amplitude is measured. In a biological system it is obviously important to distinguish whether one measures the oscillations of the concentration of species A or of species B. Sometimes it is experimentally not possible to distinguish between two species, e.g. if B represents the phosphorylated version of a protein A, so that one can only measure the amplitude of [A] + [B]. Additionally, for nonlinear networks, the amplitude distribution may be different depending on whether one measures the amplitude of the oscillations in the positive direction from the steady state or in the negative direction. If one e.g. measures the amplitude distribution of the concentration of the first species x_{1} in the positive direction, (3), i.e. = (1, 0)^{T}, can be simplified to .
We want to emphasize here that (2) allows to analytically or numerically calculate the amplitude distribution only from the knowledge of certain properties of the SDE describing the biological system without the necessity to run numerical simulations. The formula yields a good approximation for a wide class of nonlinear problems and can be used to calculate the amplitude distribution even for certain systems having more than two dimensions as we demonstrate in the second example of this paper.
2.2. Derivation of the Results
To derive the results presented in the proceeding section we analyze nonlinear SDE systems of order two as given by (1), with an asymptotically stable equilibrium point and damped oscillations in the deterministic part. We restrict our analyis to two dimensional systems, because many higher dimensional systems can be reduced to two dimensions for the purpose of analysing stochastic oscillations. In the derivation we assume w.l.o.g. that we want to compute the amplitude of the oscillations in the positive direction of x_{1}, i.e. = (1, 0)^{T}, and that the equilibrium point of the deterministic formulation ( = f(x)) of (1) is at the origin. Both assumptions can easily be satisfied for an arbitrary system by an appropriate coordinate transformation.
We establish an angular phase relationship of the state vector x(t) with respect to a fixed reference vector in the state space [12]. An oscillation period of the system (1) is defined as the time during which the angle between the state vector and the fixed reference vector changes by 2π. Further, the amplitude of an oscillation is defined as the maximal value of x_{1} during the corresponding oscillation period.
We make the following assumptions on the system (1):
1. There exists a stationary density distribution P: ℛ^{2} → ℛ^{+}: x ↦ P (x) for (1) which is sufficiently smooth in x. We demand the curvature of the level curves of this distribution to exist and not to change its sign. For the computation of the amplitude distribution P_{ A }(χ), P (x) may be computed analytically, numerically, or may be obtained by the long term limit of a measurement.
which has to be small for almost all x ∈ ^{2} having a nonnegligible probability P (x).
3. In addition we require that the average speed ν =  f(x) of (1) does not vanish for almost all x ∈ ℛ^{2} having a nonnegligible probability P (x).
Assumption 3 is obviously not satisfied in a small area around the EP where for most systems ν becomes very small. However when the system is very close to the EP, the motion of most realizations can be approximated by a random walk. Thus during the time when the system is in the small neighborhood of the EP it is not justifiable to speak about oscillations anymore until the system leaves the neighborhood again. Away from the EP,  f(x) is significantly bigger than zero for most biological relevant systems so that Assumption 3 is fulfilled.
With the help of these definitions we can formulate the following two lemmas and a theorem, which characterizes the amplitude of the stochastic oscillations.
Lemma 1 Assume the Assumptions 13 are satisfied. Let Ψ_{ f }be a realization of (1) being at the state ∈ at time t_{ 0 }. Then the amplitude of Ψ_{ f }during the current oscillation will lie with a probability of 70.8% in the set [ , + δx_{1}] and with a probability of 95.5% in the set [ , + ], with δx_{ 1 } defined by .
Lemma 2 Assume the Assumptions 13 are satisfied. Then the net flux density (x) of realizations at the state x ∈ ℛ^{2} has an absolute value proportional to the product of the probability P (x) and the average speed v(x) of (1) and is tangential to the level curve of the density distribution P at this point.
Theorem 1 (Distribution of amplitudes) Under the Assumptions 13 and if δx_{ 1 } as defined in Lemma 1 is small, the probability P_{ A }(χ) of an amplitude χ of the stochastic oscillations of (1) in the direction x_{1} is proportional to the absolute value of the net flux density   at the state ∈ satisfying = χ.
Our results for the calculation of the amplitude distribution P_{ A }(χ) presented in the previous section (2) directly follow from the application of Theorem 1. The proofs of the respective lemmas and theorem can be found in the appendix.
2.3. Remarks

To check if the assumptions necessary for the application of Theorem 1 are fulfilled, one can calculate δx_{1} according to the formula in Lemma 1. If δx_{1} is significantly smaller than for almost all ∈ having a nonnegligible probability P ( ), Theorem 1 may be applied. For a general stochastic differential equation, it is hard to exactly indicate how much smaller δx_{1} must be compared to . This highly depends on the structure of the system and also on the accuracy being required for the solution. Nevertheless, as a rule of thumb, the theorem gives good approximations if < 0.2. It is not necessary that this is true for every ∈ , but the calculated density of the amplitude may differ from the real one for the amplitudes χ ∈ for which the corresponding δx_{1} is not small enough. Nevertheless the accuracy of the amplitude distribution for other amplitudes is not affected (see proof).

It may be possible to relax the assumption that the curvature of the level curves does not change its sign globally. Intuitively it seems to be sufficient for most systems that the curvature does not change its sign only locally around the states ∈ , and that the amplitude in one oscillation period reaches its maximum in the area around . However this has to be shown for the system of interest explicitly.

If the Lie derivative L_{ f }P of the density distribution P is not small enough, the net flux density cannot be approximated to be tangential to the level curves of P anymore. Depending on the size of L_{ f }P this may lead to a significant bias in the calculation of the amplitude distribution for certain systems. In this case, the derivation of the amplitude distribution would be more involved, and also depend on the term , starting from an appropriate modification of (66) in the appendix.
2.4. The Damped Harmonic Oscillator in the Presence of Additive Gaussian White Noise
2.4.1. Equations and Properties of the System
with Γ_{ i }, i ∈ {1, 2} uncorrelated, statistical independent Gaussian white noise with zero mean and variance of one, σ ∈ ^{+} and k ∈ .
are at λ_{1,2} = k ± i, therefore the deterministic system will show damped oscillations upon perturbation from its steady state.
2.4.2. Calculation of the Density Distribution
with ⟨ x ⟩ the mean value of x and Ξ the matrix of the second moments of P, which is the solution of the equation AΞ + ΞA^{ T }+ B = 0. In system (9), ⟨ x ⟩ = 0 and Ξ may be calculated to . The state of the system is therefore Gaussian distributed with the maximum of the density being at the EP.
2.4.3. Determination of the Amplitude Distribution
For s = 10 and k = 0.01, δx_{1} is getting small compared to for ≥ 20, so our approximation should at least hold for every amplitude greater or equal to 20. The results will nevertheless show that we even get good estimations for amplitudes much smaller than 20.
This easy example was discussed to give an insight into the reasons for stochastic oscillations. In the next section, we show the practical applicability of our algorithm by predicting the amplitude distribution of a complex biochemical system and therewith give an example of a more biologically relevant application for the results of the paper.
2.5. Oscillations in the MAP Kinase Signaling Cascade
2.5.1. The Deterministic Model
For this example, we use a basic ODE model of the MAP kinase signaling cascade. The model contains a negative feedback interconnection from the last kinase MAPK to the activation of the first kinase MAPKKK [15]. We use a simpler model compared to the one in [15], and therefore sustained oscillations do not occur in the deterministic version of the model considered here.
Parameter Set for the MAP Kinase Signaling Cascade
Parameter  Value  Unit 

k _{11}  0.256  1/s 
k _{12}  0.0872  1/nM 
k _{13}  0.1144  1/nM 
k _{21}  0.4229·10^{2}  1/s 
k _{22}  0.0347  1/nM 
k _{31}  0.2049·10^{2}  1/s 
k _{32}  0.0243  1/nM 
p _{11}  0.1003  1/s 
p _{12}  0.1428  1/nM 
p _{21}  0.1218  1/s 
p _{22}  0.0794  1/nM 
p _{31}  0.1343  1/s 
p _{32}  0.0167  1/nM 
and thus the deterministic system is locally asymptotically stable.
2.5.2. The Stochastic Model
2.5.3. Transformation and Model Reduction
with z_{1} ∈ the coordinate of the fast and z_{ slow } = (z_{2}, z_{3})^{ T } the coordinates of the slow manifold, ψ_{ fast }and ^{ψ}_{ slow }= (ψ_{slow, 1}, ψ_{slow, 2})^{ T }vectors of polynomials of order two and higher, ϕ_{ fast }and ϕ_{ slow }matrices of the noise strength and Γ = (Γ_{1}, Γ_{2}, Γ_{3})^{ T }the noise vector of system (25).
Parameter Set of the Slow Manifold of the MAP Kinase Signaling Cascade
Parameter  Value  Parameter  Value  Parameter  Value 

a _{0,0}  0  a _{0,2}  0.01800  a _{4,0}  7.327·10^{8} 
a _{1,0}  0  a _{3,0}  8.476·10^{6}  a _{3,1}  8.191·10^{7} 
a _{0,1}  0  a _{2,1}  3.625·10^{5}  a _{2,2}  8.162·10^{7} 
a _{2,0}  0.003604  a _{1,2}  2.425·10^{4}  a _{1,3}  3.951·10^{6} 
a _{1,1}  0.01102  a _{0,3}  5.750·10^{4}  a _{0,4}  1.785·10^{5} 
to which Theorem 1 can be applied.
2.5.4. Calculation of the Density Distribution
As a first approach to obtain the stationary density distribution P (x), we make use of a linear approximation of the reduced system (32) around the EP. Taking this approach allows us to evaluate how well our method works with such an approximation, where the stationary density distribution can be obtained with minimal computational effort. For other systems, or if a high precision of the result is required, it could however also be necessary to solve the nonlinear FokkerPlanck equations numerically.
and Γ the vector of the disturbances Γ_{ i }, i ∈ {1, 2, 3}, of the original system (25).
2.5.5. Determination of the Amplitude Distribution
under the constraint that z_{ χ }has to lie on the slow manifold. This corresponds to that z_{ χ }has to satisfy (31), which gives us a dependency of β on α. This dependency can be obtained numerically by solving a convex optimization problem.
As can be seen in Figure 9, the simulation results fit the calculated predictions very well. However, it seems that we tend to underestimate the amplitude of the oscillations by a small amount. This can be explained as a result of Lemma 1, which states that the amplitude of a realization going through the state lies with a high probability in [ , + δx_{1}], whereas we estimate its value by the lower end of this interval (see appendix). The tendency to underestimate the amplitude seen in Figure 9 seems to be a direct result of neglecting the value of δx_{1}. Further discussion of this point can be found in the conclusions.
We also want to mention that we compare the results obtained from our method to stochastic Gillespie simulations, and not to realizations of the Langevin equation (24). The reason herefor is that Gillespie simulations better predict the behavior of biochemical networks and are thus the method of choice. However in the derivation of the amplitude distribution we approximated the stochasticity of the system using white noise terms. The results of both methods to describe the intrinsic noise of the system may under certain circumstances lead to different results, which may be a further explanation of the small bias between the measured and the theoretical predicted amplitude distribution in Figure 9. Furthermore some bias may be explained due to the linearization of the system around its steady state.
2.5.6. Numerical Approach
Sometimes it is not adequate anymore to analyze the linearized system and one has to analyze the original nonlinear one. Although there might be several special cases where the probability distribution P of a system of nonlinear stochastic differential equations is analytically computable, this is not possible in the general case. As a consequence there is only the possibility to obtain the solutions numerically. To give an example how to solve a problem with the numerical approach, we decided to analyze the same system as in the preceding section, as defined in (24), except of changing the cell volume to V = 0.017 pl, which corresponds to of the original cell volume. Because the variance of the disturbances is highly dependent on the amount of entities of the different protein species, the shape of the solution changes dramatically in a way that we cannot get good results anymore by analyzing the linearized system.
3. Conclusion
We introduce a method to determine the amplitude distribution of a wide class of linear and nonlinear stochastic systems given by (1), which display sustained stochastic oscillations. The method is applicable to systems where a stationary density distributions exists and can be computed either analytically or numerically.
The method is based on computing the flux density of realizations for the states where the tangent on the level curve of the density distribution is normal to the direction in which the oscillations are measured. We showed that under certain conditions this flux density is directly proportional to the probability of an oscillation with an adequate amplitude to occur. Our results can be used in the analysis of systems being influenced by strong internal or external noise, as we often find them in biophysical problems.
As already discussed at the end of the second example, for certain systems the calculated amplitude distribution may contain a small bias depending on the exact structure of the system and on how well the the assumption necessary for the application of our method are satisfied. For the wide class of nonlinear systems analyzed we could obtain as a result that realizations being at state will most likely have an amplitude lying in [ , + δx_{1}] as stated in Lemma 1. By requiring δx_{1} to be small we justified to approximate the amplitude by . However, if for specific systems the distribution in [ , + δx_{1}] can be calculated or estimated, it would be possible to reduce the bias for systems where δx_{1} is small but not negligible. Such an extension would lead to a more refined approximation of the amplitude distribution.
Stochastic oscillations may occur in biological systems not only as a disturbing side effect, but also in a constructive manner [21, 22]. This is not because stochastic oscillations have great benefits over more "traditional" types of oscillations like deterministic limit cycles, but due to the fact that there seems to be no reason for a preference of deterministic limit cycles. In biophysical systems, the type of oscillations we study in this paper often occurs at parameter values in the vicinity of a Hopf bifurcation in the deterministic part of the model [23]. This is important in the field of robustness analysis of biological networks [24], because stochastic oscillations can possibly improve the robustness of oscillations in a network. In this respect, our method to compute the oscillation amplitude may be helpful in order to characterize a parameter region in which either deterministic or stochastic oscillations of a comparable amplitude occur robustly.
Appendix
Proof for Lemma 1
In this section we give a short proof for Lemma 1. We therefore consider the level curve defined by of the density distribution going through the state ∈ .
After substituting the definitions made in the equations (8) and utilizing the knowledge that vanishes because of the definition of , setting
In the same way it can be shown that 95.5% of all realizations have an amplitude lower than + .
Proof for Lemma 2
Proof for Theorem 1
For small δx_{1} the integral of Ω over is approximately 1. The result of this last approximation is given by Theorem 1.
Declarations
Authors’ Affiliations
References
 Leloup JC, Goldbeter A: Proc Natl Acad Sci USA. 2003, 100 (12): 70517056. 10.1073/pnas.1132112100.View ArticleADSGoogle Scholar
 Tiedemann HB, Schneltzer E, Zeiser S, RubioAliaga I, Wurst W, Beckers J, Przemeck GKH, de Angelis MH: J Theor Biol. 2007, 248: 120129. 10.1016/j.jtbi.2007.05.014.View ArticleGoogle Scholar
 van Kampen NG: Stochastic processes in Physics and Chemistry. 2007, Amsterdam, The Netherlands: Elsevier, 3Google Scholar
 Muratov CB, VandenEijnden E: Chaos. 2008, 18: 015111101511111. 10.1063/1.2779852.View ArticleADSMathSciNetGoogle Scholar
 ElSamad H, Khammash M: Proceedings of the 45th IEEE Conference on Decision and Control. 2006, 23822387. 10.1109/CDC.2006.377181.View ArticleGoogle Scholar
 Beuter A, Glass L, Mackey MC, Titcombe MS: Nonlinear Dynamics in Physiology and Medicine, of Interdisciplinary Applied Mathematics. 2003, New York: Springer, 25:View ArticleGoogle Scholar
 Kuske R, Gordillo LF, Greenwood P: Journal of Theoretical Biology. 2007, 245: 459469. 10.1016/j.jtbi.2006.10.029.View ArticleMathSciNetGoogle Scholar
 Yu N, Kuske R, Li Y: Physical Review E. 2006, 73: 05620510562058.ADSGoogle Scholar
 Aparicio JP, Solari HG: Mathematical Biosciences. 2001, 169: 1525. 10.1016/S00255564(00)00050X.View ArticleMathSciNetMATHGoogle Scholar
 Gillespie DT: Markov Processes. 1992, San Diego, CA: Academic Press, INCMATHGoogle Scholar
 Gillespie DT: Journal of Chemical Physics. 2000, 113: 297306. 10.1063/1.481811.View ArticleADSGoogle Scholar
 Bagheri N, Stelling J, Doyle FJ: Bioinformatics. 2007, 23 (3): 358364. 10.1093/bioinformatics/btl627.View ArticleGoogle Scholar
 Gardiner CW: Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences. 2004, Springer, 3View ArticleMATHGoogle Scholar
 Pearson G, Robinson F, Gibson TB, Xu BE, Karandikar M, Berman K, Cobb MH: Endocrine Reviews. 2001, 22 (2): 153183. 10.1210/er.22.2.153.Google Scholar
 Kholodenko B: European Journal of Biochemistry. 2000, 67 (6): 15831588. 10.1046/j.14321327.2000.01197.x.View ArticleGoogle Scholar
 Hilioti Z, Sabbagh W, Paliwal S, Bergmann A, Goncalves MD, Bardwell L, Levchenko A: Curr Biol. 2008, 18 (21): 170006. 10.1016/j.cub.2008.09.027.View ArticleGoogle Scholar
 Ramsey S, Orrell D, Bolouri H: Journal of Bioinformatics and Computational Biology. 2005, 3: 415436. 10.1142/S0219720005001132.View ArticleGoogle Scholar
 Khalil HK: Nonlinear Systems. 2001, Upper Saddle River, NJ, USA: PrenticeHall, 3Google Scholar
 Jacobs SJ: Journal of Atmospheric Sciences. 1991, 48 (7): 893901. 10.1175/15200469(1991)048<0893:EOASMI>2.0.CO;2.View ArticleADSGoogle Scholar
 Sjöberg P: Numerical Solution of the FokkerPlanck Approximation of the Chemical Master Equation. PhD thesis. 2005, Uppsala UniversitetGoogle Scholar
 Rao CV, Wolf DM, Arkin AP: Nature. 2002, 420 (6912): 231237. 10.1038/nature01258.View ArticleADSGoogle Scholar
 Steuer R, Zhou C, Kurths J: Biosystems. 2003, 72 (3): 241251. 10.1016/j.biosystems.2003.07.001.View ArticleGoogle Scholar
 Waldherr S, Allgöwer F: International Journal of Systems Science. 2009, 40: 769782. 10.1080/00207720902957269.View ArticleMATHGoogle Scholar
 Kitano H: Nat Rev Genet. 2004, 5 (11): 826837. 10.1038/nrg1471.View ArticleGoogle Scholar
 Bronstein IN, Semendjajew KA: Taschenbuch der Mathematik. 1983, Frankfurt/Main, Germany: Verlag Harri Deutsch, 20MATHGoogle Scholar
Copyright
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.