Resynchronization of circadian oscillators and the east-west asymmetry of jet-lag

Experiments show that a region of the brain within the hypothalamus regulates the circadian rhythms of mammals.^1,2 This region, called the suprachiasmatic nucleus (SCN), contains of the order of 10⁴ pacemaker neuronal cells that exhibit oscillatory dynamics. The synchronization of these pacemaker neuronal cells is responsible for maintaining healthy circadian rhythms. Studies on SCN slice cultures from neonatal mice in which the coupling between pacemaker cells is suppressed show that the individual natural periods of SCN pacemaker cells are distributed with a mean slightly longer than 24 h.³ Other experiments demonstrate that the behavior of the SCN is externally affected by light induced signals transmitted from the eye and that these signals can entrain the dynamics of the SCN to a 24 h day.^3–5 

Numerous models have been formulated to study the dynamics of the SCN.^6–16 Many of these previous models rely on assumed forms for the macroscopic dynamics of the SCN,^6–9 while others model the SCN as many coupled pacemaker oscillators and then numerically compute behavior of their many oscillator systems.^10–16 Here, we start with a simple microscopic description of the coupled dynamics of individual pacemaker cells and then make use of an analytical reduction of the dynamics to a low dimensional, macroscopic description. Thus, our approach has the important advantage of simplicity and ease of interpretation provided by the availability of a relatively simple low dimensional macroscopic model, while retaining the realistic aspect of being based on the dynamics of many coupled oscillators. As we will show, this approach yields considerable insight, although perhaps at the expense of having to employ rather simplified modeling of the individual oscillators. Specifically, we consider a sinusoidally forced Kuramoto model in which the forcing models diurnally periodic external inputs (e.g., light), and we employ the technique of Refs. 17 and 18 to extract a macroscopic description from the microscopic model. We are specifically interested in using this model to study how human circadian rhythms are affected by cross-time-zone travel (i.e., the phenomenon of jet-lag) and how SCN pacemaker cells become re-entrained to a new time-zone.

We model the SCN using $N≫1$ globally coupled phase oscillators with individual phase angles θ_i (⁠$i=1,...,N$⁠). The coupling strength between pacemaker cells is assumed to be uniform for all connections and denoted by K. We model the 24 h period external light source by introducing a temporally sinusoidal external driving term. The evolution for each oscillator is given by

where ω_i is the natural frequency for oscillator i when it is decoupled from other oscillators, $σ=2π/24 (h)−1$ is the day-night frequency of the external drive, and F is its amplitude. The phase quantity p is determined by the time-zone in which the individual is located and takes values between $−π$ and π. The variables satisfying Eq. (1) may be interpreted as follows. The oscillator phases advance in time with “fully awake (noon),” for example, corresponding to the set of times $θi(t)=2πn$ with n integer. If all oscillators have a 24-h period, $ωi=σ$⁠, Eq. (1) indicates that each phase is attracted to the same value, $θi=σt+p(t)$⁠. If time t is defined to be the Greenwich time, then p = 0 corresponds to being in Greenwich. Local noon to the east of Greenwich occurs earlier than noon in Greenwich. Thus, p > 0 to the east and p < 0 to the west of Greenwich. During travel p changes with time. For simplicity, in our study we assume that all travel is very fast and approximate it as instantaneous so that we will take p(t) to change discontinuously from one constant value to another where the time of discontinuity corresponds to the time of travel. See Refs. 19 and 20 for previous work on the dynamics of Eq. (1) (not specifically in the context of jet-lag, i.e., with p(t) = 0).

With our fast travel assumption

where the fast travel takes place at time τ and the travel starts from a time-zone with phase p₁ and ends in another time-zone with phase p₂. We further suppose that the system has been in the p = p₁ state for a long time before the travel takes place, so that, by the time $t=τ$⁠, the system has evolved to its p = p₁ attractor. Thus we model recovery from jet-lag by following the solution of Eq. (1) for $t≥τ$ with p = p₂ from an initial condition on the attractor of Eq. (1) with p = p₁.

