Abstract

Many organisms use free running circadian clocks to anticipate the day night cycle. However, others organisms use simple stimulus-response 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 attractor-based ‘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.001

eLife 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 so-called ‘free-running’ clock that generates a 24-hour rhythm, and keeps ticking even in the absence of any time triggers. Others, such as certain cyanobacteria, have an ‘hourglass’ clock that is not self-sustained – 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. Free-running 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.002

Introduction

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 well-adapted (Siggia and Vergassola, 2013; Mora and Wingreen, 2010).

A striking example of regularity in environmental stimuli is the daily day-night 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 well-studied 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 non-linear dynamics to generate self-sustained 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 self-sustaining 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μm cyanobacterium with an estimated 50× 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 day-night cycle of light in fluctuating conditions. One source of fluctuations are amplitude fluctuations in the external day-night 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 dawn-dusk 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 double-edged 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 attractor-based clocks susceptible to internal fluctuations (e.g. low protein copy number [Potoyan and Wolynes, 2014]). In contrast, point attractor-based damped clocks are relatively resistant to internal fluctuations because they have no flat directions. Hence such ‘hourglass’ clocks out-perform 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 Abou-Jaoudé, 2013; Kondepudi and Prigogine, 2014; Elowitz and Leibler, 2000; Buşe et al., 2009; Potvin-Trottier et al., 2016) on clocks in cyanobacteria, plants and mammals to cell cycle and synthetic oscillators. We complement this detailed network-based study with dynamical systems theory that explains the same trade-off 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 self-sustained 24 hr rhythms even in constant dark or light conditions. A particularly simple and well-characterized free running clock is that found in S. elongatus where the clock dynamics can be reproduced by the post-translational dynamics of Kai ABC in vivo as well. Measuring the phosphorylation state at any one of several sites on KaiC reveals an orderly phosphorylation reaction with a period of 24 hr. As shown in Figure 1a, oscillations of a characteristic amplitude are sustained even in constant darkness or constant light, that is, in the absence of a periodic external drive.

Free running clocks and damped ‘hourglass’ clocks are equally good time-keepers in noiseless conditions but internal and external fluctuations reveal significant differences.

(a) Free running circadian clocks, such as the KaiABC protein clock in S. elongatus, show rhythms in both oscillating and constant light (top) or dark (bottom) conditions. (b) In contrast, damped circadian clocks, such as that in P. marinus which lacks Kai A, show rhythms only in changing light conditions and decay to a fixed state in constant conditions. (c) When subject to external noise (i.e., weather-related amplitude fluctuations in light), simulations of the free running clock show low population variance while the damped clock shows high variance. In contrast, Gillespie simulations with high internal noise due to low copy number of Kai molecule reveals that damped clocks are much more robust than free running clocks. (d) A systematic study of clock precision (i.e., mutual information between clock state and time) at fixed external noise but decreasing Kai protein copy number N reveals that free running clocks are preferred at low internal noise but damped clocks are preferable at sufficiently high internal noise.

https://doi.org/10.7554/eLife.37624.003

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

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

External fluctuations

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

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

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

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

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

Noise resistance in other biochemical clocks

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

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

A diverse range of biochemical oscillators show the trade-off between resistance to external and internal noise.

For each oscillator, the regime (green) of largest free running amplitude relative to the driving strength is most robust to external fluctuations but is most fragile to internal noise. In contrast, damped oscillations (red) are robust to internal noise and thus preferable at sufficiently high internal noise. Regimes (purple) of intermediate free running amplitude are preferred at intermediate internal noise levels. (a–g) Diverse biochemical oscillators from the literature were simulated with increasing internal noise ϵint=1/N while driven by a periodic square wave light signal with fixed strength external noise, using the external coupling and parameters specified in 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 Abou-Jaoudé, 2013; Kondepudi and Prigogine, 2014; Elowitz and Leibler, 2000; Buşe et al., 2009). Clock precision is defined as mutual information between outputs and time. The original publications identified a Hopf bifurcation parameter in these models, with free running oscillations on one side and damped oscillations on the other. Green and purple data correspond to parameter regimes with large and smaller amplitude free running oscillations relative to driving amplitude while the red data corresponds to a damped oscillator. Details in Appendix 2.

https://doi.org/10.7554/eLife.37624.004

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

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

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

Dynamical systems theory of noise resistance

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

Limit cycle clocks and point attractor clocks

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

Experiments and models of biological clocks show that external driving can be viewed as a switch between distinct day-time and night-time dynamics.

