Two strains of S. elongatus PCC 7942 were used for this study. AMC1300 is a wild-type derivative carrying a bacterial luciferase bioluminescent reporter of kaiBC expression. AMC1300 carries PkaiBC::luxAB at NS1 and PpsbAI::luxCDE at NS2, which enables the cells to produce the luciferase enzyme and the long-chain aldehyde substrate for the luminescent reaction (Chen et al., 2009). AMC408 carries a purF reporter (PpurF::luxAB at NS2, PpsbAI::luxCDE at NS1) (Nair et al., 2002). Prior to experimental measurements, all cultures were grown in test tubes in BG11M liquid medium at 30°C with shaking under cool white fluorescent bulbs (≈60 μmol photons m⁻² s⁻¹, Philips AltoII, Amsterdam, Netherlands).

To simulate different light-dark cycles, we used custom-built red LED arrays in a 96-well format (LEDs from superbrightleds.com, St. Louis, MO, cat. no. RL5-R12008; 96-well plates from Corning, Corning, NY, cat. no. 3916). The LEDs were mounted into a hollowed-out 96-well plate and the tails of the LEDs were soldered into a circuit board, where they were wired in parallel in groups of four to the analog inputs of an Arduino Mega 2560 microcontroller. Another hollow 96-well plate was glued to the bottom of the LED-carrying plate, beneath the LEDs, in order to create a light baffle and prevent light leakage between the wells. These devices were placed ≈2 mm above a black 96-well plate containing cells growing on BG11-agar, such that every well of the growth plate received illumination from a single LED (≈1.8 cm between LED and agar surface). The growing plate was sealed with transparent film, with holes punctured above each well to provide aeration. Custom Arduino scripts were written to administer appropriate light-dark schedules to cells in the growth plates. Every well received the same light intensity (1.8 V across each LED, ≈200 μmol photons m⁻² s⁻¹) in the light portion of the day. The temperature of the agar was $31.6\pm 0.8$ °C underneath LEDs that were turned on and $28.9\pm 0.3$ °C under LEDs that were off (mean ± standard deviation of 6–10 wells).

Cells were grown in test tubes until OD 0.5–0.8, as described above, and diluted to OD 0.2 immediately prior to experiment. A black 96-well plate was filled (200 μL/well) with BG11M-agar (15 g/L agar) supplemented with sodium thiosulfate (1 mM) and HEPES (20 mM, pH 8.0). When the BG11M-agar mixture cooled to room temperature, 35 μL of cells from growing culture at OD 0.2 were added to each well of the plate. The plate was sealed with transparent UniSeal (GE HealthCare Life Sciences, Pittsburgh, PA, cat no. 7704–0001), holes were punched above each well using a 26G½ needle (BD, Franklin Lakes, NJ, cat. no. 305111), and the plate was placed underneath an LED array. Bioluminescence from luciferase reporters was measured every 30 min using a TopCount scintillation counter (PerkinElmer, Boson, MA). Each well of a 96-well plate was read for 1 s per measurement.

KaiA, KaiB, and KaiC were recombinantly expressed and purified as previously described (Lin et al., 2014), although the anion exchange chromatography purification step was omitted in the preparation of KaiC used for fluorescence polarization experiments. Protein concentration was measured via Bradford Assay Kit using bovine serum albumin (BSA) as a standard (BioRad, Hercules, CA). For experiments relying on SDS-PAGE analysis of KaiC phosphorylation, KaiABC proteins were mixed in master reaction buffer (20 mM Tris [pH 8], 150 mM NaCl, 5 mM MgCl₂, 10% glycerol, 0.5 mM EDTA) supplemented with a mixture of ATP and ADP (day buffer: 2 mM ATP, night buffer: 2 mM ATP, 7.5 mM ADP). All reactions were incubated at 31°C. To mimic light-to-dark transitions, ADP was added to appropriate reactions to 7.5 mM final [ADP]. To mimic dark-to-light transitions, the reactions were passed through Zeba desalting columns (7 MW kDa cutoff, ThermoFisher, Waltham, MA) equilibrated in day buffer. Because every buffer exchange step dilutes the proteins by about 10%, the reactions were prepared at 3× standard protein concentration (10.5 μM KaiB and KaiC, 4.5 μM KaiA). KaiC phosphorylation was assayed by SDS-PAGE and quantified by densitometry, as previously described (Phong et al., 2013).

In the cases where oscillations were read out by monitoring fluorescence polarization, KaiABC proteins were mixed at 3× standard protein concentration in the master reaction buffer supplemented with a mixture of ATP and ADP (day buffer: 2.5 mM ATP, night buffer: 2.5 mM ATP, 7.5 mM ADP). All reactions were incubated at 31°C. Day-to-night transitions were mimicked by addition of 7.5 mM ADP (final), and night-to-day transitions were mimicked by passing the reactions through Zeba desalting columns twice.

For step-up perturbations shown in blue in Figure 3D, buffer exchange steps were administered every 2 hr over a 24 hr interval, and phase shifts were measured relative to an unperturbed control reaction. For every other experiment in Figure 3C–D, buffer exchange or ADP addition steps were performed every 2 hr over a 12 hr interval on two out-of-phase reactions, prepared as follows. A master mix containing proteins and appropriate nucleotides was split into two tubes, which were flash-frozen in liquid nitrogen immediately after mixing and stored at −80°C. To prepare out-of-phase reactions, one of the two tubes was thawed in a 30°C water bath 12 hr later than the other.

After all buffer exchanges were completed, the reactions were supplemented with fluorescently labeled KaiB (0.2 μM final) and transferred to the plate reader.

