Mathematical analysis of robustness of oscillations in models of the mammalian circadian clock

Circadian rhythms in a wide range of organisms are mediated by molecular mechanisms based on transcription-translation feedback. In this paper, we use bifurcation theory to explore mathematical models of genetic oscillators, based on Kim & Forger’s interpretation of the circadian clock in mammals. At the core of their models is a negative feedback loop whereby PER proteins (PER1 and PER2) bind to and inhibit their transcriptional activator, BMAL1. For oscillations to occur, the dissociation constant of the PER:BMAL1 complex, K^d, must be ≤ 0.04 nM, which is orders of magnitude smaller than a reasonable expectation of 1–10 nM for this protein complex. We relax this constraint by two modifications to Kim & Forger’s ‘single negative feedback’ (SNF) model: first, by introducing a multistep reaction chain for posttranscriptional modifications of Per mRNA and posttranslational phosphorylations of PER, and second, by replacing the first-order rate law for degradation of PER in the nucleus by a Michaelis-Menten rate law. These modifications increase the maximum allowable K^d to ~2 nM. In a third modification, we consider an alternative rate law for gene transcription to resolve an unrealistically large rate of Per2 transcription at very low concentrations of BMAL1. Additionally, we studied extensions of the SNF model to include a second negative feedback loop (involving REV-ERB) and a supplementary positive feedback loop (involving ROR). Contrary to Kim & Forger’s observations of these extended models, we find that, with our modifications, the supplementary positive feedback loop makes the oscillations more robust than observed in the models with one or two negative feedback loops. However, all three models are similarly robust when accounting for circadian rhythms (~24 h period) with K^d ≥ 1 nM. Our results provide testable predictions for future experimental studies.


Introduction
Most organisms experience perpetual day/night cycles and need to synchronize their physiological functions with this potent external rhythm of light and temperature [1]. Endogenous circadian rhythms meet this demand. These autonomous clock-like rhythms are driven by molecular mechanisms that generate oscillations of~24 h period through negative feedback on gene expression [1][2][3]. Although the genes and proteins constituting the circadian clocks in animals, plants and fungi are quite different, their essential interactions are remarkably similar. In all cases, the clock mechanism features a 'core' negative feedback loop: A activates B activates C inhibits A. In mammals, this loop consists of transcriptional regulation involving six genes: Per1/2, Cry1/2, Bmal1, and Clock [1][2][3][4]. For convenience, in this work we drop the distinction between the homologous pairs of proteins PER1/2 and CRY1/2. In this mechanism (Fig 1), the heterodimeric transcription factor BMAL1:CLOCK activates Per transcription. Per mRNA is then translated in the cytoplasm, where PER protein binds with CRY and enters the nucleus. PER:CRY then binds with BMAL1:CLOCK to block its activation of Per transcription. PER:CRY's cycle of production, nuclear entry, auto-inhibition, and subsequent degradation is widely acknowledged to be the source of circadian rhythmicity [5].
Over the past 50 years, many people have proposed mathematical models of circadian rhythms [5][6][7][8][9][10][11]. In 1965, Brian Goodwin proposed a model of periodic enzyme synthesis based on negative feedback on gene expression [12,13]. At the time, Goodwin was not attending to circadian rhythms, because nothing was known then about the negative feedback of PER on its own synthesis. But his model was picked up later by Peter Ruoff [14][15][16][17][18] to explain many characteristic features of circadian rhythms. Recently, the core negative feedback loop of Goodwin's model was extended with other feedback loops (as in Fig 1) to create more comprehensive and realistic models of circadian rhythms [19][20][21]. While all of these models have much to commend, they suffer from a technical problem with the underlying 'Goodwin oscillator'.
In his model of periodic enzyme synthesis, Goodwin assumed that the end-product of a metabolic pathway functioned as an inhibitor of expression of the gene encoding the first enzyme in the pathway. The inhibition was carried out by p molecules of the end-product binding cooperatively to the transcription factor for the gene. In this scenario the rate of transcription is given by a Hill function, a 1 K p K p þZ p , where Z = concentration of end-product, α 1 = maximum rate of transcription, and K = end-product concentration at half-maximal rate of transcription. In S1 Text, we define Goodwin's model precisely, discuss its basic problem (for the model to oscillate, p must be greater than 8, which is unreasonable), and we describe two changes to Goodwin's model that permit oscillations for smaller values of p.
One particularly interesting modification to Goodwin's model was made by Jae Kyoung Kim and Daniel Forger [20], who replaced Goodwin's view-of negative feedback by cooperative binding of a generic 'repressor' to a gene promoter-with their own model of stoichiometric binding of a repressor (PER:CRY) to an activator (BMAL1:CLOCK) of gene expression. Some characteristic features of the two models have been compared in [22,23]. In the following section, we describe the Kim-Forger model. Then, in the 'Results' section, we show that, like Goodwin's model, the Kim-Forger model also has a 'parameter fragility' problem. This analysis frames our proposals for more robust and realistic mathematical models of circadian clocks.

