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.


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: 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 Schematic oscillator designs Figure 1 Schematic oscillator designs. (Left) The repressilator is a delay oscillator here idealized as three symmetric repressive components, P 1 , P 2 , and P 3 , 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.
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. 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.

Results and discussion
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 employ-ing 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 halfmaximum 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 nonlinear 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 -α A A -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 τ 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.
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 for the elapsed time.
There is another brief transient in which A decays to 0, described using the combinations ( / )( )  . Parameter values are exactly as in Ref. [11], provided in Table 1, except with Hill coefficient = 2 to reflect cooperative dimer binding.

The second equation uses
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 At the start of the recovery phase, the complex is at its Here A ≈ 0 and it is convenient to examine the com- For α A ≠ α R , these equations give the dynamics Limiting expressions for the delay oscillator are for Limit cycle from analytical expressions (circles) and numeri-cal calculations (triangles) are displayed in three dimensions 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 param-eters 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 σ 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 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 reduc- Var Var Var  Analytic and numeric periods for β A and β R

Figure 5
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. Analytic and numeric periods for β a and β r Figure 6 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 . tion. 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 .
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, β 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 ).
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α 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 . 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, τ 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.
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 approxi-mations, 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.

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  Table 1, essentially identical to those in Ref. [11] but with Hill coefficient = 2 and concentrations reduced by 50×. 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][22][23][24][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 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
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, C 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 , 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 α R R 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 α 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 = β R -α R R -kAR + α A C. Since and R are both small, kAR 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. 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 . 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 α A K A and , and assuming that 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. 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 τ 1 and τ 2 . The first phase begins at the start of the cycle with A = R and > 0 and ends when C is at its maximum. The second phase begins when C is at its maximum and ends at the start of the next cycle:

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

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.