Emergence of Noise-Induced Oscillations in the Central Circadian Pacemaker

Bmal1 is the only known clock gene whose loss-of-function leads to complete loss of circadian rhythmicity in wheel-running behavior [7]. In addition, circadian rhythms of Per1 and Per2 mRNA expression are disrupted in the SCN of Bmal1^−/− mice [7]. To analyze the role of Bmal1 in cell- and tissue-autonomous circadian properties in central and peripheral oscillators, we measured PER2::LUC bioluminescence rhythms in various tissue explants from wild type (WT; Bmal1^+/+) and homozygous Bmal1 mutant (Bmal1^−/−) mice. Surprisingly, a persistent, yet highly variable, rhythmicity emerged in Bmal1^−/− SCN explants (Figure 1 and Video S1), whereas all peripheral tissues failed to generate PER2::LUC circadian rhythms (Figure 1A).

WT and Bmal1^−/− SCN explants show persistent PER2::LUC bioluminescence rhythms for more than 35 d in culture (Figure 1B and 1C; also see Datasets S1 and S2). Fast Fourier transform (FFT) spectrograms of the records show a consistent frequency (about one cycle per day) for the WT SCN; however, more variable frequency components are observed in the Bmal1^−/− SCN (Figure 1B and 1C, middle). Double-plotted raster plots illustrate the stable PER2::LUC rhythmicity in the WT SCN and the unstable rhythms in Bmal1^−/− SCN (Figure 1B and 1C, right). FFT analysis shows that the dominant periodicity in WT SCN is tightly clustered around 24 h, while periodicity is much more variable, ranging from about 14 to 30 h, in Bmal1^−/− SCN (Figure 1D). Because of the highly variable periods and instability of the rhythmicity from Bmal1^−/− SCN, we refer to these fluctuations as stochastic. All tissues from WT littermates, including the SCN, pituitary, liver, lung, and cornea, displayed normal circadian rhythms similar to those observed in previous studies (Figure 1A) [4],[28].

Due to the stochastic nature of the rhythms in the Bmal1^−/− SCN, period estimates based on averages of long time series do not adequately describe the cycle-to-cycle variability. Hence, we measured individual peak-to-peak intervals of the PER2::LUC expression patterns from each SCN explant (Figure 1E, left). The mean inter-peak intervals observed in the WT SCN explants were circadian at about 24 h (24.3±1.03 SD h). These values were similar to the periods estimated by Levenberg-Marquardt (LM) curve fitting of the entire time series using the LumiCycle Analysis Program (WT = 24.04±0.09 SD h, n = 23). By contrast, the Bmal1^−/− SCN showed a significantly shorter average inter-peak interval length of 17.8 h with a much greater variance (SD = 5.54 h, n = 19) than WT SCN explants.

Using the inter-peak interval data, we determined the Serial Correlation Coefficient of successive intervals, r_s (Figure 1E, right). A negative serial correlation reflects the likelihood that a long cycle will be followed by a short cycle, or vice versa, which is a characteristic feature of a functional pacemaker-driven system [39]. The average serial correlation coefficients for WT SCN explants were negative (mean r_s = −0.17±0.0606 SEM, t-test, p<0.01) as would be expected from a circadian pacemaker-driven process (Figure 1E, right). However, the coefficient for the Bmal1^−/− SCN was slightly positive and marginally significant (mean r_s = 0.07±0.0336 SEM, t-test, p<0.05). This weak serial correlation coefficient differs from the pacemaker-driven processes seen in WT SCN explants and is consistent with either an oscillator with a highly labile period or a “random walk” process [40].

