Many organisms like humans and rodents have free running clocks that show self-sustained 24 hr rhythms even in constant dark or light conditions. A particularly simple and well-characterized free running clock is that found in S. elongatus where the clock dynamics can be reproduced by the post-translational dynamics of Kai ABC in vivo as well. Measuring the phosphorylation state at any one of several sites on KaiC reveals an orderly phosphorylation reaction with a period of $24$ hr. As shown in Figure 1a, oscillations of a characteristic amplitude are sustained even in constant darkness or constant light, that is, in the absence of a periodic external drive.

Figure 1

Download asset Open asset

Not all organisms have a free-running clock; for example, many insects (Saunders, 2002) have damped ‘hourglass’ clocks that decay to a fixed point under constant light or constant dark conditions but show oscillatory dynamics under day-night cycling (see Figure 1b). In fact, a sister cyanobacterial species P. marinus has a KaiBC-protein based clock. While the details of this clock are not fully characterized, the clock lacks the KaiA-mediated negative feedback (Dufresne et al., 2003; Holtzendorff et al., 2008) loop that enables free running oscillations in S. elongatus. Consequently, in constant light or dark conditions, the clock’s state decays to a distinct day or a night state respectively (Holtzendorff et al., 2008).

Thus, both classes of clock show regular oscillations when externally driven. With cloudless day-night cycling, both systems can synchronize in phase with the external signal (i.e., ‘entrain’) and show distinct clock states at distinct times of the day. In this way, the clock state provides the rest of the cell with an estimate of the time of the day. However, while the free running clock has a natural amplitude relatively independent of the external drive, the damped clock’s amplitude is directly set by the external drive.

The day-night pattern of light on earth does not resemble the clean square wave shown in Figure 1a but is rather subject to large amplitude fluctuations during the day due to weather patterns. Such amplitude fluctuations and their spectrum have been quantified (Gu et al., 2001) and also identified as playing a critical role in several studies on the evolution and performance of circadian clocks (Domijan and Rand, 2011; Troein et al., 2009). The biological impact of such changes in light intensity in cyanobacteria have been quantified recently (Teng et al., 2013). The clock must entrain in phase to the external signal while ignoring such amplitude fluctuations.

In addition to external fluctuations, circadian clocks also deal with the intrinsically noisy nature of biochemical reactions (Swain et al., 2002). Sources of internal noise for clocks include finite copy number effects, bursty transcription, cell division and other sources (Tsimring, 2014). In particular, based on their relative sizes (Dufresne et al., 2003; Holtzendorff et al., 2008; Bryant, 2003), P. marinus is thought to have far fewer copies of the Kai clock proteins (e.g., $\sim 500$ of KaiC ) than S. elongatus ($\sim O(10000)$ copies of KaiC [Gutu et al., 2013; Kitayama et al., 2003]). Such finite numbers of molecules is known to create significant stochasticity in oscillators in the absence of an external signal (Potoyan and Wolynes, 2014).

We tested the impact of such external and internal fluctuations on the contrasting clock architectures in S. elongatus and P. marinus through simulations. We set up explicit Gillespie simulations (Gillespie, 2007) of explicit biomolecular models of the post-translational Kai clock that captures the known biochemistry (Rust et al., 2007) of S. elongatus’s clock and the putative KaiBC clock (Bryant, 2003; Holtzendorff et al., 2008) in P. marinus (Figure 1). We do not include transcriptional coupling (Zwicker et al., 2010) of the clock here and focus on the core post-translational oscillator. See Appendix 1 for details. The ATP levels in these models (Pattanayak et al., 2014) were coupled to an external square wave input of period 24 hr, representing the day-night cycle of light. To model external fluctuations, we modulated the amplitude of the input square wave over a broad range of frequencies, reflecting the broad frequency spectrum quantified by the Harvard Forest database (Moore et al., 1996). To model internal fluctuations, we varied the copy number in these Gillespie simulations.

With only external fluctuations but suppressing internal fluctuations using high copy numbers, we find that the damped oscillator develops a much larger population variance than the free running clock. In contrast, at low copy number (i.e., high internal noise) but with a noiseless external signal, we find the situation is reversed; the free running clock has significantly higher population variance. See Figure 1c.

To study this effect quantitatively, we fixed the strength of amplitude fluctuations and increased the internal noise by decreasing the copy number of all Kai proteins in the Gillespie simulation. We measured the resulting mutual information between clock state and objective time of day. (Mutual information is intuitively a measure of population variance along the most informative directions; see Appendix 4 for more.)

We see that the free running clock has higher precision than the damped clock at low internal noise (high copy number). However, as the internal noise is increased, the precision of the free running clock drops earlier and consequently, the damped oscillator has higher precision at sufficiently high internal noise (low copy number). This is shown in Figure 1d, where the precision measures the mutual information between the clock state and the time. For a fair comparison, in undriven conditions, different clocks are assumed to lose information at the same rate.

While our study here was motivated by the contrasting Kai protein-based clocks in the two cyanobacterial species S. elongatus and P. marinus, we sought to test the broader validity of our results. Hence we analyzed the internal and external noise resistance in a range of eight well-studied biochemical oscillators in the literature.

These models range from circadian clocks in numerous organisms - Neurospora (Leloup et al., 1999), Arabidopsis (Schmal et al., 2014; Locke et al., 2005), mammalian cells (Leloup and Goldbeter, 2003) - to other oscillators such as cell cycle models (Goldbeter, 1991), the Goodwin (Goodwin, 1965; Gonze and Abou-Jaoudé, 2013) oscillator, the Brusselator (Kondepudi and Prigogine, 2014) and the synthetic repressilator (Elowitz and Leibler, 2000; Buşe et al., 2009) - see Figure 2. While the internal noise properties of these oscillators in undriven conditions have been studied before (Gonze et al., 2002), here we analyzed the contrasting internal and external noise properties of these oscillators under externally driven conditions. The results are shown in Figure 2.

Figure 2

Download asset Open asset

In each case, we set all kinetic parameters to values specified in the original publications and coupled the external driving signal in the way specified in those original publications. As in the Kai clock simulations, the external signal was a square wave with amplitude fluctuations of fixed strength. Finally, we add Langevin noise to the equations to simulate internal noise; when available, we followed the finite volume prescription for rates in the original publications or related papers to set the size of Langevin noise for each reaction. Simulation and model details are in Appendix 2.

These models here all exhibit a Hopf bifurcation as kinetic parameters are tuned. The publications (Leloup et al., 1999; Schmal et al., 2014; Locke et al., 2005; Leloup and Goldbeter, 2003; Goldbeter, 1991; Goodwin, 1965; Gonze and Abou-Jaoudé, 2013; Kondepudi and Prigogine, 2014; Buşe et al., 2009) identified a parameter which when tuned leads to a Hopf bifurcation; that is, on one side of the bifurcation, we find damped oscillations while on the other side, we find free running oscillations of increasing amplitude. We picked three points along this parameter; the green and purple data correspond to free running oscillations of large and smaller natural amplitude relative to the size of the external drive. The red data corresponds to a choice of parameters on the other side of the Hopf bifurcation, that is, to damped oscillations. For the red data, we chose $\mu$ such that the relaxation timescale was comparable to the period of the external driving force, much as in the Kai model of P. marinus. The damped oscillator is a useful clock only when the relaxation timescale is comparable to the period.

In each case, we observed the same trade-off as seen in the Kai system; free running oscillations of large amplitude relative to the external drive (green) were most precise when only subject to external noise but are most fragile to internal noise. Damped oscillations in the same oscillator models are more robust and thus are preferable at sufficiently high internal noise. We find that intermediate amplitude free running oscillations show intermediate noise properties. Consequently, we can continuously trade-off resistance to internal noise for resistance to external noise by changing the amplitude of free running oscillations relative to the strength of the external drive.

We have demonstrated a trade-off between external and internal noise resistance in diverse clocks. While it is possible to trace the origin of this trade-off to specific features of each clock, here, we take a different approach based on dynamical systems theory. Dynamical systems theory has been use to make fruitful general predictions about biological clocks since Winfree’s analysis of phase singularities (Winfree, 2001). In a similar vein, we use dynamical systems theory to show this trade-off emerges due to basic features of free running and damped clock dynamics and can thus be expected to hold broadly.