We introduced a K25C mutation in KaiB using site-directed mutagenesis of the pMR0019 plasmid carrying a 6xHis-PSP-KaiB^WT construct in the pET47b(+) backbone. KaiB^K25C was expressed in BL-21 (DE3) E. coli by overnight induction with IPTG at 18°C and purified following the standard protocol (Lin et al., 2014). Labeling with 6-iodoacetofluorescein (6-IAF) was performed as previously described (Chang et al., 2012) with minor modifications. Briefly, 130 μL of KaiB^K25C stock (50–100 μM) was buffer-exchanged into labeling buffer (20 mM Tris, 1 mM TCEP, pH 7.0–7.5) using a Zeba desalting column (7 MW kDa cutoff). A freshly prepared solution of 6-IAF (Life Technologies Corp., Grand Island, NY) in DMSO was added to the protein solution in 10-fold molar excess, and the mixture was incubated overnight at 4°C in the dark. The labeling reaction was quenched by addition of DTT in approximately 10-fold molar excess relative to 6-IAF. Unincorporated dye was removed through three rounds of fivefold dilution in reaction buffer and subsequent concentration using a centrifugal filter (10 kDa cutoff, Amicon, EMD Millipore, Billerica, MA). Concentration of final fluorescein-labeled KaiB^K25C solution was determined by Bradford assay.

Oscillations in KaiABC reaction mixtures supplemented with fluorescently labeled KaiB^K25C (10.5 μM KaiB and KaiC, 4.5 μM KaiA, 0.2 μM KaiB-fluorescein) were monitored on an Infinite F500 plate reader equipped with a fluorescence polarization module (Tecan Trading AG, Switzerland). At least 30 min prior to measurement, the built-in heating module was turned on to warm the instrument and a black 384 well-plate (Greiner Bio-One, Monroe, NC, cat. no. 781900) was loaded into the plate reader. For the metabolic entrainment experiment shown in blue in Figure 2—figure supplement 1A and the step-up experiment shown in blue in Figure 3D, the plate reader temperature was 28°C; for all other experiments, the plate reader temperature was 31°C. Reactions were quickly transferred onto the pre-warmed plate (20–35 μL/well) and one to three wells were filled with master reaction buffer. The plate was sealed with a polyethylene silicone plate sealer (Nunc, ThermoFisher cat. no. 235307) and returned to the instrument. Fluorescence polarization of wells of interest (exc. 485 nm, 20 nm bandpass; em. 535 nm, 25 nm bandpass; dichroic 510 nm.) was read out every 15 min using a script created in iControl software (v. 1.12, Tecan Trading AG). Wells containing reaction buffer only were used as blanks, and the G factor was calibrated such that a solution of free fluorescein in reaction buffer produced a reading of ≈20 mP.

All oscillating trajectories were fit to sinusoids using optimization routines written in MATLAB (Mathworks, Inc.) (RRID: SCR_001622). See Computational methods for detailed fitting procedures and descriptions of simulations. Computational pipelines used for analysis and simulations have been deposited to GitHub at https://github.com/euleip/Simulations_and_analysis_pipelines_github_repo (Leypunskiy, 2017). A copy is archived at https://github.com/elifesciences-publications/Simulations_and_analysis_pipelines_github_repo).

All nonlinear fitting procedures were written in MATLAB using the lsqnonlin() or nlinfit() routines. Linear regressions were performed using the polyfit() function. Uncertainties around the best-fit slopes were evaluated using standard formulas for linear regression with known errors in dependent variables (Press et al., 1992). For in vivo experiments, these errors were computed as standard errors of the mean from replicate measurements; for in vitro measurements in Figure 2, the errors in fit phases were computed from the curvature of the cost function at the optimum using the nlpredci() function in MATLAB. Where indicated, 95% confidence intervals were also estimated by nlpredci().

Prior to fitting, bioluminescence trajectories were normalized to zero mean and unit standard deviation over the fitting interval, unless specified otherwise. In all analyses, we discarded data from the first 2.5 hr after lights-on to avoid masking effects.

The bacterial luciferase reporter system exhibits transient masking effects following dark-to-light transitions. To overcome this issue in determining the phase of gene expression in light-dark cycles of varied day length, we employed a drive-and-release strategy. In this approach (illustrated for LD 8:16, LD 12:12, and LD 16:8 in Figure 1A), cells were first entrained to a given diurnal schedule and then released into constant light for several days to determine the peak times of the entrained rhythm.

We relied on two approaches to calculate the peak times of normalized bioluminescence trajectories recorded after release into constant light: (i) we fit sinusoids to 48 hr segments of the trajectories in constant light, and (ii) we locally fit parabolas near the maxima of the trajectories. The advantage of sinusoidal fitting is that it captures phase information for the entire waveform; on the other hand, local parabolic fitting allows for a precise determination of the time of an individual peak without influence from others. In practice, we found that the estimates of the scaling between entrained phase and day length derived from the two approaches were in good agreement with each other (Table 1, Figure 1—figure supplement 2). Fitting details are described in the following two subsections.

Error bars in Figure 1C, Figure 1—figure supplement 2A, and Figure 5 represent the standard deviation of peak times calculated from technical replicates (n = 4–8). Technical replicates refer to measurements obtained from side-by-side cultures subjected to the same light-dark conditions. In rare cases, cells in individual wells died or produced noisy bioluminescent signals or trajectories that fit poorly to sinusoids or parabolas. We rejected trajectories as outliers from our analysis if they produced fits with a squared error greater than 10 (0–13 outliers per 96 wells).

To compare clock reporter dynamics during entrainment and in continuous light in Figure 1—figure supplement 1, we computed nonparametric correlations (Kendall’s τ coefficient) between the reporter signal (P_kaiBC::luxAB) measured during the ‘day’ windows in light-dark cycles (days 1–5) and during the corresponding time interval after release into continuous light (day 6). For example, in the case of LD 14:10, the normalized luminescence data recorded between 2.5 and 14 hr during each of the five entraining cycles were correlated with the luminescence dynamics measured between 2.5 and 14 hr after release into constant light. The fact that we observe Kendall’s τ > 0.8 after the second driving cycle suggests that cells are effectively entrained within three cycles and that rhythms observed during the lights-on portion of a light-dark cycle can be thought of as a fragment of a free-running rhythm.

For the %P-KaiC measurements in Figure 2C, normalized KaiC phosphorylation time courses from day 3 of each reaction were fit independently to sinusoids. Best-fit parameters were obtained by minimizing the cost function:

