Cell Reports

ArticleA Design Principle for a Posttranslational
Biochemical Oscillator
Craig C. Jolley,1,7 Koji L. Ode,3,7 and Hiroki R. Ueda1,2,3,4,5,6,*
1Laboratory for Systems Biology
2Functional Genomics Unit

RIKEN Center for Developmental Biology, 2-2-3 Minatojima-minamimachi, Chuo-ku, Kobe, Hyogo 650-0047, Japan
3Laboratory for Synthetic Biology, RIKEN Quantitative Biology Center, 2-2-3 Minatojima-minamimachi, Chuo-ku, Kobe,
Hyogo 650-0047, Japan
4Department of Biological Sciences, Graduate School of Science, Osaka University, 1-1 Machikaneyama, Toyonaka, Osaka 560-0043, Japan
5Graduate School of Frontier Biosciences, Osaka University, 1-3 Yamadaoka, Suita, Osaka 565-0871 Japan
6Department of Mathematics, Graduate School of Science, Kyoto University, Kitashirakawa Oiwake-cho, Sakyo-ku, Kyoto,
Kyoto 606-8502, Japan
7These authors contributed equally to this work.
*Correspondence: uedah-tky@umin.ac.jp
http://dx.doi.org/10.1016/j.celrep.2012.09.006

SUMMARY

Multisite phosphorylation plays an important role inbiological oscillators such as the circadian clock.
Its general role, however, has been elusive. In thistheoretical study, we show that a simple substratewith two modification sites acted upon by two
opposing enzymes (e.g., a kinase and a phosphatase)can show oscillations in its modification state. An
unbiased computational analysis of this oscillator
reveals two common characteristics: a unidirectional
modification cycle and sequestering of an enzyme bya specific modification state. These two motifs cause
a substrate to act as a coupled system in which a
unidirectional cycle generates single-molecule oscillators, whereas sequestration synchronizes the
population by limiting the available enzyme underconditions in which substrate is in excess. We also
demonstrate the conditions under which the oscillation period is temperature compensated, an important feature of the circadian clock. This theoreticalmodel will provide a framework for analyzing andsynthesizing posttranslational oscillators.
INTRODUCTION
For a wide variety of organisms, biology runs on a schedule. The
control of cellular and organismal processes by biochemicalclocks has a number of potential advantages. The circadianclock system allows organisms to anticipate ecological changes
that correlate with the day/night cycle, regulating both behavioral
patterns and metabolic fluxes. The period of circadian rhythmicity is robust and fairly unaffected by the ambient temperature(Hastings and Sweeney, 1957; Pittendrigh, 1954).
Circadian clocks are thought to be driven by cell-autonomoustranscriptional-translational oscillators (TTOs) (Dibner et al.,938 Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors

2010; Dunlap, 1999; Reppert and Weaver, 2002; Takahashi
et al., 2008; Young and Kay, 2001). In mammalian clocks, thecircadian transcriptional program is mediated by at least
three clock-controlled DNA elements, known as morning-time(E-box/E0 -box or E/E0 -box), day-time (D-box), and night-time(Rev-Erb/ROR binding element or RRE) elements (Ueda et al.,2005). The E/E0 -box-mediated transcriptional program has a critical role in the core autoregulatory loop of the mammalian circadian clock (Gekakis et al., 1998; Sato et al., 2006). In this coreloop, bHLH-PAS transcriptional activators such as BMAL1 and
CLOCK form heterodimers that directly and indirectly activate
the transcription of their target genes, including the transcriptional repressors Per and Cry. PERs and CRYs in turnform repressor complexes that physically associate with the
BMAL1/CLOCK complex to inhibit E/E0 -box-mediated transcription, closing the negative feedback loop (Griffin et al., 1999;
Kume et al., 1999; Okamura et al., 1999). This negative feedback,
together with delayed expression of Cry, is critical for selfsustained circadian clock function (Hogenesch and Ueda,
2011; Ueda et al., 2005; Ukai-Tadenuma et al., 2011).
In the mammalian clock, however, transcription and translation are not the whole story; multisite phosphorylation also playsa central role (Gallego and Virshup, 2007). PER proteins undergorobust circadian changes in phosphorylation (Lee et al., 2001).
Mutations in the potential phosphorylation sites of PERs or in
their kinases, such as casein kinase I ε/d (CKIε/d), alter the circadian period (Lowrey et al., 2000; Meng et al., 2008; Toh et al.,
2001; Xu et al., 2005, 2007). The oscillation period is also
affected by pharmacological perturbation of CKI, CKII, or
GSK3b (Chen et al., 2012; Hirota et al., 2008; Isojima et al.,2009; Meng et al., 2010; Tsuchiya et al., 2009; Walton et al.,2009). PER proteins can be phosphorylated at multiple sites bythese kinases (Camacho et al., 2001; Gallego et al., 2006a; Maieret al., 2009; Schlosser et al., 2005; Takano et al., 2000; Tsuchiya
et al., 2009; Vanselow et al., 2006). Furthermore, temperatureinsensitive phosphorylation of a PER-derived peptide by CKIε/d
has been observed in vitro (Isojima et al., 2009). In addition
to phosphorylation, the reverse reaction, dephosphorylation,

seems to play an important role in the mammalian circadianclock because the oscillation period is also affected by the pharmacological or RNA-interference-mediated perturbation ofprotein phosphatase 1 (PP1), which has been implicated in thedephosphorylation of PER (Gallego et al., 2006b; Lee et al.,2011; Schmutz et al., 2011) and PP5 (Partch et al., 2006).
The central role played by circadian posttranslational oscillators (PTOs) mediated by multisite phosphorylation is best
understood in cyanobacteria. Purified clock proteins from thecyanobacterium Synechococcus elongatus can be used toreconstitute an in vitro posttranslational circadian oscillator(Johnson et al., 2008; Nakajima et al., 2005). This system
contains three essential proteins: KaiA, KaiB, and KaiC. KaiC
exhibits both autophosphorylation and autodephosphorylation
activities; autophosphorylation is enhanced by KaiA, whereas
KaiB inhibits KaiA. When the three proteins are incubated
together with an excess of ATP, the phosphorylation level of
KaiC oscillates with an 24 hr period. The period of KaiC phosphorylation oscillations is fairly unaffected by the incubation
temperature, suggesting that the KaiC PTO is a bona fidecircadian oscillator. Studies have further shown that the phosphorylation cycle involves the sequential phosphorylation anddephosphorylation of two sites in KaiC (Nishiwaki et al., 2007;
Rust et al., 2007).
It is even possible that posttranslational mechanisms might
permit the mammalian circadian clock to operate in the absenceof regulation by gene expression. When cultured fibroblast cells
are treated with reagents that interfere with transcription, one
would expect oscillations based on a TTO to be severely
affected. Instead, circadian oscillations are robust to changes
in transcription rate (Dibner et al., 2009). In addition, cyclic
expression of Bmal1 and Clock seems not to be necessary forcircadian rhythmicity (von Gall et al., 2003). As for CRY, rhythmic,
phase-specific expression of Cry1 is critical for robust circadianoscillation (Ukai-Tadenuma et al., 2011), yet Cry1/:Cry2/cells can exhibit weak circadian oscillations under constant
expression of Cry1 (Ukai-Tadenuma et al., 2011) or a constant
supply of CRYs (Fan et al., 2007). Further evidence comesfrom red blood cells, which show a circadian variation in thecellular redox state despite the absence of transcription (O’Neill
and Reddy, 2011). While these results do not constitute direct
evidence of a posttranslational oscillator that can operate independently from transcription, they suggest that processes otherthan transcription and translation may be responsible for oscillations under some conditions.
One simple and attractive possibility is that, under someconditions, multisite phosphorylation systems can oscillate
without additional regulation. Theoretical studies have indicatedthat multisite phosphorylation systems can exhibit a large
number of stable states under very generic conditions (Markevich et al., 2004; Thomson and Gunawardena, 2009). Previous
studies of oscillatory phosphorylation systems, such as the cyanobacterial clock (Clodong et al., 2007; Rust et al., 2007; van Zon
et al., 2007), mitogen-activated protein kinase (MAPK) systems(Chickarmane et al., 2007; Liu et al., 2011; Qiao et al., 2007;
Shankaran et al., 2009), and other biochemical oscillators (Ferrell
et al., 2011; Kholodenko, 2006; Novák and Tyson, 2008) typically
employed complex regulatory schemes. Cyanobacterial clock

models (Clodong et al., 2007; Rust et al., 2007; van Zon et al.,2007) are constructed based on the autokinase and autophosphatase properties of KaiC. This unique feature makes thedetailed conclusions drawn from these models difficult to apply
to generic enzyme-substrate reactions. MAPK models are constructed based on a more general substrate-enzyme reaction
scheme, but they often involve multiple substrates sharing an
enzyme, protein synthesis, and degradation, or active feedback
regulation of enzyme activity by a substrate. It is not clear
whether such intricate regulation is required for oscillations, orwhether a more generic system would suffice.
Here we show that autonomous, robust oscillations are
possible in a system consisting of a single substrate with two
phosphorylation sites. The substrate is modified by one kinaseand one phosphatase, and no additional regulation is necessary.
In a sense, this is the simplest possible biochemical oscillator. An
extensive search of the model’s parameter space reveals two
design motifs for a simple biochemical oscillator. These features
can be readily associated with known properties of PERs, CKIε/dand other enzyme-substrate pairs described above. Under someconditions, the oscillation period is almost unaffected bysystematic changes in the enzyme reaction rates induced bytemperature differences, suggesting that temperature compensation is possible in this system. The simplicity and generalityof the model make it a potentially useful design tool for denovo biochemical oscillators, and suggest that biological posttranslational oscillation may be more common than previously
suspected.
RESULTS
Multisite Substrate Modifications Can Oscillatewith Two Opposing Enzymes
To elucidate the minimum requirements for a PTO, we began bydefining a set of components. Cyclic changes in substrate phosphorylation state require both a kinase and a phosphatase.
Single-site systems will converge to a unique steady state underrelevant conditions (see Extended Experimental Procedures;
Figure S1), and a two-site system is the simplest one inwhich oscillations might be possible (Figure 1A). The model is
described as a phosphorylation system, but the same results
could be applied to any type of posttranslational modification.
There are four distinct substrate states (Figure 1B): S00, whichis completely unphosphorylated; S01 and S10, which are singly
phosphorylated; and S11, which is doubly phosphorylated. The
same kinase (E) and phosphatase (F) act on both sites. All
reactions involve a 1:1 association of enzymes and substrates,with an intermediate enzyme-substrate complex (Figure 1C).
The reactions are described with the use of mass-action kinetics.
Although the model equations do not explicitly assume enzyme
saturation by substrate, we assumed an excess of substrateover enzyme, and enzyme-substrate binding was usuallysaturated.
We found that sustained oscillations could be observed incertain parameter regions (Figure 1D). After an initial transient,integrations that began under different initial conditions rapidly
converged to a common trajectory with the same amplitude
and period, in other words, a limit-cycle oscillation (Figure 1E).
Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors 939

Figure 1. Multisite Substrate Modifications
Can Oscillate with Two Opposing Enzymes

A
E

F

P

F

S-P

S

E

E

P
S

S

(A) A phosphatase (blue) or a kinase (green) acting
on the substrate alone will yield monotonic convergence to a steady state, as indicated by blueand green curves. A kinase and phosphataseacting in concert are required for sustained oscillations (red).
(B) PTO reaction network. The substrate can bephosphorylated on two distinct sites, and thekinase and phosphatase act on all four phosphorylation states.
(C) Enzymatic reactions involve an intermediate
enzyme-substrate complex.
(D) Time course of oscillations for an example
parameter set.
(E) Oscillations projected into the S00-S10 plane.
Different initial transients (colored arrows) converge rapidly to a single limit-cycle attractor(black).
See also Figure S1.

P
S

S

S

P

P
S

S
time
F

Km

1

K

m

3

P
8
m

K

4

,k
4

P

m

2

P

S11

S10

K

,k

,k 6

6

2

Km

Km

,k

8

S00

7

,k

3

1

,k

K

,k

7

Km

m
5

5

,k

S01

P

Concentration (x 103 M)

D

B

S00 S01 S11 S10
8

4

0

0

12

24
Time (hour)

36

48

Log10(S10) ( M)The ‘‘hit rate’’ of oscillating parameter
sets was rather low, 0.1%. The searchE
was continued until 106 oscillatingC
parameter sets had been found (Fig4
E kE
ure 2A). Because this population may
b1
k3
ES00a 1
+ P
contain several distinct types of oscilla+
kub1
S01
S00
tors, we categorized the parameter2
sets using a clustering algorithm that1accounts for the symmetry of the reaction
kub1 + k1
Km1 =network (see Extended Experimental
kb1
0
Procedures). Two major clusters were
present over a wide range of cluster
12
3
4
diameters, and we chose a diameter for
Log10(S00) ( M)which these two clusters accounted for70% of solutions (Figures 2B andIn our system, a limit-cycle oscillation is present only when the S2A). The parameter distributions for the major clusters are
equilibrium state is unstable.
shown in Figure 2C. Both clusters exhibit approximate
symmetry: fast and slow reactions or strong and weak binding
Oscillatory Systems Possess Two Design Motifs
constants tend to be located at symmetry-related positions.
Once we had represented the system as a system of dynamical For example, in cluster 1, k1, k8, Km2, and Km7 all tend to beequations, our next aim was to determine which values of the high, whereas k2, k7, Km4, and Km5 all tend to be low. The schesystem parameters would lead to sustained oscillations. The matic diagrams shown in Figure 2D will be unchanged if therate constants k1...k8 represent the number of substrate mole- diagram is rotated 180 and the kinase and phosphatase arecules converted to product molecules in a given reaction per interchanged. Symmetry will rarely be exact for individual
enzyme per minute. The binding constants Km1.Km8 are the parameter sets; the symmetric pattern emerges only when asubstrate concentrations at which the rates of reaction 1–8 reach large number of parameter sets are compared.
Based on these histograms, consensus features of each major
half of their maximum. If the binding is strong, then Km will be low
because the enzyme-substrate binding is easily saturated. This cluster can be identified (Figure 2D). Two design motifs are
system is too large and complex for straightforward analytical shared by both major clusters. The first is a unidirectional bias:
study, so we chose to identify solutions numerically through an in general, the rates of the ‘‘forward’’ (clockwise) reactions, k1,
k3, k6, and k8, are higher than the corresponding ‘‘reverse’’ reacextensive search of parameter space.
We generated solutions by selecting parameters from an tions, k5, k7, k2, and k4, respectively. The second is the presence
exponential distribution bounded to the interval 1–1,000 min1 of enzyme-sequestering steps. In cluster 1, the value of Km4 is
for k1–k8 and 0.01–1,000 mM for Km1–Km8. This range corre- usually low, meaning that S10 binds the kinase very strongly.
sponds to reasonable values typically used to model systems The kinase-catalyzed S10/S11 conversion is much slower
such as MAPK (Chickarmane et al., 2007; Liu et al., 2011; Marke- than its reverse reaction; in this sense, the binding of the kinasevich et al., 2004; Qiao et al., 2007; Shankaran et al., 2009).
by S10 is unproductive. Although the S10/S00 conversion is
940 Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors

Figure 2. Oscillatory Systems Possess Two
Design Motifs

B

A
Randomly-generated
parameter sets
(1,280,000,000)"Hit rate" = 0.094%(A) Workflow for motif identification.
(B) Results of QT clustering with a (unitless)diameter of 10.5.
(C) Parameter histograms for the two largest
clusters. Although considerable variability exists,
clear trends in rate and binding constants can be
distinguished.
(D) Schematic representation of cluster motifs.
Fast reactions are indicated by thick arrows; slow
reactions are indicated by dotted arrows. Indirect
negative regulation by enzyme sequestration isindicated by curved arrows with flat ends.
See also Figure S2.

Others (<5%)
31.7%

Parameter sets withlimit-cycle oscillations(1,201,874)

Cluster 1
59.6%
Cluster 2
8.7%

Cluster analysis
by QT algorithmC

Cluster 1

Cluster 2

