Amplitude distribution of stochastic oscillations in biochemical networks due to intrinsic noise
- Moritz Lang1,
- Steffen Waldherr1Email author and
- Frank Allgöwer1
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 long-term 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 long-term 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 non-linear 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
Keywords
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 Gaussian-like 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 multi-scale 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 Fokker-Planck 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 non-linear 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 non-linear 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 long-term 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 = (f1, f2) T , where f1 and f2 are smooth functions (∈ C∞), Σ a 2-by-2 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.
(χ) given by
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 x1 in the positive direction, (3), i.e.
= (1, 0)T, can be simplified to
.
(χ)) is directly proportional to ||
(χ)|| simplifying (2) to
with
given as above.
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 x1, 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 x1 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.
= f(x)) of (1), which is at the state
:= x(t0) at time t = t0 with the probability of the state
being P (
) =:
, will approximatively evolve tangential to the level curve of the probability distribution, thus for small Δt the probability of the state at time t + Δt won't have changed significantly (P (x(t0 + Δt)) ≈
) This assumption can be easily checked by calculating the Lie derivative of P with respect to f: 
which has to be small for almost all x ∈
2 having a non-negligible 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 non-negligible 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.
in the direction of x1 is orthogonal to the tangent on the level curve of the probability distribution in
(see Figure 1). Because of Assumption 1 there exists exactly one state
∈
for every amplitude χ satisfying
= χ. Furthermore we define the following variables, which can be calculated with the knowledge of the system (1) and the stationary probability distribution P:
Definition of
. The solid black lines are the level curves of P (x). In
∈
the unit vector
in the direction of x1 is orthogonal to the tangent on the level curve trough
.
∈
. The variables are illustrated in Figure 2.
Visualization of the variables. Visualization of the variables needed to calculate the amplitude distribution of stochastic oscillations.
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 1-3 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 [
,
+ δx1] and with a probability of 95.5% in the set [ , + ], with δx
1
defined by .
Lemma 2 Assume the Assumptions 1-3 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.
Flux through unit cell. The fluxes through the edges of an infinitesimal small unit cell around a state x ∈
2.
Theorem 1 (Distribution of amplitudes) Under the Assumptions 1-3 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 x1 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 δx1 according to the formula in Lemma 1. If δx1 is significantly smaller than
for almost all
∈
having a non-negligible probability P (
), Theorem 1 may be applied. For a general stochastic differential equation, it is hard to exactly indicate how much smaller δx1 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 δx1 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
(x) and Bi, j(x) > 0 real functions describing the influence of the system dynamics and the noise on the distribution. The diffusion term B may be calculated to
(x) is a linear function in x and ∀i, j Bi, j(x) is constant. Then (11) can be simplified to
. Following [3], the stationary density distribution for the linear Fokker-Planck equation (13) is given by
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
(χ) = (χ, 0)
T
for all χ. Because of this it is easy to calculate the necessary variables with the help of the definitions (8a)-(8e) to
For s = 10 and k = 0.01, δx1 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.
Amplitude distribution of the oscillations of the damped harmonic oscillator. Amplitude distribution of the oscillations of x2 of the damped harmonic oscillator (9) in the presence of white noise with autocorrelation σ = 10 and k = 0.01. The solid curve corresponds to the data from simulation, the dashed curve is the prediction according to Theorem 1. The simulation was run over a time scale of 100000 s which corresponds to around 16000 oscillation periods.
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
Schematic view on the MAP kinase signaling cascade. Schematic view on the MAP kinase signaling cascade. The species with the stars (*) are the phosphorylated versions of the species without stars. MAPKKK* catalyzes the phosphorylation of MAPKK, MAPKK* in turn catalyzes the phosphorylation of MAPK. A negative feedback is given by the repression of the activation of MAPKKK by MAPK*.
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 |
of the system's dynamics evaluated at the EP. The eigenvalues of A are obtained as
and thus the deterministic system is locally asymptotically stable.
2.5.2. The Stochastic Model
Oscillations of the MAP kinase signaling cascade. Plot of the oscillations of the MAP kinase signaling cascade. The plot is the result of a stochastic simulation via the Direct Gillespie approach with Dizzy.
and
with i, j ∈ {1, 2, 3}, are monotone in x and don't change significantly in an area around the EP, so that we can approximate them with their steady state values. Afterwards we can combine each function's noise terms by adding their variances
, leading to the following equations:
with
,
and
for a cell volume of V = 1 pl.
2.5.3. Transformation and Model Reduction
with z1 ∈
the coordinate of the fast and z
slow
= (z2, z3)
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).
= vz,1 and can therewith calculate the coefficients aj, k(see Table 2). The resulting manifold in the original x-coordinates is shown in Figure 7. The figure also includes an example trajectory. As can be seen, this trajectory first converges exponentially fast to an ϵ-neighborhood of the slow manifold and afterwards moves on it towards the EP.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 |
Slow Manifold of the MAP kinase signaling cascade. Slow manifold of (20) calculated up to the order of four. The black curve is a representative trajectory of the deterministic system.
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 Fokker-Planck equations numerically.
and Γ the vector of the disturbances Γ i , i ∈ {1, 2, 3}, of the original system (25).
Linear distribution of the state of the MAP kinase signaling cascade. Linear density distribution of the MAP kinase signaling cascade on the slow manifold as calculated in (37). Because we linearized the system, the distribution is Gaussian.
2.5.5. Determination of the Amplitude Distribution
not yet specified. These states are transformed into the z-coordinates by
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.
Linear amplitude distribution of the oscillations of the MAP kinase signaling cascade. Linear amplitude distribution of the oscillations of the MAP kinase signaling cascade. The solid curve corresponds to data determined by stochastic simulation of the nonlinear system (24), the dashed curve is the prediction according to the calculations in this section.
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 [
,
+ δx1], 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 δx1. 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.
Nonlinear Distribution of the state of the MAP kinase signaling cascade. Density distribution of the MAP kinase cascade on the slow manifold (31).
= (0, -1, 0)T, because then the effect of the nonlinearities is even higher than in the other direction. The reason for this is quite simple: The amount of entities of each species is not allowed to become negative. As can be seen in Figure 11 the amplitude distribution, as obtained by applying Theorem 1, quite well fits the predicted results gained from a long time simulation while the predicted values of the linearized system are far away from these values.
Nonlinear amplitude distribution of the oscillations of the MAP kinase signaling cascade. Nonlinear amplitude distribution of the oscillations of the MAP kinase signaling cascade in the negative x2-direction. The solid curve corresponds to data experimentally measured, the dashed red curve is the prediction according to the calculations in this section and the dash-dotted green curve is the outcome we get by calculating the amplitude distribution of 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 [
,
+ δx1] as stated in Lemma 1. By requiring δx1 to be small we justified to approximate the amplitude by
. However, if for specific systems the distribution in [
,
+ δx1] can be calculated or estimated, it would be possible to reduce the bias for systems where δx1 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
∈
.
at time t won't reach a state x with x1 >
+ δx1 until the next oscillation and therefore the measured amplitude will lie in [
,
+ δx1] with δx1 positive and small but yet not further determined (see Figure 2). We first calculate the Taylor series expansion for P up to the order of two:
After substituting the definitions made in the equations (8) and utilizing the knowledge that
vanishes because of the definition of
, setting
, we get
. This function is approximated with another Taylor series expansion up to order two as
must vanish independently for (44) to hold, it follows that k1 = 0 and k2 =
, giving
by the constant
by
by
. We approximate the average displacement in the direction of x1 by combination of the average displacement of this Wiener process (⟨ x2(t)⟩ = νt) with equation (45):
may have if it evolves in an area around the average displacement determined by the standard deviation
. The maximal value of Δx1(t) is time dependent and grows in the beginning due to the stochastic part of (51), but afterwards shrinks due to the deterministic part. It reaches its maximal value when the time derivatives of the mean value and of the standard deviation have the same absolute value:
of the realizations reach an amplitude greater than
+ δx1 after passing
. Because it is also possible for a realization to have its maximal value in the direction of x1 before passing
, the overall probability for a realization to have an amplitude lower than
+ δx1 is
In the same way it can be shown that 95.5% of all realizations have an amplitude lower than
+
.
Proof for Lemma 2
(x)|| follows:
Proof for Theorem 1
∈
with
= χ for every χ ∈ (0, ∞).
is a smooth simple curve in the state space which can be parameterized by
an unknown but smooth function in χ. We can get a term proportional to the average number of oscillations of all realizations per time unit
if we determine the total net flux of realizations ϕ
t
through this curve:
the normalized normal vector on
in
. We may determine the probability P
A
(χ) of the amplitude χ ∈ (0, ∞) by multiplying the flux density of realizations
with Ω (
, χ), the probability that a oscillation of a realization contributing to the flux density at
has the amplitude, and integrating over all
∈
:
we get
is directed along the x2-axis and has the absolute value P (
) ν(
). Therefore we can calculate (64) to
, χ), which may be complicated to calculate, we know a lower and an upper bound for Ω. The lower bound is given due to the fact that Ω must be zero for χ <
, because if (1) goes through the state
, it has at least the amplitude of this state. An upper bound for almost all realizations is given by Lemma 1 with
+
. We may therefore conclude that Ω is approximately zero for χ being outside of the interval [
,
+
], which simplifies (67) to
For small δx1 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): 7051-7056. 10.1073/pnas.1132112100.View ArticleADSGoogle Scholar
- Tiedemann HB, Schneltzer E, Zeiser S, Rubio-Aliaga I, Wurst W, Beckers J, Przemeck GKH, de Angelis MH: J Theor Biol. 2007, 248: 120-129. 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, Vanden-Eijnden E: Chaos. 2008, 18: 015111-1-015111-11. 10.1063/1.2779852.View ArticleADSMathSciNetGoogle Scholar
- El-Samad H, Khammash M: Proceedings of the 45th IEEE Conference on Decision and Control. 2006, 2382-2387. 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: 459-469. 10.1016/j.jtbi.2006.10.029.View ArticleMathSciNetGoogle Scholar
- Yu N, Kuske R, Li Y: Physical Review E. 2006, 73: 056205-1-056205-8.ADSGoogle Scholar
- Aparicio JP, Solari HG: Mathematical Biosciences. 2001, 169: 15-25. 10.1016/S0025-5564(00)00050-X.View ArticleMathSciNetMATHGoogle Scholar
- Gillespie DT: Markov Processes. 1992, San Diego, CA: Academic Press, INCMATHGoogle Scholar
- Gillespie DT: Journal of Chemical Physics. 2000, 113: 297-306. 10.1063/1.481811.View ArticleADSGoogle Scholar
- Bagheri N, Stelling J, Doyle FJ: Bioinformatics. 2007, 23 (3): 358-364. 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): 153-183. 10.1210/er.22.2.153.Google Scholar
- Kholodenko B: European Journal of Biochemistry. 2000, 67 (6): 1583-1588. 10.1046/j.1432-1327.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): 1700-06. 10.1016/j.cub.2008.09.027.View ArticleGoogle Scholar
- Ramsey S, Orrell D, Bolouri H: Journal of Bioinformatics and Computational Biology. 2005, 3: 415-436. 10.1142/S0219720005001132.View ArticleGoogle Scholar
- Khalil HK: Nonlinear Systems. 2001, Upper Saddle River, NJ, USA: Prentice-Hall, 3Google Scholar
- Jacobs SJ: Journal of Atmospheric Sciences. 1991, 48 (7): 893-901. 10.1175/1520-0469(1991)048<0893:EOASMI>2.0.CO;2.View ArticleADSGoogle Scholar
- Sjöberg P: Numerical Solution of the Fokker-Planck Approximation of the Chemical Master Equation. PhD thesis. 2005, Uppsala UniversitetGoogle Scholar
- Rao CV, Wolf DM, Arkin AP: Nature. 2002, 420 (6912): 231-237. 10.1038/nature01258.View ArticleADSGoogle Scholar
- Steuer R, Zhou C, Kurths J: Biosystems. 2003, 72 (3): 241-251. 10.1016/j.biosystems.2003.07.001.View ArticleGoogle Scholar
- Waldherr S, Allgöwer F: International Journal of Systems Science. 2009, 40: 769-782. 10.1080/00207720902957269.View ArticleMATHGoogle Scholar
- Kitano H: Nat Rev Genet. 2004, 5 (11): 826-837. 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.