with the period T fixed at 24 hr to match the period of the imposed metabolic cycles.

For the fluorescence polarization experiments in Figure 2C, fluorescence polarization dynamics were recorded for 48 hr after the last buffer exchange, normalized (to zero mean, unit variance) and fit to sinusoids according to the expression above, with fits to all reactions from the same experiment sharing a globally fit period T.

Step-response experiments described in this section were performed a total of three times: once using KaiC phosphorylation to read out clock phase and twice using the fluorescence polarization reporter of KaiB-KaiC interaction. The phosphorylation measurements were made using an independent preparation of proteins from the fluorescence polarization experiments.

We performed step-response experiments in Figure 3 in order to determine whether the behavior of the clock in metabolic cycles could be decomposed into a sum of phase shifts due to individual transitions. To do so, we first needed (a) to extract L and D functions from step response measurements, and then (b) to use these L and D functions in numerical simulations of clock entrainment to light-dark cycles, as described in the main text and Appendix 1.

In our experiments, we directly observed how the dynamics of the fluorescence polarization reporter or KaiC phosphorylation changed as a function of lab time (e.g. Figure 3A–B), but for our downstream simulations and analysis (e.g. Figure 4) these values had to be converted to clock phase coordinates (e.g. Figure 3C–D) consistent with how ${\displaystyle L(\theta )}$ and ${\displaystyle D(\theta )}$ are defined. Specifically, we needed to (i) convert the lab time of each step to corresponding clock phase ${\displaystyle \theta }$ (or clock time, CT), and then (ii) determine the phase of each fluorescence polarization trace or KaiC phosphorylation trajectory after the perturbation in order to compute the phase shift (${\displaystyle L(\theta )}$ and ${\displaystyle D(\theta )}$) relative to an unperturbed control.

Although these conversions are straightforward in principle, it is important to note that both (i) and (ii) rely on sinusoidal regression of phase from measured data, and that the best-fit phases in both cases are only estimates of the true values. The corresponding uncertainty in the best-fit parameters must be propagated through numerical simulations of entrainment.

Because the temporal resolutions of the KaiC phosphorylation and fluorescence polarization measurements are very different (3 hr vs. 15 min, respectively), we expected the sources of uncertainty in fitting these data to be different (see Figure 2—figure supplement 1), and we thus analyzed their errors in different ways. The sparse measurement of KaiC phosphorylation dynamics leads to relatively high uncertainty in best-fit phase and period, and we propagated the errors through our simulations using nonparametric bootstrapping of the datasets. In the case of the fluorescence polarization measurements, experiment-to-experiment variability is the major source of uncertainty. We therefore performed two replicate measurements of L and D (Figure 3C–D) and used the four possible combinations of L and D in entrainment simulations to assess the range of entrainment behavior constrained by our measurements of these step-response measurements. We describe each of these approaches in turn in the following two sections.

Periods and phases of reactions in each set were determined by global sinusoidal fitting: amplitude, offset and phase terms were fit independently for each reaction, but a single best-fit period was shared among all fits in a given step-up or step-down set. To avoid transient effects due to a metabolic pulse, only those data points which were collected at least 16 hr after a step-up or step-down perturbation were used for fitting. Mathematically, we used a non-linear least-squares optimization routine to minimize the cost functions:

and

where ${\displaystyle y_{r,i}}$ refers to the i-th data point in the r-th normalized fluorescence polarization trajectory, ${\displaystyle t_{r,i}}$ refers to the lab time when that data point was collected, and ${\displaystyle N_r}$ is the number of data points fit in the r-th reaction. Here ${\displaystyle A_r}$, ${\displaystyle {\varphi }_r}$ and ${\displaystyle C_r}$ are the best-fit parameters for the r-th reaction; ${\displaystyle T_{day}}$ and ${\displaystyle T_{night}}$ are best-fit periods for all reactions in the 100% ATP and 25% ATP datasets, respectively.

These best-fit parameters were used to determine the phase ${\displaystyle \theta }$ of each step perturbation and corresponding phase shifts ${\displaystyle L(\theta )}$ and ${\displaystyle D(\theta )}$, analogously to the way described for the KaiC phosphorylation datasets below. The step-response functions were interpolated linearly to generate smooth curves while enforcing 2π-periodicity.

Linearized step-response functions ${\displaystyle L_{lin}(\theta )}$ and ${\displaystyle D_{lin}(\theta )}$ were prepared by linear regression of ${\displaystyle L(\theta )}$ and ${\displaystyle D(\theta )}$ centered on the regions used during metabolic entrainment: D functions were linearized between 6 and 22 CT hr; ${\displaystyle L}$ functions were linearized between 18 and 34 CT hr (see Figure 3C–D). To satisfy periodicity, the step-response functions were assembled as piecewise-linear functions with the same slope everywhere but with offsets every 2π radians at ‘breakpoints,’ which were selected by manual inspection. See Figure 4—figure supplement 2C for an illustration of ${\displaystyle L(\theta )}$ and ${\displaystyle D(\theta )}$ and the corresponding ${\displaystyle L_{lin}(\theta )}$ and ${\displaystyle D_{lin}(\theta )}$.

Finally, step phases deduced from polarization trajectories were adjusted to match the phases of the KaiC phosphorylation rhythm. This conversion made use of measurements in Figure 2—figure supplement 1, which indicated that the phase of oscillation in KaiC phosphorylation lags behind the phase of the polarization reporter of KaiB-KaiC binding by approximately ${\displaystyle 2\pi /3}$.

By performing this analysis for the two step-up and two step-down datasets in Figure 3C–D, we generated two sets of {${\displaystyle L}$, ${\displaystyle L_{lin}}$, ${\displaystyle T_{day}}$, ${\displaystyle T_{night}}$} and two sets of {${\displaystyle D}$, ${\displaystyle D_{lin}}$, ${\displaystyle T_{day}}$, ${\displaystyle T_{night}}$}. The phase oscillator simulations described below require using the ${\displaystyle L}$ and ${\displaystyle D}$ functions together in each simulation. To propagate the experiment-to-experiment variability in step-response measurements through our simulations, we combined the two measurements of ${\displaystyle L}$ and two measurements of ${\displaystyle D}$ into the four possible combinations of {${\displaystyle L}$, ${\displaystyle D}$} pairs. This resulted in four sets of {${\displaystyle L}$, ${\displaystyle L_{lin}}$, ${\displaystyle D}$, ${\displaystyle D_{lin}}$, ${\displaystyle T_{day}}$, ${\displaystyle T_{night}}$}, which were used in the simulations below.