S01 (i.e., Km5 is low), inhibiting the subsequent S11/S10 step. These sequestration steps create checkpoints that ensurethat the oscillatory cycle cannot proceed
0.004
0.02
until one substrate state has been
depleted and most of the substrate population has moved into the subsequent
00state.
1
10
100
1000
10
100
1000
1
The situation in cluster 2 is only slightly
k (min-1)
k (min-1)
different. A low value of Km1 means that
0.006
0.05
S00 binds the kinase strongly. Unlike theKm1
unproductive S10-kinase complex, thisKm2complex is involved in a forward reaction
Km3
and leads to the production of S01. The
0.025
0.003
Km4
subsequent kinase-catalyzed S01/S11
Km5
reaction, however, cannot proceed until
Km6
S00 has been depleted, resulting in
Km700a checkpoint that requires nearly theKm80.01 0.1
1
10 100 1000
0.01 0.1
1
10 100 1000
entire population to accumulate as S01.
Km ( M)
Km ( M)
The low value of Km8 means that thephosphatase is sequestered by S11, andD
Cluster 1
Cluster 2
therefore most of the population has toaccumulate as S10 before the S10/S00
F P
P
S01
step can proceed.
S01
The remaining small clusters (30% of
the total) are generally similar to either
P
PE
cluster 1 or cluster 2. If all parameter
P FP
S00
S11
S00
S11
sets assigned to neither cluster 1 norcluster 2 are grouped together, theirPparameter distribution is similar to that
P E
S10
S10
of cluster 1 or cluster 2. For example, Figure S2B shows both cluster-2-typecharacteristics (high k3 and k6, low k4
fairly slow, it will proceed to completion because its reverse and k5, and low Km1 and Km8) and cluster-1-type characteristics
reaction (S00/S10) is inhibited by the kinase sequestration. (high k1 and k8, low k2 and k7, and low Km4 and Km5). To a first
The next forward reaction, S00/S01, requires the kinase and is approximation, therefore, solutions can be categorized as being
unable to proceed until S10 has been depleted by conversion similar to either cluster 1 or cluster 2. The presence of two
to S00. Once the S10/S00 reaction has gone to completion, distinct oscillatory clusters in parameter space can also benearly the entire substrate population will be in the S00 state confirmed by principal components analysis (PCA) (Figures
and the kinase will no longer be sequestered by S10. At this point, S2C and S2D).
To summarize, an extensive search of the 16-dimensional
the kinase can bind S00 and the inhibited S00/S01 reaction canproceed. In the same way, the phosphatase is sequestered by parameter space revealed two types of oscillatory parameter
Probability

0.008

k1
k2
k3
k4
k5
k6
k7
k8

Probability

Probability

Probability

0.04

Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors 941

A

B

C
2.5

or

low Km1

P

P
S11

E

low Km8

PF P
S11

S00

% oscillation

P
S01

S00

2.0

low Km5
F P
S01

1.5
1.0
0.5

P
S10

P E
S10
low Km4

or0
Biasedlow Km4
low Km5
low Km1
low Km8

-

-

+
-

+
-

+

+ +
+ +
- +

-

-

-

-

+

+

-

-

-

-

- + +
+ -

-

+ - +

+

-

-

-

-

-

-

-

-

-

+

Figure 3. Functional Significance of Two Design Motifs(A) Unidirectional modification cycle. Clockwise reactions are faster than counterclockwise ones, leading the system to prefer the ordering S00,S01,S11,S10,S00..
(B) Indirect negative regulation by enzyme sequestration.
(C) Imposing the strong binding associated with enzyme sequestration during the random parameter search results in modest gains in the oscillator discovery
rate. If this constraint is combined a clockwise rate constant bias, the discovery rate rises substantially.
See also Figure S3.

sets, illustrated in Figure 2D. Based on their parameter distributions, we have proposed two design motifs: a unidirectional biasof catalytic rate constants and indirect negative regulation
through enzyme sequestration.
Functional Significance of Two Design Motifs
The random-search results surely contain a large number ofstatistical correlations. How then can we determine whether
the design motifs described above actually promote oscillation?
One approach is to begin with a design motif and use it toconstrain the search process (Figures 3A and 3B). If an alleged
motif helps to improve the rate of oscillator discovery, then it isfunctionally significant. Enforcing a clockwise bias in the rate
constants, for example, increases the hit rate 3.4-fold (Figure 3C).
Constraining the search such that Km4 is low improves the
success rate 1.5-fold. If both constraints are combined, a 6.8fold increase is achieved. Requiring both Km4 and Km5 to be
low and maintaining the rate constant bias enforces the entire
cluster 1 pattern; this gives an oscillator discovery rate that is
23.6 times the unbiased value. Qualitatively similar results canbe obtained for cluster 2. Requiring a low value for Km1 fails toimprove the success rate, but combining constraints on Km1
and/or Km8 with the rate constant bias yields significant
improvement.
Parameter histograms similar to those in Figure 2C were also
constructed from the constrained search results (Figure S3).
Although only part of the cluster 1 pattern was enforced (clockwise bias and low Km4), the complete pattern was often present.
For example, constraining Km4 to be low resulted in a low
value for Km5, and the distributions of the catalytic rate constants
are similar to those seen in Figure 2C. Although the distribution
for Km5 is broader than the distribution imposed for Km4,
and the presence of two strong sequestration motifs is not abso942 Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors

lutely essential for oscillations, it is clear that oscillations are far
more probable with both motifs present than with only one
present.
Although the identification of design motifs from unbiased
random search results can be somewhat subjective, the ‘‘inverse
analysis’’ presented in this section shows that the motifs identified are functionally significant. Not only are unidirectional bias
and enzyme sequestration typically present in oscillating parameter sets, demanding their presence also increases the numberof oscillators found.
Period- and Amplitude-Determining Processes
Once we had identified oscillator-enriched regions of the parameter space, we aimed to determine which specific parameters
have a strong impact on system properties, in particular theoscillation period and amplitude. In Figure 4A, an example ofa symmetric cluster 1 (see Table S1) is used as the reference
parameter set for bifurcation calculations (Strogatz, 1994). The
bifurcation diagrams in Figures 4A and S4A show the equilibrium
points (either stable or unstable) of the system, as well as theperiod and amplitude of oscillations when they exist. Similarbifurcation methods can also be used to characterize the robustness of the oscillations (Figure S4B; Extended Experimental
Procedures).
The slow steps by which the enzyme-sequestering substratestates are depleted can be expected to have a large effect on
the period. The bifurcation diagram for k3 confirms this intuition.
Although the impact of k3 on the oscillation amplitude is minimal,
decreasing k3 causes the period to increase dramatically, until
oscillations terminate in a Hopf bifurcation at k3 y 7 3 104 (Figure 4A). Figure 4B illustrates the impact of k3 on period determination. When k3 is reduced by half, the transition phase from S01
to S11 (yellow-shaded) nearly doubles.

A

Figure 4. Period- and Amplitude-Determining Processes

F P
S01

k3
P

S00

S10

E

Km4
10

S00 (x103 M)

8
S00 (x103 M)

P
S11

P

6
4
2
10-4

(A) Two examples of single-parameter bifurcation
diagrams. Solid black lines denote stable fixed
points, and dashed black lines indicate an
unstable fixed point surrounded by a limit cycle.
Red vertical lines denote the parameter value inthe symmetric example set; see Table S1 for
precise values. Colored lines show oscillation
extrema, with the period indicated by the line
color; k3 has a dramatic effect on the period,
whereas Km4 primarily affects the amplitude. A
more extensive set of bifurcation diagrams isshown in Figure S4.
(B) Effects of decreasing a single rate constant.
The oscillation cycle is divided roughly into two
halves: the net conversion of S10 to S01 (via S00) is
shown on a white background, and the net
conversion of S01 to S10 (via S11) is shown ona shaded background. In the symmetric case (left),these phases are of equal length. If k3 (the rate
constant for the slow S01/S11 conversion) ishalved, the length of the S01/S11/S10 phase(shaded background) nearly doubles, indicatingthat the reaction with rate k3 is rate limiting for thisphase. See also Figure S4 and Table S1.

10-2

100
102
k3 (min-1)
10-1

8
6
4
2

104

10-2

101
103
Period (min)

105

100
102
Km4 ( M)104

B
k3= 40.8 min-1

k3= 20.4 min-1

the discrete nature of molecules may be
important and a stochastic approach
becomes necessary. Stochastic simula4000
4000
S10
S10
tions also allow us to quantify the robustS11
S11
ness of oscillations against internal noise.
00
Noise tolerance is an important requireES00a
ES00a
ment for realistic biological models.
0.8
0.8
ES00b
ES00b
The results of stochastic simulations
ES 01
ES 01
for
cluster 1 are shown in Figure 5 (see
0.4
0.4
ES10
Table S2; similar plots for cluster 2 are
ES10
shown in Figures S5A and S5B). For 600
00molecules of each type of enzyme (mainFS 01
FS 01
taining the concentrations used above for
0.8
0.8
FS 10
FS 10
a cell-like volume of 1015 L), the oscillaFS 11b
FS 11btions are extremely robust. As the system
0.4
0.4
shrinks (i.e., the total number of moleFS 11a
FS 11acules decreases), the oscillator runs
0
0
0
12
24
36
48
0
12
24
36
48
slowly relative to the deterministic limit
Time (hour)Time (hour)(Figures 5A, S5C, and S5D; see also
Extended Experimental Procedures).
Km4, the binding constant of kinase to S10, can be decreased Even when only one of each type of enzyme was present, theapparently without limit, whereas the maximum oscillatory value period length was surprisingly robust (the standard deviation ofof Km4 is limited by a supercritical Hopf bifurcation (Figure 4A). the oscillation period was 2% of the mean period). These
Tight binding of the kinase by S10 (and of the phosphatase by results suggest that our model system is fairly robust to theS01) is an important part of the design motif described above, internal noise that will prevail at the cellular scale.
and the loss of sequestration rapidly abolishes oscillations.
Coupling of Single-Molecule Oscillators Is a Design
Principle
Robustness of the Oscillation Period against Stochastic
Fluctuations
Further insight can be gained by tracking the state of individual
The mass-action approach assumes that chemical concentra- molecules. Figure 5B shows the results of stochastic simulations are continuous variables that can take on any positive tions. The top panels show the bulk population behavior, andvalue. At typical cellular volumes and concentrations, however, the middle panels show smoothed trajectories of the four
8000

S00
S01

S00
S01

Number of molecules ( M)

Number of molecules ( M)8000

Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors 943

A
SD / Mean period

Mean period (hour)

48

36

24

0
1

10
100
Number of Enzymes

0.02

0.01

0

1000

1

10
100
1000
Number of Enzymes

S00 S01 S11 S10
1

0

0

48
96
Time (hour)

1.0

S00 S01 S11 S10
5

0

0

48
96
Time (hour)

1

Molecule ID

1

100
0

48
96
Time (hour)

S00 S01 S11 S10
1

0

0

48
96
Time (hour)

Probability

0

48
96
Time (hour)

0

0

48
96
Time (hour)

1

Molecule ID

0

No sequestration

1.0

Probability

Probability

0

Molecule ID

No biasNumber of molecules (x 105)

Cluster 1

Number of molecules (x 104)

Number of molecules (x 105)

B

100
0

48
Time (hour)

96

0

48
96
Time (hour)

1

100
0

48
Time (hour)

96

Figure 5. Coupling of Single-Molecule Oscillators Is a Design Principle(A) Robustness of the oscillation period against stochastic fluctuations. Low-copy-number oscillations are slow (left) and variable (right) relative to larger systems.
In the continuous limit, the period is precisely 24 hr with zero variability.
(B) Top row: Plots of the bulk population for the cluster 1 oscillator (left), similar binding constants lacking a rate constant bias (middle), and similar rate constants
with a clockwise bias but no sequestration (right). Middle row: Smoothed traces showing the state of a single (arbitrarily chosen) substrate molecule. Bottom row:
Results for 100 individually tracked molecules, each shown as a horizontal raster line, color-coded as in the plots above. If synchronization by enzyme
sequestration is removed, the substrates still act as single-molecule oscillators but lack population-level coherence. If the unidirectional rate constant bias iseliminated, the substrates become trapped and fail to progress around the cycle. See also Figure S5 and Table S2.

phosphorylation states of an arbitrarily chosen substrate molecule. In the bottom panels of Figure 5B, each row of pixels represents a single substrate molecule, and different colors represent
the four phosphorylation states.
The middle panel of the ‘‘Cluster 1’’ column (Figure 5B) showscyclic transitions among the four phosphorylation states. Eachindividual substrate molecule acts as a noisy oscillator, synchronized with the rest of the population. In addition, the bottom
944 Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors

panel of the Cluster 1 column indicates clear bands of S11(blue) and S00 (yellow) dominance; individual substrates experience S01 (cyan) transiently during the transition into the S11band, and S10 (red) during the transition into the S00 band. Also
note that the S00/S01 and S11/S10 transitions seem fairly
sharp, whereas the S10/S00 and S01/S11 transitions are
more diffuse. This suggests that the population is coupled at
the S00/S01 and S11/S10 transitions, where they aresynchronized by enzyme sequestration. The same is true forsequestration of phosphatase by S01.
We next observed single-molecule behaviors when either
one of the two design motifs was removed. When the unidirectional motif was removed (by setting k1–k8 equal but maintaining binding affinities), the individual substrate molecules
no longer cycled through phospho-states in a defined order(the ‘‘No bias’’ column). When the sequestration motif wasremoved (by setting all binding affinities equal), phospho-statetransitions failed to synchronize (the ‘‘No sequestration’’column). Although the bulk population converged to a stable
fixed point, individual substrate molecules still transitioned
between phospho-states in a defined order (middle andbottom panels). These results confirm that (1) the design motifof unidirectional modification cycle constitutes a single-molecule oscillator, and (2) the design motif of enzyme sequestration serves to synchronize noisy single-molecule oscillators,
suggesting that coupling of single-molecule oscillators is a
design principle.
The Period of a Simple PTO Can Be Temperature
CompensatedTo qualify as a bona fide circadian clock, an 24 hr cellularrhythm must also exhibit temperature compensation, meaningthat the oscillation period is fairly robust to changes in temperature. The temperature coefficient Q10 measures the factor bywhich a process accelerates when the temperature is raisedby 10 C. Typical Q10 values for biochemical reactions range
from 2 to 3; measured values for circadian clocks range from
0.8–1.4 (Dunlap et al., 2003).
In our model system, one can calculate Q10 values for theoscillation period by assuming temperature coefficients Q10(E)and Q10(F) for the kinase and phosphatase, and rescaling therate constants in accordance with a given temperaturechange (Extended Experimental Procedures; Figures 6A andS6A–S6C). Q10(E) and Q10(F) do not uniquely determine the oscillation Q10 (Q10(cycle)); by analyzing a large number of parameter
sets, one can calculate a histogram of Q10 values (Figure 6B).
In general, larger values of Q10(E) and Q10(F) cause the peak ofQ10(cycle) to be shifted to higher values, and the width of
the Q10(cycle) distribution increases with the difference between
Q10(E) and Q10(F). We considered the case with Q10(E) = 1 andQ10(F) = 3, corresponding to a temperature-insensitive kinase(such as CKIε/d; Isojima et al., 2009) and a typical phosphatase.
The resulting histogram includes populations that are strongly
undercompensated (Q10 >> 1) as well as some that are overcompensated (Q10 < 1; Figure 6C).
Different degrees of temperature sensitivity correlate withchanges in the distribution of individual parameters (Figures
6D, S6D, and S6E). In general, undercompensated oscillators
will have higher rates for the temperature-insensitive kinase
reactions and lower rates for the temperature-sensitive phosphatase reactions, such that the rate-determining steps arestrongly affected by temperature. For well-compensated oscillators, the opposite is true (Figure 6E), and temperature compensation can be observed even in stochastic simulations (Figure 6F;
Table S3). See the Extended Experimental Procedures fora more extensive discussion of this issue.

