Optimal implementations for reliable circadian clocks

Circadian rhythms are acquired through evolution to increase the chances for survival through synchronizing with the daylight cycle. Reliable synchronization is realized through two trade-off properties: regularity to keep time precisely, and entrainability to synchronize the internal time with daylight. We found by using a phase model with multiple inputs that achieving the maximal limit of regularity and entrainability entails many inherent features of the circadian mechanism. At the molecular level, we demonstrate the role sharing of two light inputs, phase advance and delay, as is well observed in mammals. At the behavioral level, the optimal phase-response curve inevitably contains a dead zone, a time during which light pulses neither advance nor delay the clock. We reproduce the results of phase-controlling experiments entrained by two types of periodic light pulses. Our results indicate that circadian clocks are designed optimally for reliable clockwork through evolution.

Many terrestrial species, from cyanobacteria through to humans, adapt to sunlight and have acquired circadian oscillatory systems.Although their molecular implementation is species dependent [1], all exhibit a regular rhythm of 24-h that can be entrained by light input.Two fundamental properties are necessary for circadian clocks [2], i.e. regularity to keep time precisely [3][4][5] and entrainability to adjust the internal time through light stimuli [6][7][8].It is not easy to maximize entrainability and regularity simultaneously; higher entrainability implies more vulnerability to noise, whereas higher regularity signifies less flexibility to entrainment.We studied an optimal phase-response curve (PRC) problem for a one-input pathway in [9].In the present Letter, we generalize the calculations of [9] and formalize the optimal implementation problem for multiple input pathways without relying on specific oscillator models [7,10,11].We show that the simultaneous maximization of entrainability and regularity entails several inherent properties of circadian clocks.At the molecular level, we rationalize the role sharing of multiple input pathways which was reported in murine circadian clocks [12,13].At the behavioral level, we rationalize a time period during which the time is neither advanced nor delayed by light in an optimal clock [14].Our theory can also explain different gene expression patterns when entrained by two different light pulses [15].In this study we investigate the optimal implementations for multiple input pathways to achieve maximal limit of entrainability and regularity.Although all circadian clocks transmit light signals through multiple pathways [7,10,11], entrainment problems in multiple pathways have been little studied.
Circadian clocks keep time regularly under molecular noise [3][4][5].It is also entrainable by periodic light stimuli.To consider these effects, we use an N -dimensional Langevin equation with respect to x i , which is the concentration of ith molecular species: ẋi = F i (x; ρ) + Q i (x)ξ i (t) where F i (x; ρ) is the ith reaction rate (i = 1, 2, .., N ), Q i (x) is a multiplicative noise term, and ξ i (t) is white Gaussian noise with ξ i (t)ξ j (t ′ ) = 2δ ij δ(t − t ′ ).To incorporate the effect of light, we introduce a light-sensitive parameter ρ which is perturbed as ρ → ρ + dρ when stimulated by light.We quantify regularity by the temporal variance in the oscillation.From [9], the period variance is known to be 2 dθ, where T is the period of the unperturbed oscillator, U (φ) = (U 1 (φ), • • • , U N (φ)) is the infinitesimal PRC (iPRC) defined by U (φ) = ∇ x φ| x=xLC(φ) and x LC (φ) is a point on the limit-cycle trajectory at phase φ ∈ [0, 2π).Entrainability is quantified as the width of the Arnold tongue.The phase dynamics with periodic input signals are described by a tilted periodic potential.If stable points exist in the potential, the oscillator can be entrained by input signals [16].Thus we can discuss entrainability without considering noise because the existence of stable points does not depend on noise.Therefore we set Q i (x) = 0. Let p(ωt) be an input signal with angular frequency ω.For a weak input signal dρ = χp(ωt) (χ ≥ 0 is the signal strength), we obtain φ = Ω + χZ(φ)p(ωt) (φ ∈ [0, 2π) is the phase and Ω = 2π/T ) where with F i (φ; ρ) = F i (x LC (φ); ρ).Z i (φ) quantifies phase shift due to the perturbation of ρ in the ith coordinate.We hereafter refer to Z(φ) as the parametric PRC (pPRC) [9,17].The dynamics of a slow variable defined that exhibit peaks at different phase (this difference is given by ν).(c) The phase difference ν is defined by the difference between k1th and k2th molecular species.
Optimal circadian clocks are derived as optimal PRCs U i (φ) which were obtained with the variational method by maximizing the entrainability [19]) [9].We generalize the above scheme to multiple input signals.Suppose there are M input pathways and parameters ρ 1 , ρ 2 , • • • , ρ M are affected as ρ i → ρ i + dρ i under light stimuli.By introducing an auxiliary scaling parameter s i defined by dρ i = s i dρ, each parameter perturbation is reparameterized as where ρi = ρ i −s i ρ.Through this reparameterization, the multiple perturbations can be reduced to a single parameter case with respect to ρ.In this instance, ρ in Eq. ( 3) is a hypothetical parameter and has no correspondence to actual reaction rates.
Light stimuli generally affect the rate constants, but the detailed mechanisms vary between organisms [21].Our model can cover these different entrainment mechanisms.We first explain a simple case (Fig. 1(a)), and then generalize it.In Fig. 1(a), we assume that light stimuli enhance the synthesis rate of x 2 , where the dynamics of x 2 can be described by ẋ2 = ρ syn x 1 + ρ deg x 2 where ρ syn > 0 and ρ deg < 0 are the synthesis and degradation rates, respectively.Because the synthesis rate increases when stimulated by light, by taking ρ = ρ syn (ρ is the light-sensitive parameter), the rate equation can be divided into ρ-dependent and ρ-independent parts ẋ2 = F 2 (x; ρ) = ρ deg x 2 +ρx 1 = F2 (x)+ρx 1 , where F2 (x) denotes terms not including ρ in F 2 (x; ρ).We next consider a generic case with two-input pathways (M = 2).Let us assume that light stimuli affect parameters ρ 1 and ρ 2 (e.g.translation, transcription or degradation rate) which depend on k 1 th and k 2 th molecular species and affect j 1 th and j 2 th species, respectively (e.g.k 1 = 1, k 2 = 2, j 1 = 2 and j 2 = 3 in Fig. 1(b)).The rate equations F ji (x; ρ) of j 1 th and j 2 th are described as ρ-independent ρ-dependent where i = 1, 2 and Fji (x) denotes terms not including ρ i in F ji (x; ρ).From Eq. ( 1), we see that PRCs are only concerned with the derivative of F ji (x; ρ) with respect to ρ, which is specifically given by Consequently, the ρ-independent part in Eq. ( 4) plays no role in the formation of optimal PRCs.We approximate x k by the kth coordinate of x LC , denoted as x LC,k (φ) in Eq. ( 4).Because x LC,k (φ) oscillates, we model it with a sine curve.k 1 th and k 2 th molecular species may peak at different times as shown in Fig. 1(b), where the time courses of the k 1 th and k 2 th molecular species is described in the insets by solid and dashed lines, respectively.Therefore, the time course is approximated by , where u k1 and u k2 are initial phases of k 1 th and k 2 th molecular species, the amplitude α i being assumed as identical α 1 = α 2 = α (0 ≤ α ≤ 1).Let ν be the phase difference between k 1 th and k 2 th molecular species, i.e. ν = u k2 − u k1 (Fig. 1(c)).This parameter plays an important role in optimization.
Circadian clocks are entrained by sunlight whose intensity is determined by 24-h periodic solar irradiance.The solar irradiance I with respect to the zenith angle ϑ is given by I = I 0 cos ϑ, where I 0 is the irradiance at ϑ = 0, and I vanishes when the sun is below the horizon [22].Thus we approximate its time course by p(ωt) = sin(ωt) for 0 ≤ mod (ωt, 2π) < π (day) and p(ωt) = 0 for π ≤ mod (ωt, 2π) < 2π (night) with mod (x, y) = x mod y, which we call solar radiation input [9].
Our model starts from the two-input case (M = 2).The noise term was introduced only for input molecules , where q is the noise intensity), and we set u k1 = 0 (i.e.ν = u k2 ), T = 1, σ 2 T = 1, and q = 1 without loss of generality.We then need to specify additional parameters s 1 and s 2 (Eq.( 3)) that determine the strength of the input signal relative to ρ.If we are concerned with pPRCs (experimentally observed PRCs) only, the sign of s i (positive or negative) does not play any role, because s i is squared in the optimal pPRCs.We set s 1 = s 2 = 1 for simplicity; this corresponds to assuming equal weight for each light input pathway.
Figure 2(a) shows the ν dependence of maximal entrainability for α = 0.5 (solid line) and α = 1 (dashed  and (b), Z(φ), Zj 1 (φ) and Zj 2 (φ) obtained with the variational method, are plotted with thick solid lines, dashed lines and dot-dash lines, respectively.Zj 1 (φ) and Zj 2 (φ) obtained with a numerical method (NM) are plotted with circles and triangles, respectively [19].In (a), PRCs that are symmetric with respect to the horizontal axis or φ = 3π/2 are also optimal.In (b), PRCs that are symmetric with respect to the horizontal axis are also optimal.(c)-(d) Time course of the molecular species concentration for (c) α = 1 and ν = 0 and (d) α = 1 and ν = 1.47.Two species x LC,k 1 and x LC,k 2 are plotted with dashed and dot-dash lines.In (a) and (c), dashed and dot-dashed lines are indistinguishable.line).For α = 0.5, the maximum was achieved at ν = 0 where no phase difference existed.However, upon increasing α to 1, maximal entrainability was achieved at two points, ν = 0 and ν = 1.47 (Figs.3(a) and  (b)).Interestingly, optimality can be attained in the presence of phase difference.From Eq. (1), we divide the pPRC Z(φ) into contributions from two input pathways Z(φ) = Z j1 (φ) + Z j2 (φ), where Z j1 (φ) and Z j2 (φ) quantify the phase shift produced by the 1st and 2nd input pathways, respectively.Optimal Z(φ), Z j1 (φ), and Z j2 (φ) for the two-input case are plotted in Figs.3(a) and (b).Figures 3(c) and (d) are the corresponding time course of the k 1 th and k 2 th molecular species concentration.We see that optimal PRCs Z(φ) in Figs.3(a) and (b) are very similar to experimentally observed PRCs in which there is a dead zone, a time during which light neither advanced nor delayed the clock (1 φ 2 in Figs.3(a) and (b)).Intriguingly, experimental studies in different species reported the existence of the dead zone [14].Although cases with ν = 0 and ν = 1.47 achieved the same entrainability, Z(φ) for ν = 0 is asymmetric with respect to horizontal axis, which entails an asymmetric Arnold tongue.Thus for a symmetric Arnold tongue, only ν = 1.47 can achieve maximal entrainability.We calculated the ν dependence of the dead zone length L (length of null parts in PRCs [19]) in Fig. 2(b) for α = 0.5 (solid line) and α = 1.0 (dashed line).In Fig. 2(b), L quickly diminishes around ν = π, showing that a dead zone always appears in optimal PRCs except for a singular point ν = π.
The fundamental difference of ν = 1.47 from ν = 0 is the role sharing of two PRCs Z j1 and Z j2 (Fig. 3(b)).Z j1 is responsible for the phase advance and Z j2 for the delay (the positive part of Z j1 is larger than the negative part and vice versa).This effect was observed for all ν values except ν = 0, as shown below.We quantify the distance between Z j1 and Z j2 by dist(Z j1 , Z j2 ) = ´2π 0 {Z j1 (θ) − Z j2 (θ)} 2 dθ, which becomes larger when the two PRCs play more compensatory roles.The distance calculated as a function of ν for α = 0.5 (solid line) and α = 1 (dashed line) is shown in Fig. 2(c) where the distance is maximal exactly at ν = π.When a phase difference (ν > 0) exists, this role-sharing between twoinput pathways always yields a synchronization advantage.There is experimental evidence for advance and delay roles of Per 1 and Per 2, respectively, in mice [12].In this regard, Ref. [13] observed a period dependence of Per 1 and Per 2 knockout mutants on the intensity of constant light.We can reproduce this result with optimal PRCs as follows.Note that entrainability is maximal at ν = 1.47 for α = 1 (Fig. 3(b)).Consequently, we set the clock parameters to α = 1 and ν = 1.47.Under a constant light condition, the input signal is modeled by dρ = χp(t), where p(t) = 1.By integrating the phase equation φ = Ω + χZ(φ)p(t) from t = 0 to t = T where φ(t = 0) = 0, the phase at time T with input strength χ is given by φ(T ; χ) = 2π + (2π) −1 T χ ´2π 0 Z(θ)dθ.For weak χ, the period T χ , which is the period under a constant light condition, is approximated by [17,19].Assuming x j1 = [P er1] and x j2 = [P er2], we simulated Per 1 and Per 2 mutants by setting Z(φ) = Z j2 (φ) and Z(φ) = Z j1 (φ), respectively.When increasing the intensity χ of constant light, the period ratio T χ /T increases for Per 1 mutant and decreases for Per 2 mutant.This result agrees with the experimental evidence (Fig. 2 in [13]).
Our model can further suggest insights into the molecular mechanism of the clock.In hamsters, Schwartz et al. [15] reported different gene expressions of Per 1 and Per 2 when entrained by two types of periodic light pulses that have short (23.33-h) and long (24.67-h) periods.Let us reproduce Schwartz's experiment in our optimization framework with two inputs.We again set α = 1 and ν = 1.47 and assume x j1 = [P er1] and x j2 = [P er2].
Ref. [15] applied a periodic light pulse of 1-h duration, which we modeled with a periodic δ-function p(ωt) = 2πδ(mod(ωt, 2π)), (5) where a factor 2π ensures Θ(ψ) = (2π) −1 ´2π 0 p(θ − ψ)Z(θ)dθ = Z(ψ).Given the periodic light pulse (Eq.( 5)), the entrainment phase ψ st (i.e. the circadian time at which hamsters receive the light pulses) can be determined by where we used Eq. ( 2).Thus ψ st can be given as a solution of Z(ψ st ) = (ω−Ω)/χ with Z ′ (ψ st ) < 0. For the long (ω < Ω) and short (ω > Ω) pulses, (ω − Ω)/χ becomes negative and positive, respectively.This shows that the long and short pulses always act on hamsters at early (φ = 1.3 ∼ 3.2; purple in Figs.4(a) and (b)) and late (φ = 4.7 ∼ 0.47; orange in Figs.4(a) and (b)) subjective night, respectively.The effects of the light pulse on the circadian clock depend on the concentration of x k1 and x k2 at these phases (Eq.( 4)).For the long pulse, we obtain x LC,k1 < x LC,k2 (Fig. 4(b)), which indicates that the long pulse always affects the expression of x j2 whereas it influences x j1 only a little.In contrast, the short pulse affects x j1 more strongly than x j2 .Our result shows that, provided the circadian clocks are designed optimally, long and short pulses affect the expression of two different components (Per 2 and Per 1) differently.Surprisingly, our expression patterns agree with the experiments of Schwartz et al. [15].They hypothesized that light stimuli affect the transcription of Per 1 and the degradation of Per 2. In their molecular terms, k 1 th and k 2 th species in our framework correspond to Bmal 1 and Per 2, which regulate the light effect (transcription and degradation) on Per 1 and Per 2, respectively (i.e.x k1 = [Bmal1] and . The phase difference between Per 2 and Bmal 1 was experimentally determined as ν ∼ 2 [24] and close to our result (ν = 1.47).
The pPRCs hitherto discussed are intrinsic in the sense that they represent the internal clock dynamics.The intrinsic pPRCs can be observed only through the phase shift induced by short light pulses [14].Theoretically, precise measurement is possible only through δ-peaked stimuli.In experiments involving higher organisms, however, light pulses are much longer than the δ-peaked function, and observed pPRCs become different from the intrinsic ones.To study the relation between intrinsic and observed pPRCs, let us consider a squared-pulse stimu- where H(t) is the Heaviside step function, t s is onset time of the pulses and ℓ is the pulse duration.For ℓ → 0, the squared pulse reduces to a δ-function δ(t − t s ).Let Z(φ; ℓ) be an observed pPRC of Z(φ) by a light pulse with the duration ℓ.Observed and intrinsic pPRCs can be related via c µ = −ic µ ℓµΩ/{χ(1 − exp(iµΩℓ))} for µ = 0 and c µ = cµ /χ for µ = 0, where c µ and cµ are Fourier coefficients of intrinsic and observed pPRCs, respectively (Z(φ) = Nµ µ=−Nµ c µ exp(iµφ) and Z(φ; ℓ) = Nµ µ=−Nµ cµ exp(iµφ) with N µ being an expansion order).By this method, we inferred the intrinsic pPRC from an observed pPRC in human [23] (ℓ = 6.7-h)where a dead zone is seemingly nonexistent [19].The inferred pPRC (dashed line) and the observed pPRC (solid line) are shown in Fig. 4(d), where the phase (horizontal axis) represents the onset of the pulse.This result suggests that superficial pPRCs may lack a dead zone even though their innate mechanisms actually do.
We have demonstrated that key properties of circadian clocks are consequences of optimization to attain the maximal limit of entrainability and regularity.Our theory explains known experimental results such as the role sharing of two inputs and different gene expression patterns by different pulses.We also explain the superficial absence of a dead zone in human.The model can be used to reveal key molecular elements responsible for the clock.This work was supported by Grant-in-Aid for Young Scientists B (Y.H.: No. 25870171) and for Innovative Areas "Biosynthetic machinery" (M.A.) from MEXT, Japan.

Yoshihiko Hasegawa and Masanori Arita
This supplemental material describes detailed calculations introduced in the main text.Equation and figure numbers in this section are prefixed with S (e.g., Eq. (S1) or Fig. S1).Numbers without the prefix (e.g., Eq. ( 1) or Fig. 1) refer to numbers in the main text.

Optimization
We calculated the optimal phase-response curves (PRCs) using two approaches, variational and numerical methods, whose explicit procedures are described below.

Variational method
Regularity V T and entrainability E are calculated in our previous study [1].Regularity is defined by the period variance of the oscillation (higher regularity corresponds to smaller period variance), which is represented by Here N is the dimension, U i (φ) is the infinitesimal PRC (iPRC), T is the period of the unperturbed oscillator and Q i (φ) is a multiplicative noise term (see the main text).Entrainability, the extent of synchronization to a periodic input signal, is defined by the width of the Arnold tongue.Given a periodic input signal p(ωt) (ω is angular frequency), entrainability is represented by where Θ(ψ) = (2π) −1 ´2π 0 Z(ψ + θ)p(θ)dθ, ψ M = argmax ψ Θ(ψ), ψ m = argmin ψ Θ(ψ) and Z(φ) is the parametric PRC (pPRC) defined in Eq. ( 1).
The optimal PRCs, which maximize entrainability E under constant regularity V T = σ 2 T , can be calculated by the Euler-Lagrange variational method.A variational equation to be optimized is where λ is a Lagrange multiplier.With Eqs. ( 1), (S1) and (S3), Eq. (S4) is specifically given by where F i (φ; ρ) is the ith reaction rate and ρ is a light-sensitive parameter.The functional derivative of L with respect to U i yields Figure S1: Piecewise linear PRC (K = 6) as an approximation of the iPRC U i (φ), where the lth knot is V i,l (circle).We set V i,K = V i,0 for a periodic condition.
Solving Eq. (S6), we obtain the optimal PRCs U i (φ) and Z(φ) given by where the Lagrange multiplier λ is A substitution of Eq. (S8) into Eq.(S3) yields the entrainability E represented by We maximize E as functions of ∆ and δ, where ∆ = ψ M − ψ m and δ = ψ m .Substituting the obtained parameters into Eqs.(S7) and (S8), we can calculate the optimal PRCs [1].The dead zone length L is defined by the length of null parts in the PRCs.Because a dead zone emerges when ∆ = π in Eqs.(S7) and (S8), we can naturally define its length by which is plotted in Fig. 2(b) as a function of the phase difference ν.

Numerical method
Because the optimal iPRC U i (φ) and parameters (ψ M and ψ m ) are interdependent, the variational method is not trivial.Consequently, to verify the correctness of the variational method numerically, we calculated the optimal PRCs with the evolutionary algorithm (specifically, we used a differential evolution (DE) of Storn and Price [2] provided by MATHEMATICA 9).Dividing the iPRC U i (φ) into a piecewise linear function, we can reduce the variational problem to a conventional multivariate parameter optimization problem.We divide the iPRC U i (φ) into K regions (Fig. S1) and linearly connect each knot to create a piecewise linear function as an approximation to U i , where the value of the lth knot is given by V i,l for l = 0, 1, .., K − 1 (we set V i,K = V i,0 for a periodic condition): where ρ is a light sensitive parameter (Eq.(S33) corresponds to an additive case).The intrinsic pPRC of the oscillator is described by dashed line in Fig. S2(b).For test data, we artificially generated an observed pPRC which is the phase difference ∆φ caused by the squared-pulse with the duration ℓ.We applied ρ → ρ + dρ to Eq. (S33) where ρ = 0 and dρ = χp(t; t s ) (ℓ = 2 and χ = 0.1) and the phase difference ∆φ as a function of onset phase of the pulse is used as the observed pPRC which is plotted in Fig. S2(b) (solid line; the magnitude is multiplied by 10).The dot-dashed line in Fig. S2(b) shows an inferred intrinsic pPRC through Eq. (S32) which is calculated with N µ = 3 order approximation (i.e.we approximated the observed pPRC with the third order Fourier series and applied Eq. (S32)).We see excellent agreement between the inferred (dot-dashed line) and true (dashed line) intrinsic pPRCs, which verifies the reliability of the inference method.

( a )FIG. 1 :
FIG.1: Examples of (a) one-light input pathway and (b) twolight input pathway cases.The insets describe the time course of molecular species that are affected by light stimuli.In (b) the light-sensitive molecular species correspond to x1 and x2 that exhibit peaks at different phase (this difference is given by ν).(c) The phase difference ν is defined by the difference between k1th and k2th molecular species.