The ratio of oscillator frequencies in dark and light, ${\displaystyle \frac{{\omega }_D}{{\omega }_L}=\frac{T_{day}}{T_{night}}=0.93\pm 0.01}$, was determined by averaging the values of T_day and T_night from the two step-up and two step-down sets above. Best estimates for the slopes l and d in Figure 4—figure supplement 3 were computed via the following bootstrap analysis. For each L and D function in Figure 3C–D, we selected points from the regions used for linearization above (6–22 CT hr on D, 18–34 CT hr on L). We sampled these data points with replacement until we generated 500 samples containing at least three unique points. For each set of resampled points, we computed the slope of the best-fit line, for a total of 1000 samples of l and d. The crosshair in Figure 4—figure supplement 3 marks the average ± standard deviation of these values (l = 0.34 ± 0.03, d = 0.38 ± 0.05).

The phase oscillator model discussed in the main text and the Appendix 1 makes the prediction that the proportionality constant m between entrained phase and day length depends on l and d coefficients according to:

To check whether this prediction is in line with our experimental measurement of ${\displaystyle m}$, we used this formula to calculate ${\displaystyle m}$ for every pair of ${\displaystyle l}$ and ${\displaystyle d}$ samples generated in the bootstrapping procedure above, assuming ${\displaystyle {\omega }_D/{\omega }_L=0.93}$. According to this calculation, $m\left(l,d\right)=0.34\pm 0.04$ (mean ± standard deviation of the distribution), which is in agreement with the experimental measurement in Figure 2C ($m=0.38\pm 0.07$).

Recall that measuring each step-response function requires (a) conversion of the lab time of each step to corresponding clock phase ${\displaystyle \theta }$ (or clock time, CT), and then (b) determination of the phase of each KaiC phosphorylation trajectory after the perturbation in order to compute the phase shift (${\displaystyle L(\theta )}$ and ${\displaystyle D(\theta )}$) relative to an unperturbed control. In particular, uncertainties in (a) manifest as phase errors on ${\displaystyle L}$ and ${\displaystyle D}$ in Figure 3E–F; these errors are correlated across all points on the step-response function. Uncertainties in (b) manifest as phase shift errors on ${\displaystyle L}$ and ${\displaystyle D}$; these errors derive from the errors in the phase estimates of both the control and step reactions. To propagate both these sources of error, we used the following non-parametric bootstrapping strategy.

First, KaiC phosphorylation dynamics from all step-response measurements were assembled into a master dataset containing nine step-up trajectories, nine step-down trajectories, as well as two control reactions. These measurements were performed once on an independent preparation of proteins from the batch used to generate data in Figure 3C–D. This master dataset was then trimmed to include only data collected at least 16 hr after a step transition, and every trajectory was normalized. To generate bootstrapped datasets, we sampled with replacement from the entire master dataset (as opposed to resampling reactions individually) 1000 times.

Next, we used each resampled dataset to compute phase shifts in KaiC phosphorylation due to step-up and step-down transitions; in other words, each bootstrapped dataset was used to derive a bootstrapped pair of ${\displaystyle L(\theta )}$ and ${\displaystyle D(\theta )}$. For a given dataset, we globally fit all 100% ATP trajectories (nine step-up trajectories, plus the 100% ATP control reaction) such that all fits shared a best-fit period (${\displaystyle T_{day}}$), but phase, amplitude and offset terms were fit independently for each reaction, as described above for polarization datasets. Likewise, we fit all 20% ATP trajectories (nine step-down reactions, plus the 20% ATP control reaction) to obtain a globally best-fit period (${\displaystyle T_{night}}$) and independently fit phases for every reaction.

${\displaystyle L(\theta )}$ and ${\displaystyle D(\theta )}$ map the phase at which a metabolic step occurs to the resulting phase shift. The phase at which the metabolic step occurred for every reaction r was computed from the best-fit phase of the appropriate control reaction at step time ${\displaystyle t_{r,step}}$:

where ${\displaystyle {\varphi }_{con}}$ is the best-fit phase of the control reaction and ${\displaystyle T_{con}}$ is the globally-fit oscillator period in the appropriate control condition (i.e. ${\displaystyle T_{con}}$ = ${\displaystyle T_{day}}$ for step-down reactions, ${\displaystyle T_{con}}$ = ${\displaystyle T_{night}}$ for step-up reactions).

Similarly, we computed the apparent phase of each perturbed reaction at the time of each step:

where ${\displaystyle {\varphi }_r}$ is the best-fit phase of the r-th reaction and ${\displaystyle T_{pert}}$ is the globally-fit oscillator period in the perturbed condition (${\displaystyle T_{pert}}$ = ${\displaystyle T_{night}}$ for step-down reactions, ${\displaystyle T_{pert}}$ = ${\displaystyle T_{day}}$ for step-up reactions).

Finally, we defined the phase shift in response to each step perturbation as the difference in phase of the perturbed reaction and the appropriate unperturbed control.

To estimate phase shifts at other values of ${\displaystyle \theta }$, we linearly interpolated L and D between the measured values while enforcing 2π-periodicity.

For each set of L and D generated in this way, we also prepared linearized versions L_lin and D_lin. Linear fits to L and D were performed over the range of step times similar to the ones we selected for fluorescence polarization-based step-response functions as discussed above (L between 18 and 33 CT hr, D between 6 and 23 CT hr). These regions of L and D were selected by visual inspection; they are centered on phases used by the KaiABC oscillator in diurnal cycles (i.e. see arrows in Figure 3E–F) but also contain enough points (5-7) to avoid biasing the slope estimate by a single poorly fit data point.