(a) Experiments on the Kai system at distinct ATP levels corresponding to day and night conditions reveal limit cycles shifted relative to each other in a phosphorylation space for Kai (reproduced from [Leypunskiy et al., 2017[). Similar behavior (Winfree, 2001) is seen in models of diverse biochemical oscillators studied in Figure 2. (b) We build a minimal model of such driven clocks as a limit cycle of radius R whose center is shifted by a distance L between day and night. In cycling conditions (see signal in (d)), an entrained clock’s state executes a trajectory that encompasses both limit cycles as shown (bottom). (c) For damped clocks (Saunders, 2002), phenomenology suggests that the day and night limit cycle dynamics are replaced by a point attractor whose position changes between day and night. The relaxation time between the day and night attractors is comparable to 12 hours, giving rise to the trajectory shown in cycling conditions. (d) The plot shows cycling conditions of light intensities that couple to (b) and (c).

https://doi.org/10.7554/eLife.37624.005

The Kai model and indeed the diverse range of biochemical oscillators in Figure 2 show such a change in the limit cycle between day and night conditions. Here, we build a minimal model of such free-running clocks using circular day and night limit cycles of radius R in a plane. The limit cycle is defined by the dynamics τrelaxr˙=r-r3/R2,θ˙=ω about its own center; but the center of the limit cycle itself moves along the x axis in Figure 3b as (ρ(t)L,0) where ρ(t)[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 τrelax. Note that R represents the amplitude of free-running oscillations while L represents the strength or amplitude of the external driving signal.

In contrast, damped clocks are phenomenologically well-described by a day-time and a night-time point attractor with slow relaxation dynamics between them (Figure 3c), modeled as r˙=-r/τrelax,θ˙=ω about an attractor point whose location varies with current light levels as (-ρ(t)L,0). Here we assume 2τrelax24 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 time-scales as shown in Figure 4a. Much like with biochemical circuits, we subject an in silico population of dynamical system clock models to different realizations of such noisy weather patterns.

External weather-related light fluctuations are filtered out by limit cycle attractors but not by point attractors.

(a) Light intensity levels fluctuate on a range of time scales due to weather (power spectrum reproduced from Gu et al., 2001). (b) A population of limit cycle clocks of identical fixed geometry, subject to different realizations of weather conditions, show non-overlapping distributions (purple blobs) at different times of the day. Point attractor clocks form larger and more overlapping distributions. (c) A single representative dark pulse of 2.4 duration causes only a ΔΦ30 min phase lag in limit cycles since the trajectory’s deviation (purple) is fundamentally bounded by the circular attractor. In contrast, ΔΦ4 hr for the point attractor since the trajectory is in free-fall towards the blue night-time attractor. (d) The geometrically computed ΔΦ2 phase shift for a dark pulse of any fixed duration and time of occurrence (see Appendix 5) drops rapidly as (R/L)-2 for large-R/L limit cycles; this theoretical prediction agrees well with the population variance gain over a day in simulations. (e) Consequently, weakly driven limit cycles (i.e., high R/L) can tell time with high precision.

https://doi.org/10.7554/eLife.37624.006

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 σclouds2 in terms of the phase lag ΔΦ 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 ΔΦ due to the particular dark pulse shown in Figure 4c is much smaller for limit cycles (ΔΦ0.5 hr for the R,L geometry shown) than for point attractors (ΔΦ4 hr) (see Appendix 5).

In fact, this contrast in ΔΦ 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 (ΔΦ)2(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,

σ2dayσ2+σclouds2dusk(σ2+σclouds2)/s2night(σ2+σclouds2)/s2dawn(σ2+σclouds2)/s4.

Solving for steady state phase variance (σ2=(σ2+σclouds2)/s4), we obtain

(1) σlimitcycle2,extΔΦ2/(s4-1).

where we have equated σclouds2 to ΔΦ2 for a typical dark pulse Here, s2 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 s2-1L/R for large-R/L limit cycles. Plugging this and ΔΦ2(L/R)2 into Equation1, we see that σ2L/R0 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 fluctuations severely affect continuous attractors but not point attractors.

(a) We model fluctuations due to finite copy number N as Langevin noise with mean zero and standard deviation ϵint, resulting in a diffusion constant ϵint21/N for the clock state. (b) The flat direction of limit cycles cannot contain diffusion, leading to large increases ϵint2Tday in population variance of clock state during each day (and night). In contrast, point attractor dynamics have constant curvature at all times, leading to a constant population variance over time. (c) The variance drops σ2σ2/s2 at dawn and dusk for limit cycles during the off-attractor dynamics between the day and night cycles. As with external noise, the variance drop is predicted by the slope dP(θ)/dθ of the circle map between the cycles. This dawn/dusk drop goes to zero for large R/L limit cycles but variance still increases during the day and night. (d) The variance for point attractors is Dτrelax, a constant determined by the curvature τrelax-1 of the harmonic potential. (e) Thus, with only internal noise present, the precision of limit cycle clocks increases with increasing separation L/R, asymptotically approaching the performance of point attractors.

https://doi.org/10.7554/eLife.37624.007

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 ϵint1/N, where N is the overall copy number, with the ratios of different species assumed fixed (see Appendix 3). Such a Langevin approach is an approximation Gillespie, 2007 to the exact Gillespie method used in Figure 1.

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

We can understand the weakness of limit cycle attractor relative to the point attractor in terms of diffusion during day/night balanced by dawn/dusk transitions. The flat direction along the limit cycle attractor cannot contain diffusion caused by the Langevin noise during the day/night; hence the phase variance increases by σ2σ2+ϵint2Tday during a day of length Tday (and similarly at night).

Dawn and dusk times reduce the phase variance σ2σ2/s2 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 s2 entirely through geometric considerations. We define the circle map ϕ=P(θ) as relating originating points θ near dusk on the day cycle to final points on the night cycle ϕ after relaxation (experimentally characterized in Leypunskiy et al., 2017). Then s-1=dP(θ)/dθ. Figure 5c shows that this slope s-1=dP(θ)/dθ, geometrically computed in the SI, agrees with the dawn/dusk variance drop in Langevin simulations and scales as s2-1L/R for large R/L.

Thus, the population phase variance changes as

σ2Dayσ2+ϵint2TdayDusk(σ2+ϵint2Tday)/s2Night(σ2+ϵint2Tday)/s2+ϵint2TdayDawn((σ2+ϵint2Tday)/s2+ϵint2Tday)/s2.

Assuming T=Tday=Tnight and solving for steady-state phase variance (σ2=((σ2+ϵint2Tday)/s2+ϵint2Tday)/s2), we obtain

(2) σcycle2,intϵint2Ts2-1

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

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

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

σpoint2,intϵint2τrelax

in the SI, which matches Langevin simulations (Figure 5d). Since τrelaxTday to have distinct clock states through the day (Figure 3), we find σcycle2,intσpoint2,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/τrelax. But limit cycle clocks experience such ‘curved’ off-attractor dynamics only at dawn and dusk, which is offset by dephasing (Mihalcescu et al., 2004; Gonze et al., 2002) during long periods of zero curvature on the limit cycle (day/night). Hence limit cycles underperform point attractors under high internal noise.

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 ϵint21/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,

Large-R/L limit cycle attractors, which correspond to large amplitude free running clocks, outperform all other oscillators in projecting out external noise but are least robust to internal noise.

(a) Point attractors and smaller R/L limit cycles (red and purple curves) show low precision (i.e., low mutual information) but do not degrade as much as large-R/L limit cycles with increasing internal noise ϵint. Thus this simple dynamical systems model of clocks reproduces and explains the trade-off seen in the complex biochemical clocks shown in Figures 1 and 2. (b,c) Speed-precision trade-off. (b) With external noise alone, the most precise clocks (i.e., large R/L limit cycles) average over longer signal history and are thus the slowest to entrain, that is, slow to transform a population with uniform phase distribution to the steady state distribution. (c) However, with internal noise alone, there is no trade-off between speed and precision; faster entraining clocks (i.e., point attractors) are more accurate since slow clocks are exposed to more internal noise.

https://doi.org/10.7554/eLife.37624.008
(3) (L/R)optimal=ϵintϵext.

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

(4) σint2σext2Q

where Q=ϵint2ϵext2 and where σint2 is the population angular variance of the clock state due to internal noise when driven by a noiseless external signal and σext2 is the population angular variance in the absence of internal noise due to amplitude fluctuations of the external signal. Note that angular variance is a better indicator than variance because we want to know how well the system can tell time.

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

However, in holding Q fixed, our trade-off does assume that the strength of the external fluctuations ϵext – that is, the fractional size of amplitude fluctuations in the external signal – is held fixed. Similarly, we hold ϵint2, 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.

Speed-precision trade-off

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

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

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 weather-related 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 trade-off 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 Abou-Jaoudé, 2013; Kondepudi and Prigogine, 2014; Elowitz and Leibler, 2000; Buşe et al., 2009). The trade-off 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 trade-off 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 non-planar non-circular limit cycles and yet reproduce our trade-off. 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 limit-cycles 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/elifesciences-publications/elife_tradeoff_clocks).

Appendix 1

Trade-off in Kai-based clocks

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

S. elongatus clock - hexamers with collective KaiA feedback

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

Appendix 1—figure 1
Explicit biochemical KaiABC model simulated using the Gillespie algorithm.

(a) The experimentally well-characterized clock in S. elongatus consists of a negative feedback-enabled self-sustained oscillator. KaiBC complexes sequester KaiA, preventing runaway KaiC molecules from going through the cycle independently. (b) The genome of P. marinus lacks kaiA. We assume a minimal model consistent with known facts (Rust et al., 2007) about this clock; KaiC phosphorylation proceeds without KaiA and hence different KaiC hexamers can proceed independently through the cycle. (c) We combine both clocks in one model with an interpolating parameter γ that selects between an S. elongatus-like KaiA-dependent pathway and an P. marinus-like KaiA-independent pathway. All reactions shown are assumed to be first order mass-action kinetics. We simulate such a system at different overall copy numbers N using the Gillespie algorithm. (d) We find limit cycles for γ>0.9. The resulting limit cycles for γ=1,0.95 violate the simplifying assumptions used in our dynamical systems (e.g., non-circular cycles of different size); and yet our results are qualitatively validated by this model (Figure 1d from the main text).

https://doi.org/10.7554/eLife.37624.011

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 KaiA-dependent phosphorylation pathway, much like in S. elongatus, that is used during the day and driven forward by ATP.

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

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

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

However, as γ0, the KaiA-independent pathway is more active and thus the system effectively has no feedback. In fact, we find that at about γ0.82, sustained oscillations disappear (for kinetic parameters used here and reported below). Hence we chose γ=1,0.95,0 as representative of two limit cycle-based free running clocks and one point-attractor based damped clock respectively. In this way, we can view the clock dynamics of S. elongatus and P. marinus can be viewed as being on either side of the Hopf bifurcation that occurs at γ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 fATP=ATP/(ATP+ADP) during the day between the fATPday and fATPnight+(fATPday-fATPnight)/3, where fATPday, fATPnight are the ATP values during a cloudless day and night respectively. We used different day and night ATP levels for different γ that ensure that the limit cycles had periods comparable to 24 hours. For γ=1, we used ATP/ADP ratios of fATPday=0.55,fATPnight=0.45. For γ=0.95, we used fATPday=0.57,fATPnight=0.17 and for γ=0, fATPday=0.8,fATPnight=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.01hr, k+=k=2m0.04932hr1, kAon=0.2466μM1hr1, kAoff=0.02466hr1, kCC=0.2466hr1, kABC=123.30hr1,m=18. We set up Kai C and Kai A in a 1:1 stoichiometric ratio, each present at a copy number N where N was varied systematically. These rates are consistent with those measured in (Qin et al., 2010; Snijder et al., 2017; Hayashi et al., 2004).

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

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 trade-off described in the paper in these models; a summary of our results is presented in the main paper.

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

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

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

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 (PN,PC). The equations used in Gonze and Goldbeter (2006) to model this are,

dMdt=vsΩ(KIΩ)n(KIΩ)n+PNn+vmΩMKmΩ+M
dPCdt=ksMvdΩPCKdΩ+PCk1PC+k2PN
(5) dPNdt=k1PCk2PN

where νs is an mRNA transcription rate that is modulated by external signals (Leloup et al., 1999; Gonze and Goldbeter, 2006), and Ω 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): KI = 1 nM, n = 4, vm = 0.505 nM/h, Km = 0.5 nM, ks= 0.5 1/h, vd = 1.4 nM/h, Kd= 0.13 nM, k1 = 0.5 nM/h, k2 = 0.6 nM/h, and assume the volume Ω 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):

