Circadian clocks optimally adapt to sunlight for reliable synchronization

Circadian oscillation provides selection advantages through synchronization to the daylight cycle. However, a reliable clock must be designed through two conflicting properties: entrainability to synchronize internal time with periodic stimuli such as sunlight, and regularity to oscillate with a precise period. These two aspects do not easily coexist because better entrainability favors higher sensitivity, which may sacrifice the regularity. To investigate conditions for satisfying the two properties, we analytically calculated the optimal phase-response curve with a variational method. Our result indicates an existence of a dead zone, i.e., a time period during which input stimuli neither advance nor delay the clock. A dead zone appears only when input stimuli obey the time course of actual solar radiation but a simple sine curve cannot yield a dead zone. Our calculation demonstrates that every circadian clock with a dead zone is optimally adapted to the daylight cycle.


Introduction
Circadian oscillators are prevalent in organisms from bacteria to humans and serve to synchronize bodies with the environmental 24 h cycle [1,2]. Although the molecular implementation of oscillation is species-specific [3][4][5][6], every circadian clocks satisfies two requirements to achieve reliable synchronization to the environment: entrainability to synchronize internal time with periodic stimuli and regularity to oscillate with a precise period. Circadian clocks are acquired through evolution independently in bacteria, fungi, plants and animals [7]. Nonetheless, entrainability and regularity constitute major characteristics conserved in all circadian clocks [6], which strongly suggest that these two properties are essential for survival. A main source of interference with regularity is discreteness of molecular species, i.e. molecular noise [8][9][10][11][12][13]. Many studies have analysed the resistance mechanisms of circadian oscillators against noise [14][15][16][17]. Regarding entrainability, circadian clocks synchronize their internal time with the environmental cycle via sunlight, and its effect depends on the wavelength or fluence, as well as on the phase of the stimulation. However, entrainability and regularity are conflicting factors, because circadian clocks with better entrainability are sensitive not only to the periodic light stimuli, but also to the molecular noise which interferes with regularity.
Because both regularity and entrainability are important adaptive values, we expect actual circadian oscillators to optimally satisfy these two factors (figure 1). Here, we investigate the optimal phase-response curve (PRC), which is both entrainable and regular, in the phase oscillator model [18] by using the Euler-Lagrange variational method. Our main finding is the inherent existence of a dead zone in the PRC: optimality is achieved only when the PRCs have a time period during which light stimuli neither advance nor delay the clock (figure 2a). In other words, a PRC with a dead zone (figure 2a) is better adapted than those without a dead zone (figure 2b). This result is intriguing, because a dead zone, with which oscillators tend to be unaffected by stimuli (i.e. lower entrainability), achieves better entrainability. We also tested this with two types of input stimuli: a solar radiation-type input that simulated the time course of solar radiation intensity (cf. equation (2.24) and figure 4a) and a simple sinusoidal input (sine curve). Surprisingly, the dead zone in the optimal PRC emerges only for the solar radiation-type input, not for the sinusoidal input. Many experimental studies reported the existence of a dead zone in various species (figure 2c,d show experimentally observed PRCs of fruitfly [19] and mouse [20], respectively). Our results indicate that circadian oscillators in various species have adapted to solar radiation for reliable synchronization.

Phase oscillator model
Circadian oscillators basically comprise interaction between mRNAs and proteins, whose dynamics can be modelled by differential equations. A circadian oscillator of N molecular species can be represented by dx i dt ¼ F i ðxÞ ði ¼ 1; 2; . . . ; NÞ; ð2:1Þ where the N-dimensional vector x ¼ (x 1 , x 2 , . . . ,x N ) denotes the concentration of molecular species (mRNAs or proteins). The effect of noise on genetic oscillators has been a subject of considerable interest, and noise-resistant mechanisms have been extensively studied [14][15][16][17][21][22][23]. In general, the dynamics of the ith molecular concentration in a circadian oscillator subjected to molecular noise is described by the following Langevin equation (Stratonovich interpretation): where Q i (x) is an arbitrary function representing the multiplicative terms of the noise, j i (t) is white Gaussian noise with the correlation kj i ðtÞj j ðt 0 Þl ¼ 2d ij dðt À t 0 Þ (a bracket k Á l denotes expectation), and r is a model parameter.
Circadian oscillators synchronize to environmental cycles by responding to a periodic input signal (light stimuli). We let r in equation (2.2) be stimulated by the input signal: for example, r can be the degradation rate (for the sake of simplicity, we consider that the input signal affects only one parameter). We use equation (2.2) for calculating regularity and entrainability of circadian oscillators.