We extrapolated the linear approximations to L and D over the entire cycle, with a single breakpoint away from the linear region to satisfy 2π-periodicity. Breakpoints were selected by visual inspection. While interpolating between data points near the breakpoint, we assumed that L and D (or L_lin and D_lin) never generate phase shifts larger than one cycle (i.e. winding number of 0). We anticipate that this choice does not significantly affect entrained phase in most of our simulations because the regions of L and D near the breakpoint are rarely used by the oscillator during entrainment when τ <14 hr.

We repeated this procedure for each of the 1000 bootstrapped datasets, thereby generating 1000 sets of {${\displaystyle L}$, ${\displaystyle D}$, ${\displaystyle L_{lin}}$ and ${\displaystyle D_{lin}}$, ${\displaystyle T_{day}}$, ${\displaystyle T_{night}}$} that were used for subsequent analysis. Figure 3E–F shows the mean ± standard deviation of the distribution of bootstrapped L and D generated this way.

We simulated entrainment to a step-like driving cycle for a phase oscillator governed by the four combinations of {${\displaystyle L}$, ${\displaystyle L_{lin}}$, ${\displaystyle D}$, ${\displaystyle D_{lin}}$, ${\displaystyle T_{day}}$, ${\displaystyle T_{night}}$} determined from the fluorescence polarization assays and the 1000 bootstrapped sets of {${\displaystyle L}$, ${\displaystyle D}$, ${\displaystyle L_{lin}}$, ${\displaystyle D_{lin}}$, ${\displaystyle T_{day}}$, ${\displaystyle T_{night}}$} derived from KaiC phosphorylation measurements. Electing to work in units of cycles (1 cycle = 2π rad), we defined a phase variable ${\displaystyle \hat{\theta }=(\theta -\pi /2)/(2\pi )}$, such that ${\displaystyle \hat{\theta }=0}$ cycles and ${\displaystyle \hat{\theta }=1}$ cycles correspond to the trough of the KaiC phosphorylation trajectory and ${\displaystyle \hat{\theta }=0.5}$ cycles corresponds to the peak. Where circadian time (CT) is mentioned in the text, we have adopted the convention that CT = 0 and 24 hr refer to ${\displaystyle \hat{\theta }=0}$ and ${\displaystyle 1}$ cycles, respectively; CT = 12 hr refers to ${\displaystyle \hat{\theta }=0.5}$ cycles. We also defined corresponding analogs of ${\displaystyle L}$ and ${\displaystyle D}$ in units of cycles:

See Figure 4—figure supplement 2C for examples of ${\displaystyle \hat{L},\hat{D},{\hat{L}}_{lin}}$ and ${\displaystyle {\hat{D}}_{lin}}$. For each set of ${\displaystyle \left\{\hat{L},\hat{D},{\hat{L}}_{lin},{\hat{D}}_{lin},T_{day},T_{night}\right\}}$ generated in this way, we modeled oscillator entrainment to a driving cycle. As shown schematically in Figure 4A and the Appendix 1, oscillator phase increases with constant angular velocity (given by ${\omega }_L=\frac{1}{T_{day}}$ during the day and ${\displaystyle {\omega }_{\text{D}}=\frac{1}{T_{night}}}$ during the night); at dawn and dusk, oscillator phase shifts instantaneously according to ${\displaystyle \hat{L}}$ and ${\displaystyle \hat{D}}$, respectively. Therefore, oscillator phases at dusk and dawn of the n-th entraining cycle can be computed iteratively:

We simulated oscillator behavior in 24 hr light-dark cycles (driving period T_dr = 24) with day length τ lasting from 4 to 18 hr. In simulations relying on ${\displaystyle \hat{L}}$ and ${\displaystyle \hat{D}}$ determined from KaiC phosphorylation measurements, we set the initial condition ${\displaystyle {\hat{\theta }}_0}$ for a given value of τ based on our measurements of entrained phase of the KaiABC oscillator in corresponding metabolic cycles (i.e. based on values interpolated between black markers in Figure 2C). Simulations using step-response functions based on fluorescence polarization measurements were started from ${\displaystyle {\hat{\theta }}_0=0}$. We simulated 10 light-dark cycles and recorded oscillator phase after each dawn (immediately after ${\displaystyle \hat{L}\left(\hat{\theta }\right)}$ phase shift).

For each set of ${\displaystyle \hat{L}}$ and ${\displaystyle \hat{D}}$ we performed simulations for τ = 4–18 hr for both the phosphorylation-based step functions and fluorescence polarization-based step functions. Figure 4—figure supplement 1 illustrates simulations governed by one such set of ${\displaystyle \hat{L}}$ and ${\displaystyle \hat{D}}$ (derived from the magenta and blue step-response measurements shown in Figure 3C–D). We also carried out entrainment simulations for the linearized step-response functions described above (Figure 4—figure supplement 2).

We judged that the oscillator entrained stably to a given diurnal cycle if the standard deviation of clock phases at the dawns of the last five driving cycles was less than 0.01 cycles. For each value of τ, we selected those simulations where the oscillator entrained stably to the light-dark cycle, and computed the mean phase at dawn and its standard deviation (σ) for that set of simulations. These values are represented in shaded areas in Figure 4C , Figure 4—figure supplement 2A (for fluorescence polarization data) and Figure 4—figure supplement 2B (for KaiC phosphorylation data).

We were also interested in how quickly the oscillator approached entrainment in simulations for 6 ≤ τ ≤18, the range we profiled experimentally in Figure 2. To make this determination for a given value of τ, we computed how much oscillator phase at dawn varied over three successive cycles in a sliding window:

We used the first value of n for which ${\displaystyle E{C}_n<0.01}$ as a proxy for the speed of approach to entrainment.

The tables below display summary statistics for the entrainment simulations.

