1462

Biophysical Journal

Volume 107

September 2014

1462–1473

ArticleA Mammalian Circadian Clock Model Incorporating Daytime Expression
Elements
Craig C. Jolley,1 Maki Ukai-Tadenuma,1 Dimitri Perrin,1 and Hiroki R. Ueda1,2,3,4,5,*
1

Laboratory for Systems Biology, RIKEN Center for Developmental Biology, 2-2-3 Minatojima-minamimachi, Chuo-ku, Kobe, Hyogo 650-0047,
Japan; 2Laboratory for Synthetic Biology, RIKEN Quantitative Biology Center, 2-2-3 Minatojima-minamimachi, Chuo-ku, Kobe, Hyogo
650-0047, Japan; 3Department of Biological Sciences, Graduate School of Science, Osaka University, 1-1 Machikaneyama, Toyonaka, Osaka
560-0043, Japan; 4Graduate School of Frontier Biosciences, Osaka University, 1-3 Yamadaoka, Suita, Osaka 565-0871, Japan; and5
Department of Systems Pharmacology, Graduate School of Medicine, University of Tokyo, Hongo 7-3-1, Bunkyo-ku, Tokyo, 113-0033,
Japan

ABSTRACT Models of the mammalian clock have traditionally been based around two feedback loops—the self-repression of
Per/Cry by interfering with activation by BMAL/CLOCK, and the repression of Bmal/Clock by the REV-ERB proteins. Recent
experimental evidence suggests that the D-box, a transcription factor binding site associated with daytime expression, playsa larger role in clock function than has previously been understood. We present a simplified clock model that highlights therole of the D-box and illustrate an approach for finding maximum-entropy ensembles of model parameters, given experimentally
imposed constraints. Parameter variability can be mitigated using prior probability distributions derived from genome-wide
studies of cellular kinetics. Our model reproduces predictions concerning the dual regulation of Cry1 by the D-box and RevErbA/ROR response element (RRE) promoter elements and allows for ensemble-based predictions of phase response curves(PRCs). Nonphotic signals such as Neuropeptide Y (NPY) may act by promoting Cry1 expression, whereas photic signals likely
act by stimulating expression from the E/E’ box. Ensemble generation with parameter probability restraints reveals more about amodel’s behavior than a single optimal parameter set.

INTRODUCTION
Systems biology aims to develop predictive mathematicalmodels of cellular phenomena. Models can be abstracted
into two components—the model structure, which describes
the relationship between modeled quantities and specifies
rules for the dynamic updating of these quantities, and themodel parameters, which contain problem-specific quantitative information. Although the model structure is often proposed based on expert knowledge of the components of asystem and their relationships, identification of model parameters requires quantitative data. Physical scientists prefermodels with few parameters, epitomized by von Neumann’s
famous jibe about fitting an elephant with four parameters(1). Biology involves complex, heterogeneous systems—
the components of a cell are far more diverse than are theconstituents of a crystal or an atomic nucleus—and quantitative parameter identification is unavoidable.
In what has been termed sloppy behavior (2–4), a model’s
output is strongly affected by certain linear combinations ofparameters, whereas other combinations have little impact.
When a sloppy model is tuned to agree with system-level
data, the uncertainties for the estimates of individual parameters are often disturbingly large. It is only when parameters
are viewed in a systemic context that the experimental data

Submitted January 2, 2014, and accepted for publication July 2, 2014.
*Correspondence: uedah-tky@umin.ac.jp

show their true (albeit limited) ability to constrain themodel.
Given an estimated set of parameters, uncertainty quantification typically involves the calculation of confidence
intervals—finite ranges within which the true parametervalues can be localized to within a certain probability. In
many cases, the confidence interval is infinite and parameters are said to be nonidentifiable. It is often useful to distinguish between practical nonidentifiabilities that result frominsufficient or noisy data and structural nonidentifiabilities
that arise when there is an insufficient mapping of model
states to observables (5).
Confidence intervals can be estimated using the Fisher Information Matrix, which requires linearization about the
optimized solution (6). When dealing with highly nonlinear
models, it is often preferable to use a bootstrap method, inwhich parameters are repeatedly reestimated with newexperimental data to generate an ensemble of parameter
sets (7). Because new data are not always readily available,
simulated data can be generated by Monte Carlo techniques,
if the variance of the experimental data can be estimated.
Bootstrap methods are difficult to apply to models withcomplex dynamics such as oscillations because parameters
often cannot be efficiently reestimated (8). Bootstrap
methods allow the estimation of confidence intervals not
only for parameters, but also for systems-level outputs ofthe model.

Editor: Reka Albert 2014 by the Biophysical Society
0006-3495/14/09/1462/12 $2.00

http://dx.doi.org/10.1016/j.bpj.2014.07.022

Clock Model with Daytime Elements

In this study, we use Metropolis Monte Carlo (MMC)
sampling (9) to achieve the same goals as bootstrap
methods, without the need for repeated reoptimization.
The objective function for model optimization was used as
an energy function and a maximum-entropy ensemble of parameters with a noise scale similar to the experimental data
was generated. This allows us to set bounds on our parameter estimates and model predictions, and determine whichof those predictions are robust to changes in the parametervalues that are within the experimental noise.
Circadian clocks have long been a fruitful system formathematical modeling, and the mammalian clock in particular has been the subject of numerous modeling studies(10–14). Most mammalian clock models have focused on
two feedback loops. In the first, a CLOCK/BMAL1 heterodimer binds to the E-box region upstream from the Per and
Cry genes, activating their transcription. A PER/CRY heterodimer binds to the CLOCK/BMAL1 complex, abolishingits activity and indirectly repressing the expression of Per
and Cry. In the second loop, the transcription of Clockand Bmal1 are repressed by the binding of REV-ERBa
to an RRE region in their promoters, whereas RevErba
is controlled by an E-box (and therefore by CLOCK/BMAL1 and PER/CRY). Some researchers (12) havereferred to this as a positive feedback loop, emphasizingthat Bmal1 and Clock indirectly activate their own transcription by activating Per and Cry, which repress RevErba andtherefore derepress the RRE-controlled Clock and Bmal1
genes. Others (11) have modeled Bmal1 as self-inhibitory
because BMAL1 is also a direct activator of RevErba.
RevErba/ knockout mice show only a modest circadianphenotype, and most modeling studies assume that theloop involving RevErba serves only a stabilizing or regulatory role and that oscillations should persist in its absence.
Recently, RevErba/b double knockout mice have beenused to show that the (partially redundant) RevErb genes
are essential; double-knockout mice show phenotypes
similar to knockouts of core clock genes such as Per andBmal1 (15).
Recently, a different organizational scheme has been proposed (Fig. 1 A) (16–18). This approach focuses on threeclock-controlled elements (CCE), transcription factor (TF)
binding sites that give their target gene a circadian expression schedule. These three are the E/E’-box, which governs
morning expression (19–22); the D-box, which promotes
daytime expression (21,23); and the RRE, which leads toevening expression (21,24–26).
When an exogenous luciferase reporter is competitively
regulated by an E/E’-box-driven promoter and an RREdriven repressor, the result is a daytime expression phasesimilar to that obtained when the reporter is driven by aD-box. In other words, combining morning activation andnighttime repression generates a D-box. Similarly, a nighttime phase characteristic of RRE expression can be obtainedby combining a daytime promoter and a morning repressor