DISCUSSION
General Applicability of a Design Principle Composedof Two Design Motifs
Although the model developed in this study is very simple, we
can extract a design principle of wider applicability. Two design
motifs appear to be required for robust oscillations (Figure 7A).
The first is a well-defined ordering of phosphorylation states, inwhich the rate constants in the forward direction around theloop are (in general) faster than those in the reverse direction.
The other design motif is the presence of synchronizing checkpoints at which single-molecule oscillators that have progressedmore quickly must stop and wait for the others. This is accomplished by enzyme sequestration: the cycle cannot proceed
past a checkpoint until a strong-binding substrate population
has been sufficiently depleted to permit competing reactions.
The general design principles uncovered in this study can also
be applied to situations in which more than two phosphorylationsites or more than one kinase (and phosphatase) are present(Figure 7B). The principle may be also conserved when protein
degradation and synthesis are incorporated in an appropriate
position: if a substrate acts as a repressor for its own transcription and is degraded depending on phosphorylation, then sucha transcription-degradation pathway will be compatible withthe role of phosphatase reaction (i.e., the disappearance ofphosphorylated substrate will lead to the appearance of unphosphorylated substrate; Figure 7C; Extended Experimental
Procedures).
In this study we considered the behavior of an isolated PTO;
oscillators in real biological systems are responsive to externalsignals that can modify their behavior. In this context, therelationship of our model to ‘‘integrators’’ (which accumulate a
response to incoming signals) and ‘‘resonators’’ (which respondpreferentially to signals with a given frequency) may be an interesting avenue for future study (Conrad et al., 2008; Guantes andPoyatos, 2006).
Two Design Motifs Are Present in the Natural Circadian
ClocksThe design principle mentioned above can easily be related toknown properties of substrates and enzymes involved in circadian clock systems. In the cyanobacterial clock, the phosphorylation states of the clock protein KaiC have a well-defined cyclic
order (Nishiwaki et al., 2007; Rust et al., 2007), similar to the first
design motif (unidirectional modification cycle) in our model. In
higher eukaryotes, CKIε/d predominantly phosphorylates theconsensus sequence pS-x-x-S*, where pS is a phosphorylated
serine, x can be any amino acid, and S* is the target serine(Flotow et al., 1990). This motif promotes ordered phosphorylation: if the substrate contains a pS-x-x-S-x-x-S sequence,
CKIε/d will tend to phosphorylate the second serine before thethird. In mammalian PER proteins, contiguous S-x-x-S motifs
have been shown to be important for period determination(Toh et al., 2001; Vanselow et al., 2006; Xu et al., 2007).
For the second motif (enzyme sequestration), tight binding
between the enzyme and substrate is required. Tight enzymesubstrate binding is found in several enzymes involved in circadian clock systems. In mammals, CKIε/d forms a stable complex
Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors 945

B

Randomly-generated
parameter sets withlimit-cycle oscillations(5,050,984)

and Q

(F)
10

Q10(E) : Q10(F)Q10(E) : Q10(F)

1:1
1:1.5
1:2

1:3
1.5:3

0.1

Probability

Define Q

(E)
10

1.0

0.2
Probability

A

1:2.5
1:3
0
0.1

Parameter setswith successful Q10(cycle)
calculations1

10
Q10(cycle)

D

C

P

P

E

P

0.02

P

1

0
0.001 0.01 0.1

P

10

100
k1 (min-1)

S10

0.01

0
0.01

0
1000

0.1

P

S00

0.06

0.02
Probability

S10

0.03

0

E
Temperature undercompensation(Q10(E)=1, Q10(F)=3)
F P
S01

100

1000

10

100

1000

Km8
P F P
S11

S00

k8

P

Probability

P
S11

1

10

P
S01

S01

100 10001Km1 (mM)

P

1
10
Q10(cycle)

P
S11

0.02

0.03

0.01

P

S00

S11

Probability

Probability

Probability

G1
G2
G3
G4
G5

100

S01

Km1

S10

0.03

10
Q10(cycle)

P

0.06
Q10(E)=1, Q10(F)=3

1

S01

S00

0.04

2:3
2.5:3
3:3

0
0.1

100

k1

0.5

100
k8 (min-1)

1000

P
S10

0.01

0
0.01

0.1

1

10

Km8 (mM)F
S00 S01 S11 S10

P F P
S11

S00

P
S10

+10 OCE±0 OC

Temperature overcompensation(Q10(E)=1, Q10(F)=3)
F P
S01
E
S00

P

P

-10 OC
0

48

S11

P
S10

96
Time (hour)E

Figure 6. The Period of a Simple PTO Can Be Temperature Compensated(A) Workflow diagram.
(B) Histogram of Q10(cycle) values with different combinations of Q10(E) and Q10(F); Q10 for binding/unbinding reactions is assumed to be 3.0 for all calculations.
(C) Histogram of Q10(cycle) values when Q10(E) = 1, Q10(F) = 3. Parameter sets were separated into one of five groups depending on their degree of temperature
overcompensation (blue) or undercompensation (red).

946 Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors

with PERs that is important for clock function (Akashi et al., 2002;
Camacho et al., 2001; Lee et al., 2004; Vielhaber et al., 2000).
PP1 also interacts strongly with PERs (Gallego et al., 2006b). It
is not known whether these tight bindings serve to synchronize
substrates in different phases of an oscillation cycle, but bothexperimental and theoretical studies have highlighted enzyme
sequestration as a source of nonlinearity in biochemical reactions (Blüthgen et al., 2006; Legewie et al., 2006, 2007; Markevich et al., 2004).
Comparison with Other PTO Models
The design motifs uncovered in this study can also be seen inother mathematical models of PTOs. MAPK is typically modeled
as a two-site system with strictly sequential phosphorylation,such that only three states exist (Huang and Ferrell, 1996; Kholodenko, 2000; Markevich et al., 2004). Feed-forward inhibitioninteractions similar to those described in this study can generatebistability, but oscillations require larger cascades or feedback
loops. The cyanobacterial circadian clock has also beenmodeled extensively (Clodong et al., 2007; Hatakeyama andKaneko, 2012; Rust et al., 2007; van Zon et al., 2007). Thesemodels also tend to converge on design features similar toours, in that they typically employ some mechanism (often
involving conformational changes or complex formation) toensure that individual KaiC hexamers visit phosphorylation
states in a well-defined order, and the oscillations of the individual hexamers are synchronized through competition fora limited pool of KaiA. A more detailed examination of individualmodels can be found in the Extended Experimental Procedures.
Although previous studies generally focused on correctly
reproducing the behavior of a specific experimentally characterized system, here we focused on mapping out possible behaviors for a very general system. It is therefore quite striking that
our approach, which began with no constraining assumptions
about what kind of design features should promote oscillation,
should lead to fairly similar conclusions.
Strategies for Experimental Validation
Our modeling results can inform the design or discovery of PTOs
in several ways. No unusual substrate properties are necessary,
so a simple substrate, such as a short peptide, might beadequate.
The low probability of finding oscillations in our nonbiased
search may raise questions about the feasibility of discovering
PTOs experimentally. Results from the constrained-parametersearch, however, show that one can increase the probability>20-fold by requiring sequential (de)phosphorylation and tight
enzyme-substrate binding. These features are likely to be implemented (at least in part) in substrates and enzymes found innatural oscillatory systems. Accordingly, as a starting point forexperimental investigation, CKIε/d, a PER-derived peptide, andan appropriate phosphatase might be a suitable ex vivo PTO.
CONCLUSIONSThis study illustrates the potential for robust oscillations to exist
in a very simple kinase-phosphatase reaction system. Thisrequires two design motifs: a well-ordered sequence of phosphorylation states and synchronization by enzyme sequestration. Although the vast majority of possible parameter sets examined here did not result in sustained oscillations, we were able toconstruct a typical example of an oscillatory parameter set withreaction rates and binding constants that were within a reasonable range for polypeptide kinases. The ultimate test of the utilityof this proposed model, ‘‘proof by synthesis,’’ lies ahead.
EXPERIMENTAL PROCEDURES
Modeling the PTO with Two Phosphorylation Sites
The PTO can be described using a system of 14 coupled ordinary differential
equations (ODEs) describing the concentrations of four substrate phosphoforms, eight substrate-enzyme complexes, and two unbound enzymes(Figure 1B; the full system of equations is presented in Extended Experimental
Procedures). In the case of branching reactions (e.g., phosphorylation of S00 toform either S01 or S10), it is assumed that two different types of enzymesubstrate complexes are formed (ES00a or ES00b, respectively), consistent
with the presence of two distinct phosphorylation sites on the substratemolecule. Calculations were also performed with the opposite assumption(a single complex ES00 that can lead to either S01 or S10), and the results
were qualitatively similar (Figure S1).
Biased Parameter Search
When a forward bias in the rate constants was desired (Figure 3), numbers
were selected from the exponential distribution two at a time, with the larger
number being assigned to a forward (clockwise) reaction (k1, k3, k6, or k8)
and the smaller number assigned to the corresponding reverse (counterclockwise) reaction (k5, k7, k2, or k4, respectively; Figure 3A). Enzyme sequestration
was encouraged by selecting constrained binding constants from an exponential distribution over 0.01–0.05 mM rather than 0.01–1,000 mM (Figure 3B).
Defining ‘‘Typical’’ Parameter Sets
A stereotypical example was generated for each of the dominant clusters; theparameter values chosen are listed in Table S1, and oscillations for the cluster
1 example are depicted in Figure 4A. These example clusters were designed to
exhibit symmetry, so that ki = kj and Kmi = Kmj for (i,j) ˛ {(1,8),(2,7),(3,6),(4,5)}.
Rate constants were scaled to ensure an oscillation period of 24 hr.
Bifurcation Analysis
Beginning from these stereotypical parameter sets, each parameter was
systematically varied beyond the range of values over which oscillations are
possible, and concentrations at the fixed points and the oscillation extrema
were calculated. Although the starting parameter sets were symmetric, thissymmetry was not maintained during bifurcation calculations. Note that,
because of symmetry, bifurcation calculations were needed for only half ofthe parameters. For example, the diagram for k8 can be obtained by(D) Examples of parameter histograms for different oscillator groups, color-coded as in (C). The parameters k1 and k8 are located at symmetric positions in the
network (Figure 1B) and would be equivalent if the roles of kinase and phosphatase were interchanged.
(E) Schematic diagrams of the extremes of temperature sensitivity. In the undercompensated case (top) the rate-limiting steps are catalyzed by the temperaturesensitive phosphatase, whereas in the overcompensated case (bottom) the temperature-insensitive kinase steps are rate limiting.
(F) Stochastic simulation results for a particular temperature-compensated parameter set.
See also Figure S6 and Table S3.

Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors 947

A

monitor the behavior of individual substrate molecules, we modified thestochastic simulations so that in addition to a large pool of ‘‘ordinary’’ substratemolecules, the system contains a smaller number of ‘‘labeled’’ substratemolecules whose progression through the four phosphorylation states canbe tracked. The interactions of the four phospho-forms of these labeled
substrate molecules with the kinase and phosphatase are identical to their
unlabeled counterparts, but each can be followed in the stochastic simulation
as a separate molecular population, allowing us to observe how a single
substrate molecule experiences the oscillatory cycle.

F P
S01

P

S00

P
S11

P
S10

E
SUPPLEMENTAL INFORMATION
Supplemental Information includes Extended Experimental Procedures, six
figures, and three tables and can be found with this article online at http://dx.doi.org/10.1016/j.celrep.2012.09.006.

B

LICENSING INFORMATION

E’

P
S001
E

PP
S011

P PP
S111
1

S000

This is an open-access article distributed under the terms of the Creative
Commons Attribution-Noncommercial-No Derivative Works 3.0 Unported
License (CC-BY-NC-ND; http://creativecommons.org/licenses/by-nc-nd/3.0/
legalcode).
ACKNOWLEDGMENTSP
S100

We thank Drs. Tatsuo Shibata and Ryokichi Tanaka for valuable comments.
C.C.J. was supported by the RIKEN Foreign Postdoctoral Researcher
Program. Large-scale calculations used the RIKEN Integrated Cluster ofClusters. This research was supported by an intramural Grant-in-Aid fromthe RIKEN Center for Developmental Biology and the Quantitative Biology
Center (H.R.U.), the Uehara Memorial Foundation (H.R.U.), the Mitsubishi
Foundation (H.R.U.), the President’s Fund from RIKEN (H.R.U.), and a Grantin-Aid for Scientific Research on Innovative Areas (No. 3306) from the Ministryof Education, Culture, Sports, Science, and Technology, Japan (H.R.U.).

PP
S110

C
E’

P
S001
E

PP
S011

P PP
S111
1

S000

Received: June 11, 2012
Revised: August 22, 2012
Accepted: September 8, 2012
Published online: October 18, 2012
REFERENCES
Akashi, M., Tsuchiya, Y., Yoshino, T., and Nishida, E. (2002). Control of intracellular dynamics of mammalian period proteins by casein kinase I epsilon(CKIepsilon) and CKIdelta in cultured cells. Mol. Cell. Biol. 22, 1693–1703.
Blüthgen, N., Bruggeman, F.J., Legewie, S., Herzel, H., Westerhoff, H.V., andKholodenko, B.N. (2006). Effects of sequestration on signal transduction
cascades. FEBS J. 273, 895–906.

mRNA
gene S
Figure 7. Generalizations of the PTO Model(A) Simplified view of the PTO model examined in this work.
(B) The mechanisms outlined in this work could be realized in a more complex
system with a larger number of phosphorylation states and more than two
enzymes.
(C) Similar considerations could apply in a case in which the dephosphorylation of highly phosphorylated substrate is replaced by coupled processes ofdegradation and gene expression.

Camacho, F., Cilio, M., Guo, Y., Virshup, D.M., Patel, K., Khorkova, O., Styren,
S., Morse, B., Yao, Z., and Keesler, G.A. (2001). Human casein kinase Idelta
phosphorylation of human circadian clock proteins period 1 and 2. FEBS
Lett. 489, 159–165.
Chen, Z., Yoo, S.H., Park, Y.S., Kim, K.H., Wei, S., Buhr, E., Ye, Z.Y., Pan, H.L.,and Takahashi, J.S. (2012). Identification of diverse modulators of central andperipheral circadian clocks by high-throughput chemical screening. Proc.
Natl. Acad. Sci. USA 109, 101–106.
Chickarmane, V., Kholodenko, B.N., and Sauro, H.M. (2007). Oscillatory
dynamics arising from competitive inhibition and multisite phosphorylation.
J. Theor. Biol. 244, 68–76.

exchanging S00 with S11 and S01 with S10 in the k1 calculation; the bifurcation
points and period dependence will be identical.

Clodong, S., Dühring, U., Kronk, L., Wilde, A., Axmann, I., Herzel, H., andKollmann, M. (2007). Functioning and robustness of a bacterial circadianclock. Mol. Syst. Biol. 3, 90.

Stochastic Simulation
We studied the reaction networks using the stochastic simulation algorithm(Gillespie, 1977) as implemented in SPPARKS (Plimpton et al., 2009). To

Conrad, E., Mayo, A.E., Ninfa, A.J., and Forger, D.B. (2008). Rate constants
rather than biochemical mechanism determine behaviour of genetic clocks.
J. R. Soc. Interface 5(Suppl 1), S9–S15.

948 Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors

Dibner, C., Sage, D., Unser, M., Bauer, C., d’Eysmond, T., Naef, F., and Schibler, U. (2009). Circadian gene expression is resilient to large fluctuations inoverall transcription rates. EMBO J. 28, 123–134.

Lee, C., Etchegaray, J.P., Cagampang, F.R., Loudon, A.S., and Reppert, S.M.
(2001). Posttranslational mechanisms regulate the mammalian circadianclock. Cell 107, 855–867.

Dibner, C., Schibler, U., and Albrecht, U. (2010). The mammalian circadiantiming system: organization and coordination of central and peripheral clocks.
Annu. Rev. Physiol. 72, 517–549.

Lee, C., Weaver, D.R., and Reppert, S.M. (2004). Direct association between
mouse PERIOD and CKIepsilon is critical for a functioning circadian clock.
Mol. Cell. Biol. 24, 584–594.

Dunlap, J.C. (1999). Molecular bases for circadian clocks. Cell 96, 271–290.

