Analytical approximations for the amplitude and period of a relaxation oscillator

Carmen Kut; Vahid Golkhou; Joel S Bader

doi:10.1186/1752-0509-3-6

Research article
Open Access

Analytical approximations for the amplitude and period of a relaxation oscillator

Carmen Kut¹,
Vahid Golkhou^{1, 2} and
Joel S Bader^{1, 2, 3}Email author

BMC Systems Biology20093:6

https://doi.org/10.1186/1752-0509-3-6

Received: 05 September 2008
Accepted: 14 January 2009
Published: 14 January 2009

Abstract

Background

Analysis and design of complex systems benefit from mathematically tractable models, which are often derived by approximating a nonlinear system with an effective equivalent linear system. Biological oscillators with coupled positive and negative feedback loops, termed hysteresis or relaxation oscillators, are an important class of nonlinear systems and have been the subject of comprehensive computational studies. Analytical approximations have identified criteria for sustained oscillations, but have not linked the observed period and phase to compact formulas involving underlying molecular parameters.

Results

We present, to our knowledge, the first analytical expressions for the period and amplitude of a classic model for the animal circadian clock oscillator. These compact expressions are in good agreement with numerical solutions of corresponding continuous ODEs and for stochastic simulations executed at literature parameter values. The formulas are shown to be useful by permitting quick comparisons relative to a negative-feedback represillator oscillator for noise (10× less sensitive to protein decay rates), efficiency (2× more efficient), and dynamic range (30 to 60 decibel increase). The dynamic range is enhanced at its lower end by a new concentration scale defined by the crossing point of the activator and repressor, rather than from a steady-state expression level.

Conclusion

Analytical expressions for oscillator dynamics provide a physical understanding for the observations from numerical simulations and suggest additional properties not readily apparent or as yet unexplored. The methods described here may be applied to other nonlinear oscillator designs and biological circuits.

Keywords

Circadian Clock
Stochastic Simulation
Fast Subsystem
Hill Function
Bimolecular Rate Constant

Background

Analytical expressions for dynamical systems are useful for mapping underlying parameters to observed properties. For many mechanical, electrical, and atomic systems, analysis proceeds by reducing a complicated system to tractable linear system, more often than not involving coupled harmonic oscillators. The effective analytical dynamics then provide valuable intuition even when exact results are eventually obtained using numerical methods.

For dynamics of many biological networks [1], a basic tractable component is a simple switch, with effective dynamics

\dot{X} (t) = β (t) - α X (t) .

(1)

The symbol X represents the concentration of a molecule, for example a transcript, its encoded protein, or a specific protein modification state; β(t) is a time-dependent production rate, and α is a time-independent decay rate corresponding to a lifetime of α^-1. The dynamics of X can be obtained by convolution of the input β(t) with the response function e^-αt.

For gene regulatory networks, the input β(t) can often be modeled as an on-off toggle, β(t) = 0 in the repressed state and β(t) = β for full activation. This behavior arises naturally from multimeric binding of transcription factors, giving a sigmoidal Hill function as a function of transcription factor concentration.

The equilibrium behavior of this model is X = 0 for the repressed state and X = β/α for the activated state. Transients are also readily calculated. If X is itself a transcription factor, the relevant transient time is the delay between a change in the regulation of X and the subsequent change in regulation of its target genes. This introduces a second concentration, K, representing the concentration of X for half response of its target genes. For a switch from β(t) = 0 for t < 0 to β(t) = β for t ≥ 0, the transient time is α^-1 ln [1 - K/(β/α)], which approaches K/(β/α) when the threshold K is a fraction of the full response β/α. For a switch from β(t) = β for t < 0 to β = 0 for t ≥ 0, the transient time for decay from β/α to K is α^-1 ln [(β/α)/K], dominated by the protein lifetime α^-1 and more weakly dependent on the ratio of maximum to threshold concentration.

These basic equations can be used to predict the dynamics of bistable systems, such as metabolic switches [2]. Serial chains can give developmental progressions, such as the bacterial flagella gene network [3]. Negative feedback between simple switches can lead to bistable response, as observed in Delta-Notch signaling [4] and used to create a synthetic toggle switch [5].

Sustained oscillations can by created by coupling simple switches in a cycle, with each switch negatively regulating its successor, termed a repressilator [6]. The repressilator's amplitude and period can be estimated with good accuracy using the dynamics of each switch, giving an amplitude of β/α and a period of roughly 3α^-1 ln [(β/α)/K] for the three-component repressilator. Although numerical simulations are essential for a full quantitative understanding, the analytical results clearly provide intuition regarding the parameters and parameter ratios that define the oscillator behavior. For example, a stronger promoter will increase β, increasing the amplitude proportionally and increasing the period with much weaker logarithmic dependence. This insight can aid in understanding the differences between observed gene and protein circuits, and knowing which knobs to tweak when designing synthetic circuits. This direct connection also enables design of oscillators with desired period and amplitude, an important prerequisite for standardizing synthetic biology [7].

Although built in the laboratory, repressilators do not seem to be common in nature. Instead, hysteresis oscillators are thought to provide biological clocks for processes as diverse as neural signaling, basic metabolism, and development [8]. The best studied examples may be circadian clocks responsible for synchronizing living systems with day/night cycles, which are thought to have evolved independently in prokaryotes, cyanobacteria, fungi, plants, and animals [9]. A reduced model for the animal clock was introduced by Barkai and Leibler [10] (Fig. 1). Unfortunately, until now, no simple scaling rules have yet been provided oscillations arising from relaxation or hysteresis from coupled positive and negative feedback loops.

Figure 1
**Schematic oscillator designs. (Left) The repressilator is a delay oscillator here idealized as three symmetric repressive components,** P₁, P₂**, and** P₃**, with identical production rates,** β **, and decay rates,** α. (Right) The hysteresis oscillator has interlocked feedback loops involving an activator, A, and a repressor, R, which form a complex, C, with a bimolecular rate constant k. The activator and repressor have baseline production rates β_aand β_a, with total production rising to β_Aand β_Rwhen the concentration of free activator is high. The activator degrades at rate α_Awhether free or complexed; the repressor degrades at rate α_Rbut only when free.

In this oscillator, an activator protein, A, activates transcription of both itself and a repressor, R. The repressor achieves its effect by forming a complex, C, with A, that does not activate transcription. Activation of A can generate a reservoir of C that serves as a continuing source of R in the absence of A, which reduces the level of A back to baseline. Once the R molecules have degraded, A can reactivate transcription and initiate a new cycle. These dynamics exhibit hysteresis when projected onto the A-R plane.

Numerical simulations indicate that hysteresis oscillators have less noisy periods than delay oscillators; in fact, noise can actually serve to prevent the clock from falling into a stable attractor, making it more robust [11]. Intracellular communication can also improve robustness [12]. Numerical models have been used to compare clocks using transcriptional versus post-translational repression elements [13]. Clocks based on the the positive-negative feedback design have been observed computationally to be more easily tuned to desired frequencies [14].

The key parameters of the hysteresis clock model are α_Aand α_R, the decay rates of the activator and repressor proteins; β_Aand β_R, the fully activated production rates; β_a, the baseline production rate of the activator, and k is the bimolecular rate constant for $A + R \to_{}^{k} C$ . The following scaling rules and approximations are developed:

1. The maximum concentration of each component is approximately proportional to β_A/α_A, the ratio of the production rate to the decay rate of the activator.

2. A second important concentration scale is the reset point $\sqrt{β_{a} / k}$ when the activator and repressor concentrations cross, rather than a nominal steady-state baseline concentration of β_a/α.

3. The period is roughly divided into two phases, activation and recovery. The activation phase has duration (β_A/α_A)/β_R, equivalent to the time required for repressor to titrate the equilibrium concentration of activator.

4. The recovery phase has duration $α_{R}^{- 1} \ln [(β_{A} / α_{A}) / \sqrt{β_{a} / k}]$ , equivalent for the time for repressor to decay from its maximum concentration to the reset point.