1463

E4BP4

DBP
DB
BP

A

TEF
HLF
BMAL2

BMAL1

CLOCK

NPAS2

E/E'box

CRY2

D-box

CRY1
Rev-Erbα
DEC1
Rev-Erbβ
β
DEC2
RORα
PER1
RORβ
CK1ε

PER2
CK1δ

RORγ
RO
PER3
PER
R3

B

RRE
RE
E

C
A

E/E'box

D-box

2

A

B

1

C
D

5

4

B

7

6
RRE

3

9

C
8

E
D

10

E

FIGURE 1 (A) The core clockwork involves transcriptional activators(green), repressors (magenta), and protein kinases (yellow). Transcriptional
regulation involves three clock-controlled elements: the E/E’-box, theD-box, and the RRE. Gray lines indicate that a transcription factor (TF)
is regulated by a particular CCE, magenta or green lines indicate theCCE to which a given TF binds. Dotted lines denote computationally inferred relationships. (B) Simplified version, motivated by synthetic-biology
experiments. (C) Realization of the simplified clock architecture using two
activators and three repressors. Numbered boxes show regulatory interactions, color-coded by CCE. To see this figure in color, go online.

(17). Cry1 contains an E/E’-box and D-boxes in its promoter, as well as an RRE in an intron (18). When exogenous
Cry1 was driven by an engineered promoter in cultured
Cry1/:Cry2/ cells, a combination of D-box and RRE
elements recapitulated the wild-type behavior, whereas
the (widely assumed) E/E’-box/RRE combination led toaberrant circadian rhythms. CRY1 is an E/E’-box repressor
with an evening expression phase; reconstructing the E/E’box expression using D-box and RRE elements completes
the simplified network shown in Fig. 1 B and C.
This topology differs from the canonical clock structure
in several ways. Instead of a deconstructive approach inwhich system components are identified and their interactions characterized, it is informed by a synthetic approachin which interactions are added until the desired systembehavior can be replicated. As a rule, clock modelingstudies have neglected the D-box. Although the essentialityof the D-box is still unclear (27), its role in Cry1 regulation
argues for its importance.
Our study should be placed in the context of previous
clock models (10–14). The clearest difference is our inclusion of genes regulated by the D-box; previous studies
Biophysical Journal 107(6) 1462–1473

1464

Jolley et al.

have focused primarily on the Per, Cry, Bmal, and RevErb
genes, with the possible inclusion of Clock (11,13,14),
Rorg (13), Npas (14), and various circadian kinases (14).
The model by Mirsky et al. uses Rorg to capture the effectsof all of the Ror activators; this accounts for the combinatorial regulation of Rorg by the E/E’-box and the RRE, but
not the regulation of Rora and Rorb by the D-box (21).
Although recent experiments suggest dual regulation of
Cry1 by RRE and D-box elements, previous models havedescribed Cry1 as regulated either by an E-box alone(11,12) or by an E-box together with an RRE (10,13,14).
Although the precise role of the D-box is likely to remain
controversial, inclusion in modeling studies should be apart of future assessments of its function.
Another important difference is in computational methodology. Early modeling studies (10–12) used manual,
trial-and-error parameter searches. Later studies used
global search strategies such as evolutionary optimization(13) or simulated annealing (14); our work follows thismore recent trend by using differential evolution (28) for
optimization. Previous studies assessed the model’s robustness primarily by modifying individual parameters; ourmethod provides a more comprehensive picture of parameter variability by using MMC sampling, and is somewhat
in the spirit of clock modeling studies that have used bootstrap sampling (8).
Finally, previous clock models have been validated primarily based on their ability to predict knockout mutants,
including some subtle and counterintuitive effects. Doing
this successfully often requires the inclusion of multiple isoforms of the Cry, Per, and Bmal genes, to capture their partial functional redundancy (13,14). Our study has focusedinstead on a minimal model inspired by synthetic-biology
approaches; redundant isoforms are omitted and many of
the relevant knockouts have not yet been characterized.
Instead, we use our model to make predictions regarding
the phase response curves obtained by perturbations ofdifferent genes, and outline how these phase response predictions might be experimentally accessible.
MATERIALS AND METHODS
Model equations
The gene regulatory network can be represented as a system of ordinary
differential equations (ODEs). The mRNA variables a to e and the protein
variables A to E correspond to circadian clock genes as listed in Table 1.
TABLE 1 Possible assignments of model genes tomammalian clock genes
Model gene
a
b
c
d
e

Biophysical Journal 107(6) 1462–1473

Mammalian clock gene(s)Dbp/Tef/Hlf
Rora/b
Cry1
Rev-Erba/b
E4bp4

da
¼ h1 ðCÞ  hma adt

dA¼ a  hpA A
dt

db¼ g2;3 ðA; EÞ  hmb bdt

dB¼ b  hpB B
dtdc¼ g4;7 ðA; EÞ þ g5;6 ðB; DÞ  hmc c
dt

dC¼ c  hpC C
dt

dd¼ h8 ðCÞ  hmd d
dt

dD¼ d  hpD Ddt

de¼ g9;10 ðB; DÞ  hme edt

dE¼ e  hpE E
dt(1)

Details of the model equation derivation can be found in the Supporting
Material. The dependent variables a to e and A to E are unitless; the absolute concentration scales are arbitrary and we are only interested in relativechanges. The h parameters are unitless degradation rates of mRNAs andproteins. The functions h(R) and g(A,R) are used to capture the effects ofgene regulation; both are bounded to the interval [0,1], with h(R) monotonically decreasing (repression), and g(A,R) increasing in A but decreasing inR (activation and repression on the same promoter). Binding follows theHill equation; activators and repressors bind competitively to the same promoter as in the following:Hi ðRÞ ¼

1
1 þ ðci RÞn

gi;j ðA; RÞ ¼

ðci AÞn nn
1 þ ðci AÞ þ cj R
(2)Here, the ci’s are dimensionless activation constants—low ci values correspond either to weak binding or low-abundance TFs. To limit the number offree parameters, the Hill coefficient n was assumed to be equal for all
reactions.
If all dimensionless degradation rates (the h’s in Eq. 1) are multiplied bya factor gt and all of the activation constants c are multiplied by gt2, thenthe period will decrease by 1/gt, whereas the model output is otherwise
unchanged.

Knockout phenotype constraints
The circadian clock is rather resilient against gene knockouts (13,15,27),
especially at the organism level, but our simplified model lacks redundancy
and is more easily rendered arrhythmic. Rather than modifying Eq. 1, we
can obtain knockout effects by setting some of the ci variables to zero.
For example, a Dbp/Tef/Hlf triple knockout could be approximated by
setting c2 ¼ c4 ¼ 0; levels of the protein variable A will still be nonzero,but it will not have any effect on the other variables. In reality, mice withthis triple mutation are rhythmic (albeit prone to epilepsy). The simulated
knockout is always arrhythmic, because b and c are always zero in theabsence of A, and no closed feedback loops exist. Other knockout phenotypes disturb clock function in ways that significantly constrain our model.
The Cry1/2 (29) and Rev-Erba/b double knockouts (15) show complete arrhythmicity. Fig. 1 C shows that the removal of Cry1 disrupts all feedback
loops and makes oscillations impossible.
Removal of Rev-Erba/b leaves some feedback loops intact, and parameter sets with c6 ¼ c10 ¼ 0 are sometimes rhythmic. In another interesting
set of mutants, either the D-box or the RRE has been removed from Cry1(18); removal of the RRE leads to arrhythmicity, whereas D-box removal
lengthens the oscillation period. Correct phenotypes for these mutations
are not guaranteed by the model topology.