Lee, H.M., Chen, R., Kim, H., Etchegaray, J.P., Weaver, D.R., and Lee, C.
(2011). The period of the circadian oscillator is primarily determined by thebalance between casein kinase 1 and protein phosphatase 1. Proc. Natl.
Acad. Sci. USA 108, 16451–16456.

Dunlap, J.C., Loros, J.J., and Decoursey, P.J. (2003). Chronobiology: Biological Timekeeping, First Edition (Sunderland, MA: Sinauer Associates).
Fan, Y., Hida, A., Anderson, D.A., Izumo, M., and Johnson, C.H. (2007). Cyclingof CRYPTOCHROME proteins is not necessary for circadian-clock function in
mammalian fibroblasts. Curr. Biol. 17, 1091–1100.
Ferrell, J.E., Jr., Tsai, T.Y., and Yang, Q. (2011). Modeling the cell cycle: why docertain circuits oscillate? Cell 144, 874–885.
Flotow, H., Graves, P.R., Wang, A.Q., Fiol, C.J., Roeske, R.W., and Roach, P.J.
(1990). Phosphate groups as substrate determinants for casein kinase I action.
J. Biol. Chem. 265, 14264–14269.
Gallego, M., and Virshup, D.M. (2007). Post-translational modifications regulate the ticking of the circadian clock. Nat. Rev. Mol. Cell Biol. 8, 139–148.
Gallego, M., Eide, E.J., Woolf, M.F., Virshup, D.M., and Forger, D.B. (2006a).
An opposite role for tau in circadian rhythms revealed by mathematicalmodeling. Proc. Natl. Acad. Sci. USA 103, 10618–10623.
Gallego, M., Kang, H., and Virshup, D.M. (2006b). Protein phosphatase 1 regulates the stability of the circadian protein PER2. Biochem. J. 399, 169–175.
Gekakis, N., Staknis, D., Nguyen, H.B., Davis, F.C., Wilsbacher, L.D., King,
D.P., Takahashi, J.S., and Weitz, C.J. (1998). Role of the CLOCK protein in
the mammalian circadian mechanism. Science 280, 1564–1569.
Gillespie, D.T. (1977). Exact stochastic simulation of coupled chemical
reactions. J. Phys. Chem. 81, 2340–2361.
Griffin, E.A., Jr., Staknis, D., and Weitz, C.J. (1999). Light-independent role of
CRY1 and CRY2 in the mammalian circadian clock. Science 286, 768–771.
Guantes, R., and Poyatos, J.F. (2006). Dynamical principles of two-component
genetic oscillators. PLoS Comput. Biol. 2, e30.
Hastings, J.W., and Sweeney, B.M. (1957). On the mechanism of temperature
independence in a biological clock. Proc. Natl. Acad. Sci. USA 43, 804–811.

Legewie, S., Blüthgen, N., and Herzel, H. (2006). Mathematical modeling
identifies inhibitors of apoptosis as mediators of positive feedback and bistability. PLoS Comput. Biol. 2, e120.
Legewie, S., Schoeberl, B., Blüthgen, N., and Herzel, H. (2007). Competingdocking interactions can bring about bistability in the MAPK cascade.
Biophys. J. 93, 2279–2288.
Liu, P., Kevrekidis, I.G., and Shvartsman, S.Y. (2011). Substrate-dependent
control of ERK phosphorylation can lead to oscillations. Biophys. J. 101,
2572–2581.
Lowrey, P.L., Shimomura, K., Antoch, M.P., Yamazaki, S., Zemenides, P.D.,
Ralph, M.R., Menaker, M., and Takahashi, J.S. (2000). Positional syntenic
cloning and functional characterization of the mammalian circadian mutation
tau. Science 288, 483–492.
Maier, B., Wendt, S., Vanselow, J.T., Wallach, T., Reischl, S., Oehmke, S.,
Schlosser, A., and Kramer, A. (2009). A large-scale functional RNAi screen
reveals a role for CK2 in the mammalian circadian clock. Genes Dev. 23,
708–718.
Markevich, N.I., Hoek, J.B., and Kholodenko, B.N. (2004). Signaling switches
and bistability arising from multisite phosphorylation in protein kinasecascades. J. Cell Biol. 164, 353–359.
Meng, Q.J., Logunova, L., Maywood, E.S., Gallego, M., Lebiecki, J., Brown,
T.M., Sládek, M., Semikhodskii, A.S., Glossop, N.R., Piggins, H.D., et al.
(2008). Setting clock speed in mammals: the CK1 epsilon tau mutation inmice accelerates circadian pacemakers by selectively destabilizing PERIOD
proteins. Neuron 58, 78–88.

Hatakeyama, T.S., and Kaneko, K. (2012). Generic temperature compensationof biological clocks by autonomous regulation of catalyst concentration. Proc.
Natl. Acad. Sci. USA 109, 8109–8114.

Meng, Q.J., Maywood, E.S., Bechtold, D.A., Lu, W.Q., Li, J., Gibbs, J.E.,
Dupré, S.M., Chesham, J.E., Rajamohan, F., Knafels, J., et al. (2010). Entrainment of disrupted circadian behavior through inhibition of casein kinase 1(CK1) enzymes. Proc. Natl. Acad. Sci. USA 107, 15240–15245.

Hirota, T., Lewis, W.G., Liu, A.C., Lee, J.W., Schultz, P.G., and Kay, S.A. (2008).
A chemical biology approach reveals period shortening of the mammalian
circadian clock by specific inhibition of GSK-3beta. Proc. Natl. Acad. Sci.
USA 105, 20746–20751.

Nakajima, M., Imai, K., Ito, H., Nishiwaki, T., Murayama, Y., Iwasaki, H.,
Oyama, T., and Kondo, T. (2005). Reconstitution of circadian oscillation of
cyanobacterial KaiC phosphorylation in vitro. Science 308, 414–415.

Hogenesch, J.B., and Ueda, H.R. (2011). Understanding systems-level properties: timely stories from the study of clocks. Nat. Rev. Genet. 12, 407–416.

Nishiwaki, T., Satomi, Y., Kitayama, Y., Terauchi, K., Kiyohara, R., Takao, T.,
and Kondo, T. (2007). A sequential program of dual phosphorylation of KaiC
as a basis for circadian rhythm in cyanobacteria. EMBO J. 26, 4029–4037.

Huang, C.Y., and Ferrell, J.E., Jr. (1996). Ultrasensitivity in the mitogen-activated protein kinase cascade. Proc. Natl. Acad. Sci. USA 93, 10078–10083.
Isojima, Y., Nakajima, M., Ukai, H., Fujishima, H., Yamada, R.G., Masumoto,K.H., Kiuchi, R., Ishida, M., Ukai-Tadenuma, M., Minami, Y., et al. (2009).
CKIepsilon/d-dependent phosphorylation is a temperature-insensitive,
period-determining process in the mammalian circadian clock. Proc. Natl.
Acad. Sci. USA 106, 15744–15749.
Johnson, C.H., Mori, T., and Xu, Y. (2008). A cyanobacterial circadian clockwork. Curr. Biol. 18, R816–R825.
Kholodenko, B.N. (2000). Negative feedback and ultrasensitivity can bringabout oscillations in the mitogen-activated protein kinase cascades. Eur. J.
Biochem. 267, 1583–1588.
Kholodenko, B.N. (2006). Cell-signalling dynamics in time and space. Nat. Rev.
Mol. Cell Biol. 7, 165–176.
Kume, K., Zylka, M.J., Sriram, S., Shearman, L.P., Weaver, D.R., Jin, X., Maywood, E.S., Hastings, M.H., and Reppert, S.M. (1999). mCRY1 and mCRY2 areessential components of the negative limb of the circadian clock feedback
loop. Cell 98, 193–205.

Novák, B., and Tyson, J.J. (2008). Design principles of biochemical oscillators.
Nat. Rev. Mol. Cell Biol. 9, 981–991.
O’Neill, J.S., and Reddy, A.B. (2011). Circadian clocks in human red blood
cells. Nature 469, 498–503.
Okamura, H., Miyake, S., Sumi, Y., Yamaguchi, S., Yasui, A., Muijtjens, M.,
Hoeijmakers, J.H., and van der Horst, G.T. (1999). Photic induction of mPer1
and mPer2 in cry-deficient mice lacking a biological clock. Science 286,
2531–2534.
Partch, C.L., Shields, K.F., Thompson, C.L., Selby, C.P., and Sancar, A. (2006).
Posttranslational regulation of the mammalian circadian clock by cryptochrome and protein phosphatase 5. Proc. Natl. Acad. Sci. USA 103, 10467–
10472.
Pittendrigh, C.S. (1954). On temperature independence in the clock system
controlling emergence time in Drosophila. Proc. Natl. Acad. Sci. USA 40,
1018–1029.
Plimpton, S., Battaile, C., Chandross, M., Holm, L., Thompson, A., Tikare, V.,
Wagner, G., Webb, E., Zhou, X., Garcia-Cardona, C., et al. (2009). Crossing the

Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors 949

mesoscale no-man’s land via parallel kinetic Monte Carlo. Sandia Report
SAND2009-6226.
Qiao, L., Nachbar, R.B., Kevrekidis, I.G., and Shvartsman, S.Y. (2007).
Bistability and oscillations in the Huang-Ferrell model of MAPK signaling.
PLoS Comput. Biol. 3, 1819–1826.
Reppert, S.M., and Weaver, D.R. (2002). Coordination of circadian timing in
mammals. Nature 418, 935–941.
Rust, M.J., Markson, J.S., Lane, W.S., Fisher, D.S., and O’Shea, E.K. (2007).
Ordered phosphorylation governs oscillation of a three-protein circadianclock. Science 318, 809–812.
Sato, T.K., Yamada, R.G., Ukai, H., Baggs, J.E., Miraglia, L.J., Kobayashi, T.J.,
Welsh, D.K., Kay, S.A., Ueda, H.R., and Hogenesch, J.B. (2006). Feedback
repression is required for mammalian circadian clock function. Nat. Genet.
38, 312–319.

Tsuchiya, Y., Akashi, M., Matsuda, M., Goto, K., Miyata, Y., Node, K., andNishida, E. (2009). Involvement of the protein kinase CK2 in the regulation of
mammalian circadian rhythms. Sci. Signal. 2, ra26.
Ueda, H.R., Hayashi, S., Chen, W., Sano, M., Machida, M., Shigeyoshi, Y., Iino,
M., and Hashimoto, S. (2005). System-level identification of transcriptional
circuits underlying mammalian circadian clocks. Nat. Genet. 37, 187–192.
Ukai-Tadenuma, M., Yamada, R.G., Xu, H., Ripperger, J.A., Liu, A.C., andUeda, H.R. (2011). Delay in feedback repression by cryptochrome 1 is requiredfor circadian clock function. Cell 144, 268–281.
van Zon, J.S., Lubensky, D.K., Altena, P.R., and ten Wolde, P.R. (2007). An
allosteric model of circadian KaiC phosphorylation. Proc. Natl. Acad. Sci.
USA 104, 7420–7425.

Schlosser, A., Vanselow, J.T., and Kramer, A. (2005). Mapping of phosphorylation sites by a multi-protease approach with specific phosphopeptide
enrichment and NanoLC-MS/MS analysis. Anal. Chem. 77, 5243–5250.

Vanselow, K., Vanselow, J.T., Westermark, P.O., Reischl, S., Maier, B., Korte,
T., Herrmann, A., Herzel, H., Schlosser, A., and Kramer, A. (2006). Differential
effects of PER2 phosphorylation: molecular basis for the human familial
advanced sleep phase syndrome (FASPS). Genes Dev. 20, 2660–2672.

Schmutz, I., Wendt, S., Schnell, A., Kramer, A., Mansuy, I.M., and Albrecht, U.
(2011). Protein phosphatase 1 (PP1) is a post-translational regulator of the
mammalian circadian clock. PLoS ONE 6, e21325.

Vielhaber, E., Eide, E., Rivers, A., Gao, Z.H., and Virshup, D.M. (2000). Nuclear
entry of the circadian regulator mPER1 is controlled by mammalian casein
kinase I epsilon. Mol. Cell. Biol. 20, 4888–4899.

Shankaran, H., Ippolito, D.L., Chrisler, W.B., Resat, H., Bollinger, N., Opresko,
L.K., and Wiley, H.S. (2009). Rapid and sustained nuclear-cytoplasmic ERK
oscillations induced by epidermal growth factor. Mol. Syst. Biol. 5, 332.

von Gall, C., Noton, E., Lee, C., and Weaver, D.R. (2003). Light does
not degrade the constitutively expressed BMAL1 protein in the mouse suprachiasmatic nucleus. Eur. J. Neurosci. 18, 125–133.

Strogatz, S.H. (1994). Nonlinear Dynamics and Chaos: With Applications
to Physics, Biology, Chemistry, and Engineering (Reading, MA: AddisonWesley).
Takahashi, J.S., Hong, H.K., Ko, C.H., and McDearmon, E.L. (2008). The
genetics of mammalian circadian order and disorder: implications for physiology and disease. Nat. Rev. Genet. 9, 764–775.
Takano, A., Shimizu, K., Kani, S., Buijs, R.M., Okada, M., and Nagai, K. (2000).
Cloning and characterization of rat casein kinase 1epsilon. FEBS Lett. 477,
106–112.
Thomson, M., and Gunawardena, J. (2009). Unlimited multistability in multisite
phosphorylation systems. Nature 460, 274–277.
Toh, K.L., Jones, C.R., He, Y., Eide, E.J., Hinz, W.A., Virshup, D.M., Ptácek,
L.J., and Fu, Y.H. (2001). An hPer2 phosphorylation site mutation in familial
advanced sleep phase syndrome. Science 291, 1040–1043.

950 Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors

Walton, K.M., Fisher, K., Rubitski, D., Marconi, M., Meng, Q.J., Sládek, M.,
Adams, J., Bass, M., Chandrasekaran, R., Butler, T., et al. (2009). Selective
inhibition of casein kinase 1 epsilon minimally alters circadian clock period.
J. Pharmacol. Exp. Ther. 330, 430–439.
Xu, Y., Padiath, Q.S., Shapiro, R.E., Jones, C.R., Wu, S.C., Saigoh, N., Saigoh,
K., Ptácek, L.J., and Fu, Y.H. (2005). Functional consequences of a CKIdeltamutation causing familial advanced sleep phase syndrome. Nature 434,
640–644.
Xu, Y., Toh, K.L., Jones, C.R., Shin, J.Y., Fu, Y.H., and Ptácek, L.J. (2007).
Modeling of a human circadian mutation yields insights into clock regulation
by PER2. Cell 128, 59–70.
Young, M.W., and Kay, S.A. (2001). Time zones: a comparative genetics of
circadian clocks. Nat. Rev. Genet. 2, 702–715.