We take the oscillators' natural frequency distribution $g(ω)$ to be Lorentzian

where Δ and ω₀ denote the spread and the center of the natural frequency distribution, respectively. Similar to previous studies on forced Kuramoto models,^17,20 we reduce the $N≫1$ dimensional dynamics of the system to the dynamics of a complex order parameter z in a rotating frame with frequency σ

The dynamics of this order parameter for $N≫1$⁠, and $t≥τ$ is governed by the following equation:²⁰ 

where $Ω=σ−ω0, ż=dz(t)/dt$⁠, and Ω = 0 corresponds to the case where the center of the natural frequencies of the oscillators are the same as the external frequency σ. Since many experimental studies show longer-than-24-h-free-running-period of human circadian rhythms,^21–26 the average frequency of SCN pacemaker cells in typical humans is thought to have a period longer than 24 h, and thus Ω is positive.

If an individual is in a steady entrained state before travel, and at time $t=τ$ travels from a time-zone with p = p₁ to a timezone with p = p₂, then, upon his/her arrival, the individual's order parameter is $zstei(p1−p2)$⁠, where $zst$ is the entrained state in the new time-zone, i.e., $zst$ is the stable fixed point of Eq. (5), assuming that it exists. Note that Eq. (5) and $zst$ do not depend on the phase of the destination time-zone, p₂.

Because an individual in the new time-zone experiences misalignment between his/her biological time and the local time, recovery from cross-time-zone travel takes place in the following days as the order parameter gradually evolves back to the stable fixed point $zst$⁠. That is, recovery from jet-lag is described in Eq. (5) for $t≥τ$ with the initial condition

where $Δp>0$ for westward travel, $Δp<0$ for eastward travel, and traveling to a time-zone with a 1 h difference in clock setting corresponds to $|Δp|=2π/24$⁠. During recovery, in what follows, we will quantify the level of the deviation from full recovery by the distance between the order parameter and the stable fixed point, $|z(t)−zst|$⁠.

The model we have studied is highly simplified, but, as previously noted, offers the advantage of having an individual oscillator basis and reducible to a macroscopic description. Some aspects that are not included but have experimental basis are spatially structured coupling within the SCN,³ both attractive and repulsive interactions of SCN pacemaker cells,¹³ and potential influence on diurnal timing from other organs of the body.²⁷ 

The attractors of Eq. (5) have been studied in Refs. 19 and 20. It was shown that, with different sets of the parameters, K, F, Δ, and Ω, the system exhibits different dynamics separated by bifurcations.

Figures 1(a)–1(c) show dynamical phase flow portraits for the three main types of dynamics for Eq. (5). Following Antonsen et al.,¹⁹ we label these three types type A, type B, and type C, where each type occurs in its corresponding region of $(K/Δ,F/Δ,Ω/Δ)$ parameter space, and the three panels of Fig. 1 are flows computed at specific parameters sets in the appropriate region.

For type A dynamics (Fig. 1(a)), the system has one fixed point that is stable. For type B dynamics (Fig. 1(b)), there is a limit cycle that is stable and an unstable fixed point enclosed within the limit cycle. For type C dynamics (Fig. 1(c)), there are three fixed points, one stable, one unstable, and one of saddle type, and the unstable manifold of the saddle flows to the stable fixed point.

By definition (4) the magnitude of the order parameter satisfies $|z|≤1$⁠; thus, we are only interested in the dynamics within the circle $|z|=1$ which is shown as the dashed circles in Figs. 1(a)–1(c). Also from Eq. (4), larger $|z|$ corresponds to greater synchrony, and $|z|=1$ corresponds to perfect synchrony (i.e., all the θ_j are equal).

Note from Figs. 1(a) and 1(c) that the stable fixed point lies below the positive real axis. This reflects the parameter choice $Ω>0$⁠, corresponding to the typical situations for humans (i.e., a slower SCN frequency, represented by ω₀, as compared to the 24 h driving, represented by the frequency σ). Note that this negative phase of $zst$ corresponds to a human circadian response whose oscillation phase lags the forcing oscillation phase.