To determine whether the stochastic PER2::LUC rhythms from Bmal1^−/− SCN explants are cell autonomous, we studied PER2::LUC bioluminescence at the single-cell level (Figure 2). We first imaged the overt bioluminescence expression patterns from Bmal1^−/− SCN explants (Figure 2A, Video S1, Datasets S3 and S4) and analyzed bioluminescence from individual neurons (Figure 2B and 2C). The SCN cells in an intact organotypic slice were tightly synchronized and exhibited stochastic rhythms comparable to those seen in the Bmal1^−/− SCN explants using luminometry. Separately, we cultured dissociated SCN neurons and imaged bioluminescence from individual cells. In contrast to the cells in an intact organotypic slice, dissociated Bmal1^−/− neurons did not express detectable circadian rhythms (Figure 2D, 2E, 2F, Video S2, Datasets S5 and S6). A total of 243 out of 243 Bmal1^−/− SCN cells were equal to or below a threshold FFT amplitude for rhythmicity independently derived from WT SCN cells. Although Bmal1^−/− SCN cells did not generate obvious circadian oscillations, they did express fluctuating levels of PER2::LUC—likely a reflection of Per2 transcriptional and post-transcriptional noise (Figure 2F). The lack of rhythmicity in dissociated SCN cells indicates that the stochastic rhythmicity observed in Bmal1^−/− SCN explants is not a cell-autonomous property and likely arises from intercellular network interactions (to be addressed later). Thus, the integrity of the cellular anatomy of the SCN is an important factor for stochastic rhythm generation in Bmal1^−/− SCN explants.

To explore possible origins of the stochastic rhythmicity in the Bmal1^−/− SCN, we used mathematical modeling of circadian oscillators to examine the role of intercellular coupling and molecular noise. The key clock components in our mathematical model and their possible interactions within a single SCN cell are depicted in Figure 3A. This model is similar to the Forger-Peskin stochastic model [34] but with three refinements. In the current simulation, we (1) explicitly modeled the binding of CLOCK:BMAL to regulatory regions of the Per1, Per2, Cry1, and Cry2 genes; (2) explicitly modeled the interaction of BMAL with CRY1 or CRY2 and its subsequent effects on transcriptional regulation; and (3) updated the rates of degradation of mRNAs and proteins using measurements from recent experiments [41]. An explanation of these changes, as well as a list of rate constants, can be found in Figure S7 (also see Protocols S1 and S2).

We accounted for the heterogeneity of circadian period of SCN neurons by choosing all rate constants from Gaussian distributions. When no Gaussian variation is introduced to the rate constants (Figure 3B, “0%”), the standard deviation for period of the population of oscillators is 1.03 h, which reflects variability in period due to the stochastic model (i.e., no rate constant variability). As the standard deviation of the rate constants is increased, the standard deviation of the circadian period from a population of oscillators increases further (Figure 3B). When we chose the standard deviation of these Gaussians to be within 5% of their mean values, the resulting period of the model simulations had a standard deviation (SD = 1.54 h) similar to the experimental data obtained from dispersed WT SCN neurons, which have a standard deviation of 1.32 h for circadian period (Figure 3C). When these individual neurons were coupled in a network of cells (see below), the total mRNA and protein concentrations followed time courses similar to those found experimentally (Figure 3D).

This stochastic model was then used to simulate Bmal1^−/− SCN isolated cells and the Bmal1^−/− coupled SCN network, respectively. As a first step, we examined the effects of reducing the total BMAL activator concentration on the generation of rhythms in simulations of isolated cells. As the BMAL activator abundance was gradually decreased, the distribution of periods widened and the mean period shortened until a complete loss of rhythms was observed below 20% BMAL activator abundance in a population of uncoupled cells (Figure 4A and 4B). This behavior is consistent with bifurcation diagrams of single cell simulations (Figure 4C) that show the existence of a Hopf bifurcation at approximately 22% of total BMAL. Above this percentage, oscillations can be shown to emerge in our simulations. At 0% total BMAL activator abundance no element in the model could initiate transcription; thus, some level of residual BMAL activator must be present in the model in order for clock gene expression to occur.