Supplemental Information
EXTENDED EXPERIMENTAL PROCEDURES
Full System of Equations (Mass-Action Kinetics)(See ‘‘Modeling the PTO with Two Phosphorylation Sites’’ in Experimental Procedures section.)d½S00 =  ðkb1 + kb2 Þ½S00 ½E + kub1 ½ES00a  + kub2 ½ES00b  + k5 ½FS01  + k6 ½FS10 dt
d½S01 =  kb3 ½S01 ½E + kub3 ½ES01   kb5 ½S01 ½F + kub4 ½FS01  + k1 ½ES00a  + k7 ½FS11a dt
d½S10 =  kb4 ½S10 ½E + kub4 ½ES10   kb6 ½S10 ½F + kub6 ½FS10  + k2 ½ES00b  + k8 ½FS11b 
dt
d½S11 =  ðkb7 + kb8 Þ½S11 ½F + kub7 ½FS11a  + kub8 ½FS11b  + k3 ½ES01  + k4 ½ES10 
dt
d½ES00a = kb1 ½E½S00   ðkub1 + k1 Þ½ES00a 
dt
d½ES00b = kb2 ½E½S00   ðkub2 + k2 Þ½ES00b 
dt
d½ES01 = kb3 ½E½S01   ðkub3 + k3 Þ½ES01 
dt
d½ES10 = kb4 ½E½S10   ðkub4 + k4 Þ½ES10 
dt
d½FS01 = kb5 ½F½S01   ðkub5 + k5 Þ½FS01 
dt
d½FS10 = kb6 ½F½S10   ðkub6 + k6 Þ½FS10 
dt
d½FS11a = kb7 ½F½S11   ðkub7 + k7 Þ½FS11a dt
d½FS11b = kb8 ½F½S11   ðkub8 + k8 Þ½FS11b dt
d½E=  ððkb1 + kb2 Þ½S00  + kb3 ½S01  + kb4 ½S10 Þ½E + ðkub1 + k1 Þ½ES00a  + ðkub2 + k2 Þ½ES00b  + ðkub3 + k3 Þ½ES01  + ðkub4 + k4 Þ½ES10 dt
d½F=  ðkb5 ½S01  + kb6 ½S10  + ðkb7 + kb8 Þ½S11 Þ½F + ðkub5 + k5 Þ½FS01  + ðkub6 + k6 Þ½FS10  + ðkub7 + k7 Þ½FS11a  + ðkub8 + k8 Þ½FS11b 
dt[S00], [S01], [S10] and [S11] represent concentrations of the four substrate phosphorylation states, [E] and [F] are the concentrationsof free (unbound) kinase and phosphatase, and [ES00a], [ES00b], [ES01], [ES10], [FS01], [FS10], [FS11a] and [FS11b] are the concentrationsof the eight distinct enzyme-substrate complexes. The catalytic rate constants (i.e., number of enzyme-substrate complexes converted to enzyme/product pairs per unit time) are k1 through k8, binding constants (i.e., number of substrate molecules bound into
enzyme-substrate complexes per enzyme per unit time) are kb1 through kb8 and unbinding constants (i.e., number of enzymesubstrate complexes dissociating into enzyme/substrate pairs per unit time) are kub1 through kub8.
Closed, isolated systems in which all processes are reversible obey the detailed balance condition, which requires that thesystem’s time evolution converge to a unique, stable equilibrium point (Kampen, 2007). In our system, the unbinding of productfrom enzyme is irreversible, so detailed balance is violated and oscillations are possible. Using chemical reaction network theory(Feinberg, 1987a) it can easily be shown that reversible product unbinding will eliminate the possibility of oscillations. For similarreasons, replacing the kinase and phosphatase by a single bifunctional enzyme capable of catalyzing both reactions reintroduces
detailed balance and precludes oscillations.
Simplified System with Michaelis-Menten Approximation
It can easily be shown that this system of equations is not independent, and can be reformulated as a system of 11 ODEs and threealgebraic equations describing the conserved total levels of substrate, kinase, and phosphatase. The system can be simplified further
if one makes the Michaelis-Menten approximation (Michaelis and Menten, 1913), in which the binding and unbinding of substrateenzyme complexes is assumed to be in instantaneous equilibrium, so that the concentrations of the enzyme-substrate complexes[ES00a], [ES00b], [ES01], [ES10], [FS01], [FS10],[FS11a] and [FS11b] can be replaced with their steady-state expressions (Figure 1C).
Kmx =ðkubx + kx Þ
; x ˛ 1; 2; 3; 4; 5; 6; 7; 8
kbx
Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors S1

d½S00 k1
k2
k5 ½F½S01  k6 ½F½S10 =++½S00 ½E +Km1 Km2dt
Km5
Km6
d½S01  k3 ½S01 ½E k5 ½F½S01  k1 ½E½S00  k7 ½F½S11 
=

+
+
dt
Km3
Km5
Km1
Km7


d½S11 k7
k8
k3 ½E½S01  k4 ½E½S10 ++=½S11 ½F +Km7 Km8dt
Km3
Km4½S10  = Stot  ð½S00  + ½S01  + ½S11  + ðEtot  ½EÞ + ðFtot  ½FÞÞ
½E =

Etot
ð1 + ½S00 =Km1 + ½S00 =Km2 + ½S01 =Km3 + ½S10 =Km4 Þ

½F =

Ftot
ð1 + ½S01 =Km5 + ½S10 =Km6 + ½S11 =Km7 + ½S11 =Km8 Þ

The resulting equations represent the system as a system of three ODEs and three algebraic equations; it is possible to eliminate
the algebraic equations by substitution, but the resulting nonlinear expressions are unwieldy and do little to aid in understanding thesystem, so we do not present them here. The binding and unbinding rate constants have been replaced by the binding constants Km1
through Km8. Kmi is the substrate concentration at which the rate of reaction i reaches half of its maximum; if binding is strong then Kmi
will be low because the enzyme-substrate binding is easily saturated. Numerical integration indicates that this reduced system gives
identical results to the full dynamical system for physically-reasonable values of the rate constants.
Stability of the Single-Site Phosphorylation System(See ‘‘Multisite Substrate Modifications Can Oscillate with Two Opposing Enzymes’’ in Results section).
If only a single phosphorylation site is present, then the system of mass-action equations can be written as:
d½S0 
=  kb0 ½E½S0  + kub0 ½ES0  + kF ½FS1 dt
d½S1 
=  kb1 ½F½S1  + kub1 ½FS1  + kE ½ES0 
dt
d½ES0 = kb0 ½E½S0   ðkub0 + kE Þ½ES0 
dt
d½FS1 = kb1 ½F½S1   ðkub1 + kF Þ½FS1 dt
d½E=  kb0 ½E½S0  + ðkub0 + kE Þ½ES0 dt
d½F=  kb1 ½F½S1  + ðkub1 + kF Þ½FS1 
dt[S0] and [S1] represent the concentrations of the unphosphorylated and phosphorylated forms of the substrate; [ES0] and [FS1] are
the concentrations of their complexes with the kinase (E) and the phosphatase (F), respectively. If we use the conservation relations[E] = Etot – [ES0] and [F] = Ftot – [FS0], then the system can be simplified:
d½S0 
=  kb0 ðEtot  ½ES0 Þ½S0  + kub0 ½ES0  + kF ½FS1 dt
d½S1 
=  kb1 ðFtot  ½FS1 Þ½S1  + kub1 ½FS1  + kE ½ES0 
dt
d½ES0 = kb0 ðEtot  ½ES0 Þ½S0   ðkub0 + kE Þ½ES0 
dt
d½FS1 = kb1 ðFtot  ½FS1 Þ½S1   ðkub1 + kF Þ½FS1 dt
For typical conditions, we can make the Michaelis-Menten approximation and set:
½ES0  =kb0 Etot ½S0 
kub0 + kE + kb0 ½S0 ½FS1  =kb1 Ftot ½S1 
kub1 + kF + kb1 ½S1 

S2 Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors

The system can now be simplified to only two equations:d½S0 
kb0 kE Etot ½S0 
kb1 kF Ftot ½S1 =+dt
kub0 + k E + kb0 ½S0  kub1 + kF + kb1 ½S1 d½S1 
d½S0 
=dt
dt
Clearly, [S0]+[S1] is a constant of the motion, and we can write [S0] + [S1] = Seff. This allows us to further simplify to:d½S0 
kb0 kE Etot ½S0 
kb1 kF Ftot ðSeff  ½S0 Þ+=dt
kub0 + k E + kb0 ½S0  kub1 + kF + kb1 ðSeff  ½S0 Þ
This is a one-dimensional system, and therefore cannot contain any attractors other than fixed points (Strogatz, 1994). Since all
trajectories are bounded (due to the conservation relations) the only possible long-term behavior is convergence to a stable fixed
point. In other words, limit-cycle oscillations are impossible as long as the Michaelis-Menten approximation is valid.
If the Michaelis-Menten approximation is not made, it appears that much less can be said. Using the tools of chemical reaction
network theory (CRNT) (Craciun and Feinberg, 2005, 2006; Feinberg, 1987b, 1988), the system
S0 + E%ES0 /E + S1
S1 + F%FS1 /F + S0
can be shown to have a deficiency of 1 and to be injective, meaning that it exhibits a single unique physically-relevant steady state(i.e., one in which all concentrations are real and positive). A topological analysis based on CRNT is not, however, able to determine
the conditions under which this steady state might be unstable and a limit cycle might be present. By comparison, the reaction
network described in Figure 1 has a deficiency of 3 and is not injective; multistationarity cannot be ruled out on topological grounds
alone. CRNT calculations were carried out using the CRN Toolbox (http://www.che.eng.ohio-state.edu/feinberg/crnt/).

Parameter Search
The system with dual phosphorylation sites is too large and complex for straightforward analytical study, so we chose to identify solutions numerically through an extensive search of parameter space. Our underlying assumption is that the biological plausibility of an
emergent property (in this case, sustained oscillations) is proportional to the probability of observing this property in an ensemble ofsystems whose parameters are randomly selected. This is an approximation to the real situation, in which the existence of a biologicalfeature depends on its emergence via a series of random mutations, but a random search over parameter values is far more theoretically tractable than a random search in genotype space.
Solutions were generated by selecting parameters from an exponential distribution bounded to the interval 1-1000 min-1 for k1
through k8 and 0.01-1000 mM for Km1 through Km8. While the presence of oscillations can be determined by linearizing the systemof equations about the fixed point and determining its stability (see Stability Analysis below), in practice it was more computationally efficient to simply integrate the mass-action equations and determine whether they converge to a constant value or a limit
cycle. Integration used a fourth-order stiffly-stable Rosenbrock scheme (Press, 2007). Evidence of multistability was neverobserved, but multistability could not be ruled out conclusively. Integrations always began with [S00] = 104, [E] = 1, [F] = 1, andall other concentrations set to zero; any stable attractors whose basin of attraction does not include this point would not havebeen observed.

Stability Analysis
A necessary condition for the presence of sustained oscillations is that the dynamics of the system converge to a limit cycle rather
than a fixed point. The conservation relations for S, E, and F make it clear that all trajectories are bounded – if the fixed points of thesystem can be shown to be unstable, motion on a higher-dimensional attractor is the only other possibility. (While chaotic motion ona strange attractor is possible in systems with more than two dimensions, sustained aperiodic trajectories were not observed in any ofthe parameter sets examined in this study. Once initial transients have disappeared, the system is effectively two-dimensional.) Thiscan be done by finding a (real, positive) set of concentration values such that all time derivatives are equal to zero and linearizing thesystem about this point by calculating the Jacobian (Strogatz, 1994). If the real part of at least one eigenvalue of the Jacobian ispositive, then the fixed point is unstable. For some parameter sets, an unstable fixed point coexists with a fixed point to which all
trajectories converge, resulting in a false positive for instability (and hence oscillation) if the unstable solution is used to generate
the Jacobian. Rather than attempt to exhaustively determine the roots of the nonlinear system of equations (which can in general
be quite challenging) we observed that bona fide oscillatory regions in parameter space are often bounded by a Hopf bifurcation,in which two complex-conjugate eigenvalues of the Jacobian cross the imaginary axis at the same time. It is not uncommon,
however, for the (positive) real parts of these eigenvalues to diverge significantly as the system moves further from the Hopf
Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors S3

bifurcation. By requiring that two of the eigenvalues be greater than zero, oscillatory and nonoscillatory parameter sets could be
reliably distinguished.
Clustering and PCA Analysis
Parameter sets were categorized using quality threshold (QT) clustering (Heyer et al., 1999). Given a metric for determining the
distance between two points in parameter space, the QT algorithm finds the largest set of points such that all are within a distanced of each other. These are then removed from consideration, and the largest cluster within the remaining collection of points is
identified. This process continues until all points have been assigned to a cluster. This process becomes rather memory-intensive
as the size of the data set to be clustered grows. A two-step process was employed for very large data sets in which the full QT
algorithm was run on an initial seed set (10,000 points) and each of the remaining points was assigned to the closest cluster.
The results obtained by the two-step method for large data sets were qualitatively similar to those obtained by the full QT method
for smaller data sets.
The distance metric was designed to take several factors into account. First, because the parameters were exponentially distributed in the random search, the distance metric was calculated on the logarithms of parameters. This ensures that (for example) a pointwith k1 = 900 and a point with k1 = 1000 are more similar than a point with k1 = 1 and a point with k1 = 101. In Figure 1B, the essential
features of the reaction network are unchanged if S10 and S01 are exchanged. This means that, for any two points in parameter space,
two possible distances between them exist, of which the shortest was chosen. Lastly, the rate constants {ki} have units of inverse
time, and rescaling {ki} by a constant factor has no effect on the resulting dynamics except for the rescaling of the time axis. Therefore,
when each pairwise distance was calculated, rate constants were rescaled to minimize the distance metric. These measures ensurethat large distances between parameter sets will correspond to major differences in the oscillatory dynamics.
The choice of cluster diameter is somewhat arbitrary; two major clusters emerged for a range of cluster diameters and a diameter
was chosen at which these two clusters accounted for 70% of solutions (Figure 2B; Figure S2A). Each cluster was ‘‘aligned’’ byidentifying the cluster member with the shortest average distance to others in the cluster and choosing the symmetry orientationof each cluster that minimizes the distance to this reference member. The parameter distributions in aligned clusters can then be
calculated; representative histograms for the major clusters are shown in Figure 2C.
The remaining small clusters (30% of the total) are generally similar either to Cluster 1 or Cluster 2; if all parameter sets assigned
to neither Cluster 1 nor Cluster 2 are grouped together, their parameter distribution is similar to Cluster 1 or Cluster 2. For example,
Figure S2B shows both Cluster-2-type characteristics (high k3, k6, low k4, k5, low Km1, Km8) and Cluster-1-type characteristics (highk1, k8, low k2, k7, low Km4, Km5). To a first approximation, therefore, solutions can be categorized as being similar to either Cluster 1 orCluster 2.
The presence of two distinct oscillatory clusters in parameter space can also be confirmed using principal components analysis(PCA). If a covariance matrix is constructed from the logarithms of the oscillatory parameter values and normalized so that eachcolumn has zero mean and unit variance, the eigenvalues of this matrix can be used to reduce the dimensionality of parameter spaceby projecting it down onto the (orthogonal) directions that contain the largest share of the sample variance. When a random sample ofoscillatory parameter sets is projected down onto the first two principal components (which together account for about 30% of thevariance), two distinct clusters are visible, separated by an empty space (Figures S2C and S2D). Identifying the points on the PCA
projection according to the cluster to which they were assigned by the QT algorithm makes it clear that Cluster 1 and Cluster 2 arequalitatively distinct and that no other major groupings are present.
Detailed Analysis of Oscillation Mechanism(See ‘‘Period- and Amplitude-Determining Processes’’ in Results section.)
Figure 4B divides the oscillation cycle of a typical Cluster 1 oscillator into two phases: one involves the net conversion of S10 to S01(via S00) and is shown against a white background; the other involves the conversion of S01 back to S10 via S11 and is shown against
a shaded background. The system is clearly symmetric (in the left side of Figure 4B) – the two phases look identical if S10 and S00 are
exchanged with S01 and S11, respectively. In addition, the sequestration patterns of the kinase and phosphatase (lower two panels)
have identical shapes, except for a 12 hr phase shift – both enzymes spend most of their time bound to a single substrate and associate with other binding partners only during brief antiphasic episodes.
At the beginning of the S01/S11/S10 phase (colored background), S01 is highly abundant and saturates the phosphatase almost
completely, rendering it nearly inactive. Because S01 is the only highly abundant species at this point, it also binds most of the kinase,despite the fairly low affinity of the kinase for S01. Unlike the phosphatase, the kinase reacts efficiently when bound to S01, causing S01
to be converted to S11 at a respectable rate. As S11 accumulates, some of it is converted to very small amounts of S10 –initially S10 ishardly visible in the top panel (bulk protein concentrations), but the affinity of the kinase for S10 is so high that it still rapidly sequesters
most of the kinase away from S01. Once the kinase sequestration by S10 saturates (at around the middle of this phase), the depletionof S01 slows substantially. Once the levels of S01 have nearly reached zero, something interesting finally happens with the phosphatase. After spending nearly all of this phase bound unproductively to S01, it is finally free to seek other partners; the most viablecandidate is S11, now the most abundant species in the system. The association of phosphatase and S11 is, however, destined to
be short-lived, as S11 is precipitously converted to S10, completing the net conversion of S01 to S10.
S4 Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors

The S10/S00/S01 phase follows precisely the same pattern. Initially, both the kinase and the phosphatase are bound primarily by
the highly-abundant S10. The S10/S11 reaction is inefficient, so the kinase remains unproductively bound to S10 through most of thisphase. The phosphatase, however, rapidly converts S10 to S00. With the abundant but unattractive S10 out of the picture, the phosphatase can then bind almost exclusively to the trace amounts of S01 that have begun to appear, rendering it nearly-inactive for the
rest of this phase. This causes the rate of the S10/S00 conversion to slow substantially, with the slow decline of S10 effectively determining the length of this phase. Once it has disappeared completely, the kinase has nowhere to turn except for the now-abundant
S00, which is rapidly consumed to produce S01. Once S01 again becomes the dominant species in the system, the cycle is ready to
begin anew.
In the right panel, the symmetry of the system is broken by halving the rate of the S01/S11 conversion. The concentration traces in
the S10/S00/S01 phase (white background) are nearly identical to the symmetric case, indicating that this particular reaction has
a minimal impact on this phase. The S01/S11/S10 phase (colored background), on the other hand, nearly doubles in length. Close
examination shows that the shape of the transition in phosphatase binding near the end of this phase is unchanged – most of thelengthening comes from an extension of the S01/S11 conversion. This explains the large impact of this parameter on the oscillation
period – the phosphatase-catalyzed steps all wait for the sequestration by S01 to be lifted, so changes that slow down the depletion of
S01 will impact the period much more than changes that affect the more transient reactions.
Detailed Discussion of Bifurcations(See ‘‘Period- and Amplitude-Determining Processes’’ in Results section).
Determining the precise nature of bifurcations in a complex dynamical system can sometimes be challenging, and our identificationof the bifurcations illustrated in Figure 4A as Hopf bifurcations is based solely on numerical evidence. In each case, the nature of thebifurcations was determined by observing the behavior of the eigenvalues of the Jacobian at the fixed point as parameter values werevaried. The symmetric reference parameter set for Cluster 1 (see Table S1) shows two nondegenerate, positive, real eigenvalues. As
the bifurcation point is approached, these two eigenvalues approach each other and real the imaginary axis, so that they becomecomplex conjugates with positive real parts. At the bifurcation point, these two eigenvalues simultaneously cross the imaginary
axis, the typical signature of a Hopf bifurcation. Additional evidence comes from the scaling of the oscillation period with r, wherer is the distance from the bifurcation point. While the period can vary dramatically with parameter changes, it tends to remain constantas r becomes very small, another typical feature of Hopf bifurcations. Finally, the amplitudes of oscillations typically decrease as r1/2
as the bifurcation point is approached. This can be demonstrated for the Km4 bifurcation shown on the right-hand side of Figure 4A,
and for many of the other bifurcations shown in Figure S4. Unfortunately, this characteristic was observed in neither the k3 bifurcation
shown in Figure 4A nor the similar k1/k8 bifurcations for Cluster 2.
Determining the Typical Robustness of Oscillation to Parameter Changes(See ‘‘Period- and Amplitude-Determining Processes’’ in Results section).
In addition to detailed analysis of a few example cases, one can also use bifurcation analysis to examine the robustness of a larger
collection of parameter sets. We chose 1000 parameter sets each from Cluster 1 and Cluster 2 and determined the range overwhich each parameter could be varied before oscillations were lost. Figure S4B shows the average oscillatory ranges for the two
major clusters. Because we only examined bifurcations lying between 103 and 104, these are conservative estimates; Figures 4A
and S4 suggest that parameters can sometimes be varied far outside this range without affecting oscillations. On average, the parameters of Cluster 1 oscillators can tolerate about 1000-fold variation before oscillations are lost; Cluster 2 oscillators are slightly lessrobust and can tolerate roughly 800-fold variation.
Sensitivity of the PTO Period to Total Substrate Concentration(See ‘‘Period- and Amplitude-Determining Processes’’ in Results section.)
Close examination of the bifurcation diagram for Stot (Figure S4A) raises an important concern that should be addressed. Whileoscillations are lost for Stot ( 3; there is no upper bound to the oscillatory regime. For Stot[ 3, the period is almost exactly linear
in the Stot/Etot ratio (or, equivalently, Stot/Ftot). Adjusting this ratio amounts to rescaling the time axis of the simulations; all steps ofthe oscillation cycle are equally affected. When the substrate/enzyme ratio is large, enzymes will be operating near saturation,
meaning that all catalytic rates are roughly equal to the product of substrate and enzyme concentrations, and a change will affect
all reactions equally. Near the bifurcation point, enzymes are no longer saturated and this relationship breaks down.
The period of the circadian clock has been shown to be largely insensitive to the global rate of mRNA transcription (Dibner et al.,2009). Is this robustness captured in our simple model? If transcriptional repression decreases the concentrations of substrate,kinase, and phosphatase by the same factor, then the dynamics of the oscillator will be completely unchanged. Some studies(Coté and Brody, 1986), however, indicate that the circadian period may also be insensitive to a selective change in the expression
level of PER, a major phosphorylation target of CKIε/d and the most likely circadian analog to the substrate in our model system.
This difficulty can be avoided if the effective substrate/enzyme ratio is kept at a constant level regardless of the total expression
level of each substrate and enzyme. For example, subcellular compartmentalization could be dependent on the formation of a tightly
bound complex between enzymes and substrates – if a ‘‘date’’ is required for admission, then stoichiometry will be well-controlled. In
Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors S5

the circadian system, localization of PER to the nucleus requires participation in a multisubunit complex, often including CKIε/d andphosphatase (Lee et al., 2001).
Stochastic Effects on the Stability and Period of PTO(See ‘‘Robustness of the Oscillation Period against Stochastic Fluctuations’’ in Results section).
To get a sense of whether the results depicted in Figure 5 are typical of Clusters 1 and 2 more generally, we used the stochastic
simulation algorithm (SSA) to estimate what fraction of oscillators remain viable as the number of enzymes is decreased; roughly 10%
remain oscillatory even with only one of each type of enzyme present (Figure S5C). In this sense, Cluster 1 is slightly more robust than
Cluster 2, although the differences are small. Interestingly, the two clusters show opposing behaviors as molecular noise becomes
more pronounced – the period of Cluster 1 oscillators tends to increase, while that of Cluster 2 oscillators tends to decrease(Figure S5D).
Biological Significance of the Tunable Period Length
The bifurcation diagram for k3 in Figure 4A makes it clear that a wide range of periods are available in certain parameter regimes. Thisimplies the presence of broadly-applicable design features that might apply both in the circadian field and other oscillatory systems
with shorter or longer time scales. Comparison with the bifurcation diagrams shown in Figure S4 shows that, for most parameters, theperiod variation in the vicinity of the reference parameter set is rather modest – an oscillator based on this scheme is likely to bereasonably mutationally robust. Although tuning of the oscillation period is possible (and likely to be desirable) on evolutionary timescales, Figure 5A shows that the periods of individual oscillators are likely to be fairly consistent, even in the presence of stochastic
noise. These features are consistent with the current understanding of circadian rhythmicity: the (mammalian) circadian clock isrobust against temperature (Hastings and Sweeney, 1957; Pittendrigh, 1954) and global transcription rate (Dibner et al., 2009),
but changes in the activity of CKIε/d are a fragile point (Gallego and Virshup, 2007; Isojima et al., 2009; Lee et al., 2009). CKI-catalyzed
reactions likely play a period-determining role similar to that of the reactions with rates k3 and k6 in our much simpler system.
Calculating Temperature Sensitivity
For biochemical reactions, temperature sensitivity is often described using the temperature coefficient Q10:

Q10 =

R2
R1

10=

T2 T1;

where R1 and R2 are the rates of a process at temperatures T1 and T2. If Q10 = 2, the reaction rate doubles over an interval of 10
degrees. Note that, if a reaction has a constant Q10, its rate will depend exponentially on the temperature. If the reaction rate is scaled
DT=10such that R2 = aR1, the rescaling constant can be expressed as a = Q10 .
Rather than explicitly estimating activation energies, the effects of temperature were introduced by assuming Q10 values for the
kinase and the phosphatase (Q10(E) and Q10(F)) and using these to calculate rescaling factors aE and aF at a given temperature difference DT. The resulting changes in the oscillation period t(DT) could usually be fit very well with an exponential curve of the formt = t0 ebDT . The temperature coefficient of the entire cycle can then be expressed as:ðcycleÞ

Q10


=

R2
R1

10=

DT

=

 10=
DT
t1
bDT 10=
DT = e10b=et2

For a given parameter set, the rate constants were rescaled for nine points evenly spaced along a temperature range of ± 5 , theperiod was calculated for each temperature point, and an exponential fit was used to obtain Q10(cycle). Binding-unbinding reactions
had a minor impact on temperature sensitivity; we assumed a fairly high Q10 value of 3.0. Most fits were of excellent quality; those with
R2 < 0.9 were discarded (Figures S6A–S6C).
Detailed Explanation of Temperature Compensation(See ‘‘The Period of a Simple PTO Can Be Temperature Compensated’’ in Results section.)
Bona fide circadian clocks involve more than merely autonomous oscillations. A cellular rhythm is generally not considered a circadian clock unless it also exhibits temperature compensation, in which the oscillation period is roughly unchanged by the ambient
temperature. The temperature sensitivity of a system is often described using the temperature coefficient Q10, which describes
the factor by which the temperature changes as the temperature is raised by 10 C. Most biochemical reactions have a Q10 valueof 2-3, meaning that their rates increase between two and three times in response to a 10 C temperature change. For a temperature-compensated process, Q10 y 1; measured Q10 values for circadian clocks range from about 0.8 to 1.4 (Dunlap et al., 2004).
To further deepen the analogy between our simple model system and the circadian clock, we examined its potential for temperature compensation. This was done by assuming temperature coefficients Q10(E) for the kinase and Q10(F) for the phosphatase, andusing these to modify rate constants in accordance with an assumed external temperature (see Experimental Procedures; Figures
S6A and S6B).
S6 Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors

It is obvious that if Q10(E) = Q10(F) = 2, every reaction rate in the simple PTO is doubled when the temperature changes by 10 C.
In that case, the period is cut in half and temperature compensation is impossible; equal Q10 values only permit temperaturecompensation if both are equal to one. Because we are most interested in the potential behavior of systems with temperatureinsensitive kinases such as CKIε/d, we have generally assumed Q10(E) = 1. Of course it is entirely possible that Q10(F) = 1, and thiswould provide a very elegant mechanism for temperature compensation. For the sake of generality, however, we consider thecase in which Q10(F) = 3, and ask whether temperature compensation is possible even in the presence of a fairly temperature-sensitive
phosphatase.
The distribution of Q10(cycle), the temperature coefficient of the oscillation period was calculated as shown in Figure 6A: Q10(cycle)
calculations were attempted for 5 3 106 oscillating parameter sets. Many parameter sets show oscillations that are robust totemperature change, while most were unsuccessful. The systematic changes in enzyme reaction rates (induced by temperature
differences) can involve fairly significant motion in parameter space, often bringing the system outside the oscillatory regime. For
Q10(E) = 1 and Q10(F) = 3, Q10(cycle) could be calculated for 106 parameter sets. In general, larger values of Q10(E) and Q10(F) cause
the peak of Q10(cycle) to be shifted to higher values, and the width of the Q10(cycle) distribution increases as the difference between
Q10(E) and Q10(F) increases (Figures 6B and S6C). It is possible for the oscillation period to be temperature compensated or even overcompensated as long as Q10(E) s Q10(F); the probability of temperature compensation increases when values of Q10(E) or Q10(F) isclose to 1.
We next sought to determine whether temperature-compensated oscillators share characteristic parameter distributions. For thecase with Q10(E) = 1 and Q10(F) = 3 (i.e., a temperature-insensitive kinase and a temperature-sensitive phosphatase), the collection ofparameter sets was sorted according to their calculated Q10(cycle) values and separated into five groups numbered from G1 (mostovercompensated; period increases with increasing temperature) to G5 (most undercompensated; period decreases with increasing
temperature) (Figure 6C). Parameter distributions similar to those in Figure 2C could then be calculated for each group. A few representative traces are shown in Figure 6D; see also Figures S6D and S6E. The parameter distributions vary smoothly from G1 to G5 andshow a striking symmetry property: the G1 distributions are nearly identical to the G5 distributions upon an exchange of S00 and S11(reversing the roles of the kinase and phosphatase).
Figure 6E distills these parameter distributions into a schematic form. In the undercompensated case (G5), the fastest reactions inthe system are the kinase-driven forward reactions (S00/S01/S11), while the phosphatase-driven reactions (S11/S10/S00) areslower and nearly irreversible. The temperature-insensitive kinase reactions have a negligible impact on the period, while the temperature-sensitive phosphatase steps are rate-limiting, leading to overall strong temperature dependence. In addition, Km8 is typicallylower than Km6 in G5; this means that the phosphatase often will not bind to S10 until S11 has been largely depleted, introducing an
additional sequestration that favors temporal ordering of the slow phosphatase steps. Km1, however, is not typically lower than Km3,
so the fast kinase steps tend to happen concurrently. In the overcompensated G1 case, the opposite is true. The temperature-insensitive kinase reactions are slow, irreversible, and temporally-ordered, giving them a disproportionate impact on the period, while thetemperature-sensitive phosphatase steps are fast and concurrent, and have little effect on the period. It should be stressed that overcompensation is probably not due to an excessive period dependence on the temperature-insensitive kinase reactions; if the perioddepended only on these reactions then the cycle would be perfectly compensated.
Both overcompensation and extreme undercompensation (i.e., Q10(cycle) > Q10(F)) can best be understood by examining the bifurcation diagrams in Figures 4B and S4. If the system is near a bifurcation point, then small changes in parameters can result in large
period changes. The fact that 80% of Q10(cycle) calculations failed suggests that a large fraction of the randomly-discovered oscillators are near the edge of the oscillatory regime. This is a consequence of the geometry of high-dimensional space: only 19% of thearea of a two-dimensional circle is located within 10% of its edge, while 81% of the volume of a 16-dimensional ball is located within
10% of its boundary. The oscillatory parameter regime is likely to be nonspherical, further increasing the fraction of its volume locatedclose to a boundary. This intuition is supported by the fact that most of the oscillators with extreme Q10(cycle) values also showed verysmall oscillation amplitudes, as is often true near bifurcations. Most of the tails of the distribution can be attributed to low-amplitude,
marginally-stable oscillators. We would therefore expect a viable temperature-insensitive oscillator to be one in which the temperature-insensitive steps are rate-limiting (as described above) and the system is located far from the boundaries of the oscillatory
regime, so that amplitudes remain large and undercompensation is avoided.
Temperature compensation can also be observed in stochastic simulations. The Q10(cycle) values obtained from deterministic
calculations are not always valid in the stochastic case, because system size can affect the oscillation period (Figure 5A). A systematic search revealed several parameter sets for which stochastic simulations are nearly temperature-insensitive; one example isshown in Figure 6F (see Table S3). In this calculation, S00 and S10 are present only in very brief spikes, consistent with fast reactions.
Most of the time in the cycle is taken up by the slow, nearly-irreversible S01/S11 conversion, which is catalyzed by the temperatureinsensitive kinase.
In summary, precise temperature compensation may not be particularly common, but it is certainly possible. Roughly 3.4% of
oscillators show a Q10 value between 0.9 and 1.1. Most importantly, temperature compensation does not require highly unusualparameter values or modified network topologies. Instead, a range of temperature sensitivities is possible, and temperature insensitivity occupies an intermediate position in this range. To ensure temperature compensation, the temperature-insensitive steps
should be rate-limiting – this is realized through slow, irreversible reactions and a tendency toward temporal ordering through differential binding constants. Interestingly, the phosphorylation of PERs by CKIε/d is apparently a rate-determining reaction in the
Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors S7

mammalian circadian clock and is temperature insensitive when PER-derived peptide is used as a substrate in vitro (Isojima et al.,2009).
Detailed Discussion of Other Posttranslational Oscillator Models(See ‘‘Comparison with Other PTO Models’’ in the Discussion section.)
Several models have been described for bistability and oscillations in the MAPK system, and our system exhibits important similarities and differences with these. For example, Markevich et al. (2004) describe a single level of a MAPK cascade, in which a single
kinase (MAPKK) and a single phosphatase (MKP) catalyze the interconversions of MAPK enzymes with zero, one, or two phosphate
modifications. The most important difference between this scheme and the one described here is that MAPK models (following theoriginal formulation by Huang and Ferrell (Huang and Ferrell, 1996)), typically do not distinguish between the two possible singlyphosphorylated forms of MAPK, assuming instead that phosphorylation follows a strict ordering. Interestingly, Markevich et al.
describe a pair of feed-forward inhibition motifs that look strikingly similar to those presented in Figure 2. Unless a distinction ismade between the two singly-phosphorylated states, however, a directional rate bias (our first design motif) is not possible, and oscillations are not observed. Several studies have observed that, while bistability is possible within a single level of the cascade, feedback loops in a multilevel system are typically required for oscillations (Chickarmane et al., 2007; Kholodenko, 2000; Shankaran et al.,2009). Liu and coworkers (Liu et al., 2011) have shown that oscillations can be generated in a similar strictly-ordered two-site phosphorylation system with implicit feedback similar to the enzyme sequestration observed in our study, with the major difference that
their system requires a fourth component (in addition to oscillating substrate, kinase, and phosphatase) that is actively produced anddegraded (Liu et al., 2011).
Another essential difference between our study and the various MAPK models in the literature is that these models are all intended as descriptions of an experimentally-characterized system and therefore have access to parameters that have been (atleast approximately) experimentally-determined. For example, Qiao and co-workers (Qiao et al., 2007) used a Monte Carlo parameter search similar to ours to classify the dynamical behaviors available in the complete MAPK cascade system, with no explicit
feedback between the MAPK level and the MAPKKK level. By sampling 2 3 105 exponentially-distributed parameter sets, they
were able to determine the relative frequencies of single-valued, hysteretic, and oscillatory input-output relationships. In contrast
to our approach, however, in which parameter values were selected from a range covering several orders of magnitude, theirparameter values were selected to be within 25-fold of the experimentally estimated values, meaning that they characterized
the behaviors present within a reasonable neighborhood in parameter space. Our approach, in contrast, has been to sample
the parameter space as completely as possible and analyze the resulting parameter combinations to determine which regions
harbor interesting dynamical behaviors.
Our model can also be compared fruitfully to various models that exist for the in vitro reconstitution of the cyanobacterial circadianclock. Van Zon et al. (van Zon et al., 2007) describe a model in which each monomer in a KaiC hexamer can be either phosphorylatedor unphosphorylated, and the ring as a whole is in either an active state (favoring phosphorylation) or an inactive state (favoring
dephosphorylation). Highly-phosphorylated hexamers easily become inactivated and weakly-phosphorylated hexamers tend to
become activated, leading to a preferred ordering of the system’s states (our first design motif). The population of single-hexamer
oscillators is synchronized through the action of KaiA, which is assumed to bind most strongly to weakly-phosphorylated hexamers,
promoting phosphorylation (they refer to this mechanism as ‘‘differential affinity’’). This leads to a situation similar to our sequestration-based synchronization – the hexamers that are ahead of the pack are less successful in competing for KaiA and will tend to slow
down, while the laggards will tend to speed up. The supplementary information of their study formulates differential affinity in very
general terms that are easily applicable to our system.
The circadian clock model of Van Zon and co-workers also discusses a hypothesis for temperature compensation that has important similarities to ours. They begin with the assumption that the various (de)phosphorylation rates are all temperature-insensitive,
while other factors (such as binding/unbinding and enzyme/substrate affinities) may be temperature sensitive. In order for the system
as a whole to remain temperature-compensated, the temperature-insensitive steps must be rate limiting; from this observation they
infer requirements that the (un)binding and conformational change rates should be much faster than the catalytic rates and that affinities between the Kai proteins should be high enough that binding is typically saturated. A more detailed study of temperaturecompensation in the KaiABC system was presented by Hatakeyama and Kaneko, who argued that temperature compensation could
be achieved (even if individual reactions are temperature-sensitive) as long as reaction rates are limited by an enzyme that is sharedacross several reactions (Hatakeyama and Kaneko, 2012).
In our study, we assumed only that the kinase was temperature-insensitive, and examined conditions under which temperaturecompensation can still be achieved with a temperature-sensitive phosphatase. This complicates the picture substantially, but therequirement of nearly-saturated binding and fast noncatalytic processes is consistent with our results.
Another approach to the modeling of the cyanobacterial circadian clock is illustrated by Clodong and co-workers (Clodong et al.,2007), who began with a loop topology similar to the one used by van Zon et al., except that the distinction between active and inactive conformations is replaced by a distinction between KaiC hexamers without KaiB (which will tend to be phosphorylated) andKaiBC complexes which will tend to be dephosphorylated. They then used a ‘‘rewiring’’ algorithm to add additional positive or negative regulatory connections to this loop topology. The topology that generated the highest-amplitude, most robust oscillations did sousing a feed-forward inhibition synchronization scheme similar to ours, in which KaiBC complexes with low phosphorylation number
S8 Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authorsinhibit the phosphorylation of KaiC complexes. This can be accomplished if KaiA binds strongly to low-phosphorylated KaiBC
complexes, making it unavailable to promote the phosphorylation of the KaiC hexamers.
The model of Rust et al. (Rust et al., 2007) describes a two-site system similar to ours, in which the two phosphorylation states ofthe substrate (KaiC) have distinct roles. Several differences exist in the detailed implementation, however – their model lacks a separate phosphatase (relying instead on the autophosphatase activity of KaiC) and the autokinase activity of KaiC is modulated by theactivities of the regulatory proteins KaiA and KaiB. KaiA promotes the preferential phosphorylation of KaiC on T432 (KaiC-T), followedby phosphorylation on S431 (KaiC-ST); the phosphorylation on S431 to form KaiC-S is slower, and its dephosphorylation from KaiCST is inhibited by KaiA. Once low levels of KaiC-S do form, however, they recruit KaiB, which inhibits the activity of KaiA. In this way,
KaiC-S indirectly catalyzes its own formation, as well as inhibiting the formation of KaiC-T and KaiC-ST. Similarities to our model exist
in the general principles of kinetic bias (i.e., T432 is phosphorylated more quickly than S431) and synchronization by feedforward
inhibition, although the molecular details are quite different.
The various cyanobacterial clock models show differences among themselves, but the basic design principles identified in thisstudy can be discerned in each. In all cases, a kinetic bias (sometimes enforced by the presence or absence of KaiB or by an active
or inactive conformational state of KaiC) is used to give states a temporal ordering, consistent with our first design motif. In addition,
all of the models described involve some variation on the theme of synchronization by feedforward inhibition, with high concentrations of a particular phosphorylation state inhibiting the progress of reactions that are a few steps ahead in the cycle. The moststriking difference lies in the physical system – cyanobacterial clock models focus on a specific biochemical system and their parameter choices are guided by experimental evidence, while our model focuses on a very generic system and design elements werediscovered via a broad, unconstrained parameter search. Most importantly, our results demonstrate that these design principles
are not a result of specific features of the KaiABC system (such as auto(de)phosphorylation activity, multimeric arrangement of
substrates, or competitive regulation) but that they also appear in a much more general context.
More General Formulation of Design Principles(See ‘‘General Applicability of a Design Principle Composed of Two Design Motifs’’ in Discussion section.)
The general design principles uncovered in this study can also be applied to situations in which more than two phosphorylationsites are present (Figure 7B). Suppose that a protein has n phosphorylation sites. Out of the 2n possible phosphorylation states,
a subsequence of k % 2n states exists such that the rate constants in the system predispose individual substrates to visit these statesin the order S1,S2,..,Sk,S1,S2,.. This is equivalent to our first design motif. Enzyme sequestration requires that some phosphorylation
state Sj binds strongly and sequesters the enzyme that catalyzes the downstream reaction Sk/Sk+1, where k > j. This might also be
the same enzyme that catalyzes Sj/Sj+1, as in the case of our Cluster 2, where (Sj,Sk,Sk+1) corresponds to (S00,S01,S11) or(S11,S10,S00). Alternatively, it may be a different enzyme, as in the case of Cluster 1, where (Sj,Sk,Sk+1) corresponds to(S01,S11,S10) or (S10,S00,S01). It is not necessary that the system have only one type of kinase and one type of phosphatase, aslong as these conditions for enzyme competition are present. As can be seen from the examples developed in this study, oscillations
will be most robust when several overlapping modules exist. Geometric considerations (Manrai and Gunawardena, 2008;
Thomson and Gunawardena, 2009) have shown that competition for enzyme between different substrate phosphorylation states
can lead to multistability in multisite phosphorylation systems; our results lead one to imagine the presence of higher-dimensional
attractors as well.
We have described the model as a phosphorylation system because our initial interest in this system was motivated by the prevalence of multisite phosphorylation in circadian systems. However, we note that the same mathematical framework could be applied
to any type of multisite posttranslational modification such as methylation or acetylation.
Compatibility of PTO with TTO
While the system presented in this paper would be an adequate model of an in vitro experiment, realization in cellulo would require
the inclusion of transcription, translation, and degradation. Accounting for these processes would complicate the model substantially and is beyond the scope of this paper, but a few comments are in order. Figure 1B suggests that the system can be
abstracted into two processes – one is a production of S11 from S00 (with S01 as an intermediate), and the other is the disappearance of S11 and reappearance of S00. In our model, the second process is accomplished by phosphatase-catalyzed steps. This
could also easily be accomplished by a transcriptional-translational loop in which S11 is degraded and acts as a transcriptional
repressor of the gene encoding the substrate. The newly-translated substrate will be unphosphorylated, so this scheme includes
the necessary ingredients of S00 appearance, S11 disappearance, and coupling between these two processes (Figure 7C). In thiscase, the presence of the phosphatase-catalyzed pathway from S11 to S00 is superfluous, and its function can be replaced bydegradation of S11 and synthesis of S00. Thus, our PTO model can easily be generalized to a TTO model, and many of thesame considerations will apply.
This scheme is somewhat analogous to the current canonical model of the circadian clock in various organisms. Typically, thesubstrates that are phosphorylated at multiple sites and degraded upon progressive phosphorylation are also involved in the core
negative feedback loop – PER in fly and mammals, and FRQ in fungi (Dibner et al., 2010; Dunlap, 1999; Reppert and Weaver,
2002; Takahashi et al., 2008; Young and Kay, 2001). These proteins directly or indirectly repress their own transcription (Aronson
et al., 1994; Darlington et al., 1998; Kume et al., 1999), though it should be noted that the repressor activity of mammalian PER seems
Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors S9

to be relatively weak (Sangoram et al., 1998) and both mammalian PER and fungi FRQ can act as positive regulators under somecircumstances (Lee et al., 2000; Ogawa et al., 2011; Shearman et al., 2000). Multi-site phosphorylation, followed by degradation,
may also be experienced by these proteins (Baker et al., 2009; Chiu et al., 2011; Shanware et al., 2011; Tang et al., 2009; Vanselow
et al., 2006). The PTO oscillator presented in this paper is thus likely to be compatible with well-characterized TTO systems, and someof the same design principles may apply.
SUPPLEMENTAL REFERENCES
Aronson, B.D., Johnson, K.A., Loros, J.J., and Dunlap, J.C. (1994). Negative feedback defining a circadian clock: autoregulation of the clock gene frequency.
Science 263, 1578–1584.
Baker, C.L., Kettenbach, A.N., Loros, J.J., Gerber, S.A., and Dunlap, J.C. (2009). Quantitative proteomics reveals a dynamic interactome and phase-specificphosphorylation in the Neurospora circadian clock. Mol. Cell 34, 354–363.
Chickarmane, V., Kholodenko, B.N., and Sauro, H.M. (2007). Oscillatory dynamics arising from competitive inhibition and multisite phosphorylation. J. Theor.
Biol. 244, 68–76.
Chiu, J.C., Ko, H.W., and Edery, I. (2011). NEMO/NLK phosphorylates PERIOD to initiate a time-delay phosphorylation circuit that sets circadian clock speed. Cell
145, 357–370.
Clodong, S., Dühring, U., Kronk, L., Wilde, A., Axmann, I., Herzel, H., and Kollmann, M. (2007). Functioning and robustness of a bacterial circadian clock. Mol.
Syst. Biol. 3, 90.
Coté, G.G., and Brody, S. (1986). Circadian rhythms in Drosophila melanogaster: analysis of period as a function of gene dosage at the per (period) locus. J. Theor.
Biol. 121, 487–503.
Craciun, G., and Feinberg, M. (2005). Multiple equilibria in complex chemical reaction networks: I. The injectivity property. SIAM J. Appl. Math. 65, 1526–1546.
Craciun, G., and Feinberg, M. (2006). Multiple equilibria in complex chemical reaction networks: II. The species-reaction graph. SIAM J. Appl. Math. 66, 1321–
1338.
Darlington, T.K., Wager-Smith, K., Ceriani, M.F., Staknis, D., Gekakis, N., Steeves, T.D., Weitz, C.J., Takahashi, J.S., and Kay, S.A. (1998). Closing the circadian
loop: CLOCK-induced transcription of its own inhibitors per and tim. Science 280, 1599–1603.
Dibner, C., Sage, D., Unser, M., Bauer, C., d’Eysmond, T., Naef, F., and Schibler, U. (2009). Circadian gene expression is resilient to large fluctuations in overall
transcription rates. EMBO J. 28, 123–134.
Dibner, C., Schibler, U., and Albrecht, U. (2010). The mammalian circadian timing system: organization and coordination of central and peripheral clocks. Annu.
Rev. Physiol. 72, 517–549.
Dunlap, J.C. (1999). Molecular bases for circadian clocks. Cell 96, 271–290.
Dunlap, J.C., Loros, J.J., and DeCoursey, P.J. (2004). Chronobiology: Biological Timekeeping (Sunderland, MA: Sinauer Associates).
Feinberg, M. (1987a). Chemical reaction network structure and the stability of complex isothermal reactors. 1. The deficiency-zero and deficiency-one theorems.
Chem. Eng. Sci. 42, 2229–2268.
Feinberg, M. (1987b). Chemical reaction network structure and the stability of complex isothermal reactors. 1. The deficiency-zero and deficiency-one theorems.
Chem. Eng. Sci. 42, 2229–2268.
Feinberg, M. (1988). Chemical reaction network structure and the stability of complex isothermal reactors. 2. Multiple steady states for networks of deficiency
one. Chem. Eng. Sci. 43, 1–25.
Gallego, M., and Virshup, D.M. (2007). Post-translational modifications regulate the ticking of the circadian clock. Nat. Rev. Mol. Cell Biol. 8, 139–148.
Hastings, J.W., and Sweeney, B.M. (1957). On the mechanism of temperature independence in a biological clock. Proc. Natl. Acad. Sci. USA 43, 804–811.
Hatakeyama, T.S., and Kaneko, K. (2012). Generic temperature compensation of biological clocks by autonomous regulation of catalyst concentration. Proc.
Natl. Acad. Sci. USA 109, 8109–8114.
Heyer, L.J., Kruglyak, S., and Yooseph, S. (1999). Exploring expression data: identification and analysis of coexpressed genes. Genome Res. 9, 1106–1115.
Huang, C.Y., and Ferrell, J.E., Jr. (1996). Ultrasensitivity in the mitogen-activated protein kinase cascade. Proc. Natl. Acad. Sci. USA 93, 10078–10083.
Isojima, Y., Nakajima, M., Ukai, H., Fujishima, H., Yamada, R.G., Masumoto, K.H., Kiuchi, R., Ishida, M., Ukai-Tadenuma, M., Minami, Y., et al. (2009). CKIepsilon/
d-dependent phosphorylation is a temperature-insensitive, period-determining process in the mammalian circadian clock. Proc. Natl. Acad. Sci. USA 106,
15744–15749.
Kampen, N.G.v. (2007). Stochastic Processes in Physics and Chemistry, Third Edition (Amsterdam/Boston: Elsevier).
Kholodenko, B.N. (2000). Negative feedback and ultrasensitivity can bring about oscillations in the mitogen-activated protein kinase cascades. Eur. J. Biochem.
267, 1583–1588.
Kume, K., Zylka, M.J., Sriram, S., Shearman, L.P., Weaver, D.R., Jin, X., Maywood, E.S., Hastings, M.H., and Reppert, S.M. (1999). mCRY1 and mCRY2 areessential components of the negative limb of the circadian clock feedback loop. Cell 98, 193–205.
Lee, C., Etchegaray, J.P., Cagampang, F.R., Loudon, A.S., and Reppert, S.M. (2001). Posttranslational mechanisms regulate the mammalian circadian clock. Cell
107, 855–867.
Lee, H., Chen, R., Lee, Y., Yoo, S., and Lee, C. (2009). Essential roles of CKIdelta and CKIepsilon in the mammalian circadian clock. Proc. Natl. Acad. Sci. USA
106, 21359–21364.
Lee, K., Loros, J.J., and Dunlap, J.C. (2000). Interconnected feedback loops in the Neurospora circadian system. Science 289, 107–110.
Liu, P., Kevrekidis, I.G., and Shvartsman, S.Y. (2011). Substrate-dependent control of ERK phosphorylation can lead to oscillations. Biophys. J. 101, 2572–2581.
Manrai, A.K., and Gunawardena, J. (2008). The geometry of multisite phosphorylation. Biophys. J. 95, 5533–5543.

S10 Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors

Markevich, N.I., Hoek, J.B., and Kholodenko, B.N. (2004). Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades. J.
Cell Biol. 164, 353–359.
Michaelis, L., and Menten, M.L. (1913). Die Kinetik der Invertinwirkung. Biochem. Z. 49, 333–369.
Ogawa, Y., Koike, N., Kurosawa, G., Soga, T., Tomita, M., and Tei, H. (2011). Positive autoregulation delays the expression phase of mammalian clock gene Per2.
PLoS ONE 6, e18663.
Pittendrigh, C.S. (1954). On temperature independence in the clock system controlling emergence time in Drosophila. Proc. Natl. Acad. Sci. USA 40, 1018–1029.
Press, W.H. (2007). Numerical Recipes: The Art of Scientific Computing, Third Edition (Cambridge, UK/New York: Cambridge University Press).
Qiao, L., Nachbar, R.B., Kevrekidis, I.G., and Shvartsman, S.Y. (2007). Bistability and oscillations in the Huang-Ferrell model of MAPK signaling. PLoS Comput.
Biol. 3, 1819–1826.
Reppert, S.M., and Weaver, D.R. (2002). Coordination of circadian timing in mammals. Nature 418, 935–941.
Rust, M.J., Markson, J.S., Lane, W.S., Fisher, D.S., and O’Shea, E.K. (2007). Ordered phosphorylation governs oscillation of a three-protein circadian clock.
Science 318, 809–812.
Sangoram, A.M., Saez, L., Antoch, M.P., Gekakis, N., Staknis, D., Whiteley, A., Fruechte, E.M., Vitaterna, M.H., Shimomura, K., King, D.P., et al. (1998). Mammalian circadian autoregulatory loop: a timeless ortholog and mPer1 interact and negatively regulate CLOCK-BMAL1-induced transcription. Neuron 21, 1101–1113.
Shankaran, H., Ippolito, D.L., Chrisler, W.B., Resat, H., Bollinger, N., Opresko, L.K., and Wiley, H.S. (2009). Rapid and sustained nuclear-cytoplasmic ERK oscillations induced by epidermal growth factor. Mol. Syst. Biol. 5, 332.
Shanware, N.P., Hutchinson, J.A., Kim, S.H., Zhan, L., Bowler, M.J., and Tibbetts, R.S. (2011). Casein kinase 1-dependent phosphorylation of familial advanced
sleep phase syndrome-associated residues controls PERIOD 2 stability. J. Biol. Chem. 286, 12766–12774.
Shearman, L.P., Sriram, S., Weaver, D.R., Maywood, E.S., Chaves, I., Zheng, B., Kume, K., Lee, C.C., van der Horst, G.T., Hastings, M.H., and Reppert, S.M.
(2000). Interacting molecular loops in the mammalian circadian clock. Science 288, 1013–1019.
Strogatz, S.H. (1994). Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering (Reading, MA: Addison-Wesley).
Takahashi, J.S., Hong, H.K., Ko, C.H., and McDearmon, E.L. (2008). The genetics of mammalian circadian order and disorder: implications for physiology anddisease. Nat. Rev. Genet. 9, 764–775.
Tang, C.T., Li, S., Long, C., Cha, J., Huang, G., Li, L., Chen, S., and Liu, Y. (2009). Setting the pace of the Neurospora circadian clock by multiple independent FRQ
phosphorylation events. Proc. Natl. Acad. Sci. USA 106, 10722–10727.
Thomson, M., and Gunawardena, J. (2009). Unlimited multistability in multisite phosphorylation systems. Nature 460, 274–277.
van Zon, J.S., Lubensky, D.K., Altena, P.R., and ten Wolde, P.R. (2007). An allosteric model of circadian KaiC phosphorylation. Proc. Natl. Acad. Sci. USA 104,
7420–7425.
Vanselow, K., Vanselow, J.T., Westermark, P.O., Reischl, S., Maier, B., Korte, T., Herrmann, A., Herzel, H., Schlosser, A., and Kramer, A. (2006). Differential
effects of PER2 phosphorylation: molecular basis for the human familial advanced sleep phase syndrome (FASPS). Genes Dev. 20, 2660–2672.
Young, M.W., and Kay, S.A. (2001). Time zones: a comparative genetics of circadian clocks. Nat. Rev. Genet. 2, 702–715.

Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors S11

A

E
Km

m

1

,k

K

,k

E

K

m

3

1

7

,k

3

P
K

m

8

S11

P

kb1

S10

P
S01

ES00
kub1

k2

P

4

+
S00

,k
4
m

2

P

K

,k

2

,k 6

6

Km

Km

,k

8

S00

+

k1

7

Km

5

5

,k

S01

B

P

+

S10

E

Km1 = Km2
Km7 = Km8
Figure S1. PTO Model with Shared Intermediate Complexes at Branching Reactions, Related to Figure 1(A) PTO reaction network.
(B) In the model described in the main text, it was assumed that branching reactions (e.g., the conversion of S00 to either S01 or S10) involved two
distinct intermediates (e.g., ES00a for conversion to S01 and ES00b for conversion to S10). In the alternative approach shown here, branching reactions canproceed from a single complex ES00, leading to either S01 or S10. Results from calculations with this assumption are qualitatively similar to the model shown in
Figures 1B and 1C.

S12 Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors

A

C
0.25

0.8
0.6

OthersCluster 2

0.4

Cluster 1

Fraction of variance

Fraction

1

0.2

0.15
0.1
0.05

11
11
.5

10
10
.5

9
9.
5

5
8.

7.

5

8

0

0.2

0

Clustering diameter

0

5

10

15

20

Principal component

B

D

Others

0.025

0
1

10
100
k (min-1)

1000

Probability

0.03

00.01 0.1

k1
k2
k3
k4
k5
k6
k7
k8

1
10
Km ( M)

4
3
2
1

-6

-5

-4

-3

-2

-1

0
-1

0

1

2

3

4

-2
-3
-4
-5
PC 1

Km1
Km2
Km3
Km4
Km5
Km6
Km7
Km8

0.015

5

Cluster 1
Cluster 2
OthersPC 2

Probability

0.05

100 1000

Figure S2. Clustering of Oscillatory Parameter Sets, Related to Figure 2(A) Results of QT clustering with various clustering diameters.
(B) Parameter histograms obtained when all solutions not assigned to clusters 1 and 2 were pooled together.
(C) PCA of the random parameter search results. The first few principal components capture a reasonably large fraction of the variance.
(D) Projection of 10,000 parameter sets onto their first two principal components. Points are color-coded according to the clusters assigned using the QT
algorithm. PCA is able to identify cluster 1 fairly clearly, and a clear separation exists between cluster 1 and the remaining parameter sets. No other majorgroupings are apparent.

Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors S13

Probability

0.16

0.08

0
1

10

100
k (min-1)

1000

Probability

0.6

0.3

0
0.01 0.1

1 10 100 1000
Km ( M)

k1
k2
k3
k4
k5
k6
k7
k8

Km1
Km2
Km3
Km4
Km5
Km6
Km7
Km8

P
S01

S00

Biased

P

P
S11

P E
S10
low Km4

F P
S01

P

S00

P
S11

P
S10

E

Figure S3. Parameter Histogram for Oscillatory Parameter Sets Found when Clockwise Bias and Low Km4 Were Enforced, Related to Figure 3
The parameter histograms use the same coloring scheme as Figure 2C in the main text, and show clear similarities to the parameter histograms for cluster 1.

S14 Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors

A
8
6
4
2
0

10

10

1

2

10
k1 (min-1)

10

3

10

4

4
2
10-1
101
k2 (min-1)

6
4
2
-2

10

0

10
10
Km1 ( M)

2

8
6
4
2

4

10

10

10

6
4

10-3

103

S00 (x103 M)

8

8

2

10
S00 (x103 M)

S00 (x103 M)

6

10-3

10

-2

0

10
10
Km2 ( M)

2

10-1
101
k4 (min-1)

103

8
6
4
2

4

10

10-2

100
102
Km3 ( M)

104

10-1

101
103
Period (min)

105

0.8

8

S00 / Stot

S00 (x103 M)

8

S00 (x103 M)

S00 (x103 M)

S00 (x103 M)

10

6
4

0.6
0.4
0.2

2
10-2

10-1

100
101
Etot ( M)

101

102

103
Stot ( M)105

B
Cluster 1

Cluster 2

Km8

Km8

Km7

Km7

Km6

Km6

Km5

Km5

Km4

Km4
Km3

Km3
Km2

Km2

Km1

Km1

k8
k7

k8
k7

k6
k5

k6
k5

k4

k4

k3
k2

k3
k2

k1
0.1

k1
1

10

100

1000

10000

0.1

1

10

100

1000

10000

Figure S4. Effects of Parameter Variation on Oscillatory Behavior, Related to Figure 4(A) Single-parameter bifurcation diagrams. Solid black lines denote stable fixed points, and dashed black lines indicate an unstable fixed point surrounded bya limit cycle. Red vertical lines denote the parameter value in the symmetric example set; see Table S1 for precise values.
(B) Oscillation ranges for parameter sets in clusters 1 and 2. Bifurcations were sought between 103 and 104 for 1,000 parameter sets in each cluster; thecolored bars represent the average ranges over which each parameter can be varied without abolishing oscillations.

Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors S15

B

A

Number of molecules (x 105)

Cluster 2

Mean period (hour)

36

24

12

1

10

100

0

0

48
96
Time (hour)

1.0

1000

Number of enzyme molecule

0.006

0

0

48
96
Time (hour)

0.004
Molecule ID

SD / Mean period

2

Probability

0

Cluster 2
S00 S01 S11 S10

0.002
0
10
100
1000
Number of enzyme molecule

Fraction of robust oscillators

C

D

0.6
0.4

100
0

Cluster 1
Cluster 2

0.2

48Time (hour)96

0.3
Probability

11

Cluster 1
Cluster 2

0.2
0.1
0
0.5

0
1
5
10
20
Number of enzyme molecules

1
1.5
2
2.5
n=1 Period / n=5 Period

Probability

0.6

3Cluster 1
Cluster 2

0.4
0.2
0

n=5 Period / n=10 Period

Figure S5. Stochastic Simulation of a Stereotypical Cluster 2 Oscillator, Related to Figure 5(A) Effects of system size on the oscillation period. Oscillations were not observed for <20 enzyme molecules. Once they appear, the oscillations are quite robust
and the period rapidly converges on the bulk average value.
(B) Top row: Plots of the bulk population obtained in stochastic stimulations for the cluster 2 oscillator. Middle row: Smoothed traces showing the probability thata single (arbitrarily chosen) substrate molecule will occupy different phosphorylation states as a function of time. Bottom row: Results for 100 individually tracked
molecules, each shown as a horizontal raster line, color-coded as in the plots above.
(C) Fraction of parameter sets in clusters 1 and 2 that remain oscillatory as the system size is decreased. Even if only one of each type of enzyme is present,roughly 10% of parameter sets in each cluster remain robust against stochastic noise.
(D) Systematic changes in oscillator period with stochastic noise. As the system size is decreased (i.e., stochastic noise increases), cluster 1 oscillators tend toslow down (see also Figure 5A), whereas cluster 2 oscillators tend to speed up.

S16 Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors

C
Q

Probability

0.1

0.4
0.2

0

0.03

0.02

0.02

Probability

Probability

0.01

00.01 0.1 1 10 100 1000
Km2 (µM)

1000

0.01

0.01

0

Probability

0.03
0.015
0

0.03

0.06

0.015
0

1

10
100
k4 (min-1)

1000

00.01 0.1 1 10 100 1000
Km3 (µM)

10
100 1000
k6 (min-1)

1

10
100 1000
k3 (min-1)

Probability

1

1

0.03

G1

G2

G3

10
Q10(cycle)

100

G4

0.03

0.02
0.01
0
0.01 0.1 1 10 100 1000
Km6 (µM)

00.01 0.1 1 10 100 1000
Km4 (µM)

10
100 1000
k5 (min-1)

100.01 0.1 1 10 100 1000
Km7 (µM)
Probability

10100
k7 (min-1)

1

Probability

1000

0.02

0.5

0.06
Probability

Probability

Probability

Probability
Probability

10100
k2 (min-1)

0

Probability

0.02
0

1

Q10(E) : Q10(F)
1:3
1.5:3
2:3
2.5:3
3:3

0
0.1

100

0.06

0.04

0

00

10

00
0.

D
0.02

1
10
Q10(cycle)

1.0

Q10(cycle)

Q10(E) : Q10(F)

0.04

0
0.1

0

01
0.
00
1

1:
1.
5
1:
1: 2
2.
5
1:
1. 3
5:
3
2:
2. 3
5:
3

0

1:1.5
1:1.5
1:2
1:2.5
1:3

10

20

0.6

No R2 cutoff

Q10(E) : Q10(F)

0.2

1

40

No R2 cutoff

=1:3

10

60

(F)
10

0.8

0.
01

R2 of power-law ﬁt

1

80

:Q

0.
1

100

(E)
10

Probability

B% of Q10(cycle) calculations
with power-law fit R2>0.9

A

0.06
0.03
0
0.01 0.1 1 10 100 1000
Km5 (µM)

G5

E
Probability

G1

G2

0.06
0.03
0

G3

0.06

0.06

0.03

0.03
0

0
1

10 100k (min-1)

1000

1

10 100k (min-1)

G1
Probability

0.04
00.01 0.1 1 10 100 1000
Km (µM)

G5
0.06

0.03

0.03
0

0
1

1000
k1
k2

0.08

G4
0.06

10 100k (min-1)

1

k5
k6

k3
k4

G2

1000

10 100k (min-1)

1000

1

10 100k (min-1)

1000

k7
k8

G3

G4

G5

0.06
0.06
0.03

0.03
0
0.01 0.1 1 10 100 1000
Km (µM)

Km1
Km2

00.01 0.1 1 10 100 1000
Km (µM)
Km3
Km4

Km5
Km6

0.06

0.06

0.03

0.03

0
0.01 0.1 1 10 100 1000
Km (µM)00.01 0.1 1 10 100 1000
Km (µM)

Km7
Km8

Figure S6. Dependence of Temperature Sensitivity on Parameter Distributions, Related to Figure 6(A and B) Quality of the exponential fit for the calculated period at varying temperatures. Most fits are of excellent quality, regardless of the values of Q10(E) and Q10 (F)(A), whereas poorly fit parameter sets tend to show a cycle Q10 close to 1.0 (B).
(C) Histogram of cycle Q10 values when the poorly fit parameter sets (R2 < 0.9) are included.
(D) Changes in parameter distributions as a function of temperature overcompensation or undercompensation. Histogram traces are color-coded as in Figure 6C.
(E) The same histograms as in (D), organized by group rather than by parameter.

Cell Reports 2, 938–950, October 25, 2012 ª2012 The Authors S17