Email updates

Keep up to date with the latest news and content from BMC Systems Biology and BioMed Central.

Open Access Methodology article

Simulation methods with extended stability for stiff biochemical Kinetics

Pau Rué12, Jordi Villà-Freixa1* and Kevin Burrage34*

Author affiliations

1 Computational Biochemistry and Biophysics Group, Research Unit on Biomedical Informatics, IMIM/Universitat Pompeu Fabra, c/Dr. Aiguader 88, 08003, Barcelona, Catalonia, Spain

2 Departament de Física i Enginyeria Nuclear, Universitat Politècnica de Catalunya, Edifici GAIA, Rambla de Sant Nebridi s/n 08222, Terrassa, Barcelona, Spain

3 Institute for Molecular Bioscience, The University of Queensland, Brisbane, QLD, 4072, Australia

4 COMLAB and OCISB, University of Oxford, Oxford OX1 3QD, UK

For all author emails, please log on.

Citation and License

BMC Systems Biology 2010, 4:110  doi:10.1186/1752-0509-4-110


The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1752-0509/4/110


Received:19 November 2009
Accepted:11 August 2010
Published:11 August 2010

© 2010 Rué et al; licensee 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.

Abstract

Background

With increasing computer power, simulating the dynamics of complex systems in chemistry and biology is becoming increasingly routine. The modelling of individual reactions in (bio)chemical systems involves a large number of random events that can be simulated by the stochastic simulation algorithm (SSA). The key quantity is the step size, or waiting time, τ, whose value inversely depends on the size of the propensities of the different channel reactions and which needs to be re-evaluated after every firing event. Such a discrete event simulation may be extremely expensive, in particular for stiff systems where τ can be very short due to the fast kinetics of some of the channel reactions. Several alternative methods have been put forward to increase the integration step size. The so-called τ-leap approach takes a larger step size by allowing all the reactions to fire, from a Poisson or Binomial distribution, within that step. Although the expected value for the different species in the reactive system is maintained with respect to more precise methods, the variance at steady state can suffer from large errors as τ grows.

Results

In this paper we extend Poisson τ-leap methods to a general class of Runge-Kutta (RK) τ-leap methods. We show that with the proper selection of the coefficients, the variance of the extended τ-leap can be well-behaved, leading to significantly larger step sizes.

Conclusions

The benefit of adapting the extended method to the use of RK frameworks is clear in terms of speed of calculation, as the number of evaluations of the Poisson distribution is still one set per time step, as in the original τ-leap method. The approach paves the way to explore new multiscale methods to simulate (bio)chemical systems.

Background

It is by now very well known that the biochemical kinetics involving small numbers of molecules can be very different to kinetics described by the law of mass action and differential equations [1-3]. This effect is a property of the intrinsic noise of the system and is associated with the uncertainty of knowing when a reaction occurs and what that reaction is. At the molecular level such intrinsic uncertainty is, in turn, a consequence of the stochastic nature of the fluctuations of the potential energy surface for any chemical reaction in the condensed phase [4]. When considering a collection of molecules, the intrinsic noise is accentuated when some chemical species have small numbers, as is often the case in genetic regulatory models where there are small numbers of key transcription factors that can bind to a limited number of operator regions on DNA [5-15]. Kurtz [16] and Gillespie [17] realised this fact and developed discrete methods to deal with this situation. The stochastic simulation algorithm (SSA, see [18] for a review) describes the time evolution of the dynamics of the species in a well-stirred chemically reacting system as a discrete nonlinear Markov process, resulting in an exact method to sample from the probability density function described by the chemical master equation (CME). Gibson and Bruck proposed a more efficient implementation of the SSA called the next reaction method [19].

The basic idea of the SSA is that at each time point a waiting time to the next reaction and the most likely reaction to occur must be sampled from a joint probability density function leading to an appropriate update of the state vector. But if the rate constants and/or the numbers of molecules in the system are large then the waiting time (time step, τ) can be very small [18]. Because of this Gillespie [20] introduced the Poisson τ-leap method, in which all reactions are allowed to fire in a given τ with a frequency extracted from a Poisson distribution. Since then many extensions of this idea have been developed. Cao et al. [21] have considered efficient mechanisms for selecting τ and have developed implicit methods suitable for simulating stiff systems. Tian and Burrage [22] introduced a modification of Poisson τ-leap methods known as Binomial τ-leap methods that avoids the issue of obtaining negative molecular numbers from which Poisson τ-leap methods can suffer. Chatterjee et al. [23] and Auger et al. [24] have considered modifications to Binomial τ-leap methods that improve some of the implementation aspects. On the other hand, Monk [25] and Mackey [26] noted the importance of representing delays, especially when representing processes such as transcription and translation. Accordingly, Bratsun et al. [12] and Barrio et al. [27] developed a delayed version of the Stochastic Simulation Algorithm. Leier et al. [28] and Anderson [29] extended these ideas to a τ-leap setting.

Although τ-leap methods can, in some cases, substantially improve computational efficiency compared with the SSA, when there is moderate stiffness in the system the efficiencies can be quite poor. One could resort to implicit τ-leap methods but then there are considerable implementation issues and subtleties. A different approach is to explore ideas from the numerical ODE (ordinary differential equations) and numerical SDE (stochastic differential equations) communities. Thus, with ODEs it is well known that stiffness leads to a step size restriction when using explicit methods and many classes of efficient implicit methods have been designed [30]. However, in the case of moderately stiff systems explicit Runge-Kutta methods with extended stability regions along the negative real axis have proven to be especially effective [31,32]. Runge-Kutta methods are a class of one step methods which gain their efficacy by computing intermediate approximations to the solution within a step. Explicit Runge-Kutta methods with extended stability regions are based on explicit Runge-Kutta methods whose stability function is a shifted and scaled Chebyshev polynomial or some variant thereof. In the stochastic setting, there are some subtleties designing fully implicit methods due to possible unboundedness of the solution as the Wiener increment can take positive or negative values with equal likelihood [33]. Thus most methods are semi-implicit, that is implicit in the deterministic component. Abdulle and Cirilli [32] have, with some success, extended the ideas of explicit Chebyshev methods with extended stability regions to the SDE setting via their class of S-ROCK methods.