In type B dynamics, the attractor is a limit cycle and thus not fully entrained to the forcing. This case would correspond to an individual whose awake/sleep times naturally vary from day to day, even if the driving (e.g., sunlight) is the same every day. We note, however, that if type B dynamics applies with parameter values near the bifurcation point to type A dynamics, then the limit cycle is a relatively small loop around the unstable fixed point that was stable before the bifurcation. In this case, there may be little practical difference between type B dynamics and a case with a stable fixed point.

In our subsequent parameter studies, for simplicity, we will investigate parameter regions where the dynamics is of either type A or type C. Thus type B dynamics will not be further considered in this paper.

For the same number of time-zones crossed, several studies have indicated that recovery from jet-lag is typically faster following westward travel as compared with eastward travel.^28–32 A primary goal of our modeling efforts is to understand the underlying reason for this asymmetry.

Naively, we might expect that when an individual travels eastward crossing less than 12 time-zones, his/her intrinsic phase is delayed compared to the local time phase of the destination, and thus the individual re-entrains his/her intrinsic clock by phase advancing. Similarly, we might expect that an individual who travels across less than 12 time-zones westward should re-entrain his/her intrinsic clock by phase delaying. However, counter intuitively, “wrong direction” re-entrainment has been found for large enough eastward travels,^8,33,34 where, for example, an individual who travels eastward crossing 8–10 timezones re-entrains his/her intrinsic clock by phase-delaying, instead of phase advancing.

In the following, we discuss this temporal and directional asymmetry in jet-lag recovery using the dynamics of our model. We first illustrate type C dynamics in Fig. 2(a) in which we show order parameter trajectories resulting from trips crossing different numbers of time-zones. We plot four trajectories corresponding to jet-lags from traveling across 9 and 12 h westward and across 8.5 and 9.5 h eastward for the parameters given in the figure caption, where the initial condition for each trip is determined by using Eq. (6). (In Subsections III B and III C, we discuss exploration of the parameter space.) The four trajectories in Fig. 2(a) are respectively plotted as blue, cyan, green and red curves superposed on a copy of the phase space plot of Fig. 1(c) (plotted in grey). For the westward 9 h trip, the order parameter (in blue) evolves clockwise, decreasing its phase so as to realign to the stable fixed point. For the eastward 8.5 h trip, the order parameter (in green) evolves counterclockwise approaching the stable fixed point. In contrast, for the eastward 9.5 h trip, the order parameter moves clockwise. This asymmetry results from the presence of the saddle point. Specifically, although the eastward 8.5 and 9.5 h trips have close initial starting points in the order parameter space, they lie on opposite sides of the stable manifold of the saddle point, and thus they flow in two different directions. The crossover point between delays and advances of the phase lies on the stable manifolds. An eastward journey with a time-zone shift pass the crossover point results in phase-delaying recovery, while an eastward journey with a time-zone shift smaller than the crosser-over point results in phase-advancing recovery, as would typically be expected of eastward travel. We note that according to our model, an initial condition exactly on the stable manifold of the saddle would never flow to the stable fixed point and initial conditions close to the stable manifold could take a very long time to approach the vicinity of the stable fixed point (meaning recovery from jet-lag would correspondingly be very drawn-out). In a real situation, however, we do not expect arbitrarily long recovery times due to factors omitted from the model, such as day-to-day deviations in the external drive term due, for example, to light fluctuations resulting from varying weather patterns.

Analogous to Fig. 2(a) (which applies for parameters where there is type C dynamics), Fig. 2(b) shows order parameter trajectories for a set of parameters for which type A dynamics applies. In going from Fig. 2(a) to Fig. 2(b), the repeller fixed point and the saddle fixed point, along with the stable and unstable manifolds of the saddle, have annihilated. Nonetheless, the order parameter trajectories for Figs. 2(a) and 2(b) are similar. This is because of the local character of the saddle-node bifurcation and because the parameter values of Figs. 2(b) and 2(a) are not too far from the saddle-node bifurcation point. Note that due to the absence of the saddle, type A dynamics of the model cannot lead to arbitrarily long jet-lag recovery times (in contrast with type C dynamics).