Definition of regularity
Because the circadian oscillator of equation (2.2) is subjected to noise, its period varies cycle to cycle. We use the term regularity for the period variance of the oscillation (higher regularity corresponds to smaller period variance). Let us first consider the case without input signals (i.e. r is constant). As equation ( where V ¼ 2p/T is the angular frequency of the oscillation (T is a period of the oscillation). The phase f in equation (2.3) is defined only on a closed orbit of the unperturbed limitcycle oscillation. However, we can expand the definition into the entire x ¼ (x 1 , x 2 , . . . ,x N ) space, where the equiphase surface is referred to as the isochron I(f) (figure 3a). By using standard stochastic phase reduction [18], equation (2.2) can be transformed into the following Langevin equation with respect to the phase variable f (Stratonovich interpretation): iPRC U i (f) quantifies the extent of phase advance or delay when perturbed along an x i coordinate direction at phase f. The N-dimensional vector x LC (f) denotes a point on the limit-cycle trajectory at phase f, where LC stands for limit cycle. The value of iPRC U i (f) is calculated as a solution of an adjoint equation [24] or as the set of eigenvectors of a monodromy matrix in the Floquet theory [18] for arbitrary oscillators.
Let P(f; t) be the probability density function of f at time t. From equation (2.4), the Fokker-Planck equation (FPE) [25] of P(f; t) is given by With sufficiently weak noise, P(w; t) is a slowly fluctuating function of t. In such cases, F(w þ Vt) and G(w þ Vt) fluctuate much faster than P(w; t), thus these two terms can be averaged for one period while keeping P(w; t) constant ( phase averaging). In other words, we separate time scales between F(w þ Vt), G(w þ Vt) and P(w; t). By phase averaging, F(w þ Vt) vanishes because of the periodicity (use regularity optimal clock implementation entrainability higher higher feasible region infeasible region Figure 1. Illustrative relation between two trade-off properties: entrainability and regularity. There is an infeasible region with respect to entrainability and regularity (coloured area), inside which no clocks can be implemented. Actual circadian clocks are considered to optimally satisfy them and such optimal clocks lie on the edge between feasible and infeasible regions (thickdashed line).
In equation (2.4), the average period corresponds to the mean first passage time with f starting from 0 to 2p, and the period variance is the variance of the first passage time. Because direct calculation of the period variance is difficult, we approximate the period variance V T with the phase variance V f , after Kori et al. [26]. As the phase f crosses a threshold f ¼ 2p with gradient 2p/T without noise, there is a scaling relation T=ð2pÞ for sufficiently weak noise [26] (figure 3b). Consequently, the variance of the period is approximated by

Definition of entrainability
The entrainment property is an important characteristic of limit-cycle oscillators and attracts attention in systems biology [27][28][29][30][31][32]. For instance, a period mismatch in coupled oscillators is known to enhance entrainability in genetic oscillators [31]. Light stimuli affect the rate constants, i.e. the parameter r in equation (2.2) is perturbed as r þ dr by the input signal. Phase dynamics of equation (2.2) can be viewed as representing that of a tilted periodic potential (i.e. ratchet) subjected to noise. Because a synchronizable condition corresponds to the existence of stable points in the ratchet-like potential, the entrainability can be discussed without considering the noise. Consequently, in contrast to the calculation of regularity, in the evaluation of the entrainability, we consider a case without molecular noise (i.e. Q i (x) ¼ 0 in equation (2.2)). Let p(vt) be an input signal with angular frequency v. Considering a weak periodic input signal dr ¼ xp(vt), where x is the signal strength (x ! 0), and applying the phase reduction approach to equation (2.2), the time evolution of the phase variable f is given by where Z(f ) is the PRC with respect to the parameter r and corresponds to experimentally observed PRCs. In order to distinguish Z(f) from iPRC U i (f ), we will refer to Z(f ) as the parametric PRC (pPRC) [33]. Note that the definition of measured PRCs is different from pPRCs Z(f ) in a rigorous definition; the experimentally measured PRCs quantify the phase shift Df caused by light stimuli, whereas pPRCs Z(f) are normalized by the strength of perturbation, i.e. ZðfÞ ¼ @f=@r ≃Df=Dr. Therefore, the ranges of the measured  rsif.royalsocietypublishing.org J. R. Soc. Interface 11: 20131018 pPRCs have limitation 2p Df , p, whereas pPRCs Z(f ) do not. The phase reduction can yield reliable results only when the perturbed trajectory is close to the unperturbed limit-cycle trajectory (i.e. x is sufficiently small).
We next evaluate the extent of synchronization to the periodic input signal. By introducing another slow variable c ¼ f-vt in equation (2.13), we can again apply the phase-averaging procedure, which yields The oscillator of interest can synchronize to input signals when there is a stable solution of c in _ c ¼ 0 (equation (2.15)). The stable solution is an intersection point of Q(c) and -DV with dQ/dc , 0 (an empty circle in figure 3c). Then, a condition for the existence of a stable solution is ð2:17Þ where c M ¼ argmax c QðcÞ and c m ¼ argmin c QðcÞ.
We define entrainability, the extent of synchronization to the periodic input signal, by the width of the Arnold tongue, which is a domain with respect to x (signal strength) and v (signal angular frequency). The shaded region in figure 3d represents the Arnold tongue; with parameters x and v inside the Arnold tongue, the oscillator can synchronize to a periodic input signal. Because equation (2.17) constitutes a linear approximation of the Arnold tongue for sufficiently small x, the width of the Arnold tongue is given by x(Q(c M ) -Q(c m )) under the linear approximation. Thus, we define the entrainability E, or the extent of synchronization, as Intuitively, a circadian oscillator with better entrainability (i.e. larger E) can synchronize to an input signal that has a period further from that of the oscillator. The calculation above is standard in the phase reduction approach, and further details are available in Kuramoto [18].

Variational method
We use the variational method to calculate the optimal PRCs which maximize the entrainability E subject to constant variance V T ¼ s 2 T (the optimal solutions correspond to the edge in figure 1, which is described by the thick-dashed line). The constrained optimization of U i (f ) can be intuitively interpreted as maximization of weighted area (equation (2.18)), where the input being the weight, with constant area under the squared magnitude (equation (2.12)). In simple terms, the optimality is reached when the magnitude of the PRC is small during intervals when the input magnitude is small (and vice versa). In the context of neuronal oscillators, a study [34] has used the variational method to calculate the optimal PRCs for stochastic synchrony (noise-induced synchronization [35,36]).
The variational equation to be optimized is

Input signal of solar radiation model
Optimal PRCs depend on input signals, as seen in equations (2.20) and (2.21). The most common synchronizer in circadian oscillators is sunlight, for which the strength is determined by 24 h-periodic solar irradiance. The solar irradiance is calculated by I ¼ I 0 cos q and I ¼ 0 when the sun is above the horizon (0 q , p) and below the horizon ( p q , 2p), respectively, where q is the zenith angle and I 0 is the maximum irradiance [38]. It can be approximated by where ramp(x) is the ramp function defined by ramp(x) ¼ x for x ! 0 and ramp(x) ¼ 0 for x , 0. We call equation (2.24) the solar radiation input, whose plot is shown in figure 4a (the shaded region represents night). In order to show the validity of the solar radiation modelling, we compare equation (

Optimal phase-response curve of solar radiation input
Light stimuli generally affect the oscillatory dynamics multiplicatively, i.e. they act on the rate constants or transcriptional efficiency of the gene regulatory circuits [3,40]. We assume that the jth molecular species includes a parameter r as F j ðx; rÞ ¼F j ðxÞ þ rx k ; ð3:1Þ whereF j ðxÞ represents the terms that do not include r, and x k is the concentration of the kth molecular species. Here, k[f1,2, . . . ,Ng can take any value, regardless of j (both j = k and j ¼ k are allowed). For example, let figure 4c be a gene regulatory circuit of a hypothetical circadian clock, where symbols ! and s represent positive and negative regulations and x i are molecular species (see Novák & Tyson [41] for typical motifs of biochemical oscillators). Suppose x 1 and x 2 are mRNA and corresponding protein, respectively, and light stimuli increase the translational efficiency. In this case, the dynamics of light entrainment can be described by equation (3.1) with j ¼ 2, k ¼ 1 and r being the translation rate. In equation (3.1), although we can also consider an alternative case F j ðx; rÞ ¼F j ðxÞ À rx k (a negative sign), the optimal pPRCs remain unchanged under the inversion which is seen from equations (2.21) and (2.23). Consequently, we consider only the positive case to calculate the optimal PRCs (i.e. equation (3.1)). However, note that relations between iPRCs and pPRCs are affected by the inversion of the sign, and the difference matters when considering biological feasibility. When using phase reduction, the dynamics of the limit cycle are considered on the unperturbed limit-cycle trajectories x LC , and hence the points on the limit cycle can be uniquely determined by the phase f. Consequently, under the phase reduction, x k is replaced by x LC,k (f ) in equation (3.1), where x LC,k (f ) is the kth coordinate of x LC (i.e. @ r F j ðf; rÞ ¼ x LC;k ðfÞ in equation (2.20)). Here, x LC;k ðfÞ corresponds to the time course of the concentration of the kth molecular species. Because x LC,k (f ) constitutes a core clock component and is generally a smooth 2p-periodic function, we approximate it with a sinusoidal function: where u is the initial phase and a denotes the amplitude of the oscillation (figure 4d). To ensure x LC;k ðfÞ ! 0, we set 0 a 1, and a ¼ 0 recovers the additive case. Because the initial phase u does not play any role (u is offset by d in equation (2.23)) when the white Gaussian noise is additive (i.e. Q i (x) / 1), we also set u ¼ 0. The parametric approximation of equation (3.2) enables an almost closed form for the overall calculations. Although we assumed in equation (3.1) that effects of r only depend on x k , we can generalize equation (3.1) to F j ðx; rÞ ¼F j ðxÞ þ rKðxÞ, where K(x) is a nonlinear function and is assumed to be well approximated by 1asin(f þ u). By this generalization, our theory can be applied to other possible light entrainment mechanisms such as the intercellular coupling [42]. Our model needs only details about molecular species which have light input entry points but not about a whole molecular network. However, this advantage, in turn, means that we cannot specify iPRCs U i (f) of molecular species not having light input entry points. Consequently, for a noise term Q i (x), we assume that the white Gaussian noise is additive and is present only in the jth coordinate equation where q is the noise intensity and Q i (f) ¼ 0 for i = j). Figure 5a-c shows the landscape of C(D,d) as functions of D and d, and figure 5d-f expresses the optimal iPRCs U j (f) and pPRCs Z(f) for the solar radiation input (an explicit expression of C(D, d) is given in appendix A). The optimal PRC shape does not depend on the model parameters such as the period T, its variance s 2 T , or noise intensity q. These three parameters only act on the magnitude of the PRCs (i.e. the vertical scaling of the PRCs). Consequently, we normalized T ¼ 1, s 2 T ¼ 1, and q ¼ 1, as shown in figure 5. As the optimal PRCs depend on a, C(D, d) is plotted for three cases: a ¼ 0, (figure 5a), a ¼ 0.5 ( figure 5b) and a ¼ 1.0 (figure 5c), where the maximal points (D, d) yield the optimal PRCs using equations (2.20) and (2.21). The maximal parameters D and d are calculated numerically. Figure 5d-f describes the optimal iPRCs (solid line) and pPRCs (dashed line) for a ¼ 0, 0.5 and 1.0, respectively. When a ¼ 0, i.e. the input signal is additive, C(D,d) achieves a maximum for D ¼ p and arbitrary d, yielding sinusoidal PRCs as the optimal solution (figure 5d). Although the input signal p(f) is not sinusoidal, the optimal PRCs obtained using the variational method become sinusoidal. In other words, considering optimality, resonator-type oscillators have an advantage over integrator-type oscillators. For a . 0, the input signal p(f) depends on the concentration of the kth molecular species. From figure 5b, the optimal parameters for a ¼ 0.5 are (D, d) ¼ (2.31, 1.99) and (3.98, 4.30), which are different from D ¼ p (these two sets yield symmetric PRCs with respect to the horizontal axis). Figure 5e shows the optimal iPRCs U j (f) and pPRCs Z(f) for a ¼ 0.5. Interestingly, the optimal iPRCs and pPRCs for a ¼ 0. 5

Dead zone length
From the results discussed above, the optimal PRCs have a dead zone when a . 0. We next studied the length of the dead zone as a function of a (figure 6a) and improvements in the entrainability induced by the dead zone (figure 6b) for the solar radiation input. Because the dead zone, which is a null interval in PRCs, emerges when the optimal parameter is D = p, we can naturally define its length as where D is the maximum value of C(D, d). As seen in figure 6a, a dead zone clearly exists when a . 0, and the length increases with increasing a for a , 0.8. Even for a ¼ 0.1, when the oscillation amplitude of x LC,k (f ) (the concentration of a molecular species modulated by the lightsensitive parameter r; cf. figure 4d) is very small, we observe rsif.royalsocietypublishing.org J. R. Soc. Interface 11: 20131018 a dead zone with a length of L ¼ 0.475, which corresponds to about 3 h within 24 h, indicating the universality of having a dead zone in order to attain optimality. The improvement in the entrainability that is induced by a dead zone is calculated by comparing the entrainability of the optimal PRCs with that of typical sinusoidal PRCs. We consider sinusoidal functions for both the iPRC U j (f ) and pPRC Z(f ) by setting where c is the parameter to be optimized so that entrainability is maximized for each a (see appendix B for the explicit expressions). Equations (3.4) and (3.5) are scaled, so that they satisfy the constraints on the period variance (equation (2.12)). We calculated the ratios where E U and E Z represent the entrainabilities for the cases of the sinusoidal iPRC and pPRC, respectively, calculated for the solar radiation input. For the sinusoidal iPRC of equation (3.4), the entrainability is calculated with pPRC via equation (2.14). R U and R Z quantify the improvement rate of the optimal PRCs over the sinusoidal iPRC (R U ) and pPRC (R Z ). In figure 6b, the solid and dashed lines show R U and R Z , respectively, as a function of a. Both ratios monotonically increase as a increases, which shows that the optimal PRC with a dead zone exhibits better entrainability when the oscillation of x LC,k (f ) has a larger amplitude.  . In (d ), a parallel shift of the PRC is also optimal (d can be an arbitrary value). In (e), symmetric PRCs with respect to the horizontal axis are also optimal. In ( f ), symmetric PRCs with respect to the horizontal axis or f ¼ 3p/2 are also optimal (see the text).
The pPRCs correspond to experimentally observed PRCs. When the concentration of x LC,k (f ) is low, the effects of the input signal on the circadian oscillators are smaller. This is because pPRC Z(f ), which quantifies the extent of the phase shift owing to the stimulation of the parameter, depends on the concentration x LC,k (f ) (see equation (2.14)). However, even within the range f where x LC,k (f ) has smaller values, the iPRC U j (f ) contributes to an increase in the variance of the period, regardless of the concentration. From this, we see that having an iPRC with a smaller magnitude when the concentration of x LC,k (f ) is smaller results in a smaller variance, which results in a larger entrainability for a constant variance of the period. Although this qualitatively explains the benefit of a dead zone, for some input values, the optimal PRCs may not contain a dead zone for any value of a. This will be shown in the following.

Optimal phase -response curve of sinusoidal input
Because the optimal PRCs depend on input signals (equations (2.20) and (2.21)), we next consider a typical periodic input signal, a sinusoidal function (equation (2.25)). In this case,

Discussion
The existence of a dead zone optimizes both entrainability and regularity. It is rather obvious that optimization of regularity alone leads to a dead zone [43], because null response means no effect by any kind of fluctuations. Our result instead shows that optimality of both entrainability and regularity, which are in a trade-off relationship, is uniquely achieved by a dead zone. Our finding is fairly general, because a dead zone always exists in an optimal PRC unless a ¼ 0 (additive stimulation). Along with the fact that T, s T and q affect only the scaling of the optimal PRCs, when the input signal affects the dynamics multiplicatively (i.e. a . 0), the existence of a dead zone always provides a synchronization advantage. This is supported by many experimental studies of various species that report the existence of a dead zone in the PRC [1] (cf. figure 2c,d). Our general result suggests that circadian oscillators have fully adapted to solar radiation to improve synchronization. Indeed, many experimental findings imply that circadian oscillators have adapted to actual solar radiation [44]: for various animals, light-dark (LD) cycles that include a twilight period result in better entrainability than do abrupt LD cycles (on-off protocols) [44]. In this regard, another interesting problem is optimal entrainment [37] of circadian clocks by light stimuli. As two different input signals, the solar radiation and sinusoidal inputs, yield the same optimal PRCs for a ¼ 0, optimal inputs and optimal PRCs do not have one-toone correspondence. Thus, the optimal inputs are not trivial and this problem should be pursued in our future studies. The solar radiation input plays an essential role, because it yields a dead zone in the optimal PRC, whereas a sinusoidal signal does not ( figure 7). In other words, oscillators that are entrained by stimuli other than solar radiation may not exhibit a dead zone in their PRCs. This is indeed found in mammals. Mammals possess a master clock in their suprachiasmatic nucleus (SCN), which receives light stimuli via retinal  photoreceptors, and peripheral clocks in body cells [45]. The peripheral oscillators are entrained by several stimuli such as feeding and signals from the SCN through chemical pathways (e.g. hormones) [45,46]. By injection experiments of a hormone, Balsalobre et al. [47] reported that the PRCs of the peripheral oscillators in the liver do not have a dead zone.
Our result also agrees with other experimental observations. Our theory implies that a dead zone should be located where the concentration x LC,k (f) is low (0 f p in figure 4d), and that to achieve optimality, the concentration of x LC,k (f) should be maximal in the region where the PRCs exhibit a large phase shift. In Drosophila, the timeless (tim) gene is regarded as the molecular implementation of x LC,k (f). It is experimentally known that light enhances the degradation of the gene product (the TIM protein) [48,49], and the TIM protein peaks during the late evening. Figure 2c shows observations of the PRC of Drosophila against light pulses as a function time from Hall & Rosbash [19]; circles describe the experimental data, and the solid line expresses a trigonometric curve fitting (fourth order). Because the centre of the part of the PRC that can be phase shifted approximately corresponds to the peak of the concentration, as denoted above, when estimated from the PRC alone, the concentration peak of the TIM protein should occur at about 18 h. This time is also close to the experimental evidence (i.e. late evening). Therefore, our theory can be used to hypothesize further molecular behaviour affected by light stimuli.
In summary, we have constructed a model that regards circadian oscillators as a global optimization of entrainability and regularity. We have shown that our model is consistent with much experimental evidence as mentioned above. The extension and improvement of our method are possible and they are left as an area of future study. where C b ðD; dÞ ¼ C a ðÀD þ 2p; Àd þ 2pÞ. We show equation (A 1) as functions of D and d in figure 5d-f.

A.2. Sinusoidal input case
For the sinusoidal input case (equation (2.25 where we used equation (2.14).

B.2. Sinusoidal parametric phase-response curve
For the pPRC Z(f ) to be a sinusoidal function, the iPRC U j (f ) must be Equation (B 4) is normalized, so that the period variance becomes V T ¼ s 2 T . Using equation (2.14), the corresponding pPRC is a sinusoidal function: which is an explicit expression of the sinusoidal pPRC (equation (3.5)).