Here, we use the Runge-Kutta formulation to construct methods with large stability regions so that efficiencies are gained by allowing larger stepsizes. We note that this is exactly what Abdulle and Cirilli [32] do in the SDE setting, that is they use a Runge-Kutta formulation to construct methods with excellent stability properties and even though these methods are only weak order 1 they perform very well. It is noteworthy that in this work we are not using the Runge-Kutta formulation to get second order accuracy for τ-leap methods. This seems to be a difficult problem, just as it is the case for SDEs and will probably require double integrals of compensated processes to be simulated. In fact, Abdulle and Cirilli [32] also note that it is very difficult to construct weak order 2 methods with good stability properties and to our knowledge at the moment no such methods exist in the SDE setting. Note that in a stochastic setting we judge order of accuracy through two mechanisms: strong order (where trajectories are compared with the true solutions) and weak order (where moments are compared). Often a numerical method may have a higher weak order than its strong order. The Euler-Maruyama method is a case in point with weak order one and strong order a half.

Thus, in this paper, we explore a series of fully explicit multistage Runge-Kutta methods with extended stability for a fixed τ-leap stochastic simulation schema. Our methods involve the same number of Poisson evaluations per integration step as in the original τ-leap formulation but allow increasingly larger step sizes at the cost of an increasing series of deterministic evaluations in the internal stages. First we give some background on Runge-Kutta methods for ODEs and SDEs. In section Results we extend these ideas to the τ-leap methods and present a stability analysis for linear chemical kinetics, including its practical implementation. In section Numerical results we present numerical results for both the linear case and the classical stiff system described by the Schlögl reaction [34]. Finally, in section Discussion we discuss further implications of this work and, in particular, possible extensions to multiscale modelling.

Review of Runge-Kutta methods for SDEs and ODEs

Stability region for RK methods applied to ODEs

Consider the system of initial value ODEs given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M1">View MathML</a>

(1)

The class of s-stage Runge-Kutta (RK) methods for approximating the solution to (1) is given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M2">View MathML</a>

(2)

where h is the time step. This class of methods is characterised by the Butcher tableau

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M3">View MathML</a>

where bT = (β1,...,βs), w = Ae and e = (1,...,1)T. Here A is the matrix with entries αij and w is the column vector wT = (w1,...,ws)T. A Runge-Kutta method is said to be explicit if the s × s matrix A is strictly lower triangular. The method parameters are usually chosen so that a Runge-Kutta method has appropriate efficiency, order and stability characteristics. The Yi are considered to be approximations to the solution at the intermediate points tn + wih for i = 1,...s.

In a stability setting an RK method is often applied to the linear, scalar test equation

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M4">View MathML</a>

(3)

In which case it is easily seen that (2) gives rise to

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M5">View MathML</a>

Where

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M6">View MathML</a>

(4)
.

Here R(z) is the so-called stability function. This function can be extended to a linear N-dimensional equation y' = Λy in which case it becomes a matrix function of the N × N matrix Λ:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M7">View MathML</a>

(5)

where e is the unit vector, Is is the identity matrix of order s and ⊗ represents the Kronecker tensor product such that the (i, j) element of A B is aijB. Notice that, if Λ is a scalar value and taking z = hΛ, R(z) would be a scalar and take the form (4). Therefore we can refer to R seamlessly irrespective of whether the argument is a matrix or a scalar.

In the case of an explicit method, as A is a strictly lower triangular s × s matrix, its sth power is As = 0. Therefore, equation (4) can be expanded into a finite power series for A:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M8">View MathML</a>

(6)

where rj = bTAj-1e, j = 1,...,s. Hence, R(z) is a polynomial of at most degree s for any explicit method.

Since (3) is asymptotically stable for all Re [λ] < 0, the stability region of a Runge-Kutta method is defined as

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M9">View MathML</a>

(7)

Stability region for RK methods applied to SDEs

In the case of stochastic differential equations (SDEs), we consider the general m dimensional form

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M10">View MathML</a>

(8)

where W(t) = (W1(t),...,Wd(t))T is a vector of d independent Wiener processes in which an individual Wiener process has the properties

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M11">View MathML</a>

and non-overlapping Wiener increments are independent of one another. A sample of a Wiener increment W(t + h) - W(t) is simulated from a Normal random variable with mean 0 and variance h, N(0, h).

Equation (8) can arise as the limit of a discrete process through the concept of a diffusion process in which case f (t, y) will represent the mean of this process and g(t, y) is the m × d matrix such that ggT is the covariance. Equation (8) can be interpreted in several ways (see [35] for an introduction to SDEs), depending on which integral definition is used. Two such interpretations lead to Itô and Stratonovich forms of SDEs. In the Itô setting an integral is approximated by summing, over a partition, the areas of a rectangle with width the increment of the Wiener process on that subinterval and height the value of the integrand at the lefthand point of each subinterval whereas in the Stratonovich setting the integrand is evaluated at the midpoint of each interval. If (8) is interpreted in the Itô sense then the simplest numerical algorithm is given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M12">View MathML</a>

(9)

where ΔWn = (ΔW1,....ΔWd)T and ΔWi := Wi(tn + h) - Wi(tn ), i = 1,...,d are normally distributed random numbers with mean 0 and variance h. This method is known as the Euler-Maruyama method and it is known to have strong order (pathwise order) <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M13">View MathML</a> and weak order (moment order) 1.

As with the deterministic case, the quality of a numerical method can be partly characterised by its stability region associated with the scalar, linear test equation

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M14">View MathML</a>

