Mathematical Analysis of Light-sensitivity Related Challenges in Assessment of the Intrinsic Period of the Human Circadian Pacemaker

Accurate assessment of the intrinsic period of the human circadian pacemaker is essential for a quantitative understanding of how our circadian rhythms are synchronized to exposure to natural and man-made light-dark (LD) cycles. The gold standard method for assessing intrinsic period in humans is forced desynchrony (FD) which assumes that the confounding effect of lights-on assessment of intrinsic period is removed by scheduling sleep-wake and associated dim LD cycles to periods outside the range of entrainment of the circadian pacemaker. However, the observation that the mean period of free-running blind people is longer than the mean period of sighted people assessed by FD (24.50 ± 0.17 h vs 24.15 ± 0.20 h, p < 0.001) appears inconsistent with this assertion. Here, we present a mathematical analysis using a simple parametric model of the circadian pacemaker with a sinusoidal velocity response curve (VRC) describing the effect of light on the speed of the oscillator. The analysis shows that the shorter period in FD may be explained by exquisite sensitivity of the human circadian pacemaker to low light intensities and a VRC with a larger advance region than delay region. The main implication of this analysis, which generates new and testable predictions, is that current quantitative models for predicting how light exposure affects entrainment of the human circadian system may not accurately capture the effect of dim light. The mathematical analysis generates new predictions which can be tested in laboratory experiments. These findings have implications for managing healthy entrainment of human circadian clocks in societies with abundant access to light sources with powerful biological effects.


Reduction of typical Kronauer-type models to phase-only models
In the absence of zeitgebers, all Kronauer-type models can be written in the general form ẍ + μk(x) ẋ + x = 0, (S1) where time has been rescaled, such that t = ωt, where ω = 2π/τ c .(See Table S1 for the equivalent unscaled form for four popular versions of the model).For small μ, it can be shown that equation (S1) has stable almost sinusoidal solutions x ≈ r(t) cos(t + ϕ(t)) where r(t) is slowly varying and of O(1) [1] i.e. a stable limit cycle.
Defining y = ẋ, and changing to polar coordinates equation (S1) can also be written as Changing to polar coordinates using the definitions x = r cos ϕ, y = r sin ϕ gives For all the Kronauer-type models, close to the limit cycle, r and k(r cos ϕ) are O(1) and μ is small.Hence φ ≈ −1, i.e. a phase-only oscillator.
The presence of a zeitgeber introduces forcing terms to the van der Pol oscillator equations, but nevertheless, the Kronauer-type models may still be represented in the general where F represents the effect of light on the oscillator.The forms of f (x, y), g(x, y) and F for four widely used versions are given in Table S1.We note that we have changed variables from definitions in the original papers for ease of comparison between models.For example, for [2], x F ≡ −y and y F ≡ x.
Again transforming to polar coordinates, we find JFK [4] St Hilaire [5] Oscillator equation π Stimulus modulator Table S1: Comparison of four different versions of Kronauer's model.
where B is the stimulus strength.For dim-light dark cycles, the magnitude of the forcing term is small so r(t) ≈ 1, hence where and u(ϕ, 1) is the equivalent of the velocity response curve in phase only models.
For the four models listed in Table S1, the forms of u(ϕ, 1) are for Kronauer's 1998 model [3], for [4] and [5], and for [2].We note that [5] has an additional non-photic stimulus with velocity response of the form u(ϕ, 1) ≈ 1 + tanh(10 cos ϕ) sin ϕ. (S8) Across one cycle, these approximate velocity response curves have mean values of 0, 0.075, 0.0684 and 0 respectively.We note that the O(µ) terms in equation (S4) are to lowest order odd functions of ϕ and independent of r therefore do not contribute to the lowest order terms in the correction to τ F D .

Evaluating the cumulative phase response
We have, Substituting for R(ϕ) in equation (S9) using the Fourier expansion for the angular velocity response curve, equation (3) gives where f = M/T and T = ωT and we have used the fact that Similarly, for general k > 0, we have Hence, The challenge now is to evaluate the sum over k on the right hand side of equation (S11).We assume that the intrinsic period τ = 2π/ω is close to 24 h, that is, τ = T solar (1 + δ), where T solar = 24 h and |δ| ≪ 1. (For example, when τ = 24.5 h, δ = 0.020.)Then, Hence, from the map given in equation (8) of the main manuscript, For example, when T = 28 h, equation (S12) gives and when T = 20 h, equation (S12) gives FD protocols are carried out for a specified number of cycles, where the length of the cycle T and the number of cycles N are frequently chosen so that that N T is an integer number of days.For T = 28 h or T = 20 h this means that protocols typically last a multiple of six cycles.The significance of equation (S12) is that it shows that for N cycles such that N T is an integer number of days, then the FD protocol will result in ϕ 0 , ϕ 1 , . . ., ϕ N −1 where the ϕ k are approximately uniformly distributed around the circle.
For evenly spaced ϕ k = ϕ 0 + 2kπT T solar , a standard calculation can be used as follows.Consider Focussing on the sum term, and letting α = 2πT T solar , and similarly for the real parts of the sum.
FD protocols are chosen so that they are a number of beat cycles, so that N T is an integer multiple of T solar .Hence α = 2πN T T solar is an integer multiple of 2π and for any j.It follows that,

Analytical expression for the phase of entrainment for a sinusoidal velocity response curve
Equation ( 11) is Taking the tangent of equation (S28) and using the angle difference identity where Rearranging equation (S30) gives where f = M/T , A = L/ω, and n is the number of times the clock traverses ϕ = 0 during the photoperiod.Rearranging and taking the tangent of equation (S35), we obtain Applying the angle difference identity in equation (S29) to the left-hand side of equation (S36) and substituting for ϕ using equation (S34) gives where Equation (S37) can be rearranged to obtain a quadratic equation in tan(ϕ 0 /2): Since tan(ϕ 0 /2) must be real, if there are no real solutions to equation (S40), then entrainment is not possible with the given parameters.
After finding solutions to equation (S40) where they exist, we determine their stability using the condition P ′ (ϕ 0 ) < 0. In simple clock models, the derivative of the phase response curve is Once a stable solution for ϕ 0 has been found, the phase of entrainment can be calculated.