Amplitude distribution of stochastic oscillations in biochemical networks due to intrinsic noise

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

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.

Biochemical networks with stochastic dynamics can be modelled by the Langevin equation [11]. For the derivation of the theoretical results, we consider a two-dimensional Langevin equation having an equilibrium point at zero:

(1)

with the state vector x ∈ ℛ², the system's dynamics f = (f₁, f₂)^T, where f₁ and f₂ are smooth functions (∈ C^∞), Σ a 2-by-2 matrix of smooth functions and Γ = (Γ₁, Γ₂)^T, where Γ₁ and Γ₂ are uncorrelated, statistical independent Gaussian white noise with zero mean and variance of one.

Let us denote the stationary density distribution of the system (1) by P (x), and the amplitude distribution for the oscillations by P _A (χ), where χ denotes the oscillation amplitude. Under certain assumptions (see following sections) the amplitude distribution P _A (χ) then satisfies

(2)

with ν(x) = || f(x)|| the average speed at state x and (χ) given by

(3)

Hereby A is defined as A = {x ∈ ℛ²|⟨ 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₁ in the positive direction, (3), i.e. = (1, 0)^T, can be simplified to .

To normalize the probability distribution P _A of the amplitude, we have to divide equation (2) through the integral of the probability of all amplitudes:

(4)

Simple models may have additional properties allowing us to simplify (2). For linear systems the speed ν((χ)) is directly proportional to ||(χ)|| simplifying (2) to

(5)

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 x₁, 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₁ during the corresponding oscillation period.

We make the following assumptions on the system (1):

1. There exists a stationary density distribution P: ℛ² → ℛ⁺: 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.

2. The deterministic formulation ( = f(x)) of (1), which is at the state := x(t₀) at time t = t₀ 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(t₀ + Δt)) ≈ ) This assumption can be easily checked by calculating the Lie derivative of P with respect to f:

(6)

which has to be small for almost all x ∈ ² 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 ∈ ℛ² 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.

With the help of these assumptions we are able to formulate a theorem for the calculation of the amplitude distribution of stochastic oscillations. For the formulation of the theorem let us define

(7)

i.e. the unit vector in the direction of x₁ 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:

Figure 1

Definition of. The solid black lines are the level curves of P (x). In ∈ the unit vector in the direction of x₁ is orthogonal to the tangent on the level curve trough .

(8a)

(8b)

(8c)

(8d)

(8e)

The variable ν can be thought of as the deterministic speed, L is the standard deviation of the Lie-derivative of P with respect to the right hand side of (1), a the derivative of P in the direction of x₁, and κ the curvature of P with respect to x₂, each evaluated at the states ∈ . The variables are illustrated in Figure 2.

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 ₀ . Then the amplitude of Ψ _f during the current oscillation will lie with a probability of 70.8% in the set [, + δx₁] and with a probability of 95.5% in the set [ , + ], with δx ₁ defined by .

Lemma 2 Assume the Assumptions 1-3 are satisfied. Then the net flux density(x) of realizations at the state x ∈ ℛ² 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.

Lemma 2 is derived by considering the flux of realizations through an infinitesimaly small region around the state x, as illustrated in Figure 3.

Figure 3

Flux through unit cell. The fluxes through the edges of an infinitesimal small unit cell around a state x ∈ ².

Theorem 1 (Distribution of amplitudes) Under the Assumptions 1-3 and if δx ₁ as defined in Lemma 1 is small, the probability P _A (χ) of an amplitude χ of the stochastic oscillations of (1) in the direction x₁ 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₁ according to the formula in Lemma 1. If δx₁ 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 δx₁ 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₁ 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

In this section we give a motivational example for the application of Theorem 1 by discussing the behavior of the damped harmonic oscillator in the presence of additive Gaussian white noise. The oscillator satisfies the following Langevin equation:

(9a)

(9b)

with Γ _i , i ∈ {1, 2} uncorrelated, statistical independent Gaussian white noise with zero mean and variance of one, σ ∈ ⁺ and k ∈ .

It is easy to show that the deterministic system (σ = 0) is asymptotically stable for k > 0. The eigenvalues of the system matrix

(10)

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

For σ > 0 we can calculate with the help of the multivariate Fokker-Planck equation the density distribution of x _i , i ∈ {1, 2}, in the state space. The general multivariate Fokker-Planck equation (see, for example, [3, p. 193 ff.] or [13, p. 96 ff.]) is given by

(11)

