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 \dot{m} = 0 and \dot{p} = 0. After elimination and substitution, we obtain two expressions for *p** and *m**:

\begin{array}{rrr}\hfill {({p}^{\ast})}^{n+1}+{p}_{0}^{n}{p}^{\ast}-\frac{{p}_{0}^{n}}{{\mu}_{m}{\mu}_{p}}& \hfill =& \hfill 0,\end{array}

(20)

{m}^{\ast}={\mu}_{p}{p}^{\ast}.

(21)

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:

\dot{\zeta}=-{\mu}_{m}\zeta -K{\eta}_{d},

(22)

\dot{\eta}=\zeta -{\mu}_{p}\eta ,

(23)

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

\ddot{\eta}+({\mu}_{m}+{\mu}_{p})\dot{\eta}+{\mu}_{m}{\mu}_{p}\eta +K{\eta}_{d}=0.

(24)

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

\lambda =\frac{-({\mu}_{m}+{\mu}_{p})\pm \sqrt{{({\mu}_{m}+{\mu}_{p})}^{2}-4({\mu}_{m}{\mu}_{p}+K)}}{2}.

(26)

Examination of equation (26) reveals that since *μ*_{
m
}+ *μ*_{
p
}> 0 and \sqrt{{({\mu}_{m}+{\mu}_{p})}^{2}-4({\mu}_{m}{\mu}_{p}+K)}<{\mu}_{m}+{\mu}_{p}, the real part of *λ* will always be negative, implying stability. For

4 (*μ*_{
m
}*μ*_{
p
}+ *K*) > (*μ*_{
m
}+ *μ*_{
p
})^{2}

*λ* has a 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 4*K* > (*μ*_{
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

\ddot{\eta}+({\mu}_{m}+{\mu}_{p})\dot{\eta}+{\mu}_{m}{\mu}_{p}\eta =-K{\eta}_{d}+{H}_{2}{\eta}_{d}^{2}+{H}_{3}{\eta}_{d}^{3}+\dots

(30)

where the coefficients *K*, *H*_{2} and *H*_{3} are obtained by Taylor series expansion of the nonlinear term.

Assuming that the true solution is periodic, our goal is to find an analytical approximation that is valid for all *t*. The key idea is to regard the frequency *ω* as unknown in advance, and to solve for it by demanding that *η* contains no secular terms. We introduce a small parameter ϵ, and let Δ = *ϵ*^{2}*δ*, with

*η* = *ϵu*,

*τ* = *τ*_{
cr
}+ Δ = *τ*_{
cr
}+ *ϵ*^{2}*δ*.

We then stretch time by replacing the independent variable *t* by *σ*, where

*σ* = Ω*t*.

We then have

{\Omega}^{2}\frac{{d}^{2}u}{d{\sigma}^{2}}+({\mu}_{m}+{\mu}_{p})\Omega \frac{du}{d\sigma}+{\mu}_{m}{\mu}_{p}u=-K{u}_{d}+\u03f5{H}_{2}{u}_{d}^{2}+{\u03f5}^{2}{H}_{3}{u}_{d}^{3}+\dots ,

(34)

where *u*_{
d
}= *u* (*σ* - Ω*τ*). Expanding Ω in a power series in ϵ, omitting the *O* (ϵ) term for convenience, since it turns out to be zero, we have

Ω = *ω* + *ϵ*^{2}*k*_{2} + ...

Next we expand the delay term *u*_{
d
}

*u*_{
d
}= *u* (*σ* - *ωτ*_{
cr
}) - *ϵ*^{2} (*k*_{2} *τ*_{
cr
}+ *ωδ*) *u*' (*σ* - *ωτ*_{
cr
}) + *O*(*ϵ*^{3})

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

*u* (*σ*) = *u*_{0} (*σ*) + *ϵu*_{1} (*σ*) + *ϵ*^{2}*u*_{2} (*σ*) + ...

Substituting and collecting terms we find

{\omega}^{2}\frac{{d}^{2}{u}_{0}}{d{\sigma}^{2}}+({\mu}_{m}+{\mu}_{p})\omega \frac{d{u}_{0}}{d\sigma}+{\mu}_{m}{\mu}_{p}{u}_{0}=-K{u}_{0}(\sigma -\omega {\tau}_{cr})

(38)

{\omega}^{2}\frac{{d}^{2}{u}_{1}}{d{\sigma}^{2}}+({\mu}_{m}+{\mu}_{p})\omega \frac{d{u}_{1}}{d\sigma}+{\mu}_{m}{\mu}_{p}{u}_{1}=-K{u}_{1}(\sigma -\omega {\tau}_{cr})+{H}_{2}{u}_{0}^{2}(\sigma -\omega {\tau}_{cr})

(39)

\begin{array}{lll}{\omega}^{2}\frac{{d}^{2}{u}_{2}}{d{\sigma}^{2}}+({\mu}_{m}+{\mu}_{p})\omega \frac{d{u}_{2}}{d\sigma}+{\mu}_{m}{\mu}_{p}{u}_{2}\hfill & =\hfill & -K{u}_{2}(\sigma -\omega {\tau}_{cr})+{H}_{3}{u}_{0}^{3}(\sigma -\omega {\tau}_{cr})\hfill \\ +2{H}_{2}{u}_{0}(\tau -\omega {\tau}_{cr}){u}_{1}(\tau -\omega {\tau}_{cr})\hfill \\ +K({k}_{2}{\tau}_{cr}+\omega \delta ){{u}^{\prime}}_{0}(\tau -\omega {\tau}_{cr}).\hfill \end{array}

(40)

We take the solution of the *u*_{0} and *u*_{1} equations to have the form

{u}_{0}(\sigma )=\widehat{A}\mathrm{cos}\sigma ,

(41)

*u*_{1}(*σ*) = *m*_{1} sin 2*σ* + *m*_{2} sin 2*σ* + *m*_{3},

where *A* = \widehat{A}\u03f5 from (29) and (31). Next we substitute (41) and (42) into (40) to give

\begin{array}{cc}\begin{array}{r}\hfill {\omega}^{2}(-4{m}_{1}\mathrm{sin}2\sigma -4{m}_{2}\mathrm{cos}2\sigma )+\\ \hfill \omega ({\mu}_{m}+{\mu}_{p})(2{m}_{1}\mathrm{cos}2\sigma -2{m}_{2}\mathrm{sin}2\sigma )+\\ \hfill K{m}_{1}(\mathrm{sin}2\sigma \mathrm{cos}2\omega {\tau}_{cr}-\mathrm{cos}2\sigma \mathrm{sin}2\omega {\tau}_{cr})+\\ \hfill K{m}_{2}(\mathrm{cos}2\sigma \mathrm{cos}2\omega {\tau}_{cr}-\mathrm{sin}2\sigma \mathrm{sin}2\omega {\tau}_{cr})+\\ \hfill K{m}_{3}+{\mu}_{m}{\mu}_{p}({m}_{1}\mathrm{sin}2\sigma +{m}_{2}\mathrm{cos}2\sigma )\end{array}& \begin{array}{ll}=\hfill & \frac{1}{2}{H}_{2}{\widehat{A}}^{2}\mathrm{cos}2\sigma \mathrm{cos}2\omega {\tau}_{cr}+\hfill \\ \frac{1}{2}{H}_{2}{\widehat{A}}^{2}\mathrm{sin}2\sigma \mathrm{sin}2\omega {\tau}_{cr}+\frac{1}{2}{H}_{2}{\widehat{A}}^{2}.\hfill \end{array}\end{array}

Collecting sin 2*σ* terms gives

\begin{array}{cc}\begin{array}{r}\hfill -4{\omega}^{2}{m}_{1}-2({\mu}_{m}+{\mu}_{p})\omega {m}_{2}+\\ \hfill K{m}_{1}\mathrm{cos}2\omega {\tau}_{cr}+K{m}_{2}\mathrm{sin}2\omega {\tau}_{cr}+{\mu}_{m}{\mu}_{p}{m}_{1}\end{array}& \begin{array}{ll}=\hfill & \frac{1}{2}{H}_{2}{\widehat{A}}^{2}\mathrm{sin}2\omega {\tau}_{cr}.\hfill \end{array}\end{array}

Collecting cos 2*σ* terms gives

\begin{array}{cc}\begin{array}{r}\hfill -4{\omega}^{2}{m}_{2}-2({\mu}_{m}+{\mu}_{p})\omega {m}_{1}+\\ \hfill K{m}_{1}-\mathrm{sin}2\omega {\tau}_{cr}+K{m}_{2}\mathrm{cos}2\omega {\tau}_{cr}+{\mu}_{m}{\mu}_{p}{m}_{2}\end{array}& \begin{array}{ll}=\hfill & \frac{1}{2}{H}_{2}{\widehat{A}}^{2}\mathrm{cos}2\omega {\tau}_{cr}.\hfill \end{array}\end{array}

Collecting other terms gives

\begin{array}{l}K{m}_{3}+{\mu}_{m}{\mu}_{p}{m}_{3}\hfill \\ =\hfill & \frac{1}{2}{H}_{2}{\widehat{A}}^{2}.\hfill \end{array}

Letting

*α* = -4*ω*^{2} + *K* cos2*ωτ*_{
cr
}+ *μ*_{
m
}*μ*_{
p
},

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

{\widehat{H}}_{1}=\frac{1}{2}{H}_{2}{\widehat{A}}^{2}\mathrm{sin}2\omega {\tau}_{cr},

(45)

{\widehat{H}}_{2}=\frac{1}{2}{H}_{2}{\widehat{A}}^{2}\mathrm{cos}2\omega {\tau}_{cr},

(46)

{\widehat{H}}_{3}=\frac{1}{2}{H}_{2}{\widehat{A}}^{2},

(47)

we find that

{m}_{1}=\frac{\alpha {\widehat{H}}_{1}-\beta {\widehat{H}}_{2}}{{\alpha}^{2}+{\beta}^{2}},

(48)

{m}_{2}=\frac{\alpha {\widehat{H}}_{1}+\beta {\widehat{H}}_{2}}{{\alpha}^{2}-{\beta}^{2}},

(49)

{m}_{3}=\frac{{\widehat{H}}_{3}}{K+{\mu}_{m}{\mu}_{p}}.

(50)

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

\begin{array}{cc}\begin{array}{r}\hfill {k}_{2}(-{\mu}_{m}-{\mu}_{p}+K{\tau}_{cr}\mathrm{cos}\omega {\tau}_{cr})\\ \hfill +\omega \delta K\mathrm{cos}\omega {\tau}_{cr}\end{array}& \begin{array}{ll}=\hfill & \frac{3}{4}{H}_{3}{\widehat{A}}^{2}\mathrm{sin}\omega {\tau}_{cr}+\hfill \\ {H}_{2}^{2}{\widehat{A}}^{2}({m}_{1}\mathrm{cos}\omega {\tau}_{cr}+{m}_{2}\mathrm{sin}\omega {\tau}_{cr}+2{m}_{3}\mathrm{sin}\omega {\tau}_{cr})\hfill \end{array}\end{array}

and

\begin{array}{cc}\begin{array}{r}\hfill {k}_{2}(-2\omega -K{\tau}_{cr}\mathrm{sin}\omega {\tau}_{cr})\\ \hfill -\omega \delta K\mathrm{sin}\omega {\tau}_{cr}\end{array}& \begin{array}{ll}=\hfill & \frac{3}{4}{H}_{3}{\widehat{A}}^{2}\mathrm{cos}\omega {\tau}_{cr}+\hfill \\ {H}_{2}^{2}{\widehat{A}}^{2}(-{m}_{1}\mathrm{sin}\omega {\tau}_{cr}+{m}_{2}\mathrm{cos}\omega {\tau}_{cr}+2{m}_{3}\mathrm{cos}\omega {\tau}_{cr}).\hfill \end{array}\end{array}

Let

*c*_{1} = -*μ*_{
m
}- *μ*_{
p
}+ *Kτ*_{
cr
}cos*ωτ*_{
cr
},

*c*_{2} = -2*ω* - *Kτ*_{
cr
}sin *ωτ*_{
cr
},

{d}_{1}=\frac{3}{4}{H}_{3}\mathrm{sin}\omega {\tau}_{cr}+{H}_{2}^{2}({m}_{1}\mathrm{cos}\omega {\tau}_{cr}+{m}_{2}\mathrm{sin}\omega {\tau}_{cr}+2{m}_{3}\mathrm{sin}\omega {\tau}_{cr}),

(53)

{d}_{2}=\frac{3}{4}{H}_{3}\mathrm{cos}\omega {\tau}_{cr}+{H}_{2}^{2}(-{m}_{1}\mathrm{sin}\omega {\tau}_{cr}+{m}_{2}\mathrm{cos}\omega {\tau}_{cr}+2{m}_{3}\mathrm{cos}\omega {\tau}_{cr}).

(54)

Then we have

{A}^{2}=\frac{({c}_{2}\mathrm{cos}\omega {\tau}_{cr}+{c}_{1}\mathrm{sin}\omega {\tau}_{cr})}{{d}_{1}{c}_{2}-{c}_{1}{d}_{2}}\omega \Delta K,

(55)

{k}_{2}=-\frac{{d}_{2}\mathrm{cos}\omega {\tau}_{cr}+{d}_{1}\mathrm{sin}\omega {\tau}_{cr}}{{c}_{1}{d}_{2}-{c}_{2}{d}_{1}}\omega \delta K,

(56)

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

Ω = *ω* + *k*_{2}Δ,

which has the form of equation (15) (noting that *k*_{2} is negative).

### Maximum *ω*

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
}= \sqrt{c}. Rearranging equation (8) gives

{\omega}^{4}+({\mu}_{m}^{2}+{\mu}_{p}^{2}){\omega}^{2}+{\mu}_{m}^{2}{\mu}_{p}^{2}-{K}^{2}=0.

(58)

Substituting *μ*_{
m
}*μ*_{
p
}= *c* and *K* = *c*_{
K
}into the above equation gives:

{\omega}^{4}+{\omega}^{2}\left(\frac{c}{{\mu}_{p}^{2}}+{\mu}_{p}^{2}\right){\omega}^{2}+{c}^{2}-{c}_{K}^{2}=0.

(59)

Differentiating w.r.t. *μ*_{
p
}gives:

4{\omega}^{3}\frac{d\omega}{d{\mu}_{p}}+2\omega \frac{d\omega}{d{\mu}_{p}}\left(\frac{{c}^{2}}{{\mu}_{p}^{2}}+{\mu}_{p}^{2}\right)+{\omega}^{2}\left(-\frac{2{c}^{2}}{{\mu}_{p}^{3}}+2{\mu}_{p}\right)=0,

(60)

So \frac{d\omega}{d{\mu}_{p}} = 0 when

\begin{array}{r}\hfill {\omega}^{2}\left(\frac{-2{c}^{2}}{{\mu}_{p}^{3}}+2{\mu}_{p}\right)=0\\ \hfill \frac{-2{c}^{2}}{{\mu}_{p}^{3}}+2{\mu}_{p}=0\end{array}

(61)

giving

As *μ*_{
m
}*μ*_{
p
}= *c*, we have *μ*_{
m
}= *μ*_{
p
}= \sqrt{c}. The maximum value of *ω* therefore occurs when *μ*_{
m
}= *μ*_{
p
}and has the form

\omega =\sqrt{K-{\mu}^{2}}.

(63)

Writing *K* - *μ*_{2} in terms of *β* gives

K-{\mu}^{2}=\frac{n\beta -(1+\beta )}{{p}_{0}{\beta}^{1/n}{(1+\beta )}^{2}}.

(64)

Differentiating with respect to *β* gives

\frac{{\beta}^{1/n}{(1+\beta )}^{2}(n-1)-(n\beta -\beta -1)\left(\frac{1}{n}{\beta}^{1/n-1}{(1+\beta )}^{2}+2{\beta}^{1/n}(1+\beta )\right)}{{p}_{0}{\left({\beta}^{1/n}{(1+\beta )}^{2}\right)}^{2}}.

(65)

Setting this expression to zero and dividing by *β*^{1/n}(1 + *β*)/*p*_{0} gives

n+n\beta -\beta -1-(n-1)(1/n(1+\beta )+2\beta )+\frac{1}{n\beta}(1+\beta )+2=0,

(66)

which simplifies to

(*n*^{2} - 1) *β*^{2} - (*n*^{2} + 2) *β* - 1 = 0.

### Can the data tell us anything about the parameters directly?

Integrating equation (2) gives

{\displaystyle {\int}_{0}^{T}\dot{p}(t)dt}={\displaystyle {\int}_{0}^{T}m(t)dt}-{\mu}_{p}{\displaystyle {\int}_{0}^{T}p(t)dt}.

(68)

From the Fundamental Theorem of Calculus, it follows that

p(T)-p(0)={\displaystyle {\int}_{0}^{T}m(t)dt}-{\mu}_{p}{\displaystyle {\int}_{0}^{T}p(t)dt}.

(69)

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

{\mu}_{p}=\frac{{m}_{ave}}{{p}_{ave}}.

(70)

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, {\theta}_{j}^{\ast}, (*j* = 1 ... 5 in Experiment 1 and *j* = 1 ... 6 in Experiment 2) in logspace, based on the current value of {\theta}_{j}^{(i)}, using a Gaussian *proposal distribution*, q\left({\theta}_{j}^{\ast}|{\theta}_{j}^{(i)}\right) given by

q\left({\theta}_{j}^{\ast}|{\theta}_{j}^{(i)}\right)=\frac{1}{{\sigma}_{j}\sqrt{2\pi}}\mathrm{exp}\left(-\frac{{\left({\theta}_{j}^{\ast}-{\theta}_{j}^{(i)}\right)}^{2}}{2{\sigma}_{j}^{2}}\right).

(71)

*σ*_{
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*, \mathcal{A}(*θ*^{(i)}, *θ* *), given by

\mathcal{A}\left({\theta}_{j}^{(i)},{\theta}_{j}^{\ast}\right)=\mathrm{min}\left\{1,\frac{\widehat{p}\left({\theta}_{j}^{\ast}\right)q\left({\theta}_{j}^{(i)}|{\theta}_{j}^{\ast}\right)}{\widehat{p}\left({\theta}_{j}^{(i)}\right)q\left({\theta}_{j}^{\ast}|{\theta}_{j}^{(i)}\right)}\right\},

(72)

where

\widehat{p}=p(\mathcal{D}|{\theta}_{j})p({\theta}_{j}),

(73)

p(\mathcal{D}|{\theta}_{j})=p(\left\{{x}_{k}\right\}|{\theta}_{j},\sigma ,I),

(74)

and

p({\theta}_{j})=\mathcal{N}({\mu}_{j},{\sigma}_{pj}).

(75)

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

**step 1** Initialise {\theta}_{j}^{(0)}

**step 2** For *i* = 0 to *N*:

Sample *r* ~ \mathcal{N}(0, 1)

Sample {\theta}_{j}^{\ast}~q\left({\theta}_{j}^{\ast}|{\theta}_{j}^{(i)}\right) (choosing *j* at random

If r<\mathcal{A}\left({\theta}_{j}^{(i)},{\theta}_{j}^{\ast}\right)=\mathrm{min}\left\{1,\frac{\widehat{p}\left({\theta}_{j}^{\ast}\right)q\left({\theta}_{j}^{(i)}|{\theta}_{j}^{\ast}\right)}{\widehat{p}\left({\theta}_{j}^{(i)}\right)q\left({\theta}_{j}^{\ast}|{\theta}_{j}^{(i)}\right)}\right\}:

{\theta}_{j}^{(0)}={\theta}_{j}^{\ast}

else

{\theta}_{j}^{(0)}={\theta}_{j}^{(i)}