Skip to main content
  • Research article
  • Open access
  • Published:

Bifurcation analysis informs Bayesian inference in the Hes1 feedback loop



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.


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.


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.



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 2-hour 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 2-hour cycle oscillation of hes1 mRNA in a variety of cultured cells and that Hes1 protein also oscillated in a 2-hour 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.

The Hirata data comprises scaled hes1 mRNA expression levels every 30 minutes over a 12 hour period; see the discrete points in Figure 1.

Figure 1
figure 1

Experiment 2: Hirata data comprises 25 mRNA values taken every 30 minutes over a 12 hour period. Initial conditions are assumed to hold for 0 <t <τ.

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

Monk [12] suggests that the observed oscillations are due to the non-instantaneous nature of transcriptional and translational delays, and proposes a mathematical model where a delay is introduced to account for this feature. Monk's model was able to explain, via numerical simulations, the oscillation of hes1 mRNA and Hes1 protein in cultured cells observed by Hirata et al. Letting m (t) and p (t) denote the concentration of hes1 mRNA and Hes1 protein at time t, respectively, the model proposed by Monk for the Hes1 feedback loop takes the form

m ˙ = 1 1 + ( p ( t τ ) / p 0 ) n μ m m , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafmyBa0MbaiaacqGH9aqpjuaGdaWcaaqaaiabigdaXaqaaiabigdaXiabgUcaRiabcIcaOiabdchaWjabcIcaOiabdsha0jabgkHiTiabes8a0jabcMcaPiabc+caViabdchaWnaaBaaabaGaeGimaadabeaacqGGPaqkdaahaaqabeaacqWGUbGBaaaaaOGaeyOeI0IaeqiVd02aaSbaaSqaaiabd2gaTbqabaGccqWGTbqBcqGGSaalaaa@4653@
p ˙ = m μ p p , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafmiCaaNbaiaacqGH9aqpcqWGTbqBcqGHsislcqaH8oqBdaWgaaWcbaGaemiCaahabeaakiabdchaWjabcYcaSaaa@3689@