(6) dXdt=μ(x,t)+Σ(x,t)η(0,1)

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

(7) Σ(x,t)=[AB0000000ksMk1PCk2PN0000000k1PCk2PN]

where A=vsΩ(KIΩ)n(KIΩ)n+PNn and B=vmΩMKmΩ+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 dXdt=μ(x,t)+1ΩΣ(x,t)η(0,1).

As in Leloup et al. (1999) and Gonze and Goldbeter (2006), we take ν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 νs = 0.57 nM/h. Hence, in Figure 2 from the main text, we used νsDay = 0.55 nM/h,νsNight=0.05 nM/h for the point attractor (red). For the two limit cycles, we used νsDay = 0.9 nM/h,νsNight=0.6 nM/h (green), and νsDay = 0.705 nM/h,νsNight=0.695 nM/h (purple). The driving period is 18 hr, similar to the driving period of the system at v1 = 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 ML(t) and PL(t) respectively. Denoting the mRNA and protein levels of TOC1 by MT(t) and PT(t), Schmal et al. (2014) present a reduced version of the model in Locke et al. (2005) as:

dMLdt=L(t)+v1PTag1a+PTam1MLk1+ML
dPLdt=p1MLm2PLk2+PL
dMTdt=v2g2bg2b+PLbm3MTk3+MT
(8) dPNdt=p2MTm4PTk4+PT

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, g1 = 0.5 nM, g2 = 0.1 nM, m1 = 0.4 nM/h, m2 = 0.6 nM/h, m3 = 0.6 nM/h, m4 = 0.3 nM/h, k1 = 1 nM, k2 = 0.5 nM, k3 = 1 nM, k3 = 1 nM, p1 = 0.5 1/h, p2 = 0.3 1/h, v2 = 0.6 nM/h. With other parameters fixed, the system undergoes Hopf bifurcation at v1 = 0.194 nM/h We use v1 = 0.26 nM/h for limit cycles and v1 = 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 LDay = 0.05 nM/h and LNight = 0 (green data) and we set LDay = 0.01 nM/h and LNight = 0 (purple data). For point attractor, we set LDay = 0.2 nM/h and LNight = 0 (red data). The period of the driving signal is 24 hr, which is around the natural period of the system when v1 = 0.26 nM/h and L = 0.