Free running clocks are phenomenologically well-described by a limit cycle attractor, a non-linear oscillator of fixed amplitude (Winfree, 2001). While such descriptions have been used for numerous biochemical oscillators, limit cycle dynamics can be experimentally seen in molecular detail for the KaiABC clock in S. elongatus as shown in Figure 3a (reproduced from [Leypunskiy et al., 2017]). The clock follows distinct limit cycle dynamics during the day (orange data) and night (black data) (Leypunskiy et al., 2017; Pattanayak et al., 2014), with the day cycle positioned at higher phosphorylation levels due to a higher ATP/ADP ratio.

Figure 3

Download asset Open asset

The Kai model and indeed the diverse range of biochemical oscillators in Figure 2 show such a change in the limit cycle between day and night conditions. Here, we build a minimal model of such free-running clocks using circular day and night limit cycles of radius $R$ in a plane. The limit cycle is defined by the dynamics ${\tau }_{relax}\dot{r}=r-r^3/R^2,\dot{\theta }=\omega$ about its own center; but the center of the limit cycle itself moves along the x axis in Figure 3b as $(\rho (t)L,0)$ where $\rho (t)\in [0,1]$ is the normalized light intensity level at time $t$ and $L$ is a measure of the physiological changes between day and night (e.g., ATP/ADP ratio change in S. elongatus). Thus, for example in Figure 3b, the system follows the blue dynamics at night and then after dawn it relaxes to the orange day attractor on a time scale ${\tau }_{relax}$. Note that $R$ represents the amplitude of free-running oscillations while $L$ represents the strength or amplitude of the external driving signal.

In contrast, damped clocks are phenomenologically well-described by a day-time and a night-time point attractor with slow relaxation dynamics between them (Figure 3c), modeled as $\dot{r}=-r/{\tau }_{relax},\dot{\theta }=\omega$ about an attractor point whose location varies with current light levels as $(-\rho (t)L,0)$. Here we assume $2{\tau }_{relax}\sim 24$ hrs as in P. marinus (Holtzendorff et al., 2008); if relaxation were faster and completed before the day is over, the clock cannot resolve all times of the day.

Here, we will also consider a family of limit cycle clocks of varying $R/L$ to interpolate between large-$R/L$ limit cycles and point attractors (approximated by $R/L=0$).

We begin with the performance of different clocks in the presence of external intensity fluctuations. Weather patterns cause large fluctuations in the intensity of light over a wide range of time-scales as shown in Figure 4a. Much like with biochemical circuits, we subject an in silico population of dynamical system clock models to different realizations of such noisy weather patterns.

Figure 4

Download asset Open asset

When subject to weather fluctuations, we see in Figure 4b that the population variance of clock states for limit cycles at given times (purple) is fundamentally limited by the spacing between the day and night limit cycles. Point attractors develop larger overlapping population distributions at different times.

We can geometrically understand the daytime phase variance increase ${\sigma }_{clouds}^2$ in terms of the phase lag $\mathrm{\Delta }\mathrm{\Phi }$ due to a single, say $2.4$ hr dark pulse, administered during the day. Figure 4c shows that the deviation in trajectory for limit cycle clocks (purple) is fundamentally limited by the presence of a continuous attractor. In contrast, for the point attractor, the trajectory is in free fall towards the night point attractor, with no limit cycle to arrest such a fall. Consequently, the geometrically computed phase shift $\mathrm{\Delta }\mathrm{\Phi }$ due to the particular dark pulse shown in Figure 4c is much smaller for limit cycles ($\mathrm{\Delta }\mathrm{\Phi }\sim 0.5$ hr for the $R,L$ geometry shown) than for point attractors ($\mathrm{\Delta }\mathrm{\Phi }\sim 4$ hr) (see Appendix 5).

In fact, this contrast in $\mathrm{\Delta }\mathrm{\Phi }$ between limit cycles and point attractors holds for dark pulses of any duration and time of occurrence. The contrast is even greater at small $L/R$ since ${(\mathrm{\Delta }\mathrm{\Phi })}^2\sim {(L/R)}^2$ for small $L/R$, as shown geometrically in Appendix 5 and confirmed in simulations that average over random weather conditions (Figure 4d). Hence, limit cycles are more resistant to external fluctuations than point attractors.

To complete the analysis, we note that phase variance increases additively during the day and falls multiplicatively at dusk (and dawn), that is,

Solving for steady state phase variance (${\sigma }^2=({\sigma }^2+{\sigma }_{clouds}^2)/s^4$), we obtain

where we have equated ${\sigma }_{clouds}^2$ to $\mathrm{\Delta }{\mathrm{\Phi }}^2$ for a typical dark pulse Here, $s^2$ represents the variance drop during a dawn/dusk entrainment. As shown in Appendix 5 for external noise (and in Figure 5 for internal noise), this factor $s$, can be geometrically explained by the slope of the circle map relating the two cycles Leypunskiy et al., 2017; we find that $s^2-1\sim L/R$ for large-$R/L$ limit cycles. Plugging this and $\mathrm{\Delta }{\mathrm{\Phi }}^2\sim {(L/R)}^2$ into Equation1, we see that ${\sigma }^2\to L/R\to 0$ for large -$R/L$ cycles.

Figure 4e shows that the precision (i.e., mutual information between clock state and time) computed from random weather simulations agrees with this theory; clock precision drops as we interpolate from limit cycles to point attractors by changing $L$ (with equivalent results for changing $R$).

Figure 5

Download asset Open asset

Internal noise due to finite copy number effects in biochemical networks can be modeled exactly using the Gillespie method used in Figure 1. In the context of our dynamical systems model, we follow Gillespie, 2007 and add Langevin noise to all dynamical variables of the system of strength ${ϵ}_{int}\sim 1/\sqrt{N}$, where $N$ is the overall copy number, with the ratios of different species assumed fixed (see Appendix 3). Such a Langevin approach is an approximation Gillespie, 2007 to the exact Gillespie method used in Figure 1.

We simulated a population of clocks in externally noiseless day-night light cycles but with internal Langevin noise. We see in Figure 5b that limit cycle populations have significantly higher variance of clock state due to internal noise than point attractors, in contrast to Figure 4b with external noise alone.

We can understand the weakness of limit cycle attractor relative to the point attractor in terms of diffusion during day/night balanced by dawn/dusk transitions. The flat direction along the limit cycle attractor cannot contain diffusion caused by the Langevin noise during the day/night; hence the phase variance increases by ${\sigma }^2\to {\sigma }^2+{ϵ}_{int}^2{T}_{day}$ during a day of length $T_{day}$ (and similarly at night).

Dawn and dusk times reduce the phase variance ${\sigma }^2\to {\sigma }^2/s^2$ as the trajectories originating on, say, the day cycle converge on the night cycle (see Figure 5c and Leypunskiy et al., 2017; Monti and Lubensky, 2017). In fact, we can compute this variance drop $s^2$ entirely through geometric considerations. We define the circle map $\varphi =P(\theta )$ as relating originating points $\theta$ near dusk on the day cycle to final points on the night cycle $\varphi$ after relaxation (experimentally characterized in Leypunskiy et al., 2017). Then $s^{-1}=dP(\theta )/d\theta$. Figure 5c shows that this slope $s^{-1}=dP(\theta )/d\theta$, geometrically computed in the SI, agrees with the dawn/dusk variance drop in Langevin simulations and scales as $s^2-1\sim L/R$ for large $R/L$.

Thus, the population phase variance changes as

Assuming $T=T_{day}=T_{night}$ and solving for steady-state phase variance (${\sigma }^2=(({\sigma }^2+{ϵ}_{int}^2{T}_{day})/s^2+{ϵ}_{int}^2{T}_{day})/s^2$), we obtain

Consequently, as the cycles become large (large $R/L$), the dawn/dusk variance drop vanishes as $s^2-1\sim L/R\to 0$ while diffusion along the flat direction still adds ${ϵ}_{int}^{2}T$ to the variance during each day and each night; hence large-$R/L$ limit cycles have large ${\sigma }_{cycle}^{2,int}$ and thus low precision. (Unlike with external noise, internal noise introduces a diffusion length scale and hence changing $L$ and $R$ are not equivalent. To make a fair comparison, we fix $R$ and internal noise while changing $L$ in Figure 5e; see Appendix 3 for more detail about other equivalent choices).

