Biophysical clocks face a tradeoff between internal and external noise resistance
Abstract
Many organisms use free running circadian clocks to anticipate the day night cycle. However, others organisms use simple stimulusresponse strategies (‘hourglass clocks’) and it is not clear when such strategies are sufficient or even preferable to free running clocks. Here, we find that free running clocks, such as those found in the cyanobacterium Synechococcus elongatus and humans, can efficiently project out light intensity fluctuations due to weather patterns (‘external noise’) by exploiting their limit cycle attractor. However, such limit cycles are necessarily vulnerable to ‘internal noise’. Hence, at sufficiently high internal noise, point attractorbased ‘hourglass’ clocks, such as those found in a smaller cyanobacterium with low protein copy number, Prochlorococcus marinus, can outperform free running clocks. By interpolating between these two regimes in a diverse range of oscillators drawn from across biology, we demonstrate biochemical clock architectures that are best suited to different relative strengths of external and internal noise.
https://doi.org/10.7554/eLife.37624.001eLife digest
The daily rising and setting of the sun is perhaps the most predictable pattern on Earth. Many organisms, from ancient bacteria to animals and plants, have evolved internal biological clocks to anticipate specific events such as dusk and dawn. However, biological clocks also need to continue working when faced with irregularities – both arising from within the organism and from external factors, such as a passing cloud that darkens the sky.
Some organisms, including humans, have a socalled ‘freerunning’ clock that generates a 24hour rhythm, and keeps ticking even in the absence of any time triggers. Others, such as certain cyanobacteria, have an ‘hourglass’ clock that is not selfsustained – rather, these clocks show a simple response to the sunrise (or sunset) that would gradually perish without another sunset (or sunrise).
So far, it has been unclear why organisms have different kinds of clocks and if one type of clock is better suited for some conditions than others. Here, Pittayakanchit, Lu et al. analyzed and compared mathematical models of clocks in a variety of organisms, from cyanobacteria and fungi to plants and animals.
The results revealed that internal and external irregularities put opposing pressures on biological clocks. Freerunning clocks are more precise and more robust to external fluctuations, but more susceptible to internal ones. In contrast, hourglass clocks can remain accurate when internal irregularities are high but can be disturbed by external ones.
Biological clocks affect the health of the entire organism and faulty clocks are implicated in numerous diseases. The study of Pittayakanchit, Lu et al. showed that the optimal architecture of a biological clock depends on the balance of irregularities in the external and internal environment of an organism. A next step will be to understand whether an organism can change its clock architecture while the environment changes. A better understanding of how biological clocks are regulated may help us find ways to tune faulty clocks to account for both the external environment and the internal state of an organism.
https://doi.org/10.7554/eLife.37624.002Introduction
Extracting information from a noisy external signal is fundamental to the survival of organisms in dynamic environments (Bowsher and Swain, 2014). From yeast anticipating the length of starvation (Mitchell et al., 2015) and bacteria estimating sugar availability (Tu et al., 2008), to dictyostelium counting cAMP pulses (Cai et al., 2014), organisms must often infer properties of the environment that are masked by noisy irregular fluctuations in order to be welladapted (Siggia and Vergassola, 2013; Mora and Wingreen, 2010).
A striking example of regularity in environmental stimuli is the daily daynight cycle of light on earth; organisms from all kingdoms of life use circadian clocks to synchronize  or ‘entrain’  in phase to these $24$hour periodic signals in order to anticipate and prepare for future changes in light (Winfree, 2001). The most remarkable and wellstudied examples of clocks are free running circadian clocks, found in organisms ranging from the cyanobacterium S. elongatus to insects, plants and humans. Such clocks use nonlinear dynamics to generate selfsustained $24$hr rhythms of a preferred amplitude even in the absence of external driving. Many salient properties have been ascribed to such free running internal rhythms (Troein et al., 2009; Winfree, 2001).
However, several organisms have only damped clocks or ‘hourglass clocks’; their response to daily changes in light is not a selfsustaining oscillation, but rather a physiological program that decays to a steady state over a day. For example, some strains of P. marinus, a smaller $0.5\mu m$ cyanobacterium with an estimated $50\times $ smaller protein copy number than S. elongatus (Bryant, 2003; Gutu et al., 2013; Holtzendorff et al., 2008; Dufresne et al., 2003; Kitayama et al., 2003), appear to have such a damped ‘hourglass’ clock, despite the clock being constituted from Kai proteins similar to those in S. elgonatus.
The potential benefits and drawbacks of these timing systems are not immediately obvious. In particular, it is unclear when an ‘hourglass’ clock might be sufficient or even preferred over free running clocks.
Here, we compare such classes of clocks when driven by the daynight cycle of light in fluctuating conditions. One source of fluctuations are amplitude fluctuations in the external daynight signal due to weather patterns (Gu et al., 2001) or other environmental disturbances. Phase entrainment to such fluctuating environmental signals is a challenge because while amplitude fluctuations are uninformative of phase, an entrainment mechanism looking for dawndusk transitions might conflate such amplitude fluctuations with true variations in phase. Biomolecular clocks also face an internal source of fluctuations (Lestas et al., 2010), due to various causes like finite copy number effects (Tsimring, 2014), bursty transcription, interactions with the cell cycle and cell division (Teng et al., 2013). It is clear that the inability to deal with either of these fluctuations will lead to poor phase entrainment, with a host of associated fitness costs in cyanobacteria (Woelfle et al., 2004), plants, rodents and humans (Evans and Davidson, 2013). However, it is not clear what kinds of clock architecture are best at dealing with internal and external fluctuations and whether these demands are compatible.
We find that free running clocks, based on limit cycle attractors, are a doubleedged sword when subject to such internally and externally fluctuating conditions. The flat direction along such continuous limit cycle attractors can selectively project out external amplitude fluctuations while retaining phase information. However, the flat direction along the attractor makes these continuous attractorbased clocks susceptible to internal fluctuations (e.g. low protein copy number [Potoyan and Wolynes, 2014]). In contrast, point attractorbased damped clocks are relatively resistant to internal fluctuations because they have no flat directions. Hence such ‘hourglass’ clocks outperform free running clocks at sufficiently high internal noise.
We first demonstrate our results in diverse biochemical oscillators, drawn from the literature (Leloup et al., 1999; Schmal et al., 2014; Locke et al., 2005; Leloup and Goldbeter, 2003; Goldbeter, 1991; Goodwin, 1965; Gonze and AbouJaoudé, 2013; Kondepudi and Prigogine, 2014; Elowitz and Leibler, 2000; Buşe et al., 2009; PotvinTrottier et al., 2016) on clocks in cyanobacteria, plants and mammals to cell cycle and synthetic oscillators. We complement this detailed networkbased study with dynamical systems theory that explains the same tradeoff in terms of the broad features common to the diverse models studied here. In all cases, our approach involves systematically deforming the dynamics to interpolate between free running and ‘hourglass’ clocks and using information theoretic measures to quantify clock performance in the presence of fluctuations.
By continuously interpolating between these clock architectures, our work predicts that a survey of clock systems in different environmental niches will reveal that clock architecture vary systematically with the relative strength of external and internal fluctuations (Laughlin, 1981). Further, our work suggests intriguing forward evolution experiments in the lab where the same structured external environment can nevertheless result in distinct regulatory systems, depending on the size of internal fluctuations. Finally, the existence of ‘hourglass’ clocks are easier to overlook experimentally than free running oscillations. Hence our theoretical demonstration that ‘hourglass’ clocks have functional benefits in specific conditions highlights the importance of experiments that specifically look for such damped clocks. More broadly, our work highlights the need to experimentally probe regulatory strategies by varying different kinds of noise independently when possible, since the strategies to deal with different kinds of noise are not equivalent and can be in conflict.
Results
Free running clocks and damped ‘hourglass’ clocks
Many organisms like humans and rodents have free running clocks that show selfsustained 24 hr rhythms even in constant dark or light conditions. A particularly simple and wellcharacterized free running clock is that found in S. elongatus where the clock dynamics can be reproduced by the posttranslational 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.
Not all organisms have a freerunning 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 daynight cycling (see Figure 1b). In fact, a sister cyanobacterial species P. marinus has a KaiBCprotein based clock. While the details of this clock are not fully characterized, the clock lacks the KaiAmediated 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 daynight 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.
External fluctuations
The daynight 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.
Internal 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).
Noise resistance of Kaibased clocks
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 posttranslational 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 posttranslational 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 daynight 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.
Noise resistance in other biochemical clocks
While our study here was motivated by the contrasting Kai proteinbased 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 wellstudied 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 AbouJaoudé, 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.
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 AbouJaoudé, 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 tradeoff 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 tradeoff 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.
Dynamical systems theory of noise resistance
We have demonstrated a tradeoff between external and internal noise resistance in diverse clocks. While it is possible to trace the origin of this tradeoff 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 tradeoff emerges due to basic features of free running and damped clock dynamics and can thus be expected to hold broadly.
Limit cycle clocks and point attractor clocks
Free running clocks are phenomenologically welldescribed by a limit cycle attractor, a nonlinear 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.
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 freerunning 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 freerunning oscillations while $L$ represents the strength or amplitude of the external driving signal.
In contrast, damped clocks are phenomenologically welldescribed by a daytime and a nighttime 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$).
External noise
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 timescales 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.
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$).
Internal noise
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 ${\u03f5}_{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 daynight 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}+{\u03f5}_{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 steadystate phase variance (${\sigma}^{2}=(({\sigma}^{2}+{\u03f5}_{int}^{2}{T}_{day})/{s}^{2}+{\u03f5}_{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 ${\u03f5}_{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 daynight 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’ offattractor 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.
Combination of external and 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 ${\u03f5}_{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,
In the SI, also we show that, under certain simplifying assumptions, Equations 1 and 2 can be combined to give an explicit tradeoff relationship,
where $Q={\u03f5}_{int}^{2}{\u03f5}_{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 tradeoff explicit and also clarifies which parameters are varied and which parameters are held fixed in this tradeoff. 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 tradeoff does assume that the strength of the external fluctuations ${\u03f5}_{ext}$ – that is, the fractional size of amplitude fluctuations in the external signal – is held fixed. Similarly, we hold ${\u03f5}_{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.
Speedprecision tradeoff
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 speedprecision tradeoff 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}={\u03f5}_{int}/{\u03f5}_{ext}$  will have the highest precision.
Discussion
Free running circadian clocks are a remarkable result of evolution in a changing but predictable environment and are thought to provide numerous benefits (Troein et al., 2009). Here, we showed that the limit cycle attractor underlying such a clock is able to effectively project out weatherrelated amplitude changes that are perpendicular to the flat direction of the attractor. Similar roles for the flat direction of continuous attractors in projecting out external (or input) fluctuations have been explored in neuroscience (Burak and Fiete, 2012); Seung (1996), for example, for head and eye motor control and spatial navigation. However, we also found that the same flat direction becomes a vulnerability with internal fluctuations since such fluctuations cannot be restricted to be perpendicular to the attractor.
We confirmed the tradeoff between resistance to external and internal noise in diverse models of biochemical clocks and oscillators, using parameters inferred from experimental data by the original publications (Leloup et al., 1999; Schmal et al., 2014; Locke et al., 2005; Leloup and Goldbeter, 2003; Goldbeter, 1991; Goodwin, 1965; Gonze and AbouJaoudé, 2013; Kondepudi and Prigogine, 2014; Elowitz and Leibler, 2000; Buşe et al., 2009). The tradeoff in each of these models can be given explanations that are specific to those models; for example, one can identify specific bottlenecks for external and internal noise in these models (Cheong et al., 2011). However, we have provided an alternative kind of analysis based on the geometry of the dynamical systems involved. Such an explanation misses aspects specific to each clock  for example, how specific biologically tuneable parameters in each model affect internal and external noise resistance. However, the dynamical systems picture has the advantage in that it identifies the common origin of the tradeoff across these systems. Such a dynamical systems picture has been fruitful in making other general but falsifiable predictions in biology (Gan and O'Shea, 2017; Leypunskiy et al., 2017; Corson and Siggia, 2017), going back to Winfree’s phase singularity (Winfree, 2001).
Our dynamical systems theory shows that the critical parameter for noise resistance is the strength of the external driving relative to the amplitude of free running oscillations, captured by the geometric ratio $L/R$ in our analysis. Weak driving provides resistance to external noise while strong driving provides resistance to internal noise. While our dynamical systems theory involve planar circular limit cycles, the models in Figure 2 have complex nonplanar noncircular limit cycles and yet reproduce our tradeoff. Finally, while the internal noise discussed here is set by finite copy number, this dependence is not essential to the results here. Any source of disturbance (e.g., bursty transcription) that perturbs the phase of the oscillator in constant light conditions is equivalent to internal noise. Similarly, external noise can involve any kind of fluctuation (e.g., multiplicative fluctuations, phase fluctuations) of the external signal that does not result in a persistent phase shift of the external signal itself.
Our work suggests that the damped oscillators are not merely poor cousins of the remarkable free running oscillators found for example, in S. elongatus. At the low protein copy numbers such as those found in P. marinus, damped clocks keep time more reliably than free running clocks. (Low copy number has been linked to a trend towards reduced genome size and cell size in P. marinus [Bryant, 2003].) In addition to the noisy internal environment of P. marinus, the external environment might also play a role in selecting a damped clock; P. marinus is typically found in the open ocean, where the external environment may be more regular than the fresh water habitat of S. elongatus. In addition to P. marinus, damped oscillators are found elsewhere in biology, often in specific physiological conditions (Saunders, 2002; Kidd et al., 2015). Understanding the benefits and drawbacks of such damped oscillators in different conditions is critical since such oscillations are easily overlooked experimentally, in comparison to free running oscillations.
While numerous upstream and downstream considerations can modify (Rand et al., 2004; Domijan and Rand, 2011) the ultimate biological impact of fluctuations, we find that the core oscillator’s geometry in itself can continuously trade off protection against external fluctuations for protection against internal fluctuations in the diverse range of models studied here.
Note added in proofs: The study of Monti et al. (2018, in press) independently arrived at the conclusion that free running clocks based on limitcycles are more robust to external noise. Experiments in Chew et al. (2018, in press) suggest that the free running clock in S. elongatus is less robust to internal noise than the hourglass clock in P. marinus.
Materials and methods
We incorporated most of our methods in Results and Discussion. For details beyond those presented in Results, please see Appendices. Code to simulate the systems can be found at https://github.com/WeerapatP/elife_tradeoff_clocks (Pittayakanchit, 2018; copy archived at https://github.com/elifesciencespublications/elife_tradeoff_clocks).
Appendix 1
Tradeoff in Kaibased clocks
We demonstrate our tradeoff using Gillespie simulations of an explicit biomolecular KaiABC model of the posttranslational clocks in S. elongatus and P. marinus.
S. elongatus clock  hexamers with collective KaiA feedback
The S. elongatus clock has been wellcharacterized 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 KaiAdependent which allows for feedback that enables collective coherent oscillations in a cell. After complete phosphorylation of KaiAC complexes (usually by the end of the day), KaiC forms a KaiBC complex which then dephosphorylates in an ordered manner. Crucially, the KaiBC 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).
P. marinus model  independent hexamers
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.
Hybrid model
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 KaiAdependent phosphorylation pathway, much like in S. elongatus, that is used during the day and driven forward by ATP.
But to also include P. marinuslike 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. elongatuslike KaiA dependent pathway is accessible. When $\gamma =0$, only the P. marinuslike 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 KaiBbased 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 freerunning limit cycle oscillations in S. elongatus.
However, as $\gamma \to 0$, the KaiAindependent 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 cyclebased free running clocks and one pointattractor 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$.
Gillespie simulations
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: $dt=0.01\phantom{\rule{thickmathspace}{0ex}}hr,$ ${k}_{+}={k}_{}=2m\cdot 0.04932\phantom{\rule{thickmathspace}{0ex}}h{r}^{1},$ ${k}_{Aon}=0.2466\phantom{\rule{thickmathspace}{0ex}}\mu {M}^{1}h{r}^{1},$ ${k}_{Aoff}=0.02466\phantom{\rule{thickmathspace}{0ex}}h{r}^{1},$ ${k}_{C\to C\ast}=0.2466\phantom{\rule{thickmathspace}{0ex}}h{r}^{1},$ ${k}_{ABC}=123.30\phantom{\rule{thickmathspace}{0ex}}h{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 lightdark 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 uppathways and $d$ is the net phosphorylation state of KaiC in the KaiBbound ‘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 timeofday was binned with bins of size $0.5$ hrs.
Phase portrait
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.
Appendix 2
Other oscillators
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 tradeoff 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 attractorbased ‘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 daynight 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 ${\u03f5}_{ext}$, defined as the noisetosignal ratio of the amplitude fluctuations in the external signal, fixed. We varied internal noise ${\u03f5}_{int}$ along the $x$ axis of plots in Figure 2 of the main paper. Here, ${\u03f5}_{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.
Neurospora and Drosophila circadian clocks by Goldbeter
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
Arabidopsis circadian clock by Millar et al
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.
B.3 Mammalian PerCry circadian clock by Leloup et al
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.
B.4 cdc2cyclin cell cycle by Goldbeter
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 selfsustained 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.
B.5 Goodwin oscillator
One of the earliest models of biochemical oscillators was proposed by Goodwin (1965) (later corrected). We use the simplest widelystudied version of such a Goodwin oscillator (Gonze and AbouJaoudé, 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 $\alpha =100$ are taken to the output of the clock for computing MI.
B.6 Repressilator
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 nonlinearity 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 nondriven 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 $s(t)=0$ are taken to the output of the clock for computing $MI$.
B.7 Brusselator
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$.
Appendix 3
Dynamical Systems
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}{\u03f5}_{int}^{2}$ is the diffusion constant, d$W$ is a Wiener process, and here we assume that the internal noise is a homogeneous whitenoise in the 2dimensional 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 xaxis 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 StuartLandau 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 daynight 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 daynight cycle.
Simulations
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 ${\u03f5}_{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 ${\u03f5}_{int}^{2}$ is a measure of phase diffusion, independent of limit cycle size. In other words, ${\u03f5}_{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 ${\u03f5}_{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 ${\u03f5}_{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 steadystate 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 speederror 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 nighttime 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.
Optimal Dynamical system and tradeoff
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 stronglydriven limit cycles,
For a small $L/R$, we had found that $\mathrm{\Delta}\mathrm{\Phi}\sim {\u03f5}_{ext}L/R$ where ${\u03f5}_{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$ (‘StuartLandau’) regime, the above equations reduce to,
The population variance when both noises are present is approximately given by $\text{max}({\u03f5}_{int}^{2}R/L,{\u03f5}_{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 tradeoff relationship,
The tradeoff relationship above clarifies which parameters are held fixed and which ones are varied in our tradeoff. If $Q$ is held fixed, this tradeoff 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 tradeoff 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 ${\u03f5}_{ext}$ and internal ${\u03f5}_{int}$ noise fixed. In all the models studied here, ${\u03f5}_{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 noisetosignal 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 signaltonoise ratio of the external signal itself.
Analogously, the strength of internal noise ${\u03f5}_{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 ${\u03f5}_{int}$.
For some purposes, it may make sense to hold the dimensionful diffusion constant ${D}_{int}={\u03f5}_{int}^{2}{R}^{2}$ fixed while making comparisons. In this case, in addition to the tradeoff 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 ${\u03f5}_{int}^{2}$ for clock phase). In this case, it is insightful to rewrite $Q={\u03f5}_{ext}^{2}{\u03f5}_{int}^{2}={\u03f5}_{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 tradeoff should be understood as one at fixed power.
Appendix 4
Supplementary Methods
Modeling external noise (weather fluctuations) We generate a square wave of period $24$ hours to model the daynight 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 daynight 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 tradeoff 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 applesapples comparison between different clocks, we define the strength of the internal noise ${\u03f5}_{int}$ to be the phase diffusion constant in undriven condition. See Appendix 3 on Dynamical systems simulations for more.
Measures of clock timetelling quality
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 timetelling 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 $\overrightarrow{c}$ and time $t$ as,
Here $H(T)$ is a constant, independent of the clock mechanism. Thus, $MI$ depends entirely on the entropy of the distribution $p(tc)$ of real times given clock state $c$, averaged over different clock states,
Consider a clock whose statespace is two dimensional with a periodic xaxis as shown in Appendix 4—figure 1. Further, assume that the distribution $p(\overrightarrow{c}t)$ of clock states at a given time is supported on a rectangle of size $a}_{x}\times {a}_{y$ as shown in Appendix 4—figure 1 and that the clock states move along the xaxis at a uniform velocity $u$. This situation implies that
So,
Since $MI(C;T)=H(T)H(TC)$, $MI$ depends on $\mathrm{log}{a}_{x}$ and is independent of $a}_{y$, meaning that only the spread in the direction of motion $a}_{x$ affects the mutual information. Consequently, to understand the quality of timetelling at different times of the day, we project the population variance of $p(\overrightarrow{c}t)$ to the direction of the instantaneous velocity of the center of mass of $p(\overrightarrow{c}t)$.
CramerRao bounds
CramerRao (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 $lo{g}_{2}50=5.64$ bits. As shown in the main paper, as $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 $L/R\to 0$ (assuming no internal noise) when computed with the same number of temporal bins ($50$ in this case).
Appendix 5
Circle Map  Dark pulse phase shift
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 $R/L$, the size $R$ of limit cycles relative to their separation $L$. In particular, we will show geometrically that the gain in phase variance during the day $\sigma}^{2}\Rightarrow {\sigma}^{2}+{\sigma}_{clouds}^{2$ scales as $(L/R{)}^{2}$, in perfect agree with stochastic weather simulations.
To compute the scaling relationship of $\sigma}_{clouds}^{2$, we compute the phase shift $\mathrm{\Delta}\mathrm{\Phi}$ caused by a single dark pulse with width $\tau$ on the limit cycles with angular speed $\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 $\theta$ on the day cycle. A dark pulse of length $\tau$ administered just then will change the dynamics to that of the night cycle. This clock has state $\varphi =P(\theta )$ with respect to the night cycle and will evolve for a time $\tau$ according to the night cycle dynamics, reaching a new state $\varphi +\omega \tau$, at a radial position determined by $R,L$. At the end of the dark pulse, we use the nightday circle map, $\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 $R$ and $L$, as shown in Appendix 5— figure 2. We can write each mapping using simple trigonometry:
and
Notice the mapping $Q$ only differs from $P$ by changing $L$ to $L$. We also include the diagram showing the transition due to dark pulse in Appendix 5— figure 2. The process ‘1’ corresponds to $\varphi =P(\theta )$, ‘2’ corresponds to the rotation on the night cycle $\varphi \to \varphi +\omega \tau$, and ‘3’ corresponds to the transition back to the day cycle ${\theta}^{\ast}=Q(\varphi +\omega \tau )$. Combining this three processes, we write $\theta}^{\ast$ as ${\theta}^{\ast}(\theta ,\tau ,L/R)$ and expand it in the limit that $L/R\Rightarrow 0$ to obtain that
where $\mathrm{\Delta}\mathrm{\Phi}={\theta}^{\ast}(\theta +\omega \tau )$ because $\theta +\omega \tau$ is the phase of the clock if it did not experience the dark pulse.
This expression $\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 $\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 $\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 $\tau$, there exists a time of administration such that $\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 counterintuitive 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 $(L/R{)}^{2}$ in the limit that $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 ($\theta =\pi /2$), around noon ($\theta =\omega \tau /2$ with small $\omega \tau$), and in the evening ($\theta =\pi /2$). Starting the day with variance $\sigma}_{0}^{2$, by the end of the day the variance becomes
Thus, the variance gained due to fluctuation, $\sigma}^{2}{\sigma}_{0}^{2}={\sigma}_{clouds}^{2$, is proportional to $(L/R{)}^{2}$. This simple derivation may not rigorously reflect the correct constant in front of $(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 $L/R$. Appendix 5— figure 2d shows that averaging $\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.
Circle Map  Step Response Curve
In our main paper, we claim that the variance of the clock state across a population drops $\sigma}^{2}\Rightarrow {\sigma}^{2}/{s}^{2$ at dusk where ${s}^{2}1\sim L/R$ as $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 $\varphi ={P}_{T}(\theta )$ to be the phase on the night cycle that a clock evolves to, after time a time $T$, if the lights were suddenly turned off when the clock is at state $\theta$ on the day cycle. See Appendix 5—figure 1a,b. In principle, with complex relaxation dynamics between the limit cycles, ${P}_{T}(\theta )$ could show complex dependence on $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, $T$ only causes an overall shift in $\varphi ={P}_{T}(\theta )$; that is, we can write ${P}_{T}(\theta )=P(\theta )+\omega T$ where $\omega$ is the angular frequency of the clock. In what follows, we will be interested in the derivative of ${\mathrm{\partial}}_{\theta}{P}_{T}(\theta )$; hence we will work with $P(\theta )$ instead of ${P}_{T}(\theta )$.
This circle map, $\varphi =P(\theta )$, is important since it determines whether two differing daytime 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 $\theta}_{0$, ${\theta}_{0}+\mathrm{\Delta}\theta$ at dusk. After dusk, these two clocks will be mapped to $P({\theta}_{0})$ and $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 $\mathrm{\Delta}\theta$ to $\mathrm{\Delta}\varphi$ where,
By a similar argument, if the phase variance of clock states across a population is $\sigma}^{2$ before dusk, it will be reduced by,
This expression is valid in the regime where the population variance $\sigma}^{2$ is small enough to linearize the circle map $P(\theta )$. Similar considerations apply to the dawn transition between the night and day cycle as well. Both circle maps were recently experimentally characterized for S. elongatus in Leypunskiy et al. (2017).
In our simple theoretical model where clock frequency does not change with amplitude (i.e. the radial coordinate), we can easily compute $P(\theta )$ from geometry. In Appendix 5—figure 1, we draw a diagram of the transition from a particle on the day cycle at the phase $\theta$ to the night cycle at the phase $\varphi$. By trigonometry, we write
and derive
where $\theta$ corresponds to the angle on the day cycle at dusk, which is at $\pi /2$ in Appendix 5—figure 1a. This equation implies that as the day and night limit cycle gets closer, the geometric focusing effect $s$ converges to one. This asymptotic behavior is intuitive because if $L=0$, meaning no transition, then the variance should remain the same ($s=1$, so $\sigma}^{2}\to {\sigma}^{2}/{1}^{2$ at the transition).
Remarkably, our geometric derivation of ${s}^{2}1$ matches the variance drop $\sigma}^{2}\to {\sigma}^{2}/{s}^{2$ seen in stochastic simulations of weather conditions; see Appendix 5—figure 1b. The variance gain during the day is the result of the fluctuation of sunlight, simulated as random dark pulses of random intervals, amplitude and time of delivery. Such variance is accumulated during the day and the drop over dusk time is measured (green Xs).
Appendix 5—figure 1e shows the variance drop seen in simulations with internal noise in Langevin simulations. While the cause of variance increase during the day is different (finite copy number effects), the variance drop at dusk agrees well with the geometric computation of $s}^{2$ and thus with the external noise simulations as well. In both cases, the simulations and geometric theory show that ${s}^{2}1\sim L/R$ as $L/R\Rightarrow 0$.
Appendix 6
Langevin model of finite copy number fluctuations
Chemical reactions that occur in the bulk of a homogeneous solution can be described by a set of ordinary differential equations. However, within a single cell the copy number of molecule is limited and thus the reaction carries internal noise from the stochastic fluctuations. Gillespie showed that chemical reactions under finite copy number can be approximated by a Langevin dynamics using the following argument (Gillespie, 2007),
Consider an elementary reaction
with the forward rate constant ${k}_{+}$, during each infinitesimal time $\delta t$, the probability of the occurrence of this reaction follows a Poisson distribution whose mean and variance both equal to ${R}_{+}\delta t={k}_{+}\cdot {N}_{A}\cdot {N}_{B}\cdot \delta t$. Integration over a larger time step, the Poisson distribution can be approximated into a Gaussian form, resulting in Langevin dynamics,
where $W$ is a standard Wiener process of mean $0$ and autocorrelation function $\u27e8W({t}_{1})W({t}_{2})\u27e9=\delta ({t}_{1}{t}_{2})$.
To describe a chemical reaction network, the Langevin equation for each species consists of contributions to the noise from each reaction where the species is involved. Now consider adding another reaction
with the rate constant ${k}_{}$, then the Langevin equation for species A becomes,
where $R}_{+}={k}_{+}\cdot {N}_{A}\cdot {N}_{B$ and $R}_{}={k}_{}\cdot {N}_{C}\cdot {N}_{D$ respectively denote the number rates of the forward and the backward reaction; $\mathrm{d}{W}_{1}$ and $\mathrm{d}{W}_{2}$ are identical independent standard Wiener processes.
To fully determine the effect of the noise using the Langevin dynamics for a chemical reaction network, one needs to consider all of the reactions corresponding to the species of interest; the noise term usually becomes timedependent and multiplicative. To simplify the description of internal noise in our phenomenological model of limit cycle/point attractor, we take a first order approximation that the diffusion coefficient in the reaction coordinate space is homogeneous in both space and time. (See similar treatments of another biological system in Potoyan and Wolynes (2014). In contrast, our explicit KaiABC simulations, as well as numerical simulations on the other types of biooscillators, presented later, do not make this simplifying assumption of homogeneous diffusion.) This allows us to write a 2dimension phenomenological stochastic differential equation
where the $f(\overrightarrow{z},t)$ denotes the deterministic dynamics driven by daynight cycles and the diffusion constant $D$ is assumed to be inversely proportional to the total number of KaiC molecules within the cell. For limit cycles of radius $R$, we set $D\sim {R}^{2}{\u03f5}_{int}^{2}$. Then, $\u03f5}_{int}^{2$ is the diffusion constant for the phase of the oscillator. We hold $\u03f5}_{int}^{2$ fixed while changing $R$ to make a fair comparison across systems of different size.
Population variance
For the cell to carry out a reliable computation, the population variance from the internal noise needs to be reduced. Such noise reduction comes from the dynamics of the attractor. In the limit cycle attractor mechanism, the internal noise reduction is performed only along the radial axis but not along the flat attractor direction.
In contrast, the point attractor mechanism is able to limit population variance due to internal noise in all directions due to the effective ‘curvature’ of the dynamics. Here we analytically estimate the steadystate population variance for a point attractor mechanism. The population variance is together determined by the diffusive term $\sqrt{2D}\cdot \mathrm{d}\overrightarrow{W}$, and the noise reduction effect from the restoring force of the point attractor’s harmonic well. During each infinitesimal time $\delta t$, the internal noise increase the variance by
In contrast, the overdamped deterministic motion within a harmonic well provides a focusing effect that reduces the variance exponentially with time. To quantify this focusing effect, consider a 1d overdamped dynamics of a particle within a harmonic energy well of $V(r)=k\cdot {r}^{2}$. The solution to the equation of motion is $r(t)={r}_{0}\cdot {e}^{2kt}$, with initial position $r(0)={r}_{0}$. Consider an ensemble of points with a mean initial position $\mu}_{0$ and a initial variance of $\sigma}_{0}^{2$, one can solve the dynamics of the mean as
and the dynamics of the variance as
Thus, per infinitesimal time $\delta t$, the geometric focusing effect of the energy well of the point attractor reduces the population variance by
where $g={e}^{4k\delta t}$.
Under the competition between the spreading effect from the internal noise and the geometrical focusing effect from the deterministic dynamics, the population variance reaches a steady value solved by
and by taking the limit of $\delta t$ goes to $0$, we have ${\sigma}_{st}^{2}=D/2k$.
Data availability
The code and data used in the simulations are available via GitHub at https://github.com/WeerapatP/elife_tradeoff_clocks
References

Environmental sensing, information transfer, and cellular decisionmakingCurrent Opinion in Biotechnology 28:149–155.https://doi.org/10.1016/j.copbio.2014.04.010

Existence of limit cycles in the repressilator equationsInternational Journal of Bifurcation and Chaos 19:4097–4106.https://doi.org/10.1142/S0218127409025237

The free energy cost of accurate biochemical oscillationsNature Physics 11:772–778.https://doi.org/10.1038/nphys3412

Health consequences of circadian disruption in humans and animal modelsProgress in Molecular Biology and Translational Science 119:283–323.https://doi.org/10.1016/B9780123969712.000105

Stochastic simulation of chemical kineticsAnnual Review of Physical Chemistry 58:35–55.https://doi.org/10.1146/annurev.physchem.58.032806.104637

Circadian rhythms and molecular noiseChaos: An Interdisciplinary Journal of Nonlinear Science 16:026110.https://doi.org/10.1063/1.2211767

Robustness of circadian rhythms with respect to molecular noiseProceedings of the National Academy of Sciences 99:673–678.https://doi.org/10.1073/pnas.022628299

Oscillatory behavior in enzymatic control processesAdvances in Enzyme Regulation 3:425–437.https://doi.org/10.1016/00652571(65)900671

Cloud modulation of surface solar irradiance at a pasture site in southern brazilAgricultural and Forest Meteorology 106:117–129.https://doi.org/10.1016/S01681923(00)002094

Stoichiometric interactions between cyanobacterial clock proteins KaiA and KaiCBiochemical and Biophysical Research Communications 316:195–202.https://doi.org/10.1016/j.bbrc.2004.02.034

BookModern Thermodynamics: From Heat Engines to Dissipative StructuresJohn Wiley & Sons.

A simple coding procedure enhances a neuron’s information capacityZeitschrift Fur Naturforschung. Section C, Biosciences 36:910–912.

Limit cycle models for circadian rhythms based on transcriptional regulation in Drosophila and NeurosporaJournal of Biological Rhythms 14:433–448.https://doi.org/10.1177/074873099129000948

Modelling genetic networks with noisy and varied experimental data: the circadian clock in Arabidopsis thalianaJournal of Theoretical Biology 234:383–393.https://doi.org/10.1016/j.jtbi.2004.11.038

Seasonal variation in radiative and turbulent exchange at a deciduous forest in central massachusettsJournal of Applied Meteorology 35:122–134.https://doi.org/10.1175/15200450(1996)035<0122:SVIRAT>2.0.CO;2

Limits of sensing temporal concentration changes by single cellsPhysical Review Letters 104:248101.https://doi.org/10.1103/PhysRevLett.104.248101

Design principles underlying circadian clocksJournal of the Royal Society Interface 1:119–130.https://doi.org/10.1098/rsif.2004.0014

Global entrainment of transcriptional systems to periodic inputsPLoS Computational Biology 6:e1000739.https://doi.org/10.1371/journal.pcbi.1000739

Modeling and simulating the Arabidopsis thaliana circadian clock using XPPAUTOMethods in Molecular Biology 1158:337–358.https://doi.org/10.1007/9781493907007_23

Weather and seasons together demand complex biological clocksCurrent Biology 19:1961–1964.https://doi.org/10.1016/j.cub.2009.09.024

Noise in biologyReports on Progress in Physics 77:026601.https://doi.org/10.1088/00344885/77/2/026601

BookThe Geometry of Biological TimeSpringer Science & Business Media.https://doi.org/10.1007/9781475734843

The adaptive value of circadian clocks: an experimental assessment in cyanobacteriaCurrent Biology : CB 14:1481–1486.https://doi.org/10.1016/j.cub.2004.08.023

Robust circadian clocks from coupled proteinmodification and transcriptiontranslation cyclesProceedings of the National Academy of Sciences 107:22540–22545.https://doi.org/10.1073/pnas.1007613107
Article and author information
Author details
Funding
Simons Foundation
 Arvind Murugan
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Aaron Dinner, John Hopfield, Eugene Leypunskiy, Charles Matthews, Brian Moths, Thomas Witten, and the Rust and Murugan labs for fruitful discussions.
Version history
 Received: April 17, 2018
 Accepted: June 23, 2018
 Accepted Manuscript published: July 10, 2018 (version 1)
 Version of Record published: July 25, 2018 (version 2)
Copyright
© 2018, Pittayakanchit et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 3,651
 views

 605
 downloads

 33
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Cell Biology
 Computational and Systems Biology
Bats have unique characteristics compared to other mammals, including increased longevity and higher resistance to cancer and infectious disease. While previous studies have analyzed the metabolic requirements for flight, it is still unclear how bat metabolism supports these unique features, and no study has integrated metabolomics, transcriptomics, and proteomics to characterize bat metabolism. In this work, we performed a multiomics data analysis using a computational model of metabolic fluxes to identify fundamental differences in central metabolism between primary lung fibroblast cell lines from the black flying fox fruit bat (Pteropus alecto) and human. Bat cells showed higher expression levels of Complex I components of electron transport chain (ETC), but, remarkably, a lower rate of oxygen consumption. Computational modeling interpreted these results as indicating that Complex II activity may be low or reversed, similar to an ischemic state. An ischemiclike state of bats was also supported by decreased levels of central metabolites and increased ratios of succinate to fumarate in bat cells. Ischemic states tend to produce reactive oxygen species (ROS), which would be incompatible with the longevity of bats. However, bat cells had higher antioxidant reservoirs (higher total glutathione and higher ratio of NADPH to NADP) despite higher mitochondrial ROS levels. In addition, bat cells were more resistant to glucose deprivation and had increased resistance to ferroptosis, one of the characteristics of which is oxidative stress. Thus, our studies revealed distinct differences in the ETC regulation and metabolic stress responses between human and bat cells.

 Computational and Systems Biology
 Neuroscience
Normal aging leads to myelin alterations in the rhesus monkey dorsolateral prefrontal cortex (dlPFC), which are positively correlated with degree of cognitive impairment. It is hypothesized that remyelination with shorter and thinner myelin sheaths partially compensates for myelin degradation, but computational modeling has not yet explored these two phenomena together systematically. Here, we used a twopronged modeling approach to determine how agerelated myelin changes affect a core cognitive function: spatial working memory. First, we built a multicompartment pyramidal neuron model fit to monkey dlPFC empirical data, with an axon including myelinated segments having paranodes, juxtaparanodes, internodes, and tight junctions. This model was used to quantify conduction velocity (CV) changes and action potential (AP) failures after demyelination and subsequent remyelination. Next, we incorporated the single neuron results into a spiking neural network model of working memory. While complete remyelination nearly recovered axonal transmission and network function to unperturbed levels, our models predict that biologically plausible levels of myelin dystrophy, if uncompensated by other factors, can account for substantial working memory impairment with aging. The present computational study unites empirical data from ultrastructure up to behavior during normal aging, and has broader implications for many demyelinating conditions, such as multiple sclerosis or schizophrenia.