In Bmal1 null mutants, other functional E-box activator proteins may be present. The most likely candidate for “residual BMAL” activity in Bmal1 mutants is BMAL2 (also known as MOP9; RefSeq: Arntl2NM_172309↗), a paralog of BMAL1, which is regionally co-expressed in the SCN and forms transcriptionally active complexes with CLOCK [7]. BMAL1 and BMAL2 show similar sensitivity to CRY1-mediated transcriptional repression (Figure S1; for methods, see Text S1), which suggests that the core negative feedback loop of the circadian oscillator could remain functional in the presence of BMAL2. Using quantitative PCR methods that normalize the amplification of Bmal1 and Bmal2 RNA, we estimate that Bmal2 mRNA levels are approximately 10% of the level of Bmal1 and that Bmal2 expression levels are unaffected by Bmal1 loss-of-function (i.e., no compensation of Bmal2 in Bmal1^−/− mutants) (Figure S2; for methods, see Text S1). Indeed recent work has shown that Bmal2 can substitute for Bmal1 if ectopically driven at high levels in Bmal1^−/− mice [42]. Since BMAL1 and BMAL2 appear to be functionally equivalent, we use the term “BMAL” in our simulations to represent the sum total of BMAL1 and BMAL2.

On the basis of the expression analysis, we chose the 10% BMAL activator concentration to simulate the Bmal1^−/− SCN neuron because Bmal2 represents approximately 10% of the wild-type level of Bmal1; and in Bmal1^−/− cells, the remaining Bmal2 would contribute about 10% of the total BMAL1 and BMAL2 activity. Our simulated and experimental isolated cell rhythms show remarkable agreement at the cell autonomous level. In Figure 4D, we show representative experimental records from dissociated WT SCN cells (first row) along with their associated spectrograms (second row). Similarly, the third and fourth rows in parallel display the computational simulations for the isolated WT SCN neurons with their associated spectrograms. Comparable traces are shown for isolated Bmal1^−/− cells in Figure 4E. Thus, the 10% BMAL activator simulations faithfully capture the stochastic behavior of PER2::LUC expression seen in dissociated Bmal1^−/− SCN neurons.

Although there are many mathematical models of the SCN that incorporate coupling of a population of circadian oscillators [16],[19]–[21], very few models have studied coupling in a population of oscillators in which the individual oscillators are stochastic in nature [26],[27]. The results in this paper represent the most detailed stochastic simulation of coupled SCN oscillators to date. To model a population of stochastic oscillators in the SCN, we utilized a group of 100 “cell autonomous” stochastic oscillators as described in the previous section and coupled the population of oscillators (see Figures 5, S3, and S4). For the coupling pathway, we used a model of the VIP signaling pathway in the SCN, whereby a rhythmic coupling agent (CA) is released by one SCN cell and affects neighboring SCN cells by activating an adenosine 3′,5′-monophosphate (cAMP) signaling pathway, which ultimately activates cAMP response element binding (CREB) protein and cAMP response elements (CRE) on the Per1 and Per2 promoters to induce PER1 and PER2 proteins. In our proposed model, the intercellular coupling is all-to-all. We assumed that CRE activation of Per1 and Per2 could occur only when repressors were not bound as suggested in previous experiments [43]. Thus, this is the equivalent of allowing CRY to repress CREB-mediated PER production.

When coupling is included in a population of WT single-cell models, this coupled cellular network was able to produce coherent ∼24-h rhythms (Figure 5B, top, and 5C, bottom) as seen in WT SCN explants and in other coupled oscillator models of the SCN [16],[19]–[21]. Remarkably, as seen in experiments, rhythmicity also emerged when a population of nonrhythmic Bmal1^−/− SCN simulated neurons was coupled together (Figure 5C, top and middle). To compare these rhythms to those found experimentally, we calculated the peak-to-peak intervals of the rhythms from the SCN network with 100%, 20%, and 10% BMAL activator concentrations (Figure 5B). The distribution of these inter-peak intervals is similar to experimental measurements from WT and Bmal1^−/− explant rhythms (Figure 5B; compare with Figure 1E). The observed decrease in period in the simulated SCN explants with lower activator concentrations is similar to that observed in simulated single cells, implying that the short period of the Bmal1^−/− explant may be attributable to a shortened intracellular, single-cell periodicity. However, given the general absence of rhythms in individual uncoupled Bmal1^−/− cells, it is interesting that coupling of these intrinsically nonrhythmic cells (10% BMAL) can generate stochastic rhythms.

