 Research article
 Open Access
Bifurcation analysis informs Bayesian inference in the Hes1 feedback loop
 Catherine F Higham^{1}Email author
https://doi.org/10.1186/17520509312
© Higham; licensee BioMed Central Ltd. 2009
Received: 24 June 2008
Accepted: 26 January 2009
Published: 26 January 2009
Abstract
Background
Ordinary differential equations (ODEs) are an important tool for describing the dynamics of biological systems. However, for ODE models to be useful, their parameters must first be calibrated. Parameter estimation, that is, finding parameter values given experimental data, is an inference problem that can be treated systematically through a Bayesian framework.
A Markov chain Monte Carlo approach can then be used to sample from the appropriate posterior probability distributions, provided that suitable prior distributions can be found for the unknown parameter values. Choosing these priors is therefore a vital first step in the inference process. We study here a negative feedback loop in gene regulation where an ODE incorporating a time delay has been proposed as a realistic model and where experimental data is available. Our aim is to show that a priori mathematical analysis can be exploited in the choice of priors.
Results
By focussing on the onset of oscillatory behaviour through a Hopf Bifurcation, we derive a range of analytical expressions and constraints that link the model parameters to the observed dynamics of the system. Computational tests on both simulated and experimental data emphasise the usefulness of this analysis.
Conclusion
Mathematical analysis not only gives insights into the possible dynamical behaviour of gene expression models, but can also be used to inform the choice of priors when parameters are inferred from experimental data in a Bayesian setting.
Keywords
 Markov Chain Monte Carlo
 Hopf Bifurcation
 Markov Chain Monte Carlo Method
 Inference Problem
 Extrinsic Noise
Background
Aims
Mathematical models can help biologists understand the mechanisms and dynamics behind their experimental observations (Tomlin et al. [1]). The most widely used approach to modelling the dynamics of a genetic network is to employ systems of ordinary differential equations (ODEs) (Voit [2] and de Jong [3]). These models have biological parameters, some of which can be measured experimentally and some of which cannot. Parameter estimation, that is, recovering unknown parameters from experimental data, is an important step towards obtaining a good model that can not only explain observed results but can also be used for prediction and "what if" scenarios.
The inference of parameters, from real biological data, within a Bayesian framework is a relatively new, albeit currently very active area. Bayesian inference and Markov chain Monte Carlo (MCMC) methods have been recently advocated for the estimation of model parameters from biological systems described by ODEs (Rogers et al. [4]). This paper aims to show that a general technique, using dynamical systems analysis to inform the choice of priors, can add value. The negative feedback loop studied here is a key feature of many complex biological networks and is referred to as a motif by Alon [5]. Hence our work will be of interest whenever larger networks are dealt with by a modular 'divide and conquer' approach. A Bayesian setting links the quantity that we are interested in, the probability that our parameters take certain values given the data, to two quantities that we can assign, the probability that we would have observed the measured data if the parameters took those values and our prior biological knowledge or ignorance about these parameters (Sivia [6]). Whereas traditional parameter estimation methods are deterministic and point valued (for example, COPASI [7]), these new methods use Bayes' Theorem to assign probabilities to parameter values and can handle noise inherently.
Using a Bayesian approach to parameter estimation for nonlinear dynamical systems throws up several challenges, and these typically increase when time delays are included. Key issues are

dimensionality: models may involve several undetermined parameters,

identification: different parameter combinations may produce similar dynamics, for example, increasing a production rate may be almost equivalent to decreasing a decay rate,