B.3 Mammalian Per-Cry 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 DXi for species Xi. 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 vsP (a transcriptional rate) to be the light input function. We use vsPDay = 1.09 nM/h and vsPNight = 1.07 nM/h for the purple limit cycle data in Figure 2, vsPDay = 1.15 nM/h and vsPNight = 1.07 nM/h for the green limit cycle data. For the point attactor data (red), we set vsPDay = 1.5 nM/h and vsPNight = 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 vsP = 1.07 nM/h.

B.4 cdc2-cyclin 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 self-sustained nature of cell cycles.

dCdt=viΩkdCvdXΩCKdΩ+C
dMdt=v1CKCΩ+C1MK1+(1M)V2MK2+M
(9) dXdt=v3M(1X)K3+(1X)V4XK4+X

where Ω 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 Ω 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): Ki=0.1 (i=1–4), VM1=0.5min-1, V2=0.167min-1, VM3=0.2min-1, V4=0.1min-1, vd=0.1μMmin-1, KC=0.3μM, Kd=0.02μM, kd=3.33×10-3min-1.

Goldbeter (Goldbeter, 1991) suggested that vi is modulated by external signals. So, we use viDay = 0.0106 μMmin-1 and viNight = 0.0105 μMmin-1 for small L/R limit cycle, viDay = 0.0111 μMmin-1 and viNight = 0.0105 μMmin-1 for large L/R limit cycle, and viDay = 0.009 μMmin-1 and viNight = 0 for point attractor. The bifurcation from point attractor to limit cycle happen around vi = 0.01 μMmin-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 widely-studied version of such a Goodwin oscillator (Gonze and Abou-Jaoudé, 2013; Woller et al., 2014),

dXdt=α(t)1+ZnX
dYdt=XY
(10) dZdt=YZ

When n=9, the limit cycles disappear at a Hopf bifurcation found at α7. As is commonly done (Woller et al., 2014), we couple the external signal to the bifurcation parameter α(t). We use αDay=120,αNight=80 for the green limit cycle in Figure 2c of the main paper, αDay=108,αNight=92 for purple limit cycle data, and αDay=2.5,α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 α=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 non-linearity spread equally amongst all three reactions, lowering the requisite Hill coefficient of any one reaction to a biochemically realistic n=3. Repressilator circuits are not usually driven by an external signal, except in a few theoretical analyses (e.g., Russo et al., 2010; Schultz et al., 2013). We use the simplest version of these, with the driving signal modulating the transcription rate of only one of the proteins