(10)

The solutions of (10) in the Itô and Stratonovich cases are, respectively,

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M15">View MathML</a>

In the later case, the solution is mean square stable <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M16">View MathML</a> if Re [a] + Re [b2] ≤ 0.

A very general class of stochastic Runge-Kutta methods [36] was constructed for the solution of (8) which, when applied to the scalar test SDE (10) produces

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M17">View MathML</a>

where R is a multinomial in p and q if the method is explicit and where p = ha, <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M18">View MathML</a>. Analogous to the deterministic case, the mean square stability region of a method is defined as

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M19">View MathML</a>

In the case of the Euler-Maruyama method

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M20">View MathML</a>

and in the (p, q) plane, with p, q ∈ ℝ, the stability region is a circle of radius 1 centered in (-1,0).

Results

The τ-leap Runge-Kutta framework with bounded variance and extended stability domain

As stated in the Background section, the SSA describes the time evolution of a vector of integer numbers of molecules in the presence of intrinsic noise. More formally, suppose that there are N chemical species S1,...,SN undergoing m chemical reactions. Let Xi(t), i = 1,...,N denote the number of molecules of species Si and X(t) = (X1(t),...,XN(t))T. Now any set of chemical reactions is uniquely characterised by two sets of quantities. These are the update (stoichiometric) vectors ν1,...,νm for each of the m reactions and the propensity functions a1(X(t)),...,am(X(t)), which are proportional to the probabilities of each of the reactions occurring. For example, given the reaction

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M21">View MathML</a>

then X(t) = (A(t), B(t), C(t))T, ν1 = (-1, -1, 1)T and a1(X(t)) = cA(t)B(t).

Given X(t) at time t, the SSA determines a waiting time τ to the next reaction assuming an exponential waiting time distribution <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M22">View MathML</a>, where <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M23">View MathML</a>, and then selects the most likely reaction, say k, based on the relative sizes of a1(X(t)),...,am(X(t)). The state vector is then updated as

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M24">View MathML</a>

and the algorithm repeats.

Since a typical stepsize (waiting time) is of the size 1/a0(X(t)), this can be very small if some of the rate constants are large and/or some species have large numbers of molecules. Accordingly τ-leap methods attempt to take a larger step size in which all the reactions can occur based on a certain frequency. This can be written as

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M25">View MathML</a>

(11)

Gillespie [20] chose the number of Rj reactions per step, Kj , as coming from a Poisson distribution with mean τaj(Xn), that is

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M26">View MathML</a>

(12)

Using the so-called compensated process given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M27">View MathML</a>

(13)

which satisfies <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M28">View MathML</a>[L (τ, x)] = 0 and <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M28">View MathML</a> [L (τ, x)2] = τx, equation (11) can be restated as

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M29">View MathML</a>

(14)

Where <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M30">View MathML</a>.

As noted by Gillespie [20] and Tian and Burrage [22], and as a consequence of the Law of Large Numbers, as → ∞, L(τ, x) converges to a normal random variable with zero mean and variance τx, N(0, τx), and this can be considered as a sample <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M31">View MathML</a> of <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M32">View MathML</a>. Substituting this into (14) gives

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M33">View MathML</a>

(15)

This is precisely the Euler-Maruyama method applied to the SDE

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M34">View MathML</a>

(16)

Thus in the continuous limit the Poisson τ-leap method can be viewed as the Euler-Maruyama method applied to a form of the Chemical Langevin Equation. Indeed Li [37] has shown that the Poisson τ-leap method has mean square strong order <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M35">View MathML</a> and weak order 1 and this is consistent with the previous remarks. In addition, equation (16) is a particular case of the general SDE

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M36">View MathML</a>

These relationships naturally lead to the introduction of the class of Runge-Kutta τ-leap methods which bears a relationship, similar to the one discussed above, to the general class of Stochastic Runge-Kutta methods for solving SDEs [36]. This general class of explicit s-stage Runge-Kutta τ-leap methods takes the form

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M37">View MathML</a>

(17)

where L(τ, x) is given by (13) and <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M38">View MathML</a> represents the drift or expected stepchange. As our focus is explicit methods, the matrix A is strictly lower diagonal. We note that (17) requires the same number of samples of Poisson random variables per step as the Poisson τ-leap method.

The Poisson τ-leap method given by (11) and (12) is equivalent to (17) with

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M39">View MathML</a>

Indeed any Runge-Kutta method for solving an ODE can be incorporated into this framework. We also note that other methods proposed in the literature can be put into this framework. For example, the midpoint method of Gillespie [20] can be represented with s = 2, bT = (0, 1), w = (0, 0.5)T and where the row-wise entries of A are 0, 0, 0.5, 0.

The linear case

As in the case of stability settings in the ODE and SDE regimes, we analyse (17) when applied to linear kinetics, which in this case are described by sets of unimolecular reactions. A general set of m unimolecular reactions can be described by m propensity functions given by the following linear functions

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M40">View MathML</a>

(18)

where x is the state vector of dimension N and <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M41">View MathML</a>, j = 1,...,m are m vectors of dimension N defining the propensities. A more convenient way to describe this linear kinetics system is by using the N × N matrix W

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M42">View MathML</a>

so that now the drift or expected step-change can be represented as

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M43">View MathML</a>

If the Runge-Kutta method for ODEs underlying a Runge-Kutta τ-leap method (17) has stability function given by (4), then when the latter is applied to (18) we show (Additional file 1) that

Additional file 1. Supplementary material for "Simulation methods with extended stability for stiff biochemical kinetics". Technical results, coefficients for the optimal stability polynomials and notes on the mitogen-activated protein kinase (MAPK) cascade simulation results.

Format: PDF Size: 310KB Download file

This file can be viewed with: Adobe Acrobat ReaderOpen Data

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M44">View MathML</a>