local maxima: the likelihood function may have many locally optimal values.
MCMC methods [8] can go some way to addressing these difficult issues. In this work we show that a priori mathematical analysis can also be an effective tool. Our focus is on oscillatory time series, and we use the standard and widely applicable tools of Hopf bifurcation analysis and Lindstedt's method [9] in order to inform the choice of priors, thereby simplifying and focussing the parameter estimation process in order to find a particular parameter regime. More precisely, we follow the approach of Verdugo and Rand [10] in studying the case where oscillations arise via bifurcation through the delay parameter. Verdugo and Rand show that this leads to biologically realistic parameter values. We show here that inference with real data further supports this viewpoint.
The setting for our work is the biologically important instance where the expression of a gene is down regulated by its protein product. This arises, for example, with the p53 tumor suppressor protein whose intracellular activity is regulated through a feedback loop involving its transcriptional target [11]. We focus here on the case of the delayed Hes1 feedback loop featuring hes1 mRNA and Hes1 protein, where both mathematical models and quantitative experimental data are available in [12] and [13].
Biological Data
It has been observed that mRNAs for Notch signalling molecules such as the bHLH factor Hes1 oscillate with 2hour cycles during somite segmentation, although the molecular mechanism of such oscillations remains to be fully determined. Hirata et al. [13] investigated the oscillations of mRNAs for Notch signalling molecules by examining the time course of hes1 mRNA and its Hes1 protein in detail. They observed that a single serum treatment induced a 2hour cycle oscillation of hes1 mRNA in a variety of cultured cells and that Hes1 protein also oscillated in a 2hour cycle after the single serum treatment with a delay of about 15 minutes relative to the hes1 mRNA oscillation. The half lives of hes1 mRNA and Hes1 protein were measured and the proteases for Hes1 protein degradation were identified.
Hirata et al. confirmed experimentally that degradation of Hes1 protein is required for hes1 mRNA increase and that de novo production of the protein is required for reduction of hes1 mRNA. These facts together support the theory that Hes1 is an essential component of a two hour cycle clock. They showed that the same mechanism applies to hes1 mRNA oscillation in the presomitic mesoderm.
A system featuring negative feedback and described by only two ODEs (one for each species, in this case hes1 mRNA and Hes1 protein), can be shown never to generate sustained oscillations [14]. Hirata et al. proposed that the observed oscillations could be modelled by introducing a Hes1 interacting factor as a third molecular species, but there is no experimental evidence for such an interacting factor. Monk [12] showed that observed oscillatory behaviour can be explained by incorporating a time delay.
Model for Hes1 Feedback Loop
where μ_{ m }and μ_{ p }are the rates of degradation of mRNA and protein, respectively, p_{0} is the normalised repression threshold, n is a hill coefficient and τ is a constant delay caused by transcription and translation. In this system, the Hes1 gene transcribes hes1 mRNA which passes from the nucleus to the cytoplasm. There it is translated into Hes1 protein. An unusual feature is that the Hes1 protein binds to the gene promoter and represses the transcription of hes1 mRNA. Following Monk, the equations have been rewritten so that the delay appears in the equation describing the regulation of transcription and not in the equation describing protein translation or synthesis. Throughout the rest of this paper, we will refer to τ as a transcriptional delay. The equations have also been rescaled so that the dynamics of the model are determined by five parameters. Using an ad hoc parameter fitting procedure, Monk showed that this model reproduced the broad features of the oscillation of hes1 mRNA and Hes1 protein in cultured cells observed by Hirata et al. [13].
Thus Monk [12] argues that the observed oscillatory behaviour is best accounted for by the introduction of a delay parameter that acts as a proxy for the many subprocesses that make up transcription and translation, rather than by introducing an unknown third agent [13]. By performing mathematical analysis of the model (1)–(2), we will show that the convincing but ad hoc parameter fitting exercise in [12] can be extended to a fully Bayesian setting.
Mathematical Analysis
Verdugo and Rand [10] gave some mathematical analysis of Monk's model [12] for the Hes1 feedback loop. In particular, they derived closed form approximations for the amplitude and frequency of oscillation, where oscillatory behaviour is assumed to arise through Hopf bifurcation in the delay parameter. The analysis in [10] applies to the case where the decay rates of hes1 mRNA and Hes1 protein, key components of the feedback, are equal, that is, μ_{ m }= μ_{ p }in (1)–(2). In this work, we study the more realistic case where the decay rates are allowed to be different, also focussing on oscillatory behaviour. By allowing these degradation rates to differ, we expose some interesting features of the nonlinear system under consideration; namely that oscillations of this type only occur when the difference between degradation rates is small. The novelty of this work lies in using bifurcation tools to aid what turns out to be a very difficult inference problem. Details of the mathematical analysis can be found in the methods section. Our results are of interest in their own right as a means to understand how the system dynamics are driven by the model parameters, but our main aim here is to inform the corresponding Bayesian inference problem.
Experiments
In order to illustrate the use of mathematically informed priors when inferring model parameters we conducted two experiments. Experiment 1 uses synthetic data, computer generated from the mathematical equations describing the Monk model with known parameters. Experiment 2 looks at the published data of Hirata et al. Data is noisy and we have noted the error bars in the Hirata data when making assumptions about this type of extrinsic noise. Bayesian methods are very useful when dealing with this type of experimental noise because we can include noise as a hyper parameter, treat it as a nuisance parameter and integrate it out. Intrinsic noise will also come into play at the molecular level [15] and there is recent work on the subject of inferring parameters in stochastic models [16]. However, we do not believe that there is enough data to make it feasible to fit a stochastic model that would take into account intrinsic noise. We therefore use the ODEbased model in line with other authors [12] and [15]. Consequently, we are dealing with extrinsic and experimental noise. To make the synthetic data in Experiment 1 more realistic, Gaussian noise is added, based on the absolute size of the error bars of Hirata with a mean of zero and a standard deviation of 0.2. In our example, we have used larger initial conditions so the relative noise is lower. However, investigations with a higher level of noise to match the relative values from Hirata give broadly similar conclusions. The justification for adding Gaussian noise to make the data more realistic lies with the Central Limit Theorem, which shows that under appropriate conditions the overall effect of extrinsic noise sources will be normal [4]. An advantage of using synthetic data is that the parameters to be inferred are known and so the accuracy of the method to infer them can be evaluated. For Experiment 2, two of the parameters values, μ_{ m }and μ_{ p }, have been measured independently. So mathematically informed priors are complemented with biological priors where available.
In Experiment 2 (Hirata Data) there is data for both species m and p and so in Experiment 1 (Synthetic Data) we generate data for both species m and p as well. We will show that some aspects of the data for p inform the mathematical analysis, however the actual inference of parameters for both m and p is based on the data for m only. This is because protein expression levels are harder to measure accurately in practice so it is useful to establish what can be inferred from the more accurately measured species alone. With the Hirata data, there is also the added complication that the protein levels measured have been scaled in an undisclosed manner [15].
The location in time of the sampling points is clearly an important issue. For species that oscillate, it may be difficult to know in advance which time points will provide most information. In this work the real data is already available and hence our initial test was chosen to match those details. Generally, for oscillatory behaviour, we would expect that the key constraint on the sample times is that the sampling frequency should be sufficiently large.
Methodology
We have some data and a proposed model. We use a Bayesian framework to infer the parameters underlying the model. In a Bayesian framework, the posterior pdfs are given by the prior pdfs modified by the likelihood of the data given those parameters. More precisely the posterior pdfs are proportional to the product of the likelihood and the priors. The posterior pdf for a particular parameter depends on the values of all the other parameters and are shown for individual parameters by integrating over the other parameter values. Hence the posterior pdfs that we show are marginal distributions for each parameter and not distributions for fixed values of the remaining parameters. The resulting equations for the posterior pdfs are intractable so we use MCMC methods to sample from the posterior distributions. This section spells out further details of this methodology which differs slightly for Experiment 1 and Experiment 2.
Bayesian Inference
Bayesian methods have several advantages over other approaches to parameter estimation [17]: use of background information, the ability to include uncertainity in all parameter values and the ease of making inferences about some parameters irrespective of the values of others. The development of highspeed computing means that the potential of Bayesian methods is being realised, but there are still many computational difficulties. In this section, we summarise the key features of Bayesian inference and highlight the technical decisions which have to be made. For more information about Bayesian analysis we refer readers to [6] and reviews such as [17].
Our initial state of knowledge (or ignorance) about the model parameters is encapsulated in the prior probability. The priors are then modified by the experimental measurements through the likelihood function to give us the posterior probability, which represents our new state of knowledge about the value of the parameters. A likelihood function needs to be chosen and after using a normal distribution for Experiment 1, we chose to base the likelihood function for Experiment 2 on the chisquared distribution. The chisquared distribution, χ^{2}, has longer tails than the normal distribution and so encourages the inclusion of outliers in the data, which is important for the estimation of model parameters which influence in particular the amplitude of the oscillations. This proved to be an important consideration in Experiment 2, possibly because of experimental noise levels.
For Experiment 1, we therefore suppose that the probability of the k th datum having a value x_{ k }is given by N (f_{ k }, σ) where f_{ k }= f(θ, t_{ k }) is the true value of the function of parameters θ = {p_{0}, n, μ_{ m }, μ_{ p }, τ} of interest, and the variance, σ, accounts for the error in its measurement. In our case, f_{ k }is the mRNA level m (p_{0}, n, μ_{ m }, μ_{ p }, τ) at time t_{ k }from the delay differential equation model (1)–(2).
For Experiment 2, we suppose that the probability of the k th datum having a value x_{ k }is given by χ^{2} (f_{ k }, v) where the degree of freedom v, accounts for the error in its measurement. Given that the Hirata data has been scaled, it will also be necessary to infer an additional scaling parameter, which we refer to as k_{ s }, so θ = {p_{0}, n, μ_{ m }, μ_{ p }, τ, k_{ s }}.
For both experiments, we assume that the value of σ and v is known. Our inference about the value of θ is expressed, therefore, by the posterior pdf p (θ  {x_{ k }}, I) where I denotes all other background information; to help us calculate it, we use Bayes' Theorem [6]
p (θ  {xk}, v, I) ∝ p ({x_{ k }}  θ, v, I) × p (θ, v, I).
The relation in (3) is expressed using proportionality, because the term p (dataI) has been omitted from the denominator in the right hand side. This is fine for data analysis problems, such as this one, involving parameter estimation, since the missing denominator is simply a normalisation constant not depending explicitly on the parameters. However if we were to consider model selection, this term would play a crucial role, and is thus given the special name of evidence.
According to (3), if we multiply the prior stated in (5) by the likelihood function (4), we obtain the posterior pdf. Our approach here will be to deal with the posterior pdf numerically, using MCMC.
MCMC Method
Experiment 1 (Synthetic Data).
θ  sp1  sp2  sp3  σ  acc rate  $\widehat{R}$ 