This preview of the full results is accurate in the limit that production and decay rates are faster for the activator than the repressor. More accurate (but slightly more complicated) expressions are derived in the Methods. The strategy is to separate the oscillator into fast and slow subsystems that depend on the oscillator phase: during activation, A and C are fast and R is slow, and during recovery, C and R are fast and A is slow. This strategy is distinct from treatments that identify the activator as the fast subsystem throughout the entire cycle [11].

The Results show that the analytical results are accurate over a wide range of parameter space, compared with the numerical solutions to corresponding ODEs and also a stochastic simulation at the original literature values [10]. More detailed comparisons across parameters are done using ODEs alone, as the focus of this work is obtaining tractable analytics for a nonlinear system rather than investigating important known differences between ODEs and stochastic systems for small particle counts [15], or for systems where stochastic dynamics are essential for generating oscillations [16].

The value of the analytical expressions is then demonstrated through a comparison of operating characteristics: the noise, quantified as the variance of the period due to variance in production and decay rates; the cost or inverse efficiency, defined as the rate of protein production averaged over a cycle; and the dynamic range, quantified in decibels as the log-ratio of the concentration of the activating component at its maximum and minimum values. We conclude with a physical interpretation of the clock formulas and use these formulas to interpret results of computational studies.

Results and discussion

Hysteresis oscillator model

The protein concentrations [A], [R], and [C] of the activator, repressor, and complex are in units of molecules/cell and are denoted A, R, and C when the meaning is clear by context. The terms A and R refer only to free molecules and do not include those contained in complex C. The corresponding mRNA concentrations for activator and repressor are A_mand R_m. Continuous concentrations are assumed throughout. The mathematical model is

\begin{matrix} {\dot{A}}_{m} = {β^{'}}_{a} + ({β^{'}}_{A} - {β^{'}}_{a}) A^{n} / [A^{n} + K_{A}^{n}] - {α^{'}}_{A} A_{m} \\ \dot{A} = {β^{″}}_{A} A_{m} - α_{A} A - k A R \\ {\dot{R}}_{m} = {β^{'}}_{r} + ({β^{'}}_{R} - {β^{'}}_{r}) A^{n} / [A^{n} + K_{R}^{n}] - {α^{'}}_{R} R_{m} \\ \dot{R} = {β^{″}}_{R} R_{m} - α_{R} R - k A R + α_{A} C \\ \dot{C} = k A R - α_{A} C . \end{matrix}

(2)

The parameters ${β^{'}}_{a}$ and ${β^{'}}_{r}$ are baseline transcriptional rates for A_mand R_m; ${β^{'}}_{A}$ and ${β^{'}}_{R}$ are the fully activated transcriptional rates. Transcriptional activation is represented by Hill functions with half-response at A = K_Afor the activator and A = K_Rfor the repressor. The same Hill exponent n is used for both activator and repressor. This exponent is related to the number of activator proteins that form a transcriptional complex, and cooperative binding can result in Hill coefficients larger than 1. Although n = 1 was used in the original model [10], transcriptional activation in the metazoan clock is thought to be due to (hetero)dimers [9]. If binding is cooperative, n = 2 may be more appropriate.

Because mRNA decay rates ${β^{'}}_{A}$ and ${β^{'}}_{R}$ are fast compared to protein decay rates, min^-1 compared to hr^-1, mRNA transients are brief compared to protein response. Taking the limit of fast mRNA response is equivalent to employing steady-state approximations for mRNA levels, ${\dot{A}}_{m}$ ≈ 0 and ${\dot{R}}_{m}$ ≈ 0 (see Methods).

The steady-state approximation incorporates mRNA dynamics implicitly through effective protein synthesis rates,

\begin{array}{l} β_{a} & = & {β^{″}}_{A} {β^{'}}_{a} / {α^{'}}_{A} \\ β_{A} & = & {β^{″}}_{A} {β^{'}}_{A} / {α^{'}}_{A} \end{array}

(3)

\begin{array}{l} β_{r} & = & {β^{″}}_{R} {β^{'}}_{r} / {α^{'}}_{R} \\ β_{R} & = & {β^{″}}_{R} {β^{'}}_{R} / {α^{'}}_{R}, \end{array}

(4)

and effective equations

\begin{array}{l} \dot{A} & = & β_{a} + (β_{A} - β_{a}) A^{n} / [A^{n} + K_{A}^{n}] - α_{A} A - k A R \\ \dot{R} & = & β_{r} + (β_{R} - β_{r}) A^{n} / [A^{n} + K_{R}^{n}] - α_{R} R - k A R + α_{A} C \\ \dot{C} & = & k A R - α_{A} C . \end{array}

(5)

For notational clarity, subscripted values of t are used to denote times relative to the start of the cycle, and subscripted values of τ refer to time intervals.

Central parameter values, presented in Table 1, were based on Ref. [11]. Trajectories generated using exactly these values, except with the Hill coefficient changed from 1 to 2 to reflect cooperative binding of transcription factors as dimers, are similar whether from the analytical approximate solutions to Eq. 5 derived below, the continuous-time ODE solution, or the corresponding stochastic dynamics for discrete particle numbers (Fig. 2).

Table 1

Circadian clock parameter values

Parameter	Units	Vilar et al. (Ref. [11])	This Work Default Value^a	This Work Range
α _A	hr^-1	1	1	10^-1 .. 10^0.5
β _a	(mol/cell) hr^-1	250	5	1 .. 10
β _A	(mol/cell) hr^-1	2500	50	10 .. 10²
K _A	(mol/cell)	50	1	1 .. 10^0.3
α _R	hr^-1	0.2	0.2	10^-1 .. 10^0.5
β _r	(mol/cell) hr^-1	0.1	0.0	0.001 .. 0.01
β _R	(mol/cell) hr^-1	500	10	10^0.5 .. 10^1.5
K _R	(mol/cell)	100	2	10^0.2..10^0.5
n	Hill exponent	1	2	1 .. 10
k	(mol/cell)^-1 hr^-1	2	100	20 .. 200

^aDefault values are based on Ref. [11] but with concentrations rescaled by a factor of 50, corresponding to half activation of A at one molecule per cell. The two values in bold, β_rand n, differ from Ref. [11].

Figure 2
**Trajectories of activator, repressor, and complex concentrations are displayed for dynamics from analytical expressions (top panel), continuous ODEs (middle), and stochastic simulations (bottom)**. Parameter values are exactly as in Ref. [11], provided in Table 1, except with Hill coefficient = 2 to reflect cooperative dimer binding.

Unlike other studies that focus on the differences between deterministic ODEs and stochastic simulations, our aim is to develop analytical expressions that reproduce ODE behavior. For this purpose, we scaled the concentration parameters so that the concentration K_Aof A giving half-maximum self-activation occurs at 1 molecule per cell, compared with 50 molecules per cell from Ref. [11]. This rescaling gives a bimolecular rate constant k of 100 (mol/cell)^-1 hr^-1 for a pair of molecules, close to a first-principles estimate of the diffusion-limited bimolecular collision rate within the nuclear volume (see Methods). In addition to the Hill coefficient of 2 mentioned above, we also set the baseline production rate of repressor, β_r, to zero, as opposed to the value of 0.002 that would be obtained from concentration scaling. Results using β_r= 0 and β_r= 0.002 are virtually identical (see Results, Comparison with ODE solutions).

Analytical expressions for period and phase

Here we provide an overview of the method based on two assumptions:

1. Bimolecular collisions are fast compared to protein synthesis and decay rates.

2. Fast collisions between activator and repressor molecules means that the A + R → C reaction effectively goes to completion, with either A ≈ 0 or R ≈ 0 at all times.

The second assumption permits A, R, and C to be calculated from pairwise combinations that eliminate the non-linear bimolecular term. The Methods uses an expansion of the Hill functions to derive tractable dynamical equations, with summary analytical expressions provided in Table 2. The main results are sketched here using a simplified logic approximation, replacing the Hill functions with step functions [1]. These results are accurate when A passes quickly through its threshold value, which occurs for much of parameter space.