Note that Equation 2 is invalid for strictly undriven clocks (i.e., $s=1$); such clocks show a variance that increases indefinitely with time. Here, we focus on driven clocks, which always settle to the finite variance given by Equation 2.

In contrast, for point attractors, the population variance stays constant during the day-night cycle and is shown to be

in the SI, which matches Langevin simulations (Figure 5d). Since ${\tau }_{relax}\sim T_{day}$ to have distinct clock states through the day (Figure 3), we find ${\sigma }_{cycle}^{2,int}\ge {\sigma }_{point}^{2,int}$.

In summary, in both cases, population variance is reduced by the geometric ‘curvature’ of the dynamics, that is, convergence of nearby trajectories. Point attractor trajectories experience a constant curvature of $1/{\tau }_{relax}$. But limit cycle clocks experience such ‘curved’ off-attractor dynamics only at dawn and dusk, which is offset by dephasing (Mihalcescu et al., 2004; Gonze et al., 2002) during long periods of zero curvature on the limit cycle (day/night). Hence limit cycles underperform point attractors under high internal noise.

We now subject the clock systems to both internal and external noise at the same time. We find results (see Figure 6a) that parallel those for explicit molecular models of biochemical oscillators studied in Figure 2. Large-$R/L$ limit cycles outperform other clocks in filtering out external noise when internal noise is low, but their precision degrades more rapidly than other clocks as internal noise ${ϵ}_{int}^2\sim 1/N$ is increased. Point attractors have poor precision with only external noise but do not significantly degrade with internal noise and outperform all other clocks at high internal noise. At comparable strengths of internal and external noise, limit cycles with an intermediate value of $R/L$ are most precise. In the SI, we show that the optimal geometry is set by the ratio of internal and external noise strength,

Figure 6

Download asset Open asset

In the SI, also we show that, under certain simplifying assumptions, Equations 1 and 2 can be combined to give an explicit trade-off relationship,

where $Q={ϵ}_{int}^2{ϵ}_{ext}^2$ and where ${\sigma }_{int}^2$ is the population angular variance of the clock state due to internal noise when driven by a noiseless external signal and ${\sigma }_{ext}^2$ is the population angular variance in the absence of internal noise due to amplitude fluctuations of the external signal. Note that angular variance is a better indicator than variance because we want to know how well the system can tell time.

Equation 4 makes our trade-off explicit and also clarifies which parameters are varied and which parameters are held fixed in this trade-off. As long as $Q$ is held fixed, we allow all other parameters to vary – for example, the overall strength of the external drive $L$, the size of the cycle $R$, and as discussed in the SI, all other parameters characterizing the normal form of limit cycles near a Hopf bifurcation.

However, in holding $Q$ fixed, our trade-off does assume that the strength of the external fluctuations ${ϵ}_{ext}$ – that is, the fractional size of amplitude fluctuations in the external signal – is held fixed. Similarly, we hold ${ϵ}_{int}^2$, the phase diffusion constant, fixed – that is, we are comparing clocks that would show the same population phase variance (in units of radians) over the same time in undriven conditions. See Appendix 3 for alternative comparisons and other details.

Another measure of clock quality is the entrainment speed, that is, the time taken to reach steady state population variance, starting from a population uniformly distributed in clock phase. In Figure 6b, we see that with external noise only, the most precise clocks (i.e., small-$L/R$ limit cycles) are the slowest to entrain because they retain a longer history of the external signal, allowing them to average out external noise better.

But strikingly, such a speed-precision trade-off is absent if internal noise is high. In Figure 6c, only internal noise is present and the external signal has no fluctuations. We see that clocks most robust to internal noise are also the fastest to entrain. Intuitively, the phase of slow entraining clocks is affected by the cumulative effect of internal fluctuations over a longer period of time. With both external and internal noise present, clocks with intermediate entraining speed - that is, intermediate ${(L/R)}_{optimal}={ϵ}_{int}/{ϵ}_{ext}$ - will have the highest precision.

We demonstrate our trade-off using Gillespie simulations of an explicit biomolecular KaiABC model of the post-translational clocks in S. elongatus and P. marinus.

The S. elongatus clock has been well-characterized experimentally (Bryant, 2003; Gutu et al., 2013; Dufresne et al., 2003; Kitayama et al., 2003 - see Appendix 1—figure 1a). The clock is fundamentally based on the ordered phosphorylation and dephosphorylation of KaiC (Rust et al., 2007). Phosphorylation of KaiC is KaiA-dependent which allows for feedback that enables collective coherent oscillations in a cell. After complete phosphorylation of KaiA-C complexes (usually by the end of the day), KaiC forms a KaiB-C complex which then dephosphorylates in an ordered manner. Crucially, the KaiB-C complex also sequesters KaiA in a KaiABC complex, reducing the pool of available KaiA for phosphorylation of other KaiC hexamers. This negative feedback enables coherent oscillations of the population of KaiC molecules in a single cell (Rust et al., 2007).

Appendix 1—figure 1

Download asset Open asset

P. marinus lacks the kaiA gene but possesses and expresses kaiB and kaiC. While the details of the protein clock dynamics are not as fully known as with S. elongatus, gene expression shows cycling in cycling conditions but decays in constant conditions (Holtzendorff et al., 2008). A conservative model, consistent with all these known facts about P. marinus, is shown in Appendix 1—figure 1b; without KaiA feedback, different hexamer units phosphorylate independently and settle to a hyperphosphorylated state at the end of the day. At night, they dephosphorylate along a distinct pathway (homologous to that used by S. elongatus but without KaiA) and reach a hypophosphorylated state by dawn.

We created the following hybrid model that includes S. elongatus and P. marinus models as different limits. In our model, shown in Appendix 1—figure 1c, KaiC has a KaiA-dependent phosphorylation pathway, much like in S. elongatus, that is used during the day and driven forward by ATP.

But to also include P. marinus-like behavior in the model, we allow for a second parallel phosphorylation pathway for KaiC that is independent of KaiA. The relative access of these two pathways is controlled by a parameter $\gamma$. When $\gamma =1$, only the S. elongatus-like KaiA dependent pathway is accessible. When $\gamma =0$, only the P. marinus-like KaiA independent pathway is accessible. Collectively, we call these states along these phosphorylation pathways, the UP states of KaiC - phosphorylation are going UP along these pathways which are usually used during the day.

After maximum phosphorylation (usually at dusk), KaiA unbinds (if present) and a KaiB-based dephosphorylation pathway takes over (common to both systems). We call these states the DOWN states of KaiC.

Critically, KaiA is assumed to be sequestered through the formation of KaiABC complexes during this dephosphorylation stage. In S. elongatus, reduced KaiA availability prevents other KaiC hexamers from proceeding independently through the UP stage while most of the population is in the DOWN state. Such negative feedback is critical in maintaining free-running limit cycle oscillations in S. elongatus.

However, as $\gamma \to 0$, the KaiA-independent pathway is more active and thus the system effectively has no feedback. In fact, we find that at about $\gamma \approx 0.82$, sustained oscillations disappear (for kinetic parameters used here and reported below). Hence we chose $\gamma =1,0.95,0$ as representative of two limit cycle-based free running clocks and one point-attractor based damped clock respectively. In this way, we can view the clock dynamics of S. elongatus and P. marinus can be viewed as being on either side of the Hopf bifurcation that occurs at $\gamma \approx 0.82$.

We ran explicit Gillespie simulations corresponding to the deterministic equations above at different overall copy number $N$ with fixed stoichiometric ratios of the molecules KaiA,B,C.

We simulated external input noise by varying the ATP levels during the day. External noise in these simulations were implemented by changing ATP levels in the following way: we fluctuated the ATP levels $f_{ATP}=ATP/(ATP+ADP)$ during the day between the $f_{ATP}^{day}$ and $f_{ATP}^{night}+(f_{ATP}^{day}-f_{ATP}^{night})/3$, where $f_{ATP}^{day}$, $f_{ATP}^{night}$ are the ATP values during a cloudless day and night respectively. We used different day and night ATP levels for different $\gamma$ that ensure that the limit cycles had periods comparable to $24$ hours. For $\gamma =1$, we used ATP/ADP ratios of $f_{ATP}^{day}=0.55,f_{ATP}^{night}=0.45$. For $\gamma =0.95$, we used $f_{ATP}^{day}=0.57,f_{ATP}^{night}=0.17$ and for $\gamma =0$, $f_{ATP}^{day}=0.8,f_{ATP}^{night}=0.2$. The corresponding limit cycles and point attractors are shown in Appendix 1—figure 1d.

