Analytical approximations for the amplitude and period of a relaxation oscillator

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.

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

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.

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 β _a and β _a , with total production rising to β _A and β _R when the concentration of free activator is high. The activator degrades at rate α _A whether free or complexed; the repressor degrades at rate α _R but only when free.

Full size image

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 α _A and α _R , the decay rates of the activator and repressor proteins; β _A and β _R , the fully activated production rates; β _a , the baseline production rate of the activator, and k is the bimolecular rate constant for A+R\underset{}{\overset{k}{\to }}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{{\beta }_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 {\alpha }_R^{-1}\ln [({\beta }_A/{\alpha }_A)/\sqrt{{\beta }_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.

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 _m and R _m . Continuous concentrations are assumed throughout. The mathematical model is

The parameters {{\beta }^{\prime }}_a and {{\beta }^{\prime }}_r are baseline transcriptional rates for A _m and R _m ; {{\beta }^{\prime }}_A and {{\beta }^{\prime }}_R are the fully activated transcriptional rates. Transcriptional activation is represented by Hill functions with half-response at A = K _A for the activator and A = K _R for 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 {{\beta }^{\prime }}_A and {{\beta }^{\prime }}_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,

and effective equations

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

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.

Full size image

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 _A of 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.

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 - α _A A - kAR, and the nullcline \dot{A} = 0 crosses A = R at the value \sqrt{{\beta }_a/k+{({\alpha }_A/2k)}^2}. In the limit of a fast bimolecular reaction, k>{\alpha }_A^2/4{\beta }_a, the crossing point is at A = R = \sqrt{{\beta }_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,

where Θ(u) is the unit step function. Starting from A = R = \sqrt{{\beta }_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

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

t = (β _A - α _A K _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

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

The total time for the activation phase is therefore

τ₁ = (β _A - α _A K _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

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

which continue until the concentration of R dips below A to trigger a new cycle. The concentrations cross at \sqrt{{\beta }_a/k}, and the duration of the recovery phase is

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

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.

Full size image

Analytic and numeric periods for α _A and α _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 α _A and 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.

Full size image

Analytic and numeric periods for β _A and β _R . Other parameter values are from Table 1, except that β _a is 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.

Full size image

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

Full size image

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 /α _R ranging 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

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

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

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

For purposes of comparison, we assume that \text{Var}({\beta }_1)/{\beta }_1^2={\sigma }^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 β _A and β _R during 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 /β _A is no greater than 0.5 suggests these costs:

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

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.

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 β _R to titrate an equilibrium concentration of activator molecules, β _A /α _A . The second phase has duration approximately equal to {\alpha }_R^{-1}\ln [({\beta }_A/{\alpha }_A)/\sqrt{{\beta }_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{{\beta }_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 _A or K _R , rather than the small value \sqrt{{\beta }_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.

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

Full size image

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.

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 _m with dynamics

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 \tilde{f}(s) is the Laplace transform of f(t),

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

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

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:

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 - α _A C ≈ 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 /α _A with transient time on the order of {\alpha }_A^{-1}. As A is small, the concentration C at time 0 is approximately β _a /α _A . Using \dot{C} ≈ 0 and A = R,

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

During this phase, A is much smaller than K _A and K _R , \dot{R} is small, and α _R R and β _r are small relative to activated transcription of R, yielding the approximate dynamics

Recovery of A implies \dot{A} > 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

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

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 \dot{R} = β _R - α _R R - kAR + α _A C. Since \dot{R} and R are both small, kAR ≈ β _R + α _A C. 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

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

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,

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 = α _A C, giving \dot{C} ≈ 0 as well. The dynamics for A in this regime are

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

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 α _A A = 0 yield kAR = β _a , and

The trajectories are

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

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{{\beta }_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 _X is calculated as

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 \Delta X=\sqrt{{\beta }_a/k}-({\beta }_r/{\alpha }_R), and assuming that β _r /α <<β _A τ _X ,

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

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,

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

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 )τ₁/τ_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 _i represents 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,

with dynamics

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

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.

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.

Affiliations

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

   • Carmen Kut
   • , Vahid Golkhou
   •  & Joel S Bader

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

   • Vahid Golkhou
   •  & Joel S Bader

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

   • Joel S Bader

Corresponding author

Correspondence to Joel S Bader.

Authors' contributions

CK and VG assisted with the mathematical analysis and performed the numerical simulations. JSB conceived the modeling strategy, designed the study, and performed some of the simulations. All three participated in writing the manuscript.

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

Reprints and Permissions