While only three sets of simulations are shown in Figure 5B and 5C, we systematically varied the coupling strength and BMAL concentration to determine the amplitude of rhythms and synchrony of the network in Figure S3. This shows that the stochastic rhythmicity can be seen over a wide range of values for the coupling strength. In addition, we note that the relative amplitudes of the Bmal1^−/− simulations are reduced when compared with those found experimentally. Figure S3 shows that the amplitude of these rhythms depends greatly on the coupling strength and mechanism. In particular as the coupling becomes nonlinear, the amplitude increases. Including other nonlinear aspects of the coupling (e.g., electrical activity) could also increase the amplitude. In addition, there are likely to be other transcriptional inputs to the Per and Cry genes that could affect amplitude and that were not included in the model.

Since the emergent rhythms in the Bmal1^−/− SCN explants were stochastic, we hypothesized that the molecular noise, inherent in all biochemical reactions and included in our current model through our use of the Gillespie stochastic simulation algorithm, could play an important role. To test this hypothesis, we developed a deterministic version of our proposed model to understand the implications of the noise. The deterministic model represents a version of our proposed model without any noise. Each cell in the deterministic model had its parameters chosen from a distribution as was done in the stochastic model. Once the parameter was chosen it remained fixed for the rest of the simulation. For example, the bifurcation analysis in Figure 4C used the deterministic model with its mean parameters. The stochastic and deterministic models produced equivalent results when the number of each molecule species (Figure 5D) in the stochastic model was increased. This equivalence serves as a check on the accuracy of both the stochastic and deterministic simulation approaches. The output of 100 simulated cells is shown using both the stochastic and deterministic models in Figure 5C and 5E. The total BMAL concentration was varied from 10% to 20% and finally 100% of the activator level in WT SCN, and the ensemble averages from each simulation are plotted in yellow. Strikingly, when a population of cells using the deterministic model was coupled, rhythmicity could not be recovered at 10% BMAL activator levels and only damped rhythms could be seen at 20% BMAL activator levels (Figure 5E).

The discrepancy between the stochastic and deterministic models clearly indicates that noise is a necessary but not sufficient condition for the rhythmicity observed in the Bmal1^−/− SCN explants. We reach this conclusion because the deterministic model fails to show any sustained rhythmicity at the 10% or 20% level of BMAL activator concentration. Our analysis indicates that noise alone could not restore the rhythmicity in individual cells, nor is the coupling mechanism alone sufficient to induce oscillations in a population of cells without noise. Thus, what is necessary to induce rhythms in a population is a joint requirement for both molecular noise and intercellular coupling.

The coupled stochastic model of the Bmal1^−/− SCN predicts that both noise and coupling are required to generate stochastic rhythms. To test the role of coupling in SCN explants, we used a variety of agents previously reported to uncouple SCN neurons. We first used tetrodotoxin (TTX) treatment on Bmal1^−/− SCN rhythms (Figure 6). TTX prevents action potentials by selectively and reversibly blocking voltage-dependent Na⁺ channels. TTX application has been shown to desynchronize neurons within an intact SCN [44], suggesting that action potentials and/or consequent neuronal transmission are required for maintaining SCN synchrony. Under continuous TTX administration, WT SCN tissue showed persistent PER2::LUC rhythms. However, the amplitude of the rhythm gradually diminished cycle by cycle (damping), which is attributable to intercellular desynchrony and a reduction in amplitude at the individual oscillator level [44]. In contrast, the Bmal1^−/− SCN exhibited an immediate loss of stochastic rhythmicity of PER2::LUC output when treated with TTX. The PER2::LUC rhythms returned immediately upon removal of TTX in both WT and Bmal1^−/− SCN explants (Figures 6B and S5). To determine whether the loss of stochastic rhythmicity was due to desynchrony of individual cellular rhythms, we used bioluminescence imaging to measure single-cell behavior during TTX treatment (Figures 6C and S5). Uncoupling cells by TTX treatment within an organotypic Bmal1^−/− SCN slice caused a loss of stochastic rhythmicity at the single-cell level with PER2::LUC patterns similar to those of SCN neurons in dispersed culture. Thus, TTX treatment abolished stochastic oscillations by preventing the sustainment of rhythmicity at the single-cell level.