We used the following kinetic parameters in all simulations: ${\displaystyle dt=0.01hr,}$ ${\displaystyle k_{+}=k_{-}=2m\cdot 0.04932h{r}^{-1},}$ ${\displaystyle k_{Aon}=0.2466\mu M^{-1}h{r}^{-1},}$ ${\displaystyle k_{Aoff}=0.02466h{r}^{-1},}$ ${\displaystyle k_{C\to C\ast }=0.2466h{r}^{-1},}$ ${\displaystyle k_{ABC}=123.30h{r}^{-1},m=18}$. We set up Kai C and Kai A in a $1:1$ stoichiometric ratio, each present at a copy number $N$ where $N$ was varied systematically. These rates are consistent with those measured in (Qin et al., 2010; Snijder et al., 2017; Hayashi et al., 2004).

Much like with Langevin simulations of dynamical systems performed in this paper, we run the Gillespie simulation until equilibration of the population. However, the system appears to reach the equilibrium state much faster (only over five light-dark cycles of 12 hr: 12 hr). We extracted one day of such a trajectory on day six and repeated the simulation 100–400 times. We repeat 400 times when the copy number is low ($<1200$) since the spread will be big and we found that the probability distribution is not smooth. We run only 100 times for the high copy number ($>1200$). Pooling together these trajectories, we computed the mutual information between clock state (i.e., $(u,d)$ where $u$ is the net phosphorylation state of KaiC in the up-pathways and $d$ is the net phosphorylation state of KaiC in the KaiB-bound ‘down’ pathways in Appendix 1—figure 1c ) and time of day. The $(u,d)$ space was binned using bins of fixed size of dimension $(0.05,0.05)$ while the 24 hr time-of-day was binned with bins of size $0.5$ hrs.

With these choices, we see in Appendix 1—figure 1d that this model has limit cycles of different position during the day and night. The corresponding experimental data, reproduced from Leypunskiy et al. (2017), are presented in the main paper.

Here, we study the effect of internal and external noise on a diverse array of biochemical oscillator models from the literature in the parameter regimes described in the original publications. We confirm the same trade-off described in the paper in these models; a summary of our results is presented in the main paper.

In all of the models described below, we set all parameters to values used in the original or cited papers with only two exceptions: (a) the parameter identified as coupling to external signals in these publications is varied over time as a square wave with amplitude fluctuations added, (b) the parameter designated by the relevant original publication as controlling the distance from the Hopf bifurcation was used to simulate a point attractor-based ‘hourglass’ oscillator (red lines in Figure 2 of the main paper) and limit cycles of different free running oscillation amplitude (green and purple lines in Figure 2 of the main paper). This latter parameter roughly corresponds to $R$, the size of limit cycle, while the amplitude of square wave coupled to the former parameter corresponds to $L$, the separation of the limit cycles, in our dynamical systems theory, that is, the separation of the ‘day’ and ‘night’ limit cycles. (In several papers, these two are the same parameter, in which case the day-night difference reflects $L$ while the mean value reflects $R$.) Finally, we add Langevin noise to the equations to simulate internal noise; when available, we followed the finite volume prescription for rates in these papers to set the size of Langevin noise for each reaction.

We keep the strength of external noise ${ϵ}_{ext}$, defined as the noise-to-signal ratio of the amplitude fluctuations in the external signal, fixed. We varied internal noise ${ϵ}_{int}$ along the $x$ axis of plots in Figure 2 of the main paper. Here, ${ϵ}_{int}$ is defined as the phase diffusion constant of a clock in undriven conditions (see how we define internal noise in the section on Neurospora and Drosophilia below); this normalization, which depends on the Hopf bifurcation parameter in (b) above, allows us to make a fair comparison between different clocks since they develop the same phase variance over the same time in undriven conditions.

As seen in Figure 2 of the main paper, these diverse models agree with the trends found in our analysis of dynamical systems and with simulations of the KaiABC system, showing that our results are not tied to any particular molecular model.

The circadian clock in Neurospora has been modeled (Leloup et al., 1999) as arising from interactions between mRNA ($M$) and a protein that can shuttle in and out of a nucleus ($P_N,P_C$). The equations used in Gonze and Goldbeter (2006) to model this are,

where ${\nu }_s$ is an mRNA transcription rate that is modulated by external signals (Leloup et al., 1999; Gonze and Goldbeter, 2006), and $\mathrm{\Omega }$ is the volume of the system which in turn determines the strength of stochastic noise. (A model with very similar equations has also been suggested as a model of the Drosophila circadian clock [Leloup et al., 1999].)

We use the same parameters used in the Ref. (Leloup et al., 1999; Gonze and Goldbeter, 2006): $K_I$ = 1 nM, $n$ = 4, $v_m$ = 0.505 nM/h, $K_m$ = 0.5 nM, $k_s$= 0.5 1/h, $v_d$ = 1.4 nM/h, $K_d$= 0.13 nM, $k_1$ = 0.5 nM/h, $k_2$ = 0.6 nM/h, and assume the volume $\mathrm{\Omega }$ dependence of these parameters to be exactly as used in Gonze and Goldbeter (2006). We add internal stochasticity by adding Langevin noise with a diffusion matrix (Gonze and Goldbeter, 2006):

where $\mu (x,t)$ is the RHS of Equation 5, $\eta (0,1)$ is a vector whose entries are independent standardized Gaussian noise (mean 0, variance 1), and

where $A=v_s\mathrm{\Omega }\frac{{(K_I\mathrm{\Omega })}^n}{{(K_I\mathrm{\Omega })}^n+P_N^n}$ and $B=v_m\mathrm{\Omega }\frac{M}{K_m\mathrm{\Omega }+M}$. This is how internal noise get added into other oscillators models as well. However, for the system of equations that use concentration instead of the number of molecules, the equation has to be modified to $\frac{dX}{dt}=\mu (x,t)+\frac{1}{\sqrt{\mathrm{\Omega }}}\mathrm{\Sigma }(x,t)\eta (0,1)$.

As in Leloup et al. (1999) and Gonze and Goldbeter (2006), we take ${\nu }_s$ to be modulated by the external signal (light). As shown in Leloup et al. (1999) and Gonze and Goldbeter (2006), a Hopf bifurcation occurs at ${\nu }_s$ = 0.57 nM/h. Hence, in Figure 2 from the main text, we used ${\nu }_s^{Day}$ = 0.55 nM/h,${\nu }_s^{Night}$=0.05 nM/h for the point attractor (red). For the two limit cycles, we used ${\nu }_s^{Day}$ = 0.9 nM/h,${\nu }_s^{Night}$=0.6 nM/h (green), and ${\nu }_s^{Day}$ = 0.705 nM/h,${\nu }_s^{Night}$=0.695 nM/h (purple). The driving period is 18 hr, similar to the driving period of the system at $v_1$ = 0.7 nM/h

A model of the circadian clock in Arabidopsis thaliana was introduced in Locke et al. (2005). While many biologically important features have been added in the years since then, the original model was based on a single negative feedback loop and involves two transcription factors (LHY and CCA1) that inhibit their activator $TOC1$. A reduced model with the same phenomenology was presented in Schmal et al. (2014), in which LHY and CCA1 are combined into one variable, representing their mRNA and protein levels by $M_L(t)$ and $P_L(t)$ respectively. Denoting the mRNA and protein levels of $TOC1$ by $M_T(t)$ and $P_T(t)$, Schmal et al. (2014) present a reduced version of the model in Locke et al. (2005) as:

where $L(t)$ is a light input function, and other parameters except the variables specified on the left hand sides are constant.