(19)

where R is the multidimensional version of (4) given by (5). Note that this is a natural generalization of the deterministic case when a Runge-Kutta method is applied to the problem y' = Λy giving Xn = R(hΛ)Xn-1. Thus with fixed stepsize τ

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M45">View MathML</a>

(20)

Therefore, boundedness in the mean requires that the spectral radius, ρ, of R(τW) satisfies

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M46">View MathML</a>

In order to analyse the framework (17) from the perspective of both mean and variance behaviour we consider the reversible isomerisation reaction with fixed total number of molecules given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M47">View MathML</a>

(21)

as the linear scalar test equation. It is easy to see that this system is a analogous to (3) for ODEs and (10) to SDEs with constant nonzero term. The system is chosen to have constant nonzero term in order to compare its variance, which otherwise would fade to zero, to the variance given by the framework methods (17). In this case

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M48">View MathML</a>

(22)
.

For this set of reactions, the Chemical Master Equation (which describes the probability density function associated with the evolving Markov process X) can be solved analytically [18,38]. In particular, it can be shown that the stationary state <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M49">View MathML</a> has a probability density function (PDF) that follows a binomial distribution with

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M50">View MathML</a>

Where

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M51">View MathML</a>

and T = X1(t) + X2(t) is the (fixed) total number of molecules in the system. Thus from the properties of the binomial distribution with e = (1, 1)T

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M52">View MathML</a>

(23)

In the case of non-negative coefficients in the underlying RK method and for constant τ one can show (see details in Additional file 1) that if (17) is applied to (21) with constant τ such that |R(z)| < 1, z = -τ(k1 + k2), then in the limit as n → ∞ the mean vector converges to the theoretical mean, that is

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M53">View MathML</a>

Note that with the constraint |R(z)| < 1, z = -τ(k1 + k2) then the spectral radius of R(τW) is less than or equal to 1, and as there is only one eigenvalue equal to one hence we have boundedness of the mean.

Furthermore, if Var [X] denotes the variance of the new method at steady state (X1 and X2 have the same variance) and if R2(z) ≠ 1, z = -τ(k1 + k2), then (see details in Additional File 1)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M54">View MathML</a>

where

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M55">View MathML</a>

(24)

We call this the relative variance at the stationary state associated to R.

Let us consider some particular cases of this result:

Poisson τ-leap For this method R(z) = 1 + z and <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M56','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M56">View MathML</a>. Thus, the equilibrium variance doubles at z = -1, it rises fourfold at z = -1.5 and is unbounded at z = -2.

Two stage methods with α21 ≠ 0 For the family of explicit two-stage methods with α21 ≠ 0

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M57">View MathML</a>

the stability function is R(z) = 1 + z + γz2, where γ = β2α21 and the variance behaviour is determined by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M58">View MathML</a>

In this case we have one free parameter of the method, γ, which allows us to control both the stability function R and the relative variance at steady state. We might be interested in setting γ to a value that both allows large time-steps to be used (by maximising the region (-l, 0] for which z fulfils |R(z)| < 1) and keeps the relative variance, ψ(z) close to one. In the case γ <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M59">View MathML</a>, ψ grows as z becomes more negative. More interesting is the case γ > <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M59">View MathML</a>, where the maximum and minimum of ψ occur for <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M60">View MathML</a>, respectively and in this case

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M61','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M61">View MathML</a>

Constraining ψ to be around 1 with a certain fixed tolerance ϵ, |ψ(z) -1| < ϵ, for a range z ∈ (-l, 0] to be maximised is achieved with

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M62">View MathML</a>

and with a stability region (-l, 0] with

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M63">View MathML</a>

For instance, for 0.5 < ψ(z) < 1.5, setting γ = 0.20096 gives a maximum stability region of (-3.68026, 0] and thus the method

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M64','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M64">View MathML</a>

This is the methodology we propose in the following section for the derivation of particular Runge-Kutta methods with s steps. Note that if we required the same limitation on the variance with the standard Poisson τ-leap method we could only take <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M65','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M65">View MathML</a>. Thus with the two stage method we can take a stepsize almost six times as large.

Implicit midpoint rule For the implicit midpoint rule

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M66','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M66">View MathML</a>

(25)

This was first shown by Cao et al. [38]. In fact only those Runge-Kutta methods that have a stability function given by (25) can preserve the variance exactly for linear problems. These methods include the implicit midpoint and trapezoidal rules and have to be implicit.

Methods with bounded variance and extended stability domain

For the general case of s stages we require ψ(z) to be as close to 1 as possible for as large a range of z as possible, this is, for as large a range of z fulfilling the stability condition |R(z)| < 1). We proceed by first showing that if we consider a bound on the relative variance, ψ, around one, we automatically fulfil the stability conditions for a certain range. In this sense, let ϵ ≥ 0 (and ϵ < 1), we impose the constraint

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M67','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M67">View MathML</a>

(26)

and optimise the value of ls, ϵ such that the range for which this holds is (-ls, ϵ, 0].

Noticing from (24) that

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M68">View MathML</a>

(27)

inequality (26) can be restated in terms of R(z)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M69','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M69">View MathML</a>

(28)

with z ∈ (-ls ϵ, 0]. Hence, we can translate constraints in the relative variance into constraints in the stability function. Since we are interested in constructing explicit methods we can ask how we can make ψ(z) close to 1 in an explicit framework for which we already know the stability function is a polynomial of at most degree s (equation (6))

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M70','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M70">View MathML</a>

Thus, similar to the case s = 2 in which we had one free parameter, γ, to optimise, if we assume r1 = bTe = 1 then we have s - 1 parameters, r2,...,rs, we can optimise. In this case, though, the search of the optimal set of parameters has to be performed with numerical optimisation methods rather than analytically. The problem of finding optimal sets of parameters can be stated as a nonlinear program, NLP, and thus its solution approximated numerically (see details in Additional file 1).