Table 2

Analytical results for oscillator period and phase

Times	Hysteresis^a	Delay^b
τ₁ (A > R)	n = 1: (β_A/α_Aβ_R) + K_A/β_A+(K_A/β_A) ln(K_A/ $\sqrt{β_{a} / k}$ )	α^-1 ln [(β₁ - β₀)/(β₁ - αK)]
	n ≠ 1: (β_A/α_Aβ_R) + K_A/β_A+ [β_A(n- 1)]^-1 $K_{A}^{n}$ (β_a/k)^1-n
τ₂ (A <R)	α_A= α_R: α^-1 ln [(β_A/α)/ $\sqrt{β_{a} / k}$ ]	α^-1 ln [(β₁ - αK - β₀)/(αK - β₀)]
	α_A≠ α_R: $α_{\min}^{- 1}$ ln [(β_A/\|α_A- α_R\|)/ $\sqrt{β_{a} / k}$ ]
τ _tot	τ₁ + τ₂	3(τ₁ + τ₂)
τ_tot, limiting	α^-1 [β_A/β_R+ ln(β_A/α $\sqrt{β_{a} / k}$ )]	3α^-1 ln [β₁/αK]
Concentrations
A _max	(β_A/α_A)- (β_R/α_A) ln(β_A/β_R)	(β₁/α) - K
R _max	α_A= α_R: [(β_A- β_a)/α-K_A]/e
	$α_{A} \neq α_{R} : (β_{A} / α_{A}) {(α_{R} / α_{A})}^{α_{R} / α_{A}}$
C _max	β_A/α_A
A _min	β_a/(α_A+ kR_max)	β₀/α
A_min, limiting	β_a/(kβ_A/αe)	β₀/α

^aLimiting expressions for the hysteresis oscillator are for α_A= α_R= α, β_A> β_R, β_A> αK_A.

^bLimiting expressions for the delay oscillator are for β₁/α > K > β₀/α.

The start of the cycle, t = 0, is defined to occur when A and R are both low with A increasing just past R. From the dynamics of A, $\dot{A}$ ≈ β_a- α_AA - kAR, and the nullcline $\dot{A}$ = 0 crosses A = R at the value $\sqrt{β_{a} / k + {(α_{A} / 2 k)}^{2}}$ . In the limit of a fast bimolecular reaction, $k > α_{A}^{2} / 4 β_{a}$ , the crossing point is at A = R = $\sqrt{β_{a} / k}$ . The first phase of dynamics ('activation phase') ends when A has risen to its maximum and then returned to a low value, with C high and R still small. The end of this first phase, with duration τ₁, is defined when C is at its maximum. In the second phase ('recovery phase'), C declines to a baseline value and R rises and falls. The end of the second phase, with duration τ₂, is defined when R has just crossed below A.

During the activation phase,

\begin{array}{l} (d / d t) (A + C) & \approx & β_{a} + (β_{A} - β_{a}) Θ (A - K_{A}) - α_{A} (A + C) \\ (d / d t) (R + C) & \approx & β_{r} + (β_{R} - β_{r}) Θ (A - K_{R}), \end{array}

(6)

where Θ(u) is the unit step function. Starting from A = R =

\sqrt{β_{a} / k}

, A rises quickly to K_A. The time required is approximately K_A/β_a(see Methods for a more accurate expression). The subsequent dynamics are

\begin{array}{l} A + C & \approx & β_{A} / α_{A} + [(β_{a} - β_{A}) / α_{A}] \exp (- α_{A} t) \\ R + C & \approx & β_{R} [t - (K_{R} - K_{A}) / β_{A}], \end{array}

(7)

neglecting β_rrelative to β_R. Since R ≈ 0 in the activation phase, R + C ≈ C, and (A + C) - (R + C) > K_Ais required to maintain activation. Assuming that A + C is close to its asymptotic value of β_A/α_Aat this point implies R + C = β_A/α_A- K_A, or

t = (β_A- α_AK_A)/α_Aβ_R+ (K_R- K_A)/β_A

for the elapsed time.

There is another brief transient in which A decays to 0, described using the combinations

\begin{array}{l} (d / d t) (C + R) & \approx & - α_{R} R \approx 0 \approx \dot{C} \\ (d / d t) (A - R) & \approx & β_{a} - α_{A} (A + C) - α_{R} R \approx β_{a} - β_{A} \approx \dot{A} . \end{array}

(9)

The second equation uses C ≈ β_A/α_Aduring this interval. The time for this transient is approximately K_A/(β_A- β₀).

The total time for the activation phase is therefore

τ₁ = (β_A- α_AK_A)/α_Aβ_R+ (K_R- K_A)/β_A+ K_A/(β_A- β_a) ≈ (β_A/α_A)/β_R

This value can be rationalized as the amount of time required for enough repressor to be produced, at rate β_R, to neutralize the total amount of activator both free and in complex, β_A/α_A.

At the start of the recovery phase, the complex is at its maximum concentration of approximately (β_A/α_A) - K_A≈ β_A/α_A. Here A ≈ 0 and it is convenient to examine the combinations C + A and R – A with

\begin{matrix} (d / d t) (C + A) = β_{a} - α_{A} (C + A) \approx β_{a} - α_{A} C \approx \dot{C} \\ (d / d t) (R - A) = - β_{a} + α_{A} (C + A) - α_{R} R \approx \dot{R} . \end{matrix}

(11)

For α_A≠ α_R, these equations give the dynamics

\begin{array}{l} C (t) & = & β_{a} / α_{a} + [(β_{A} - β_{a}) / α_{A}] \exp (- α_{A} t) \\ R (t) & = & [β_{A} / (α_{A} - α_{R})] \cdot [\exp (- α_{R} t) - \exp (- α_{A} t)], \end{array}

(12)

which continue until the concentration of R dips below A to trigger a new cycle. The concentrations cross at

\sqrt{β_{a} / k}

, and the duration of the recovery phase is

τ_{2} = {[\min (α_{A}, α_{R})]}^{- 1} \ln [(β_{A} / | α_{A} - α_{R} |) / \sqrt{β_{a} / k}] .

(13)

Comparison with the ODE solutions

A three-dimensional visualization of the dynamics (Fig. 3) demonstrates that the analytical expressions are in excellent agreement with numerical results from ODEs. For a more complete examination of agreement across parameter space, parameters representing decay rates (Fig. 4), activated production rates (Fig. 5), and baseline production rates (Fig. 6) were scanned over an order of magnitude. The period and the time for the individual phases are compared with numerical ODE results using 4D contour plots [17, 18].

Figure 3
**Limit cycle from analytical expressions (circles) and numerical calculations (triangles) are displayed in three dimensions**. Parameter values are from Table 1. Points are spaced at 200 equal time increments along the trajectory, and colors of points represent time progression in ROYGBIV order. The total period is 27.4 hr from the numerical ODE solution and 32.0 hr from the analytical expression.

Figure 4
**Analytic and numeric periods for** α_Aand α_R. Other parameter values are set to defaults from Table 1. The background is colored according to periods calculated from numerical solutions of the ODE model, and black contour lines are from analytical expressions in Table 2. **(A)** Full period, τ_tot. **(B)** Activation phase, τ₁. **(C)** Recovery phase, τ₂. In this and subsequent figures, contour lines are estimated from values calculated on a grid indicated by light background lines. A finer grid would yield smoother contours. The general agreement of the colored contours (ODEs) and the thick contour lines (analytic) indicates an accurate model for the dependence of the period on specific parameters. For example, the time τ₁ for the activation phase depends only on α_Aand not α_R, while the recovery phase depends on both. As the degradation rates decrease, the period increases rapidly. Decreasing the degradation rates further can eliminate oscillations.