Using the coupled stochastic model of the SCN, we mimicked the experiments that blocked coupling via TTX and then restored the coupling by removing TTX (Figure 6D). Our simulation results reproduce the experimental data for both the WT and Bmal1^−/− SCNs; these results demonstrate that abolition of coupling alone in the model can account for the observed patterns of PER2::LUC bioluminescence from Bmal1^−/− SCN when blocked with TTX.

Similar findings were observed when other pharmacological agents were used to uncouple the SCN neurons (Figure S6). For example, pertussis toxin treatment, which blocks G_i/o-mediated signal transduction, and bicuculline, which blocks GABA_A-receptors, both abolished PER2::LUC stochastic rhythmicity—these agents have previously been shown to uncouple SCN cells [45],[46].

A prominent pathway for coupling SCN neurons involves the VIP receptor-signaling pathway that activates cAMP and CREB-mediated induction of Per1 and Per2 transcription [37],[38],[43],[47],[48]. To test the role of the cAMP signaling pathway in coupling in SCN explants, we used MDL-12,330A (MDL), a potent inhibitor of adenylyl cyclase, and H-89, an inhibitor of cAMP-activated protein kinase (protein kinase A; PKA) (Figure 7). Previous studies have shown that MDL reduces cAMP concentrations to basal levels in SCN cells [47], and H-89 prevents VIP- or calcium-induced circadian transcription in SCN cells [49]. Similar to the effect of TTX treatment, WT SCN explants showed damped rhythms with prolonged exposure to these inhibitors of cAMP-mediated signal transduction (Figure 7B and 7C, left), and the Bmal1^−/− SCN explants showed an immediate loss of stochastic rhythmicity (Figure 7B and 7C, right). Thus, the cAMP pathway appears to be necessary for the generation of stochastic oscillations.

In the circadian clock mechanism, the half-life of the PER proteins contributes to the determination of circadian period, which is under the control of CK1ε/δ-mediated PER phosphorylation [50]–[52]. Modeling studies have also suggested that the Per2 feedback loop, in particular, may have a dominant role in setting the period of circadian oscillation [53]. The VIP/cAMP signaling pathway ultimately converges on Per induction, and thus the periodicity of stochastic rhythms in Bmal1-null mutant SCN may also depend upon PER. To test a role for PER proteins in determining the period of stochastic rhythms, we examined the effects of agents that inhibit CK1ε/δ to lengthen circadian period. We treated SCN explants with a potent protein kinase inhibitor, SP600125, which has been shown to lengthen circadian period and to inhibit CK1ε kinase activity [51]. SP600125 treatment lengthened WT SCN circadian period to over 30 h (Figure 8A, top and 8B, left) as shown previously. Interestingly, SP600125 also significantly lengthened the period (average inter-peak intervals) of the stochastic rhythms from Bmal1^−/− SCN explants (Figure 8A, bottom and 8B, right). Modeling experiments in which the phosphorylation rate of PER by CK1 was reduced also lengthened the periodicity of WT and Bmal1-mutant SCN simulations in a similar fashion (Figure 8C). Thus, the kinetics of the quasi-circadian periodicity that emerges from the SCN network are regulated by a PER-dependent process in which the phosphorylation rate of PER by CK1 can determine the periodicity of the emergent network oscillation.