To simulate this system, we use the parameters used in Schmal et al. (2014): a = 2, b = 2, $g_1$ = 0.5 nM, $g_2$ = 0.1 nM, $m_1$ = 0.4 nM/h, $m_2$ = 0.6 nM/h, $m_3$ = 0.6 nM/h, $m_4$ = 0.3 nM/h, $k_1$ = 1 nM, $k_2$ = 0.5 nM, $k_3$ = 1 nM, $k_3$ = 1 nM, $p_1$ = 0.5 1/h, $p_2$ = 0.3 1/h, $v_2$ = 0.6 nM/h. With other parameters fixed, the system undergoes Hopf bifurcation at $v_1$ = 0.194 nM/h We use $v_1$ = 0.26 nM/h for limit cycles and $v_1$ = 0.05 nM/h for point attractor. $L(t)$ is a light input function. For the two limit cycles in Figure 2 in the main text, we set $L^{Day}$ = 0.05 nM/h and $L^{Night}$ = 0 (green data) and we set $L^{Day}$ = 0.01 nM/h and $L^{Night}$ = 0 (purple data). For point attractor, we set $L^{Day}$ = 0.2 nM/h and $L^{Night}$ = 0 (red data). The period of the driving signal is 24 hr, which is around the natural period of the system when $v_1$ = 0.26 nM/h and $L$ = 0.

The circadian clock in mammalian cells was modeled in detail by Leloup and Goldbeter (2003), using $19$ equations representing the interactions between $Per,Cry$ and other genes. We simulate this entire system explicitly with the parameter values specified in the original publication (Leloup and Goldbeter, 2003). To introduce Langevin noise, we use a simplified diagonal diffusion matrix with entry $\sqrt{D{X}_i}$ for species $X_i$. We do not reproduce these $19$ equations or parameter values used from Leloup and Goldbeter 2003) here in interest of space; the only modification we made is to introduce Langevin noise to each of the $19$ equations.

Leloup and Goldbeter (2003) identified parameter $v_{sP}$ (a transcriptional rate) to be the light input function. We use $v_{sP}^{Day}$ = 1.09 nM/h and $v_{sP}^{Night}$ = 1.07 nM/h for the purple limit cycle data in Figure 2, $v_{sP}^{Day}$ = 1.15 nM/h and $v_{sP}^{Night}$ = 1.07 nM/h for the green limit cycle data. For the point attactor data (red), we set $v_{sP}^{Day}$ = 1.5 nM/h and $v_{sP}^{Night}$ = 0. In addition, Leloup and Goldbeter (2003) identified parameters $KAC,vmB$ as controlling the distance from the Hopf bifurcation. For the point attractor, we used KAC = 0.4 nM, and vmB = 0.9 nM/h (also used in Leloup and Goldbeter, 2003). For the limit cycles, we used KAC = 0.6 nM, and vmB = 0.8 nM/h, which lies on the other side of the Hopf bifurcation. The period of the input signal is at 21.5 hr, which is around the natural period of the system when $v_{sP}$ = 1.07 nM/h.

A classic model of the cell cycle was proposed by Goldbeter (1991). While many additional details have been added on since then, the model captures the essential mechanism behind the self-sustained nature of cell cycles.

where $\mathrm{\Omega }$ is the size of the system and other parameters are constants. The three variables are the cyclin concentration $C$, the fraction of active cdc2 kinase $M$, and the fraction of active cyclin protease $X$. For $C$, the internal noise is proportional to the square root of the rates, but for $M$ and $X$, it is proportional to the square root of the rates divided by $\mathrm{\Omega }$ because they are fractions and not concentrations (following the prescription in Gonze and Goldbeter, 2006 for a similar model). Parameter values were taken from Goldbeter (1991): $K_i=0.1$ ($i$=1–4), $V_{M1}=0.5{\text{min}}^{-1}$, $V_2=0.167{\text{min}}^{-}1$, $V_{M3}=0.2{\text{min}}^{-1}$, $V_4=0.1{\text{min}}^{-1}$, $v_d=0.1\mu M{\text{min}}^{-1}$, $K_C=0.3\mu M$, $K_d=0.02\mu M$, $k_d=3.33\times {10}^{-3}{\text{min}}^{-1}$.

Goldbeter (Goldbeter, 1991) suggested that $v_i$ is modulated by external signals. So, we use $v_i^{Day}$ = 0.0106 $\mu M{\text{min}}^{-1}$ and $v_i^{Night}$ = 0.0105 $\mu M{\text{min}}^{-1}$ for small L/R limit cycle, $v_i^{Day}$ = 0.0111 $\mu M{\text{min}}^{-1}$ and $v_i^{Night}$ = 0.0105 $\mu M{\text{min}}^{-1}$ for large L/R limit cycle, and $v_i^{Day}$ = 0.009 $\mu M{\text{min}}^{-1}$ and $v_i^{Night}$ = 0 for point attractor. The bifurcation from point attractor to limit cycle happen around $v_i$ = 0.01 $\mu M{\text{min}}^{-1}$. The period of the driving signal is 35 min.

One of the earliest models of biochemical oscillators was proposed by Goodwin (1965) (later corrected). We use the simplest widely-studied version of such a Goodwin oscillator (Gonze and Abou-Jaoudé, 2013; Woller et al., 2014),

When $n=9$, the limit cycles disappear at a Hopf bifurcation found at $\alpha \approx 7$. As is commonly done (Woller et al., 2014), we couple the external signal to the bifurcation parameter $\alpha (t)$. We use ${\alpha }^{Day}=120,{\alpha }^{Night}=80$ for the green limit cycle in Figure 2c of the main paper, ${\alpha }^{Day}=108,{\alpha }^{Night}=92$ for purple limit cycle data, and ${\alpha }^{Day}=2.5,{\alpha }^{Night}=1$ for the red point attractor data. The input signal has a period of 4, which is roughly the natural period of the limit cycle at ${\displaystyle \alpha =100}$ are taken to the output of the clock for computing MI.

The repressilator is a model of an early synthetic biology system (Elowitz and Leibler, 2000) that demonstrated oscillations in a synthetically wired gene regulatory circuit. While resembling the Goodwin oscillator in topology, the network has the total non-linearity spread equally amongst all three reactions, lowering the requisite Hill coefficient of any one reaction to a biochemically realistic $n=3$. Repressilator circuits are not usually driven by an external signal, except in a few theoretical analyses (e.g., Russo et al., 2010; Schultz et al., 2013). We use the simplest version of these, with the driving signal modulating the transcription rate of only one of the proteins

where $\alpha$ is the bifurcation parameter and $s(t)$ is the variable coupled to the input signal. In a non-driven repressilator, $s(t)=0$. When $n=3$, the Hopf bifurcation occurs at $\alpha \approx 2.5$, so for limit cycles, we use $\alpha =5.2$ and for point attractor we use $\alpha =1.9$.

We use $s^{Day}=0.7/5.2,s^{Night}=-0.7/5.2$ for the green limit cycle in Figure 2g in the main paper, $s^{Day}=4.8/5.2,s^{Night}=-1.7/5.2$ for the purple limit cycle data, and $s^{Day}=0.5/1.9,s^{Night}=-1.9/1.9$ for the point attractor (red). The input signal has a period of 4, which is roughly the natural period of the limit cycle at $\alpha =5.2$ and ${\displaystyle s(t)=0}$ are taken to the output of the clock for computing $MI$.

The Brusselator is a model of autocatalytic reactions that show limit cycle oscillations. This model has been extensively studied over the years; while the explicit biochemical reactions can be found in Kondepudi and Prigogine, 2014), these reactions are modeled by the ODEs:

where $b$ has been identified as a bifurcation parameter (Kondepudi and Prigogine, 2014). Most studies do not consider driven Brusselator models; we follow the driving prescriptions of the Goodwin model and couple the external light to the bifurcation parameter, converting the constant $b$ into $b(t)$. The bifurcation point is at $b=2$. For $b<2$, we have a point attractor and for $b>2$ we have a limit cycle. We use $b^{Day}=2.25$ and $b^{Night}=2.2$ for the purple limit cycle data in Figure 2f from main text. We use $b^{Day}=2.8$ and $b^{Night}=2.2$ for the green limit cycle data. Lastly, we use $b^{Day}=1.8$ and $b^{Night}=0.5$ for the point attractor (red). The signal has a period of 6.4 which is around the natural period of the system when $b=2.2$.

To complement our study of detailed biochemical implementations of such systems, we study two kinds of dynamical systems in this paper; limit cycles and point attractors. The minimal model of limit cycles and point attractors is given by the ‘normal form’ near a Hopf bifurcation:

For $\alpha >0$, the above equation describes a circular limit cycle of radius $R$ and frequency $\omega$. This equation undergoes a Hopf bifurcation at $\alpha =0$, where the limit cycle shrinks to zero and resulting in point attractor for $\alpha <0$. The ‘normal form’ can be seen as the universal simple form – for example, circular limit cycles of radius $R$ – that any limit cycle and point attractor will reduce to in the neighborhood of a Hopf bifurcation. We add White noise in the Cartesian space representation of the Dynamical equations to represent the internal noise as follows:

where $D\sim R^2{ϵ}_{int}^2$ is the diffusion constant, d$W$ is a Wiener process, and here we assume that the internal noise is a homogeneous white-noise in the 2-dimensional space. (Similar assumptions are made in reference Potoyan and Wolynes, 2014)

While we assume this simple form here as a minimal model, we do not assume that the oscillator is weakly driven. Instead, based on experimental observations of the Kai clock (Leypunskiy et al., 2017) and models of numerous other clocks (Winfree, 2001), we assume that the origin of the limit cycle or point attractor equations above moves by a finite amount $L$ as the external light signal switches between day and night values. In fact, we move the origin along the x-axis as a function of time as $(Ls(t),0)$ where $s(t)$ is the external light signal, assumed to be of amplitude $1$. Thus we are assuming a simple circular form of limit cycles and point attractors but do not restrict to weak driving. (In the limit of weak driving, that is, small $L/R$, our model can be shown to reduce to the universal Stuart-Landau model of weakly driven oscillators as a special case.)

In Equation 14, ${\tau }_{\text{relax}}\sim \frac{1}{|\alpha |}$ is the relaxation time for perturbations away from the limit cycle or point attractor. For limit cycles, perturbations away from the limit cycle tend to decay fast relative to the period $2\pi /\omega$, typically on the order of hours (Leypunskiy et al., 2017).

In contrast, the point attractor in damped ‘hourglass’ clocks P. marinus needs to have relaxation ${\tau }_{\text{relax}}\sim \frac{1}{|\alpha |}\sim 2\pi /\omega$ comparable to the period of the day-night cycle itself. As explained in the main paper, if relaxation were much faster, the clock state would decay to a fixed point before the end of the day (or night) and thus not show distinct states at distinct times of the day-night cycle.

For both limit cycles and point attractors, we simulate a population of clocks, each represented by a particle in the given dynamical system, subject to external and/or internal noise.

We use $\alpha =5$ for limit cycle system and $\alpha =-5$ for point attractor system where $\omega =2\pi$ in these units. For point attractors, we set $R=1000L$, where $L$ is the separation of the day and night attractor. In such a limit, the point attractors are quadratic potentials with linear restoring forces since $\frac{r^3}{R^2}$ is small. The center of the cycle and point attractors during the day are assumed to be at $(-L,0)$ and at $(0,0)$ at night; or more generally at $(-Ls(t),0)$ were $s(t)$ the light signal (assumed to be of amplitude $1$).

To simulate external noise, we use a square wave signal $s(t)$ of amplitude $1$ with amplitude fluctuations set by ${ϵ}_{ext}$. As explained in Appendix 4 on external noise, we take the fluctuations in $s(t)$ to have a correlation time of $2.4$ hours.

To simulate internal noise, we add Langevin noise to Equations 14 as described in Appendix 6 on Langevin noise. As explained in Appendix 6, our measure of internal noise ${ϵ}_{int}^2$ is a measure of phase diffusion, independent of limit cycle size. In other words, ${ϵ}_{int}^2$ is a measure of the population phase variance (i.e., variance in $\theta$) developed by limit cycles of any size in undriven conditions in a given period of time.

To interpolate between limit cycles and point attractors, we systematically change $L$ holding $R$ fixed. For limit cycle simulations, changing $L$ and $R$ are equivalent. To see this for external noise simulations, note that $L/R$ is the only dimensionless parameter. For internal noise, our definition of ${ϵ}_{int}^2$ above as the phase diffusion in undriven conditions for limit cycles of any size, ensures that changing $L$ and $R$ are equivalent for limit cycles.

For point attractors, we set the separation $L$ be the diameter of the limit cycles simulated in the same plots to keep the size of the resolvable chemical spaces roughly comparable for limit cycles and point attractors. While this precise choice is arbitrary to some extent, note that the point attractor results for external noise do not depend on this parameter at all since $L$ is the only relevant length scale in external noise simulations. The separation $L$ does affect the clock precision with internal noise (of fixed absolute strength ${ϵ}_{int}$) but the mutual information changes only logarithmically with $L$.

We evolve our dynamical system using the Euler method with time step $dt=0.001$ days until the value of mutual information from one day to the next does not change by more than 2–3% - that is. the system has reached steady state. Reaching steady-state usually takes around 200 days, but if the ratio of $L/R$ is smaller than $0.1$, then we may need to run the simulation until day 500 to reach an equilibrium (See speed-error tradeoff in Figure 6b and c in the main text).

For limit cycles, we initialize the population of ${10}^4$ particles by uniformly distributing them along the perimeter of the night cycle. In the point attractor system, we initialize a population of ${10}^5$ at the night-time point attractor.

We use a larger population with point attractors since the particles tend to be distributed over a larger area of the dynamical system. Note that we bin the population by position to compute mutual information between position in the 2d state space and time. Doing so reliably requires a smooth distribution after binning. For limit cycles, the particles usually stay close to attractor and thus provide sufficient count in each bin. However, for the point attractor, the population is usually spread over the entire 2d area between the two point attractors. Therefore, we need ${10}^5$ particles to get an accurate value of mutual information of point attractor system.

To find the optimal dynamical system geometry that operates with best accuracy when both internal noise and external noises are present, recall that we derived the following equations for strongly-driven limit cycles,

For a small $L/R$, we had found that $\mathrm{\Delta }\mathrm{\Phi }\sim {ϵ}_{ext}L/R$ where ${ϵ}_{ext}^2$ is a measure of the variance of light during the day. Further, we showed that $s^2-1\sim L/R$ in this limit. Hence, in the small $L/R$ (‘Stuart-Landau’) regime, the above equations reduce to,

The population variance when both noises are present is approximately given by $\text{max}({ϵ}_{int}^{2}R/L,{ϵ}_{ext}^{2}L/R)$. This variance is minimized when the two terms are equal, giving

which defines the optimal geometry of the dynamical system for given strength of internal and external noise.

In contrast, by taking the product of the equations above, we find the trade-off relationship,

The trade-off relationship above clarifies which parameters are held fixed and which ones are varied in our trade-off. If $Q$ is held fixed, this trade-off relationship holds under variations of all the parameters of the normal form of limit cycles (Equation 14). (While $L/R$ allows us to navigate the trade-off by increasing one of ${\sigma }_{ext}^2,{\sigma }_{int}^2$ and decreasing the other, other parameters such as the relaxation time leave both ${\sigma }_{ext}^2,{\sigma }_{int}^2$ relatively unaffected.)

Holding $Q$ fixed does involve holding the strength of external ${ϵ}_{ext}$ and internal ${ϵ}_{int}$ noise fixed. In all the models studied here, ${ϵ}_{ext}$ is simply defined as the size of the amplitude fluctuations in the external signal relative to the amplitude of the external signal itself - that is, the noise-to-signal ratio of the external signal - with no reference to the clock dynamics. Changing $L/R$ and other parameters can strengthen or weaken the coupling of this noisy external signal to the clock but do not affect the signal-to-noise ratio of the external signal itself.

Analogously, the strength of internal noise ${ϵ}_{int}^2$ is defined as the diffusion constant for the phase of a clock (e.g., in radians${}^2$) in the absence of an external driving signal. As discussed in Appendix 3, this definition ensures that limit cycles of different sizes develop the same phase variance over the same time when subject to the same ${ϵ}_{int}$.