Figure 1 shows the stability function and relative variance function for the Poisson τ-leap, and optimal methods for s = 3 and s = 5 under the constraints |ψ(z)-1| < 0.1, 0.25 and 0.5 and Table 1 summarises the numerical values for these conditions.

Table 1. Stability regions for methods with bounded relative variance and optimal stability

thumbnailFigure 1. Stability and relative variance for the different methods. Stability and relative variance functions for the Poisson τ-leap method (solid line) and RK τ-leap methods with optimal stability regions and bounded relative variance (ψ) with 3 stages (dotted line) and 5 stages (dashed line). Regions fulfilling the bounds on ψ are shown in grey. Square dots correspond to relative variances computed from 106 simulations each. (a), (b): Relative variance bounded by 0.1. (c), (d): Relative variance bounded by 0.25. (e), (f): Relative variance bounded by 0.5.

Efficient methods with bounded variance and extended stability

Runge-Kutta methods with a given stability polynomial R(z) are not unique. This is because the stability polynomial only reflects the application of a Runge-Kutta method to a linear problem. Nonlinear problems require many additional order conditions to be satisfied in order for a method to have a certain order of accuracy. Thus many different methods can have the same stability polynomial. Furthermore, we have already seen that the relative variance ψ does not directly depend on A but on R(z) thus making all methods with the same stability function behave identically in terms of stationary variance for linear problems. In order to distinguish between methods with the same stability function we would have to consider more complicated nonlinear chemistry and this is beyond the scope of this work. However, we have an explicit way of constructing an efficient method that has a given stability polynomial (i.e. to find values for b and A of the Butcher tableau, see details in Additional file 1). Furthermore, the tableaus build in this way are such that βs = 1, βj = 0, j = 0,...,s - 1 and A has all its elements set to zero except those on the first subdiagonal. These Runge-Kutta schemes obtained in this way are very natural, can be regarded as fixed point iterations and allow the following efficient reformulation of (17)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M71','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M71">View MathML</a>

(29)

It is thus clear that these methods are computationally more efficient than the general case as they only require s-1 evaluations of the expected step-change f(·) instead of the s(s - 1)/2 required in the general framework (17). A collection of methods have been implemented in a branch of the ByoDyn package, v.5.0 [39].

Numerical results

Reversible isomerisation

We compare the new Runge-Kutta framework to the Poisson τ-leap to solve three systems of chemical reactions. The first is the reversible isomerisation test problem in (21) for which we have already developed theoretical results. Numerical simulation of the number of molecules for each of the two components in the system was carried out using the different methods discussed in the previous section with k1 = k2 = 10 (z = -20τ) and X(0) = (100, 100)T. We sampled 106 trajectories for each of the methods and for different fixed τ values. Figure 2 shows a comparison between the true probability density function (PDF) and the histograms of X1 obtained from the different methods and some of the values of τ. Note that the Poisson τ-leap method becomes unstable for τ > 0.1 and so does RK τ-leap with three stages for τ > 0.4. Figure 1 shows that the stationary variances obtained by the simulations are in exact accordance with the theoretical values derived in the previous section.

thumbnailFigure 2. Histogram of X1 in the Reversible isomerisation reaction. Histogram of X1 in the Reversible isomerisation reaction (106 samples used) solved by the SSA (grey background), Poisson τ-leap (dashed line), and optimal RK τ-leap methods with bounded relative variance. (a) τ = 0.05 (z = -1), Optimal RK τ-leap s = 3, ϵ = 0.1 (solid line) and s = 5, ϵ = 0.1 ("+" marks), (b) τ = 0.4 (z = -8), Optimal RK τ-leap s = 3, ϵ = 0.5 (solid line) and s = 5, ϵ = 0.1 ("+" marks), Poisson τ-leap is unstable for this time step. (c) τ = 0.6 (z = -12), Optimal RK τ-leap s = 5, ϵ = 0.5 ("+" marks), Poisson τ-leap and Optimal RK τ-leap s = 3 are unstable for this time step.

Schlögl reaction

We also consider Schlögl's autocatalytic reaction system [34,40] to illustrate the accuracy of the presented framework, developed for the linear case, for nonlinear systems. We use here the same set of parameters as Cao et al. [38] for which this system presents a bimodal PDF for the species X in the stationary state. We have also considered that the non-autocatalytic species are buffered (assuming they are constant) hence reducing the system to a scalar problem (see Table 2). We have again performed 106 simulations for each method and τ value. Figure 3 shows histograms computed by the SSA, Poisson τ-leap and the methods with s = 3, 5. Visual inspection of the plots shows a consistent improvement over the original τ-leap method by means of the multistage RK methods developed here. A more precise comparison of the plots is given in Figure 4, which shows the estimated Kullback-Leibler divergences between the exact PDF (PE) and the PDFs of each of these methods (PM), given by:

Table 2. Details of the Schlögl reaction system

thumbnailFigure 3. Histogram of X in the Schlögl reaction. Histogram of X in the Schlögl reaction (106 samples used) solved by the SSA (grey background), Poisson τ-leap (dashed line), and Optimal RK τ-leap methods. (a) τ = 0.4, Optimal RK τ-leap s = 3, ϵ = 0.1 (solid line) and s = 5, ϵ = 0.1 ("+" marks), (b) τ = 0.8, Optimal RK τ-leap s = 3, ϵ = 0.5 (solid line) and s = 5, ϵ = 0.5 ("+" marks).

thumbnailFigure 4. Kullback-Leibler divergence for the Schlögl reaction. Kullback-Leibler divergence between the exact stationary distribution of X in the Schlögl reaction (estimated by 106 samples solved by SSA) and the approximate stationary distributions obtained with the Poisson τ-leap (black), Optimal RK s = 3; ϵ = 0.5 (grey lines) and Optimal RK s = 5, ϵ = 0.5 (white). Bars are shown only for the stable method and τ settings. Asterisks denote methods that have a rate of failure above 10-3.

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M76','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M76">View MathML</a>