with (x) and B_{i, 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

(12)

with Σ the matrix of the noise terms of (1). We call (11) linear if ∀ _i (x) is a linear function in x and ∀i, j B_{i, j}(x) is constant. Then (11) can be simplified to

(13)

For the system (9), A_{i, j}are given by the system matrix A from (10), and B_{i, j}are the elements of the matrix B = σ²I₂, where I₂ is the second order identity matrix. To get the stationary density distribution, we set . Following [3], the stationary density distribution for the linear Fokker-Planck equation (13) is given by

(14)

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

To get the amplitude distribution P _A (χ) of the oscillations of the harmonic oscillator (9), we apply Theorem 1 after checking the assumptions made in the theorem. Notice that due to the easy structure of this first example and its symmetry, we are always measuring the oscillations in a direction tangential to a principal axis of the distribution. Therefore (χ) = (χ, 0)^Tfor all χ. Because of this it is easy to calculate the necessary variables with the help of the definitions (8a)-(8e) to

(15a)

(15b)

(15c)

(15d)

After some calculations we obtain

(16)

For s = 10 and k = 0.01, δx₁ 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.

Due to the linearity of the oscillator (9) and the symmetry of the harmonic oscillator, we can determine the amplitude distribution by utilizing formula (5) as

(17)

For s = 10 and k = 0.01 the predicted amplitude distribution P _A is plotted over the amplitude χ in Figure 4. For comparison we also displayed the amplitude distribution of the stochastic oscillations as measured in a long term simulation run.

Figure 4

Amplitude distribution of the oscillations of the damped harmonic oscillator. Amplitude distribution of the oscillations of x₂ 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

In the following section we apply the algorithm developed in this paper to a realistic example from biochemical signal transduction. A frequent module in many eucaryotic cells from yeast to mammals is the mitogen activated protein (MAP) kinase signaling cascade. MAPK cascades are typically activated by extracellular stimuli such as growth factors, and regulate the activity of various genes, thereby provoking a cellular response to the applied stimulus. Many important cellular functions such as differentiation, proliferation and death are controlled by MAPK cascades [14]. MAPK cascades consist of three layers of kinases, where each kinase phosphorylates and thereby activates the kinase on the next layer, as shown in Figure 5. The kinases are named MAPKKK, MAPKK, and MAPK, in the order of activation, and the phosphorylated, active forms are denoted with a star. The basic structure may be complemented by additional feedback interconnections, giving rise to deterministic limit cycle oscillations [15]. More recently, oscillations in the MAPK cascade have been determined experimentally in yeast cells, where the oscillatory activity of the MAPK controls periodic changes in cell shape during the mating process [16].

Figure 5

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.

Denote x₁, x₂ and x₃ the concentrations of MAPKKK*, MAPKK* and MAPK*, respectively. Using three conservation relations

(18a)

(18b)

(18c)

we get the deterministic description of the system using mass balancing with Michaelis-Menten reaction kinetics:

(19a)

(19b)

(19c)

with the reaction rates given by

(20a)

(20b)

(20c)

(20d)

(20e)

(20f)

The parameter sets for the phosphorylation (k _ij ) and for the dephosphorylation (p _kl ) are given in Table 1.

Table 1 Parameter Set for the MAP Kinase Signaling Cascade

The biochemically relevant equilibrium point x₀ of the system is computed as

(21a)

(21b)

(21c)

We can determine the local stability of this EP by analyzing the Jacobian A = of the system's dynamics evaluated at the EP. The eigenvalues of A are obtained as

(22a)

(22b)

and thus the deterministic system is locally asymptotically stable.

2.5.2. The Stochastic Model

For validation purposes and to get a first insight into the system's dynamics, we did stochastic simulations of the biochemical network model of the MAPK cascade as described by (19). We used the stochastic simulation software Dizzy by the Institute for Systems Biology [17], which is able to do stochastic simulations of models not only with mass action kinetics, but also with arbitrary kinetics such as Michaelis-Menten kinetics, which are used in our model. For the simulations we assumed a cell volume of V = 1 pl, which corresponds to the absolute number of molecules of each enzyme of

(23a)

(23b)

(23c)

where N _A is the Avogadro constant. We provide a typical plot of the oscillations of x₁, x₂ and x₃ against the time in Figure 6. As can be seen from the plot, the oscillations don't vanish as the deterministic calculations suggested, but are sustained oscillations with a determined amplitude range. In the following we will show that these oscillations can be analyzed with the help of Theorem 1.

Figure 6

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.

As discussed in [11] it is possible, under certain assumptions, to describe a system of biochemical reactions by a set of stochastic differential equations, giving rise to the so called Langevin approximation. These assumptions are, roughly speaking, fulfilled if the amount of entities of each species does not get too small, like in our system. The Langevin approximation of the MAPK cascade model (Figure 5) is given by

(24a)

(24b)

(24c)

with Γ_{i, f}and Γ_{j, b}, i, j ∈ {1, 2, 3} uncorrelated, statistical independent Gaussian white noise with zero mean and a variance of one, N _A the Avogadro constant and V the cell volume. This system corresponds to the deterministic one with an additional noise part for every reaction added, with a standard deviation corresponding to the square root of the magnitude of the reaction rate divided by the cell volume. The factors of the white noises, 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:

(25a)

(25b)

(25c)

with , and for a cell volume of V = 1 pl.

2.5.3. Transformation and Model Reduction

The theoretical results in this paper have been developped for two-dimensional systems only. In higher order systems stochastic oscillations may appear on a two dimensional manifold in the state space. Such a system can be reduced to an order two system by computing the slow manifold, making use of a time scale separation [18]. Following this approach, we first transform the system (25) to

(26a)

(26b)

with z₁ ∈ the coordinate of the fast and z _slow = (z₂, z₃)^T the coordinates of the slow manifold, ψ _fast and ^ψ _slow = (ψ_{slow, 1}, ψ_{slow, 2})^Tvectors of polynomials of order two and higher, ϕ _fast and ϕ _slow matrices of the noise strength and Γ = (Γ₁, Γ₂, Γ₃)^Tthe noise vector of system (25).

The transformation needs to be done in such a way that the absolute values of the real parts of the eigenvalues of A _fast are much larger than these of A _slow . To this end, we compute the eigenvectors corresponding to the eigenvalues (22a)-(22b) as

(27a)

(27b)

The two eigenvectors with non-zero imaginary parts, e₂ and e₃, are both eigenvectors corresponding to eigenvalues having real parts with small absolute values. Therefore it is straightforward to define the desired transformation as

(28)

with the transformation matrix T = (e₁, e₂ + e₃, (e₂ - e₃)·i). We can get the reaction rates of the transformed system by calculating

(29)

Because all of the eigenvalues (22a)-(22b) have non-zero real parts, the Hartman-Grobman theorem states that there exists a local transformation z = H(ζ) such that the system in ζ-coordinates obeys the differential equation

(30)

Because in the ζ-system there exists a slow manifold spanned by the eigenvectors corresponding to the two eigenvalues with real parts having low absolute values, in the z-system there exists a slow manifold, too (see [19]). We now search the function z₁ = m(z₂, z₃) which describes the dependence of the states on the manifold. We approximate this function with a truncated Taylor series expansion:

(31)

We substitute (31) in = v_z,1 and can therewith calculate the coefficients a_{j, 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.

Table 2 Parameter Set of the Slow Manifold of the MAP Kinase Signaling Cascade

Figure 7

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.

The slow manifold is attractive enough so that realizations of the stochastic system (25) will stay close to it. Due to this, it is sufficient to just take the oscillations on the slow manifold into account and therewith simplify the problem to two dimensions. We substitute the formula of the slow manifold (31) in the differential equations for z₂ and z₃ (29) and get a two dimensional reduced description of the system

(32a)

(32b)

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.

The linear approximation of the reduced system (32) is given by

(33)

with the system matrix

(34)

and Γ the vector of the disturbances Γ _i , i ∈ {1, 2, 3}, of the original system (25).

From the coordinate transformation T applied to the original system, Σ_η is determined by

(35)

where diag(σ _i ) is a diagonal matrix with the diagonal elements σ _i being the standard deviations of the stochastic terms in the system (25). For the system (33) with V = 1 pl, Σ_η is obtained as

(36)

Following the same approach as in the previous example, the stationary density distribution is obtained as

(37)

with

(38)

for a cell with a volume of V = 1 pl (see Figure 8).

Figure 8

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

With the preliminary work of the preceding sections we are now able to calculate the amplitude distribution according to Theorem 1. We consider the states x_χ = (χ, α, β)^Twith a, β ∈ not yet specified. These states are transformed into the z-coordinates by

(39)

The probability of the amplitude χ can now be calculated with the formula

(40)

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.

The resulting amplitude distributions in each of the three original coordinate directions x _i , i ∈ {1, 2, 3} are shown in Figure 9. For comparison, the figure also contains measured amplitude distributions, which are the results of a long term simulation via the Direct Gillespie approach with Dizzy with a simulation over 900000 steps.

Figure 9

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 [, + δx₁], 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₁. 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.

To calculate the density distribution of the nonlinear system we have to solve the nonlinear Fokker-Planck equation (11). We therefore applied the algorithm developed in [20]. The shape of the density distribution, which is non-Gaussian, is shown in Figure 10.

Figure 10

Nonlinear Distribution of the state of the MAP kinase signaling cascade. Density distribution of the MAP kinase cascade on the slow manifold (31).

From the density distribution P (x), we can compute the amplitude distribution P _A (χ) with the method proposed in this paper. In this example, we analyze the amplitude distribution P _A (χ) of x₂ for peak amplitudes below the equilibrium point values, i.e. = (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.

Figure 11

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 x₂-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 [, + δx₁] as stated in Lemma 1. By requiring δx₁ to be small we justified to approximate the amplitude by . However, if for specific systems the distribution in [, + δx₁] can be calculated or estimated, it would be possible to reduce the bias for systems where δx₁ 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 ∈ .

In the following we argue that there exists a δx₁ such that for a high probability an arbitrarily chosen realization of (1) being at state at time t won't reach a state x with x₁ > + δx₁ until the next oscillation and therefore the measured amplitude will lie in [, + δx₁] with δx₁ 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:

(41)

After substituting the definitions made in the equations (8) and utilizing the knowledge that vanishes because of the definition of , setting

, we get

(42)

If we stay on a level curve, it must hold that δP = 0 and δx₁ has to be a function of δx₂ locally around . This function is approximated with another Taylor series expansion up to order two as

(43)

with k₁, k₂ yet unknown constants. By substituting (43) in (42) we get

(44)

Because the coefficients of δx₂ and must vanish independently for (44) to hold, it follows that k₁ = 0 and k₂ = , giving

(45)

This means that we may approximate a level curve of the density distribution locally by a parabola. For small deviations and small times t it is possible to make some approximations for (1). First we may approximate Σ (x) in an area around by the constant

(46)

We can rewrite the norm of the first row elements of by

(47)

which can be easily validated by the definitions of L and a given in the equations (8). We furthermore approximate the second element of f(x) in an area around by

(48)

Putting everything together we get a one dimensional Wiener process for the movement of (1) in the direction of x₂ for small times Δt:

(49)

with Δx(t) = x(t) - . We approximate the average displacement in the direction of x₁ by combination of the average displacement of this Wiener process (⟨ x₂(t)⟩ = νt) with equation (45):

(50)

Because the stochastic part of (49) may be neglected (see Assumption 4), we can get a linear one dimensional stochastic equation for the movement of (1) in the direction of x₁ for small times t:

(51)

It can be shown [10, p. 129 ff.] that the variance of (51) evolves as

(52)

We may think of the solution of (51) as a growing Gaussian distribution moving along the trajectory of the deterministic system (see Figure 2). We now determine the maximal value of Δx₁(t) that a realization of (1) starting at may have if it evolves in an area around the average displacement determined by the standard deviation . The maximal value of Δx₁(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:

(53)

(54)

We can now calculate δx₁ to

(55)

(56)

For a Gaussian distribution it is true that 68.3% of all realizations stay in an area around the average displacement determined by the standard deviation. This means that only of the realizations reach an amplitude greater than + δx₁ after passing . Because it is also possible for a realization to have its maximal value in the direction of x₁ before passing , the overall probability for a realization to have an amplitude lower than + δx₁ is

(57)

In the same way it can be shown that 95.5% of all realizations have an amplitude lower than + .

Proof for Lemma 2

Assume a small unit cell with the edge lengths dr and dφ with its center being x and the level curve of P through x given and going through the center of the edges with the lengths dr (see Figure 3). If dr and dφ are small, we may assume the density distribution P and the values f(x) and Σ(x) of (1) as constant inside the unit cell. Due to the small size of the Lie-derivative of P along f, the movement normal to the level curve is dominated by the stochastic part of (1). The net fluxes ϕ₁ and ϕ₃ trough the edges 1 and 3 are vanishing and therefore the overall flux is tangential to the level curve. The net flux ϕ₂ (ϕ₄), trough the edges 2 (the edge 4) can be determined by integration (see, for example, [[25], p. 622]):

(58)

(59)

(60)

By letting dr and dφ go to zero the absolute value of the net flux density ||(x)|| follows:

(61)

Proof for Theorem 1

Because the density distribution P is smooth and the curvature of its level curves does not change its sign, there exists exactly one state ∈ with = χ for every χ ∈ (0, ∞). is a smooth simple curve in the state space which can be parameterized by

(62)

with 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:

(63)

with 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 ∈ :

(64)

From the condition we get

(65)

Because of Lemma 2 we know that the flux density is directed along the x₂-axis and has the absolute value P () ν(). Therefore we can calculate (64) to

(66)

(67)

Although we do not know the exact shape of Ω (, χ), 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

(68)

If δx₁ is small, we can assume P and ν being approximately constant in the integration interval, which simplifies (68) to

(69)

For small δx₁ the integral of Ω over is approximately 1. The result of this last approximation is given by Theorem 1.