We now attempt to estimate a reasonable set of the parameters (K, F, Ω, and Δ) that might be thought to correspond to a typical healthy person. However, we recognize that very substantial variation of these parameters is indicated by the observation that the severity of jet-lag has large variation from one individual to another. Thus, we will use our putative “typical” parameter set as a reference point from which we will study the effect of parameter deviations. Our choice of a typical parameter set is based on the following:

1. Pacemaker cells in the SCN should be in a fairly synchronized state, which suggests that the stable fixed point of the order parameter has a fairly large amplitude (i.e., $|z|<1$ should not be small compared with one).

2. Empirical studies (Refs. 22 and 35) suggest that an individual confined to a dark enclosure may have substantial (but noisy) regularity in his sleep-awake cycle. Thus we surmise that the system can synchronize when F = 0.

3. No direct measurements of SCN oscillators are available for human, and thus any estimation of the oscillator distribution must be based on indirect deductions. Observations of blind people²¹ and of an individual in a dark enclosure,²² as well as other earlier experimental results,^23,24 suggested that typical circadian rhythms have a period very roughly in the range 24.5–26 h when there is no external drive. However, it has been argued²⁵ that such observations are not indicative of the natural period of the relevant biological oscillators, and, using a different methodology, Ref. 25 has estimated an average period of roughly 24.2 h. However, another, more recent study,²⁶ estimated a period of 24.5 h. Here, we will take 24.5 h as value for the average period.³⁶ 

4. The speed for an individual to recover from jet-lag is different for eastward and westward travels.³² For example, it is sometimes claimed that, on average, an individual approximately retrains his/her phase by 1 h per day after an eastward travel and retrains his/her phase by 1.5 h per day after a westward travel. Also, 9-h eastward travel typically corresponds to the border between jet-lag recovery by phase advancing or delaying,^37,38 and thus is considered most disruptive.

Based on the above and on numerical simulations of our model (to incorporate the information from point 4), we choose

as our reference parameter set. (For this set type A dynamics applies.) The corresponding recovery curves from several travels eastward and westward are shown in Fig. 3, where the vertical axis is $|z(t)−zst|$⁠, which we regard as a quantitative measure of the severity of jet-lag at time t. We take recovery from jet-lag to have occurred when $|z(t)−zst|$ drops below 0.2 (the horizontal black line). We see that for this parameter set there is a strong asymmetry in recovery, with westward recovery being significantly faster, and that the longest recovery time occurs for 9 h time-zone eastward travel.

As previously mentioned, the effects of cross-time-zone travel vary substantially from individual to individual. With this motivation, we investigate the parameter dependence of the model dynamics and its implications for recovery from jet-lag.

Dividing Eq. (5) by Δ, we see that keeping $K/Δ, F/Δ$ and $Ω/Δ$ constant while varying the value of Δ is equivalent to changing the time scale of the dynamics. Thus, in the following, we keep Δ constant at our reference value, Eq. (7), and study how the dynamics change when K, F, and Ω are individually varied away from the reference parameter set (Eqs. (8)–(10)).

The parameter K, which represents the oscillator coupling strength, controls the degree of synchronization for fixing values of the external drive, F. In Fig. 4(a), we show the recovery time in a range of K for several cross-time-zone journeys with the other two parameters at their reference values (⁠$F=3.5Δ$ and $Ω=1.4Δ$⁠). For $K/Δ$ between 2.0 and about 8.0, there is only one fixed point that is stable (type A dynamics). The magnitude of the stable fixed point increases as K increases (not shown), corresponding to increasing strength of the degree of synchronization. When K exceeds $K≈8.0Δ$⁠, two other fixed points, one unstable and one saddle, appear through a saddle-node bifurcation (type C dynamics). The large recovery time for eastward 9-h travel is due to the newly born saddle point (see Figs. 1 and 2) that is very close to the eastward 9-h travel initial point. In simulations with other parameter sets (not shown), a sharp, narrow spike in the recovery time can occur in plots like those in Fig. 4(a) when, as $K/Δ$ is varied, the stable manifold of the saddle point crosses the initial point of the jet-lag simulation, causing the recovery time to go to infinity.