Kim & Forger's model
In 2012, Kim & Forger [20] presented a model of the negative feedback loop generating autonomous circadian rhythms in mammalian cells (Fig 2A). The Kim-Forger (KF) ODEs are: Kim-Forger SNF Model.
ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi SNF stands for 'single negative feedback' (i.e., the core negative feedback loop involving PER:CRY inhibition of BMAL1:CLOCK). As originally written, the KF model has three A T is the total concentration of BMAL1:CLOCK in the nucleus, and b A free is the concentration of 'free' BMAL1:CLOCK (i.e., not bound to PER:CRY) in the nucleus. (The 'hat' on each variable indicates a concentration in nanomole/liter; and b t is time in hours.) The fac- A T is the probability that BMAL1:CLOCK is not bound to its repressor, PER:CRY. By expressing the rate of transcription of Per mRNA to be proportional to b A free = b A T , Kim & Forger are implicitly assuming that the total number of BMAL1:CLOCK dimers is large enough to saturate the E-boxes on Per genes, and that PER:CRY binds equally well to BMAL1:CLOCK dimers that are either bound or not bound to an E-box. Eq (4) is derived by solving the condition for equilibrium binding of BMAL1:CLOCK (A) and PER:CRY (P) to form an inactive The b a's and b b's are rate constants with appropriate units of concentration and time. It is commonplace in these models to , because this condition is most conducive to oscillations. First of all, we cast the equations on the left side of (1)-(3) into dimensionless form on the right side by defining dimensionless concentrations,   (1) In addition to the SNF model, Kim & Forger proposed two extended models, in which the core negative feedback loop involving PER and BMAL1 is supplemented with either an additional negative feedback from REV-ERB on transcription of the Bmal1 gene (called the NNF model, Fig 2B) or an additional positive feedback from ROR on transcription of the Bmal1 gene (called the PNF model, Fig 2C). Evidences for these interactions are found in references [24][25][26][27][28]. The ODEs of Kim & Forger's 'NNF' and 'PNF' models are presented in S2 Text.
Notice that, in the SNF model, nonlinearity in the transcription term is due to tight stoichiometric binding between PER and BMAL1, not (as in Goodwin's equations) to cooperative participation of nuclear PER in the regulation of Per gene expression. Consequently, the SNF model circumvents the unreasonable cooperativity constraint (p > 8) of Goodwin's model. While the SNF model appears to oscillate robustly and avoid Goodwin's unrealistic constraint (p > 8), the SNF model has an unrealistic constraint of its own. To elaborate, we derive an equation for oscillations to arise in the SNF model.

Locus of Hopf bifurcations in Kim & Forger's SNF model
The condition for a Hopf bifurcation in Eqs (1)-(4) is A T À P ss À 1 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ðA T À P ss À 1Þ 2 þ 4A T where P ss , the steady-state solution of Eqs (1)-(4), satisfies the equation ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi Solving Eqs (5) and (6) simultaneously, we find that and from Eq (6) we derive ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Substituting (8) into (7), we find, after a little algebra, the condition for a Hopf bifurcation: Solving the quadratic Eq (9), we obtain α as a function of A T , as plotted in Fig 3A. We must locate a wild-type (WT) cell somewhere within the oscillatory domain, far enough from the HB locus so that mutant cells overexpressing or under-expressing BMAL1 and PER are still rhythmic. To this end, we propose the following 'five-point criterion' for choosing the values of α and A T for a WT cell: if the point ða WT ; A WT T Þ locates a WT cell on the bifurcation diagram, then the points 1 2 should also lie within the oscillatory domain. We introduce this constraint because: Bmal1 +/− and Clock +/− cells, i.e., a WT ; 1 2 A WT T À � are rhythmic [29,30]; Per1 −/− Per2 +/+ and Per1 +/+ Per2 −/− cells, i.e., We can estimate b K d from the fact that there is a maximum of~30,000 molecules of PER in a mammalian cell [36]; hence, 3 � 10 cytoplasmic PER is transported into the nucleus, so the rate at which PER molecules are lost from must equal the rate at which PER molecules are gained in the nucleus, b a 3 b P c � V N , assuming that there is not significant degradation of PER in the cytoplasm (for an order-of-magnitude estimation, this is a reasonable simplifying assumption). In this case, From the simulation in Fig 3B, we find that P tot = P+P c �2500 at the peak of its oscillation, and from BioNumbers [37], we find that the volume of a typical mam- A T � 40 nM, and the total number of BMAL1 molecules in a nucleus of volume~500 fL would be~12,000. The observed number of BMAL1 molecules in a cell is~25,000 [36], which is not too far off, considering that some fraction of BMAL1 molecules may not localize to the nucleus or act as functional transcription factors.

Is b
K d ¼ 0:04 nM a reasonable estimate of the affinity of PER for BMAL1? We expect the time constant for dissociation of the PER:BMAL1 complex to be on the order of minutes (i.e., b k unbind > 10 À 3 s À 1 ), because, if dissociation of the complex were slower, then the negative feedback of PER on BMAL1 would react sluggishly to changes in nuclear PER concentration. Furthermore, Eq (4) implies equilibrium of PER-BMAL1 binding and would not hold with a much slower dissociation constant. With this estimate of the dissociation rate constant, the binding constant for the complex would have to be However, protein-protein binding rate constants are typically on the order of 10 6 M −1 s −1 [38]. So, we estimate that a physically realistic, minimum value for the dissociation constant of the PER:BMAL1 complex is b K d;min � 1 nM, and we conclude that the dissociation constant used in the SNF model is unrealistically small by at least 25-fold.
Fribourgh et al. [39] recently studied the docking of PER2:CRY1/2 to the core PAS domain of BMAL1:CLOCK and measured b K d � 400 nM. This estimate of b K d is likely too large because the authors used only partial protein sequences. Since the true value of b K d is likely between 1 and 100 nM, we will take b K d ¼ 10 nM as our benchmark. To summarize, we find that circadian oscillations in KF's original SNF model require a value of the PER:BMAL1 dissociation constant, b K d � 0.04 nM (or smaller), that is 250-fold less than a realistic estimate of b K d;est = 10 nM, and 25-fold less than the minimum value, b K d;min = 1 nM. In this work we consider some realistic changes to the SNF model that increase the maximum permissible value of b K d for oscillations. In the process, we come up with some other surprising reassessments of the KF model and its extensions.

Longer feedback loop and saturating PER degradation increase the oscillatory robustness of the Kim-Forger SNF model
Our primary goal in modifying KF models is to alleviate the unreasonable constraint on b K d , the dissociation constant of the PER:BMAL1 complex. To this end, we consider two changes to the SNF model: first, increasing the number of dynamical species in the PER-BMAL1 negative feedback loop, and second, introducing a Michaelis-Menten rate law for the degradation of nuclear PER. These same changes are known to increase the robustness of Goodwin's model (as explained in S1 Text).
Longer feedback loop. In the SNF model, there is only one intermediate (P c ) between Per mRNA (M) and nuclear PER protein (P). However, the primary gene transcript must be processed and exported to the cytoplasm, where it is translated into nascent PER protein. PER protein must be phosphorylated multiple times (PER has 10-20 phosphorylation sites [40,41]) and bound to CRY before it is transported into the nucleus. These steps insert a considerable time lag between Per gene transcription and the negative feedback on BMAL1 activity. To account for this time delay, we replace P c in the SNF model by a sequence of species, P 1 , . . ., P J (note that the first few intermediates could be mRNA species), to obtain the modified ODEs: where N = J+2 is the total length of the negative feedback loop, and A free is still given by Eq (4). This change lengthens the time between Per mRNA transcription and the negative feedback signal generated by nuclear PER and consequently increases the oscillatory potential of the negative feedback loop [42]. The longer feedback loop changes the condition for a Hopf bifurcation to arise in ODEs (11)- (14): the number '8' on the left-hand-side of Eq (5) is replaced by the number Following the same derivation as before, we find that Eq (9) determines α as a function of A T at the Hopf bifurcation, provided that the conservation law for nuclear transport mentioned before, with b a 3 replaced by b a Jþ2 ), we rewrite the relation above as 3 � 10 4 (When there is no chance of misunderstanding, we write P tot for max t P tot ðtÞ.) So, in this case we might estimate that b K d = 100 nM/540 � 0.2 nM. However, 'P tot ' includes Per mRNA species as well as PER protein species. So a better estimate of P tot might be '300', in which case b K d � 0.33 nM, which is still 30-fold smaller than our estimate of b K d;est = 10 nM for the binding of PER to BMAL1. Furthermore, in this case, we estimate b A T = 13 nM (4000 molecules in a nucleus of volume 500 fL), which is perhaps too small compared to the observed number of~25,000 BMAL1 molecules.
Saturating degradation of nuclear PER. PER is degraded by proteasomes after it is polyubiquitinated by the E3 ligase β-TrCP [43]. Because the rate of this enzyme-mediated reaction likely saturates at large substrate concentration, it is reasonable to replace the linear kinetics for degradation of nuclear PER in Eq (3) by a Michaelis-Menten rate law [43], β max and K m are dimensionless parameters; in particular, This change also has the potential to increase the oscillatory robustness of the model. Intuitively, the upper limit to the rate of PER degradation introduced by the Michaelis-Menten rate law causes nuclear PER concentration to react sluggishly to changes in the rate of Per mRNA production, which is another sort of 'lag' in the negative feedback loop.
To keep track of these changes, we introduce the notation SNF (0DN), where D denotes the PER degradation rate law (L for linear or M for Michaelian), and N denotes the number of species in the negative feedback loop. For example, the original KF model is denoted SNF (0L3). The significance of the '0' will become evident shortly.
For the case of 'saturating degradation,' we still scale all concentrations with respect to b K d , but we can no longer derive a closed-form algebraic equation for the locus of Hopf bifurcations. Instead, for N = 8, we searched the four-dimensional parameter space (α, A T , β max , K m ) for oscillations with the smallest value of max t P tot ðtÞ, subject to the constraints that K m > 1 and that the model gives a reasonable domain of oscillations in the (α, A T ) plane (i.e., large enough to satisfy the five-point criterion). We found (see S3 Text) several different combinations of β max and K m that could satisfy these criteria with similar values of max t P tot ðtÞ, suggesting that the model is robust with respect to these criteria. A typical combination is β max = 3.  experimental value, 1 nM < b K d < 10 nM, and the estimate of the total number of BMAL1 molecules per nucleus is acceptable, considering our uncertainty about the localization of BMAL1.
For further information, see S3 Text for notable patterns in the optimization results for SNF (0M8).
A disturbing property of this SNF (0M8) model is that oscillations persist even as A T ! 0, which is clearly impossible because there can be no expression of the Per gene when BMAL1  concentration is zero. The problem, of course, is that the rate law for Per transcription (rate / A free /A T ) is valid only if BMAL1 saturates Per E-boxes, which clearly cannot be true as A T ! 0. To get around this problem, we propose a revised rate law for Per gene transcription.

