Analytical approximations for the amplitude and period of a relaxation oscillator
- Carmen Kut^{1},
- Vahid Golkhou^{1, 2} and
- Joel S Bader^{1, 2, 3}Email author
https://doi.org/10.1186/1752-0509-3-6
© Kut et al; licensee BioMed Central Ltd. 2009
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
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.
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].
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}\mathrm{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.
Results and discussion
Hysteresis oscillator model
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).
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.
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^{2} |
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 |
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.
Analytical results for oscillator period and phase
Times | Hysteresis^{ a } | Delay^{ b } |
---|---|---|
τ_{1} (A > R) | n = 1: (β_{ A }/α_{ A }β_{ R }) + K_{ A }/β_{ A }+(K_{ A }/β_{ A }) ln(K_{ A }/$\sqrt{{\beta}_{a}/k}$) | α^{-1} ln [(β_{1} - β_{0})/(β_{1} - αK)] |
n ≠ 1: (β_{ A }/α_{ A }β_{ R }) + K_{ A }/β_{ A }+ [β_{ A }(n- 1)]^{-1} ${K}_{A}^{n}$(β_{ a }/k)^{1-n} | ||
τ_{2} (A <R) | α_{ A }= α_{ R }: α^{-1} ln [(β_{ A }/α)/$\sqrt{{\beta}_{a}/k}$] | α^{-1} ln [(β_{1} - αK - β_{0})/(αK - β_{0})] |
α_{ A }≠ α_{ R }: ${\alpha}_{\mathrm{min}}^{-1}$ ln [(β_{ A }/|α_{ A }- α_{ R }|)/$\sqrt{{\beta}_{a}/k}$] | ||
τ _{tot} | τ_{1} + τ_{2} | 3(τ_{1} + τ_{2}) |
τ_{tot}, limiting | α^{-1} [β_{ A }/β_{ R }+ ln(β_{ A }/α $\sqrt{{\beta}_{a}/k}$)] | 3α^{-1} ln [β_{1}/αK] |
Concentrations | ||
A _{max} | (β_{ A }/α_{ A })- (β_{ R }/α_{ A }) ln(β_{ A }/β_{ R }) | (β_{1}/α) - K |
R _{max} | α_{ A }= α_{ R }: [(β_{ A }- β_{ a })/α-K_{ A }]/e | |
${\alpha}_{A}\ne {\alpha}_{R}:({\beta}_{A}/{\alpha}_{A}){({\alpha}_{R}/{\alpha}_{A})}^{{\alpha}_{R}/{\alpha}_{A}}$ | ||
C _{max} | β_{ A }/α_{ A } | |
A _{min} | β_{ a }/(α_{ A }+ kR_{max}) | β_{0}/α |
A_{min}, limiting | β_{ a }/(kβ_{ A }/αe) | β_{0}/α |
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 τ_{1}, 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 τ_{2}, is defined when R has just crossed below A.
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.
The second equation uses C ≈ β_{ A }/α_{ A }during this interval. The time for this transient is approximately K_{ A }/(β_{ A }- β_{0}).
The total time for the activation phase is therefore
τ_{1} = (β_{ 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 }.
Comparison with the ODE solutions
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
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 σ^{2}, 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 })τ_{1}/τ_{tot}. For the delay oscillator, the synthesis rate during phase 1 is β_{0} + 2β_{1}, and the synthesis rate during phase 2 is 2β_{0} + β_{1}, with power requirement β_{0} + β_{1} + (β_{1}τ_{1} + β_{0}τ_{2})/τ_{tot}.
If the activated production rates are similar, β_{1} ≈ β_{ A }, then the hysteresis oscillator will have a cost of (3/2)τ_{1}/τ_{tot} relative to the delay oscillator, giving a cost advantage when τ_{1}/τ_{tot} < 2/3. Using typical parameters, the activation phase is faster than the recovery phase, with τ_{1}/τ_{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_{10}(A_{max}/A_{min}).
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α^{2} ≈ 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 β_{ R }to titrate an equilibrium concentration of activator molecules, β_{ A }/α_{ A }. The second phase has duration approximately equal to ${\alpha}_{R}^{-1}\mathrm{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, τ_{2}, 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.
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_{1} and D_{2} and effective radius h is
k = 4π(D_{1} + D_{2})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^{2}/sec. Using a nuclear diameter of roughly 7 μ m (volume = 1.8 × 10^{-16} m^{3}) 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
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
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 }
defines the crossing value X.
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 }.
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 }.
from Eq. 33. The complex will be seen to have obtained its maximum value C_{max} at this time.
The repressor makes a round trip
in the case that α_{ A }= α_{ R }= α.
in the case that α_{ A }= α_{ R }= α.
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
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.
Energy cost
The energy cost of sustaining oscillations is estimated as the number of protein molecules synthesized per cycle, giving equal weight to activator and repressor proteins. During the first phase, activator proteins are synthesized at rate β_{ A }and repressor proteins at rate β_{ R }. During the second phase, synthesis rates are β_{ a }and β_{ r }. The cost per cycle is
Cost = (β_{ A }+ β_{ R })τ_{1}/τ_{tot} + (β_{ a }+ β_{ r })τ_{2}/τ_{tot}.
Repressilator dynamics
The repressilator is a synthetic oscillator constructed from an odd cycle of negative feedback loops [6].
Again using a steady-state approximation for mRNA levels, the dynamics for the standard three-component symmetric continuous repressilator are
where P_{ i }represents one of i = 1 to 3 components, and P_{i-1}= P_{3} when i = 1.
Using the same approximations as for the hysteresis oscillator, the repressilator cycle starts when P_{3} falls just below K, with P_{1} at its baseline level β_{0}/α and P_{2} roughly equal to (β_{1}/α) - K. In this regime,
These dynamics continue until P_{1}(t) = K, defined to occur after an interval τ_{1},
τ_{1} = α^{-1} ln [(β_{1} - β_{0})/(β_{1} - α K)] ≈ K/β_{1},
These dynamics continue until time τ_{1} + τ_{2} when P_{2}(τ_{1} + τ_{2}) = K,
τ_{2} = α^{-1} ln [(β_{1} - αK - β_{0})/(αK - β_{0})] ≈ α^{-1} ln(β_{1}/αK).
This pattern repeats three times for a total period τ_{tot} of
τ_{tot} = 3(τ_{1} + τ_{2}) = (3/α) ln [(β_{1} - β_{0})(β_{1} - αK - β_{0})/(β_{1} - αK)(αK - β_{0})].
Energy cost
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.
Authors’ Affiliations
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.