dXdt=α(1+s(t))1+YnX
dYdt=α1+ZnY
(11) dZdt=α1+XnZ

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

We use sDay=0.7/5.2,sNight=-0.7/5.2 for the green limit cycle in Figure 2g in the main paper, sDay=4.8/5.2,sNight=-1.7/5.2 for the purple limit cycle data, and sDay=0.5/1.9,sNight=-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 α=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:

dXdt=1(1+b(t))X+X2Y
(12) dYdt=b(t)XX2Y

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 bDay=2.25 and bNight=2.2 for the purple limit cycle data in Figure 2f from main text. We use bDay=2.8 and bNight=2.2 for the green limit cycle data. Lastly, we use bDay=1.8 and bNight=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:

(13) r˙=α(rr3R2)
(14) θ˙=ω

For α>0, the above equation describes a circular limit cycle of radius R and frequency ω. This equation undergoes a Hopf bifurcation at α=0, where the limit cycle shrinks to zero and resulting in point attractor for α<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:

(15) dx=(α(1-x2+y2R2)x-ωy)dt+2DdW
(16) dy=(α(1-x2+y2R2)y+ωx)dt+2DdW

where DR2ϵint2 is the diffusion constant, dW is a Wiener process, and here we assume that the internal noise is a homogeneous white-noise in the 2-dimensional space. (Similar assumptions are made in reference Potoyan and Wolynes, 2014)

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

In Equation 14, τrelax1|α| 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π/ω, 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 τrelax1|α|2π/ω comparable to the period of the day-night cycle itself. As explained in the main paper, if relaxation were much faster, the clock state would decay to a fixed point before the end of the day (or night) and thus not show distinct states at distinct times of the day-night cycle.

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 α=5 for limit cycle system and α=-5 for point attractor system where ω=2π 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 r3R2 is small. The center of the cycle and point attractors during the day are assumed to be at (-L,0) and at (0,0) at night; or more generally at (-Ls(t),0) were s(t) the light signal (assumed to be of amplitude 1).

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

To simulate internal noise, we add Langevin noise to Equations 14 as described in Appendix 6 on Langevin noise. As explained in Appendix 6, our measure of internal noise ϵint2 is a measure of phase diffusion, independent of limit cycle size. In other words, ϵint2 is a measure of the population phase variance (i.e., variance in θ) 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 ϵint2 above as the phase diffusion in undriven conditions for limit cycles of any size, ensures that changing L and R are equivalent for limit cycles.

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

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

For limit cycles, we initialize the population of 104 particles by uniformly distributing them along the perimeter of the night cycle. In the point attractor system, we initialize a population of 105 at the night-time point attractor.

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

Optimal Dynamical system and trade-off

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

(17) σint2ϵint2Ts2-1
(18) σext2ΔΦ2s4-1

For a small L/R, we had found that ΔΦϵextL/R where ϵext2 is a measure of the variance of light during the day. Further, we showed that s2-1L/R in this limit. Hence, in the small L/R (‘Stuart-Landau’) regime, the above equations reduce to,

(19) σint2ϵint2R/L
(20) σext2ϵext2L/R

The population variance when both noises are present is approximately given by max(ϵint2R/L,ϵext2L/R). This variance is minimized when the two terms are equal, giving

(LR)optimalϵintϵext,

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

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

(21) σext2σint2Qϵext2ϵint2

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

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

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

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

Appendix 4

Supplementary Methods

Modeling external noise (weather fluctuations) We generate a square wave of period 24 hours to model the day-night cycle of light on Earth with the day length of 12 hr. However, such a square wave is modulated by weather fluctuations, for example, periods of reduced intensity due to passing clouds during the daytime. We model such fluctuating intensity as follows. We assume each weather condition lasts a random interval of time drawn from an exponential distribution of mean 2.4 hrs (1/10 of a day). During a given weather condition, we set the intensity of light to a random value, drawn uniformly from [1-noiseext,1] where noiseext 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 noiseext to 1. However, in our simulations for eight different models of biological clocks, noiseext ranges from 0.5 to 1 because when noiseext 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 ρ[0,1], we switch the dynamics to an alternative limit cycle (or point attractor) at a fractional distance ρ 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 ρ[0,1], we follow dynamics due to a limit cycle located at (-ρL,0). We follow the same rules for the point attractor.

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

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

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

Measures of clock time-telling 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 c and the time t,

(22) MI(C;T)=cC,tTp(c,t)log2(p(c,t)p(c)p(t))