Figure 5
**Analytic and numeric periods for** β_Aand β_R. Other parameter values are from Table 1, except that β_ais set to 0.1β_A. Plots are as in Fig. 4. The dependence on production rates is primarily in the activation phase, not the recovery phase.

Figure 6
**Analytic and numeric periods for** β_aand β_r. Plots are as in Fig. 4. The period depends primarily on β_aand is insensitive to β_r.

Decay rates have a strong influence on the period (Fig. 4). The activation phase depends almost entirely on α_A, with virtually no dependence on α_R; the recovery phase depends on the smaller of the two. The parameter space searched is symmetric about α_A= α_R, with α_A/α_Rranging from 1/30 to 30. Robust oscillations are still observed even when α_R>> α_A, demonstrating that oscillations do not require that A is the faster subsystem.

The dependence of the full period on the activated production rates scales roughly as β_A/β_R(Fig. 5). Most of this dependence arises from the activation phase. For the recovery phase, both the analytical and the numerical estimates suggest very little effect. This result is consistent with a low level of activator during the recovery phase. The baseline production rate of the activator does affect the time for the recovery phase, as well as the activation phase (Fig. 6). The production rate of the repressor is generally taken to be low in clock models, and over two orders of magnitude has virtually no effect on the dynamics.

Delay oscillator

The first well-known engineered biological clock was a delay oscillator termed the repressilator [6]. Equations for a standard simplified continuous, symmetric three-component model are presented (see Methods), again using the approximation that mRNA levels decay faster than protein levels. Repressilator dynamics, periods, and amplitudes are reviewed in the Methods and included in Table 2.

In the comparisons that follow for noise (defined by variance in the period), efficiency, and dynamic range, it is necessary to introduce a correspondence between parameters of the delay oscillator and hysteresis oscillator. We assume that production rates, variance in production rates, and decay rates are equivalent in the two systems.

Variance of the period

The noise in the oscillator period is analyzed here through the variance, Var(τ_tot), providing an analytical route to sensitivity analysis of robustness [19]. Assuming that production events of A and R are correlated but uncorrelated with decay events, and that both A and R decay at the same rate α, this variance is

\begin{array}{l} Var (τ_{tot}) & = & {(\partial τ_{tot} / \partial β_{A})}^{2} Var (β_{A}) + {(\partial τ_{tot} / \partial β_{B})}^{2} Var (β_{B}) + \\ (\partial τ_{tot} / \partial β_{A}) (\partial τ_{tot} / \partial β_{B}) Cov (β_{A}, β_{B}) + {(\partial τ_{tot} / \partial α)}^{2} Var (α) . \end{array}

(14)

Using the limiting form for the hysteresis oscillator period, Table 2,

\begin{array}{l} Var (τ_{tot}) & = & α^{- 2} [{(1 + β_{A} / β_{R})}^{2} Var (β_{A}) / β_{A}^{2} + {(β_{A} / β_{R})}^{2} Var (β_{R}) / β_{R}^{2} \\ - 2 (β_{A} / β_{R}) (1 + β_{A} / β_{R}) Cov (β_{A}, β_{R}) / β_{A} β_{R}] \\ + α^{- 2} Var (β_{a}) / 4 β_{a}^{2} + {[τ_{tot} + α^{- 1}]}^{2} Var (α) / α^{2} . \end{array}

(15)

To simplify this expression, we designated the correlation between β_Aand β_Ras r and assume that the coefficient of variation is roughly uniform for each component,

\begin{array}{l} Var (β_{A}) / β_{A}^{2} \approx Var (β_{R}) / β_{R}^{2} \approx Var (β_{a}) / β_{a}^{2} \approx σ^{2} : \\ \begin{array}{l} Var (τ_{tot}) & = & (σ^{2} / α^{2}) [(5 / 4) + 2 (β_{A} / β_{R}) (1 + β_{A} / β_{R}) (1 - r)] \\ + {[τ_{tot} + α^{- 1}]}^{2} Var (α) / α^{2} \end{array} \end{array}

(16)

The corresponding variance for the limiting form of the delay oscillator period is

Var (τ_{tot}) = (9 / α^{2}) Var (β_{1}) / β_{1}^{2} + {[τ_{tot} + 3 α^{- 1}]}^{2} Var (α) / α^{2}

(17)

For purposes of comparison, we assume that $Var (β_{1}) / β_{1}^{2} = σ^{2}$ as well, and that the decay rate α is also the same for the hysteresis and delay designs.

Consider separately the variance due to transcriptional noise, scaling as σ², and variance due to decay noise, scaling as Var(α). The variance due to transcriptional noise will be smaller for the hysteresis oscillator than the delay oscillator when

r ≥ 1 – 31/[8(β_A/β_R)(1 + β_A/β_R)].

The parameters in Table 1 suggest β_A/β_R≈ 5, and noise reduction for the hysteresis oscillator requires a rather large correlation in transcription rates, r ≥ 7/8. A smaller ratio gives a smaller correlation required for noise reduction. For example, β_A/β_R= 2 gives r ≥ 1/3. Correlation arises naturally because production rates of both the activator and the repressor depend on the concentration of free activator. For the true biological system, correlation would also arise from fluctuations in the concentrations of polymerase and ribosomes.

Variation in the period due to decay noise is always larger for the delay oscillator compared to the hysteresis oscillator. Assuming a period of 24 hours and a protein lifetime of 0.5 to 5 hours, the delay oscillator has 8% to 80% higher variance in its period than the hysteresis oscillator.

Efficiency

The efficiency of an oscillator is defined here as the inverse of its power requirements, where power is the rate of protein production averaged over a period. For the hysteresis oscillator, activator and represser molecules are produced at rates β_Aand β_Rduring the activation phase, and are not produced during the recovery phase, yielding the power requirement (β_A+ β_R)τ₁/τ_tot. For the delay oscillator, the synthesis rate during phase 1 is β₀ + 2β₁, and the synthesis rate during phase 2 is 2β₀ + β₁, with power requirement β₀ + β₁ + (β₁τ₁ + β₀τ₂)/τ_tot.

Assuming that baseline synthesis rates are small relative to activated synthesis rates and β_R/β_Ais no greater than 0.5 suggests these costs:

\begin{matrix} Hysteresis : & (3 / 2) β_{A} τ_{1} / τ_{tot} \\ Delay : & β_{1} . \end{matrix}

(19)

If the activated production rates are similar, β₁ ≈ β_A, then the hysteresis oscillator will have a cost of (3/2)τ₁/τ_tot relative to the delay oscillator, giving a cost advantage when τ₁/τ_tot < 2/3. Using typical parameters, the activation phase is faster than the recovery phase, with τ₁/τ_tot ≈ 1/3. This gives the hysteresis oscillator a two-fold cost reduction, or equivalently a two-fold efficiency increase relative to the delay oscillator.

Dynamic range

To be functional, an oscillator must couple to other biological components. The most straightforward coupling is for the activator molecule to serve as a transcription factor for output elements. These elements may have varying binding affinities for the activator, and it may therefore be advantageous for the activator to have a large dynamic range during a cycle. The dynamic range is quantified as decibels (dB) as 10 log₁₀(A_max/A_min).

Using the limiting forms from Table 2, the dynamic ranges of the oscillators are

\begin{matrix} Hysteresis : & 10 \log_{10} (β_{A} / β_{a}) + 10 \log_{10} (k β_{A} / e α^{2}) \\ Delay : & 10 \log_{10} (β_{1} / β_{0}) \end{matrix}

(20)

For the delay oscillator, the ratio of activated to baseline production rates is typically a factor of 10 to 100, yielding a 10 to 20 dB dynamic range. The hysteresis oscillator has a similar contribution of 10 dB from the ratio β_A/β_a. The hysteresis oscillator has an additional contribution, however, because the minimum concentration of A is much lower than the conventional steady-state baseline β_a/α. Again using typical values, kβ_A/eα² ≈ 7000 to 8000, and the effect is a boost of about 40 dB to the dynamic range.

Conclusion