In the large majority of our simulations, we found that the phase oscillator entrained stably within three light-dark cycles. For simulations derived from KaiC phosphorylation datasets, we found that the oscillator either entrained to the driving cycle quickly or not at all. Indeed, when we restricted our analysis only to those simulations that were judged as entrained within eight cycles, over 96% entrained within three light-dark cycles. Generally, the oscillator entrained readily for day lengths shorter than 14 hr, but often failed to entrain for longer day lengths. We determined that this occurs because for τ >14 dawn phases often sample the $\widehat{L}$ function near the breakpoint of the curve (near 18 CT hr in Figure 3F), leading to erratic responses to driving cues (i.e. lack of entrainment) or disagreement between simulations and experiment (Figure 4—figure supplement 2B). Relatedly, we believe that the better agreement with experiment achieved in simulations using the step functions derived from the fluorescence polarization data than from the KaiC phosphorylation data reflects the better temporal resolution of the breakpoint (2 hr using the polarization approach vs. 4 hr for KaiC phosphorylation).

In simulations of entrainment to light-dark cycles of varying period in Figure 3—figure supplement 2, we relied on a single set of {${\displaystyle L}$, ${\displaystyle D}$, ${\displaystyle L_{lin}}$, ${\displaystyle D_{lin}}$, ${\displaystyle T_{day}}$, ${\displaystyle T_{night}}$} based on the step-response measurements shown in blue and magenta in Figure 3C–D. We simulated entrainment to driving periods from 6 to 48 hr, in increments of 0.0025 hr. For each driving period T_drive, we subjected the phase oscillator to 1000 cycles with equal day and night durations (τ − T_drive/2) and plotted the phases attained by the oscillator at the end of nighttime (immediately preceding the action of L) at the last 950 cycles.

We used ${\displaystyle {\hat{L}}_{lin}\text{and}{\hat{D}}_{lin}}$ derived from the same step-response measurements as in Figure 4—figure supplement 1 to simulate response of a phase oscillator to 12 hr dark pulses administered throughout the circadian cycle (Figure 5—figure supplement 1). Phase evolution of the oscillator was simulated explicitly for 120 hr using a timestep of dt = 0.01 hr:

For phase resetting and wedge experiments (Figure 5), clock phase was estimated by sinusoidal regression of normalized bioluminescence data collected 36–48 hr after the end of the applied dark pulse. Phase shifts in response to dark pulses were computed as differences in average peak times between perturbed wells (${\displaystyle t_{pk,DP}}$) and controls (${\displaystyle t_{pk,LL}}$): ${\displaystyle \mathrm{\Delta }{t}_{pk}={\overline{t}}_{pk,DP}-{\overline{t}}_{pk,LL},}$ where overbars indicate averages over replicate wells (n = 4–8). Clock phases at which the dark pulses were applied were determined from the average fit phase and period of the unperturbed (control) wells (${\overline{\varphi }}_{LL}$ and ${\overline{T}}_{LL})$) according to ${\displaystyle {\hat{\theta }}_t=(t/{\overline{T}}_{LL})-}$${\displaystyle ({\overline{\varphi }}_{LL}+0.5\pi )/(2\pi )}$, where ${\displaystyle {\hat{\theta }}_t}$ is measured in cycles and ${\displaystyle \hat{\theta }=0}$ corresponds to the minimum of an oscillatory trajectory.

In the Appendix 1, we show that in the regime where the circadian clock is well-approximated by a phase oscillator governed by linear L and D step-response functions, the slopes of seasonal entrainment and phase resetting of the clock can be described by a model with two free parameters ${\beta }_1$ and ${\beta }_2$: for phase resetting, ${\displaystyle \mathrm{\Delta }{t}_{pk}={\hat{\theta }}_t/({\omega }_L{\beta }_2)+\delta \left(1+{\beta }_1/{\beta }_2\right)+C_1}$ and for entrainment $t_{pk}=\tau \left(1-{\beta }_1-{\beta }_2\right)+T{\beta }_1+C_2$. Here, δ is dark pulse duration (in hours), τ is day length (in hours), ${\displaystyle {\hat{\theta }}_t}$ is the clock phase at time t (in cycles), ${\omega }_L$ is the clock frequency in the light (in units of cycles/hour), T is the driving period (in hours), and $C_1$ and $C_2$ are constants that do not depend on ${\displaystyle {\hat{\theta }}_t}$, $\delta$, $\tau$, or $T$. For phase resetting experiments, ${\displaystyle {\hat{\theta }}_t}$ and ${\omega }_L$ were estimated based on the average of the best sinusoidal fits to unperturbed wells in each experiment (${\overline{\varphi }}_{LL}$ and ${\displaystyle {\overline{T}}_{LL}}$). In particular, we set ${\displaystyle {\hat{\theta }}_t={\omega }_{L}t-}$${\displaystyle ({\overline{\varphi }}_{LL}+0.5\pi )/(2\pi )}$. In the formula for ${\displaystyle \mathrm{\Delta }{t}_{pk},}$ the term ${\displaystyle {\hat{\theta }}_t/({\omega }_L{\beta }_2)}$ thus simplifies to ${\displaystyle \left(t/{\beta }_2\right)-({\overline{\varphi }}_{LL}+0.5\pi )/(2\pi {\beta }_2)}$, and the term to the right of the minus sign was incorporated into the constant term C₁.

In the global fits of all datasets in Figure 5A–C, we varied ${\beta }_1$ and ${\beta }_2$ to simultaneously fit ${\displaystyle \mathrm{\Delta }{t}_{pk}}$ to our phase-resetting and wedge data and $t_{pk}$ to our seasonal entrainment data. The constant terms were allowed to vary as follows:

• in Figure 5C, a single C₂ term was fit for all curves, referred to as C_entrainment below;

• in Figure 5B, a single C₁ intercept was fit for both datasets (t_DP = 36 and t_DP = 39 hr), referred to as C_wedge below;

• in Figure 5A, the points before and after the breakpoint were fit using the same slope, but varying constant (C₁) terms. The breakpoint was selected by visual inspection. Below, the intercepts to the left and right of the breakpoint are referred to as C_{PRC_left} and C_{PRC_right}, respectively.

In total, only two parameters (${\beta }_1$ and ${\beta }_2)$) were used to determine the slopes of all curves, and six parameters were used for the entire global fit of nine linear segments, which minimized the cost function:

