Formula display:

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.

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

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.

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

The parameters and are baseline transcriptional rates for A_m and R_m; and 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 and 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, ≈ 0 and ≈ 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).

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, ≈ β_a - α_AA - kAR, and the nullcline = 0 crosses A = R at the value . In the limit of a fast bimolecular reaction, , the crossing point is at A = R = . 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 = , 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

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

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

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

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.

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

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

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

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 = 100 molecules/hr. In comparison, the maximum production rates of A and R are expected to be only 10 to 20 molecules/hr for typical parameters, and decay rates are slower still.

With this timescale separation, the bimolecular reaction effectively goes to completion: all available molecules of A and R are paired into complexes, until either A ≈ 0 or R ≈ 0. Depending on whether A or R is limiting, this leads to the approximation that ≈ 0 or ≈ 0. These approximations lead to analytical expressions, as detailed below.

Recovery from time 0 to full activation at time t_K

The start of the cycle, time 0, is defined when A and R are small, with A just crossing above R, and = kAR - α_AC ≈ 0. In this regime, (d/dt)(C + A) = β_a - α_A(A + C) + (β_A - β_a)Aⁿ/[Aⁿ + ], which is approximately equal to β_a - α_A(A + C) for A <<K_A. At the end of the previous cycle, the sum A + C therefore approached an equilibrium value of β_a/α_A with transient time on the order of . As A is small, the concentration C at time 0 is approximately β_a/α_A. Using ≈ 0 and A = R,

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, is small, and α_RR and β_r are small relative to activated transcription of R, yielding the approximate dynamics

Recovery of A implies > 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_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 α_AA 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 = β_R - α_RR - kAR + α_AC. Since and R are both small, kAR ≈ β_R + α_AC. The dynamics for C are therefore = β_R, or

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 ≈ 0 and R ≈ 0 yield kAR = α_AC, giving ≈ 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 are small, while R and C change more rapidly. The approximations = 0 with α_AA = 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 = 0, is A(t) ≈ β_a/[α_A + kR(t)]. According to this expression, A has its minimum value when R is at its maximum. These dynamics continue until R = X = , at which point R ≈ A and the cycle restarts. For convenience, let α_min = min(α_A, α_R) and Δα = |α_A - α_R|. Again ignoring the small contribution from β_r, the crossing time t_X is calculated as

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 - α_AK_A and , 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 > 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

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

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

with dynamics

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

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

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

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.

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.

References