P _{0}  85  95  91  0.5  26  1.0009 
n  4  5  5.3  0.7  34  1.0001 
μ _{ m }  0.024  0.028  0.022  0.1  43  1.0004 
μ _{ p }  0.036  0.034  0.033  0.1  33  1.0000 
τ  19  19.2  18.7  0.01  20  1.0000 
Experiment 2 (Hirata Data).
θ  sp1  sp2  sp3  σ  acc rate  $\widehat{R}$ 

P _{0}  99  103  103  0.55  29  1.0000 
n  5  4.8  5.3  0.7  34  1.0003 
μ _{ m }  0.028  0.03  0.03  0.3  28  1.0001 
μ _{ p }  0.028  0.03  0.031  0.5  31  1.0002 
τ  18  19  18.5  0.5  27  1.0001 
k _{ s }  2.5  2.2  2  0.2  27  1.0000 
Results and discussion
Mathematical Analysis
We summarise the main analytical results here. Further details about the derivations of the formulae can be found in the Methods section.
The result (7) characterises the ratio between mRNA and protein amplitudes, but not their absolute values. To proceed we use Lindstedt's method which is a technique for uniformly approximating periodic solutions to ODEs when regular perturbation approaches fail. We regard the frequency as unknown in advance and solve by demanding that an appropriate series expansion contains no secular terms. The result is closed form approximate expressions for the amplitude and frequency of oscillation expressed in terms of the model parameters. For the derivation of these formulae and more details, see the Methods section. In summary, our approximation to A^{2} and the observed frequency – may be written in the form
A^{2} = f_{ A }(n, p_{0}, μ_{ m }, μ_{ p }) Δ,
Ω = ω  f_{Ω}(n, p_{0}, μ_{ m }, μ_{ p }) Δ,
where f_{ A }and f_{Ω} are positive functions defined in terms of the model parameters and Δ = τ  τ_{ cr }. As Δ increases we would expect A and B to increase proportionally to $\sqrt{\Delta}$ and Ω to decrease.
The frequency of oscillation, ω, takes a maximum that can be explicitly defined in terms of the model parameters. We may show that for any μ_{ m }μ_{ p }= c, a unique maximum value of ω occurs when μ_{ m }= μ_{ p }= $\sqrt{c}$. We may then find the value of c and hence μ_{ m }and μ_{ p }for which the maximum of ω occurs. We can use this information combined with equation (15) to deduce that ω > Ω. This result allows us to make an informed judgement about ω given Ω, which can be used in the setting of priors.
Biology meaningful solutions must also have ω > 0 and τ_{ cr }> 0. Imposing these conditions in equation (9) gives
ω^{2} > μ_{ m }μ_{ p },
ωμ_{ m }+ ωμ_{ p }<K,
The data can tell us something about the parameters directly. Given sufficient data points from complete periodic cycles, the model parameter μ_{ p }can be approximated by the ratio of the average of m to the average of p. Further details can be found in the Methods section; see equation (70).
We now show how our results can be used to inform the choice of priors in our two inference experiments.
Experiment 1 (synthetic data)
Experiment 1 (Synthetic Data).
θ  mean  std dev 