Clock Model with Daytime Elements

1465

To ensure correct knockout phenotypes, our model was reintegrated threetimes: once with c6 ¼ c10 ¼ 0 (Rev-Erb knockout), once with c5 ¼c6 ¼
0 (Cry1DRRE), and once with c4 ¼ c7 ¼ 0 (Cry1DD-box). If any of these
did not produce the correct oscillation phenotype (arrhythmic, arrhythmic,
and long-period, respectively), the objective function (see below) was set to~10300 (the maximum double-precision floating point value), effectively
removing these parameter sets from consideration.

To estimate physically reasonable mRNA-protein phase lags, an mRNAm and a protein p were modeled using the following:

Parameter search

These equations can be solved analytically (see Supporting Material) toobtain an asymptotic protein oscillation amplitude of

Initially, parameters were chosen randomly from the probability distributions described in Table S1 and Fig. S1 in the Supporting Material. The distribution of mRNA decay constants was derived from the data set of
Sharova et al. (30), who measured the decay rates of 1825 TF mRNAs ina mouse cell culture. The rate distributions for mRNA transcription, as
well as protein translation and degradation, were drawn from Schwanhäusser et al. (31), who used parallel metabolic pulse labeling to estimate
these rates for more than 5000 genes in mammalian cells. The distribution
for TF-DNA binding constants was chosen to have significant density over
the range from 0.002 to 10 nM (32), but to have a density near zero forvalues above 1000 nM, a region dominated by nonspecific binding (33).
The nondimensional parameters c and h above were obtained from these
distributions as described in the Supporting Material. The pairs (c1, c8),
(c2, c4), (c3, c7), (c5, c9), and (c6, c10) were constrained to be equal in
the initial search, because they describe the binding constants for a single
TF on similar promoters. For each randomly chosen parameter set, thefixed point was identified using a multidimensional Newton search (6)and the stability determined from the eigenvalues of the JacobianJik ðtÞ ¼ vfi ðtÞ=vyj at the fixed point (34).
In this initial search, ~ 1.86% of parameter sets showed limit-cycle oscillations. When knockout constraints were included in the initial search, thesuccess rate decreased to 0.015%. Because the experimental parameter distributions were measured for a broad range of cellular genes, rather than specifically for oscillating transcription factors, one can expect to see somesystematic differences between the distribution used in the search and the
distribution of parameters in successful oscillators. Typical oscillation periods were around 200 h; parameters were rescaled as described above togive 24-h periods. In general, this led to degradation rates that were fasterin the oscillating parameter sets than for cellular genes and proteins generally. To compensate for this difference, we calculated new parameter distributions from the search results that were used in scoring (Fig. S1).