The parameter F controls the magnitude of the external drive which is meant to capture the effects of factors such as the sensitivity of an individual's response to sunlight, cloudiness, geographical latitude, and seasonality. In Fig. 4(b), we show the recovery time in a range of F for several cross-time-zone journeys with the other two parameters at their reference values (⁠$K=4.5Δ$ and $Ω=1.4Δ$⁠). In the parameter range illustrated, the dynamics are always type A. Not surprisingly, the recovery times for all seven cross-time-zone journeys decrease with increasing F. For increasingly large values of F, the stable fixed point approaches the positive real z-axis (not shown), reflecting decrease in the difference between the external drive phase and the internal overall phase. The decrease of recovery time with increasing F is consistent with the idea that recovery from jet-lag is facilitated by being outdoors in bright daytime sunlight (effectively increasing F).

The parameter $Ω=σ−ω0$ is the difference between the external drive frequency σ (corresponding to a 24 h day) and the average of the pacemaker cells' natural frequencies ω₀. In Fig. 4(c) we show the recovery time in a range of Ω from $−3.7Δ$ to $3.7Δ$ for several cross-time-zone journeys with the other two parameters at their reference values (⁠$K=4.5Δ$ and $F=3.5Δ$⁠). From Fig. 4(c), we notice that the eastward and the westward recovery time curves (dashed and solid curves of the same color) are reflections of each other about the $Ω/Δ=0$ axis. This can be seen from the fact that the complex conjugate of Eq. (5) with Ω replaced by $−Ω$ yields an equation for $z*(t)$ that is the same as the original equation (Eq. (5)) for z(t).

In this paper, we present a model of circadian dynamics that is based on the microscopic dynamics of individual pacemaker cells in the suprachiasmatic nucleus (SCN). Our model employs an extremely simple phase oscillator description of the individual SCN cells. This description, while not reflecting the true complexity of the mechanisms governing the individual cells (see, e.g., Refs. 10 and 11), has the advantage that in the large N limit (⁠$N=$ the number of pacemaker cells), it allows an exact reduction of the high dimensional microscopic dynamics to that of a relatively simple low dimensional system for the evolution of the macroscopic system state. Using this low dimensional macroscopic description, we model the dynamics of resynchronization of the SCN following rapid cross-time-zone travel, i.e., jet-lag. We believe that this approach has resulted in increased understanding of the east/west jet-lag asymmetry and of the origin of potential differences in the responses of different individuals to cross-time-zone travel.

While our approach provides insights into the key features of jet-lag from a dynamical systems perspective, a number of avenues remain for further exploration. For example, our basic uncoupled phase oscillator (⁠$dθi/dt=ωi$⁠) could be replaced by other, potentially more realistic, phase oscillator descriptions,^39,40 while still allowing a reduction to low dimensional dynamics. Also, because many features concerning the coupling of SCN pacemaker cells are unknown, as well as for the sake of simplicity, we have assumed an all-to-all network topology for the connection between SCN pacemaker cells. We note, however, that experience with other networked systems suggests that heterogeneous connectivity patterns may play an important role in the resynchronization of circadian oscillators in individuals experiencing jet-lag. Finally, we have implemented rapid cross-time-zone travel as a discontinuous phase shift in a sinusoidal external drive term, but more complex external drive functions may be appropriate and could even be used, for example, to explore the effectiveness of different strategies for combating jet-lag.

This work was supported by the Army Research Office under Grant No. W911NF-12-1-0101, and by the National Science Foundation under Grant No. PHY-1461089.