P _{0}  90  10 
n  5  1 
μ _{ m }  0.025  0.001 
μ _{ p }  0.035  0.001 
τ  18.5  1 
Results
MCMC provides samples and in order to recover the posterior pdfs these samples are binned, smoothed and normalised. We can also simulate the dynamic activity associated with each of these samples using the associated parameter values. The variance of this dynamic activity at each time point indicates the likely range of activity that could have produced the data. The shaded area in Figure 3 shows this range (5–95 percentile) for Experiment 1. The fit of this area in terms of position and width shows how well we have recovered the model parameters from the data.
Experiment 2 (biological data)
Experiment 2 (Hirata Data).
θ  mean  std dev 

P _{0}  100  10 
n  5  1 
μ _{ m }  0.0288  0.002 
μ _{ p }  0.0311  0.004 
τ  19  1 
k _{ s }  2.2  0.1 
Results
The shaded area in Figure 1 shows the fit of the parameter estimates within the 5–95 percentile for Experiment 2. The early data points are covered by the prediction. Compared to Monk's results, there is broad agreement, however we predict a slightly longer delay (around 18.9 compared to Monk's 18.5) and slightly different values for μ_{ m }and μ_{ p }(around 0.0285 and 0.0305, respectively compared to Monk's 0.03 for both). We also inferred a different scaling factor which contributed significantly to the improved match to the Hirata et al. data that can be seen in Figure 1.
Conclusion
We have shown that mathematical analysis can give insights into the possible dynamical behaviour of gene expression models. We derived a range of analytical expressions and constraints that link the model parameters to the observed dynamics of the system. The basis of the analysis is the assumption that oscillatory behaviour arises through Hopf bifurcation in the delay parameter. It is then possible to say for which range of parameters and to what extent this behaviour will occur. For the problem of inferring parameter values from data, we can use this analysis to inform the choice of priors. For these types of problems, where different parameter values can result in similar activity, identification of the parameters through the data can be ambiguous. MCMC methods can go some way towards resolving these difficulties but we showed here that mathematically informed priors can add value too. On the experimental data from Hirata et al., a full Bayesian MCMC computation improved on the parameter fit published by Monk, predicting nonequal values for mRNA and protein decay rates and a longer time delay for transcription.
Methods
We describe here in more detail the mathematical and computational aspects of this work. We generalise the approach of Verdugo and Rand to the more biologically realistic case where the degradation rates of mRNA and protein are not assumed to be equal. The analysis reveals that we would expect to see oscillations when the difference between μ_{ m }and μ_{ p }is small which justifies, generally, the assumption that the degradation rates are almost equal.
Stability of Equilibrium
We now look for a solution of the form η = e^{ λt }, which leads to
λ^{2} + (μ_{ m }+ μ_{ p }) λ + μ_{ m }μ_{ p }+ Ke^{λτ}= 0.
Examination of equation (26) reveals that since μ_{ m }+ μ_{ p }> 0 and $\sqrt{{({\mu}_{m}+{\mu}_{p})}^{2}4({\mu}_{m}{\mu}_{p}+K)}<{\mu}_{m}+{\mu}_{p}$, the real part of λ will always be negative, implying stability. For
4 (μ_{ m }μ_{ p }+ K) > (μ_{ m }+ μ_{ p })^{2}
λ has a nonzero imaginary part as well as a negative real part and the equilibrium point (m*, p*) will be a stable spiral. When μ_{ m }= μ_{ p }, the inequality (27) becomes K > 0 which always holds. However when μ_{ m }≠ μ_{ p }this inequality becomes 4K > (μ_{ m } μ_{ p })^{2}. In other words, we only have a stable spiral when μ_{ m } μ_{ p } is sufficiently small. For larger values of μ_{ m } μ_{ p } we will have two negative real roots for λ, resulting in a fixed stable point with no oscillations. We conclude that, overall, the nondelay model always has a linearly stable fixed point.
Following the approach in [10], we now assume, by continuity, that as τ increases the roots λ will cross the imaginary axis (this will not be at the origin because μ_{ m }> 0 and μ_{ p }> 0 and therefore the real and imaginary parts of λ cannot both be zero) at some critical value τ = τ_{ cr }, and for τ > τ_{ cr }the steady state (m*, p*) will lose stability giving rise to a Hopf bifurcation. We thus assume that for τ = τ_{ cr }the system (22)–(23) will exhibit a pair of pure imaginary eigenvalues ± ω i corresponding to the solution
ζ(t) = B cos (ωt + ϕ),
η(t) = A cos ωt,
where A and B are the amplitudes of the η (t) and ζ (t) oscillations, and where ϕ is a phase angle. More generally, for values of delay τ close to τ_{ cr }, the nonlinear system is expected to relax to a periodic solution which can be written in the approximate form of equations (28) and (29).
Substituting equations (28) and (29) into equations (22) and (23) and matching time dependent terms results in the explicit formulae (6)–(11), which generalise those in [10].
Lindstedt's method
where the coefficients K, H_{2} and H_{3} are obtained by Taylor series expansion of the nonlinear term.
Assuming that the true solution is periodic, our goal is to find an analytical approximation that is valid for all t. The key idea is to regard the frequency ω as unknown in advance, and to solve for it by demanding that η contains no secular terms. We introduce a small parameter ϵ, and let Δ = ϵ^{2}δ, with
η = ϵu,
τ = τ_{ cr }+ Δ = τ_{ cr }+ ϵ^{2}δ.
We then stretch time by replacing the independent variable t by σ, where
σ = Ωt.
where u_{ d }= u (σ  Ωτ). Expanding Ω in a power series in ϵ, omitting the O (ϵ) term for convenience, since it turns out to be zero, we have
Ω = ω + ϵ^{2}k_{2} + ...
Next we expand the delay term u_{ d }
u_{ d }= u (σ  ωτ_{ cr })  ϵ^{2} (k_{2} τ_{ cr }+ ωδ) u' (σ  ωτ_{ cr }) + O(ϵ^{3})
Finally we expand u (t) in a power series in ϵ to obtain
u (σ) = u_{0} (σ) + ϵu_{1} (σ) + ϵ^{2}u_{2} (σ) + ...
u_{1}(σ) = m_{1} sin 2σ + m_{2} sin 2σ + m_{3},
Letting
α = 4ω^{2} + K cos2ωτ_{ cr }+ μ_{ m }μ_{ p },
Let
c_{1} = μ_{ m } μ_{ p }+ Kτ_{ cr }cosωτ_{ cr },
Multiplying (56) by ϵ^{2} and substituting back into (35) gives
Ω = ω + k_{2}Δ,
which has the form of equation (15) (noting that k_{2} is negative).
Maximum ω
which simplifies to
(n^{2}  1) β^{2}  (n^{2} + 2) β  1 = 0.
Can the data tell us anything about the parameters directly?
This allows us, from appropriate data sets, to infer the parameter μ_{ p }by averaging data from m and p. This, of course, presupposes that the period T can be estimated accurately and that there is enough data to form good approximations to the averages.
We remark that μ_{ m }cannot be approximated this way because of the nonlinearity in (1).
MCMC (Metropolis Hastings) Method
σ_{ j }values for the proposal rates are chosen in order that the acceptance rates are between 20% and 35%.
A pseudocode outline for the MCMC Metropolis Hastings algorithm is as follows:
step 1 Initialise ${\theta}_{j}^{(0)}$
step 2 For i = 0 to N:
Sample r ~ $\mathcal{N}$(0, 1)
Sample ${\theta}_{j}^{\ast}~q\left({\theta}_{j}^{\ast}{\theta}_{j}^{(i)}\right)$ (choosing j at random
Declarations
Acknowledgements
CFH gratefully acknowledges the support of a Daphne Jackson Fellowship funded by the Leverhulme Trust, based in the Department of Computing Science at the University of Glasgow. The author benefitted from useful discussions with Mark Girolami, Raya Khanin, Simon Rogers, Dirk Husmeier and Neil Lawrence and from helpful reviewers' comments.
Authors’ Affiliations
References
 Tomlin C, Axelrod J: Biology by numbers: mathematical modelling in developmental biology. Nature Reviews Genetics. 2007, 8: 331340. 10.1038/nrg2098View ArticlePubMedGoogle Scholar
 Voit E: Computational Analysis of Biochemical Systems. 2000, Cambridge University PressGoogle Scholar
 de Jong H: Modeling and simulation of genetic regulatory systems: a literature review. J Comput Biol. 2002, 9: 67103. 10.1089/10665270252833208View ArticlePubMedGoogle Scholar
 Rogers S, Khanin R, Girolami M: Bayesian modelbased inference of transcription factor activity. BMC Bioinformatics. 2006, 8 Suppl 2: S2Google Scholar
 Alon U: An Introduction to Systems Biology: Design Principles of Biological Circuits. 2007, Chapman and Hall/CRCGoogle Scholar
 Sivia D: Data Analysis: A Bayesian Tutorial. 2006, Oxford University Press, 2Google 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: 306774. 10.1093/bioinformatics/btl485View ArticlePubMedGoogle Scholar
 Andrieu C, de Freitas N, Doucet A, Jordan MI: An introduction to MCMC for machine learning. Machine Learning. 2003, 50: 543. 10.1023/A:1020281327116.View ArticleGoogle Scholar
 Strogatz SH: Nonlinear Dynamics and Chaos. 1994, Addison WesleyGoogle Scholar
 Verdugo A, Rand R: Hopf bifurcation in a DDE model of gene expression. Communications in Nonlinear Science and Numerical Simulation. 2008, 13: 235242. 10.1016/j.cnsns.2006.05.001.View ArticleGoogle Scholar
 BarOr R, Maya R, Segel L, Alon U, Levine A, Oren M: Generation of oscillations by the p53Mdm2 feedback loop: a theoretical and experimental study. Proceedings of the National Academy of Sciences USA. 2000, 97: 1125011255. 10.1073/pnas.210171597.View ArticleGoogle Scholar
 Monk NA: Oscillatory Expression of Hes1, p53, and NFkB Driven by Transcriptional Time Delays. Current Biology. 2003, 13: 14091413. 10.1016/S09609822(03)004949View ArticlePubMedGoogle Scholar
 Hirata H, Yoshiura S, Ohtsuka T, Bessho Y, Harada T, Yoshikawa K, Kageyama R: Oscillatory Expression of the bHLH Factor Hes1 Regulated by a Negative Feedback Loop. Science. 2002, 298:Google Scholar
 Bernard S, Čajavec B, PujoMenjouet L, Mackey MC, Herzel H: Modeling transcriptional feedback loops: the role of Gro/TLE1 in Hes1 oscillations. Philos Transact A Math Phys Eng Sci. 2006, 364 (1842): 115570. 10.1098/rsta.2006.1761View ArticleGoogle Scholar
 Barrio M, Burrage K, Leier A, Tian T: Oscillatory Regulation of Hes1: Discrete Stochastic Delay Modelling and Simulation. PLoS Computational biology. 2006, 2: 10171030. 10.1371/journal.pcbi.0020117.View ArticleGoogle Scholar
 Wilkinson DJ: Stochastic Modelling for Systems Biology. 2006, Chapman and Hall/CRCGoogle Scholar
 Beaumont MA, Rannala B: The Bayesian revolution in genetics. Nature Reviews Genetics. 2004, 5: 251261. 10.1038/nrg1318View ArticlePubMedGoogle Scholar
 Gelman A, Carlin JB, Stern HS, Rubin DB: . Bayesian Data Analysis. 2004, Chapman and Hall CRCGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.