(30)

The MAPK cascade

Finally, we have tested the performance of our methods on a larger system of chemical reactions with stiffness due to different reaction time scales and species amounts ranging over several orders of magnitude. For this purpose we considered the Huang and Ferrell model for the mitogen-activated protein kinase (MAPK) cascade [41]. This model is available from the BioModels database [42] and consists of 22 species interacting through 30 reaction channels. The set of parameters used here (see Additional file 1 for details) renders the model stiff and with species amounts ranging from none up to 3·105 molecules. With the chosen initial conditions the system undergoes a transient change and finally settles down into a stationary state at around t = 150 minutes. We have simulated the model using SSA (Gillespie's Direct Method), the Poisson τ-leap and the RK methods presented here. To produce fair comparisons, all methods have been rewritten in ANSI C using the Mersenne twister [43] pseudorandom number generator from the GNU Scientific Library. The GNU C Compiler was used to compile the sources with the -O2 optimisation flag. The algorithms were run on an Intel(R) Core(TM)2 Duo Processor E8500 at 3.16 GHz and 6 MB cache. We have run the system to a final time T = 200. Simulations run with SSA took 61, 841 ± 74 seconds.

We have compared the methods in two distinct situations. First we have run them with the same time step τ = 5·10-5. In this case, the Poisson τ-leap method took 51.7 ± 0.4 seconds while Optimal RK τ-leap methods with s = 3 and s = 5 took 86.1 ± 0.4 seconds and 113.9 ± 0.3 seconds respectively. Hence, at the same time step the RK methods are approximately 66% and 120% slower than the Poisson τ-leap due to the multiple evaluations of the propensity functions per step. However, there is an important difference in the results. The relative variance at the steady state is 1.3 (see Additional file 1) for the Poisson τ-leap while for both RK τ-leap methods with s = 3 and s = 5 (ϵ = 0.1) it is less than 1.04.

Then we have compared these methods when run at their respective maximum time steps such that the relative variance at the stationary state is bounded to 1.1 (estimated from the simulations). The maximum time steps allowed with this constraint were: τ = 2·10-5 for the Poisson τ-leap, τ = 3.5·10-4 for the RK τ-leap (3, 0.1) and τ = 9.5·10-4 for the Optimal RK τ-leap (5, 0.1). With this setting, the runtimes obtained were: 111.9 ± 0.7 seconds for the Poisson τ-leap, 15.7 ± 0.06 seconds for the RK τ-leap (3, 0.1) and 7.8 ± 0.02 seconds for the Optimal RK τ-leap (5, 0.1). Thus, in this case the Poisson τ-leap approximately 7.1 and 14.3 times slower than the RK methods, respectively.

Discussion

Biochemical kinetics typically deals with multiscale problems, in which several scales of time, space and concentrations, simultaneously affect the dynamical behaviour of the system. Thus, the systems biology community is deeply interested in the development of methods that lead to a multiscale view of biochemical systems. As a first step in this workflow, we have presented here a new set of methods that considerably expands the classical τ-leap implementation, from a stability perspective. The importance of the results shown here embraces not only the increase in computational speed for stochastic simulations, a key element for the understanding of the intrinsically noisy biological systems, but more importantly, a way to deal with fast reactions in multiscale settings. The methods developed here have been demonstrated for a first example of stiff system, the classical Schlögl autocatalytic reaction, and can be straightforwardly incorporated into hybrid SSASDE-ODE frameworks.

We see from Table 1 that if we require a bound on the equilibrium variance of 0.1 then the Poisson τ-leap method must take <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/110/mathml/M77','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/110/mathml/M77">View MathML</a> while for the RK methods the bounds on |z| are approximately 4 and 10, respectively with s = 3, 5. This is a very considerable improvement and all the more striking given that the same number of Poisson random variables are simulated per step in all cases.

Initially we had hoped that an approach via Chebyshev methods using ideas from ODEs and SDEs applied to the discrete cases would have been fruitful. It turns out that while such methods have good mean behaviour, the variance behaviour is poor. This is because the variance growth function satisfies (24) and an s-stage Chebyshev method would have s - 1 poles and zeros due to the oscillations in the stability function. Similar issues arise even in the damped forms of the Chebyshev formulation. This means that our optimisation approach is the only way of getting good bounds on ψ(z).

Our results on the nonlinear bimodal Schlögl problem show that the RK methods still behave appropriately even on nonlinear problems. For example, from Figure 3 we see that the Poisson τ-leap method is not very accurate with τ = 0.4 and quite poor in picking up the second peak with τ = 0.8. On the other hand the RK methods match the peak quite well, albeit with a slight shift in that peak. Furthermore, numerical results from the MAPK cascade simulations show that our methods can run an order of magnitude faster than the Poisson τ-leap and still give the same accuracy in the results.

Finally, we note that we could extend our RK methods to allow more than one set of Poisson random variables to be simulated per step. We imagine that this would allow even bigger stepsizes but at the cost of taking more simulation time in that the additional Poisson sampling is expensive. We emphasise that although our analysis of these new methods has been given for unimolecular reactions, the simulations of the nonlinear Schlögl reaction and the MAPK cascade indicate that these methods have a more general applicability and we will consider nonlinear analysis via Taylor series expansions in future work.

Authors' contributions

KB, PR and JVF designed the research. KB and PR developed the algorithms and PR implemented them and performed and analysed the simulations. KB, PR and JVF wrote the manuscript. All authors have read and approved the final version of the manuscript.

Acknowledgements

PR would like to thank Marta Dies for helpful discussions. PR acknowledges Obra Social "la Caixa" for funding through the Graduate Fellowship program. Support from Spanish MCINN grant CTQ2008-00755/BQU and from EC-funded projects BioBridge (FP6-2005-LIFESCIHEALTH-7 037909), QosCosGrid (FP6-2005-IST-5 033883) and VPH (FP7-2007-IST-223920) is highly appreciated. JVF participates in the COM-BIOMED network.

References

  1. Turner TE, Schnell S, Burrage K: Stochastic approaches for modelling in vivo reactions. [http:/ / www.sciencedirect.com/ science/ article/ B73G2-4CS4GV4-1/ 2/ f17f5571a06a80aaaaa53238eed83faf] webcite

    Computational Biology and Chemistry 2004, 28(3):165-178. PubMed Abstract | Publisher Full Text OpenURL

  2. Klipp E, Herwig R, Kowald A, Wierling C, Lehrach H: Systems biology in practice. Wiley-VCH Weinheim; 2005. Publisher Full Text OpenURL

  3. Wilkinson D: Stochastic Modelling for Systems Biology. CRC Press; 2006. OpenURL

  4. Villà J, Warshel A: Energetics and Dynamics of Enzymatic Reactions.

    J Phys Chem B 2001, 105:7887-907. Publisher Full Text OpenURL

  5. McAdams H, Arkin A: Stochastic mechanisms in gene expression.

    Proceedings of the National Academy of Sciences 1997, 94(3):814-819. Publisher Full Text OpenURL

  6. Hasty J, Pradines J, Dolnik M, Collins J: Noise-based switches and amplifiers for gene expression.

    Proceedings of the National Academy of Sciences 2000, 97(5):2075-2080. Publisher Full Text OpenURL

  7. Thattai M, van Oudenaarden A: Intrinsic noise in gene regulatory networks.

    Proceedings of the National Academy of Sciences 2001, 151588598. OpenURL

  8. Ozbudak E, Thattai M, Kurtser I, Grossman A, van Oudenaarden A: Regulation of noise in the expression of a single gene.

    Nature Genetics 2002, 31:69-73. PubMed Abstract | Publisher Full Text OpenURL

  9. Isaacs F, Hasty J, Cantor C, Collins J: Prediction and measurement of an autoregulatory genetic module.

    Proceedings of the National Academy of Sciences 2003, 100(13):7714-7719. Publisher Full Text OpenURL

  10. Thattai M, van Oudenaarden A: Stochastic Gene Expression in Fluctuating Environments. [http://www.genetics.org/cgi/content/abstract/167/1/523] webcite

    Genetics 2004, 167:523-530. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  11. Tian T, Burrage K: Bistability and switching in the lysis/lysogeny genetic regulatory network of bacteriophage [lambda]. [http:/ / www.sciencedirect.com/ science/ article/ B6WMD-4B76G7M-1/ 2/ 64146757d7c19ee4acd0247e0d997cb5] webcite

    Journal of Theoretical Biology 2004, 227(2):229-237. PubMed Abstract | Publisher Full Text OpenURL

  12. Bratsun D, Volfson D, Tsimring LS, Hasty J: Delay-induced stochastic oscillations in gene regulation. [http://www.pnas.org/content/102/41/14593.abstract] webcite

    Proceedings of the National Academy of Sciences of the United States of America 2005, 102(41):14593-14598. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  13. Kaern M, Elston T, Blake W, Collins J, et al.: Stochasticity in gene expression: from theories to phenotypes.

    Nat Rev Genet 2005, 6(6):451-464. PubMed Abstract | Publisher Full Text OpenURL

  14. Komili S, Silver P: Coupling and coordination in gene expression processes: a systems biology view.

    Nature Reviews Genetics 2008, 9:38. PubMed Abstract | Publisher Full Text OpenURL

  15. Halley J, Winkler D, Burden F: Toward a Rosetta stone for the stem cell genome: Stochastic gene expression, network architecture, and external influences.

    Stem Cell Research 2008, 1(3):157-168. PubMed Abstract | Publisher Full Text OpenURL

  16. Kurtz TG: The Relationship between Stochastic and Deterministic Models for Chemical Reactions. [http://link.aip.org/link/?JCP/57/2976/1] webcite

    The Journal of Chemical Physics 1972, 57(7):2976-2978. Publisher Full Text OpenURL

  17. Gillespie DT: Exact stochastic simulation of coupled chemical reactions. [http://pubs.acs.org/doi/abs/10.1021/j100540a008] webcite

    The Journal of Physical Chemistry 1977, 81(25):2340-2361. Publisher Full Text OpenURL

  18. Gillespie DT: Stochastic simulation of chemical kinetics. [http://www.ncbi.nlm.nih.gov/pubmed/17037977] webcite

    Annual review of physical chemistry 2007, 58:35-55.

    [10.1146/annurev. physchem.58.032806.104637]

    PubMed Abstract | Publisher Full Text OpenURL

  19. Gibson MA, Bruck J: Efficient Exact Stochastic Simulation of Chemical Systems with Many Species and Many Channels. [http://pubs.acs.org/doi/abs/10.1021/jp993732q] webcite

    The Journal of Physical Chemistry A 2000, 104(9):1876-1889. Publisher Full Text OpenURL

  20. Gillespie DT: Approximate accelerated stochastic simulation of chemically reacting systems. [http://link.aip.org/link/?JCP/115/1716/1] webcite

    The Journal of Chemical Physics 2001, 115(4):1716-1733. Publisher Full Text OpenURL

  21. Cao Y, Gillespie DT, Petzold LR: Efficient step size selection for the tau-leaping simulation method. [http:/ / scitation.aip.org/ getabs/ servlet/ GetabsServlet?prog=normaln\&id=JCPS A6000124000004044109000001\&idtype= cvips\&gifs=yes] webcite

    The Journal of Chemical Physics 2006., 124(4) PubMed Abstract | Publisher Full Text OpenURL

  22. Tian T, Burrage K: Binomial leap methods for simulating stochastic chemical kinetics.

    The Journal of Chemical Physics 2004, 121(21):10356-10364. PubMed Abstract | Publisher Full Text OpenURL

  23. Chatterjee A, Vlachos DG, Katsoulakis MA: Binomial distribution based tau-leap accelerated stochastic simulation. [http://link.aip.org/link/?JCP/122/024112/1] webcite

    The Journal of Chemical Physics 2005, 122(2):024112. PubMed Abstract | Publisher Full Text OpenURL

  24. Auger A, Chatelain P, Koumoutsakos P: R-leaping: Accelerating the stochastic simulation algorithm by reaction leaps. [http://link.aip.org/link/?JCP/125/084103/1] webcite

    The Journal of Chemical Physics 2006, 125(8):084103. PubMed Abstract | Publisher Full Text OpenURL

  25. Monk NA: Oscillatory expression of Hes1, p53, and NF-kappaB driven by transcriptional time delays. [http://view.ncbi.nlm.nih.gov/pubmed/12932324] webcite

    Curr Biol 2003, 13(16):1409-1413. PubMed Abstract | Publisher Full Text OpenURL

  26. Yildirim N, Mackey MC: Feedback Regulation in the Lactose Operon: A Mathematical Modeling Study and Comparison with Experimental Data. [http://www.biophysj.org/cgi/content/abstract/84/5/2841] webcite

    Biophys J 2003, 84(5):2841-2851. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  27. Barrio M, Burrage K, Leier A, Tian T: Oscillatory Regulation of Hes1: Discrete Stochastic Delay Modelling and Simulation.

    PLoS Comput Biol 2006, 2(9):e117. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  28. Leier A, Marquez-Lago TT, Burrage K: Generalized binomial tau-leap method for biochemical kinetics incorporating both delay and intrinsic noise. [http://link.aip.org/link/?JCP/128/205107/1] webcite

    The Journal of Chemical Physics 2008, 128(20):205107. PubMed Abstract | Publisher Full Text OpenURL

  29. Anderson DF: A modified next reaction method for simulating chemical systems with time dependent propensities and delays. [http://link.aip.org/link/?JCP/127/214107/1] webcite

    The Journal of Chemical Physics 2007, 127(21):214107. PubMed Abstract | Publisher Full Text OpenURL

  30. Hairer E, Norsett SP, Wanner G: Solving Ordinary, Differential Equations II. Stiff and Differential-Algebraic Problems. Volume 2. Springer-Verlag, Second Revised Editio edition; 2002.

    [Index]

    OpenURL

  31. Abdulle A, Medovikov A: Second order Chebyshev methods based on orthogonal polynomials.

    Numerische Mathematik 2001, 90:1-18. Publisher Full Text OpenURL

  32. Abdulle A, Cirilli S: S-ROCK: Chebyshev Methods for Stiff Stochastic Differential Equations. [http://portal.acm.org/citation.cfm?id=1350482&jmp=cit&coll=GUIDE&dl=] webcite

    SIAM J Sci Comput 2008, 997-1014. Publisher Full Text OpenURL

  33. Hernandez D, Spigler R: Convergence and stability of implicit runge-kutta methods for systems with multiplicative noise.

    BIT Numerical Mathematics 1993, 33(4):654-669. Publisher Full Text OpenURL

  34. Schlögl F: Chemical reaction models for non-equilibrium phase transitions.

    Zeitschrift für Physik A Hadrons and Nuclei 1972, 253(2):147-161. OpenURL

  35. Øksendal B: [http:/ / www.amazon.ca/ exec/ obidos/ redirect?tag=citeulike09-20\&amp;pa th=ASIN/ 3540047581] webcite

    Stochastic Differential Equations: An Introduction with Applications (Universitext). Springer; 2005. OpenURL

  36. Burrage K, Burrage PM: High strong order explicit Runge-Kutta methods for stochastic ordinary differential equations.

    Applied Numer Maths 1996, 22:81-101. Publisher Full Text OpenURL

  37. Li T: Analysis of Explicit Tau-Leaping Schemes for Simulating Chemically Reacting Systems. [http://link.aip.org/link/?MMS/6/417/1] webcite

    Multiscale Modeling and Simulation 2007, 6(2):417-436. Publisher Full Text OpenURL

  38. Cao Y, Petzold LR, Rathinam M, Gillespie DT: The numerical stability of leaping methods for stochastic simulation of chemically reacting systems.

    J Chem Phys 2004, 121(24):12169-12178. PubMed Abstract | Publisher Full Text OpenURL

  39. de Lomana ALG, Gómez-Garrido A, Sportouch D, Villà-Freixa J: Optimal Experimental Design in the Modelling of Pattern Formation. [http://www.springerlink.com/content/kk7774170666m254/] webcite

    LNCS 2008, 5101:610-619. OpenURL

  40. van Kampen NG: Stochastic Processes in Physics and Chemistry. Elsevier; 2007. OpenURL

  41. Huang CY, Ferrell JE: Ultrasensitivity in the mitogen-activated protein kinase cascade. [http://www.pnas.org/content/93/19/10078.abstract] webcite

    Proceedings of the National Academy of Sciences of the United States of America 1996, 93(19):10078-10083. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  42. Novère NL, Bornstein BJ, Broicher A, Courtot M, Donizelli M, Dharuri H, Li L, Sauro HM, Schilstra MJ, Shapiro BE, Snoep JL, Hucka M: BioModels Database: a free, centralized database of curated, published, quantitative kinetic models of biochemical and cellular systems.

    Nucleic Acids Research 2006, (34 Database):689-691. Publisher Full Text OpenURL

  43. Matsumoto M, Nishimura T: Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator.

    ACM Transactions on Modeling and Computer Simulation 1998, 8:3-3. Publisher Full Text OpenURL