All animal care and experimental treatments were approved and performed in strict accordance with Northwestern University and University of Texas Southwestern Medical Center guidelines for animal care and use. Mice were bred from PER2::LUC (http://jaxmice.jax.org/strain/006852.html↗) homozygous parents [4] mated with Bmal1^+/− mice [4],[7]. From the first generation of pups, mice carrying a copy of the Per2::luc (luc/+) and heterozygous for the mutation in Bmal1 were then crossed to obtain Bmal1^+/+, Bmal1^+/−, and Bmal1^−/− mutants harboring the PER2::LUC reporter.

At approximately 8–10 wk of age, mice were placed in individual running wheel cages, and activity was recorded using the ClockLab data collection system (Actimetrics, Wilmette, IL) [75]. After 2 wk in LD 12∶12, the mice were released into constant darkness (DD) for an additional 4 wk. Animals were then returned to LD 12∶12 for at least 2 wk before their tissues were harvested for bioluminescence experiments.

Animals were euthanized by cervical dislocation between ZT 11–13. The tissues were removed immediately and put in Hank balanced salt solution (HBSS; with 10 mM HEPES, 25 units/ml penicillin, and 25 µg/ml streptomycin) on ice. Brain slices containing the SCN were sectioned at 300 µm using a vibratome followed by scalpel dissection of the SCN, resulting in a piece of tissue about 1 mm×1 mm in size. The peripheral tissues were dissected into pieces approximately 1 mm³ in size, with the exception of the pituitary, which was cultured as a whole. Each dissected tissue was cultured on a Millicell culture membrane (PICMORG50, Millipore) with 1.2 ml DMEM medium (Cellgro), supplemented with 10 mM HEPES (pH 7.2), 2% B27 (Invitrogen), 25 units/ml penicillin, 25 µg/ml streptomycin, and 0.1 mM luciferin (Promega) as described previously [76]. Medium changes were performed by lifting the Millicell culture membrane and placing it into a new culture dish prepared with fresh medium. For peripheral tissues, forskolin (10 µM) was administered for ∼30 min to synchronize the cells before placement into fresh medium. For uncoupling experiments, TTX (1 µM), BIC (200 µM), PTX, (5 ng/ml), MDL (1 µM), and H-89 (10 µM) were used. To test the role of PER2 in intercellular coupling and period determination, the protein kinase inhibitor, SP600125 (25 µM), was used. Drugs were added during medium change and were left undisturbed until replacement with fresh medium. All reagents were purchased from Sigma-Aldrich.

Explant cultures were maintained at 36 °C either in an incubator or in a temperature-controlled room. The bioluminescence was continuously monitored with LumiCycle photomultiplier tube (PMT) detector systems (Actimetrics, Wilmette, IL). Dark counts from the PMT were measured with luciferin-containing medium alone and subtracted from overall bioluminescence. Bioluminescence was measured immediately upon placement in culture, continuing without interruption for >7 d.

Bioluminescence analyses were performed using the LumiCycle Analysis Program (Actimetrics, Wilmette, IL). Raw data were baseline fitted. Then peak-to-peak durations (inter-peak intervals) were measured by manually identifying individual peaks. Serial correlation coefficient estimates were calculated for 7 to 10 cycle epochs using methods described previously [39].

Bmal1^−/− SCN slices used for imaging were dissected from 4–10-d-old pups, cut by tissue chopper (Stoelting) to a thickness of 400 µm, and cultured on Millicell-CM membrane inserts (PICMORG50). SCN slices were maintained in culture for 2–3 wk before imaging. For TTX experiments, adult SCN explants were dissected as described in the previous section.

For preparation of dispersed SCN neurons, we used 2–5-d-old Bmal1^−/− pups (or WT controls from the same heterozygous breeding line). Cylindrical punches of unilateral SCN were made from 400 µm coronal sections using a 20-gauge needle. Punches from 2–6 mice were pooled in each preparation, and the experiment was performed twice for each genotype. Cells were dissociated using papain and cultured as previously described [1], except that medium contained 5% FBS instead of rat serum. SCN neurons were maintained in culture for 2–7 wk before imaging.

Bioluminescence imaging was performed as previously reported [3],[28],[77]. Just before imaging, medium was changed to fresh explant medium containing 1 mM luciferin. Culture dishes were sealed and placed on the stage of an inverted microscope (Olympus IX71) in a dark room. A heated lucite chamber around the microscope stage (Solent Scientific, UK) kept the cells at a constant 36°C. For experiments presented in Figure 2 and Videos S1 and S2, images were collected using an Olympus 4× XLFLUOR (NA 0.28) or UPlanSApo (NA 0.16) objective and transmitted to a CCD camera (Spectral Instruments SI800, Tucson, AZ) cooled to −92°C. For dispersed neurons, signal-to-noise ratio was improved by 2×2 binning of pixels. Images of 29.8 min exposure duration were collected at 30 min intervals for 7 d. SCN neuron viability was assessed by cell morphology and stability of average daily bioluminescence, and no differences were observed between WT and Bmal1^−/− cells. For TTX experiments, adult SCN slice images were collected with a dual microchannel plate intensified gallium arsenide phosphide XR/MEGA-10Z CCD camera (Stanford Photonics, Inc., Palo Alto, CA) cooled to −20°C.

Bioluminescence was analyzed using MetaMorph (Molecular Devices) as previously described [3],[28],[77]. For TTX experiments, CellCycle single-cell analysis software (Actimetrics, Wilmette, IL) was used. Cells that were clearly discriminable from adjacent cells and that remained bioluminescent for the entire experiment were selected for analysis. Bioluminescence time series were first imported into LumiCycle Analysis v. 2.31 (Actimetrics). A linear baseline was subtracted from raw data (polynomial order = 1). Due to high initial transients of luminescence in some cases, the first 12 h of data were excluded from analysis in all cells. Maximum spectral power value was determined using FFT-Relative Power, which was defined as relative spectral power density at the peak within the range of 0–36 h, i.e., proportion of total variance within a 0.14 cycles/day window centered at highest point within the range of 0–36 h. A cell was considered to show significant circadian periodicity when the spectral analysis indicated a peak in the circadian range (20–36 h) large enough such that a 0.14 cycles/day window centered on the peak accounted for at least 10% of the total variance in the record (FFT power spectrum, Blackman-Harris windowing, peak amplitude ≥0.1) as described previously [78]. This is the same criterion used in our previous study of mutant SCN neurons [28], chosen so as to include all clearly rhythmic WT cells. Comparisons significant by ANOVA (p<0.05) were further explored by pairwise t-tests. Period was defined as the period of the best-fit sine wave as estimated by a Levenberg-Marquardt algorithm.

Simulations used a revised version of the Forger-Peskin stochastic model of the mammalian circadian clock [14],[34]. We explicitly modeled the binding and unbinding of CRY1 or CRY2 (and any other bound proteins) with CLOCK:BMAL (which has the same biochemical properties as the CLOCK:BMAL1 or CLOCK:BMAL2 complexes, which were treated equivalently). Based on experimental evidence [28],[79], we assumed that the concentration of these protein complexes was constant (1,600 molecules per cell for the WT case). The original Forger-Peskin model [14] contained a binding event of the CRY proteins to activator on the promoter. This was replaced by a binding event of CLOCK:BMAL to the promoter. Consistent with fitting of experimental data [41] we assumed that there was just one active E-box on the promoters of Per1, Per2, Cry1, and Cry2. Reaction rates can be found in Figure S7.

We also modeled a CRE element where a factor (e.g., CREB) can bind and activate transcription of Per1 and Per2 or a secreted CA. The “coupling” factor is equivalent to AVP or VIP and mediates signaling between SCN neurons (see Figure 5). Since medium changes can stimulate transcription, presumably through this CRE element, we assumed that there was a constant low level of this factor in the network (one molecule per cell). Removing this constitutive activation had little effect on our simulations (unpublished data).

Due to the revised transcription mechanism, new rates of transcription were chosen to give qualitative agreement of the uncoupled model with protein concentration levels determined experimentally by Lee et al. (2001) [11]. Other parameters were fit to published experimental data on the mRNA and protein time courses used in previous models [14],[34]. Degradation rates for Per1/2 and Cry1/2 mRNAs and CRY1 and PER2 proteins were obtained directly from published experimental data [41]. This yielded a model with a slightly short period (22.9 h) that was scaled to give a 24-h period.

Simulations were performed using the Gillespie Algorithm as in Forger and Peskin (2005) [34] or MATLAB's ode15s with a maximum time step of 0.05 h. In the limit of large numbers of molecules, these simulations matched, as did a version of the model coded in Mathematica. Simulations typically lasted 500 h; the first 50–100 h were discarded to remove the effects of initial transients. Using MATLAB, spectrograms of the experimental and simulated single cell and slice data were produced. Specifically, the “spectrogram” command was used, with window size 128, overlap set to 100, and sampling set to 48 cycles per day. Simulated data were collected and analyzed in the same manner as the experimental data (Figure 4D, 4E). Network simulations contained 100 cells with mean field coupling. We then calculated the normalized amplitude and the order parameter defined by Garcia-Ojalvo and colleagues [80]. Codes written in the C programming language are available as Protocol S1.