A more realistic rate law for Per transcription does not affect the robustness of the SNF model
We propose to replace the KF expression for the rate of Per gene transcription (Rate Law 0) by a revised Rate Law 1 that is more realistic for small A T (see S4 Text): For rate law 1, the maximum rate of transcription is a A T K A þA T , which depends on how strongly BMAL1:CLOCK binds to the E-box, as characterized by the dimensionless dissocia- and also, when A T becomes small, the transcription rate is proportional to A free /K A (not A free /A T ). Rate law 0 applies to the case in which binding between PER: CRY and BMAL1:CLOCK is independent of the binding between BMAL1:CLOCK and E-box, and BMAL1:CLOCK complexes saturate Per E-boxes. Rate law 1 relaxes the assumption of saturation of Per E-boxes by BMAL1:CLOCK.
Modified Kim-Forger SNF equations. Taking all of the aforementioned changes into account, we have (see S5 Text): ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where ( ( In the notation SNF(TDN), T denotes the Per transcription rate law (0 or 1), and N = J + 2 = total length of the negative feedback loop.

Models of form SNF (1LN) can be analyzed exactly as SNF(0LN), and the condition for a
Hopf bifurcation is Eq (9), F�α 2 −C(A T )�α+O(A T ) = 0, where, Solving this quadratic equation for α as a function of A T , we plot the locus of Hopf bifurcations for N = 3 and 8 in Fig 6A and 6B. Clearly, this change in rate law makes little difference in the robustness of oscillations for 1�K A �20. As K A increases further, the bifurcation locus moves 'up' and P tot increases, so the estimated value of b K d gets smaller. For SNF (1MN) we have no closed-form algebraic equation for the locus of Hopf bifurcations, so as before, we set N = 8 and searched the five-dimensional parameter space (α, A T , β max , K m , K A ) for oscillations with the smallest value of max t P tot ðtÞ, subject to the constraints a 2 ½10 À 2 ; 10 3 �; A T 2 ½10 À 2 ; 10 2 �; b max 2 ½10 À 2 ; 10 3 �; K m 2 ½1; 10 2 �; K A 2 ½1; 10 2 �; and that the amplitude of oscillation of P tot (t) be larger than 0.5, where ampl ¼ maxÀ min maxþmin . The amplitude constraint is to select for 'robust' oscillations. A summary of these calculations is provided in S3 Text. Briefly, we found~1000 parameter sets with max t P tot ðtÞ = 32.3 ± 3.5 and Period = 20.7 ± 1.1. Then we checked for parameter sets that satisfy the 'five-point criterion'. Results of a typical parameter set (β max = 5, K m = 5.5 and K A = 20) are illustrated in Fig 7A and  7B. In this case, max t P tot ðtÞ = 75, which is larger than '33' because the WT cell ða WT ¼ 50; A WT T ¼ 20Þ must be centered in the oscillatory domain of Fig 7A. Discounting P tot for mRNA species, we estimate max t P tot ðtÞ � 40. Hence, b K d = 2.5 nM and b A T = 50 nM (~15,000 molecules of BMAL1 per nucleus). We conclude that, although rate law 1 is more accurate than rate law 0 for values of A T � K A , it does not improve significantly on our estimates of b K d and b A T .

Adding a negative feedback loop involving REV-ERB does not increase the robustness of circadian oscillations
The dimensional equations on the left-hand-side are cast into dimensionless form with the same definitions used in Eqs (1) To compare NNF to SNF, we must adopt some constraints on the new parameters. First of all, from Narumi et al. [36], we find that the maximum number of REV-ERB molecules during the circadian rhythm in mouse liver cells is 50,000. If all molecules are confined to a nucleus of 500 fL, then 50;000 molec 500 fL For b K V to be greater than, say, 10 nM, we will constrain V MAX so that max t VðtÞ < 10. We continue to insist that the relative amplitude of P tot (t) be > 0.5, and, in addition, we constrain A MAX so that the relative amplitude of A(t) over an oscillation is > 0.2 [36]. We also require that max t A T ðtÞ= max t P tot ðtÞ be as close to 1 as possible (total numbers of BMAL1 and PER proteins are close [36]). Subject to these constraints, we search over the available parameter space a 2 ½10 À 2 ; 10 3 �; A MAX 2 ½10 À 2 ; 10 3 �; b max 2 ½10 À 2 ; 10 3 �; K m 2 ½1; 10 2 �; K A 2 ½1; 10 2 �; V max 2 ½10 À 2 ; 10 2 �; d 2 ½10 À 2 ; 10 2 � to minimize the objective function max t P tot ðtÞ; i.e., to maximize the value of b K d . A summary of these calculations is provided in S3 Text. Briefly, we found~1000 parameter sets with max t P tot ðtÞ � 91.7 ± 24.5 and Period = 22.7 ± 2.1. Again, after checking for parameter sets that satisfy the 'five-point criterion', we plot results of a typical parameter set (β max = 4.5, K m = 2.

An additional positive feedback loop involving ROR increases the robustness of circadian oscillations at a cost
Next, we explore Kim & Forger's PNF model, with similarly modified rate laws for gene transcription: Eq (15-1) for the rates of transcription of Per and Ror genes, and b In the latter rate law, b K R is the dissociation constant for ROR binding to the RORE promoter and ε is the fractional reduction in The dimensional equations on the left-hand-side are cast into dimensionless form with the same definitions used earlier, plus To compare PNF to SNF, we adopt the following constraints. First of all, since the maximum number of ROR molecules during the circadian rhythm in mouse liver cells is 25,000 [36], we estimate that 25;000 molec 500 fL RðtÞ, and consequently we constrain R MAX so that max t RðtÞ < 5. We continue to insist that the relative amplitudes of P tot (t) be > 0.5 and of A(t) be > 0.2, and that max t A T ðtÞ= max t P tot ðtÞ be as close to 1 as possible. Subject to these constraints, we search over the available parameter space a 2 ½10 À 2 ; 10 3 �; A MAX 2 ½10 À 2 ; 10 3 �; b max 2 ½10 À 2 ; 10 3 �; K m 2 ½1; 10 2 �; K A 2 ½1; 10 2 �; R max 2 ½10 À 2 ; 10 2 �; d 2 ½10 À 2 ; 10 2 �; ε 2 ½10 À 4 ; 10 À 1 � to minimize the objective function max t P tot ðtÞ. A summary of these calculations is provided in the S3 Text. Briefly, we found~1000 parameter sets with max t P tot ðtÞ � 3.6 ± 2.2 and Period = 43.9 ± 25.4. Again, after checking for parameter sets that satisfy the 'five-point criterion', we plot results of a typical parameter set (β max = 1.85, K m = 8.5, K A = 35, R MAX = 6.2, ε = 0.0003 and δ = 8) in Fig 9. For the WT cell ða WT ¼ 10; A WT MAX ¼ 15Þ; max t P tot ðtÞ = 12 and avgð b A T Þ = 3. Discounting for mRNA species, max t P tot ðtÞ = 7; hence, b K d = 15 nM and avgð b A T Þ � 45 nM (i.e., 13,000 molecules of BMAL1 per nucleus), which are quite reasonable estimates. Clearly, the PNF (1M8) exhibits more robust oscillations than the SNF(1M8) and NNF(1M8) models and is consistent with b K d � 10 nM, but the oscillation waveforms are unbelievable. For about half of the oscillatory period, R (t), A T (t) and M(t) are � 0, which is inconsistent with observations [36].

The three 1M8 models are about equally robust with respect to circadian oscillations
In Fig 10 we redraw the bifurcation plots for the '1M8' models with colors to indicate oscillatory periods. (For each model, we choose a value of b b 1 , as indicated in the legend, to convert from dimensionless period τ to a period of~24 h given the wild-type parameter values, specified in the legend.) For SNF and NNF models the oscillatory period varies over a range of 22-25 h. The PNF model, though very robust in terms of oscillatory potential, is restricted in exhibiting circadian oscillations (say, 23-26 h) to two regions of expression of BMAL1 (parameter A MAX ) and PER (parameter α); namely, a broad band around the diagonal α + 16�A MAX � 8000, which is clearly seen in Fig 10C, and a triangular region α + A MAX < 35 (for 5 < A MAX < 35, 5 < α < 35) seen in Fig 10D. Our WT simulation (Fig 9) is found in the 'triangular region. ' Oscillations in the 'broad band' (results not shown), although they are consistent with 1 nM < b K d < 10 nM, predict average concentrations of nuclear BMAL1 that are much too small (< 1 nM = 300 molecules per nucleus). Indeed, max t A T ðtÞ is so small despite A MAX being very large, because the concentration of ROR in this region is very small (max t RðtÞ < 0:003) and dA T /dt�A MAX (ε+R)<0.003�A MAX . We note, in passing, that robust oscillations and a broad distribution of periods, as we see in the PNF model, is a common feature of models that combine positive and negative feedback loops [42,44].
We can quantify the 'robustness' of the 1M8 models by measuring the area of 'circadian oscillations' in the (A T , α) or (A MAX , α) plane. To standardize the area, we measure A T (A MAX ) and α as multiples of their WT values, as given in the figure legend. In these units, the areas are NNF:SNF:PNF = 6.6:4:2. These ratios are probably not significantly different in terms of what might possibly be measured experimentally.
In S2 Fig we explore the dependence of oscillatory period on fold-changes in BMAL1 expression (i.e., co-expression of BMAL1 and CLOCK in experiments). SNF (0LN) models are quite insensitive to fold-changes in BMAL1 expression: the change in period is~1 h across the range of oscillations. SNF (0M8), SNF (1M8) and NNF (1M8) models are more sensitive, with a change of 5-8 h across the range. Apparently, the saturating rate law for nuclear PER degradation is responsible for the increased sensitivity to variable expression of BMAL1. PNF (1M8) exhibits very long periods of oscillation for large values of A MAX , as noted, and is comparably sensitive (ΔT � 8 h) over a restricted range of BMAL1 overexpression (up to 2.5 x WT level). When Lee et al. constitutively overexpressed BMAL1:CLOCK dimers (~four-fold) in MEF cells, they observed rhythms of reduced amplitude but normal 24 h period (see their Fig  3C) [34]. This observation is within the limits predicted by models with linear degradation of nuclear PER but is not consistent with the assumption of saturating degradation. S3 Fig shows the trends with respect to PER overexpression (e.g., co-expression of PER2 and CRY1 in experiments). SNF (0L8) is quite insensitive (ΔT � 1 h), whereas SNF (0M8) and SNF(1M8) are more sensitive (ΔT � 4.5 h). NNF (1M8) and PNF (1M8) are even more sensitive to overexpression of PER (8-10 h). These trends are subject to experimental investigation by overexpression of BMAL1:CLOCK and PER2:CRY1.

Discussion
The Kim-Forger (KF) models of mammalian circadian rhythms (called SNF, NNF and PNF) are appealing in many respects, but they rely on an unrealistic requirement for robust oscillations, namely that the equilibrium dissociation constant of the PER:CRY::BMAL1:CLOCK complex must be b K d < 0.04 nM, which is 250-fold smaller than a reasonable value for the dissociation of the PER:CRY::BMAL1:CLOCK complex. This difficulty can be ameliorated by lengthening the core negative feedback loop between Per mRNA transcription and PER:CRY inactivation of BMAL1:CLOCK (the transcription factor driving Per expression), and/or by implementing a Michaelis-Menten rate law for the degradation of nuclear PER. The KF models were further modified by introducing an alternative rate law for BMAL1:CLOCK-mediated transcription of clock genes (Per, Rev-erb and Ror) to correct a problem at low expression of the Bmal1 gene, and by providing more accurate rate laws for the effects of REV-ERB and ROR on Bmal1 expression.
With these modifications, we find (Fig 7) that the SNF (1M8) model can exhibit oscillations for b K d � 2 nM. From biophysical considerations and in vitro measurements, we estimate that b K d;est � 10 nM, so the model constraint is not too far off from expectations. For b K d ¼ 2 nM, the SNF (1M8) model oscillates over a 14-fold range of total BMAL1 concentrations, 10 nM < b A T < 140 nM. For b A T = 30 nM, the corresponding number of BMAL1 molecules in a nucleus of volume 500 fL would be 30 nM � 500 fL ð Þ 6�10 14 1 nmol À � 10 À 15 L 1 fL À � ¼ 9; 000, which is about onethird the observed number (~25,000) of BMAL1 molecules in a mammalian cell [36]. If the remaining BMAL1 molecules are dispersed through the cytoplasm of volume 5000 fL, the cytoplasmic concentration of BMAL1 would be about one-tenth the nuclear concentration, which is not unreasonable for a 'nuclear' protein such as BMAL1. Furthermore, the model focuses on BMAL1:CLOCK complexes that bind E-boxes to regulate gene expression. BMAL1 in this form may account for only a fraction of total BMAL1, if BMAL1, like PER, undergoes multistep post-translational modifications. Indeed, both BMAL1 and CLOCK are known to be phosphorylated at multiple sites, which affects their stability and nuclear accumulation, as well as activity of the BMAL1:CLOCK complex [45][46][47].
Replacing the linear rate law for nuclear PER degradation by a Michaelis-Menten rate law causes a dramatic change in the sensitivity of oscillation to the expression levels of Per2 and Bmal1 (compare Figs 4A and 5A). Models with linear PER degradation, e.g., SNF (0L8), predict that oscillations are possible over an ever broadening range of rates of Per and Bmal1 expression; e.g., 1:8 < aÀ 50 A T À 10 < 75 (approximately) in Fig 4A. For a comparable model with Michaelis-Menten degradation, SNF (0M8), the oscillatory domain is bounded by maximal rates of expression: α < 45 and A T < 45 in Fig 5A. These contrasting results provide a testable prediction for future experimental exploration. By overexpressing Per/Cry genes and/or Bmal1/Clock genes under control of their normal (regulated) promoters (i.e., by manipulating α and A T ), one could test whether nuclear PER degradation operates in a saturated (Michaelis-Menten) or unsaturated (linear) manner, which would be difficult to measure directly in vivo.
In the same experiment, by measuring the dependence of oscillation period on α and A T (i.e., on fold-changes in expression of Per/Cry and Bmal1/Clock), one could investigate a second property of our models (S2 and S3 Figs) that period length is much more sensitive to α and to A T in models with saturated degradation than in models with linear degradation of nuclear PER.
The single negative feedback loop (SNF), whereby PER inhibits its own synthesis, can be supplemented with an auxiliary positive feedback from ROR (PNF) or a second negative feedback from REV-ERB (NNF) on the synthesis of BMAL1. For their versions of these three models (0L3 versions), Kim & Forger observed a 'robustness trend' NNF > SNF > PNF, in terms of the size of the oscillatory domain in parameter space. For our versions of these models, we find that SNF and NNF have similar oscillatory domains, while PNF is much more robust. However, if 'robustness' is defined as the size of the domain of circadian oscillations (22-26 h) in parameter space (fold-changes in expression of PER:CRY and BMAL1:CLOCK complexes), then the SNF(1M8), NNF(1M8) and PNF(1M8) models are nearly equally robust.
Our models could be employed in the future to explore other features of the mammalian circadian clock. For instance, following the lead of Kim and colleagues [48,49], we could address our models to the circadian clock's temperature-compensation and/or phase-shifting properties. Adding these key features may answer some remaining questions about the behaviors of these models. Another question that could be addressed with these models is the function of an anti-sense transcript of the PER2 gene [50]. Furthermore, these models could be applied in chronopharmacology and chronotherapy studies [51]. One such application would be modeling PER2's interaction with the tumor suppressor protein p53 in stressed (e.g., DNA damage) cells compared to unstressed cells [52,53].

Materials and methods
The 'wiring diagrams' (molecular mechanisms) of our models (see Fig 2) were converted into nonlinear ODEs, as described in the main text. The ODEs were solved by standard numerical algorithms, as implemented in MATLAB and XPP-AUTO. Software codes are provided in S1 and S2 Codes for XPP and MATLAB, respectively. Bifurcation diagrams were computed using XPP-Auto, which may be downloaded from www.math.pitt.edu/~bard/xpp/xpp.html. To optimize the parameters of SNF (0M8), SNF (1M8), NNF (1M8) and PNF (1M8) models, we used MATLAB's simulated annealing method ('simulannealbnd') within physiologically reasonable ranges. The parameter ranges and optimization criteria we used for each model and the corresponding cost functions are provided in S3 Text.