where each ${\displaystyle {\hat{y}}_i}$ represents the average measurement of phase shift or peak time and ${\displaystyle {\sigma }_i}$ represents the standard error of that measurement. The reduced chi-squared value of the fit was ${\chi }_{\nu }^2=5.67$. The best-fit coefficients determined in the fit are summarized in the table below.

To simulate entrainment to light-dark cycles of different day length in the geometric resetting framework in Figure 6, we modeled an oscillator running along the daytime limit cycle (centered at 0, radius 1) in the light and the nighttime cycle (centered at X, radius R) in the dark. We set the angular frequency to be ${\displaystyle \frac{2\pi \mathrm{r}\mathrm{a}\mathrm{d}}{24\mathrm{h}\mathrm{r}}}$ in both light and dark. As described in the Appendix 1, dusk and dawn transitions were modeled as radial jumps from one cycle to the nearest point on the other cycle.

We considered values of R and X spanning six orders of magnitude and studied entrainment to 24 hr cycles with day length lasting from 6 to 18 hr. For each pair of R and X, we simulated oscillator dynamics in 30 light-dark cycles of a given day length. We judged that an oscillator failed to entrain to a given diurnal schedule if oscillator phases after 29 and 30 cycles were more than π/180 radians apart, or if simulations starting from different initial phases (θ₀ = π/4 and 5π/4) reached phases over π/180 radians apart after 30 light-dark cycles.

For every simulation that passed the entrainment criteria above, we computed the best linear fit and slope m of oscillator phase dependence on day length τ. We assessed goodness of fit by computing the mean fit error, defined as the average absolute value of the deviations between the linear fit and simulation results. If the mean fit error was greater than 10% of the deviation between maximum and minimum phases to which the oscillator entrained for this range of τ, the phase dependence on day length was judged to be non-linear.

The results of the simulations described immediately above strongly suggest that the relative geometry of day and night orbits determines the scaling of entrained phase with day length. However, those simulations are based on idealized infinitely-attracting circular limit cycles with constant angular frequency, assumptions which are likely to be violated for biological clocks. For example, the KaiC phosphorylation limit cycles in Figure 6F are somewhat elliptical. We were therefore interested in understanding whether the relative geometry of the limit cycles would be the dominant determinant of oscillator entrainment if our assumptions were relaxed. To this end, we explored entrainment in limit cycle models where these features—(i) orbit attraction strength, (ii) orbit ellipticity, and (iii) variation in angular frequency throughout the cycle—could be treated explicitly.

Phase of transcriptional dynamics of dawn and dusk genes in Figure 6—figure supplement 3 were derived from a sinusoidal fit to the transcriptional profiles of dawn and dusk genes identified in microarray time courses by Vijayan et al. (2009). Briefly, normalized time course data was downloaded from the GEO depository, and time series for transcripts annotated as dawn genes were averaged to obtain the average dawn gene transcriptional profile. The average dusk gene trajectory was obtained analogously. Average dawn and dusk waveforms were fit to sinusoids with 24 hr period, and the resulting phase estimates were used to define the phase shift between the orange and gray curves in Figure 6—figure supplement 3.

In Figure 6—figure supplement 3, nightfall in LD 12:12 for all values of m coincides with the 12 hr dark pulses administered during the initial synchronization in the Vijayan et al. experiment. In simulations of other day lengths (τ), we assumed that the only effect of diurnal cycling is to adjust the phase of the circadian transcriptional program relative to dawn and dusk, without affecting the shape or relative timing of dawn and dusk gene transcriptional waveforms. The bias in allocation of daytime resources between dawn and dusk genes was then computed according to the expression:

where $I\left(gene\right)$ is the integrated RNA signal for that gene over the daytime hours.

Circadian rhythms are characterized by the ability to maintain self-sustaining ≈24 hr oscillations in constant conditions and to adjust the phase of those oscillations in response to environmental signals (e.g. transitions between light and dark). To represent the first feature, we model a circadian clock as an oscillator that runs along a ‘light’ limit cycle with speed ${\omega }_L$ when lights are on, and along a ‘dark’ limit cycle with speed ${\omega }_D$ when lights are off (Appendix 1—figure 1A). Because the oscillator is always on one of the limit cycles, its state can be fully described by a single variable, the phase $\theta \left(t\right)$:

We define frequencies in units of cycles per hour (${\displaystyle 1\mathrm{c}\mathrm{y}\mathrm{c}\mathrm{l}\mathrm{e}=2\pi \mathrm{r}\mathrm{a}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}\mathrm{s}}$), so that the duration of one cycle in constant conditions is $T=1/\omega$ hours, and phase $\theta \left(t\right)$ is defined between 0 and 1.

We assume that the oscillator transitions between the two limit cycles immediately at dark-to-light (‘dawn’) and light-to-dark (‘dusk’) transitions, with a change in phase. We define $L\left(\theta \right)$ and $D\left(\theta \right)$ as phase-shifts accompanying clock responses to lights on and lights off cues, respectively:

Though in principle $L\left(\theta \right)$ and $D\left(\theta \right)$ could be arbitrary periodic functions, our measurements in Figure 3 suggest that $L\left(\theta \right)$ and $D\left(\theta \right)$ are approximately linear over the range of clock phases when dawn and dusk occurred in our experiments. Henceforth we consequently assume that ${\displaystyle L}$ and ${\displaystyle D}$ are linear functions of ${\displaystyle \theta }$:

Given this framework, we can now determine the change in clock state over any interval simply by adding the changes in phase at the transitions and accumulated during the free runs. We now explicitly compute the change in clock state in response to a dark pulse and the state of stable entrainment for light-dark cycles of arbitrary period.

To probe the response of an oscillator, one can subject it to a perturbation of duration $\delta$ and compare the resulting phase to that of an unperturbed oscillator. Plotting the results of such an assay as a function of the phase ${\theta }_t$ at which the perturbation is delivered yields a phase response curve (PRC; Appendix 1—figure 1C, Figure 5A). Alternatively, varying the duration of the perturbation while delivering it at a fixed phase corresponds to a ‘wedge’ experiment (Figure 5B).