For some purposes, it may make sense to hold the dimensionful diffusion constant $D_{int}={ϵ}_{int}^2{R}^2$ fixed while making comparisons. In this case, in addition to the trade-off effect discussed in this paper, large limit cycles are given an additional robustness to internal noise, trivially by virtue of their size, since the diffusion constant $D_{int}$ in chemical space is held fixed (instead of the dimensionless diffusion constant ${ϵ}_{int}^2$ for clock phase). In this case, it is insightful to re-write $Q={ϵ}_{ext}^2{ϵ}_{int}^2={ϵ}_{ext}^2{D}_{int}/R^2$ and recognize that $P=Q^{-1}$ is a measure of the power needed to maintain free running clock oscillations Cao et al., 2015 - larger cycles cost more energy per cycle to maintain. Thus, in this case, our trade-off should be understood as one at fixed power.

Modeling external noise (weather fluctuations) We generate a square wave of period $24$ hours to model the day-night cycle of light on Earth with the day length of 12 hr. However, such a square wave is modulated by weather fluctuations, for example, periods of reduced intensity due to passing clouds during the daytime. We model such fluctuating intensity as follows. We assume each weather condition lasts a random interval of time drawn from an exponential distribution of mean $2.4$ hrs (1/10 of a day). During a given weather condition, we set the intensity of light to a random value, drawn uniformly from $[1-{\text{noise}}_{\text{ext}},1]$ where ${\text{noise}}_{\text{ext}}$ is the strength of the external noise: $0$ means no external noise and $1$ means full external noise. This random value will range from 0 to 1 where one represents the maximum intensity during the day. (At night, the intensity is held at zero with no fluctuations). In the simulation of our limit cycle model, we set ${\text{noise}}_{\text{ext}}$ to 1. However, in our simulations for eight different models of biological clocks, ${\text{noise}}_{\text{ext}}$ ranges from 0.5 to 1 because when ${\text{noise}}_{\text{ext}}$ is too high, the system may not get entrained due to the difference in the natural and driving frequencies.

This noisy external signal is coupled to the diverse range of systems studied here in different ways as described in the respective sections. For each system, we simulate a population of organisms where each individual is subject to a different realization of the weather conditions described above. (a) In the Kai clock, the light signal is taken to affect the cellular ATP levels. (b) In the other eight diverse oscillators in the main paper, we coupled the light signal to the parameter specified as coupled to external signals in the original publications. (c) For the dynamical systems model, we assume that the position of the limit cycle is moved by the light signal. When the light intensity is reduced during the day to a value $\rho \in [0,1]$, we switch the dynamics to an alternative limit cycle (or point attractor) at a fractional distance $\rho$ between the ideal day and night cycles. For example, assume the night cycle is centered at $(0,0)$ and the day cycle is centered at $(-L,0)$. During a weather condition with intensity $\rho \in [0,1]$, we follow dynamics due to a limit cycle located at $(-\rho L,0)$. We follow the same rules for the point attractor.

Note that we have used a square wave to approximate the natural cycle of light on earth. The square wave also allows for an intuitive derivation of Equations 1, 2 by dividing up the day-night cycle into four parts: diffusion during the day and during the night, contracting variance during dawn and dusk. For other waveforms, such a clear separation is not possible and all these processes occur concurrently. However, Equations 1, 2 are expected to still hold up to $O(1)$ prefactors. Numerically, we tested sinusoidal inputs and verified our trade-off relationship.

Modeling internal noise The internal noise represents any source of stochasticity intrinsic to a single cell that would exist even in constant conditions. Such noise could be due to finite copy numbers of molecules, bursty of transcription etc.

We model internal noise in the Kai clock using explicit Gillespie simulations at finite copy number $N$ as described in the section on Kai clocks. For the diverse other biochemical clocks studied here, we add Langevin noise to the dynamical equations, following the prescriptions laid out in the original publications when available. In the dynamical systems models, we model internal noise by adding Langevin noise to the dynamical equations as described in the section on Langevin noise. Each individual particle in our simulation is subject to an independent random realizations of such Langevin noise. In order to ensure an apples-apples comparison between different clocks, we define the strength of the internal noise ${ϵ}_{int}$ to be the phase diffusion constant in undriven condition. See Appendix 3 on Dynamical systems simulations for more.

We develop and use two distinct measures of performance of noisy clocks driven by noisy inputs.

Mutual information: The performance of the clock is quantified by the mutual information between the clock state $\overrightarrow{c}$ and the time $t$,

for all $\overrightarrow{c}$ in the set of available positions $C$ and all $t$ in the available time bins $T$. (In the dynamical systems model, $\overrightarrow{c}$ represents the position in the 2d $r,t$ plane. For the explicit KaiABC biomolecular model, $\overrightarrow{c}$ represents the phosphorylation state of KaiC.) We simulate a population of clocks, where each clock is subject to a different realization of input signals, representing different weather conditions and also subject to different realizations of internal Langevin noise (or Gillespie fluctuations). We then collect the trajectories of each clock on the last day of the simulations and calculate the probability distribution $p(\overrightarrow{c}|t)$ of clock states at a given (objective) time $t\in [0,24]$ hrs of the last day in the simulation. The probability function $p(\overrightarrow{c})$ is calculated by accumulating the distribution of $p(\overrightarrow{c}|t)$ over time $t\in [0,24]$ hrs of the last day. The position $\overrightarrow{c}$ and time $t$ are binned into different bins depending on their values. We start the minimum and maximum values of the bins to the minimum and maximum values of the variables. The bin size in the time dimension is 0.48 hr or 28.8 min, while The bin size in the x and y dimensions are both $0.01$.

We refer to this mutual information measure as ‘Precision’ in Appendix 1—figure 1d, 2 and 6a from the main text.

Population variance along direction of motion: Mutual information is a good indicator of how well the clock encodes information about time. However, it is calculated for the entire day. Often, we want to see how the time-telling ability of a clock changes during the day (e.g., day vs night or before and after dusk). Hence we develop a new measure, closely related to mutual information, but can be computed at specific times of day.

Intuitively, the mutual information quantifies how much the population distributions of clock states at different times $t$ overlap. If these distributions are not overlapping, the clock state is a good readout of the time $t$. Such distributions are shown in Figures 4b and 5b (purple) in the main text.

We argue that only the spread of the clock distribution along the direction of motion of the clock in state space affects mutual information. The spread of the distribution in orthogonal directions does not affect mutual information as much.

To see this, we write mutual information between clock state ${\displaystyle \overrightarrow{c}}$ and time ${\displaystyle t}$ as,

Here $H(T)$ is a constant, independent of the clock mechanism. Thus, ${\displaystyle MI}$ depends entirely on the entropy of the distribution ${\displaystyle p(t|c)}$ of real times given clock state ${\displaystyle c}$, averaged over different clock states,

Consider a clock whose state-space is two dimensional with a periodic x-axis as shown in Appendix 4—figure 1. Further, assume that the distribution ${\displaystyle p(\overrightarrow{c}|t)}$ of clock states at a given time is supported on a rectangle of size ${\displaystyle a_x\times a_y}$ as shown in Appendix 4—figure 1 and that the clock states move along the x-axis at a uniform velocity ${\displaystyle u}$. This situation implies that

Appendix 4—figure 1

Download asset Open asset

So,

Since ${\displaystyle MI(C;T)=H(T)-H(T|C)}$, ${\displaystyle MI}$ depends on ${\displaystyle -\log a_x}$ and is independent of ${\displaystyle a_y}$, meaning that only the spread in the direction of motion ${\displaystyle a_x}$ affects the mutual information. Consequently, to understand the quality of time-telling at different times of the day, we project the population variance of ${\displaystyle p(\overrightarrow{c}|t)}$ to the direction of the instantaneous velocity of the center of mass of ${\displaystyle p(\overrightarrow{c}|t)}$.

Cramer-Rao (CR) bounds quantify the total available information about phase in a given length of history of the signal. Any estimator working with that length of history must necessarily have higher variance (i.e., higher error) than the Cramer- lower bound corresponding to that length of history. In the limit of infinitely long histories, the CR bound in this context corresponds to zero error; with any finite binning in time, the upper bound on MI is simply set by the number of bins in time. In our case, this bound is given by ${\displaystyle lo{g}_{2}50=5.64}$ bits. As shown in the main paper, as ${\displaystyle L/R\to 0}$, limit cycles process longer and longer histories of the external signal. Consequently, the mutual information for such cycles approaches the upper bound in the limit ${\displaystyle L/R\to 0}$ (assuming no internal noise) when computed with the same number of temporal bins (${\displaystyle 50}$ in this case).