for all c in the set of available positions C and all t in the available time bins T. (In the dynamical systems model, c represents the position in the 2d r,t plane. For the explicit KaiABC biomolecular model, 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(c|t) of clock states at a given (objective) time t[0,24] hrs of the last day in the simulation. The probability function p(c) is calculated by accumulating the distribution of p(c|t) over time t[0,24] hrs of the last day. The position c and time t are binned into different bins depending on their values. We start the minimum and maximum values of the bins to the minimum and maximum values of the variables. The bin size in the time dimension is 0.48 hr or 28.8 min, while The bin size in the x and y dimensions are both 0.01.

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

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

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

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

To see this, we write mutual information between clock state c and time t as,

(23) MI(C;T)=H(T)H(T|C).

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

(24)H(T|C)=p(c)dcH(T|c)(25)=p(c)dc[dtp(t|c)logp(t|c)]

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

Appendix 4—figure 1
Mutual information MI(c,t) between clock state c and time t is only affected by the variance of the clock state distribution p(c|t) at a given time t along the direction of motion and not orthogonal to it.

In this toy example, we assume the distribution p(c|t) to be supported on a rectangle of size ax and ay in a 2d clock state space. The clock state moves at a speed u in the x-direction. Time telling quality is affected by how much the population at different times overlap with each other. Consequently, clocks with large ax and small ay (bottom) have lower mutual information MI(c,t) relative to clocks with small ax and large ay (top). Consequently, we use the population variance along the direction of motion as an instantaneous measure of time-telling ability in the paper.

https://doi.org/10.7554/eLife.37624.015
p(t|c)={0for |cxut|>axu2axfor |cxut|ax

So,

H(T|C)=p(c)dct=(cxax)/u(ax+cx)/udtu2axlog(u2ax)=log(2axu)

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

Cramer-Rao bounds

Cramer-Rao (CR) bounds quantify the total available information about phase in a given length of history of the signal. Any estimator working with that length of history must necessarily have higher variance (i.e., higher error) than the Cramer- lower bound corresponding to that length of history. In the limit of infinitely long histories, the CR bound in this context corresponds to zero error; with any finite binning in time, the upper bound on MI is simply set by the number of bins in time. In our case, this bound is given by log250=5.64 bits. As shown in the main paper, as L/R0, 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/R0 (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 σ2σ2+σclouds2 scales as (L/R)2, in perfect agree with stochastic weather simulations.

To compute the scaling relationship of σclouds2, we compute the phase shift ΔΦ caused by a single dark pulse with width τ on the limit cycles with angular speed ω (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 θ on the day cycle. A dark pulse of length τ administered just then will change the dynamics to that of the night cycle. This clock has state ϕ=P(θ) with respect to the night cycle and will evolve for a time τ according to the night cycle dynamics, reaching a new state ϕ+ωτ, at a radial position determined by R,L. At the end of the dark pulse, we use the night-day circle map, θ=Q(ϕ), 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:

Appendix 5—figure 1
The population variance of clock states is reduced by dusk and can be computed geometrically.

(a) A population of clocks near state θ on the day cycle is mapped to the neighborhood of state ϕ on the night cycle by the dusk transition. We define ϕ=P(θ) to be the map relating the clock state θ on the day cycle just before dusk to its eventual position ϕ on the night cycle after dusk (assumed greater than the relaxation time). (b) This map can be analytically computed for circles of size R with centers separated by length L. (c) For a given R/L = 2 , we obtain P(θ) shown here. Since θ=π/2 corresponds to the dusk time of the entrained trajectory, the slope s1=dP/dθ at θ=π/2 determines the change in population variance of clock states at dusk. (d,e) The variance drop s2 at dusk, defined as σ2σ2/s2 at dusk, seen in both the external (averaging over weather) and internal noise (averaging over Langevin noise) simulations agree well with the geometrically computed s(R/L), especially at large R/L. We find that s21L/R for large-R/L limit cycles.

https://doi.org/10.7554/eLife.37624.017
Appendix 5—figure 2
Increase in population variance due to random weather conditions can be estimated from the phase shifts ΔΦ due to dark pulses (i.e., the Phase Response Curve).

(a) A single dark pulse administered during the day shifts the phase of a clock (purple) relative to a clock that experiences no such dark pulse (black). (b) We can compute the phase shift ΔΦ due to such a dark pulse geometrically by computing the deviation in trajectory. Assuming a dark pulse of length τ, the clock evolves for a time τ according to the night cycle dynamics. At the end of such a pulse, we switch back to the day limit cycle and compute the resulting phase shift ΔΦ. (c) The resulting phase shift ΔΦ due to a pulse of length τ=2.4 hrs, depends on the time θ when it is administered but is generally smaller for larger R/L. (d) We find that ΔΦ2 for a specific τ=2.4 hrs dark pulse administered at the same time (8 AM) falls as (L/R)2 for large-R/L limit cycles. This trend matches the variance gain σclouds2 seen in stochastic simulations that average over random weather conditions (pulses of different length, intensity and time of application). The broken brown curve shows a theoretical prediction for such an average ΔΦ2, obtained by sampling the curve shown in (c) at different points of application and differing intensity. Despite the presence of a variance-reducing zero around mid-day in (c), σclouds2 drops as (L/R)2, much as ΔΦ2 for any particular pulse. (Brown theory curves translated together using one fitting parameter).

https://doi.org/10.7554/eLife.37624.018
(26) ϕ=P(θ)=arctan(L+RsinθRcosθ)

and

(27) θ=Q(ϕ)=arctan(L+Rsin(ϕ+ωτ)Rcos(ϕ+ωτ)).

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 ϕ=P(θ), ‘2’ corresponds to the rotation on the night cycle ϕϕ+ωτ, and ‘3’ corresponds to the transition back to the day cycle θ=Q(ϕ+ωτ). Combining this three processes, we write θ as θ(θ,τ,L/R) and expand it in the limit that L/R0 to obtain that

(28) ΔΦ=LR(cos(θ+ωτ)cos(θ))+𝒪(LR)2

where ΔΦ=θ(θ+ωτ) because θ+ωτ is the phase of the clock if it did not experience the dark pulse.

This expression ΔΦ 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 |ΔΦ|2(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 ΔΦ(θ) 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 τ, there exists a time of administration such that ΔΦ=0! In fact, such a dark pulse has an entraining effect, reducing the population variance. We leave experimental and theoretical investigation of the counter-intuitive effects of such specially timed dark pulses to future work.

Here, we show that even if we include such dark pulses with an entraining effect, the variance gained at the end of the day is still proportional to (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 (θ=π/2), around noon (θ=ωτ/2 with small ωτ), and in the evening (θ=π/2). Starting the day with variance σ02, by the end of the day the variance becomes

(29)σ2=σ02+(ΔΦ)θ=π/22(1+(dΔΦdθ)θ=ωτ/2)2+(ΔΦ)θ=π/22(30)σ02+(LRsinωτ)2(1+2LRsin(ωτ2))2+(LRsinωτ)2
(31) σ2σ02+2(LRsinωτ)2+𝒪(LR)2.

Thus, the variance gained due to fluctuation, σ2σ02=σclouds2, 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 ΔΦ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 σ2σ2/s2 at dusk where s21L/R as L/R0. Data from Langevin simulations was presented. Here we will derive this result using a simple geometric argument about circle maps.

We define ϕ=PT(θ) 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 θ on the day cycle. See Appendix 5—figure 1a,b. In principle, with complex relaxation dynamics between the limit cycles, PT(θ) 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 ϕ=PT(θ); that is, we can write PT(θ)=P(θ)+ωT where ω is the angular frequency of the clock. In what follows, we will be interested in the derivative of θPT(θ); hence we will work with P(θ) instead of PT(θ).

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

(32) ΔϕΔθdP(θ)dθ|θ=θ0

By a similar argument, if the phase variance of clock states across a population is σ2 before dusk, it will be reduced by,

(33) σ2duskσ2(dP(θ)dθ|θ=θ0)2

This expression is valid in the regime where the population variance σ2 is small enough to linearize the circle map P(θ). 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(θ) from geometry. In Appendix 5—figure 1, we draw a diagram of the transition from a particle on the day cycle at the phase θ to the night cycle at the phase ϕ. By trigonometry, we write

(34) ϕ=P(θ)=arctan(L+RsinθRcosθ),

and derive

(35)s21=(dP(θ)dθ)21(36)=L(2L3+7LR23LR2cos(2θ)+4R(2L2+R2)sinθ)2R2(R+Lsinθ)2(37)=2sin(θ)LR+𝒪(LR)2,

where θ corresponds to the angle on the day cycle at dusk, which is at π/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 σ2σ2/12 at the transition).

Remarkably, our geometric derivation of s21 matches the variance drop σ2σ2/s2 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 s2 and thus with the external noise simulations as well. In both cases, the simulations and geometric theory show that s21L/R as L/R0.

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

(38) A+B>C+D

with the forward rate constant k+, during each infinitesimal time δt, the probability of the occurrence of this reaction follows a Poisson distribution whose mean and variance both equal to R+δt=k+NANBδt. Integration over a larger time step, the Poisson distribution can be approximated into a Gaussian form, resulting in Langevin dynamics,

(39) dNA=k+NANBdt+R+dW

where W is a standard Wiener process of mean 0 and autocorrelation function W(t1)W(t2)=δ(t1t2).

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

(40) C+D>A+B

with the rate constant k-, then the Langevin equation for species A becomes,

(41) dNA=k+NANBdt+kNCNDdt+R+dW1+RdW2

where R+=k+NANB and R=kNCND respectively denote the number rates of the forward and the backward reaction; dW1 and dW2 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 time-dependent 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 bio-oscillators, presented later, do not make this simplifying assumption of homogeneous diffusion.) This allows us to write a 2-dimension phenomenological stochastic differential equation

(42) dz=f(z,t)dt+2DdW

where the f(z,t) denotes the deterministic dynamics driven by day-night cycles and the diffusion constant D is assumed to be inversely proportional to the total number of Kai-C molecules within the cell. For limit cycles of radius R, we set DR2ϵint2. Then, ϵint2 is the diffusion constant for the phase of the oscillator. We hold ϵint2 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 steady-state population variance for a point attractor mechanism. The population variance is together determined by the diffusive term 2DdW, and the noise reduction effect from the restoring force of the point attractor’s harmonic well. During each infinitesimal time δt, the internal noise increase the variance by

(43) σ2(t+δt)=σ2(t)+2Dδt.

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 1-d overdamped dynamics of a particle within a harmonic energy well of V(r)=kr2. The solution to the equation of motion is r(t)=r0e2kt, with initial position r(0)=r0. Consider an ensemble of points with a mean initial position μ0 and a initial variance of σ02, one can solve the dynamics of the mean as

(44) μ(t)=μ0e2kt

and the dynamics of the variance as

(45) σ2(t)=σ02e4kt

Thus, per infinitesimal time δt, the geometric focusing effect of the energy well of the point attractor reduces the population variance by

(46) σ2(t+δt)=σ2(t)/g

where g=e4kδ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

(47) σst2=σst2+2Dδtg=σst2+2Dδte4kδt

and by taking the limit of δt goes to 0, we have σst2=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

    1. Gonze D
    2. Goldbeter A
    (2006) Circadian rhythms and molecular noise
    Chaos: An Interdisciplinary Journal of Nonlinear Science 16:026110.
    https://doi.org/10.1063/1.2211767
  1. Book
    1. Kondepudi D
    2. Prigogine I
    (2014)
    Modern Thermodynamics: From Heat Engines to Dissipative Structures
    John Wiley & Sons.
    1. Laughlin S
    (1981)
    A simple coding procedure enhances a neuron’s information capacity
    Zeitschrift Fur Naturforschung. Section C, Biosciences 36:910–912.
  2. Book
    1. Saunders DS
    (2002)
    Insect Clocks (Third Edition)
    Elsevier.

Article and author information

Author details

  1. Weerapat Pittayakanchit

    1. Department of Physics, University of Chicago, Chicago, United States
    2. The James Franck Institute, University of Chicago, Chicago, United States
    Contribution
    Conceptualization, Software, Formal analysis, Visualization, Methodology, Writing—original draft, Writing—review and editing
    Contributed equally with
    Zhiyue Lu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7940-3184
  2. Zhiyue Lu

    1. Department of Physics, University of Chicago, Chicago, United States
    2. The James Franck Institute, University of Chicago, Chicago, United States
    Contribution
    Conceptualization, Software, Formal analysis, Visualization, Methodology, Writing—original draft, Writing—review and editing
    Contributed equally with
    Weerapat Pittayakanchit
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0216-4346
  3. Justin Chew

    Medical Scientist Training Program, Pritzker School of Medicine, University of Chicago, Chicago, United States
    Contribution
    Conceptualization, Software, Validation, Investigation, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4749-547X
  4. Michael J Rust

    1. Department of Physics, University of Chicago, Chicago, United States
    2. The James Franck Institute, University of Chicago, Chicago, United States
    3. Department of Molecular Genetics and Cell Biology, University of Chicago, Chicago, United States
    Contribution
    Conceptualization, Investigation, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7207-4020
  5. Arvind Murugan

    1. Department of Physics, University of Chicago, Chicago, United States
    2. The James Franck Institute, University of Chicago, Chicago, United States
    Contribution
    Conceptualization, Formal analysis, Supervision, Funding acquisition, Visualization, Methodology, Writing—original draft, Project administration, Writing—review and editing
    For correspondence
    amurugan@uchicago.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5464-917X

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

  1. Received: April 17, 2018
  2. Accepted: June 23, 2018
  3. Accepted Manuscript published: July 10, 2018 (version 1)
  4. 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,551
    Page views
  • 589
    Downloads
  • 27
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

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)

  1. Weerapat Pittayakanchit
  2. Zhiyue Lu
  3. Justin Chew
  4. Michael J Rust
  5. Arvind Murugan
(2018)
Biophysical clocks face a trade-off between internal and external noise resistance
eLife 7:e37624.
https://doi.org/10.7554/eLife.37624

Share this article

https://doi.org/10.7554/eLife.37624

Further reading

    1. Cancer Biology
    2. Computational and Systems Biology
    Bingrui Li, Fernanda G Kugeratski, Raghu Kalluri
    Research Article

    Non-invasive early cancer diagnosis remains challenging due to the low sensitivity and specificity of current diagnostic approaches. Exosomes are membrane-bound nanovesicles secreted by all cells that contain DNA, RNA, and proteins that are representative of the parent cells. This property, along with the abundance of exosomes in biological fluids makes them compelling candidates as biomarkers. However, a rapid and flexible exosome-based diagnostic method to distinguish human cancers across cancer types in diverse biological fluids is yet to be defined. Here, we describe a novel machine learning-based computational method to distinguish cancers using a panel of proteins associated with exosomes. Employing datasets of exosome proteins from human cell lines, tissue, plasma, serum, and urine samples from a variety of cancers, we identify Clathrin Heavy Chain (CLTC), Ezrin, (EZR), Talin-1 (TLN1), Adenylyl cyclase-associated protein 1 (CAP1), and Moesin (MSN) as highly abundant universal biomarkers for exosomes and define three panels of pan-cancer exosome proteins that distinguish cancer exosomes from other exosomes and aid in classifying cancer subtypes employing random forest models. All the models using proteins from plasma, serum, or urine-derived exosomes yield AUROC scores higher than 0.91 and demonstrate superior performance compared to Support Vector Machine, K Nearest Neighbor Classifier and Gaussian Naive Bayes. This study provides a reliable protein biomarker signature associated with cancer exosomes with scalable machine learning capability for a sensitive and specific non-invasive method of cancer diagnosis.

    1. Computational and Systems Biology
    2. Immunology and Inflammation
    Alain Pulfer, Diego Ulisse Pizzagalli ... Santiago Fernandez Gonzalez
    Tools and Resources

    Intravital microscopy has revolutionized live-cell imaging by allowing the study of spatial–temporal cell dynamics in living animals. However, the complexity of the data generated by this technology has limited the development of effective computational tools to identify and quantify cell processes. Amongst them, apoptosis is a crucial form of regulated cell death involved in tissue homeostasis and host defense. Live-cell imaging enabled the study of apoptosis at the cellular level, enhancing our understanding of its spatial–temporal regulation. However, at present, no computational method can deliver robust detection of apoptosis in microscopy timelapses. To overcome this limitation, we developed ADeS, a deep learning-based apoptosis detection system that employs the principle of activity recognition. We trained ADeS on extensive datasets containing more than 10,000 apoptotic instances collected both in vitro and in vivo, achieving a classification accuracy above 98% and outperforming state-of-the-art solutions. ADeS is the first method capable of detecting the location and duration of multiple apoptotic events in full microscopy timelapses, surpassing human performance in the same task. We demonstrated the effectiveness and robustness of ADeS across various imaging modalities, cell types, and staining techniques. Finally, we employed ADeS to quantify cell survival in vitro and tissue damage in mice, demonstrating its potential application in toxicity assays, treatment evaluation, and inflammatory dynamics. Our findings suggest that ADeS is a valuable tool for the accurate detection and quantification of apoptosis in live-cell imaging and, in particular, intravital microscopy data, providing insights into the complex spatial–temporal regulation of this process.