To determine the phase response, we need to determine the phase evolution as the system transitions to dark, free runs in the dark, and then transitions to light. Note that because this model treats the transitions between the light and dark cycles as instantaneous, we expect that it will only be valid for dark pulses that are long compared to the true relaxation time of the clock. Mathematically,

where the prime denotes the perturbed oscillator, and the + and – in the subscripts distinguish the phase after and before the transition from dark-to-light at $t+\delta$. That is,

Substituting the linear forms above for $D$ and $L$ and grouping like terms,

Subtracting the unperturbed phase ${\theta }_{t+\delta }={\theta }_t+\delta {\omega }_L$,

Thus we obtain a formula for the change in phase in response to a dark pulse, with linear dependences on the phase of the perturbation, ${\theta }_t$, and its duration, $\delta$.

We can similarly describe the phase of a clock stably entrained to light-dark cycles of varying day length. To this end, consider a phase oscillator subjected to light-dark cycles of period $T$ and day length $\tau$ (Appendix 1—figure 1D). In this case, rather than a phase difference, we need to determine ${\theta }_{\tau }$ such that ${\theta }_{T+\tau }={\theta }_{\tau }+1$. If this condition is met, the oscillator is stably entrained to a cycle of period T. In other words, we enforce periodicity.

Proceeding analogously to above, we determine the evolution over a full driving period, starting at dusk. This comprises a transition to the dark at time $\tau$, free run in the dark for duration $T-\tau$, a transition to the light at time $T$, and free run in the light for duration $\tau$. Mathematically,

Again, substituting the linear forms above for $D$ and $L$ and grouping like terms,

Setting ${\theta }_{T+\tau }={\theta }_{\tau }+1$ and solving for ${\theta }_{\tau }$,

Again, we obtain a formula with linear dependence, here on the day length $\tau$ and the driving period $T$.

The expressions above describe phases, but we do not measure phase directly. Rather, we measure the time of peak expression of a reporter. Thus it is necessary to convert from phase to the time of peak expression. To this end, we arbitrarily assign the peak expression to occur at $\theta =0.5$ (midway through the cycle). Because phase advance and delay correspond to earlier and later peak times respectively, the time of the peak shifts with opposite sign to the phase. In each case above, we are interested when the peak (or $\theta =0.5$) occurs relative to the last dawn (i.e. immediately after the action of L). We can solve for this peak time:

where $x=t+\delta$ (end of dark pulse) in the case of phase resetting and $x=0$ (last dawn) in the case of entrainment. To solve for ${\displaystyle t_{pk}}$ in entrained conditions, we make use of the fact that ${\theta }_0={\theta }_{\tau }-{\omega }_L\tau$, allowing us to substitute the expression above for ${\theta }_{\tau }$. For phase resetting, we must apply this formula to both the perturbed and unperturbed phase and then subtract:

Furthermore, we note that certain combinations of the frequencies ${\omega }_L$ and ${\omega }_D$ and the linear sensitivities $l$ and $d$ occur repeatedly. We thus define

Then, for phase resetting

and for entrainment

where $C_1$ and $C_2$ are constants that do not depend on ${\theta }_t$, $\delta$, $\tau$, or $T$.

In the global fits for all datasets in Figure 5A–C, we varied ${\beta }_1$ and ${\beta }_2$ to simultaneously fit ${\displaystyle t_{pk}}$ to our phase-resetting and wedge data and ${\displaystyle t_{pk}}$ to our seasonal entrainment data. Thus, the slopes of all the linear fits in Figure 5A–C were derived from the same values of ${\beta }_1$ and ${\beta }_2$. See Computational Methods for further details.

In our experimental measurements, we observed L and D step-response functions with linear regions with slopes < 1. Here we present a simple dynamical systems picture where a tunable, linear dependence of phase shift magnitude on oscillator phase arises from the geometry of how a strongly attracting limit cycle is deformed by a changing external input.

We suppose that the oscillator is described by a circular limit cycle in the plane with unit radius during the day. During the night, we assume that the limit cycle remains in the plane, but is offset a distance X from the daytime cycle and has an altered radius R (Figure 6A). At the moment of the light-to-dark transition, the state of the system is on the daytime limit cycle—it is then attracted to the nighttime limit cycle. The step-response function ${\displaystyle D(\theta )}$ describes the shift between the new phase angle (measured relative to the center of the nighttime limit cycle) and the original phase ${\displaystyle \theta }$ (Figure 6A).

If we place the origin at the center of the nighttime limit cycle, the Cartesian coordinates of the system state immediately before the light-to-dark transition, when the system is still on the daytime orbit, are:

If we assume that the fixed point giving rise to the nighttime orbit is very strongly attracting, then the system approaches the limit cycle much faster than any circulation occurs. In this limit, the light-to-dark transition leads to a purely radial jump from the daytime orbit towards the center of the nighttime orbit (Figure 6A–B). Under these conditions, the new phase ${\displaystyle {\theta }_{night}}$ on the nighttime cycle is simply the angle measured relative to the center of the new limit cycle:

To develop intuition about how this mapping between phase on day and night orbits is affected by their geometric arrangement, we can perform a Taylor expansion for angles near ${\displaystyle \theta }$ = 0 at dusk, corresponding to conditions close to LD 12:12:

Likewise, it can be shown that for transitions from nighttime to daytime cycles, near θ = π , the phase immediately after transition to the daytime cycle at dawn can be given by:

These expressions show that, in the linear portion of the step response function, angles are compressed by a factor tunable by the relative geometry of the light and dark limit cycles. In this way, the geometric arrangement of the limit cycles sets the slopes of ${\displaystyle L(\theta )}$ and ${\displaystyle D(\theta )}$. Numerical simulations in Figure 6 and Figure 6—figure supplement 1 were carried out using trigonometric formulas as above to convert phases between the daytime and nighttime limit cycles. Details of simulations in Figure 6—figure supplement 2 are described in Computational Methods.