Biochemical fluctuations, optimisation and the linear noise approximation
© Pahle et al; licensee BioMed Central Ltd. 2012
Received: 10 August 2011
Accepted: 1 March 2012
Published: 17 July 2012
Stochastic fluctuations in molecular numbers have been in many cases shown to be crucial for the understanding of biochemical systems. However, the systematic study of these fluctuations is severely hindered by the high computational demand of stochastic simulation algorithms. This is particularly problematic when, as is often the case, some or many model parameters are not well known. Here, we propose a solution to this problem, namely a combination of the linear noise approximation with optimisation methods. The linear noise approximation is used to efficiently estimate the covariances of particle numbers in the system. Combining it with optimisation methods in a closed-loop to find extrema of covariances within a possibly high-dimensional parameter space allows us to answer various questions. Examples are, what is the lowest amplitude of stochastic fluctuations possible within given parameter ranges? Or, which specific changes of parameter values lead to the increase of the correlation between certain chemical species? Unlike stochastic simulation methods, this has no requirement for small numbers of molecules and thus can be applied to cases where stochastic simulation is prohibitive.
We implemented our strategy in the software COPASI and show its applicability on two different models of mitogen-activated kinases (MAPK) signalling -- one generic model of extracellular signal-regulated kinases (ERK) and one model of signalling via p38 MAPK. Using our method we were able to quickly find local maxima of covariances between particle numbers in the ERK model depending on the activities of phospho-MKKK and its corresponding phosphatase. With the p38 MAPK model our method was able to efficiently find conditions under which the coefficient of variation of the output of the signalling system, namely the particle number of Hsp27, could be minimised. We also investigated correlations between the two parallel signalling branches (MKK3 and MKK6) in this model.
Our strategy is a practical method for the efficient investigation of fluctuations in biochemical models even when some or many of the model parameters have not yet been fully characterised.
Random fluctuations in discrete molecular numbers can have significant impact, both detrimental and constructive, on the functioning of biochemical systems [1, 2]. Systems that contain only relatively small numbers of particles of a certain chemical species, such as in signal transduction or gene expression, are particularly prone to this intrinsic noise. Here, the underlying discreteness of the system and stochastic timing of reactive events can lead to fluctuations in species abundances of high relative amplitude. Even when particle numbers are high, stochastic effects can significantly affect the dynamic behaviour of certain biochemical networks .
Biochemical systems have evolved to be robust against molecular fluctuations by attenuation, or even to exploit them (see [4, 5] for examples). Therefore, these fluctuations should be considered whenever quantitative and dynamic models are devised to describe biochemical systems. Different mathematical formalisms have been developed to allow stochastic modelling and to explicitly take into account random fluctuations. Such systems are usually modelled by a continuous-time Markov process which follows the chemical master equation. The chemical master equation describes the time evolution of the system state probability distribution, i.e. how probable it is that a chemical species in the system will have specific particle numbers at a specific point in time. Both analytic and numerical solutions of this chemical master equation are difficult to obtain for most biologically relevant systems. Even though there exist methods to numerically solve the master equation  these are only feasible for relatively simple systems. A popular substitute is to apply Gillespie's stochastic simulation algorithm  to calculate single trajectories of the system's dynamics. By calculating very many of such (random) instances one can then approximate the trajectory of the probability density function of each chemical species and calculate relevant time-dependent statistics, such as the mean value or covariances. However, the stochastic simulation of single trajectories alone can be computationally demanding. The calculation of very many of them quickly becomes impracticable even when accelerated approximate stochastic simulation methods  are employed.
For a quick characterisation of the fluctuations in a biochemical system there exists an alternative, namely the linear noise approximation (LNA; see, e.g., [9–11]). This approximate method is based on van Kampen's system-size expansion of the chemical master equation [12–14]. The LNA estimates the variances of the species abundances and the covariances between them. Even though, theoretically, the LNA is only locally valid in the vicinity of macroscopic steady states or other system trajectories, in practical terms, it often gives good results even when the behaviour of the stochastic model and the behaviour of the corresponding deterministic model are quite different . The LNA is particularly interesting because it is independent of computationally demanding stochastic simulations but, instead, only uses information about the stoichiometries in the system and the macroscopic reaction rates -- therefore it can be calculated very quickly. Other approaches have also been proposed for the estimation of steady state noise. For instance, in , analytical estimates of the fluctuations are found using error growth techniques. These are based on ideas from nonlinear dynamics and do not begin from a master equation. This is in contrast to the work presented here, where the molecular basis of the model is central, and where the nature of the fluctuations can be explicitly calculated. There are, of course, many studies of fluctuations in biochemical systems. For instance, in  the authors use data from time series to infer the values of the model parameters. This is in some sense the converse of our approach.
Often, in practice, one or more of the parameters of a model, such as reaction rates or initial concentrations, cannot be exactly determined. For instance, such parameters might only be known to lie within a certain range or nothing might be known about them at all. This uncertainty about parameters can translate into uncertainty about the system behaviour when it has high sensitivity towards those parameters. This is also true for molecular fluctuations in the system since their expected amplitude and other properties depend on parameter values. If only one or two parameters are unknown it is possible to exhaustively scan this parameter space using a regular grid or other techniques to probe how the model is affected by variations in values of those parameters. However, this approach is not feasible if the number of unknown parameters is large since the hyper-volume of the parameter search space increases exponentially with the number of uncertain parameters, and consequently so does the computational time.
In this article we introduce a different strategy to study random fluctuations in biochemical models with parameters that are not well characterised. Our approach combines the LNA with optimisation methods to search the unknown parameter space for parameter values that lead to extrema in covariance estimates. This can dramatically reduce the required computation time compared to exhaustive searches with stochastic simulations, thereby permitting types of studies of stochastic fluctuations that were not possible before. We will show a relevant biological example of a search for conditions that minimise the noise in the output of a p38 MAPK signalling system. Scanning the parameter space and using stochastic simulation is clearly impossible here because this would take more than 2.4 · 1017 years. Our method, in contrast, was able to find these conditions in 25 min. Therefore, the strategy we are proposing makes it possible to gain biological insight about the noise structure of relevant biological systems even if these systems are big and the parameters are not well defined.
Global optimisation methods have been shown to be effective in finding good extrema estimates of dynamic properties of biochemical network models even in high-dimensional search spaces [18, 19]. The strategy proposed here is similar to an earlier one successfully applied to the search for extreme values of sensitivities .
The application of this strategy passes through a closed loop containing the automatic calculation of a steady state, the LNA method and one optimisation algorithm; alternatively the method is also appropriate to use with parameter scanning or sampling algorithms instead of the optimisation. We implemented this strategy in the software COPASI [21, 22], which already contains optimisation, scanning and sampling algorithms. We demonstrate the application of this new strategy on two different models of mitogen-activated kinase (MAPK) signalling pathways, namely a model of extracellular signal-regulated kinases (ERK) by Kholodenko  and a model of p38 MAPK by Hendriks et al. .
Implementation of the method in COPASI
Our implementation allows arbitrary objective functions to be optimised. For instance, LNA estimates of covariances of different chemical species, as well as other model observables, can be combined into a complex objective function. This allows the calculation of various quantities of interest, for instance, Fano factors  or coefficients of variation (CV), as shown below. In terms of parameter search, our implementation can use a large variety of numerical optimisation algorithms, both local and global, that are accessible in COPASI -- gradient-based, particle swarm , simulated annealing, evolutionary algorithms and others . This is particularly important since the performance of global optimisation algorithms has been shown to be problem-dependent, and no single one is guaranteed to converge to a global optimum for all problems .
Application of the method on MAPK signalling systems
Signalling through mitogen-activated protein kinases (MAPK) is involved in a broad range of cellular processes, such as proliferation, differentiation, stress responses and apoptosis. Therefore it is also implicated in a variety of diseases like cancer, stroke or diabetes . As such, it has been the object of a number of computational modelling studies that helped elucidate dynamic properties of the system, such as amplification of signals, noise reduction or switching behaviour .
There exist different specific MAPK signalling pathways with different functions, for example ERK1/2, p38 or JNK, with different topologies and characteristics. However, in most cases the basic structure is that of a three-tier cascade. Here, the MAPKs on the output level, such as ERK1/2 or p38, phosphorylate transcription factors or other proteins to trigger specific cellular responses. The MAPKs are, in turn, activated via phosphorylation by other protein kinases, so-called MAP2K (or MKK) that are themselves activated by MAP3K (or MKKK) further upstream.
Fluctuations in a model of ultrasensitivity in ERK MAP kinase signalling
It is interesting to see how the magnitude of fluctuations changes with the reaction parameters. As an example, we used our LNA implementation in COPASI in combination with a parameter scan to investigate how changes in the reaction parameter v2 affect the variance of MKKK (MAPK kinase kinase). Values of v2 were scanned within a certain range and the LNA automatically calculated for each value of v2. In the model, this parameter corresponds to the V max of phospho-MKKK dephosphorylation and so refers to the activity of MKKK-phosphatase.
Fluctuations in a model of p38 MAPK signalling
The so-called p38 mitogen-activated protein kinases (p38 MAPK) are responsive to proinflammatory cytokines and stress factors . One prominent signal are lipopolysaccharides (LPS), which are components in the cell wall of bacteria. Their presence indicates a bacterial infection and triggers a strong immune response in animals. The MAPK of this pathway, p38, can, inter alia, activate MAP kinase-activated protein kinase 2 (MK2). One substrate of MK2 is the heat shock protein 27 (Hsp27) and the concentration of the active/phosphorylated form of Hsp27 is regularly used to estimate the activity of the p38 MAPK signalling pathway. The level of Hsp27 will also represent the main signalling output in the model.
The model we use for this study was developed in Hendriks et al. 2008 . Its structure is shown in the additional file 1. The original model included the rapid inactivation of a (TAK1:TAB1:TAB2) complex. This was represented by a degradation reaction which, after an initial peak, led to an abrogation of p38 MAPK activity. For this study we removed this degradation reaction which allows the system to reach a steady state of sustained p38 MAPK signalling depending on the amount of LPS. We also reformulated the model in such a way that it no longer contained three compartments (medium, cytosol and nucleus) but, instead, uses a single reference volume for all species, including the nuclear ones, corresponding to the volume of the whole cell. This was needed because the current implementation of the LNA can only handle models with one compartment. In  the model was fitted to experimental measurements, and in the following we will use the set of parameters which showed the best fit.
As mentioned above, random fluctuations in signalling systems are particularly interesting to study, since here copy numbers of the different species are often low. For instance, MKK3 and MKK6 are typically present in the order of only ten thousand particles per cell. This could lead to pronounced fluctuations which hamper reliable information transfer through this signalling pathway. But perhaps there are conditions (parameter values) for which these fluctuations are minimised, which is what we want to investigate.
By contrast, phospho-Hsp27, the endpoint of the modelled signalling pathway, shows a decrease in its variance with increasing stimulation (Figure 5, panel B).
An interesting property of the p38 MAPK pathway is the existence of two parallel signalling branches, through MKK3 and MKK6, that both can phosphorylate p38 MAPK. Therefore, we were interested in whether fluctuations in the MKK3 branch correlate with fluctuations in the MKK6 branch. First, we scanned the estimated covariance of phospho-MKK3 and phospho-MKK6 over a range of stimulus strengths. We found that the fluctuations in the two branches seem to be mostly uncorrelated (the LNA actually estimates a very weak anti-correlation for higher initial concentrations of LPS, data not shown), an indication that the largest part of the fluctuations does not originate from the common upstream part of the two branches but rather from within the branches themselves.
We now wanted to investigate how the parameters in the system influence this anti-/correlation. Therefore, we searched for extreme values of the LNA-estimated correlation coefficient of phospho-MKK3 and phospho-MKK6 within a fixed, but large, range of all parameter values.
We therefore ran the LNA in combination with the particle swarm optimisation algorithm of COPASI, using the correlation coefficient as the objective function for maximisation. In addition we set constraints on the number of steady state particle numbers in the system. Both phospho-MKK3 and phospho-MKK6 particle numbers were allowed to change only 4-fold, i.e. within 50% - 200% of their original values. All other species' particle numbers were allowed to change 100-fold, i.e. within 10% - 1000% of their original values. The reasons for this were, firstly, that we did not primarily want to change the steady state of the system but rather only wanted to affect the fluctuations around the steady state. Secondly, if particle numbers are not constrained the optimisation often converges towards degenerate cases where one or both of the steady state particle numbers are very close to zero, i.e. the lower limit -- a situation where the LNA estimation can have large errors due to its assumption of Gaussian fluctuations.
We used a particle swarm optimisation  method with a swarm size of 50. The parameters to vary were all 29 reaction rates of the first 20 reactions in the model (see ), which includes all receptor (complex)-related reactions, both MKK3 and MKK6 branches, and the phosphorylation and dephosphorylation of p38. The parameters were allowed to change 100-fold, e.g. within 10% - 1000% of their original values. With these settings our method was able to find conditions where the estimated correlation between phospho-MKK3 and phospho-MKK6 was larger than 0.95 with a computation time of roughly 70 min.
Optimisation of the coefficient of variation of phospho-Hsp27 particle numbers
complex + MKK6 ↔ complex_MKK6
complex_MKK6 → complex + MKK6P
MKK6_phosphatase + MKK6P ↔ Ppase_MKK6P
Ppase_MKK6P → MKK6_phosphatase + MKK6
complex + MKK3 ↔ complex_MKK3
complex_MKK3 → complex + MKK3P
MKK3_phosphatase + MKK3P ↔ Ppase_MKK3P
Ppase_MKK3P → MKK3_phosphatase + MKK3
MKK6P + p38 ↔ MKK6P_p38
MKK6P_p38 → MKK6P + p38P
MKK3P + p38 ↔ MKK3P_p38
MKK3P_p38 → MKK3P + p38P
p38_phosphatase + p38P ↔ Ppase_38P
Ppase_p38P → p38_phosphatase + p38
We would like to note here that a (naive) comprehensive search for optima using a regular grid approach and stochastic simulations of the system in this particular case would have taken a prohibitively long computation time. Assuming that, within the given limits, we only look at ten different values per parameter we would have 10<no.parameters>= 1021 sample points. For each point we would need to carry out a stochastic simulation that, including the calculation to allow the system to settle down to a steady state, takes approximately 7700 s on a typical desktop computer (for a simulated time of 10000 s). Neglecting the time needed to calculate the actual statistics on the simulated time series this would lead to a computation time of more than 1021 · 7700 s ≈ 2.4 · 1017 years. And this would only explore ten values of each parameter (i.e. it would be at a low resolution.) In contrast our method, using the linear noise approximation in combination with numerical optimisation, took 25 min to converge. This clearly shows the utility of the method we propose here: it makes tractable to calculate many phenomena that otherwise would be computationally prohibitive. Finally, although an approximation had to be adopted, it is typically so good that this has very little impact on the accuracy of the method.
Our contribution with this work is two-fold. First, we implemented the linear noise approximation in the freely available software COPASI, and thus made it accessible to a large group of users. Secondly, we showed how the LNA in combination with multi-dimensional parameter scans or with global numerical optimisation methods is appropriate to quickly characterise the influence of parameters on intrinsic fluctuations in biochemical models even when there is considerable uncertainty about a number of parameters. We showed, with realistic biochemical signalling models, that using this approach one is able to explore parameter space such that conditions can be found for which there is minimal, or maximal, noise. It is also possible to search for conditions where specific model variables are highly (or poorly) correlated. This new method thus provides a new and important way to explore the universe of behaviours displayed by models. Given the importance of noise and fluctuations in intracellular biochemistry, this method is therefore of great value for the study of those systems.
In the recent article by Komorowski et al.  a related method is proposed. There, the linear noise approximation is used to calculate Fisher information matrices for stochastic models, primarily to inform experimental design, e.g. by examining the information content of different experimental samples. Our approach, on the other hand, focuses on exploring the model independently of any physical measurements. Therefore, the two approaches are complementary.
In certain cases, however, care should be taken when using the LNA. This is due to the assumption that the fluctuations are Gaussian in nature. Problems can arise if the system is close to a boundary. For example, if the number of molecules for a particular species is very close to zero the probability distribution for the fluctuations becomes 'squashed' (which the LNA does not take into account), to satisfy the requirement that the probability to have a negative number of molecules present is zero. Boundaries can also arise due to conservation relations, which are discussed in the Methods section, as these add constraints to the system. When using the LNA in combination with one of the optimisation algorithms in COPASI, such systems near boundaries are sometimes found, especially when the user wishes to minimise a covariance, as we found when studying the p38 MAPK model. This is because the fluctuations can be very small when the system is close to a boundary, which can give the impression that the fluctuations of two different species are uncorrelated, which may not be the case away from the boundary. In these cases, adding constraints to the particle numbers (as we did when studying the p38 MAPK model) helps to keep the system away from these states. The current implementation of the LNA in COPASI is only able to consider models in which the reactions all occur within one compartment. As many biochemical models involve multiple compartments we hope to extend our work, so that in future it will be possible to use the LNA to study a wider range of models.
Michaelis-Menten reaction mechanism
v1 = k1
S + E → SE
v2 = k2 · [S] · [E]
SE → S + E
v3 = k3 · [SE]
SE → P + E
v4 = k4 · [SE]
v5 = k5 · [P]
A substrate, S, is converted to a product, P, via an enzyme E. The substrate and enzyme form a complex, SE. A constant flux of S molecules is supplied to the system and P molecules are able to leave the system. In our notation above, Y1 = S, Y2 = E, Y3 = P and Y4 = SE. Also, r11 = 1, r21 = 0 and so on. The total number of enzyme molecules, i.e. the number of free enzyme molecules plus the enzyme molecules bound in a complex, is fixed. If the number of SE molecules decreases by one, the number of E molecules increases by one. This conservation relation means that there is only one independent variable here, not two. In general, if the system contains Λ conservation relations, then the dimension of the system can be reduced from to . It is necessary to reduce the size of the system in this way to facilitate the linear algebra to be done later. In the Michaelis-Menten system above, and Λ = 1, so K = 3.
(The ODEs for the species that have been eliminated can be found by using the conservation relations.) However, the large system size limit is inappropriate for many systems of interest, in particular when the molecular populations are low (and the volume is small, as in most cells), then the discrete nature of the molecules has important consequences. In these cases a stochastic description is required.
In all the investigations we will carry out in this paper, we will be interested in fluctuations about the stationary state. In terms of the deterministic dynamics Eq. (2), the solution x(t) will be replaced by its fixed point value x*, and so the A and B matrices will be independent of time.
The equation can be solved for C numerically by employing the Bartels-Stewart algorithm . Here, A is transformed to lower real Schur form and A T is transformed to real Schur form. This allows elements of the transformed matrix C to be solved for successively. The solution of C is found by reversing the original transformation. It is important that the conservation relations are used to reduce the dimension of the system before the Bartels-Stewart algorithm is applied. The matrices for the 'unreduced' system will contain linearly dependent columns or, equivalently, zero eigenvalues. When this is the case, the solution to the equation is no longer unique .
Therefore our implementation of the LNA first automatically determines existing conservation relations (also known as conserved moieties) and reduces the system from to K independent chemical species. Then the LNA is applied to the reduced system and the corresponding covariance matrix is calculated. In the last step the covariance matrix for the full system is recovered as follows.
where the c j and α jk are constants.
with Cred. the K × K covariance matrix of the reduced system.
Covariance matrix C for the Michaelis-Menten reaction system.
As just mentioned, we implemented the LNA described above in the software COPASI [21, 22]. COPASI is a widely used software for the analysis and simulation of biochemical networks. It lets the users access sophisticated mathematical methods, such as deterministic, stochastic and hybrid simulation, metabolic control analysis, sensitivity analysis, optimisation and parameter fitting, to study their models. COPASI also allows closed-loop applications of parameter scanning, sampling and optimisation with one of the other analyses, for example sensitivity analysis or the linear noise approximation. Models can be conveniently imported and exported using the Systems Biology Markup Language (SBML) . COPASI is an open source software and is freely available under the Artistic license version 2.0 at .
Briefly, our LNA implementation in COPASI first detects dependent species (conservation relations) and carries out the corresponding reduction of the system, if needed. Then an automatic search for a steady state of the system is started. If a steady state has been found the Lyapunov equation Eq. (10) for the reduced system is solved at this steady state using the Bartels-Stewart algorithm. Finally, the covariance matrix for the full system is recovered as described above.
In addition, before the LNA is carried out COPASI automatically checks the model according to a number of criteria that preclude a direct calculation of the LNA. For instance, if there are reversible reactions present in the model COPASI will notify the user that all reversible reactions have to be split into irreversible reactions before the LNA can be applied. There exists a tool in COPASI which can do this in an automated way for a large class of models.
Optimisation is a general modelling tool with a wide application to the solution of diverse problems. Essentially, if something can be specified as a maximum or minimum of some function, optimisation will be the way to solve such a problem. In biochemical network modelling the most common application is parameter estimation; another one is the design of genetically engineered pathways (commonly known as metabolic engineering) where one seeks to maximise a flux, titre or a yield of a biotransformation . Other popular applications are those where a specific parameter set is sought that produces a desired behaviour of the model. This is the basis of a method that tabulates all maxima and minima of each parameter sensitivity towards a specific model variable , as a means of approaching global sensitivity analysis. All of these applications require that a simulator be integrated with an optimisation algorithm in a closed loop.
There are many different numeric algorithms for searching minima (or maxima) of functions: the traditional gradient-based methods, direct search that use geometric heuristics, population-based algorithms like evolutionary algorithms and particle swarm [27, 28], and stochastic searches like simulated annealing. Often problems are complex in that the objective function is not convex and can have several local minima yet one seeks one of the global minima. For such problems it is necessary to employ algorithms that are not trapped in local minima, as are the gradient-based algorithms. Empirical evidence shows that population-based and stochastic search algorithms are commonly the most efficient at finding global minima. In our experience with biochemical networks this is usually achieved by evolutionary algorithms  or the particle swarm algorithm. All of these algorithms are available in our COPASI implementation.
JP and PM thank the UK's BBSRC (grant BB/F018398/1) and BBSRC/EPSRC (grant BB/C008219/1), the USA's NIH (grant GM080219), the German Federal Ministry of Education and Research (BMBF), the Klaus Tschira Foundation, and the EU FP7 project UniCellSys, grant no. 201142. JDC thanks EPSRC for the award of a PhD studentship. This is a contribution from the Manchester Centre for Integrative Systems Biology.
- McAdams HH, Arkin A: It's a noisy business! Genetic regulation at the nanomolar scale. Trends Genet 1999,15(2):65-69. 10.1016/S0168-9525(98)01659-XView ArticleGoogle Scholar
- Elowitz MB, Levine AJ, Siggia ED, Swain PS: Stochastic gene expression in a single cell. Science 2002, 297: 1183-1186. 10.1126/science.1070919View ArticleGoogle Scholar
- Kummer U, Krajnc B, Pahle J, Green AK, Dixon CJ, Marhl M: Transition from stochastic to deterministic behavior in calcium oscillations. Biophys J 2005,89(3):1603-1611. 10.1529/biophysj.104.057216View ArticleGoogle Scholar
- Rao CV, Wolf DW, Arkin AP: Control, exploitation and tolerance of intracellular noise. Nature 2002, 420: 231-237. 10.1038/nature01258View ArticleGoogle Scholar
- Eldar A, Elowitz MB: Functional roles for noise in genetic circuits. Nature 2010, 467: 167-173. 10.1038/nature09326View ArticleGoogle Scholar
- Munsky B, Khammash M: The finite state projection approach for the analysis of stochastic noise in gene networks. IEEE Trans Automat Contr 2008, 53: 201-214.View ArticleGoogle Scholar
- Gillespie DT: A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J Comp Phys 1976, 22: 403-434. 10.1016/0021-9991(76)90041-3View ArticleGoogle Scholar
- Pahle J: Biochemical simulations: stochastic, approximate stochastic and hybrid approaches. Brief Bioinform 2009, 10: 53-64.View ArticleGoogle Scholar
- Elf J, Ehrenberg M: Fast evaluation of fluctuations in biochemical networks with the linear noise approximation. Genome Res 2003, 13: 2475-2484. 10.1101/gr.1196503View ArticleGoogle Scholar
- Tomioka R, Kimura H, Kobayashi TJ, Aihara K: Multivariate analysis of noise in genetic regulatory networks. J Theor Biol 2004,229(4):501-521. 10.1016/j.jtbi.2004.04.034View ArticleGoogle Scholar
- McKane AJ, Nagy JD, Newman TJ, Stefani MO: Amplified biochemical oscillations in cellular systems. J Stat Phys 2007,128(1/2):165-191.View ArticleGoogle Scholar
- van Kampen N: A power series expansion of the master equation. Can J Phys 1961,39(4):551-567. 10.1139/p61-056View ArticleGoogle Scholar
- van Kampen N: The expansion of the master equation. Adv Chem Phys 1976, 34: 245-309.Google Scholar
- van Kampen N: Stochastic Processes in Physics and Chemistry. Amsterdam: Elsevier; 1981.Google Scholar
- Hayot F, Jayaprakash C: The linear noise approximation for molecular fluctuations within cells. Phys Biol 2004,1(4):205-210. 10.1088/1478-3967/1/4/002View ArticleGoogle Scholar
- Orrell D, Ramsey S, de Atauri P, Bolouri H: A method for estimating stochastic noise in large genetic regulatory networks. Bioinformatics 2005,21(2):208-217. 10.1093/bioinformatics/bth479View ArticleGoogle Scholar
- Munsky B, Trinh B, Khammash M: Listening to the noise: random fluctuations reveal gene network parameters. Mol Syst Biol 2009, 5: 318.View ArticleGoogle Scholar
- Mendes P, Kell DB: Non-linear optimization of biochemical pathways: applications to metabolic engineering and parameter estimation. Bioinformatics 1998,14(10):869-883. 10.1093/bioinformatics/14.10.869View ArticleGoogle Scholar
- Moles CG, Mendes P, Banga JR: Parameter estimation in biochemical pathways: a comparison of global optimization methods. Genome Res 2003,13(11):2467-2474. 10.1101/gr.1262503View ArticleGoogle Scholar
- Sahle S, Mendes P, Hoops S, Kummer U: A new strategy for assessing sensitivities in biochemical models. Phil Trans R Soc A 2008,366(1880):3619-3631. 10.1098/rsta.2008.0108View ArticleGoogle Scholar
- Hoops S, Sahle S, Gauges R, Lee C, Pahle J, Simus N, Singhal M, Xu L, Mendes P, Kummer U: COPASI - a COmplex PAthway SImulator. Bioinformatics 2006,22(24):3067-3074. 10.1093/bioinformatics/btl485View ArticleGoogle Scholar
- Kholodenko BN: Negative feedback and ultrasensitivity can bring about oscillations in the mitogen-activated protein kinase cascades. Eur J Biochem 2000, 267: 1583-1588. 10.1046/j.1432-1327.2000.01197.xView ArticleGoogle Scholar
- Hendriks BS, Hua F, Chabot JR: Analysis of mechanistic pathway models in drug discovery: p38 pathway. Biotechnol Prog 2008, 24: 96-109. 10.1021/bp070084gView ArticleGoogle Scholar
- Hucka M, Finney A, Sauro HM, Bolouri H, Doyle JC, Kitano H, the rest of the SBML Forum, Arkin AP, Bornstein BJ, Bray D, Cornish-Bowden A, Cuellar AA, Dronov S, Gilles ED, Ginkel M, Gor V, Goryanin II, Hedley WJ, Hodgman TC, Hofmeyr JH, Hunter PJ, Juty NS, Kasberger JL, Kremling A, Kummer U, Le Novère N, Loew LM, Lucio D, Mendes P, Minch E, Mjolsness ED, Nakayama Y, Nelson MR, Nielsen PF, Sakurada T, Schaff JC, Shapiro BE, Shimizu TS, Spence HD, Stelling J, Takahashi K, Tomita M, Wagner J, Wang J: The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics 2003,19(4):524-531. 10.1093/bioinformatics/btg015View ArticleGoogle Scholar
- Fano U: Ionization yield of radiations. II. The fluctuations of the number of ions. Phys Rev 1947, 72: 26-29. 10.1103/PhysRev.72.26View ArticleGoogle Scholar
- Kennedy J, Eberhart R: Particle Swarm Optimization. Proceedings of the Fourth IEEE International Conference on Neural Networks, Perth, Australia 1995, 4: 1942-1948.View ArticleGoogle Scholar
- Fogel DB, Fogel LJ, Atmar JW: Meta-evolutionary programming. In 25th Asilomar Conference on Signals, Systems and Computers. Volume 1. IEEE Computer Society, Asilomar; 1992:540-545.Google Scholar
- Wolpert DH, Macready WG: No free lunch theorems for optimization. IEEE Trans Evolutionary Comput 1997, 1: 67-82. 10.1109/4235.585893View ArticleGoogle Scholar
- Lawrence MC, Jivan A, Shao C, Duan L, Goad D, Zaganjor E, Osborne J, McGlynn K, Stippec S, Earnest S, Chen W, Cobb MH: The roles of MAPKs in disease. Cell Res 2008,18(4):436-442. 10.1038/cr.2008.37View ArticleGoogle Scholar
- Kholodenko BN, Birtwistle MR: Four-dimensional dynamics of MAPK information-processing systems. Wiley Interdisciplinary Rev: Syst Biol Med 2009, 1: 28-44. 10.1002/wsbm.16Google Scholar
- Chang L, Karin M: Mammalian MAP kinase signalling cascades. Nature 2001, 410: 37-40. 10.1038/35065000View ArticleGoogle Scholar
- Komorowski M, Costa M, Rand D, Stumpf M: Sensitivity, robustness, and identifiability in stochastic chemical kinetics models. PNAS 2011,108(21):8645-8650. 10.1073/pnas.1015814108View ArticleGoogle Scholar
- Grima R: An effective rate equation approach to reaction kinetics in small volumes: Theory and application to biochemical reactions in nonequilibrium steady-state conditions. J Chem Phys 2010, 133: 035101. 10.1063/1.3454685View ArticleGoogle Scholar
- Bartels RH, Stewart GW: Solution of the matrix equation AX + XB = C [F4]. Comm ACM 1972,15(9):820-826. 10.1145/361573.361582View ArticleGoogle Scholar
- Reder C: Metabolic control theory: a structural approach. J Theor Biol 1988,135(2):175-201. 10.1016/S0022-5193(88)80073-0View ArticleGoogle Scholar
- Funahashi A, Morohashi M, Kitano H, Tanimura N: CellDesigner: a process diagram editor for gene-regulatory and biochemical networks. BIOSILICO 2003,1(5):159-162. 10.1016/S1478-5382(03)02370-9View ArticleGoogle Scholar
- Funahashi A, Matsuoka Y, Jouraku A, Morohashi M, Kikuchi N, Kitano H: CellDesigner 3.5: a versatile modeling tool for biochemical networks. Proc of the IEEE 2008,96(8):1254-1265.View ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.