where μ m and μ p are the rates of degradation of mRNA and protein, respectively, p0 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 sub-processes 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.


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 ODE-based 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.


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 high-speed 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 chi-squared distribution. The chi-squared 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 θ = {p0, 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 (p0, 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 θ = {p0, 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 (data|I) 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.

As in [4], we assume that the data are independent, so that the measurement of one datum does not influence what we can infer about the outcome of another (when given the values of θ). The likelihood function is then given by the product of the probabilities of obtaining the N individual data:

p ( { x k } | θ , I ) = k = 1 N p ( x k | θ , I ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiCaaNaeiikaGIaei4EaSNaemiEaG3aaSbaaSqaaiabdUgaRbqabaGccqGG9bqFcqGG8baFcqaH4oqCcqGGSaalcqWGjbqscqGGPaqkcqGH9aqpdaqeWbqaaiabdchaWjabcIcaOiabdIha4naaBaaaleaacqWGRbWAaeqaaOGaeiiFaWNaeqiUdeNaeiilaWIaemysaKKaeiykaKIaeiOla4caleaacqWGRbWAcqGH9aqpcqaIXaqmaeaacqWGobGta0Gaey4dIunaaaa@4E6C@

We assign a Gaussian pdf as prior for each parameter θ j , with mean μ j and variance σ j , where j = 1:5 in Experiment 1 and j = 1:6 in Experiment 2:

p ( θ j | σ , I ) = p ( θ j | I ) = 1 σ j 2 π exp ( ( θ j μ j ) 2 2 σ j 2 ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiCaaNaeiikaGIaeqiUde3aaSbaaSqaaiabdQgaQbqabaGccqGG8baFcqaHdpWCcqGGSaalcqWGjbqscqGGPaqkcqGH9aqpcqWGWbaCcqGGOaakcqaH4oqCdaWgaaWcbaGaemOAaOgabeaakiabcYha8jabdMeajjabcMcaPiabg2da9KqbaoaalaaabaGaeGymaedabaGaeq4Wdm3aaSbaaeaacqWGQbGAaeqaamaakaaabaGaeGOmaiJaeqiWdahabeaaaaGccyGGLbqzcqGG4baEcqGGWbaCdaqadaqaaiabgkHiTKqbaoaalaaabaGaeiikaGIaeqiUde3aaSbaaeaacqWGQbGAaeqaaiabgkHiTiabeY7aTnaaBaaabaGaemOAaOgabeaacqGGPaqkdaahaaqabeaacqaIYaGmaaaabaGaeGOmaiJaeq4Wdm3aa0baaeaacqWGQbGAaeaacqaIYaGmaaaaaaGccaGLOaGaayzkaaGaeiOla4caaa@61DF@

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

MCMC techniques are used to sample from distributions and more generally to solve optimisation and integration problems in large dimensional spaces (Andrieu et al. [8]). Here, we use MCMC to sample from the posterior pdf defined in (3). The approach is, essentially, to compute candidate values for our unknown parameters by constructing a Markov chain with an appropriate invariant measure. The candidate values are chosen randomly using a proposal function and are either accepted or rejected depending on an acceptance probability. The frequency at which a particular parameter is successfully sampled is indicative of its posterior probability. Once we have enough samples from the solution space, they can be binned, smoothed and normalised to indicate the relative posterior likelihood of a given parameter. Sampling can start from any particular parameter value. The proposal function is chosen so that samples are drawn widely and performance can be monitored by the acceptance rate. With Markov chains, after a burn in period, the memory of the initial starting point is lost. At this point sampling will be exclusively from the required distribution. We need to establish at which point this has occurred. Gelman et al. [18] propose setting off parallel sampling chains and monitoring their convergence with a statistic R ^ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafmOuaiLbaKaaaaa@2D12@ . We have followed this convention. More details can be found in the Methods section and Tables 1 (Experiment 1) and 2 (Experiment 2).

Table 1 Experiment 1 (Synthetic Data).
Table 2 Experiment 2 (Hirata Data).

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.

Following [10] we can derive an expression for the equilibrium state of m and p (m* and p*, respectively, see equations (21)–(20)) of model (1)–(2) and show that by increasing the delay the system moves from a stable state to an unstable state. At a critical value τ cr of the time delay a Hopf Bifurcation occurs. For values of delay τ close to τ cr , the nonlinear system is expected to exhibit a periodic solution which can be expressed for both the mRNA and protein in terms of their amplitudes, A and B, respectively, their oscillatory frequency ω (at τ cr ) and their phase difference ϕ. Algebraic manipulation of these expressions results in explicit formulae linking the model's biological parameters which generalise those in [10]. Specifically, at the Hopf bifurcation we have

tan ϕ = ω μ p , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGagiiDaqNaeiyyaeMaeiOBa4Maeqy1dyMaeyypa0tcfa4aaSaaaeaacqaHjpWDaeaacqaH8oqBdaWgaaqaaiabdchaWbqabaaaaOGaeiilaWcaaa@39A6@
B A = ω 2 + μ p 2 , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqWGcbGqaeaacqWGbbqqaaGccqGH9aqpdaGcaaqaaiabeM8a3naaCaaaleqabaGaeGOmaidaaOGaey4kaSIaeqiVd02aa0baaSqaaiabdchaWbqaaiabikdaYaaaaeqaaOGaeiilaWcaaa@38F9@
w 2 = ( μ m 2 + μ p 2 ) ± ( μ m 2 + μ p 2 ) 2 4 ( μ m 2 μ p 2 K 2 ) 2 MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaem4DaC3aaWbaaSqabeaacqaIYaGmaaGccqGH9aqpcqGHsisldaWcaaqaaiabcIcaOiabeY7aTnaaDaaaleaacqWGTbqBaeaacqaIYaGmaaGccqGHRaWkcqaH8oqBdaqhaaWcbaGaemiCaahabaGaeGOmaidaaOGaeiykaKIaeyySae7aaOaaaeaacqGGOaakcqaH8oqBdaqhaaWcbaGaemyBa0gabaGaeGOmaidaaOGaey4kaSIaeqiVd02aa0baaSqaaiabdchaWbqaaiabikdaYaaakiabcMcaPmaaCaaaleqabaGaeGOmaidaaOGaeyOeI0IaeGinaqJaeiikaGIaeqiVd02aa0baaSqaaiabd2gaTbqaaiabikdaYaaaaeqaaOGaeqiVd02aa0baaSqaaiabdchaWbqaaiabikdaYaaakiabgkHiTiabdUealnaaCaaaleqabaGaeGOmaidaaOGaeiykaKcabaGaeGOmaidaaaaa@5A6F@
τ c r = 1 ω arctan ( ω μ p + ω μ m ω 2 μ m μ p ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGH9aqpjuaGdaWcaaqaaiabigdaXaqaaiabeM8a3baakiGbcggaHjabckhaYjabcogaJjabcsha0jabcggaHjabc6gaUnaabmaajuaGbaWaaSaaaeaacqaHjpWDcqaH8oqBdaWgaaqaaiabdchaWbqabaGaeyOeI0IaeqyYdCNaeqiVd02aaSbaaeaacqWGTbqBaeqaaaqaaiabeM8a3naaCaaabeqaaiabikdaYaaacqGHsislcqaH8oqBdaWgaaqaaiabd2gaTbqabaGaeqiVd02aaSbaaeaacqWGWbaCaeqaaaaaaOGaayjkaiaawMcaaiabcYcaSaaa@55C3@
sin ω τ c r = ω μ p + ω μ m K , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGagi4CamNaeiyAaKMaeiOBa4MaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGH9aqpjuaGdaWcaaqaaiabeM8a3jabeY7aTnaaBaaabaGaemiCaahabeaacqGHRaWkcqaHjpWDcqaH8oqBdaWgaaqaaiabd2gaTbqabaaabaGaem4saSeaaOGaeiilaWcaaa@4578@
cos ω τ c r = ω 2 μ m μ p K , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGagi4yamMaei4Ba8Maei4CamNaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGH9aqpjuaGdaWcaaqaaiabeM8a3naaCaaabeqaaiabikdaYaaacqGHsislcqaHjpWDdaWgaaqaaiabd2gaTbqabaGaeqiVd02aaSbaaeaacqWGWbaCaeqaaaqaaiabdUealbaakiabcYcaSaaa@44D7@


K = n β p ( 1 + β ) 2 , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaem4saSKaeyypa0tcfa4aaSaaaeaacqWGUbGBcqaHYoGyaeaacqWGWbaCdaahaaqabeaacqGHxiIkaaGaeiikaGIaeGymaeJaey4kaSIaeqOSdiMaeiykaKYaaWbaaeqabaGaeGOmaidaaaaakiabcYcaSaaa@3B89@


β = ( p p 0 ) n . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqOSdiMaeyypa0ZaaeWaaKqbagaadaWcaaqaaiabdchaWnaaCaaabeqaaiabgEHiQaaaaeaacqWGWbaCdaWgaaqaaiabicdaWaqabaaaaaGccaGLOaGaayzkaaWaaWbaaSqabeaacqWGUbGBaaGccqGGUaGlaaa@386D@

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 A2 and the observed frequency – may be written in the form

A2 = f A (n, p0, μ m , μ p ) Δ,

Ω = ω - fΩ(n, p0, μ m , μ p ) Δ,

and substituting (14) into (7) gives an approximation for B2

B 2 = ( ω 2 + μ p 2 ) f A ( n , p 0 , μ m , μ p ) Δ , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOqai0aaWbaaSqabeaacqaIYaGmaaGccqGH9aqpcqGGOaakcqaHjpWDdaahaaWcbeqaaiabikdaYaaakiabgUcaRiabeY7aTnaaDaaaleaacqWGWbaCaeaacqaIYaGmaaGccqGGPaqkcqWGMbGzdaWgaaWcbaGaemyqaeeabeaakiabcIcaOiabd6gaUjabcYcaSiabdchaWnaaBaaaleaacqaIWaamaeqaaOGaeiilaWIaeqiVd02aaSbaaSqaaiabd2gaTbqabaGccqGGSaalcqaH8oqBdaWgaaWcbaGaemiCaahabeaakiabcMcaPiabfs5aejabcYcaSaaa@4CF5@

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 Δ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaOaaaeaacqqHuoaraSqabaaaaa@2D56@ 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 = c MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaOaaaeaacqWGJbWyaSqabaaaaa@2D3F@ . 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,

from which we deduce that

μ m μ p < 1 p 0 ( 2 n 2 ) 1 / n ( n n 2 ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqiVd02aaSbaaSqaaiabd2gaTbqabaGccqaH8oqBdaWgaaWcbaGaemiCaahabeaakiabgYda8KqbaoaalaaabaGaeGymaedabaGaemiCaa3aaSbaaeaacqaIWaamaeqaamaabmaabaWaaSaaaeaacqaIYaGmaeaacqWGUbGBcqGHsislcqaIYaGmaaaacaGLOaGaayzkaaWaaWbaaeqabaGaeGymaeJaei4la8IaemOBa4gaamaabmaabaWaaSaaaeaacqWGUbGBaeaacqWGUbGBcqGHsislcqaIYaGmaaaacaGLOaGaayzkaaaaaiabc6caUaaa@4823@

This is a key result because we have defined a computable, finite region of the μ m - μ p plane in which oscillations can occur and in which we can search for oscillatory solutions; see Figure 2.

Figure 2
figure 2

Constrained region of the μ m - μ p plane where the inequality (19) predicts that oscillations may occur for arbitary τ. (In this picture n = 5 and p0 = 90).

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)

The synthetic data in Figure 3 comprises 49 mRNA values. Details of the parameter values used with the Monk equations and the added noise can be found in the caption of Figure 3. The corresponding 49 protein values were also generated and feature in Figure 4 where they are considered for the μ p prior. The protein values, as for Experiment 2, are not used in the inference problem. The following section describes how the mathematical analysis is used to find priors for the model parameters.

Figure 3
figure 3

Experiment 1: Synthetic data for the inference experiment comprising 49 mRNA values and 49 protein values (not shown here as they are not used directly in the inference problem). Data comes from the model with the following parameter values: p0 = 90, n = 5, μ m = 0.025, μ p = 0.035 and τ = 19.5. Initial conditions are 3 for the mRNA and 100 for the protein. Independent Gaussian noise (mean = 0 and standard deviation = 0.2) is then added to make the data more realistic.

Figure 4
figure 4

Experiment 1: We use the observed data to estimate μ p = m a v e p a v e = 4.9783 144.8241 = 0.0344 MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqiVd02aaSbaaSqaaiabdchaWbqabaGccqGH9aqpjuaGdaWcaaqaaiabd2gaTnaaBaaabaGaemyyaeMaemODayNaemyzaugabeaaaeaacqWGWbaCdaWgaaqaaiabdggaHjabdAha2jabdwgaLbqabaaaaOGaeyypa0tcfa4aaSaaaeaacqaI0aancqGGUaGlcqaI5aqocqaI3aWncqaI4aaocqaIZaWmaeaacqaIXaqmcqaI0aancqaI0aancqGGUaGlcqaI4aaocqaIYaGmcqaI0aancqaIXaqmaaGaeyypa0JaeGimaaJaeiOla4IaeGimaaJaeG4mamJaeGinaqJaeGinaqdaaa@51B4@ . We use the data points representing m max = 7.136, m min = 3.092, p max = 171.9 and p min = 110.3 to calculate B A = 7.136 3.092 171.9 110.3 = 0.0656 MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqWGcbGqaeaacqWGbbqqaaGccqGH9aqpjuaGdaWcaaqaaiabiEda3iabc6caUiabigdaXiabiodaZiabiAda2iabgkHiTiabiodaZiabc6caUiabicdaWiabiMda5iabikdaYaqaaiabigdaXiabiEda3iabigdaXiabc6caUiabiMda5iabgkHiTiabigdaXiabigdaXiabicdaWiabc6caUiabiodaZaaakiabg2da9iabicdaWiabc6caUiabicdaWiabiAda2iabiwda1iabiAda2aaa@4BAB@ .

We can estimate μ p as m a v e p a v e MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqWGTbqBdaWgaaqaaiabdggaHjabdAha2jabdwgaLbqabaaabaGaemiCaa3aaSbaaeaacqWGHbqycqWG2bGDcqWGLbqzaeqaaaaaaaa@37A7@ from equation (70). In computing m ave and p ave , care needs to be taken to discard the transient data and to select data points from a complete number of cycles. In our example, as the oscillatory period appears to be about 2 hours, we take the last 20 data points which form approximately 5 cycles, see Figure 4. This gives an estimate for μ p of 0.0344. We estimate the ratio of the amplitudes B A MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqWGcbGqaeaacqWGbbqqaaaaaa@2E8B@ using Bm max - m min and Ap max - p min , see Figure 4. Substituting μ p ≈ 0.0344 and B A MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqWGcbGqaeaacqWGbbqqaaaaaa@2E8B@ ≈ 0.0656 into (7) gives an estimate for ω of 0.0559. We now use the fact that A, B and ω are functions of the unknowns p0, n and μ m in order to solve numerically (7), (8), (14) and (15). This gives p0 ≈ 90, n ≈ 5 and μ p ≈ 0.026. We also require a prior for τ. Lindstedt's method tells us that the closeness in value of τ and τ cr depends on the difference between ω and the observed frequency Ω. Our observed period is about 2 hours giving an observed frequency of 0.052. Substituting our estimates into (9) and (15) we obtain τ cr ≈ 17.6567 and τ ≈ 19.9508. See Table 3 for a summary of our priors for the synthetic experiment. The variances for the priors are based, in the case of μ m and μ p , on the decay rate measurements taken by Hirata et al., and in the case of the other variables on biological assumptions combined with the mathematical analysis.

Table 3 Experiment 1 (Synthetic Data).


Figure 5 shows the prior pdfs and the corresponding posterior pdfs for each of the parameters inferred in Experiment 1. The parameter values used in the model to simulate the data are shown on each x-axis with a dot. The peak of each curve indicates which parameter values are most probable, and a sharper peak indicates more certainty. For all parameters in Experiment 1 the peak lies over the input value to within visual accuracy. We are also interested in the extent to which the posterior pdfs, which express the probability of each parameter given the data, modify our prior beliefs. Here we used Gaussian priors with mean based on the point estimates outlined above and small variance. The posterior pdfs for p0, n, μ m and μ p provide only slight modifications of the prior pdfs suggesting that the a priori mathematical analysis has allowed us to select these priors very effectively. The fact that these mathematically informed priors, which are based on broad features of the data, are reasonably close to the posterior pdfs suggests that the bifurcation analysis is relevant and consistent with the data. The parameter τ shows the most change from prior to posterior. We can show from (16) that the amplitude of mRNA is sensitive to τ, and this experiment (and Experiment 2 below) suggests that the likelihood function is also strongly affected by τ. In Figure 6 we rerun the same experiment using uniform priors over the point estimates. We see more modification of the posterior pdfs, relative to the prior pdfs but, as shown by the flatness of the pdf curves, there is increased uncertainity in the prediction of parameter values. We ran several more experiments with different priors (not shown here) and conclude that when we apply more general priors we run into trouble with this hard inference problem and get stuck in suboptimal posteriors. This highlights the importance of exploiting whatever information can be gleamed from the structure of the model.

Figure 5
figure 5

Experiment 1: Prior pdfs, posterior pdfs and parameter value used in model.

Figure 6
figure 6

Experiment 1: Binned posterior samples compared to mean ie uniform prior.

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)

Hirata et al. measured the half lives of hes1 mRNA and Hes1 protein and estimated that hes1 mRNA is degraded with a half life of 24.1 ± 1.7 minutes and that Hes1 protein is degraded with a half life of 22.3 ± 3.1 minutes. These half lifes are close but not identical. Expressed as decay rates, the rate for hes1 mRNA is 0.0288 ± 0.002 and the rate for Hes1 protein is 0.0311 ± 0.004. We used these rates to give biologically informed priors for the parameters μ m and μ p . Using Lindstedt's method we obtained an expression for the observed oscillatory frequency Ω from equation (15). As k2Δ < 0 we have ω > Ω. In other words, we have a lower bound for ω. Hirata et al. observed an oscillatory period of around 2 hours which corresponds to an oscillatory frequency (in mins) of 2 π 120 MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqaIYaGmcqaHapaCaeaacqaIXaqmcqaIYaGmcqaIWaamaaaaaa@31F2@ . Using this observation, we can argue that ω > 0.0524. Substituting ω = 0.0524, μ m = 0.0288 and μ p = 0.0311 into equation (8) gives K > 0.036. The value of K in equation (12) can be shown, via numerical computation, to be sensitive to the value of n, whereas the values of m* and p* are more sensitive to p0 than n. Given m* (or p* as m* = μ p p* from equation (21)) and the estimate for K, we can obtain point estimates from equations (21) and (12) of n ≈ 5 and p0 ≈ 100. We also require a prior for τ. Substituting our estimates for ω, μ m and μ p into equation (9) gives τ cr < 19.8 and we use this estimate of τ cr , which is close to τ, for choosing our prior for τ. We also need a scaling factor k s because the data has been scaled by a undisclosed constant, and this will also affect the initial condition. We set k s between 2 and 3, the initial condition of p to be 100 and m to be k s . See Table 4 for a summary of the priors.

Table 4 Experiment 2 (Hirata Data).


Figure 7 shows the prior pdfs and the corresponding posterior pdfs for each of the parameters inferred in Experiment 2. Once again, the customised priors are very close to the final pdfs. In this experiment, we chose a χ2 distribution rather than a normal distribution for the likelihood, in order to encourage inclusion of any outlying data points.

Figure 7
figure 7

Experiment 2: Prior and posterior pdfs.

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.


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 non-equal values for mRNA and protein decay rates and a longer time delay for transcription.


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

Following the approach of Verdugo and Rand [10], equilibrium points for the system (1) and (2) are found by setting m ˙ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafmyBa0Mbaiaaaaa@2D41@ = 0 and p ˙ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafmiCaaNbaiaaaaa@2D47@ = 0. After elimination and substitution, we obtain two expressions for p* and m*:

( p ) n + 1 + p 0 n p p 0 n μ m μ p = 0 , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeGabiWaaaqaaiabcIcaOiabdchaWnaaCaaaleqabaGaey4fIOcaaOGaeiykaKYaaWbaaSqabeaacqWGUbGBcqGHRaWkcqaIXaqmaaGccqGHRaWkcqWGWbaCdaqhaaWcbaGaeGimaadabaGaemOBa4gaaOGaemiCaa3aaWbaaSqabeaacqGHxiIkaaGccqGHsisljuaGdaWcaaqaaiabdchaWnaaDaaabaGaeGimaadabaGaemOBa4gaaaqaaiabeY7aTnaaBaaabaGaemyBa0gabeaacqaH8oqBdaWgaaqaaiabdchaWbqabaaaaaGcbaGaeyypa0dabaGaeGimaaJaeiilaWcabaaabaaabaaaaaaa@4A0B@
m = μ p p . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemyBa02aaWbaaSqabeaacqGHxiIkaaGccqGH9aqpcqaH8oqBdaWgaaWcbaGaemiCaahabeaakiabdchaWnaaCaaaleqabaGaey4fIOcaaOGaeiOla4caaa@367A@

To find out whether m* and p* are stable, we linearize about these points and define ζ and η to be deviations from the equilibrium: ζ (t) = m (t) - m*, η (t) = p (t) - p*, and η d = η (t - τ). This results in the linear system:

ζ ˙ = μ m ζ K η d , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafqOTdONbaiaacqGH9aqpcqGHsislcqaH8oqBdaWgaaWcbaGaemyBa0gabeaakiabeA7a6jabgkHiTiabdUealjabeE7aOnaaBaaaleaacqWGKbazaeqaaOGaeiilaWcaaa@3B07@
η ˙ = ζ μ p η , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafq4TdGMbaiaacqGH9aqpcqaH2oGEcqGHsislcqaH8oqBdaWgaaWcbaGaemiCaahabeaakiabeE7aOjabcYcaSaaa@3769@

where K and β are given in (12) and (13). Equations (22) and (23) can be written as the second order DDE:

η ¨ + ( μ m + μ p ) η ˙ + μ m μ p η + K η d = 0. MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafq4TdGMbamaacqGHRaWkcqGGOaakcqaH8oqBdaWgaaWcbaGaemyBa0gabeaakiabgUcaRiabeY7aTnaaBaaaleaacqWGWbaCaeqaaOGaeiykaKIafq4TdGMbaiaacqGHRaWkcqaH8oqBdaWgaaWcbaGaemyBa0gabeaakiabeY7aTnaaBaaaleaacqWGWbaCaeqaaOGaeq4TdGMaey4kaSIaem4saSKaeq4TdG2aaSbaaSqaaiabdsgaKbqabaGccqGH9aqpcqaIWaamcqGGUaGlaaa@4AE6@

We now look for a solution of the form η = eλt, which leads to

λ2 + (μ m + μ p ) λ + μ m μ p + Ke-λτ= 0.

First, we focus on the non-delay case, τ = 0. Here we have

λ = ( μ m + μ p ) ± ( μ m + μ p ) 2 4 ( μ m μ p + K ) 2 . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeq4UdWMaeyypa0tcfa4aaSaaaeaacqGHsislcqGGOaakcqaH8oqBdaWgaaqaaiabd2gaTbqabaGaey4kaSIaeqiVd02aaSbaaeaacqWGWbaCaeqaaiabcMcaPiabgglaXoaakaaabaGaeiikaGIaeqiVd02aaSbaaeaacqWGTbqBaeqaaiabgUcaRiabeY7aTnaaBaaabaGaemiCaahabeaacqGGPaqkdaahaaqabeaacqaIYaGmaaGaeyOeI0IaeGinaqJaeiikaGIaeqiVd02aaSbaaeaacqWGTbqBaeqaaiabeY7aTnaaBaaabaGaemiCaahabeaacqGHRaWkcqWGlbWscqGGPaqkaeqaaaqaaiabikdaYaaacqGGUaGlaaa@537C@

Examination of equation (26) reveals that since μ m + μ p > 0 and ( μ m + μ p ) 2 4 ( μ m μ p + K ) < μ m + μ p MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaOaaaeaacqGGOaakcqaH8oqBdaWgaaqaaiabd2gaTbqabaGaey4kaSIaeqiVd02aaSbaaeaacqWGWbaCaeqaaiabcMcaPmaaCaaabeqaaiabikdaYaaacqGHsislcqaI0aancqGGOaakcqaH8oqBdaWgaaqaaiabd2gaTbqabaGaeqiVd02aaSbaaeaacqWGWbaCaeqaaiabgUcaRiabdUealjabcMcaPaqabaGaeyipaWJaeqiVd02aaSbaaeaacqWGTbqBaeqaaiabgUcaRiabeY7aTnaaBaaabaGaemiCaahabeaaaaa@4B05@ , the real part of λ will always be negative, implying stability. For

4 (μ m μ p + K) > (μ m + μ p )2

λ has a non-zero 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 non-delay 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

Changing the first order system (22)–(23) into a second order DDE results in

η ¨ + ( μ m + μ p ) η ˙ + μ m μ p η = K η d + H 2 η d 2 + H 3 η d 3 + MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafq4TdGMbamaacqGHRaWkcqGGOaakcqaH8oqBdaWgaaWcbaGaemyBa0gabeaakiabgUcaRiabeY7aTnaaBaaaleaacqWGWbaCaeqaaOGaeiykaKIafq4TdGMbaiaacqGHRaWkcqaH8oqBdaWgaaWcbaGaemyBa0gabeaakiabeY7aTnaaBaaaleaacqWGWbaCaeqaaOGaeq4TdGMaeyypa0JaeyOeI0Iaem4saSKaeq4TdG2aaSbaaSqaaiabdsgaKbqabaGccqGHRaWkcqWGibasdaWgaaWcbaGaeGOmaidabeaakiabeE7aOnaaDaaaleaacqWGKbazaeaacqaIYaGmaaGccqGHRaWkcqWGibasdaWgaaWcbaGaeG4mamdabeaakiabeE7aOnaaDaaaleaacqWGKbazaeaacqaIZaWmaaGccqGHRaWkcqWIMaYsaaa@59B9@

where the coefficients K, H2 and H3 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.

We then have

Ω 2 d 2 u d σ 2 + ( μ m + μ p ) Ω d u d σ + μ m μ p u = K u d + ϵ H 2 u d 2 + ϵ 2 H 3 u d 3 + , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeuyQdC1aaWbaaSqabeaacqaIYaGmaaqcfa4aaSaaaeaacqWGKbazdaahaaqabeaacqaIYaGmaaGaemyDauhabaGaemizaqMaeq4Wdm3aaWbaaeqabaGaeGOmaidaaaaakiabgUcaRiabcIcaOiabeY7aTnaaBaaaleaacqWGTbqBaeqaaOGaey4kaSIaeqiVd02aaSbaaSqaaiabdchaWbqabaGccqGGPaqkcqqHPoWvjuaGdaWcaaqaaiabdsgaKjabdwha1bqaaiabdsgaKjabeo8aZbaakiabgUcaRiabeY7aTnaaBaaaleaacqWGTbqBaeqaaOGaeqiVd02aaSbaaSqaaiabdchaWbqabaGccqWG1bqDcqGH9aqpcqGHsislcqWGlbWscqWG1bqDdaWgaaWcbaGaemizaqgabeaakiabg2da9mbvLH1qn1uy0Hws0fgBPngaryWyT1wAXadaiqaacqWF1pGScqWGibasdaWgaaWcbaGaeGOmaidabeaakiabdwha1naaDaaaleaacqWGKbazaeaacqaIYaGmaaGccqGHRaWkcqWF1pGSdaahaaWcbeqaaiabikdaYaaakiabdIeainaaBaaaleaacqaIZaWmaeqaaOGaemyDau3aa0baaSqaaiabdsgaKbqaaiabiodaZaaakiabgUcaRiablAciljabcYcaSaaa@7716@

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

Ω = ω + ϵ2k2 + ...

Next we expand the delay term u d

u d = u (σ - ωτ cr ) - ϵ2 (k2 τ cr + ωδ) u' (σ - ωτ cr ) + O(ϵ3)

Finally we expand u (t) in a power series in ϵ to obtain

u (σ) = u0 (σ) + ϵu1 (σ) + ϵ2u2 (σ) + ...

Substituting and collecting terms we find

ω 2 d 2 u 0 d σ 2 + ( μ m + μ p ) ω d u 0 d σ + μ m μ p u 0 = K u 0 ( σ ω τ c r ) MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqyYdC3aaWbaaSqabeaacqaIYaGmaaqcfa4aaSaaaeaacqWGKbazdaahaaqabeaacqaIYaGmaaGaemyDau3aaSbaaeaacqaIWaamaeqaaaqaaiabdsgaKjabeo8aZnaaCaaabeqaaiabikdaYaaaaaGccqGHRaWkcqGGOaakcqaH8oqBdaWgaaWcbaGaemyBa0gabeaakiabgUcaRiabeY7aTnaaBaaaleaacqWGWbaCaeqaaOGaeiykaKIaeqyYdCxcfa4aaSaaaeaacqWGKbazcqWG1bqDdaWgaaqaaiabicdaWaqabaaabaGaemizaqMaeq4WdmhaaOGaey4kaSIaeqiVd02aaSbaaSqaaiabd2gaTbqabaGccqaH8oqBdaWgaaWcbaGaemiCaahabeaakiabdwha1naaBaaaleaacqaIWaamaeqaaOGaeyypa0JaeyOeI0Iaem4saSKaemyDau3aaSbaaSqaaiabicdaWaqabaGccqGGOaakcqaHdpWCcqGHsislcqaHjpWDcqaHepaDdaWgaaWcbaGaem4yamMaemOCaihabeaakiabcMcaPaaa@66E8@
ω 2 d 2 u 1 d σ 2 + ( μ m + μ p ) ω d u 1 d σ + μ m μ p u 1 = K u 1 ( σ ω τ c r ) + H 2 u 0 2 ( σ ω τ c r ) MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqyYdC3aaWbaaSqabeaacqaIYaGmaaqcfa4aaSaaaeaacqWGKbazdaahaaqabeaacqaIYaGmaaGaemyDau3aaSbaaeaacqaIXaqmaeqaaaqaaiabdsgaKjabeo8aZnaaCaaabeqaaiabikdaYaaaaaGccqGHRaWkcqGGOaakcqaH8oqBdaWgaaWcbaGaemyBa0gabeaakiabgUcaRiabeY7aTnaaBaaaleaacqWGWbaCaeqaaOGaeiykaKIaeqyYdCxcfa4aaSaaaeaacqWGKbazcqWG1bqDdaWgaaqaaiabigdaXaqabaaabaGaemizaqMaeq4WdmhaaOGaey4kaSIaeqiVd02aaSbaaSqaaiabd2gaTbqabaGccqaH8oqBdaWgaaWcbaGaemiCaahabeaakiabdwha1naaBaaaleaacqaIXaqmaeqaaOGaeyypa0JaeyOeI0Iaem4saSKaemyDau3aaSbaaSqaaiabigdaXaqabaGccqGGOaakcqaHdpWCcqGHsislcqaHjpWDcqaHepaDdaWgaaWcbaGaem4yamMaemOCaihabeaakiabcMcaPiabgUcaRiabdIeainaaBaaaleaacqaIYaGmaeqaaOGaemyDau3aa0baaSqaaiabicdaWaqaaiabikdaYaaakiabcIcaOiabeo8aZjabgkHiTiabeM8a3jabes8a0naaBaaaleaacqWGJbWycqWGYbGCaeqaaOGaeiykaKcaaa@7883@
ω 2 d 2 u 2 d σ 2 + ( μ m + μ p ) ω d u 2 d σ + μ m μ p u 2 = K u 2 ( σ ω τ c r ) + H 3 u 0 3 ( σ ω τ c r ) + 2 H 2 u 0 ( τ ω τ c r ) u 1 ( τ ω τ c r ) + K ( k 2 τ c r + ω δ ) u 0 ( τ ω τ c r ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeaabmWaaaqaaiabeM8a3naaCaaaleqabaGaeGOmaidaaKqbaoaalaaabaGaemizaq2aaWbaaeqabaGaeGOmaidaaiabdwha1naaBaaabaGaeGOmaidabeaaaeaacqWGKbazcqaHdpWCdaahaaqabeaacqaIYaGmaaaaaOGaey4kaSIaeiikaGIaeqiVd02aaSbaaSqaaiabd2gaTbqabaGccqGHRaWkcqaH8oqBdaWgaaWcbaGaemiCaahabeaakiabcMcaPiabeM8a3LqbaoaalaaabaGaemizaqMaemyDau3aaSbaaeaacqaIYaGmaeqaaaqaaiabdsgaKjabeo8aZbaakiabgUcaRiabeY7aTnaaBaaaleaacqWGTbqBaeqaaOGaeqiVd02aaSbaaSqaaiabdchaWbqabaGccqWG1bqDdaWgaaWcbaGaeGOmaidabeaaaOqaaiabg2da9aqaaiabgkHiTiabdUealjabdwha1naaBaaaleaacqaIYaGmaeqaaOGaeiikaGIaeq4WdmNaeyOeI0IaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGGPaqkcqGHRaWkcqWGibasdaWgaaWcbaGaeG4mamdabeaakiabdwha1naaDaaaleaacqaIWaamaeaacqaIZaWmaaGccqGGOaakcqaHdpWCcqGHsislcqaHjpWDcqaHepaDdaWgaaWcbaGaem4yamMaemOCaihabeaakiabcMcaPaqaaaqaaaqaaiabgUcaRiabikdaYiabdIeainaaBaaaleaacqaIYaGmaeqaaOGaemyDau3aaSbaaSqaaiabicdaWaqabaGccqGGOaakcqaHepaDcqGHsislcqaHjpWDcqaHepaDdaWgaaWcbaGaem4yamMaemOCaihabeaakiabcMcaPiabdwha1naaBaaaleaacqaIXaqmaeqaaOGaeiikaGIaeqiXdqNaeyOeI0IaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGGPaqkaeaaaeaaaeaacqGHRaWkcqWGlbWscqGGOaakcqWGRbWAdaWgaaWcbaGaeGOmaidabeaakiabes8a0naaBaaaleaacqWGJbWycqWGYbGCaeqaaOGaey4kaSIaeqyYdCNaeqiTdqMaeiykaKIafmyDauNbauaadaWgaaWcbaGaeGimaadabeaakiabcIcaOiabes8a0jabgkHiTiabeM8a3jabes8a0naaBaaaleaacqWGJbWycqWGYbGCaeqaaOGaeiykaKIaeiOla4caaaaa@B56E@

We take the solution of the u0 and u1 equations to have the form

u 0 ( σ ) = A ^ cos σ , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemyDau3aaSbaaSqaaiabicdaWaqabaGccqGGOaakcqaHdpWCcqGGPaqkcqGH9aqpcuWGbbqqgaqcaiGbcogaJjabc+gaVjabcohaZjabeo8aZjabcYcaSaaa@3B17@

u1(σ) = m1 sin 2σ + m2 sin 2σ + m3,

where A = A ^ ϵ MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafmyqaeKbaKaatqvzynutnfgDOLeDHXwAJbqegmwBTLwmWaaceaGae8x9dilaaa@36DD@ from (29) and (31). Next we substitute (41) and (42) into (40) to give

ω 2 ( 4 m 1 sin 2 σ 4 m 2 cos 2 σ ) + ω ( μ m + μ p ) ( 2 m 1 cos 2 σ 2 m 2 sin 2 σ ) + K m 1 ( sin 2 σ cos 2 ω τ c r cos 2 σ sin 2 ω τ c r ) + K m 2 ( cos 2 σ cos 2 ω τ c r sin 2 σ sin 2 ω τ c r ) + K m 3 + μ m μ p ( m 1 sin 2 σ + m 2 cos 2 σ ) = 1 2 H 2 A ^ 2 cos 2 σ cos 2 ω τ c r + 1 2 H 2 A ^ 2 sin 2 σ sin 2 ω τ c r + 1 2 H 2 A ^ 2 . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeqabeGaaaqaauaabiqaheaaaaqaaiabeM8a3naaCaaaleqabaGaeGOmaidaaOGaeiikaGIaeyOeI0IaeGinaqJaemyBa02aaSbaaSqaaiabigdaXaqabaGccyGGZbWCcqGGPbqAcqGGUbGBcqaIYaGmcqaHdpWCcqGHsislcqaI0aancqWGTbqBdaWgaaWcbaGaeGOmaidabeaakiGbcogaJjabc+gaVjabcohaZjabikdaYiabeo8aZjabcMcaPiabgUcaRaqaaiabeM8a3jabcIcaOiabeY7aTnaaBaaaleaacqWGTbqBaeqaaOGaey4kaSIaeqiVd02aaSbaaSqaaiabdchaWbqabaGccqGGPaqkcqGGOaakcqaIYaGmcqWGTbqBdaWgaaWcbaGaeGymaedabeaakiGbcogaJjabc+gaVjabcohaZjabikdaYiabeo8aZjabgkHiTiabikdaYiabd2gaTnaaBaaaleaacqaIYaGmaeqaaOGagi4CamNaeiyAaKMaeiOBa4MaeGOmaiJaeq4WdmNaeiykaKIaey4kaScabaGaem4saSKaemyBa02aaSbaaSqaaiabigdaXaqabaGccqGGOaakcyGGZbWCcqGGPbqAcqGGUbGBcqaIYaGmcqaHdpWCcyGGJbWycqGGVbWBcqGGZbWCcqaIYaGmcqaHjpWDcqaHepaDdaWgaaWcbaGaem4yamMaemOCaihabeaakiabgkHiTiGbcogaJjabc+gaVjabcohaZjabikdaYiabeo8aZjGbcohaZjabcMgaPjabc6gaUjabikdaYiabeM8a3jabes8a0naaBaaaleaacqWGJbWycqWGYbGCaeqaaOGaeiykaKIaey4kaScabaGaem4saSKaemyBa02aaSbaaSqaaiabikdaYaqabaGccqGGOaakcyGGJbWycqGGVbWBcqGGZbWCcqaIYaGmcqaHdpWCcyGGJbWycqGGVbWBcqGGZbWCcqaIYaGmcqaHjpWDcqaHepaDdaWgaaWcbaGaem4yamMaemOCaihabeaakiabgkHiTiGbcohaZjabcMgaPjabc6gaUjabikdaYiabeo8aZjGbcohaZjabcMgaPjabc6gaUjabikdaYiabeM8a3jabes8a0naaBaaaleaacqWGJbWycqWGYbGCaeqaaOGaeiykaKIaey4kaScabaGaem4saSKaemyBa02aaSbaaSqaaiabiodaZaqabaGccqGHRaWkcqaH8oqBdaWgaaWcbaGaemyBa0gabeaakiabeY7aTnaaBaaaleaacqWGWbaCaeqaaOGaeiikaGIaemyBa02aaSbaaSqaaiabigdaXaqabaGccyGGZbWCcqGGPbqAcqGGUbGBcqaIYaGmcqaHdpWCcqGHRaWkcqWGTbqBdaWgaaWcbaGaeGOmaidabeaakiGbcogaJjabc+gaVjabcohaZjabikdaYiabeo8aZjabcMcaPaqaaaqaaaaaaeaafaqaaeWbcaaaaeaaaeaaaeaaaeaaaeaaaeaaaeaaaeaaaeaaaeaaaeaacqGH9aqpaeaajuaGdaWcaaqaaiabigdaXaqaaiabikdaYaaakiabdIeainaaBaaaleaacqaIYaGmaeqaaOGafmyqaeKbaKaadaahaaWcbeqaaiabikdaYaaakiGbcogaJjabc+gaVjabcohaZjabikdaYiabeo8aZjGbcogaJjabc+gaVjabcohaZjabikdaYiabeM8a3jabes8a0naaBaaaleaacqWGJbWycqWGYbGCaeqaaOGaey4kaScabaaabaqcfa4aaSaaaeaacqaIXaqmaeaacqaIYaGmaaGccqWGibasdaWgaaWcbaGaeGOmaidabeaakiqbdgeabzaajaWaaWbaaSqabeaacqaIYaGmaaGccyGGZbWCcqGGPbqAcqGGUbGBcqaIYaGmcqaHdpWCcyGGZbWCcqGGPbqAcqGGUbGBcqaIYaGmcqaHjpWDcqaHepaDdaWgaaWcbaGaem4yamMaemOCaihabeaakiabgUcaRKqbaoaalaaabaGaeGymaedabaGaeGOmaidaaOGaemisaG0aaSbaaSqaaiabikdaYaqabaGccuWGbbqqgaqcamaaCaaaleqabaGaeGOmaidaaOGaeiOla4caaaaaaaa@2273@

Collecting sin 2σ terms gives

4 ω 2 m 1 2 ( μ m + μ p ) ω m 2 + K m 1 cos 2 ω τ c r + K m 2 sin 2 ω τ c r + μ m μ p m 1 = 1 2 H 2 A ^ 2 sin 2 ω τ c r . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeqabeGaaaqaauaabiqadeaaaeaacqGHsislcqaI0aancqaHjpWDdaahaaWcbeqaaiabikdaYaaakiabd2gaTnaaBaaaleaacqaIXaqmaeqaaOGaeyOeI0IaeGOmaiJaeiikaGIaeqiVd02aaSbaaSqaaiabd2gaTbqabaGccqGHRaWkcqaH8oqBdaWgaaWcbaGaemiCaahabeaakiabcMcaPiabeM8a3jabd2gaTnaaBaaaleaacqaIYaGmaeqaaOGaey4kaScabaGaem4saSKaemyBa02aaSbaaSqaaiabigdaXaqabaGccyGGJbWycqGGVbWBcqGGZbWCcqaIYaGmcqaHjpWDcqaHepaDdaWgaaWcbaGaem4yamMaemOCaihabeaakiabgUcaRiabdUealjabd2gaTnaaBaaaleaacqaIYaGmaeqaaOGagi4CamNaeiyAaKMaeiOBa4MaeGOmaiJaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGHRaWkcqaH8oqBdaWgaaWcbaGaemyBa0gabeaakiabeY7aTnaaBaaaleaacqWGWbaCaeqaaOGaemyBa02aaSbaaSqaaiabigdaXaqabaaakeaaaaaabaqbaeaabmGaaaqaaaqaaaqaaaqaaaqaaiabg2da9aqaaKqbaoaalaaabaGaeGymaedabaGaeGOmaidaaOGaemisaG0aaSbaaSqaaiabikdaYaqabaGccuWGbbqqgaqcamaaCaaaleqabaGaeGOmaidaaOGagi4CamNaeiyAaKMaeiOBa4MaeGOmaiJaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGGUaGlaaaaaaaa@8223@

Collecting cos 2σ terms gives

4 ω 2 m 2 2 ( μ m + μ p ) ω m 1 + K m 1 sin 2 ω τ c r + K m 2 cos 2 ω τ c r + μ m μ p m 2 = 1 2 H 2 A ^ 2 cos 2 ω τ c r . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeqabeGaaaqaauaabiqadeaaaeaacqGHsislcqaI0aancqaHjpWDdaahaaWcbeqaaiabikdaYaaakiabd2gaTnaaBaaaleaacqaIYaGmaeqaaOGaeyOeI0IaeGOmaiJaeiikaGIaeqiVd02aaSbaaSqaaiabd2gaTbqabaGccqGHRaWkcqaH8oqBdaWgaaWcbaGaemiCaahabeaakiabcMcaPiabeM8a3jabd2gaTnaaBaaaleaacqaIXaqmaeqaaOGaey4kaScabaGaem4saSKaemyBa02aaSbaaSqaaiabigdaXaqabaGccqGHsislcyGGZbWCcqGGPbqAcqGGUbGBcqaIYaGmcqaHjpWDcqaHepaDdaWgaaWcbaGaem4yamMaemOCaihabeaakiabgUcaRiabdUealjabd2gaTnaaBaaaleaacqaIYaGmaeqaaOGagi4yamMaei4Ba8Maei4CamNaeGOmaiJaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGHRaWkcqaH8oqBdaWgaaWcbaGaemyBa0gabeaakiabeY7aTnaaBaaaleaacqWGWbaCaeqaaOGaemyBa02aaSbaaSqaaiabikdaYaqabaaakeaaaaaabaqbaeaabmGaaaqaaaqaaaqaaaqaaaqaaiabg2da9aqaaKqbaoaalaaabaGaeGymaedabaGaeGOmaidaaOGaemisaG0aaSbaaSqaaiabikdaYaqabaGccuWGbbqqgaqcamaaCaaaleqabaGaeGOmaidaaOGagi4yamMaei4Ba8Maei4CamNaeGOmaiJaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGGUaGlaaaaaaaa@8308@

Collecting other terms gives

K m 3 + μ m μ p m 3 = 1 2 H 2 A ^ 2 . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeaabiWaaaqaaiabdUealjabd2gaTnaaBaaaleaacqaIZaWmaeqaaOGaey4kaSIaeqiVd02aaSbaaSqaaiabd2gaTbqabaGccqaH8oqBdaWgaaWcbaGaemiCaahabeaakiabd2gaTnaaBaaaleaacqaIZaWmaeqaaaGcbaaabaaabaaabaGaeyypa0dabaqcfa4aaSaaaeaacqaIXaqmaeaacqaIYaGmaaGccqWGibasdaWgaaWcbaGaeGOmaidabeaakiqbdgeabzaajaWaaWbaaSqabeaacqaIYaGmaaGccqGGUaGlaaaaaa@42ED@


α = -4ω2 + K cos2ωτ cr + μ m μ p ,

β = -2ω(μ m + μ p ) + K sin 2ωτ cr ,

H ^ 1 = 1 2 H 2 A ^ 2 sin 2 ω τ c r , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafmisaGKbaKaadaWgaaWcbaGaeGymaedabeaakiabg2da9KqbaoaalaaabaGaeGymaedabaGaeGOmaidaaOGaemisaG0aaSbaaSqaaiabikdaYaqabaGccuWGbbqqgaqcamaaCaaaleqabaGaeGOmaidaaOGagi4CamNaeiyAaKMaeiOBa4MaeGOmaiJaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGGSaalaaa@430B@
H ^ 2 = 1 2 H 2 A ^ 2 cos 2 ω τ c r , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafmisaGKbaKaadaWgaaWcbaGaeGOmaidabeaakiabg2da9KqbaoaalaaabaGaeGymaedabaGaeGOmaidaaOGaemisaG0aaSbaaSqaaiabikdaYaqabaGccuWGbbqqgaqcamaaCaaaleqabaGaeGOmaidaaOGagi4yamMaei4Ba8Maei4CamNaeGOmaiJaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGGSaalaaa@4303@
H ^ 3 = 1 2 H 2 A ^ 2 , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafmisaGKbaKaadaWgaaWcbaGaeG4mamdabeaakiabg2da9KqbaoaalaaabaGaeGymaedabaGaeGOmaidaaOGaemisaG0aaSbaaSqaaiabikdaYaqabaGccuWGbbqqgaqcamaaCaaaleqabaGaeGOmaidaaOGaeiilaWcaaa@376B@

we find that

m 1 = α H ^ 1 β H ^ 2 α 2 + β 2 , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemyBa02aaSbaaSqaaiabigdaXaqabaGccqGH9aqpjuaGdaWcaaqaaiabeg7aHjqbdIeaizaajaWaaSbaaeaacqaIXaqmaeqaaiabgkHiTiabek7aIjqbdIeaizaajaWaaSbaaeaacqaIYaGmaeqaaaqaaiabeg7aHnaaCaaabeqaaiabikdaYaaacqGHRaWkcqaHYoGydaahaaqabeaacqaIYaGmaaaaaOGaeiilaWcaaa@4027@
m 2 = α H ^ 1 + β H ^ 2 α 2 β 2 , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemyBa02aaSbaaSqaaiabikdaYaqabaGccqGH9aqpjuaGdaWcaaqaaiabeg7aHjqbdIeaizaajaWaaSbaaeaacqaIXaqmaeqaaiabgUcaRiabek7aIjqbdIeaizaajaWaaSbaaeaacqaIYaGmaeqaaaqaaiabeg7aHnaaCaaabeqaaiabikdaYaaacqGHsislcqaHYoGydaahaaqabeaacqaIYaGmaaaaaOGaeiilaWcaaa@4029@
m 3 = H ^ 3 K + μ m μ p . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemyBa02aaSbaaSqaaiabiodaZaqabaGccqGH9aqpjuaGdaWcaaqaaiqbdIeaizaajaWaaSbaaeaacqaIZaWmaeqaaaqaaiabdUealjabgUcaRiabeY7aTnaaBaaabaGaemyBa0gabeaacqaH8oqBdaWgaaqaaiabdchaWbqabaaaaOGaeiOla4caaa@3BFB@

Now we substitute (41) and (42) into (38) and equate to zero the coefficients of the resonant terms sin σ and cos σ. This gives

k 2 ( μ m μ p + K τ c r cos ω τ c r ) + ω δ K cos ω τ c r = 3 4 H 3 A ^ 2 sin ω τ c r + H 2 2 A ^ 2 ( m 1 cos ω τ c r + m 2 sin ω τ c r + 2 m 3 sin ω τ c r ) MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeqabeGaaaqaauaabiqadeaaaeaacqWGRbWAdaWgaaWcbaGaeGOmaidabeaakiabcIcaOiabgkHiTiabeY7aTnaaBaaaleaacqWGTbqBaeqaaOGaeyOeI0IaeqiVd02aaSbaaSqaaiabdchaWbqabaGccqGHRaWkcqWGlbWscqaHepaDdaWgaaWcbaGaem4yamMaemOCaihabeaakiGbcogaJjabc+gaVjabcohaZjabeM8a3jabes8a0naaBaaaleaacqWGJbWycqWGYbGCaeqaaOGaeiykaKcabaGaey4kaSIaeqyYdCNaeqiTdqMaem4saSKagi4yamMaei4Ba8Maei4CamNaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaaakeaaaaaabaqbaeaabmGaaaqaaaqaaaqaaiabg2da9aqaaKqbaoaalaaabaGaeG4mamdabaGaeGinaqdaaOGaemisaG0aaSbaaSqaaiabiodaZaqabaGccuWGbbqqgaqcamaaCaaaleqabaGaeGOmaidaaOGagi4CamNaeiyAaKMaeiOBa4MaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGHRaWkaeaaaeaacqWGibasdaqhaaWcbaGaeGOmaidabaGaeGOmaidaaOGafmyqaeKbaKaadaahaaWcbeqaaiabikdaYaaakiabcIcaOiabd2gaTnaaBaaaleaacqaIXaqmaeqaaOGagi4yamMaei4Ba8Maei4CamNaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGHRaWkcqWGTbqBdaWgaaWcbaGaeGOmaidabeaakiGbcohaZjabcMgaPjabc6gaUjabeM8a3jabes8a0naaBaaaleaacqWGJbWycqWGYbGCaeqaaOGaey4kaSIaeGOmaiJaemyBa02aaSbaaSqaaiabiodaZaqabaGccyGGZbWCcqGGPbqAcqGGUbGBcqaHjpWDcqaHepaDdaWgaaWcbaGaem4yamMaemOCaihabeaakiabcMcaPaaaaaaaaa@9FC4@


k 2 ( 2 ω K τ c r sin ω τ c r ) ω δ K sin ω τ c r = 3 4 H 3 A ^ 2 cos ω τ c r + H 2 2 A ^ 2 ( m 1 sin ω τ c r + m 2 cos ω τ c r + 2 m 3 cos ω τ c r ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeqabeGaaaqaauaabiqadeaaaeaacqWGRbWAdaWgaaWcbaGaeGOmaidabeaakiabcIcaOiabgkHiTiabikdaYiabeM8a3jabgkHiTiabdUealjabes8a0naaBaaaleaacqWGJbWycqWGYbGCaeqaaOGagi4CamNaeiyAaKMaeiOBa4MaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGGPaqkaeaacqGHsislcqaHjpWDcqaH0oazcqWGlbWscyGGZbWCcqGGPbqAcqGGUbGBcqaHjpWDcqaHepaDdaWgaaWcbaGaem4yamMaemOCaihabeaaaOqaaaaaaeaafaqaaeWacaaabaaabaaabaGaeyypa0dabaqcfa4aaSaaaeaacqaIZaWmaeaacqaI0aanaaGccqWGibasdaWgaaWcbaGaeG4mamdabeaakiqbdgeabzaajaWaaWbaaSqabeaacqaIYaGmaaGccyGGJbWycqGGVbWBcqGGZbWCcqaHjpWDcqaHepaDdaWgaaWcbaGaem4yamMaemOCaihabeaakiabgUcaRaqaaaqaaiabdIeainaaDaaaleaacqaIYaGmaeaacqaIYaGmaaGccuWGbbqqgaqcamaaCaaaleqabaGaeGOmaidaaOGaeiikaGIaeyOeI0IaemyBa02aaSbaaSqaaiabigdaXaqabaGccyGGZbWCcqGGPbqAcqGGUbGBcqaHjpWDcqaHepaDdaWgaaWcbaGaem4yamMaemOCaihabeaakiabgUcaRiabd2gaTnaaBaaaleaacqaIYaGmaeqaaOGagi4yamMaei4Ba8Maei4CamNaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGHRaWkcqaIYaGmcqWGTbqBdaWgaaWcbaGaeG4mamdabeaakiGbcogaJjabc+gaVjabcohaZjabeM8a3jabes8a0naaBaaaleaacqWGJbWycqWGYbGCaeqaaOGaeiykaKIaeiOla4caaaaaaaa@9CD9@


c1 = -μ m - μ p + cr cosωτ cr ,

c2 = -2ω - cr sin ωτ cr ,

d 1 = 3 4 H 3 sin ω τ c r + H 2 2 ( m 1 cos ω τ c r + m 2 sin ω τ c r + 2 m 3 sin ω τ c r ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemizaq2aaSbaaSqaaiabigdaXaqabaGccqGH9aqpjuaGdaWcaaqaaiabiodaZaqaaiabisda0aaakiabdIeainaaBaaaleaacqaIZaWmaeqaaOGagi4CamNaeiyAaKMaeiOBa4MaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGHRaWkcqWGibasdaqhaaWcbaGaeGOmaidabaGaeGOmaidaaOGaeiikaGIaemyBa02aaSbaaSqaaiabigdaXaqabaGccyGGJbWycqGGVbWBcqGGZbWCcqaHjpWDcqaHepaDdaWgaaWcbaGaem4yamMaemOCaihabeaakiabgUcaRiabd2gaTnaaBaaaleaacqaIYaGmaeqaaOGagi4CamNaeiyAaKMaeiOBa4MaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGHRaWkcqaIYaGmcqWGTbqBdaWgaaWcbaGaeG4mamdabeaakiGbcohaZjabcMgaPjabc6gaUjabeM8a3jabes8a0naaBaaaleaacqWGJbWycqWGYbGCaeqaaOGaeiykaKIaeiilaWcaaa@7032@
d 2 = 3 4 H 3 cos ω τ c r + H 2 2 ( m 1 sin ω τ c r + m 2 cos ω τ c r + 2 m 3 cos ω τ c r ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemizaq2aaSbaaSqaaiabikdaYaqabaGccqGH9aqpjuaGdaWcaaqaaiabiodaZaqaaiabisda0aaakiabdIeainaaBaaaleaacqaIZaWmaeqaaOGagi4yamMaei4Ba8Maei4CamNaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGHRaWkcqWGibasdaqhaaWcbaGaeGOmaidabaGaeGOmaidaaOGaeiikaGIaeyOeI0IaemyBa02aaSbaaSqaaiabigdaXaqabaGccyGGZbWCcqGGPbqAcqGGUbGBcqaHjpWDcqaHepaDdaWgaaWcbaGaem4yamMaemOCaihabeaakiabgUcaRiabd2gaTnaaBaaaleaacqaIYaGmaeqaaOGagi4yamMaei4Ba8Maei4CamNaeqyYdCNaeqiXdq3aaSbaaSqaaiabdogaJjabdkhaYbqabaGccqGHRaWkcqaIYaGmcqWGTbqBdaWgaaWcbaGaeG4mamdabeaakiGbcogaJjabc+gaVjabcohaZjabeM8a3jabes8a0naaBaaaleaacqWGJbWycqWGYbGCaeqaaOGaeiykaKIaeiOla4caaa@7111@

Then we have

A 2 = ( c 2 cos ω τ c r + c 1 sin ω τ c r ) d 1 c 2 c 1 d 2 ω Δ K , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemyqae0aaWbaaSqabeaacqaIYaGmaaGccqGH9aqpjuaGdaWcaaqaaiabcIcaOiabdogaJnaaBaaabaGaeGOmaidabeaacyGGJbWycqGGVbWBcqGGZbWCcqaHjpWDcqaHepaDdaWgaaqaaiabdogaJjabdkhaYbqabaGaey4kaSIaem4yam2aaSbaaeaacqaIXaqmaeqaaiGbcohaZjabcMgaPjabc6gaUjabeM8a3jabes8a0naaBaaabaGaem4yamMaemOCaihabeaacqGGPaqkaeaacqWGKbazdaWgaaqaaiabigdaXaqabaGaem4yam2aaSbaaeaacqaIYaGmaeqaaiabgkHiTiabdogaJnaaBaaabaGaeGymaedabeaacqWGKbazdaWgaaqaaiabikdaYaqabaaaaOGaeqyYdCNaeuiLdqKaem4saSKaeiilaWcaaa@5C32@
k 2 = d 2 cos ω τ c r + d 1 sin ω τ c r c 1 d 2 c 2 d 1 ω δ K , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaem4AaS2aaSbaaSqaaiabikdaYaqabaGccqGH9aqpcqGHsisljuaGdaWcaaqaaiabdsgaKnaaBaaabaGaeGOmaidabeaacyGGJbWycqGGVbWBcqGGZbWCcqaHjpWDcqaHepaDdaWgaaqaaiabdogaJjabdkhaYbqabaGaey4kaSIaemizaq2aaSbaaeaacqaIXaqmaeqaaiGbcohaZjabcMgaPjabc6gaUjabeM8a3jabes8a0naaBaaabaGaem4yamMaemOCaihabeaaaeaacqWGJbWydaWgaaqaaiabigdaXaqabaGaemizaq2aaSbaaeaacqaIYaGmaeqaaiabgkHiTiabdogaJnaaBaaabaGaeGOmaidabeaacqWGKbazdaWgaaqaaiabigdaXaqabaaaaOGaeqyYdCNaeqiTdqMaem4saSKaeiilaWcaaa@5C03@

Multiplying (56) by ϵ2 and substituting back into (35) gives

Ω = ω + k2Δ,

which has the form of equation (15) (noting that k2 is negative).

Maximum ω

We have seen that K takes the same value for any μ m μ p = c where c is a constant. We are interested in investigating ω along μ m μ p = c. We will show that ω has a unique maximum along this line and the maximum value of ω occurs when μ m = μ p = c MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaOaaaeaacqWGJbWyaSqabaaaaa@2D3F@ . Rearranging equation (8) gives

ω 4 + ( μ m 2 + μ p 2 ) ω 2 + μ m 2 μ p 2 K 2 = 0. MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqyYdC3aaWbaaSqabeaacqaI0aanaaGccqGHRaWkcqGGOaakcqaH8oqBdaqhaaWcbaGaemyBa0gabaGaeGOmaidaaOGaey4kaSIaeqiVd02aa0baaSqaaiabdchaWbqaaiabikdaYaaakiabcMcaPiabeM8a3naaCaaaleqabaGaeGOmaidaaOGaey4kaSIaeqiVd02aa0baaSqaaiabd2gaTbqaaiabikdaYaaakiabeY7aTnaaDaaaleaacqWGWbaCaeaacqaIYaGmaaGccqGHsislcqWGlbWsdaahaaWcbeqaaiabikdaYaaakiabg2da9iabicdaWiabc6caUaaa@4D8C@

Substituting μ m μ p = c and K = c K into the above equation gives:

ω 4 + ω 2 ( c μ p 2 + μ p 2 ) ω 2 + c 2 c K 2 = 0. MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqyYdC3aaWbaaSqabeaacqaI0aanaaGccqGHRaWkcqaHjpWDdaahaaWcbeqaaiabikdaYaaakmaabmaabaqcfa4aaSaaaeaacqWGJbWyaeaacqaH8oqBdaqhaaqaaiabdchaWbqaaiabikdaYaaaaaGccqGHRaWkcqaH8oqBdaqhaaWcbaGaemiCaahabaGaeGOmaidaaaGccaGLOaGaayzkaaGaeqyYdC3aaWbaaSqabeaacqaIYaGmaaGccqGHRaWkcqWGJbWydaahaaWcbeqaaiabikdaYaaakiabgkHiTiabdogaJnaaDaaaleaacqWGlbWsaeaacqaIYaGmaaGccqGH9aqpcqaIWaamcqGGUaGlaaa@4D7E@

Differentiating w.r.t. μ p gives:

4 ω 3 d ω d μ p + 2 ω d ω d μ p ( c 2 μ p 2 + μ p 2 ) + ω 2 ( 2 c 2 μ p 3 + 2 μ p ) = 0 , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeGinaqJaeqyYdC3aaWbaaSqabeaacqaIZaWmaaqcfa4aaSaaaeaacqWGKbazcqaHjpWDaeaacqWGKbazcqaH8oqBdaWgaaqaaiabdchaWbqabaaaaOGaey4kaSIaeGOmaiJaeqyYdCxcfa4aaSaaaeaacqWGKbazcqaHjpWDaeaacqWGKbazcqaH8oqBdaWgaaqaaiabdchaWbqabaaaaOWaaeWaaeaajuaGdaWcaaqaaiabdogaJnaaCaaabeqaaiabikdaYaaaaeaacqaH8oqBdaqhaaqaaiabdchaWbqaaiabikdaYaaaaaGccqGHRaWkcqaH8oqBdaqhaaWcbaGaemiCaahabaGaeGOmaidaaaGccaGLOaGaayzkaaGaey4kaSIaeqyYdC3aaWbaaSqabeaacqaIYaGmaaGcdaqadaqaaKqbakabgkHiTmaalaaabaGaeGOmaiJaem4yam2aaWbaaeqabaGaeGOmaidaaaqaaiabeY7aTnaaDaaabaGaemiCaahabaGaeG4mamdaaaaakiabgUcaRiabikdaYiabeY7aTnaaBaaaleaacqWGWbaCaeqaaaGccaGLOaGaayzkaaGaeyypa0JaeGimaaJaeiilaWcaaa@68C4@

So d ω d μ p MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqWGKbazcqaHjpWDaeaacqWGKbazcqaH8oqBdaWgaaqaaiabdchaWbqabaaaaaaa@3422@ = 0 when

ω 2 ( 2 c 2 μ p 3 + 2 μ p ) = 0 2 c 2 μ p 3 + 2 μ p = 0 MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeGabiqaaaqaaiabeM8a3naaCaaaleqabaGaeGOmaidaaOWaaeWaaeaajuaGdaWcaaqaaiabgkHiTiabikdaYiabdogaJnaaCaaabeqaaiabikdaYaaaaeaacqaH8oqBdaqhaaqaaiabdchaWbqaaiabiodaZaaaaaGccqGHRaWkcqaIYaGmcqaH8oqBdaWgaaWcbaGaemiCaahabeaaaOGaayjkaiaawMcaaiabg2da9iabicdaWaqaaKqbaoaalaaabaGaeyOeI0IaeGOmaiJaem4yam2aaWbaaeqabaGaeGOmaidaaaqaaiabeY7aTnaaDaaabaGaemiCaahabaGaeG4mamdaaaaakiabgUcaRiabikdaYiabeY7aTnaaBaaaleaacqWGWbaCaeqaaOGaeyypa0JaeGimaadaaaaa@5128@


μ = c . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqiVd0Maeyypa0ZaaOaaaeaacqWGJbWyaSqabaGccqGGUaGlaaa@3137@

As μ m μ p = c, we have μ m = μ p = c MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaOaaaeaacqWGJbWyaSqabaaaaa@2D3F@ . The maximum value of ω therefore occurs when μ m = μ p and has the form

ω = K μ 2 . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqyYdCNaeyypa0ZaaOaaaeaacqWGlbWscqGHsislcqaH8oqBdaahaaWcbeqaaiabikdaYaaaaeqaaOGaeiOla4caaa@34D5@

Writing K - μ2 in terms of β gives

K μ 2 = n β ( 1 + β ) p 0 β 1 / n ( 1 + β ) 2 . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaem4saSKaeyOeI0IaeqiVd02aaWbaaSqabeaacqaIYaGmaaGccqGH9aqpjuaGdaWcaaqaaiabd6gaUjabek7aIjabgkHiTiabcIcaOiabigdaXiabgUcaRiabek7aIjabcMcaPaqaaiabdchaWnaaBaaabaGaeGimaadabeaacqaHYoGydaahaaqabeaacqaIXaqmcqGGVaWlcqWGUbGBaaGaeiikaGIaeGymaeJaey4kaSIaeqOSdiMaeiykaKYaaWbaaeqabaGaeGOmaidaaaaakiabc6caUaaa@4A67@

Differentiating with respect to β gives

β 1 / n ( 1 + β ) 2 ( n 1 ) ( n β β 1 ) ( 1 n β 1 / n 1 ( 1 + β ) 2 + 2 β 1 / n ( 1 + β ) ) p 0 ( β 1 / n ( 1 + β ) 2 ) 2 . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqaHYoGydaahaaqabeaacqaIXaqmcqGGVaWlcqWGUbGBaaGaeiikaGIaeGymaeJaey4kaSIaeqOSdiMaeiykaKYaaWbaaeqabaGaeGOmaidaaiabcIcaOiabd6gaUjabgkHiTiabigdaXiabcMcaPiabgkHiTiabcIcaOiabd6gaUjabek7aIjabgkHiTiabek7aIjabgkHiTiabigdaXiabcMcaPmaabmaabaWaaSaaaeaacqaIXaqmaeaacqWGUbGBaaGaeqOSdi2aaWbaaeqabaGaeGymaeJaei4la8IaemOBa4MaeyOeI0IaeGymaedaaiabcIcaOiabigdaXiabgUcaRiabek7aIjabcMcaPmaaCaaabeqaaiabikdaYaaacqGHRaWkcqaIYaGmcqaHYoGydaahaaqabeaacqaIXaqmcqGGVaWlcqWGUbGBaaGaeiikaGIaeGymaeJaey4kaSIaeqOSdiMaeiykaKcacaGLOaGaayzkaaaabaGaemiCaa3aaSbaaeaacqaIWaamaeqaamaabmaabaGaeqOSdi2aaWbaaeqabaGaeGymaeJaei4la8IaemOBa4gaaiabcIcaOiabigdaXiabgUcaRiabek7aIjabcMcaPmaaCaaabeqaaiabikdaYaaaaiaawIcacaGLPaaadaahaaqabeaacqaIYaGmaaaaaiabc6caUaaa@7525@

Setting this expression to zero and dividing by β1/n(1 + β)/p0 gives

n + n β β 1 ( n 1 ) ( 1 / n ( 1 + β ) + 2 β ) + 1 n β ( 1 + β ) + 2 = 0 , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOBa4Maey4kaSIaemOBa4MaeqOSdiMaeyOeI0IaeqOSdiMaeyOeI0IaeGymaeJaeyOeI0IaeiikaGIaemOBa4MaeyOeI0IaeGymaeJaeiykaKIaeiikaGIaeGymaeJaei4la8IaemOBa4MaeiikaGIaeGymaeJaey4kaSIaeqOSdiMaeiykaKIaey4kaSIaeGOmaiJaeqOSdiMaeiykaKscfa4aaSaaaeaacqaIXaqmaeaacqWGUbGBcqaHYoGyaaGccqGGOaakcqaIXaqmcqGHRaWkcqaHYoGycqGGPaqkcqGHRaWkcqaIYaGmcqGH9aqpcqaIWaamcqGGSaalaaa@57AE@

which simplifies to

(n2 - 1) β2 - (n2 + 2) β - 1 = 0.

Can the data tell us anything about the parameters directly?

Integrating equation (2) gives

0 T p ˙ ( t ) d t = 0 T m ( t ) d t μ p 0 T p ( t ) d t . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaa8qmaeaacuWGWbaCgaGaaiabcIcaOiabdsha0jabcMcaPiabdsgaKjabdsha0bWcbaGaeGimaadabaGaemivaqfaniabgUIiYdGccqGH9aqpdaWdXaqaaiabd2gaTjabcIcaOiabdsha0jabcMcaPiabdsgaKjabdsha0bWcbaGaeGimaadabaGaemivaqfaniabgUIiYdGccqGHsislcqaH8oqBdaWgaaWcbaGaemiCaahabeaakmaapedabaGaemiCaaNaeiikaGIaemiDaqNaeiykaKIaemizaqMaemiDaqhaleaacqaIWaamaeaacqWGubava0Gaey4kIipakiabc6caUaaa@5532@

From the Fundamental Theorem of Calculus, it follows that

p ( T ) p ( 0 ) = 0 T m ( t ) d t μ p 0 T p ( t ) d t . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiCaaNaeiikaGIaemivaqLaeiykaKIaeyOeI0IaemiCaaNaeiikaGIaeGimaaJaeiykaKIaeyypa0Zaa8qmaeaacqWGTbqBcqGGOaakcqWG0baDcqGGPaqkcqWGKbazcqWG0baDaSqaaiabicdaWaqaaiabdsfaubqdcqGHRiI8aOGaeyOeI0IaeqiVd02aaSbaaSqaaiabdchaWbqabaGcdaWdXaqaaiabdchaWjabcIcaOiabdsha0jabcMcaPiabdsgaKjabdsha0bWcbaGaeGimaadabaGaemivaqfaniabgUIiYdGccqGGUaGlaaa@52CB@

If p is periodic and T is a multiple of the period then the left-hand side of this equation vanishes and the right-hand side may be written m ave - μ p p ave . This gives us the result

μ p = m a v e p a v e . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqiVd02aaSbaaSqaaiabdchaWbqabaGccqGH9aqpjuaGdaWcaaqaaiabd2gaTnaaBaaabaGaemyyaeMaemODayNaemyzaugabeaaaeaacqWGWbaCdaWgaaqaaiabdggaHjabdAha2jabdwgaLbqabaaaaOGaeiOla4caaa@3D3E@

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

Our approach was to choose candidate values for our unknown parameters, θ j MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqiUde3aa0baaSqaaiabdQgaQbqaaiabgEHiQaaaaaa@3004@ , (j = 1 ... 5 in Experiment 1 and j = 1 ... 6 in Experiment 2) in logspace, based on the current value of θ j ( i ) MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqiUde3aa0baaSqaaiabdQgaQbqaaiabcIcaOiabdMgaPjabcMcaPaaaaaa@3222@ , using a Gaussian proposal distribution, q ( θ j | θ j ( i ) ) MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemyCae3aaeWaaeaacqaH4oqCdaqhaaWcbaGaemOAaOgabaGaey4fIOcaaOGaeiiFaWNaeqiUde3aa0baaSqaaiabdQgaQbqaaiabcIcaOiabdMgaPjabcMcaPaaaaOGaayjkaiaawMcaaaaa@3AD9@ given by

q ( θ j | θ j ( i ) ) = 1 σ j 2 π exp ( ( θ j θ j ( i ) ) 2 2 σ j 2 ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemyCae3aaeWaaeaacqaH4oqCdaqhaaWcbaGaemOAaOgabaGaey4fIOcaaOGaeiiFaWNaeqiUde3aa0baaSqaaiabdQgaQbqaaiabcIcaOiabdMgaPjabcMcaPaaaaOGaayjkaiaawMcaaiabg2da9KqbaoaalaaabaGaeGymaedabaGaeq4Wdm3aaSbaaeaacqWGQbGAaeqaamaakaaabaGaeGOmaiJaeqiWdahabeaaaaGccyGGLbqzcqGG4baEcqGGWbaCdaqadaqaaiabgkHiTKqbaoaalaaabaWaaeWaaeaacqaH4oqCdaqhaaqaaiabdQgaQbqaaiabgEHiQaaacqGHsislcqaH4oqCdaqhaaqaaiabdQgaQbqaaiabcIcaOiabdMgaPjabcMcaPaaaaiaawIcacaGLPaaadaahaaqabeaacqaIYaGmaaaabaGaeGOmaiJaeq4Wdm3aa0baaeaacqWGQbGAaeaacqaIYaGmaaaaaaGccaGLOaGaayzkaaGaeiOla4caaa@5F11@

σ j values for the proposal rates are chosen in order that the acceptance rates are between 20% and 35%.

The candidate values are accepted or rejected with an acceptance probability, A MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWefv3ySLgznfgDOfdaryqr1ngBPrginfgDObYtUvgaiqaacqWFaeFqaaa@3715@ (θ(i), θ *), given by

A ( θ j ( i ) , θ j ) = min { 1 , p ^ ( θ j ) q ( θ j ( i ) | θ j ) p ^ ( θ j ( i ) ) q ( θ j | θ j ( i ) ) } , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWefv3ySLgznfgDOfdaryqr1ngBPrginfgDObYtUvgaiqaacqWFaeFqdaqadaqaaiabeI7aXnaaDaaaleaacqWGQbGAaeaacqGGOaakcqWGPbqAcqGGPaqkaaGccqGGSaalcqaH4oqCdaqhaaWcbaGaemOAaOgabaGaey4fIOcaaaGccaGLOaGaayzkaaGaeyypa0JagiyBa0MaeiyAaKMaeiOBa42aaiWaaeaacqaIXaqmcqGGSaaljuaGdaWcaaqaaiqbdchaWzaajaWaaeWaaeaacqaH4oqCdaqhaaqaaiabdQgaQbqaaiabgEHiQaaaaiaawIcacaGLPaaacqWGXbqCdaqadaqaaiabeI7aXnaaDaaabaGaemOAaOgabaGaeiikaGIaemyAaKMaeiykaKcaaiabcYha8jabeI7aXnaaDaaabaGaemOAaOgabaGaey4fIOcaaaGaayjkaiaawMcaaaqaaiqbdchaWzaajaWaaeWaaeaacqaH4oqCdaqhaaqaaiabdQgaQbqaaiabcIcaOiabdMgaPjabcMcaPaaaaiaawIcacaGLPaaacqWGXbqCdaqadaqaaiabeI7aXnaaDaaabaGaemOAaOgabaGaey4fIOcaaiabcYha8jabeI7aXnaaDaaabaGaemOAaOgabaGaeiikaGIaemyAaKMaeiykaKcaaaGaayjkaiaawMcaaaaaaOGaay5Eaiaaw2haaiabcYcaSaaa@7D2B@


p ^ = p ( D | θ j ) p ( θ j ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafmiCaaNbaKaacqGH9aqpcqWGWbaCcqGGOaaktuuDJXwAK1uy0HwmaeHbfv3ySLgzG0uy0Hgip5wzaGabaiab=nq8ejabcYha8jabeI7aXnaaBaaaleaacqWGQbGAaeqaaOGaeiykaKIaemiCaaNaeiikaGIaeqiUde3aaSbaaSqaaiabdQgaQbqabaGccqGGPaqkcqGGSaalaaa@4910@
p ( D | θ j ) = p ( { x k } | θ j , σ , I ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiCaaNaeiikaGYefv3ySLgznfgDOfdaryqr1ngBPrginfgDObYtUvgaiqaacqWFdeprcqGG8baFcqaH4oqCdaWgaaWcbaGaemOAaOgabeaakiabcMcaPiabg2da9iabdchaWjabcIcaOiabcUha7jabdIha4naaBaaaleaacqWGRbWAaeqaaOGaeiyFa0NaeiiFaWNaeqiUde3aaSbaaSqaaiabdQgaQbqabaGccqGGSaalcqaHdpWCcqGGSaalcqWGjbqscqGGPaqkcqGGSaalaaa@53C3@


p ( θ j ) = N ( μ j , σ p j ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiCaaNaeiikaGIaeqiUde3aaSbaaSqaaiabdQgaQbqabaGccqGGPaqkcqGH9aqptuuDJXwAK1uy0HwmaeHbfv3ySLgzG0uy0Hgip5wzaGabaiab=1q8ojabcIcaOiabeY7aTnaaBaaaleaacqWGQbGAaeqaaOGaeiilaWIaeq4Wdm3aaSbaaSqaaiabdchaWjabdQgaQbqabaGccqGGPaqkcqGGUaGlaaa@4A65@

A pseudo-code outline for the MCMC Metropolis Hastings algorithm is as follows:

step 1 Initialise θ j ( 0 ) MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqiUde3aa0baaSqaaiabdQgaQbqaaiabcIcaOiabicdaWiabcMcaPaaaaaa@31B5@

step 2 For i = 0 to N:

Sample r ~ N MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWefv3ySLgznfgDOfdaryqr1ngBPrginfgDObYtUvgaiqaacqWFneVtaaa@372F@ (0, 1)

Sample θ j ~ q ( θ j | θ j ( i ) ) MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqiUde3aa0baaSqaaiabdQgaQbqaaiabgEHiQaaakiabc6ha+jabdghaXnaabmaabaGaeqiUde3aa0baaSqaaiabdQgaQbqaaiabgEHiQaaakiabcYha8jabeI7aXnaaDaaaleaacqWGQbGAaeaacqGGOaakcqWGPbqAcqGGPaqkaaaakiaawIcacaGLPaaaaaa@4096@ (choosing j at random

If r < A ( θ j ( i ) , θ j ) = min { 1 , p ^ ( θ j ) q ( θ j ( i ) | θ j ) p ^ ( θ j ( i ) ) q ( θ j | θ j ( i ) ) } MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOCaiNaeyipaWdcbaGae8xqaeKae8NCai3aaeWaaeaacqaH4oqCdaqhaaWcbaGaemOAaOgabaGaeiikaGIaemyAaKMaeiykaKcaaOGaeiilaWIaeqiUde3aa0baaSqaaiabdQgaQbqaaiabgEHiQaaaaOGaayjkaiaawMcaaiabg2da9iGbc2gaTjabcMgaPjabc6gaUnaacmaabaGaeGymaeJaeiilaWscfa4aaSaaaeaacuWGWbaCgaqcamaabmaabaGaeqiUde3aa0baaeaacqWGQbGAaeaacqGHxiIkaaaacaGLOaGaayzkaaGaemyCae3aaeWaaeaacqaH4oqCdaqhaaqaaiabdQgaQbqaaiabcIcaOiabdMgaPjabcMcaPaaacqGG8baFcqaH4oqCdaqhaaqaaiabdQgaQbqaaiabgEHiQaaaaiaawIcacaGLPaaaaeaacuWGWbaCgaqcamaabmaabaGaeqiUde3aa0baaeaacqWGQbGAaeaacqGGOaakcqWGPbqAcqGGPaqkaaaacaGLOaGaayzkaaGaemyCae3aaeWaaeaacqaH4oqCdaqhaaqaaiabdQgaQbqaaiabgEHiQaaacqGG8baFcqaH4oqCdaqhaaqaaiabdQgaQbqaaiabcIcaOiabdMgaPjabcMcaPaaaaiaawIcacaGLPaaaaaaakiaawUhacaGL9baaaaa@75A7@ :

θ j ( 0 ) = θ j MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqiUde3aa0baaSqaaiabdQgaQbqaaiabcIcaOiabicdaWiabcMcaPaaakiabg2da9iabeI7aXnaaDaaaleaacqWGQbGAaeaacqGHxiIkaaaaaa@36F4@


θ j ( 0 ) = θ j ( i ) MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqiUde3aa0baaSqaaiabdQgaQbqaaiabcIcaOiabicdaWiabcMcaPaaakiabg2da9iabeI7aXnaaDaaaleaacqWGQbGAaeaacqGGOaakcqWGPbqAcqGGPaqkaaaaaa@3912@


  1. Tomlin C, Axelrod J: Biology by numbers: mathematical modelling in developmental biology. Nature Reviews Genetics. 2007, 8: 331-340. 10.1038/nrg2098

    Article  CAS  PubMed  Google Scholar 

  2. Voit E: Computational Analysis of Biochemical Systems. 2000, Cambridge University Press

    Google Scholar 

  3. de Jong H: Modeling and simulation of genetic regulatory systems: a literature review. J Comput Biol. 2002, 9: 67-103. 10.1089/10665270252833208

    Article  CAS  PubMed  Google Scholar 

  4. Rogers S, Khanin R, Girolami M: Bayesian model-based inference of transcription factor activity. BMC Bioinformatics. 2006, 8 Suppl 2: S2-

    Google Scholar 

  5. Alon U: An Introduction to Systems Biology: Design Principles of Biological Circuits. 2007, Chapman and Hall/CRC

    Google Scholar 

  6. Sivia D: Data Analysis: A Bayesian Tutorial. 2006, Oxford University Press, 2

    Google Scholar 

  7. 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: 3067-74. 10.1093/bioinformatics/btl485

    Article  CAS  PubMed  Google Scholar 

  8. Andrieu C, de Freitas N, Doucet A, Jordan MI: An introduction to MCMC for machine learning. Machine Learning. 2003, 50: 5-43. 10.1023/A:1020281327116.

    Article  Google Scholar 

  9. Strogatz SH: Nonlinear Dynamics and Chaos. 1994, Addison Wesley

    Google Scholar 

  10. Verdugo A, Rand R: Hopf bifurcation in a DDE model of gene expression. Communications in Nonlinear Science and Numerical Simulation. 2008, 13: 235-242. 10.1016/j.cnsns.2006.05.001.

    Article  Google Scholar 

  11. Bar-Or R, Maya R, Segel L, Alon U, Levine A, Oren M: Generation of oscillations by the p53-Mdm2 feedback loop: a theoretical and experimental study. Proceedings of the National Academy of Sciences USA. 2000, 97: 11250-11255. 10.1073/pnas.210171597.

    Article  CAS  Google Scholar 

  12. Monk NA: Oscillatory Expression of Hes1, p53, and NF-kB Driven by Transcriptional Time Delays. Current Biology. 2003, 13: 1409-1413. 10.1016/S0960-9822(03)00494-9

    Article  CAS  PubMed  Google Scholar 

  13. 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 

  14. Bernard S, Čajavec B, Pujo-Menjouet 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): 1155-70. 10.1098/rsta.2006.1761

    Article  CAS  Google Scholar 

  15. Barrio M, Burrage K, Leier A, Tian T: Oscillatory Regulation of Hes1: Discrete Stochastic Delay Modelling and Simulation. PLoS Computational biology. 2006, 2: 1017-1030. 10.1371/journal.pcbi.0020117.

    Article  CAS  Google Scholar 

  16. Wilkinson DJ: Stochastic Modelling for Systems Biology. 2006, Chapman and Hall/CRC

    Google Scholar 

  17. Beaumont MA, Rannala B: The Bayesian revolution in genetics. Nature Reviews Genetics. 2004, 5: 251-261. 10.1038/nrg1318

    Article  CAS  PubMed  Google Scholar 

  18. Gelman A, Carlin JB, Stern HS, Rubin DB: . Bayesian Data Analysis. 2004, Chapman and Hall CRC

    Google Scholar 

Download references


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.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Catherine F Higham.

Authors’ original submitted files for images

Rights and permissions

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 (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Higham, C.F. Bifurcation analysis informs Bayesian inference in the Hes1 feedback loop. BMC Syst Biol 3, 12 (2009).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: