Parameter values
All concentration parameters have units (number of molecules per cell), and all time units are (hours) or (per hour) for rates. A summary of default parameter values for the core circadian oscillator, based on Ref. [11], is provided (Table 1). The default parameters are generally consistent with full production response with a few activator proteins per cell, activated mRNA production rates of about one per five minutes, mRNA lifetimes of about 10 minutes, and translation rates of 5 proteins per transcript per hour.
The parameter K
A
, with units of molecules per cell, provides a concentration scale. Concentration and production-rate parameters from Ref. [11] were reduced by a factor of 50 to give K
A
= 1, corresponding to half activation at one transcription factor protein per cell. The product of the bimolecular rate constant k and a concentration yields a rate. Since concentrations were reduced by a factor of 50, the bimolecular rate constant was increased by a factor of 50 to yield an identical effective rate, maintaining exactly the same balance between bimolecular collision rates and protein decay rates.
Other than the concentration scaling, two other changes were made relative to Ref. [11]. First, the baseline repressor production rate was set to 0 mol/hr, compared to 0.1 mol/hr previously, or 0.002 mol/hr after the 50× concentration scaling. Second, the Hill exponent was set to 2 rather than 1, reflecting that the transcription factors in this system may bind cooperatively as dimers rather than monomers [9].
Several other parameter sets for more complete models of the circadian clock in Drosophila and mammalian systems are available [21–25]. These parameter sets are quite different, and rather than attempting to reconcile them we scanned parameters over a range of values.
A first-principles route to estimating the value of the bimolecular rate constant is to calculate a diffusion-limited reaction rate for proteins. The bimolecular rate constant for two proteins with diffusion constants D1 and D2 and effective radius h is
k = 4π(D1 + D2)h/V
where V corresponds to a spherical volume in which diffusion occurs [26], which we take as the nuclear diameter. Typical diffusion constants for proteins are in the range of 10-6 to 10-7 cm2/sec. Using a nuclear diameter of roughly 7 μ m (volume = 1.8 × 10-16 m3) and an effective protein radius of 1 nm gives a bimolecular rate constant k ~150 to 165/(molecule-hr) when A and R have units (molecules)/(nuclear volume).
We followed the Barkai-Leibler model in assuming that the activator molecule can be degraded while free and in a complex, while the repressor can only be degraded when free and is protected in the complex. We also assumed, as in the original model, that the activator decay rate is the same for free and bound molecules.
Justification of steady-state approximations for mRNA
Consider a protein X with mRNA X
m
with dynamics
(22)
where β'(t) is an mRNA synthesis rate that has explicit time dependence, for example due to fluctuating concentrations of a transcription factor. The protein synthesis rate, β"X
m
, has time dependence through the mRNA concentration X
m
with a constant coefficient β". Using the notation that (s) is the Laplace transform of f(t),
(23)
The transfer function relating protein output to mRNA production input is β"(s + α')-1(s + α")-1. When mRNA decay rates are fast compared to protein decay rates, the ratio α"/α' can be used as a small expansion parameter, yielding
(24)
The notation O(...) indicates asymptotic order. Transforming back to the time domain yields
(25)
with effective initial conditions
Xeff(0) = X(0) + (β"/α')X
m
(0) + O(α"/α').
In the model presented here, protein and mRNA concentrations are both close to 0 at time 0, and the offset to X(0) due to mRNA may be ignored.
ODE calculations
Solutions to ODEs were obtained with R [27] using the interface to lsoda [28, 29].
Stochastic simulations
Stochastic simulations were performed using the Gillespie stochastic simulation algorithm [30, 31] as implemented in R by the GillespieSSA package [32]. Simulations with a maximum particle count below 1000 used the direct method. Simulations with a maximum count of 1000 or more were accelerated using the binomial tau-leap method [33, 34]. Six equations define the stochastic dynamics:
(27)
Separating fast and slow subsystems can improve the performance of stochastic simulations [35].
Analytical approximations
The analytical approximations are motivated by the timescale separation between fast diffusion-limited bimolecular collisions between A and R, and the remaining slow timescales in the system. When the concentrations A and R are each 1 molecule/cell, for example, and k is the diffusion-limited rate constant 100/(molecule × hr), the formation rate of C is = 100 molecules/hr. In comparison, the maximum production rates of A and R are expected to be only 10 to 20 molecules/hr for typical parameters, and decay rates are slower still.
With this timescale separation, the bimolecular reaction effectively goes to completion: all available molecules of A and R are paired into complexes, until either A ≈ 0 or R ≈ 0. Depending on whether A or R is limiting, this leads to the approximation that ≈ 0 or ≈ 0. These approximations lead to analytical expressions, as detailed below.
Recovery from time 0 to full activation at time t
K
The start of the cycle, time 0, is defined when A and R are small, with A just crossing above R, and = kAR - α
A
C ≈ 0. In this regime, (d/dt)(C + A) = β
a
- α
A
(A + C) + (β
A
- β
a
)An/[An+ ], which is approximately equal to β
a
- α
A
(A + C) for A <<K
A
. At the end of the previous cycle, the sum A + C therefore approached an equilibrium value of β
a
/α
A
with transient time on the order of . As A is small, the concentration C at time 0 is approximately β
a
/α
A
. Using ≈ 0 and A = R,
(28)
defines the crossing value X.
As A becomes activated but is still below K
A
and K
R
, the Hill functions are approximated to yield the equations
(29)
During this phase, A is much smaller than K
A
and K
R
, is small, and α
R
R and β
r
are small relative to activated transcription of R, yielding the approximate dynamics
(30)
Recovery of A implies > 0 and β
X
> 0. The decay term α
A
A is assumed to be small compared to the production rate, which will be seen in the results to be a good approximation.
The dynamics of Eq. 30 continue until the turning point of the Hill function, which occurs at A ≈ K
A
.
Terming t
K
as the time when A = K
A
in this regime, Eq. 30 is integrated implicitly to yield
(31)
In the case that n = 1, it is possible to retain the term α
A
A in the analytical solution, and the corresponding equations are
(32)
With typical parameters, β
X
≈ β
A
/K
A
> α
A
, and α
A
can be ignored as with the n > 1 case. This approximation ignores the very short transient in which A is between K
A
and K
R
. As will be seen from the results, very little error is made.
The activator makes a round trip
The next analytical regime begins with A = K
A
and rising rapidly, C ≈ β
a
/α
A
, and R ≈ 0. Here full transcriptional activation is a good approximation, A >> K
A
and A >> K
R
, and the Hill functions become 1. The dynamics for R are = β
R
- α
R
R - kAR + α
A
C. Since and R are both small, kAR ≈ β
R
+ α
A
C. The dynamics for C are therefore = β
R
, or
C(t) = (β
a
/α
A
) + β
R
(t - t
K
).
An interpretation of this equation is that every new molecule of R, produced at rate β
R
, is converted immediately to a complex, adding to the baseline amount of complex at time t
K
.
The trajectory for A using Eq. 33 is
(34)
The time when A reaches its maximum and the maximum value are
(35)
These equations continue to hold until A returns to value K
A
. The time t'
K
when this occurs is obtained from Eq. 34 using the approximation that the exponential transient has decayed and that β
r
is negligible,
(36)
from Eq. 33. The complex will be seen to have obtained its maximum value Cmax at this time.
Once A falls below the threshold K
A
, it continues to fall rapidly and the Hill functions are approximately 0. The conditions that ≈ 0 and R ≈ 0 yield kAR = α
A
C, giving ≈ 0 as well. The dynamics for A in this regime are
(37)
These dynamics continue until A = 0, at which point C is still equal to Cmax. Defining this time as t
C
,
(38)
The repressor makes a round trip
In this phase of the cycle, A becomes the slow subsystem because A and hence are small, while R and C change more rapidly. The approximations = 0 with α
A
A = 0 yield kAR = β
a
, and
(39)
The trajectories are
(40)
(41)
in the case that α
A
= α
R
= α.
Ignoring the small contribution from the baseline production rate β
r
, the time t
R
for the maximum of R and its maximum value are
(42)
in the case that α
A
= α
R
= α.
An improved approximation for A in this region, again from = 0, is A(t) ≈ β
a
/[α
A
+ kR(t)]. According to this expression, A has its minimum value when R is at its maximum. These dynamics continue until R = X = , at which point R ≈ A and the cycle restarts. For convenience, let αmin = min(α
A
, α
R
) and Δα = |α
A
- α
R
|. Again ignoring the small contribution from β
r
, the crossing time t
X
is calculated as
(43)
The final term in t
X
is a perturbative estimate of the correction due to the faster transient. When α
A
= α
R
= α, an iterative solution is used for τ
X
= t
X
- t
C
. Denoting Δβ = β
A
- β
a
- α
A
K
A
and , and assuming that β
r
/α <<β
A
τ
X
,
(44)
where in the last equation the term τ
X
on the right-hand side has been replaced by its leading order value α-1 ln [Δβ/α ΔX] to terminate the expansion.
Time intervals
For convenience, time intervals are denoted with symbol τ as
(45)
Simplified formulas for the intervals are obtained by using leading-order contributions based on biologically reasonable assumptions that threshold values (K
A
and K
R
) are small compared to the activated level of A, and that the fully activated production rate for A, β
A
is much larger than the baseline rate β
a
or the fully activated production rate β
R
of repressor, and that β
r
is 0.
In these limits,
(46)
Both τ
K
and τ
C
are fast, scaling as inverse production rates, while τ'
K
and τ
X
are slow, scaling as inverse decay rates. To compare with numeric simulation results, we break the entire cycle into two phases with times τ1 and τ2. The first phase begins at the start of the cycle with A = R and > 0 and ends when C is at its maximum. The second phase begins when C is at its maximum and ends at the start of the next cycle:
(47)
Energy cost
The energy cost of sustaining oscillations is estimated as the number of protein molecules synthesized per cycle, giving equal weight to activator and repressor proteins. During the first phase, activator proteins are synthesized at rate β
A
and repressor proteins at rate β
R
. During the second phase, synthesis rates are β
a
and β
r
. The cost per cycle is
Cost = (β
A
+ β
R
)τ1/τtot + (β
a
+ β
r
)τ2/τtot.
Repressilator dynamics
The repressilator is a synthetic oscillator constructed from an odd cycle of negative feedback loops [6].
Again using a steady-state approximation for mRNA levels, the dynamics for the standard three-component symmetric continuous repressilator are
where P
i
represents one of i = 1 to 3 components, and Pi-1= P3 when i = 1.
Using the same approximations as for the hysteresis oscillator, the repressilator cycle starts when P3 falls just below K, with P1 at its baseline level β0/α and P2 roughly equal to (β1/α) - K. In this regime,
These dynamics continue until P1(t) = K, defined to occur after an interval τ1,
τ1 = α-1 ln [(β1 - β0)/(β1 - α K)] ≈ K/β1,
where the approximation holds when β1/α >> K > β0/α. At this point,
(52)
with dynamics
(53)
These dynamics continue until time τ1 + τ2 when P2(τ1 + τ2) = K,
τ2 = α-1 ln [(β1 - αK - β0)/(αK - β0)] ≈ α-1 ln(β1/αK).
This pattern repeats three times for a total period τtot of
τtot = 3(τ1 + τ2) = (3/α) ln [(β1 - β0)(β1 - αK - β0)/(β1 - αK)(αK - β0)].
Energy cost
During the τ+ phase, components 1 and 2 are synthesized at rate β1, while component 3 is synthesized at rate β0. During the τ _ phase, only component 1 is synthesized at the higher rate. The energy consumption, in terms of proteins per cycle, is
(56)
Implementation and availability
The analytical and numerical equations are provided as R code as supplementary material released under the GNU Lesser General Public License version 3 [see Additional file 1]. Source code is also available from the author at http://www.baderzone.org.