This work provides a physical interpretation for the period and dynamic range of a model for a hysteresis oscillator. The period has a first phase whose duration, approximately (β_A/α_A)/β_R, can be interpreted as the time required to synthesize sufficient repressor molecules at rate β_Rto titrate an equilibrium concentration of activator molecules, β_A/α_A. The second phase has duration approximately equal to $α_{R}^{- 1} \ln [(β_{A} / α_{A}) / \sqrt{β_{a} / k}]$ . This has the familiar form of a protein lifetime, α^-1, multiplied by the log-ratio of an initial concentration to a final concentration. The initial concentration, β_A/α_A, corresponds to the same equilibrium concentration as before. The final concentration, $\sqrt{β_{a} / k}$ , is the value when the activator and repressor concentrations cross.

It is intriguing that the critical step triggering the start of a new phase is the crossing of the activator and repressor concentrations. In the context of gene expression, predictors based on crossing of mRNA abundances have been remarkably powerful for classification problems [20]. The results generated here for a particular nonlinear system show that such crossings can be important in marking the transition between distinct states.

The signal that oscillations are not supported is that the time for the second phase of the cycle, τ₂, becomes long. This is apparent in the contour plots showing the time for each phase. For example, in Fig. 4, when the decay constants become small, the period becomes rapidly larger. In this region, both A and R are small, and C is close to the value β_a/α_A. An expansion of the dynamical equations in this region could provide and analytically tractable expression for stability analysis.

Our results agree with numerical simulations that have found the hysteresis design to be more robust with respect to noise, but permit the ability to ascribe variance independently to production and decay sources. The hysteresis oscillator period is estimated to be roughly ten times more robust to fluctuations in decay rates. Reduction of noise from production rates requires positive correlations of at least 35% between activator and repressor production fluctuations. This observation is interesting because the hysteresis model explicitly couples the synthesis of these components during a single phase of the cycle, whereas the delay oscillator produces components continuously throughout the cycle. Same-time production of activator and repressor molecules should naturally introduce correlations in production rates because both depend on the same fluctuating concentration of activator molecules for transcriptional activation.

Our results also provide an intuitive explanation for a recent observation that coupled positive-negative feedback oscillators can cover a wider frequency spectrum than pure negative feedback oscillators [14]. For the period of the negative feedback oscillator considered here, the protein production rate appears inside a logarithm, giving it only weak influence on the period. For the hysteresis oscillator, however, protein production rates appear to linear order, with a much stronger ability to influence the period. Moreover, activator and repressor production rates appear as a ratio, permitting greater leverage for changing the period.

The analytical expressions also permit easy examination of other oscillator properties. Sustaining oscillations requires a cost that can be measured in the biomolecules that must be synthesized and then degraded over the course of a period. We estimate that the energetic cost of the hysteresis oscillator is about half that of delay oscillator. Our results suggest that efficiency may be an important oscillator property; to our knowledge, it has yet to be studied in computational models.

The output of an oscillator should have a large dynamic range to maximize its ability to couple to output systems. While the dynamic range of a delay oscillator is 10 to 20 dB, the dynamic range of a hysteresis oscillator with similar transcription and decay rates is 50 to 60 dB, an impressive gain. The large dynamic range arises from a state where all concentrations are close to 0. This does not necessarily make stochastic behavior important, however: the relevant count is the number of particles required to activate transcription, K_Aor K_R, rather than the small value

\sqrt{β_{a} / k}

marking the start of a new cycle. As seen in Fig. 2 with K_A= 50, trajectories from ODEs, analytical approximations, and stochastic simulations are quite similar. When stochastic simulations are run using parameters scaled by a factor of 50 to give K_A= 1, however, the period is shorter and more variable than continuous ODEs or the analytical formulas (Fig. 7). For the analytical and ODE trajectories, the only difference is a 50× scaling in the output amplitudes.

Figure 7
**Trajectories of activator, repressor, and complex concentrations are displayed for dynamics from analytical expressions (top panel), continuous ODEs (middle), and stochastic simulations (bottom)**. Parameter values are the default values from Table 1, essentially identical to those in Ref. [11] but with Hill coefficient = 2 and concentrations reduced by 50×.

In summary, these results provide a direct connection between parameters and observed properties of a circadian clock model. By showing how period and amplitude scale with parameters, these results help explain results observed in numerical simulations and suggest oscillator efficiency as an area where additional computational analysis may be valuable. The analysis strategy is to convert a nonlinear system into a series of linear systems connected at their boundaries, with a key transition marked by the crossing of activator and repressor concentrations. While the analysis here is for a particular clock model, the analysis strategy is general and should be applicable to other nonlinear biological systems.

Methods

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 D₁ and D₂ and effective radius h is

k = 4π(D₁ + D₂)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 cm²/sec. Using a nuclear diameter of roughly 7 μ m (volume = 1.8 × 10^-16 m³) 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_mwith dynamics

\begin{matrix} {\dot{X}}_{m} (t) = β^{'} (t) - α^{'} X_{m} \\ \dot{X} (t) = β^{″} X_{m} - α^{″} X, \end{matrix}