During the daytime, sunlight intensity fluctuates because of cloud cover and we have referred to these fluctuations as external noise. In our simulations, we subject each individual in a population to a different realization of these weather conditions and compute the resulting population variation of clock state. Such variation limits the ability of the cell to read out the objective time from the clock state.

Here, we relate the population phase variance caused by random cloud cover in our dynamical systems model to the geometrically computed Phase Response Curve (PRC) due to a single dark pulse administered during the day. Using this geometric method, we will find that the ability of limit cycle to withstand external intensity fluctuations increases with ${\displaystyle R/L}$, the size ${\displaystyle R}$ of limit cycles relative to their separation ${\displaystyle L}$. In particular, we will show geometrically that the gain in phase variance during the day ${\displaystyle {\sigma }^2\Rightarrow {\sigma }^2+{\sigma }_{clouds}^2}$ scales as ${\displaystyle (L/R{)}^2}$, in perfect agree with stochastic weather simulations.

To compute the scaling relationship of ${\displaystyle {\sigma }_{clouds}^2}$, we compute the phase shift ${\displaystyle \mathrm{\Delta }\mathrm{\Phi }}$ caused by a single dark pulse with width ${\displaystyle \tau }$ on the limit cycles with angular speed ${\displaystyle \omega }$ (i.e., the Phase Response Curve (PRC) corresponding to such a dark pulse). Appendix 5— figure 2a shows an example of a dark pulse in the signal and how it affects the trajectory. Consider a clock at state ${\displaystyle \theta }$ on the day cycle. A dark pulse of length ${\displaystyle \tau }$ administered just then will change the dynamics to that of the night cycle. This clock has state ${\displaystyle \varphi =P(\theta )}$ with respect to the night cycle and will evolve for a time ${\displaystyle \tau }$ according to the night cycle dynamics, reaching a new state ${\displaystyle \varphi +\omega \tau }$, at a radial position determined by ${\displaystyle R,L}$. At the end of the dark pulse, we use the night-day circle map, ${\displaystyle \theta =Q(\varphi )}$, to find the clock state back on the day cycle. Note that all these shifts depend on the limit cycle geometry, that is, on ${\displaystyle R}$ and ${\displaystyle L}$, as shown in Appendix 5— figure 2. We can write each mapping using simple trigonometry:

Appendix 5—figure 1

Download asset Open asset

Appendix 5—figure 2

Download asset Open asset

and

Notice the mapping ${\displaystyle Q}$ only differs from ${\displaystyle P}$ by changing ${\displaystyle L}$ to ${\displaystyle -L}$. We also include the diagram showing the transition due to dark pulse in Appendix 5— figure 2. The process ‘1’ corresponds to ${\displaystyle \varphi =P(\theta )}$, ‘2’ corresponds to the rotation on the night cycle ${\displaystyle \varphi \to \varphi +\omega \tau }$, and ‘3’ corresponds to the transition back to the day cycle ${\displaystyle {\theta }^{\ast }=Q(\varphi +\omega \tau )}$. Combining this three processes, we write ${\displaystyle {\theta }^{\ast }}$ as ${\displaystyle {\theta }^{\ast }(\theta ,\tau ,L/R)}$ and expand it in the limit that ${\displaystyle L/R\Rightarrow 0}$ to obtain that

where ${\displaystyle \mathrm{\Delta }\mathrm{\Phi }={\theta }^{\ast }-(\theta +\omega \tau )}$ because ${\displaystyle \theta +\omega \tau }$ is the phase of the clock if it did not experience the dark pulse.

This expression ${\displaystyle \mathrm{\Delta }\mathrm{\Phi }}$ indicates the amount of phase shifted that the cloud causes. With different clocks experiencing different weather conditions, the variance gained among the population due to the fluctuation of sunlight grows like ${\displaystyle |\mathrm{\Delta }\mathrm{\Phi }{|}^2\sim (L/R{)}^2}$. We see good agreement between stochastic weather simulations and this geometric computation as shown in Appendix 5— figure 2d.

In this calculation, we focused on dark pulses administered at a fixed generic time (8 AM in Appendix 5— figure 2d). However, the PRC ${\displaystyle \mathrm{\Delta }\mathrm{\Phi }(\theta )}$ for dark pulses has a zero at a specific time of the day (see Appendix 5— figure 2c). That is, for each dark pulse of width ${\displaystyle \tau }$, there exists a time of administration such that ${\displaystyle \mathrm{\Delta }\mathrm{\Phi }=0}$! In fact, such a dark pulse has an entraining effect, reducing the population variance. We leave experimental and theoretical investigation of the counter-intuitive effects of such specially timed dark pulses to future work.

Here, we show that even if we include such dark pulses with an entraining effect, the variance gained at the end of the day is still proportional to ${\displaystyle (L/R{)}^2}$ in the limit that ${\displaystyle L/R}$ goes to zero. To simplify our derivation but retain the essence of what dark pulses do during the daytime, let’s us consider dark pulses coming at three times: in the morning (${\displaystyle \theta =-\pi /2}$), around noon (${\displaystyle \theta =-\omega \tau /2}$ with small ${\displaystyle \omega \tau }$), and in the evening (${\displaystyle \theta =\pi /2}$). Starting the day with variance ${\displaystyle {\sigma }_0^2}$, by the end of the day the variance becomes

Thus, the variance gained due to fluctuation, ${\displaystyle {\sigma }^2-{\sigma }_0^2={\sigma }_{clouds}^2}$, is proportional to ${\displaystyle (L/R{)}^2}$. This simple derivation may not rigorously reflect the correct constant in front of ${\displaystyle (L/R{)}^2}$ term, but the full rigorous derivation, concerning the dark pulses coming randomly at random time during the day, should yield the same power law dependent on ${\displaystyle L/R}$. Appendix 5— figure 2d shows that averaging ${\displaystyle \mathrm{\Delta }{\mathrm{\Phi }}^2}$ over pulses administered at different times numerically (dashed line) results in the same power law as for single pulses and as seen in stochastic weather simulations.

In our main paper, we claim that the variance of the clock state across a population drops ${\displaystyle {\sigma }^2\Rightarrow {\sigma }^2/s^2}$ at dusk where ${\displaystyle s^2-1\sim L/R}$ as ${\displaystyle L/R\Rightarrow 0}$. Data from Langevin simulations was presented. Here we will derive this result using a simple geometric argument about circle maps.

We define ${\displaystyle \varphi =P_T(\theta )}$ to be the phase on the night cycle that a clock evolves to, after time a time ${\displaystyle T}$, if the lights were suddenly turned off when the clock is at state ${\displaystyle \theta }$ on the day cycle. See Appendix 5—figure 1a,b. In principle, with complex relaxation dynamics between the limit cycles, ${\displaystyle P_T(\theta )}$ could show complex dependence on ${\displaystyle T}$. However, we work in a simplified model where the angular frequency of the clock is independent of the amplitude of oscillations. In this limit, ${\displaystyle T}$ only causes an overall shift in ${\displaystyle \varphi =P_T(\theta )}$; that is, we can write ${\displaystyle P_T(\theta )=P(\theta )+\omega T}$ where ${\displaystyle \omega }$ is the angular frequency of the clock. In what follows, we will be interested in the derivative of ${\displaystyle {\mathrm{\partial }}_{\theta }{P}_T(\theta )}$; hence we will work with ${\displaystyle P(\theta )}$ instead of ${\displaystyle P_T(\theta )}$.

This circle map, ${\displaystyle \varphi =P(\theta )}$, is important since it determines whether two differing day-time clock states are brought closer or taken further at dusk and thus determines the rate of entrainment of a population to the external signal. Consider two organisms that have nearby but distinct clock states ${\displaystyle {\theta }_0}$, ${\displaystyle {\theta }_0+\mathrm{\Delta }\theta }$ at dusk. After dusk, these two clocks will be mapped to ${\displaystyle P({\theta }_0)}$ and ${\displaystyle P({\theta }_0+\mathrm{\Delta }\theta )\approx P({\theta }_0)+\mathrm{\Delta }\theta dP(\theta )d\theta {|}_{\theta ={\theta }_0}}$ respectively. Thus, dusk changes the difference between the clock states from ${\displaystyle \mathrm{\Delta }\theta }$ to ${\displaystyle \mathrm{\Delta }\varphi }$ where,