Estimating protein oscillation phases
Model genes were identified with real circadian genes (see Table 1) for
which qPCR data were available for mRNA abundances in the mouse
Suprachiasmatic nucleus (SCN) (21). The available protein data were farless quantitative; published results (22,35–38) were available only for proteins sampled from the mouse liver, typically at 4-h time intervals, and(roughly) quantified by Western blotting. Our approach was to estimate
amplitudes from published data and use the amplitude data to derive otherfeatures of the protein abundance curves.
Protein amplitudes were obtained either by calculating band densities
from published blot images using ImageJ (39), or (when band densities
were quantified in the original publication) converting published plots
into a data table format using Engauge (http://digitizer.sourceforge.net/).
For DBP and REV-ERBa, oscillations were so large that a reliable cosine
fit was impossible; these were assigned a conservative oscillation amplitude of 0.8. The amplitudes obtained by curve fitting for CRY1 (¼ 0.521)
and E4BP4 (¼ 0.461) were usable. No oscillation amplitude could beobtained for (ROR)a because, in contrast to the SCN, it does not oscillate in the mouse liver. A conservative value of 0.25 was assigned, equal
to the smallest oscillation amplitude found for any circadian protein(CLOCK).



dm
1
2pt
¼
cos
þ 1 bm  am m
dt
2
24

(3)

dp
¼ bp m  ap p
dt

am ap
rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
:



2
2
2
2
a p þ u 0 u0 þ a m(4)

The phase lag between the mRNA and protein oscillations is obtained using
the following:

tan1



ap sin 4 þ u0 cos 4
4
ap cos 4  u0 sin 44 ¼ tan1

u0
am

(5)

Note that both quantities depend only on the degradation rates am and ap,
and the driving frequency u0 ¼ 2p/24. When values for am and ap weredrawn at random from the experimentally derived distributions in Table
S1, a roughly linear relationship was found between phase lag and protein
amplitude (Fig. S2). This allows us to estimate physically reasonable phaselags based on protein amplitudes. For DBP, CRY1, REV-ERBa, andE4BP4, the resulting phase lags were consistent with the available Western
blot data. The higher-amplitude proteins (DBP and REV-ERBa) were assigned the shortest delays (1.7 h) whereas E4BP4 was given a longer delay(3.4 h). For protein RORa, Western blot data was unavailable and the phaselag was set to 4.5 h, consistent with the relatively low amplitude assumed
above. These phase lags were then added to the (more accurate) mRNA
phases to obtain the target phases depicted in Fig. S3.

Definition of the objective function
Parameter sets were scored by integrating the dynamical equations until they
converged to a limit cycle, then measuring the total squared deviation between
the experimentally derived data points and the model abundances (Figs. 2,
Fig. S4). For the protein components, cosine curves were calculated using
the amplitudes and phases described above, with the same sampling density(12 points over 48 h) as the qPCR data. The published qPCR data (and theestimated protein data) had been normalized to have a mean of 1.0. Because
the model equations are nondimensional, the absolute oscillation baseline isarbitrary and model outputs were similarly normalized before comparison.
The score function should be insensitive to initial conditions, so the score isminimized relative to an arbitrary phase shift 4 as in the following:X

S ¼ minfX!

x˛fa;b;c;d;e;A;B;C;D;Eg

 modelðx; ti þ fÞÞðexpðx; ti Þ

i

(6)

2

In addition to measuring the quality of the model output, one can assess theplausibility of the model parameters themselves. The previous probabilityof each individual parameter value pj was calculated using the probability
distributions Pj obtained from the initial parameter search (see Table S2).

J ¼ ln

Y
j

 
Pj pj

!

þ ln

Y

 
max Pj!

(7)

jBiophysical Journal 107(6) 1462–1473

1466

Jolley et al.

1RUP FRQF



'ES D



no P



P



Data

















5RUD E















&U\ F




















5HY(UEα G









(ES H







FIGURE 2 Ensemble-based predictions ofmRNA levels. Parameters were chosen by minimizing the total squared deviation from the data
points shown, either without (‘‘no P’’) or with(‘‘P’’) probability constraints on parameter values.
Bold lines show the output for the optimized
parameter values; shaded areas contain 80% of
the traces from Monte Carlo ensembles. Similarplots for the protein components can be found in
Fig. S4. To see this figure in color, go online.




































7LPH KIn some of the calculations described below, this term was added to thescore function defined above.

SP ¼ S þ gJ

(8)Generation of phase response curves
Phase response curves (PRCs) were generated following the method of
Kramer et al. (41). Briefly, a general set of coupled ODEs with a parameter
vector a can be written as in the following:

The weighting term g controls the relative balance between the two components. The value chosen (0.1) was sufficient to constrain the distribution ofmodel parameters without severely compromising the quality of the fit(Fig. 2).

Parameter set optimization
Parameters were optimized using a differential evolution (DE) algorithm(28). 104 oscillating parameter sets were randomly generated from the
distributions in Table S1 and rescaled to have 24-h periods; the 256 highest-scoring sets were used as the starting population. DE converged on anear-optimal parameter set; the result was further refined by simplex minimization (40).

Generation of model ensembles
Once a global optimum was identified, an ensemble of nearby points could
be generated using the MMC algorithm (9). The fitted qPCR data comprised
five transcripts, sampled every four hours over a 48-h period; an error scale
was established by calculating the median absolute deviation between themodel results and these 60 data points. An MMC sampling temperature
was chosen that gave roughly the same median absolute deviation between
the optimized model and the members of the ensemble. The scoring
schemes with fewer constraints (i.e., omitting the parameter-probability
and knockout terms) were able to fit the experimental data more closely
and therefore required lower-temperature MMC sampling to generate an
ensemble with the same median error.
Biophysical Journal 107(6) 1462–1473

dy¼ fðy; a; tÞdt(9)We consider a to include not only the parameters that enter explicitly into f
but also the initial conditions y(0). Using the Jacobian, we can calculate theadjoint Green’s function matrix as in the following:d y 0
K ðt ; tÞ þ Ky ðt0 ; tÞJðt0 Þ ¼ 0;
dt0

Ky ðt; tÞ ¼ I; t0 %t(10)Once the system has converged to the limit cycle, its state can be completely
described by specifying its phase along the limit cycle trajectory. A differential perturbation vyk to the state variable yk, delivered at time t’, will leadto a differential delay or advance vt that can be calculated as in thefollowing:

Qk ðt0 Þ ¼

vt¼ lim Kjky ðt0 ; tÞt/Nvyk ðt0 Þ

dyj ðtÞ
dt(11)

In practice, this calculation requires integrating Eq. 9 (in our case, Eq. 1)
forward from a defined starting point on the limit cycle to time t and saving
values of J at all of the integration time points. The t/N limit is satisfied
as long as the forward integration has adequately converged to the limit
cycle. J must then be interpolated (e.g., using a cubic spline) to give values
at arbitrary time points for the integration of Eq. 9 backward to time t’. This
backward integration only needs to be performed once; one can save theClock Model with Daytime Elements

1467

values of K at each integration time point and interpolate between them toobtain Qk(t’) at arbitrary times.
It is also possible (albeit less precise) to calculate PRCs by a method
resembling an experimental protocol, in which state variables or model
parameters are perturbed at a defined point during the integration and theresulting phase shift measured. This method produces curves with different
magnitudes, but roughly the same shapes.

RESULTS
Ensemble parameter distributions
Fig. 3 shows the distribution of parameters obtained underfour different scoring conditions: both with and without
the probability constraints (‘‘P’’) and the knockout constraints (‘‘KO’’). For each set of conditions, 6400 parameter
sets were sampled from the MMC ensemble and the pooled
set of 25,600 parameter sets, along with the four optimized
solutions, was projected onto their first two principal components and visualized using minimum-area convex hulls(Fig. 3 A). The ensembles with probability constraints
occupy a smaller area than those without them, illustrating


A













3&





3&


B











3&






QR3
3



QR.2QR3
QR.23



3&

FIGURE 3 (A) Parameter values were drawn from a Monte Carloensemble and projected onto their first two principal coordinates. The solidlines show a minimum-area convex hull containing 90% of points. (B) Minimum-area convex hulls for parameter sets sampled under four different
scoring conditions: both with and without the parameter value probability
constraints (‘‘P’’ versus ‘‘no P’’) and either with or without (‘‘no KO’’)
the knockout constraints. Locations of optima are shown by round dots.
The probability constraint generally results in a smaller sampling region,and removal of the knockout constraints permits the exploration of otherwise forbidden regions of parameter space. To see this figure in color, go
online.

the constraints’ success in limiting access to unrealistic
parameter values. The calculation with neither probability
nor knockout constraints (‘‘no KO, no P’’) occupies a smallerarea than the version with knockout constraints only (‘‘noP’’); this is likely because the ‘‘no KO, no P’’ calculation
could provide a closer fit to the experimental data and wastherefore sampled at a lower temperature. The ‘‘no KO’’sampling also ventures into regions of parameter spacethat appear to be inaccessible when the knockout constraints
are present. Fig. S5 shows the degree to which knockout
constraints are violated when they are not enforced during
optimization and sampling. When the probability constraints
are present but the knockout constraints are not (‘‘no KO, P’’
in Fig. 3 B), the arrhythmic RevErb knockout phenotype isalmost never reproduced correctly. This appears to result
from specific changes in the parameter distribution.
Fig. 4 presents another perspective on parameter variability—here individual parameter values are comparedacross the same four sampling schemes. Again, the two ensembles using the probability constraint show less variationin parameter values than the two that omitted it.
One unanticipated effect of the probability constraint termis that strong regulatory interactions (i.e., large c values)
were compensated by weakening other interactions, often
to the point of insignificance. The net effect was a ‘‘pruning’’of the regulatory network that identified interactions that
may not be important for the correct functioning of themodel. For example, c3 and c7, the parameters controlling
E4BP4’s interactions with Rora and Cry1, were lower thanthe other c parameters. When these were set to zero (i.e.,an E4bp4 knockout), rhythmicity was nearly always retained. This suggests that E4bp4 may function primarily asa circadian output gene, rather than as a member of thecore clock network. Similar effects are observed for c9
and c10, the parameters controlling the regulation of E4bp4.
Similar pruning effects allow us to clearly see the effects
of knockout constraints. The clearest effect of the knockout
constraint (i.e., the difference between ‘‘P’’ and ‘‘no KO, P’’)is in the values of c4 and c6, the parameters that control regulation of Cry1 by DBP (via the D-box) and REV-ERBa (via
the RRE), both of which are repressed by CRY1 via the E/E’box. Fig. S5 suggests that the major difference between these
two ensembles is that the knockout-unconstrained ensemble
fails to reproduce the RevErb knockout phenotype. In otherwords, a correct RevErb knockout phenotype is correlatedwith the regulation of Cry1 by E/E’-box-controlled genes.
If other clock genes with as-yet-uncharacterized knockout
phenotypes (such as Rora or E4bp4) are shown to be essential, a similar approach will allow us to identify the regulatory interactions underlying their essentiality.
Experimental predictions
Ensemble-based parameterization allows us to quantify uncertainty in our predictions of experimental results. Because
Biophysical Journal 107(6) 1462–1473

1468

Jolley et al.
mRNA degradationProtein degradation












ηmb

ηmc

ηmd

ηpB

ηpC

ηpD

ηpE

ηme





10

ηpA



ηma



Hill coefficient (n)8Activation constants








6

4





2

0



χ1

χ2

χ3

χ4

χ5

χ6

χ7

χ8

χ9

χ10




FIGURE 4 Parameter distributions for the four different scoring scenarios; data series are named as in Fig. 3. The probability constraint term suppresses
overly fast mRNA turnover and keeps activation constants confined within a narrow range but has little effect on protein degradation rates or the Hillcoefficient. To see this figure in color, go online.

the model was parameterized with system-level data, it will
be most useful to make a prediction of system-level
behavior, determining how the clock as a complete system
will respond to a perturbation.
PRCs are used in chronobiology to describe the effect ofexternal signals (light, chemicals, etc.) on the phase of theclock (42). A PRC is obtained by plotting the phase shift
as a function of the preperturbation phase. For example, under constant-dark conditions, an early-morning light pulse
will shift the clock to an earlier time (dawn is perceived tohave come early), whereas an evening pulse produces a positive phase shift (dusk is perceived to have come late). Under constant-light conditions, a dark-pulse PRC can also beconstructed by briefly turning the lights off at defined times.
Although our model does not explicitly incorporate theeffects of external stimuli, phase-response curves can becalculated for differential perturbations to the state variables. In cultured SCN slices, phase shifts have been
observed as a result of chemical stimulation (43–46), but
the connection between the signal-response cascades andthe core clockwork is unclear, especially for the dark-pulse
response. For example, external stimulation could lead totranscription activation, targeted protein degradation, orother downstream effects. If the phase response of state variable perturbations is similar to that of experimental perturbations, one can begin to search for a causal link between
external perturbations and internal changes.
Biophysical Journal 107(6) 1462–1473

Fig. 5 shows the PRCs calculated for perturbations to thefive mRNA components of the model (see Fig. S6 for protein
components). The PRCs calculated for the parameterizations
with and without the probability constraints have similar
shapes, but sometimes wildly different magnitudes. The
shape of the PRC curves can be compared with ‘‘lighttype’’ PRCs such as N-Methyl-D-aspartate (NMDA) (47)or ‘‘dark-type’’ PRCs such as NPY (48). Such comparisonsshould be made cautiously, however; the experimental data
come from organotypic slice cultures in the case of NMDA
or intact mice in the case of NPY and may reflect higherlevel effects that are not observed in the single-cell clock.
The PRCs for Dbp and RevErba show some similarities
to NMDA. Light stimulation has been known to promote
Per1/2 transcription (49–51); all four of these genes haveexpression schedules that are controlled by an E/E’-box.
Stronger similarities can be seen between the dark-type
PRC for NPY and the calculated PRC for Cry1. The mechanism by which NPY affects the core clockwork is still unclear; Cry1 may be a fruitful target for future investigation.
The NMDA data of Asai et al. (47) used cultured SCN slices for which circadian time is more difficult to interpret thanin an intact rodent. They defined CT0 as the luminescencepeak of luciferase driven by a Per1 promoter. When SCN tissues are synchronized to an external zeitgeber, the luminescence peak coincides with the middle of the light period (52).
All of our experimental protein/mRNA phases (as well as theClock Model with Daytime Elements

1469

A related test of our model’s predictive power is suggested by the rescue of Cry1/:Cry2/ knockout cells using exogenous Cry1 driven by a synthetic promoter (18).
When the phase of Cry1 expression was altered by usingsynthetic promoters with varying numbers of D-box and
RRE elements, it was found that weakening of the D-box
regulation correlates with a later peak of CRY1 expression(relative to the period) and an overall lengthening of thecircadian period. We can mimic this effect by rescaling c4
and c7, which control the binding of DBP and E4BP4 tothe promoter region of Cry1. For the optimized parameter
set, a decrease in these two binding constants led to a diminished role for the D-box in the regulation of Cry1, leading inturn to a phase delay of CRY1 expression and an increase in
the oscillation period (Fig. 6), consistent with the experimental trend.
When eight parameter sets were randomly chosen from
the MMC ensemble, seven show the expected positive
trend. Period/phase curves were calculated for modified
versions of 6400 MC-sampled parameter sets. The resulting distribution of points is mostly in the positive quadrant(lengthened period and delayed CRY1 phase), suggestingan overall trend that agrees qualitatively with experiments.
It may be possible to generalize: systems-level outputs
become less consistent as they are less similar to the typeof behavior used to parameterize the model. Our score

NPY PRC) were drawn from intact rodents subjected to alight period lasting from CT0-12. The luminescence peakshould be at roughly CT6; the data points in Fig. 5 are shiftedby 6 h relative to the presentation in Asai et al.
The PRC for E4bp4 is more enigmatic. In the probabilityconstrained calculation, perturbations to E4bp4 levels havevirtually no effect on the phase. When the model is optimized without probability constraints, the effects of perturbations to E4bp4 levels are huge. As was mentioned above,
Fig. 4 suggests that E4bp4 may be a circadian output, rather
than an integral part of the core feedback loops. In this case,
we would expect direct perturbations to have no effect on
the core clockwork.
The PRC shapes can be interpreted using Fig. 2, whichshows mRNA concentrations as a function of circadiantime. The expression level of Cry1 peaks around lights-off(CT12) and has a minimum near lights-on (CT0); theextrema for CRY1 protein will be about 2 h later. If Cry1levels are increased during the daytime, when its levels
are rising, the clock will speed up as shown by the positive
daytime amplitude in Fig. 5. Conversely, an increase in Cry1
during the night will prolong its decline and slow down theclock, consistent with the negative nighttime amplitude. In
general, PRCs cross the x axis in the negative direction
near the abundance peak, and they cross in the positive direction near the abundance minimum.




'ES D









 









no P

35& QR 3

35& 3





NMDA PRCNMDA data
NPY PRC






P





NPY data

5RUα E

&U\ F

















































 






















5HY(UEα G










(ES H











FIGURE 5 Phase response curves for the fivemRNA components in the model. PRCs were
calculated both with (‘‘P’’) and without (‘‘no P’’)
probability constraints, and shaded regions contain
80% of the curves from a Monte Carlo ensemble.
The left y axis shows the scale for the ‘‘P’’
ensemble; the right axis shows the scale for the‘‘no P’’ ensemble. The calculated PRC for an increase in Cry1 expression is qualitatively similar
to a measured PRC for neuropeptide Y (NPY),
whereas the calculated PRCs for Dbp and RevErba
share some features in common with the PRC for
NMDA. Similar plots for the protein components
can be found in Fig. S6. To see this figure in color,
go online.

 













 


















Biophysical Journal 107(6) 1462–1473

1470

Jolley et al.

often be less able to reproduce systems-level data than asloppy model parameterized at the systems level.
This is not to say that individual parameter values are irrelevant to systems-level predictions. For example, because concentration scales are arbitrary in our dimensionless model,mRNA or protein degradation rates can increase significantly,
as long as they are offset by increases in the c parameters that
control the protein’s ability to regulate other genes. Phase
response calculations, however, require taking a derivativeof the phase with respect to a perturbation to a concentration
variable—the absolute value of these dimensionless variables
determines the scale of the PRC. The extreme variability of
the c parameters related to E4bp4 (c3 and c7 as inputs, c9
and c10 as inputs; see Fig. 4) is therefore reflected in the
wildly differing scales of E4bp4 PRCs in Fig. 5.
Our solution to this problem is facilitated by recent
genome-scale surveys of the synthesis and degradation rates
of proteins and mRNAs. Even if the components of a modelcannot be identified with any of the measured species, large
data sets can be used to construct previous probability distributions for model parameters. Fig. 2 indicates that the presence of probability constraints results in only minor changes
to the model outputs, whereas Figs. 3 and 4 show dramatic
decreases in parameter variability.

&5<GHOD\ &7



2SWLPL]HG
















FIGURE 6 Effects of weakening the D-box regulation of Cry1 The x axis
indicates the increase in period length (beyond the original ~24 h), and the yaxis shows the increase in the CRY1 phase (from the original CT 14),
normalized to circadian time. For the parameter set optimized with probability constraints, weaker D-box regulation results in a longer oscillation
period and a later CRY1 peak (relative to the period length). When eight
random samples were drawn from a Monte Carlo ensemble, the same positive trend was observed in seven of them. Shaded regions show convex
hulls containing the period/delay curves for parameters drawn from the
MC ensemble; the outermost hull contains 90% of points. To see this figure
in color, go online.

function largely measured the shape of the limit cycle. The
abundance oscillations shown in Fig. 2 show behavior on
the limit cycle and are quite consistent. PRCs (Fig. 5)
show the effects of small perturbations from the limit cycle
and are still fairly consistent. The Cry1/2 rescue experiments perturb the system significantly from the wild-type
limit cycle, and only qualitative predictions can be madewith any confidence.
Cleaning up sloppy models
Model sloppiness poses a challenge for systems biology—if
a model is parameterized to systems-level data, estimates
for the values of individual parameters are often unreliable.
Because the physical environment of the cell differs fromthe dilute solutions used in most in vitro assays, experimental estimates of parameter values will often be of limited
utility. Furthermore, whereas the ‘‘sloppy’’ eigendirections
are severely underconstrained by systems-level data, the‘‘stiff’’ eigendirections are often more tightly constrained
than they would be by direct experimental measurement(3). A model built with directly measured parameters will
TABLE 2
MutantDbp/Tef/Hlf
Rora/b
Cry1/2
RevErba/b
E4bp4

Knockout predictions
Our simplified model eliminates much of the redundancy of
the real clock system and is less resilient to gene knockouts
than the actual system. Knockouts of Cry1/2 are arrhythmic
because of the network topology (29); if C is removed fromthe diagram in Fig. 1 C, no closed negative feedback loops
persist. The Rev-Erba knockout is required to be arrhythmic
by our scoring scheme (15), although many parameter sets
showed damped oscillations to a stable fixed point. Other
possible knockouts have not been as well-characterized
experimentally and are discussed below and in Table 2.
Dbp/Tef/Hlf

If these D-box activators are knocked out, Rora also will not
be synthesized; because A and B are activators required forthe transcription of Cry1, this precludes oscillations. The
concentrations of Rora/b, Cry1, and E4bp4 should decrease
dramatically. If this mutation does not affect these levels atall (particularly Rora/b, which contain D-boxes), then the

Effects of modeled knockout mutations
Predicted rhythmicity

Predicted expression changes

Experimental
observation

Reference(s)Arrhythmic due to model topology
Damped oscillations
Arrhythmic due to model topology
Arrhythmic due to constraint; sometimes shows
damped oscillations
RhythmicRor, Cry1, E4bp4 all decrease dramaticallyE4bp4 decreases
Dbp/Tef/Hlf, Ror, RevErb, E4bp4 all increase
Dbp/Tef/Hlf decreases; others do not change much

Rhythmic
Rhythmic?
Arrhythmic
Arrhythmic(27)(54,55)(29)(15)

Biophysical Journal 107(6) 1462–1473

No effect

Clock Model with Daytime Elements

system may contain an additional D-box activator (either
oscillatory or constitutive) that is compensating for theloss of Dbp/Tef/Hlf. If instead, Rora/b levels do decrease
and the phenotype is similar to the removal of D-box regulation from Cry1 (Fig. 6) with a later Cry1 phase and alonger period, this suggests that D-box mediated expression
is affected but RRE-mediated expression is not. In such acase, it is possible that the RRE is being activated by apathway other than the D-box-mediated ROR proteins.
Rora/bKnocking out the Ror genes does not abolish Cry1 (which isalso activated by DBP/TEF/HLF), but it does abolishE4bp4. The remaining three genes do not form any closed
feedback loops, so sustained oscillations are impossible.
Members of the MMC ensemble with probability andknockout constraints generally show damped oscillations;in ~ 30% of cases this damping is fairly slow (i.e., morethan a few oscillation periods are observed) and population-level oscillations may still be possible, similar to
Per1 and Cry1 single-knockout mutants (53). Levels of
E4bp4 decrease dramatically, but changes to Cry1 aresmaller. This is probably because the Ror genes have incoherent effects on Cry1, both activating it directly and repressing it indirectly via E4bp4. If arrhythmicity and a decreasein E4bp4 are not observed, then RRE expression may be
activated by a different pathway. Experimentally, the Roraknockout shows a shortened circadian period (54), whereas
the Rorb knockout shows a lengthened one (55); the doubleknockout has (to our knowledge) not yet been characterized.
E4bp4

When parameter sets are drawn from the MMC ensemble
with probability and knockout constraints, simulated knockouts of E4bp4 have no effect on the period, the fit quality, or
the correctness of other knockout phenotypes (RevErba,
Cry1DD, Cry1DRRE). This is consistent with our interpretation of Figs. 4 and 5—E4bp4 may function primarily as an
output gene, rather than a core component of the clock.

DISCUSSION
A new view of the mammalian clock
Models that focus on the central Per-Cry/Bmal1-Clockfeedback loop have been successful at reproducing several
aspects of the circadian clock, such as gene knockout phenotypes (10,12–14) and entrainment to light stimuli (10,11).
The recent finding that Cry1 is best regarded as being coregulated by the RRE and the D-box (rather than the E/E’box) will require a different topological organization for
future models. Although our model is less detailed thanmany in the literature, it more correctly describes the transcriptional regulation of Cry1 (18). Our decision not toinclude Per1 and Per2 in the model was motivated primarily

1471

by the current incomplete state of knowledge about thedetailed biochemical mechanisms of these two importantcomponents, in particular the unusually long delay between
their mRNA and protein expression peaks (22,35,36). As
more reliable experimental characterizations of these genes
and their products become available, it will be possible toincorporate them into our modeling framework.
Our results suggest that the current understanding of the
D-box and its role is far from complete. When the modelis optimized with parameter constraints, the parameters governing the effects of E4BP4 (a D-box repressor) on Cry1 andRora take on extremely low values (Fig. 4), suggesting thatthe observed expression patterns for these genes can be fit
without any significant help from E4bp4. The inconsistent,
extremely low-amplitude PRC for E4bp4 (Fig. 5) likewiseindicates that E4bp4 may function primarily as a circadian
output, rather than as a central part of the clockwork.
At the same time, experiments show that (at least in
cultured cells) D-box regulation is essential for the correct
timing of Cry1 expression (18). It is possible that the effectsof the D-box on the core clockwork are mediated solely by
the D-box activators Dbp/Tef/Hlf or that the essentiality of
E4bp4 depends on interactions that were not included in
our model. In any case, careful experimental studies of Dbox-binding transcription factors and their targets are likely
to yield valuable insights into clock function.
Experimental perturbations
PRCs are determined experimentally by perturbing theclock system using a stimulus—light, neurotransmitters,
etc.—that has an indirect effect on the core clock system,
mediated by signaling pathways. A more direct (andpossibly feasible) approach would be transient overexpression of a clock gene in cultured cells. Simple genetic circuits
such as an incoherent feed-forward loop can generate transient pulses of gene expression in response to an externalsignal (56,57). If an exogenous clock gene driven by thistype of circuit were transfected into a cell line containing
a bioluminescent circadian reporter, one could observe theeffects of an overexpression pulse on the clock’s phase.
This type of protocol avoids the problems of redundancy.
For example, our model gene b (a D-box driven activatorof the RRE) could correspond to Rora or Rorb, and otherredundant genes may yet be discovered. Overexpressionof a single redundant gene, however, may be sufficient toproduce the expected phase shift. This contrasts withknockout experiments in which all redundant versions of agene must be removed before an effect is observed.
Mechanism of chemically induced
phase resetting
Previous models of the mammalian clock have focused on
the light-stimulated induction of Per1/2 transcription as
Biophysical Journal 107(6) 1462–1473

1472

Jolley et al.
11. Leloup, J. C., and A. Goldbeter. 2003. Toward a detailed computationalmodel for the mammalian circadian clock. Proc. Natl. Acad. Sci. USA.
100:7051–7056.

the sole external input to the clock (10,11,58,59). Although
this is undoubtedly important, it is likely that other pathways
exist, particularly for nonphotic signals such as NPY.
Because it focuses on a different feedback loop structure
than what has traditionally been employed in clock models,
our model provides an ideal tool with which to search forother possible modes of perturbation.
Comparison of our calculated PRCs with experimentallyderived PRCs (Fig. 5), suggests that photic signals such as
NMDA may act by stimulating the synthesis of E/E’-box
controlled genes such as Dbp or RevErba, consistent withthe light stimulation of E/E’-box-regulated Per. In addition,
nonphotic signals such as NPY may act on the core clockwork by increasing the expression of Cry1.

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

SUPPORTING MATERIAL

17. Ukai-Tadenuma, M., T. Kasukawa, and H. R. Ueda. 2008. Proof-bysynthesis of the transcriptional logic of mammalian circadian clocks.
Nat. Cell Biol. 10:1154–1163.

Six figures, two tables, and supporting text are available at http://www.
biophysj.org/biophysj/supplemental/S0006-3495(14)00740-1

18. Ukai-Tadenuma, M., R. G. Yamada, ., H. R. Ueda. 2011. Delay infeedback repression by cryptochrome 1 is required for circadian clockfunction. Cell. 144:268–281.

C. C. J. and D. P. are supported by the RIKEN Foreign Postdoctoral
Researcher Program. Calculations used the RIKEN Integrated Cluster ofClusters. This research was supported by an intramural Grant-in-Aid
from the 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 Grant-in-Aid for Scientific Research on Innovative Areas(No. 3306) from the Ministry of Education, Culture, Sports, Science, andTechnology, Japan (H. R. U.).
The authors would also like to thank Rikuhiro Yamada for helpful
discussions.

REFERENCES
1. Dyson, F. 2004. A meeting with Enrico Fermi. Nature. 427:297–298.
2. Brown, K. S., and J. P. Sethna. 2003. Statistical mechanical approaches
to models with many poorly known parameters. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 68:021904.
3. Gutenkunst, R. N., J. J. Waterfall, ., J. P. Sethna. 2007. Universallysloppy parameter sensitivities in systems biology models. PLOS Comput. Biol. 3:1871–1878.
4. Daniels, B. C., Y. J. Chen, ., C. R. Myers. 2008. Sloppiness, robustness, and evolvability in systems biology. Curr. Opin. Biotechnol.
19:389–395.
5. Raue, A., C. Kreutz, ., J. Timmer. 2009. Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics. 25:1923–1929.
6. Press, W. H. 2007. Numerical Recipes: The Art of Scientific
Computing. Cambridge University Press, New York.
7. Joshi, M., A. Seidel-Morgenstern, and A. Kremling. 2006. Exploiting
the bootstrap method for quantifying parameter confidence intervals
in dynamical systems. Metab. Eng. 8:447–455.

12. Becker-Weimann, S., J. Wolf, ., A. Kramer. 2004. Modeling feedback
loops of the mammalian circadian oscillator. Biophys. J. 87:3023–
3034.
13. Mirsky, H. P., A. C. Liu, ., F. J. Doyle, 3rd. 2009. A model of the cellautonomous mammalian circadian clock. Proc. Natl. Acad. Sci. USA.
106:11107–11112.
14. Kim, J. K., and D. B. Forger. 2012. A mechanism for robust circadiantimekeeping via stoichiometric balance. Mol. Syst. Biol. 8:630.
15. Cho, H., X. Zhao, ., R. M. Evans. 2012. Regulation of circadianbehaviour and metabolism by REV-ERB-a and REV-ERB-b. Nature.
485:123–127.

19. Hogenesch, J. B., W. K. Chan, ., C. A. Bradfield. 1997. Characterization of a subset of the basic-helix-loop-helix-PAS superfamily that interacts with components of the dioxin signaling pathway. J. Biol. Chem.
272:8581–8593.
20. Gekakis, N., D. Staknis, ., C. J. Weitz. 1998. Role of the CLOCK protein in the mammalian circadian mechanism. Science. 280:1564–1569.
21. Ueda, H. R., S. Hayashi, ., S. Hashimoto. 2005. System-level identification of transcriptional circuits underlying mammalian circadianclocks. Nat. Genet. 37:187–192.
22. Yoo, S. H., C. H. Ko, ., J. S. Takahashi. 2005. A noncanonical E-box
enhancer drives mouse Period2 circadian oscillations in vivo. Proc.
Natl. Acad. Sci. USA. 102:2608–2613.
23. Falvey, E., L. Marcacci, and U. Schibler. 1996. DNA-binding specificity of PAR and C/EBP leucine zipper proteins: a single amino acidsubstitution in the C/EBP DNA-binding domain confers PAR-likespecificity to C/EBP. Biol. Chem. 377:797–809.
24. Harding, H. P., and M. A. Lazar. 1993. The orphan receptor Rev-ErbA
alpha activates transcription via a novel response element. Mol. Cell.
Biol. 13:3113–3121.
25. Preitner, N., F. Damiola, ., U. Schibler. 2002. The orphan nuclear receptor REV-ERBalpha controls circadian transcription within the positive limb of the mammalian circadian oscillator. Cell. 110:251–260.
26. Ueda, H. R., W. Chen, ., S. Hashimoto. 2002. A transcription factor
response element for gene expression during circadian night. Nature.
418:534–539.
27. Gachon, F., P. Fonjallaz, ., U. Schibler. 2004. The loss of circadian
PAR bZip transcription factors results in epilepsy. Genes Dev.
18:1397–1412.
28. Storn, R., and K. Price. 1997. Differential evolution—a simple and efficient heuristic for global optimization over continuous spaces. J. Glob.
Optim. 11:341–359.
29. van der Horst, G. T., M. Muijtjens, ., A. Yasui. 1999. Mammalian
Cry1 and Cry2 are essential for maintenance of circadian rhythms. Nature. 398:627–630.

8. St John, P. C., and F. J. Doyle, 3rd. 2013. Estimating confidence intervals in predicted responses for oscillatory biological models. BMC
Syst. Biol. 7:71.

30. Sharova, L. V., A. A. Sharov, ., M. S. Ko. 2009. Database for mRNAhalf-life of 19 977 genes obtained by DNA microarray analysis ofpluripotent and differentiating mouse embryonic stem cells. DNA
Res. 16:45–58.

9. Metropolis, N., and S. Ulam. 1949. The Monte Carlo method. J. Am.
Stat. Assoc. 44:335–341.

31. Schwanhäusser, B., D. Busse, ., M. Selbach. 2011. Global quantification of mammalian gene expression control. Nature. 473:337–342.

10. Forger, D. B., and C. S. Peskin. 2003. A detailed predictive model ofthe mammalian circadian clock. Proc. Natl. Acad. Sci. USA.
100:14806–14811.

32. Nalefski, E. A., E. Nebelitsky, ., S. R. Gullans. 2006. Single-molecule
detection of transcription factor binding to DNA in real time: specificity,
equilibrium, and kinetic parameters. Biochemistry. 45:13794–13806.

Biophysical Journal 107(6) 1462–1473

Clock Model with Daytime Elements

1473

33. Kolesov, G., Z. Wunderlich, ., L. A. Mirny. 2007. How gene order isinfluenced by the biophysics of transcription regulation. Proc. Natl.
Acad. Sci. USA. 104:13948–13953.

48. Huhman, K. L., and H. E. Albers. 1994. Neuropeptide Y microinjected
into the suprachiasmatic region phase shifts circadian rhythms in constant darkness. Peptides. 15:1475–1478.

34. Strogatz, S. H. 1994. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. Addison-Wesley, Reading, MA.
35. Debruyne, J. P., E. Noton, ., S. M. Reppert. 2006. A clock shock:mouse CLOCK is not required for circadian oscillator function.
Neuron. 50:465–477.

49. Albrecht, U., Z. S. Sun, ., C. C. Lee. 1997. A differential response of
two putative mammalian circadian regulators, mper1 and mper2, tolight. Cell. 91:1055–1064.

36. Lee, C., J. P. Etchegaray, ., S. M. Reppert. 2001. Posttranslational
mechanisms regulate the mammalian circadian clock. Cell. 107:
855–867.
37. Mitsui, S., S. Yamaguchi, ., H. Okamura. 2001. Antagonistic role of
E4BP4 and PAR proteins in the circadian oscillatory mechanism.
Genes Dev. 15:995–1006.
38. Stratmann, M., F. Stadler, ., J. A. Ripperger. 2010. Flexible phaseadjustment of circadian albumin D site-binding protein (DBP) gene
expression by CRYPTOCHROME1. Genes Dev. 24:1317–1328.
39. Schneider, C. A., W. S. Rasband, and K. W. Eliceiri. 2012. NIH Image
to ImageJ: 25 years of image analysis. Nat. Methods. 9:671–675.
40. Nelder, J. A., and R. Mead. 1965. A simplex-method for function minimization. Comput. J. 7:308–313.
41. Kramer, M. A., H. Rabitz, and J. M. Calo. 1984. Sensitivity analysis ofoscillatory systems. Appl. Math. Model. 8:328–340.
42. Dunlap, J. C., J. J. Loros, and P. J. DeCoursey. 2004. Chronobiology:
Biological Timekeeping. Sinauer Associates, Sunderland, MA.
43. Prosser, R. A., R. R. Dean, ., J. D. Miller. 1993. Serotonin and the
mammalian circadian system: I. In vitro phase shifts by serotonergic
agonists and antagonists. J. Biol. Rhythms. 8:1–16.
44. Hamada, T., S. Shibata, ., S. Watanabe. 1993. Effect of somatostatinon circadian rhythms of firing and 2-deoxyglucose uptake in rat suprachiasmatic slices. Am. J. Physiol. 265:R1199–R1204.
45. Trachsel, L., H. C. Heller, and J. D. Miller. 1995. Nicotine phase-advances the circadian neuronal activity rhythm in rat suprachiasmatic
nuclei explants. Neuroscience. 65:797–803.
46. Biello, S. M., D. A. Golombek, and M. E. Harrington. 1997. Neuropeptide Y and glutamate block each other’s phase shifts in the suprachiasmatic nucleus in vitro. Neuroscience. 77:1049–1057.
47. Asai, M., S. Yamaguchi, ., H. Okamura. 2001. Visualization of mPer1
transcription in vitro: NMDA induces a rapid phase shift of mPer1 genein cultured SCN. Curr. Biol. 11:1524–1527.

50. Shearman, L. P., M. J. Zylka, ., S. M. Reppert. 1997. Two period homologs: circadian expression and photic regulation in the suprachiasmatic nuclei. Neuron. 19:1261–1269.
51. Shigeyoshi, Y., K. Taguchi, ., H. Okamura. 1997. Light-induced
resetting of a mammalian circadian clock is associated with rapid induction of the mPer1 transcript. Cell. 91:1043–1053.
52. Inagaki, N., S. Honma, ., K. Honma. 2007. Separate oscillating cellgroups in mouse suprachiasmatic nucleus couple photoperiodically tothe onset and end of daily activity. Proc. Natl. Acad. Sci. USA.
104:7664–7669.
53. Liu, A. C., D. K. Welsh, ., S. A. Kay. 2007. Intercellular coupling
confers robustness against mutations in the SCN circadian clock
network. Cell. 129:605–616.
54. Sato, T. K., S. Panda, ., J. B. Hogenesch. 2004. A functional genomics
strategy reveals Rora as a component of the mammalian circadianclock. Neuron. 43:527–537.
55. André, E., F. Conquet, ., M. Becker-André. 1998. Disruption of retinoid-related orphan receptor beta changes circadian behavior, causes
retinal degeneration and leads to vacillans phenotype in mice.
EMBO J. 17:3867–3877.
56. Mangan, S., and U. Alon. 2003. Structure and function of the feed-forward loop network motif. Proc. Natl. Acad. Sci. USA. 100:11980–
11985.
57. Basu, S., R. Mehreja, ., R. Weiss. 2004. Spatiotemporal control ofgene expression with pulse-generating networks. Proc. Natl. Acad.
Sci. USA. 101:6355–6360.
58. Okamura, H., S. Miyake, ., G. T. J. van der Horst. 1999. Photic induction of mPer1 and mPer2 in cry-deficient mice lacking a biologicalclock. Science. 286:2531–2534.
59. Zylka, M. J., L. P. Shearman, ., S. M. Reppert. 1998. Three period homologs in mammals: differential light responses in the suprachiasmatic
circadian clock and oscillating transcripts outside of brain. Neuron.
20:1103–1110.

Biophysical Journal 107(6) 1462–1473