(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_mwith a constant coefficient β". Using the notation that

\tilde{f}

(s) is the Laplace transform of f(t),

\begin{matrix} {\tilde{X}}_{m} (s) = {(s + α^{'})}^{- 1} [X_{m} (0) + {\tilde{β}}^{'} (s)] \\ \tilde{X} (s) = {(s + α^{″})}^{- 1} X (0) + β^{″} {(s + α^{'})}^{- 1} {(s + α^{″})}^{- 1} [X_{m} (0) + {\tilde{β}}^{'} (x)] \end{matrix}

(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

\tilde{X} (s) = {(s + α^{″})}^{- 1} X (0) + (β^{″} / α^{'}) {(s + α^{″})}^{- 1} [1 - (s / α^{'}) + O {(s / α^{'})}^{2}] [X_{m} (0) = {\tilde{β}}^{'} (x)] .

(24)

The notation O(...) indicates asymptotic order. Transforming back to the time domain yields

\dot{X} (t) = (β^{″} / α^{'}) β^{'} (t) - α^{″} X (t) + O (α^{″} / α^{'})

(25)

with effective initial conditions

X_eff(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:

\begin{matrix} A & \overset{β_{a} + (β_{A} - β_{a}) {(A / K_{A})}^{n} / [1 + {(A / K_{A})}^{n}]}{\to} & A + 1 \\ A & \overset{α_{A} A}{\to} & A - 1 \\ R & \overset{β_{r} + (β_{R} - β_{r}) {(A / K_{R})}^{n} / [1 + {(A / K_{R})}^{n}]}{\to} & R + 1 \\ R & \overset{α_{R} R}{\to} & R - 1 \\ (A, R, C) & \overset{k A R}{\to} & (A - 1, R - 1, C + 1) \\ (R, C) & \overset{α_{A} C}{\to} & (R + 1, C - 1) \end{matrix}

(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 $\dot{C}$ = 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 $\dot{A}$ ≈ 0 or $\dot{R}$ ≈ 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

\dot{C}

= kAR - α_AC ≈ 0. In this regime, (d/dt)(C + A) = β_a- α_A(A + C) + (β_A- β_a)Aⁿ/[Aⁿ+

K_{A}^{n}

], 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/α_Awith transient time on the order of

α_{A}^{- 1}

. As A is small, the concentration C at time 0 is approximately β_a/α_A. Using

\dot{C}

≈ 0 and A = R,

X \equiv A (0) = R (0) = \sqrt{β_{a} / k}

(28)

defines the crossing value X.

As A becomes activated but is still below K_Aand K_R, the Hill functions are approximated to yield the equations

\begin{array}{l} \dot{A} & = & β_{a} + (β_{A} - β_{a}) {(A / K_{A})}^{n} - α_{A} A - k A R \\ \dot{R} & = & β_{r} + (β_{R} - β_{r}) {(A / K_{R})}^{n} - α_{R} R - k A R + α_{A} C \\ \dot{C} & = & k A R - α_{A} C . \end{array}

(29)

During this phase, A is much smaller than K_Aand K_R,

\dot{R}

is small, and α_RR and β_rare small relative to activated transcription of R, yielding the approximate dynamics

\begin{matrix} k A R = (β_{R} - β_{r}) {(A / K_{R})}^{n} + α_{A} C \\ \begin{matrix} \approx (β_{R} - β_{r}) {(A / K_{R})}^{n} + β_{a}, & and \end{matrix} \\ \dot{A} = [(β_{A} - β_{a}) / K_{A}^{n} - (β_{R} - β_{r}) / K_{R}^{n}] A^{n} - α_{A} A \\ \begin{matrix} \approx B_{X} A^{n}, & with \end{matrix} \\ β_{X} \equiv (β_{A} - β_{a}) / K_{A}^{n} - (β_{R} - β_{r}) / K_{R}^{n} . \end{matrix}

(30)

Recovery of A implies $\dot{A}$ > 0 and β_X> 0. The decay term α_AA 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_Kas the time when A = K_Ain this regime, Eq. 30 is integrated implicitly to yield

\begin{matrix} t (A) = \int_{X}^{A} d A {[β_{X} A^{n}]}^{- 1} = {[β_{X} (n - 1)]}^{- 1} [X^{- n + 1} - A^{- n + 1}], \\ t_{K} = {[B_{X} (n - 1)]}^{- 1} [X^{1 - n} - K_{A}^{1 - n}], \\ A (t) = {[X^{1 - n} - (n - 1) β_{X} t]}^{1 / (1 - n)} . \end{matrix}

(31)

In the case that n = 1, it is possible to retain the term α_AA in the analytical solution, and the corresponding equations are

\begin{matrix} t_{K} = {(β_{X} - α_{A})}^{- 1} \ln (K_{A} / X) \\ A (t) = X \exp [(β_{X} - α_{A}) t] . \end{matrix}

(32)

With typical parameters, β_X≈ β_A/K_A> α_A, and α_Acan be ignored as with the n > 1 case. This approximation ignores the very short transient in which A is between K_Aand 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_Aand rising rapidly, C ≈ β_a/α_A, and R ≈ 0. Here full transcriptional activation is a good approximation, A >> K_Aand A >> K_R, and the Hill functions become 1. The dynamics for R are $\dot{R}$ = β_R- α_RR - kAR + α_AC. Since $\dot{R}$ and R are both small, kAR ≈ β_R+ α_AC. The dynamics for C are therefore $\dot{C}$ = β_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

\begin{matrix} \dot{A} = β_{A} - α_{A} A - β_{R} - α_{A} C \\ = (β_{A} - β_{a} - β_{R}) - α_{A} A - α_{A} β_{R} (t - t_{K}) \\ A (t) = K_{A} \exp [- α_{A} (t - t_{K})] + [(β_{A} - β_{a}) / α_{A}] \cdot {1 - \exp [- α_{A} (t - t_{K})]} - β_{R} (t - t_{K}) . \end{matrix}

(34)

The time when A reaches its maximum and the maximum value are

\begin{matrix} t_{A} = t_{K} + (1 / α_{A}) \ln [(β_{A} - β_{a} - α_{A} K_{A}) / β_{n}] \\ A_{\max} = (β_{A} - β_{a} - β_{R}) / α_{A} - (β_{R} / α_{A}) \ln [(β_{A} - β_{a} - α_{A} K_{A}) / β_{R}] . \end{matrix}

(35)

These equations continue to hold until A returns to value K_A. The time t'_Kwhen this occurs is obtained from Eq. 34 using the approximation that the exponential transient has decayed and that β_ris negligible,

\begin{matrix} {t^{'}}_{K} = t_{K} + (β_{A} - β_{a} - α_{A} K_{A}) / α_{A} β_{R}, & and \\ C ({t^{'}}_{K}) = (β_{A} / α_{A}) - K_{A} \equiv C_{\max} \end{matrix}

(36)

from Eq. 33. The complex will be seen to have obtained its maximum value C_max 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

\dot{R}

≈ 0 and R ≈ 0 yield kAR = α_AC, giving

\dot{C}

≈ 0 as well. The dynamics for A in this regime are

\begin{matrix} \dot{A} = β_{a} - α_{A} A - β_{A} + α_{A} K_{A} \\ A (t) = K_{A} + [(β_{A} - β_{a}) / α_{A}] \cdot {[\exp (- α_{A} (t - {t^{'}}_{k})] - 1} . \end{matrix}

(37)

These dynamics continue until A = 0, at which point C is still equal to C_max. Defining this time as t_C,

t_{C} = {t^{'}}_{K} - (1 / α_{A}) \ln [1 - α_{A} K_{A} / (β_{A} - β_{a})] \approx {t^{'}}_{K} + K_{A} / (β_{A} - β_{a}) .

(38)

The repressor makes a round trip

In this phase of the cycle, A becomes the slow subsystem because A and hence

\dot{A}

are small, while R and C change more rapidly. The approximations

\dot{A}

= 0 with α_AA = 0 yield kAR = β_a, and

\begin{array}{l} \dot{C} = β_{a} - α_{A} C \\ \dot{R} = β_{r} + α_{A} C - α_{R} R - β_{a} . \end{array}

(39)

The trajectories are

\begin{array}{l} C (t) & = & β_{a} / α_{A} + [(β_{A} - β_{a}) / α_{A} - K_{A}] \exp [- α_{A} (t - t_{C})] \\ R (t) & = & [(β_{A} - β_{a} - α_{A} K_{A}) / (α_{A} - α_{R})] \cdot {\exp [- α_{R} (t - t_{C})] - \exp [- α_{A} (t - t_{C})]} \\ + (β_{r} / α_{R}) \cdot {1 - \exp [- α_{R} (t - t_{C})]}, or \end{array}

(40)

\begin{array}{l} R (t) & = & (β_{A} - β_{a} - α_{A} K_{A}) \cdot (t - t_{C}) \cdot \exp [- α (t - t_{C})] \\ + (β_{r} / α) \cdot {1 - \exp [- α (t - t_{C})]} \end{array}

(41)

in the case that α_A= α_R= α.

Ignoring the small contribution from the baseline production rate β_r, the time t_Rfor the maximum of R and its maximum value are

\begin{matrix} t_{R} = t_{C} + [1 / (α_{A} - α_{R})] \ln (α_{A} / α_{R}), \\ R_{\max} = [(β_{A} - β_{a}) / α_{A} - K_{A}] \cdot {(α_{R} / α_{A})}^{α_{R} / (α_{A} - α_{R})}; or \\ t_{R} = t_{C} + α^{- 1}, \\ R_{\max} = [(β_{A} - β_{a}) / α_{A} - K_{A}] \cdot e^{- 1} \end{matrix}

(42)

in the case that α_A= α_R= α.

An improved approximation for A in this region, again from

\dot{A}

= 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 =

\sqrt{β_{a} / k}

, 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_Xis calculated as

\begin{array}{l} t_{X} & = & t_{C} + α_{\min}^{- 1} \ln {[(β_{A} - β_{a} - α_{A} K_{A}) / Δ α \sqrt{β_{a} / k}]} \\ + α_{\min}^{- 1} \ln {1 - {[Δ α \sqrt{β_{a} / k} / (β_{A} - β_{a} - α_{A} K_{A})]}^{Δ α / α_{\min}}} . \end{array}

(43)

The final term in t_Xis 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- α_AK_Aand

Δ X = \sqrt{β_{a} / k} - (β_{r} / α_{R})

, and assuming that β_r/α <<β_Aτ_X,

\begin{array}{l} τ_{X} & = & α^{- 1} \ln [(Δ β / Δ X) τ_{X}] \\ = & α^{- 1} \ln [(Δ β / Δ X) α^{- 1} \ln [(Δ β / Δ X) τ_{X}] \\ = & α^{- 1} \ln [Δ β / α Δ X] + α^{- 1} \ln {\ln [(Δ β / Δ X) τ_{X}]} \\ \approx & α^{- 1} \ln [Δ β / α Δ X] + α^{- 1} \ln {\ln [(Δ β / α Δ X) \ln (Δ β / α Δ X)]}, \end{array}

(44)

where in the last equation the term τ_Xon 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

\begin{matrix} τ_{K} = t_{K} \\ {τ^{'}}_{K} = {t^{'}}_{K} - t_{K} \\ τ_{C} = t_{C} - {t^{'}}_{K} \\ τ_{X} = t_{X} - t_{C} . \\ τ_{tot} = t_{X} = τ_{K} + {τ^{'}}_{K} + τ_{C} + τ_{X} . \end{matrix}

(45)

Simplified formulas for the intervals are obtained by using leading-order contributions based on biologically reasonable assumptions that threshold values (K_Aand K_R) are small compared to the activated level of A, and that the fully activated production rate for A, β_Ais much larger than the baseline rate β_aor the fully activated production rate β_Rof repressor, and that β_ris 0.

In these limits,

\begin{matrix} τ_{K} = {[β_{A} (n - 1)]}^{- 1} K_{A}^{n} {(β_{a} / k)}^{1 - n}, or \\ = (K_{A} / β_{A}) \ln (K_{A} / \sqrt{β_{a} / k}) if n = 1 \\ {τ^{'}}_{K} = α_{A}^{- 1} (β_{A} / β_{R}) \\ τ_{C} = K_{A} / β_{A} \\ τ_{X} = α_{\min}^{- 1} \ln [β_{A} / | α_{A} - α_{R} | \sqrt{β_{a} / k}] \\ = α^{- 1} \ln [β_{A} / α \sqrt{β_{a} / k}] if α_{A} = α_{R} . \end{matrix}

(46)

Both τ_Kand τ_Care fast, scaling as inverse production rates, while τ'_Kand τ_Xare slow, scaling as inverse decay rates. To compare with numeric simulation results, we break the entire cycle into two phases with times τ₁ and τ₂. The first phase begins at the start of the cycle with A = R and

\dot{A}

> 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:

\begin{array}{l} τ_{1} & = & τ_{K} + {τ^{'}}_{K} + τ_{C} \approx {τ^{'}}_{K} \approx α_{A}^{- 1} β_{A} / β_{R} \\ τ_{2} & = & τ_{X} \approx α_{\min}^{- 1} \ln [β_{A} / | α_{A} - α_{R} | \sqrt{β_{a} / k}] \end{array}

(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 β_Aand repressor proteins at rate β_R. During the second phase, synthesis rates are β_aand β_r. The cost per cycle is

Cost = (β_A+ β_R)τ₁/τ_tot + (β_a+ β_r)τ₂/τ_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_irepresents one of i = 1 to 3 components, and P_i-1= P₃ when i = 1.

Using the same approximations as for the hysteresis oscillator, the repressilator cycle starts when P₃ falls just below K, with P₁ at its baseline level β₀/α and P₂ roughly equal to (β₁/α) - K. In this regime,

These dynamics continue until P₁(t) = K, defined to occur after an interval τ₁,

τ₁ = α^-1 ln [(β₁ - β₀)/(β₁ - α K)] ≈ K/β₁,

where the approximation holds when β₁/α >> K > β₀/α. At this point,

\begin{array}{l} P_{1} (τ_{1}) & = & K \\ P_{2} (τ_{1}) & = & (β_{1} / α) - K (β_{1} - α K) / (β_{1} - β_{0}) \approx (β_{1} / α) - K \\ P_{3} (τ_{1}) & = & K (β_{1} - α K) / (β_{1} - β_{0}) + (β_{0} / α) (α K - β_{0}) / (β_{1} - β_{0}) \approx K \end{array}

(52)

with dynamics

\begin{array}{l} {\dot{P}}_{1} & = & β_{1} - α P_{1} \\ {\dot{P}}_{2} & = & β_{0} - α P_{2} \\ {\dot{P}}_{3} & = & β_{0} - α P_{3} . \end{array}

(53)

These dynamics continue until time τ₁ + τ₂ when P₂(τ₁ + τ₂) = K,

τ₂ = α^-1 ln [(β₁ - αK - β₀)/(αK - β₀)] ≈ α^-1 ln(β₁/αK).

This pattern repeats three times for a total period τ_tot of

τ_tot = 3(τ₁ + τ₂) = (3/α) ln [(β₁ - β₀)(β₁ - αK - β₀)/(β₁ - αK)(αK - β₀)].

Energy cost

During the τ₊ phase, components 1 and 2 are synthesized at rate β₁, while component 3 is synthesized at rate β₀. During the τ _ phase, only component 1 is synthesized at the higher rate. The energy consumption, in terms of proteins per cycle, is

\begin{array}{l} Cost & = & β_{0} + β_{1} + (β_{1} τ_{1} + β_{0} τ_{2}) / τ_{tot} \\ \approx & β_{1} + 2 β_{0} + K / τ_{2} . \end{array}

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

Declarations

Acknowledgements

This material is based upon work supported by the National Science Foundation under Grant No. NSF CAREER 0546446. JSB acknowledges additional support from the Whitaker Foundation, helpful comments from students in the Fall 2007 semester of JHU 580.420, and assistance from teaching assistants Mr. Hailiang Huang, Mr. Scott Patterson, Ms. Yan Qi, and Mr. Yasir Suhail.

Electronic supplementary material

12918_2008_274_MOESM1_ESM.r Additional file 1: clocksim.r. Source code for analytical and numerical simulations and plots. (R 24 KB)

Below are the links to the authors’ original submitted files for images.

12918_2008_274_MOESM2_ESM.ppt Authors’ original file for figure 1

12918_2008_274_MOESM3_ESM.tiff Authors’ original file for figure 2

12918_2008_274_MOESM4_ESM.pdf Authors’ original file for figure 3

12918_2008_274_MOESM5_ESM.eps Authors’ original file for figure 4

12918_2008_274_MOESM6_ESM.eps Authors’ original file for figure 5

12918_2008_274_MOESM7_ESM.eps Authors’ original file for figure 6

12918_2008_274_MOESM8_ESM.pdf Authors’ original file for figure 7

Authors’ Affiliations

(1)

Department of Biomedical Engineering, Johns Hopkins University School of Medicine, Baltimore, MD 21205, USA

(2)

High-Throughput Biology Center, Johns Hopkins University School of Medicine, Baltimore, MD 21205, USA

(3)

Department of Biomedical Engineering, Johns Hopkins University School of Medicine, Baltimore, MD 21218, USA

References

Alon U: An introduction to systems biology : design principles of biological circuits. 2007, Chapman & Hall/CRCGoogle Scholar
Rosenfeld N, Alon U: Response delays and the structure of transcription networks. J Mol Biol. 2003, 329 (4): 645-654. 10.1016/S0022-2836(03)00506-0View ArticlePubMedGoogle Scholar
Kalir S, McClure J, Pabbaraju K, Southward C, Ronen M, Leibler S, Surette MG, Alon U: Ordering genes in a flagella pathway by analysis of expression kinetics from living bacteria. Science. 2001, 292 (5524): 2080-2083. 10.1126/science.1058758View ArticlePubMedGoogle Scholar
Collier JR, Monk NA, Maini PK, Lewis JH: Pattern formation by lateral inhibition with feedback: a mathematical model of delta-notch intercellular signalling. J Theor Biol. 1996, 183 (4): 429-446. 10.1006/jtbi.1996.0233View ArticlePubMedGoogle Scholar
Gardner TS, Cantor CR, Collins JJ: Construction of a genetic toggle switch in Escherichia coli. Nature. 2000, 403 (6767): 339-342. 10.1038/35002131View ArticlePubMedGoogle Scholar
Elowitz MB, Leibler S: A synthetic oscillatory network of transcriptional regulators. Nature. 2000, 403 (6767): 335-338. 10.1038/35002125View ArticlePubMedGoogle Scholar
Shetty RP, Endy D, Knight TF: Engineering BioBrick vectors from BioBrick parts. J Biol Eng. 2008, 2: 5- 10.1186/1754-1611-2-5PubMed CentralView ArticlePubMedGoogle Scholar
Goldbeter A: Biochemical oscillations and cellular rhythms : the molecular bases of periodic and chaotic behaviour. 1996, Cambridge University PressView ArticleGoogle Scholar
Young MW, Kay SA: Time zones: a comparative genetics of circadian clocks. Nat Rev Genet. 2001, 2 (9): 702-715. 10.1038/35088576View ArticlePubMedGoogle Scholar
Barkai N, Leibler S: Circadian clocks limited by noise. Nature. 2000, 403 (6767): 267-268.PubMedGoogle Scholar
Vilar JMG, Kueh HY, Barkai N, Leibler S: Mechanisms of noise-resistance in genetic oscillators. Proc Natl Acad Sci USA. 2002, 99 (9): 5988-5992. 10.1073/pnas.092133899PubMed CentralView ArticlePubMedGoogle Scholar
Liu AC, Welsh DK, Ko CH, Tran HG, Zhang EE, Priest AA, Buhr ED, Singer O, Meeker K, Verma IM, Doyle FJ, Takahashi JS, Kay SA: Intercellular coupling confers robustness against mutations in the SCN circadian clock network. Cell. 2007, 129 (3): 605-616. 10.1016/j.cell.2007.02.047PubMed CentralView ArticlePubMedGoogle Scholar
Guantes R, Poyatos JF: Dynamical principles of two-component genetic oscillators. PLoS Comput Biol. 2006, 2 (3): e30- 10.1371/journal.pcbi.0020030PubMed CentralView ArticlePubMedGoogle Scholar
Tsai TYC, Choi YS, Ma W, Pomerening JR, Tang C, Ferrell JE: Robust, tunable biological oscillations from interlinked positive and negative feedback loops. Science. 2008, 321 (5885): 126-129. 10.1126/science.1156951PubMed CentralView ArticlePubMedGoogle Scholar
McAdams HH, Arkin A: Stochastic mechanisms in gene expression. Proc Natl Acad Sci USA. 1997, 94 (3): 814-819. 10.1073/pnas.94.3.814PubMed CentralView ArticlePubMedGoogle Scholar
Kim J, Heslop-Harrison P, Postlethwaite I, Bates DG: Stochastic noise and synchronisation during dictyostelium aggregation make cAMP oscillations robust. PLoS Comput Biol. 2007, 3 (11): e218- 10.1371/journal.pcbi.0030218PubMed CentralView ArticlePubMedGoogle Scholar
Goodyear CP: Spawning stock biomass per recruit in fisheries management: foundation and current use. Canadian Special Publication of Fisheries and Aquatic Sciences. 1993, 120: 67-81.Google Scholar
Prager MH: 4D Contour Plots. http://addictedtor.free.fr/graphiques/RGraphGallery.php?graph=90
Stelling J, Gilles ED, Doyle FJ: Robustness properties of circadian clock architectures. Proc Natl Acad Sci USA. 2004, 101 (36): 13210-13215. 10.1073/pnas.0401463101PubMed CentralView ArticlePubMedGoogle Scholar
Tan AC, Naiman DQ, Xu L, Winslow RL, Geman D: Simple decision rules for classifying human cancers from gene expression profiles. Bioinformatics. 2005, 21 (20): 3896-3904. 10.1093/bioinformatics/bti631PubMed CentralView ArticlePubMedGoogle Scholar
Smolen P, Baxter DA, Byrne JH: A reduced model clarifies the role of feedback loops and time delays in the Drosophila circadian oscillator. Biophys J. 2002, 83 (5): 2349-2359. 10.1016/S0006-3495(02)75249-1PubMed CentralView ArticlePubMedGoogle Scholar
Leloup JC, Goldbeter A: Toward a detailed computational model for the mammalian circadian clock. Proc Natl Acad Sci USA. 2003, 100 (12): 7051-7056. 10.1073/pnas.1132112100PubMed CentralView ArticlePubMedGoogle Scholar
Ruoff P, Christensen MK, Sharma VK: PER/TIM-mediated amplification, gene dosage effects and temperature compensation in an interlocking-feedback loop model of the Drosophila circadian clock. J Theor Biol. 2005, 237: 41-57. 10.1016/j.jtbi.2005.03.030View ArticlePubMedGoogle Scholar
Leise TL, Moin EE: A mathematical model of the Drosophila circadian clock with emphasis on posttranslational mechanisms. J Theor Biol. 2007, 248: 48-63. 10.1016/j.jtbi.2007.04.013View ArticlePubMedGoogle Scholar
Xie Z, Kulasiri D: Modelling of circadian rhythms in Drosophila incorporating the interlocked PER/TIM and VRI/PDP1 feedback loops. J Theor Biol. 2007, 245 (2): 290-304. 10.1016/j.jtbi.2006.10.028View ArticlePubMedGoogle Scholar
Berg HC: Random Walks in Biology. 1983, Princeton University PressGoogle Scholar
R Development Core Team: R: A Language and Environment for Statistical Computing. 2006, R Foundation for Statistical Computing, Vienna, Austria, http://www.R-project.orgGoogle Scholar
Hindmarsh A: ODEPACK, A Systematized Collection of ODE Solvers. Scientific Computing: Applications of Mathematics and Computing to the Physical Sciences, IMACS Transactions on Scientific Computing. 1983, 1: 55-64. Stepleman Rea, Amsterdam, Netherlands; New York, U.S.A.: North-HollandGoogle Scholar
Petzold L: Automatic Selection of Methods for Solving Stiff and Nonstiff Systems of Ordinary Differential Equations. SIAM Journal on Scientific and Statistical Computing. 1983, 4: 136-148. 10.1137/0904010.View ArticleGoogle Scholar
Gillespie DT: A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions. Journal of Computational Physics. 1976, 22: 403-10.1016/0021-9991(76)90041-3.View ArticleGoogle Scholar
Gillespie DT: Exact stochastic simulation of coupled chemical reactions. Journal of Physical Chemistry. 1977, 81 (25): 2340-2361. 10.1021/j100540a008.View ArticleGoogle Scholar
Pineda-Krch M: GillespieSSA: Implementing the Gillespie Stochastic Simulation Algorithm in R. Journal of Statistical Software. 2008, 25 (12): 1-18. http://www.jstatsoft.org/v25/i12View ArticleGoogle Scholar
Tian T, Burrage K: Binomial leap methods for simulating stochastic chemical kinetics. J Chem Phys. 2004, 121 (21): 10356-10364. http://dx.doi.org/10.1063/1.1810475 10.1063/1.1810475View ArticlePubMedGoogle Scholar
Chatterjee A, Vlachos DG, Katsoulakis MA: Binomial distribution based tau-leap accelerated stochastic simulation. J Chem Phys. 2005, 122 (2): 024112- 10.1063/1.1833357View ArticlePubMedGoogle Scholar
Munsky B, Khammash M: The finite state projection algorithm for the solution of the chemical master equation. J Chem Phys. 2006, 124 (4): 044104- 10.1063/1.2145882View ArticlePubMedGoogle Scholar

Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Analytical approximations for the amplitude and period of a relaxation oscillator

Background

Results

Conclusion

Hysteresis oscillator model

Analytical expressions for period and phase

Comparison with the ODE solutions

Delay oscillator

Variance of the period

Efficiency

Dynamic range

Parameter values

Justification of steady-state approximations for mRNA

ODE calculations

Stochastic simulations

Analytical approximations

Recovery from time 0 to full activation at time t K

The activator makes a round trip

The repressor makes a round trip

Time intervals

Energy cost

Repressilator dynamics

Energy cost

Implementation and availability

Acknowledgements

BMC Systems Biology

Contact us

Recovery from time 0 to full activation at time t_K