Measuring Relative Coupling Strength in Circadian Systems

Abstract

Modern imaging techniques allow the monitoring of circadian rhythms of single cells. Coupling between these single cellular circadian oscillators can generate coherent periodic signals on the tissue level that subsequently orchestrate physiological outputs. The strength of coupling in such systems of oscillators is often unclear. In particular, effects on coupling strength by varying cell densities, by knockouts, and by inhibitor applications are debated.

In this study, we suggest to quantify the relative coupling strength via analyzing period, phase, and amplitude distributions in ensembles of individual circadian oscillators. Simulations of different oscillator networks show that period and phase distributions become narrower with increasing coupling strength. Moreover, amplitudes can increase due to resonance effects. Variances of periods and phases decay monotonically with coupling strength, and can serve therefore as measures of relative coupling strength.

Our theoretical predictions are confirmed by studying recently published experimental data from PERIOD2 expression in slices of the suprachiasmatic nucleus during and after the application of tetrodotoxin (TTX). On analyzing the corresponding period, phase, and amplitude distributions, we can show that treatment with TTX can be associated with a reduced coupling strength in the system of coupled oscillators. Analysis of an oscillator network derived directly from the data confirms our conclusions.

We suggest that our approach is also applicable to quantify coupling in fibroblast cultures and hepatocyte networks, and for social synchronization of circadian rhythmicity in rodents, flies, and bees.

Keywords

circadian clock oscillator coupling synchronization suprachiasmatic nucleus

Introduction

Networks of coupled oscillators are ubiquitous in living nature and physical systems (Pikovsky et al., 2001; Strogatz, 2003). Therein, even weak coupling can lead to a synchronization or pattern formation (Glass and Mackey, 1988; Winfree, 1990). Often, oscillator networks constitute incoherent states (i.e., oscillators have independent periods and random phases) or synchronized states due to overcritical coupling strength (Kuramoto, 2003).

In chronobiology, dispersed fibroblasts represent an incoherent state (Nagoshi et al., 2004), whereas the suprachiasmatic nucleus (SCN) exhibits complete synchronization in wild-type animals under standard conditions (Yamaguchi et al., 2003). The observation of frequency locking in the synchronized state allows a binary classification associated with undercritical and overcritical coupling.

There are, however, chronobiological questions requiring a more detailed quantification of coupling strength: do increasing cell densities strengthen coupling between fibroblasts even if no global synchronization is achieved (Noguchi et al., 2013)? Do inhibitors of neuronal communication weaken coupling within the SCN despite persisting synchronization (Ono et al., 2016)? Is there a dependence of coupling strength on the number of cohabitating animals in systems with social synchronization (Paul et al., 2015)?

In all these cases, modern recording techniques allow to monitor rhythms of individual oscillators. In this work, we suggest to study the corresponding distributions of individual periods, phases, and amplitudes to quantify coupling strength beyond a binary classification into incoherent and coherent states. It will be shown that period distributions are helpful indicators in the case of undercritical coupling, whereas narrowing of phase distribution and amplitude expansions can serve as markers of coupling strength in synchronized states.

Circadian systems constitute networks of coupled oscillators on different levels. Gene-regulatory feedback loops generate oscillations in almost every cell. Synchronization of these noisy oscillations can lead to precise self-sustained circadian rhythms on the organismic level.

In mammals, synchronized neurons in the SCN orchestrate physiological rhythms of the whole body with an astonishingly high precision. The network of SCN neurons illustrates that single-cell oscillations can be quite different from network oscillations. Synchrony, precision, and large amplitudes arise as emergent properties due to coupling (Herzog et al., 2004, 2017; Welsh et al., 2010; Mohawk and Takahashi, 2011).

There is also evidence that other cell types such as fibroblasts or hepatocytes can be coupled. It has been shown, for example, that cell densities of fibroblast cultures affect rhythmicity (O’Neill and Hastings, 2008; Noguchi et al., 2013). In many cases, however, the coupling strength is not sufficient to allow complete synchronization.

Moreover, there are reports on social synchronization of circadian rhythmicity in weakly interacting mice (Paul et al., 2015), hamsters (Paul et al., 2014), flies (Levine et al., 2002), and bees (Fuchikawa et al., 2016). Here, tactile contacts, pheromones, and common environmental factors are discussed as coupling mechanisms (Bloch et al., 2013).

In all these examples of coupled circadian oscillators, the type and strength of coupling are debated. Close to the synchronization threshold, the Kuramoto theory provides a sophisticated framework to quantify synchronization (Kuramoto, 1975). For example, the synchronization index R, which is tantamount to the global phase coherence, is widely used to characterize the onset of synchronization (see, e.g., Strogatz, 2000).

The situation is less clear for overcritical coupling as found in the SCN or undercritical coupling as in dispersed fibroblasts. In these cases, wide ranges of coupling strength lead to practically the same phenotype: complete synchronization or desynchronization, respectively.

We present in this article a conceptional framework of how measured data can be exploited to estimate coupling strength between oscillators. We illustrate the concept using well-studied networks of amplitude phase oscillators (Bodenstein et al., 2012; Gu et al., 2013) and an established experimental protocol to pharmacologically reduce circadian coupling. Ensembles of coupled amplitude phase oscillators are systematically analyzed to extract features that are associated with increasing coupling strength. These features are based on distributions of periods, phases, and amplitudes. We show that narrowing of period and phase distributions as well as amplitude expansion are related to an increasing coupling strength. Although an explicit calculation of coupling strength from ensembles of oscillators is currently impossible, a quantitative comparison between situations in which putative coupling mechanisms are differentially manipulated can be performed. In the case of SCN networks, for example, the effects of GABA inhibitors could be tested; and, in fibroblasts, cultures of low and high cell densities can be compared.

As a first step, we analyze simulated networks to test the following hypotheses that have been previously described in the dynamical systems literature:

For undercritical coupling (desynchronized oscillators), the period distribution becomes narrower with increasing coupling strength (Liu et al., 1997). It is known from engineering that frequencies can approach each other due to coupling (frequency pulling) (Wiener, 1961). Furthermore, when increasing coupling, some of the oscillators might constitute frequency-locked clusters (Hoppensteadt and Keener, 1982).

For overcritical coupling (complete synchronization), the phase distribution becomes narrower with increasing coupling strength. Even in the case of complete frequency locking, phases of individual oscillators differ due to varying intrinsic frequencies or amplitudes (Strogatz, 2000). On a population level, these different phases are termed “chronotypes” (Roenneberg et al., 2003). It is plausible that for strong synchronizing signals, the phase dispersion shrinks.

Coupling can induce amplitude expansion (“resonances”). Whereas rigid limit cycles can be described successfully as phase oscillators (Ashwin et al., 2016), amplitudes of weak oscillators near the Hopf bifurcation can be amplified by synchronizing agents. This implies that an observation of resonant behavior also provides information on the underlying oscillator type.

After studying simulated networks, we will analyze data from brain slices representing the network of SCN neurons. We exploit local fits with amplitude-phase models to quantify period, phase, and amplitude distributions for 3 different situations: 1) synchronized SCN mouse explants, 2) application of TTX to reduce spike-associated coupling, and 3) washout of TTX to restore coupling. It is reasonable to assume that these situations represent strong, weak, and medium coupling, respectively (Yamaguchi et al., 2003; Baba et al., 2008; Webb et al., 2009, 2012; Enoki et al., 2012; Evans et al., 2013; Taylor et al., 2017).

Materials and Methods

Numerical Simulations

Numerical solutions of Equations (11)–(12), have been obtained for a total integration time of $400 \times 24 h$ at ∆ $t = 0.25 h$ time steps by using the odeint function of SCIentificPYthon. Initial conditions of the $N = 1000$ oscillators have been chosen to be equally distributed among the unit circle, that is, the radial component was set to $r_{i} (t = 0) = 1$ , whereas $θ_{i} (t = 0)$ has been drawn from uniformly distributed random numbers between $[0, 2 π]$ . The initial condition was chosen once and has then been fixed throughout all simulations.

Analysis of Simulated Time Series

The instantaneous phase $θ (t)$ and amplitude $A (t)$ of an oscillatory time series $s (t)$ are estimated by means of an analytic signal approach as described in more detail in, for example, Pikovsky et al. (2001). The analytic signal is defined as $z (t) = s (t) + i ℋ (s (t))$ , where $ℋ (s (t))$ denotes the Hilbert transform. Along these lines, the instantaneous phase and amplitude are defined by $θ (t) : = atan2 (ℑ (z (t)), ℜ (z (t)))$ and $A (t) : = \sqrt{ℜ {(z (t))}^{2} + ℑ {(z (t))}^{2}}$ , respectively, where $ℜ (z (t))$ and $ℑ (z (t))$ denote the real and imaginary parts of $z (t)$ . Throughout this article, Hilbert transformations of discrete time series are performed numerically by the SCIentificPYthon function hilbert.

A further analysis of the instantaneous phases and amplitudes was restricted to a time interval between $t \in T = [150 \times 24 h,250 \times 24 h]$ to ignore transient dynamics of the system and to minimize edge effects during the numerical estimation of the Hilbert transform (Pikovsky et al., 2001). Based on the definition $\frac{d θ_{i} (t)}{d t} = ω_{i} (t) = \frac{2 π}{τ_{i} (t)}$ of the instantaneous frequency of oscillator $i$ , we calculate the time-averaged period (as found, e.g., in Figures 2 and 3A) in the interval $T$ via $〈 τ_{i}^{HT} 〉 = \frac{2 π}{a_{θ}}$ , where $a_{θ}$ is the slope of a linear fit to the (unwrapped) instantaneous phase $θ_{i} (t)$ in the interval $T$ . The average amplitude of an oscillator $i$ (as depicted, e.g., in Figure 3C) is simply determined by the time average of its instantaneous amplitude $A_{i} (t)$ over the time interval $T$ .

Analysis of Experimental Time Series

Previously published data from PER2::LUC recordings in coronal SCN slices, as published in Abel et al. (2016), are analyzed in the Experimental Data from the Suprachiasmatic Nucleus section. First, time series corresponding to bioluminescence recordings of single SCN neurons are detrended by means of a Hodrick–Prescott filter. To this end, we use the hpfilter function of the statsmodels Python module and choose a smoothing parameter of $λ = 0.05 {(\frac{24 h}{Δ t})}^{4}$ , where ∆ $t$ is the experimental sampling rate of the time series, as suggested in St. John and Doyle (2015). Second, we independently fit under all 3 experimental conditions the stochastic linearly damped oscillator model:

\frac{d x}{d t} = - γ x - \frac{2 π}{τ} y + ξ_{x}

\frac{d y}{d t} = - γ y + \frac{2 π}{τ} x + ξ_{y}

by means of an autocorrelation approach as described in Westermark et al. (2009). Here, $γ$ and $τ$ denote the radial relaxation rate and the intrinsic free-running period of the oscillator, respectively, whereas $ξ_{j}$ for $j \in {x, y}$ are independent, delta-correlated white-noise terms with zero mean, that is, $〈 ξ_{i} (t) ξ_{i} (s) 〉 = 2 D δ (t - s)$ and $〈 ξ_{i} (t) 〉 = 0$ . In essence, the autocorrelation function

C (t) = \frac{D}{γ} e^{- γ t} \cos (\frac{2 π}{τ} t),

which can be analytically obtained from Equations (1)–(2), is fitted to the numerically obtained autocorrelation function of the experimental time series of interests, obtained by means of the acovf function of the statsmodels Python module. By this means, we can estimate the oscillator parameters $γ$ , $τ$ , and the noise strength $D$ for a given experimental time series. The amplitude of the resulting oscillations as plotted in Figure 6C can be obtained via $A = \sqrt{D / γ}$ ; see Westermark et al. (2009) for details.

Simulations of the stochastic ordinary differential Equations (1)–(2), as depicted in Figure 5B–D, have been obtained by means of Euler’s method for stochastic equations using XPP-AUTO (Ermentrout, 2002; Schmal et al., 2014).

Fitting a Linearized Hopf Oscillator

In addition to the method described in the previous section, we also fit a (stochastic) linearized Hopf oscillator

\frac{d r}{d t} = γ (A - r) + ξ_{r}

\frac{d θ}{d t} = \frac{2 π}{τ} + ξ_{θ},

with $〈 ξ_{i} (t) ξ_{i} (s) 〉 = 2 D_{i} δ (t - s)$ for $j \in {r, θ}$ , to time-series data under all 3 conditions separately, using an autocorrelation approach as described in Westermark et al. (2009). As before, the corresponding autocorrelation function

C (t) = \frac{1}{2} (A^{2} + \frac{D_{r}}{γ} e^{- γ t}) e^{- D_{θ} t} \cos (\frac{2 π}{τ} t)

is fitted to the numerically obtained autocorrelation function of the experimental time series of interest. Because fitting Equation (6) is an ill-conditioned problem due to the summation of 2 exponentials, we set $γ = 0.1 h^{- 1}$ throughout all fits.

Networks of Linearized Hopf Oscillators

In the Network Model Based on Fitted Parameters section, we study mean-field coupled networks of linearized Hopf oscillators. The underlying equations read as

\frac{d x_{i}}{d t} = γ_{i} \frac{x_{i}}{r_{i}} (A_{i} - r_{i}) - \frac{2 π}{τ_{i}} y_{i} + M_{i}

\frac{d y_{i}}{d t} = γ_{i} \frac{y_{i}}{r_{i}} (A_{i} - r_{i}) + \frac{2 π}{τ_{i}} x_{i}

with $r_{i} = \sqrt{x_{i}^{2} + y_{i}^{2}}$ and $M = \frac{K}{N} Σ_{i = 1}^{N} x_{i} (t)$ .

Random Networks

Oscillator dynamics on a random network are investigated as follows: we assume that the probability of oscillator $j$ , influencing oscillator $i$ , in a network of $N$ oscillators is given by $p$ . Next, the mean-field term $M_{i}$ in Equations (7) and (11) is replaced by $M_{i}^{R} = k_{i, in}^{- 1} \sum_{j = 1}^{N} K_{i j} x_{j} (t)$ , where $K_{i j}$ is assumed to be $K$ if oscillator $j$ influences oscillator $i$ , and $0$ otherwise. $k_{i, in}$ denotes the number of oscillators that have an impact on oscillator $i$ (i.e., the in-degree of oscillator $i$ ). In summary, the network topology is a directed random graph, including self-loops; see Supplementary Figures S16A–C and S17A–C for examples.

Results

Analysis of Coupled Amplitude-Phase Oscillators

Model Equations

In complex systems of coupled oscillators, the emergent behavior at the network level is influenced by both the properties of the individual oscillators and the coupling topology and strength between the oscillatory units. At an abstract level, the oscillatory entities can be of diverse nature. In the case of “social synchronization,” each oscillatory unit describes the dynamics at the level of the organism. In the case of the SCN or fibroblast cultures, the individual oscillatory units correspond to cell-autonomous circadian clocks. The intrinsic dynamical properties of the single oscillators can be conveniently parameterized by means of the generic phase-amplitude model:

\frac{d r_{i}}{d t} = γ_{i} r_{i} (A_{i} - r_{i})

\frac{d θ_{i}}{d t} = \frac{2 π}{τ_{i}},

which is considered a simple description of a stable limit cycle in a 2-dimensional (2D) plane and is commonly termed the Poincaré oscillator (Glass and Mackey, 1988). Equation (9) determines the dynamical evolution of the radial coordinate $r_{i} (t)$ (i.e., the time-dependent distance from the origin), whereas Equation (10) determines the change of the angular coordinate $θ_{i} (t)$ ; see Figure S1 for an illustration. The parameters $τ_{i}$ , $A_{i}$ , and $γ_{i}$ denote the free-running period, amplitude, and amplitude relaxation rate of oscillator $i$ , respectively. In contrast to contextual (i.e., organism or even cell-type-specific) models, such conceptual or phenomenological oscillator models allow for a straightforward analysis and interpretation of how specific entrainment or synchronization behaviors rely on fundamental clock properties such as amplitudes or free-running periods (Beersma, 2005; Roenneberg et al., 2008; Gonze, 2011).

What kind of emergent network properties, not inherent to the isolated oscillatory units, arise due to mutual coupling, and how do these properties depend on the coupling strength? To test this, we consider an ensemble of $N = 1000$ mutually coupled Poincaré oscillators and quantify the emergent properties during mutual coupling by numerical simulations. For the sake of simplicity, we assume that all oscillators interact with each other through a mean field, as proposed in previously published models of the SCN network (Gonze et al., 2005; Locke et al., 2008). This kind of coupling tacitly implies a relatively fast diffusion of coupling agents (e.g., neuropeptides in the case of SCN neurons) compared to the ~24 h time scale of circadian free-running periods and an equally weighted contribution of each oscillator to the mean field (Ullner et al., 2009). The network dynamics in the presence of coupling can then be given by

\frac{d x_{i}}{d t} = γ_{i} x_{i} (A_{i} - r_{i}) - \frac{2 π}{τ_{i}} y_{i} + M_{i}

\frac{d y_{i}}{d t} = γ_{i} y_{i} (A_{i} - r_{i}) + \frac{2 π}{τ_{i}} x_{i}

where we have rewritten Equations (9)–(10) in Cartesian coordinates and made, without loss of generality, the assumption that the mean field $M_{i} = \frac{K}{N}$ $Σ_{i = 1}^{N} x_{i} (t) = M$ additively couples solely to the x-coordinate, see Equation (11), as previously proposed (Gu et al., 2013; Schmal et al., 2015). Parameter $K$ denotes the strength of the coupling between the mean field and the single oscillatory units.

To get an intuition about the impact of increasing coupling strength, we investigate the dynamics of system (11)–(12) for 3 representative values of $K$ , as illustrated in Figure 1. In the corresponding simulations, we assume that all individual oscillators have an identical amplitude $A_{i} = 1$ and an amplitude relaxation rate of $γ_{i} = 0.1 h^{- 1}$ , as suggested by Westermark et al. (2009), Ramkisoensing et al. (2014), and Tokuda et al. (2015). The intrinsic free-running periods $τ_{i}$ are chosen from a normal distribution with mean $μ_{τ} = 24 h$ and a standard deviation (SD) of $σ_{τ} = 2 h$ , as suggested by experiments with dispersed (i.e., presumably uncoupled) SCN neurons (Welsh et al., 1995; Liu et al., 1997; Honma et al., 1998, 2004).

In the case of $K = 0.04$ , depicted in Figure 1A, no particular order can be observed, and all oscillators seem to run at or close to their own intrinsic frequency. In oscillator theory, such a (macro-)state is commonly referred to as the “incoherent state” (Strogatz and Mirollo, 1991). If we subsequently increase the coupling (e.g., to $K = 0.07$ ), order emerges: a huge fraction of the oscillators appears to run at the same pace, thereby leading to a dramatic increase of the mean-field oscillation amplitude. In this “partially synchronized state” or “mixed state,” a cluster of synchronized oscillators coexists with a nonsynchronized set of oscillators. Finally, large enough coupling (e.g., $K = 0.1$ ) leads to the emergence of a “fully synchronized state,” in which all oscillators are locked to the mean field; see Figure 1C. As expected from previous analysis (Abraham et al., 2010), an increasing modulation of the individual oscillators’ amplitude can be observed when increasing coupling strength $K$ (i.e., the oscillators show a resonance effect). In addition, it can be noticed that, within the synchronized cluster, oscillators with smaller free running periods $τ_{i}$ tend to phase-lead with respect to oscillators having larger values of $τ_{i}$ . This situation is analogous to the occurrence of chronotypes (i.e., the phase of entrainment of an organism with respect to an external zeitgeber). If synchronized to an external stimulus, the phase is analogously correlated with its internal free-running period (Aschoff and Wever, 1966; Pittendrigh and Daan, 1976; Granada et al., 2013; Schmal et al., 2015; Woelders et al., 2017).

Emergent network properties are fingerprints of coupling

The examples of the previous section show that certain dynamical properties such as amplitudes and periods of the individual oscillators are subject to change if the mutual coupling between the oscillators is varied. To study changes in amplitudes, phase organization, and periods in more detail, we investigate the oscillatory properties of every single oscillator for all 3 coupling strengths ( $K = 0.04$ , $K = 0.07$ , and $K = 0.1$ ) numerically. To this end, instantaneous phases and amplitudes are determined by means of a Hilbert transformation, as described in further detail in the Materials and Methods section. An average period $〈 τ_{i}^{HT} 〉$ is estimated for each oscillator by fitting a straight line to the (unwrapped) phase $θ_{i}^{HT} (t)$ and defining the inverse slope of this line times $2 π$ as $〈 τ_{i}^{HT} 〉$ . By this means, a single, time-independent value for the average period of each oscillator in the coupled state is determined, even if the dynamics show nonperiodic behavior such as beating; see Supplementary Figure S2 for an example analysis.

A direct comparison of the intrinsic free-running periods $τ_{i}$ in the uncoupled state and the corresponding average periods $〈 τ_{i}^{HT} 〉$ in the coupled state reveals a typical dependency between the size of the synchronized cluster and the coupling strength $K$ . In the case of small coupling strengths K, the average period $〈 τ_{i}^{HT} 〉$ of a given oscillator is barely affected by the mean field coupling, as one can see in Figure 2 for $K = 0.04$ (blue dots) (i.e., no frequency locking occurs). When a certain critical coupling strength is reached, a set of oscillators whose free-running periods are close to the mean of the ensemble period starts to form a frequency- or period-locked cluster, in which all oscillators share a common oscillation period, that is, they form a frequency plateau in the $(〈 τ_{i}^{HT} 〉, τ_{i})$ diagram; see Figure 2 for $K = 0.07$ (green dots). Due to the nonvanishing oscillation amplitude of the mean field in such cases (remember Figure 1B), even oscillators that are not locked to the synchronized cluster experience a change of their average period toward the period of the synchronized cluster, a phenomenon commonly known as “frequency pulling” (Wiener, 1961). For high coupling strengths K, the ensemble forms a giant synchronized cluster, in which all oscillators tick at a common pace; see Figure 2 for $K = 0.1$ (red dots). As a consequence, the variance in the distribution of the average individual oscillator periods $〈 τ_{i}^{HT} 〉$ decreases with increasing coupling; see Figure 3A.

As noted in the previous section, we can observe a re-organization of the oscillation phases with varying coupling strength. Although (instantaneous) oscillator phases $θ_{i}^{HT}$ are equally distributed among their whole co-domain $- π \leq θ_{i}^{HT} < π$ in the case of small coupling strength such as $K = 0.04$ , the variance of their distribution successively decreases with increasing $K$ ; compare Figure 3B. Moreover, we find an increase of the mean value in the distribution of average amplitudes $〈 A_{i}^{HT} 〉$ with increasing coupling due to resonance effects; see Figure 3C.

Finally, in Figure 4, we systematically explore the above-described findings for a broader range of coupling strength $K$ . Based on synchronization properties, we can divide the dynamics into 3 different regimes. In the case of undercritical coupling strength K, which approximately lies in the range of 0 ≤ $K ≲ 0.05$ , the oscillators do not synchronize to each other, and the SD in the distribution of periods is constantly close to the value of the free-running periods $σ_{τ} = 2 h$ ; see Figure 4A. In this incoherent state, phases spread all over the unit circle due to variations in the intrinsic periods, which leads to high values of the (circular) variance in the phase distribution and, hence, a constantly low value of the classical order parameter R; see Figure 4B. After reaching a critical coupling strength $K_{c} \approx 0.05$ , the fraction of synchronized oscillators increases, leading to a decreasing variance in the period and phase distributions and an increasing global phase coherence $R$ ; see Figure 4A and 4B. For overcritical coupling strength $K ≳ 0.1$ , all oscillators form a fully synchronized cluster such that the variance in the distribution of periods drops to zero; see Figure 4A. The variance in the phase distribution and the order parameter $R$ are close to zero and one, respectively, and slowly saturate to these values for a further increase in K; see Figure 4B. Overall, it can be noticed that the evolution of the global phase coherence $R$ for increasing $K$ closely resembles the behavior as known from the phase oscillator description in the paradigmatic Kuramoto model (Kuramoto, 1975; Strogatz, 2000; Sonnenschein and Schimansky-Geier, 2013); see Supplementary Text 1.1 and Supplementary Figure S3 for further details. Pure phase oscillator models, however, per definitionem exclude an analysis of amplitude effects and are thus insufficient in this study, as we will see in the following paragraphs.

Similar to the oscillation period, the amplitude of the oscillators experiences no change in the regime of undercritical coupling; see Figure 4C. As soon as the coupling exceeds the threshold $K_{c} \approx 0.05$ , synchronization with the mean field leads to an amplitude expansion with increasing $K$ . From symmetry considerations, we can estimate the resonance effect on amplitudes as

μ_{A} (K) = A + \frac{K}{4 γ} (1 + \sqrt{1 - {(\frac{2 Δ ω}{K})}^{2}})

in the case of overcritical coupling (i.e., for all K ≥ $2 Δ ω$ ); see Supplementary Text 1.2 for a detailed derivation. Here, $µ_{A} (K)$ denotes the mean of the distribution of individual oscillator amplitudes in the network for a given $K \geq 2 Δ ω$ , whereas $A$ and $γ$ are the intrinsic amplitude and radial relaxation rate of the individual oscillators. Finally, $Δ ω$ equals $\frac{2 π}{τ + σ_{τ}} - \frac{2 π}{τ - σ_{τ}}$ , where $τ$ and $σ_{τ}$ are the mean and SD of the intrinsic free-running periods. A comparison of the numerical results (Fig. 4C, bold line) and the analytical approximation (Fig. 4C, dashed line) reveals a good agreement in the investigated ranges of coupling strength $K$ (see also Suppl. Figs. S4–S7).

Dependence on Single-Cell Properties

From Equation (13), it becomes clear that the amplitude not only increases with increasing coupling strength $K$ but also depends on the intrinsic properties of the individual single-cell oscillators (i.e., the radial relaxation rate $γ$ as well as the spread of the free-running periods $σ_{τ}$ ). Such dependence on intrinsic oscillator parameters can also be observed with respect to other emergent network properties.

The critical coupling strength $K_{c}$ increases with an increasing spread $σ_{τ}$ of the free-running period distribution (see Suppl. Fig. S5A), similar to what can be observed in the Kuramoto model; compare Supplementary Text 1.1 and Figure S3B. Contrarily, $σ_{τ}$ has little effect on the amplitude expansion for high coupling strength $K ≫ K_{c}$ , because Equation (13) converges to $μ_{A} (K) \approx A + \frac{K}{2 γ}$ for $K ≫ 2 Δ ω$ ; see also Supplementary Figure S5B.

For fixed values of $A_{i} = A$ and $σ_{τ} = 2 h$ , an increasing radial relaxation rate $γ_{i} = γ$ leads to a mild increase in the critical coupling strength $K_{c}$ (see Suppl. Fig. S6A) while having a dramatic effect on the amplitude expansion (see Suppl. Fig. S6B). Along these lines, rigid oscillators (e.g., for $γ_{i} = 1 h^{- 1}$ ) exhibit almost no variations in amplitude with coupling, whreeas nonrigid oscillators (e.g., for $γ_{i} = 0.01 h^{- 1}$ ) exhibit a considerable increase in amplitude; see Supplementary Figure S6B. This is tantamount to what has been observed in theoretical entrainment studies for nonrigid and rigid oscillators; see, for example, Abraham et al. (2010).

Finally, as also can be deduced from Equation (13), we can observe that an increasing intrinsic amplitude $A_{i} = A$ of the individual oscillators has virtually no effect on the critical coupling strength $K_{c}$ ; see Supplementary Figure S7A. Furthermore, an increasing intrinsic oscillator amplitude proportionally increases the oscillation amplitude at the network level while having little effect on the amplitude expansion with coupling; compare Supplementary Figure S7B.

In conclusion, different oscillator properties at the single-cell level can produce differential responses with coupling. Thus, the behavior at the network level for different coupling regimes might provide information on the underlying oscillator type.

Experimental Data from the Suprachiasmatic Nucleus

Studying simulated networks of mean-field coupled oscillators resulted in the following theoretical predictions: increasing coupling can lead to frequency locking, a narrowing period distribution, a decreasing variance of the phase distribution, and an increasing amplitude due to resonance effects. We will use these predictions to interpret data from slices of SCN tissue ex vivo. It has been shown that, even though single SCN neurons are able to generate circadian rhythms, the remarkable precision of the SCN as a tissue crucially depends on heterogeneous intercellular couplings (Herzog et al., 2004). Here, we focus on previously published data from PER2::LUC recordings in 5 coronal SCN slices of neonatal to 7-day-old mice, in which spike-associated couplings have been altered by pharmacological treatments with tetrodotoxin (TTX) (Abel et al., 2016). Therein, after 4 days of bioluminescence recordings following the slice preparations, 2.5 µM TTX was applied and remained in the medium for 6 days. Subsequently, TTX was washed out from the culture medium, and recordings were maintained for at least 8 more days; see Abel et al. (2016) and Taylor et al. (2017) for details. In Abel et al. (2016), single SCN neurons were automatically identified and tracked in the 2D bioluminescence time-lapse recordings. Ensemble average bioluminescence intensities (black line) and SDs (gray-shaded area) for each time point of the recordings of one SCN slice are plotted in Figure 5A. Time-series data that correspond to the other SCN slice preparations are plotted in Supplementary Figure S8. As already noted in Abel et al. (2016) and Taylor et al. (2017), application of TTX leads to a reduction in synchrony and a reduced amplitude of the SCN cells that can be partially reversed by washing out TTX. Here, we estimate oscillation properties of these single-cell time series separately under all 3 conditions (i.e., before TTX application, during TTX application, and after washing TTX from the medium). To this end, we fit a (stochastic) linearly damped oscillator model with 3 parameters to the data separately under all these conditions, using an autocorrelation function approach as described previously (Westermark et al., 2009; Tokuda et al., 2015); see the Materials and Methods section for further details. By this method, we determine for each neuronal cell its free-running period, amplitude, radial relaxation rate, and fluctuation or noise strength independently under every condition. In Figure 5B–D, representative simulated time traces for the estimated parameters under each condition are shown.

We can rediscover the characteristic patterns of different coupling regimes, as known from our simulation studies, when comparing period, phase, and amplitude distributions under all 3 experimental conditions; see Figure 6 for example distributions from the first SCN slice. Before application of TTX, the distribution of the single-cell periods has a mean of $μ_{τ} \approx 26.18 h$ and an empirical SD of $s_{τ} \approx 0.29 h$ (Fig. 6A, red). During the application of TTX, the ensemble mean of the period drops to a value of $μ_{τ} \approx 24.91 h$ , accompanied by a broadening of its distribution to $s_{τ} \approx 2.10 h$ (Fig. 6A, blue), lying close to the value that was previously published for distributions of dispersed SCN neurons (Herzog et al., 2004; Welsh et al., 1995; Honma et al., 1998). Finally, washing out TTX from the medium partially reverses this effect: the distribution’s mean rises to a value of $μ_{τ} \approx 25.46 h$ , whereas its empirical SD decreases to $s_{τ} \approx 0.35 h$ , lying close to the value before application of the neurotoxin (Fig. 6A, green). Analogous changes of periods due to coupling changes have been observed in previous network simulations (Gonze et al., 2005; Gu et al., 2009, 2016b). Similarly, representative phase distributions show a narrow, broad, and intermediate spread before applying, during the application of, and after washing out TTX, respectively; see Figure 6B. Lastly, the single-cell oscillations in bioluminescence intensities experience, on average, an approximately 12-fold decline in their amplitudes during application of TTX, which is again partly recovered by the washing procedure; see Figure 6C.

In summary, we observe an increase of the empirical SD in the distributions of single-cell periods and phases as well as a decreasing mean of the amplitude distribution with application of TTX compared to the control condition or after washing out TTX. In analogy to our simulation results and previous studies (Gu et al., 2016a; Beersma et al., 2017), it is therefore plausible to assume that the application of TTX corresponds to a dynamical regime with relatively weak intercellular couplings compared to the other conditions. Note that application of TTX also might affect single-cell oscillators. Such additional effects of TTX could contribute to the period and amplitude changes.

Quite similar results were obtained in the 4 other SCN slice preparations, as shown in Supplementary Figures S9–S12. A summary of the statistical properties in the distributions of periods, phases, and amplitudes among all 5 SCN slice preparations can be found in Supplementary Tables S1–S3, respectively.

Network Model Based on Fitted Parameters

Inspired by our previous results, we next sought to investigate a network model, based on amplitude-phase oscillators that were fitted to the experimental single-cell data. To this end, we fit a “linearized Hopf oscillator” with 5 parameters to the time-series data separately under all 3 conditions, using an autocorrelation approach as performed in Westermark et al. (2009) and Tokuda et al. (2015); see the Materials and Methods section for details on the model and its corresponding autocorrelation function. Because the fit of such an autocorrelation function is an ill-conditioned problem due to the summation of 2 exponentials (see Equation (6)), we set $γ = 0.1 h^{- 1}$ throughout all fits, which is tantamount to the value used for the simulations shown in Figures 1 –4. It can be noted that the periods and amplitudes obtained by this method are highly correlated with those obtained by fitting the linearly damped model, although amplitudes are slightly but systematically estimated higher; see Supplementary Figure S13.

Our fitting procedure has been applied to the time-series data obtained under the application of TTX. Thus, we can investigate an oscillator network closely related to experimental data. We observe a qualitatively similar response of this biologically motivated network to changes in the coupling strength $K$ in comparison to mean-field coupled Poincaré oscillators. The variances of period and phase distributions decrease while amplitudes expand with increasing coupling strength K; see Figure 7. The amplitude expansion with increasing coupling strength $K$ is model dependent (compare Fig. 4C, Fig. 7C, Suppl. Fig. S14C, and Suppl. Fig. S15C), which is confirmed by our analytical approximations; see Supplementary Text S1.2.

Finally, we investigate the behavior of oscillator networks under different random network topologies. Here, we assume that the probability that a certain oscillator has an influence on any other oscillator is given by a small number p; see the Materials and Methods section for details. In summary, we find that the results depicted in Figures 4 and 7 are robust for a broad range of network topologies as long as the connectivity between the oscillators is high enough. This holds true for coupled Poincaré oscillators (Suppl. Fig. S16) as well as coupled linearized Hopf oscillators that are based on fitted parameters (Suppl. Fig. S17).

Discussion

The circadian clock can be regarded as a system of coupled oscillators on the levels of cells, peripheral organs, and organisms. In many cases, for example the wild-type SCN, the individual oscillators are perfectly synchronized. In the case of undercritical coupling (e.g., for dispersed fibroblasts), no synchronization is observed. For varying conditions such as cell densities (Noguchi et al., 2013), population size (Bloch et al., 2013), knockouts (Ono et al., 2013), developmental stages (Tokuda et al., 2015), inhibitor applications (Herzog et al., 2017), or entrainment schedules (Buijink et al., 2016), changes of coupling strength can be expected. As long as no transitions between coherent and incoherent states are induced, however, no obvious qualitative changes occur.

A direct measurement of the absolute coupling strength is currently unavailable. This would be best estimated by monitoring the response of a single oscillatory entity with respect to its network embedding after an isolated perturbation of this entity. Here, we emphasize to indirectly quantify the relative coupling strength between different experimental conditions by analyzing distributions of periods, phases, and amplitudes. Network simulations show that for undercritical coupling, period distributions become narrower due to partial synchronization and frequency pulling. For overcritical coupling, phase distributions are informative because large coupling induces narrow phase distributions. This effect resembles reduced spreads of chronotypes for strong zeitgebers (Roenneberg et al., 2003; Granada et al., 2013). Furthermore, amplitudes might increase with coupling due to resonance effects. We provide analytical calculations for the amplitude expansion that agrees well with our numerical simulations.

We confirmed the theoretical predictions from oscillator theory using data from recent experiments with SCN slices (Abel et al., 2016). The untreated or restored SCNs exhibit larger amplitudes and narrow period and phase distributions compared to data with TTX treatment. According to our theoretical results, this is indicative of a reduction of coupling with application of TTX, which is in accordance with previous interpretations (Yamaguchi et al., 2003; Abel et al., 2016).

There are numerous additional applications of coupling strength quantifications in SCN preparations. For example, it has been suggested that coupling is changing in early development (Ono et al., 2016). Moreover, the relative strength of coupling via gap junctions and coupling factors such as vasoactive intestinal peptide (VIP), γ-aminobutyric acid (GABA), and arginine vasopressin (AVP) is discussed (Herzog et al., 2017). Our study encourages closer inspections of the corresponding period, phase, and amplitude distributions in forthcoming studies. Importantly, all these properties continuously change in a dose-dependent fashion. Often, a binary classification (synchronized or not synchronized) is used. Our approach allows to study dose-dependent changes of coupling strength.

Furthermore, circadian clocks have been found in virtually all mammalian tissues (Dibner et al., 2010). Thus, the quantification of coupling between cells in peripheral organs has been raised (Abraham et al., 2010; Guenthner et al., 2014). There are indications that hepatocytes are coupled, leading to a rather strong liver clock (Saini et al., 2013; Guenthner et al., 2014). In dispersed fibroblasts, mathematical modeling predicts coupling, even though no full synchronization is achieved (Rougemont and Naef, 2007). Indeed, Noguchi et al. (2013) have shown that increasing cell densities improve rhythmicity accompanied by increasing amplitudes. According to our concept, these observations indicate coupling.

In plant circadian systems, early studies suggested weak intercellular couplings of cell-autonomous clocks (Thain et al., 2000), but recent data show that circadian coupling allows synchrony among shoot apex clocks as quantified by phase variances (Takahashi et al., 2015). Furthermore, the statistical analysis of individual periods and phases indicates coupling in social groups of rodents (Paul et al., 2014, 2015) and honeybees (Fuchikawa et al., 2016).

In addition to the method described in this article, there are also alternative indirect methods to quantify coupling. For example, monitoring relaxation to equilibria after a perturbation provides information on coupling strength. In the SCN, local perturbations might be induced by optogenetic tools (Jones et al., 2015). In flies, short light pulses induced transient desynchrony, and the relaxation of the neuronal network revealed the kinetics of retuning (Roberts et al., 2015).

Another strategy is to study the kinetics of SCN plasticity on different protocols of light–dark cycles. Regional phase differences within the SCN network can be induced by very long photoperiods (e.g., 20 h of light and 4 h of darkness) (Evans et al., 2013), or by entrainment with unnaturally short and long zeitgeber periods such as 22 h or 26 h (Azzi et al., 2017). The relaxation dynamics of the induced phase differences can be exploited to explore coupling mechanisms.

Mathematically, a network of N uncoupled oscillators with different individual frequencies represents an N-dimensional torus (Grebogi et al., 1985). Contrarily, complete synchronization constitutes a 1D limit cycle in phase space. Thus, the attractor dimension (Farmer et al., 1983) could be regarded as a measure of coupling strength. Estimations of attractor dimensions, however, require long, stationary time series rarely available in chronobiology. Still, the observations of multiple period clusters (Ono et al., 2016), of “relative coordination” (Holst, 1939), and of “splitting” (Gwinner, 1974; de la Iglesia et al., 2000) point to loss of synchrony due to reduced coupling.

Our simulations in Figures 1 –4 were based on amplitude-phase oscillators with global coupling. Our data-based networks with global (Figure 7) or random coupling (Suppl. Fig. S16–S17) do not address spatial structures. For more local coupling schemes or spatial clusters of differential oscillator properties, complex spatiotemporal patterns of periods, phases, and amplitudes can be expected (Bernard et al., 2007; Myung et al., 2012; DeWoskin et al., 2015). In such situations, measures of spatial coherence such as mutual information (Abel et al., 2016) or Moran’s I (Schmal et al., 2017) can complement the analysis of overall distributions discussed in this article.

Our approach assumes that periods, phases, and amplitudes can be measured for individual oscillators. In many cases, only overall signals are available for tissue slices or cell cultures (Maier et al., 2009). In such cases, only indirect cues provide information on coupling strength. It has been shown theoretically and experimentally (Abraham et al., 2010) that coupling can establish strong oscillators (Bordyugov et al., 2015), that is, robust oscillations with limited sensitivity to perturbations and narrow entrainment ranges. Indeed, in vertebrates with strongly coupled pacemaker cells, the entrainment ranges are much narrower than in invertebrates (Aschoff and Pohl, 1978). A narrow entrainment range implies also a high sensitivity of the entrainment phase with respect to period variations (Hoffmann, 1969; Schmal et al., 2015), which is known also for human chronotypes (Granada et al., 2013; Wright et al., 2005).

In summary, we stress that a quantification of coupling strength is essential to understand emergent properties of circadian systems. Oscillator theory predicts that period, phase, and amplitude distributions are useful indicators of coupling strength. We analyze several oscillator networks and find that variances of periods and phases shrink monotonically with increasing coupling strength. The increase of amplitudes appears to be model dependent as confirmed by numerical studies and analytical calculations. The theoretically predicted effects are illustrated by experimental data from SCN slices treated with TTX. We emphasize that our approach is also applicable to peripheral tissues and social networks.

Footnotes

Acknowledgements

The authors gratefully thank Dr Bharath Ananthasubramaniam, Dr Grigory Bordyugov, Dr Jihwan Myung, and Dr Isao T. Tokuda for stimulating discussions, as well as Dr Daniel Granados-Fuentes and Thomas Wang for experimental support.

This work has been supported by the Deutsche Forschungsgemeinschaft (grant no. BO 3612/2-1) to CS, the Bundesministerium für Bildung und Forschung (grant no. 01GQ1503) to HH, and the National Institute of Neurological Disorders and Stroke (grant no. 095367-01) to EH. CS acknowledges support from the Joachim Herz Stiftung.

Author Contributions

CS, EH, and HH conceived and designed the study. CS performed the simulations. EH contributed the experimental data. CS, EH, and HH wrote the article. All authors gave final approval for publication.

Conflict of Interest Statement

The author(s) have no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Notes

References

Abel

Meeker

Granados-Fuentes

John

PCS

Wang

Bales

Doyle

Herzog

Petzold

(2016) Functional network inference of the suprachiasmatic nucleus. Proc Natl Acad Sci U S A 113(16):4512–4517.

Abraham

Granada

Westermark

Heine

Kramer

Herzel

(2010) Coupling governs entrainment range of circadian clocks. Mol Syst Biol 6:438.

Aschoff

Pohl

(1978) Phase relations between a circadian rhythm and its zeitgeber within the range of entrainment. Die Naturwissenschaften 65(2):80–84.

Aschoff

Wever

(1966) Circadian period and phase-angle difference in chaffinches (Fringilla coelebs L.). Comparative Biochemistry and Physiology 18(2):397–404.

Ashwin

Coombes

Nicks

(2016) Mathematical frameworks for oscillatory network dynamics in neuroscience. J Math Neurosci 6:2.

Azzi

Evans

Leise

Myung

Takumi

Davidson

Brown

(2017) Network dynamics mediate circadian clock plasticity. Neuron 93(2):441–450.

Baba

Ono

Honma

(2008) A TTX-sensitive local circuit is involved in the expression of PK2 and BDNF circadian rhythms in the mouse suprachiasmatic nucleus. Eur J Neurosci 27(4):909–916.

Beersma

DGM

(2005) Why and how do we model circadian rhythms? J Biol Rhythms 20(4):304–313.

Beersma

DGM

Gargar

Daan

(2017) Plasticity in the period of the circadian pacemaker induced by phase dispersion of its constituent cellular clocks. J Biol Rhythms 32(3):237–245.

10.

Bernard

Gonze

Cajavec

Herzel

Kramer

(2007) Synchronization-induced rhythmicity of circadian oscillators in the suprachiasmatic nucleus. PLoS Comput Biol 3(4):e68.

11.

Bloch

Herzog

Levine

Schwartz

(2013) Socially synchronized circadian oscillators. Proc R Soc Lond B Biol Sci 280(1765):20130035.

12.

Bodenstein

Gosak

Schuster

Marhl

Perc

(2012) Modeling the seasonal adaptation of circadian clocks by changes in the network structure of the suprachiasmatic nucleus. PLoS Comput Biol 8(9):e1002697.

13.

Bordyugov

Abraham

Granada

Rose

Imkeller

Kramer

Herzel

(2015) Tuning the phase of circadian entrainment. J R Soc Interface 12(108):20150282.

14.

Buijink

Almog

Wit

Roethler

Engberink

AHOO

Meijer

Garlaschelli

Rohling

JHT

Michel

(2016) Evidence for weakened intercellular coupling in the mammalian circadian clock under long photoperiod. PLOS ONE 11(12):e0168954.

15.

de la Iglesia

Meyer

Carpino

Schwartz

(2000) Antiphase oscillation of the left and right suprachiasmatic nuclei. Science 290(5492):799–801.

16.

DeWoskin

Myung

Belle

MDC

Piggins

Takumi

Forger

(2015) Distinct roles for GABA across multiple timescales in mammalian circadian timekeeping. Proc Natl Acad Sci U S A 112(29):3911–3919.

17.

Dibner

Schibler

Albrecht

(2010) The mammalian circadian timing system: organization and coordination of central and peripheral clocks. Annu Rev Physiol 72:517–549.

18.

Enoki

Kuroda

Ono

Hasan

Ueda

Honma

(2012) Topological specificity and hierarchical network of the circadian calcium rhythm in the suprachiasmatic nucleus. Proc Natl Acad Sci U S A 109(52):21498–21503.

19.

Ermentrout

(2002) Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Researchers and Students. SIAM 2002, Philadelphia, USA.

20.

Evans

Leise

Castanon-Cervantes

Davidson

(2013) Dynamic interactions mediated by nonredundant signaling mechanisms couple circadian clock neurons. Neuron 80(4):973–983.

21.

Farmer

Ott

Yorke

(1983) The dimension of chaotic attractors. Physica D 7(1):153–180.

22.

Fuchikawa

Eban-Rothschild

Nagari

Shemesh

Bloch

(2016) Potent social synchronization can override photic entrainment of circadian rhythms. Nat Commun 7:11662.

23.

Glass

Mackey

(1988) From Clocks to Chaos: The Rhythms of Life. Princeton, NJ: Princeton University Press.

24.

Gonze

(2011) Modeling circadian clocks: from equations to oscillations. Cent Eur J Biol 6(5):699–711.

25.

Gonze

Bernard

Waltermann

Kramer

Herzel

(2005) Spontaneous synchronization of coupled circadian oscillators. Biophys J 89(1):120–129.

26.

Granada

Bordyugov

Kramer

Herzel

(2013) Human chronotypes from a theoretical perspective. PLoS One 8(3):e59464.

27.

Grebogi

Ott

Yorke

(1985) Attractors on an N-torus: quasiperiodicity versus chaos. Physica D 15(3):354–373.

28.

Liang

Yang

Rohling

JHT

(2016a) Heterogeneity induces rhythms of weakly coupled circadian neurons. Scientific Reports 6:21412.

29.

Rohling

JHT

Liang

Yang

(2016b) Impact of dispersed coupling strength on the free running periods of circadian rhythms. Phys Rev E 93(3):032414.

30.

Wang

Liu

(2009) Free-running period of neurons in the suprachiasmatic nucleus: its dependence on the distribution of neuronal coupling strengths. Phys Rev E 80(3):030904.

31.

Guenthner

Luitje

Pyle

Molyneux

Leise

Harrington

(2014) Circadian rhythms of PER2::LUC in individual primary mouse hepatocytes and cultures. PLoS One 9:e87573.

32.

Gwinner

(1974) Testosterone induces “splitting” of circadian locomotor activity rhythms in birds. Science 185(4145):72–74.

33.

Herzog

Aton

Numano

Sakaki

Tei

(2004) Temporal precision in the mammalian circadian system: a reliable clock from less reliable neurons. J Biol Rhythms 19(1):35–46.

34.

Herzog

Hermanstyne

Smyllie

Hastings

(2017) Regulating the suprachiasmatic nucleus (SCN) circadian clockwork: interplay between cell-autonomous and circuit-level mechanisms. Cold Spring Harb Perspect Biol 9(1):a027706.

35.

Hoffmann

(1969) Zum Einfluß der Zeitgeberstärke auf die Phasenlage der synchronisierten circadianen Periodik. Zeitschrift für vergleichende Physiologie 62(1):93–110.

36.

Holst

(1939) Die relative Koordination. Ergebnisse der Physiologie, biologischen Chemie und experimentellen Pharmakologie 42(1):228–306.

37.

Honma

Nakamura

Shirakawa

Honma

(2004) Diversity in the circadian periods of single neurons of the rat suprachiasmatic nucleus depends on nuclear structure and intrinsic period. Neurosci Lett 358(3):173–176.

38.

Honma

Shirakawa

Katsuno

Namihira

Honma

(1998) Circadian periods of single suprachiasmatic neurons in rats. Neurosci Lett 250(3):157–160.

39.

Hoppensteadt

Keener

(1982) Phase locking of biological clocks. J Math Biol 15(3):339–349.

40.

Jones

Tackenberg

McMahon

(2015) Manipulating circadian clock neuron firing rate resets molecular circadian rhythms and behavior. Nat Neurosci 18(3):373–375.

41.

Kuramoto

(1975) Self-entrainment of a population of coupled non-linear oscillators. In: Araki

(ed.) International Symposium on Mathematical Problems in Theoretical Physics. Lecture Notes in Physics no. 39. Berlin: Springer. pp. 420–422.

42.

Kuramoto

(2003) Chemical Oscillations, Waves, and Turbulence. New York: Dover.

43.

Levine

Funes

Dowse

Hall

(2002) Resetting the circadian clock by social experience in Drosophila melanogaster. Science 298(5600):2010–2012.

44.

Liu

Weaver

Strogatz

Reppert

(1997) Cellular construction of a circadian clock: period determination in the suprachiasmatic nuclei. Cell 91:855–860.

45.

Locke

Westermark

Kramer

Herzel

(2008) Global parameter search reveals design principles of the mammalian circadian clock. BMC Syst Biol 2(1):22.

46.

Maier

Wendt

Vanselow

Wallach

Reischl

Oehmke

Schlosser

Kramer

(2009) A large-scale functional RNAi screen reveals a role for CK2 in the mammalian circadian clock. Genes Dev 23(6):708–718.

47.

Mohawk

Takahashi

(2011) Cell autonomy and synchrony of suprachiasmatic nucleus circadian oscillators. Trends Neurosci 34(7):349–358.

48.

Myung

Hong

Hatanaka

Nakajima

Schutter

Takumi

(2012) Period coding of Bmal1 oscillators in the suprachiasmatic nucleus. J Neurosci 32(26):8900–8918.

49.

Nagoshi

Saini

Bauer

Laroche

Naef

Schibler

(2004) Circadian gene expression in individual fibroblasts: cell-autonomous and self-sustained oscillators pass time to daughter cells. Cell 119(5):693–705.

50.

Noguchi

Wang

Welsh

(2013) Fibroblast PER2 circadian rhythmicity depends on cell density. J Biol Rhythms 28(3):183–192.

51.

O’Neill

Hastings

(2008) Increased coherence of circadian rhythms in mature fibroblast cultures. J Biol Rhythms 23(6):483–488.

52.

Ono

Honma

(2013) Cryptochromes are critical for the development of coherent circadian rhythms in the mouse suprachiasmatic nucleus. Nat Commun 4:1666.

53.

Ono

Honma

(2016) Differential roles of AVP and VIP signaling in the postnatal changes of neural networks for coherent circadian rhythms in the SCN. Sci Adv 2(9):e1600960.

54.

Paul

Indic

Schwartz

(2014) Social forces can impact the circadian clocks of cohabiting hamsters. Proc R Soc Lond B Biol Sci 281(1779):20132535.

55.

Paul

Indic

Schwartz

(2015) Social synchronization of circadian rhythmicity in female mice depends on the number of cohabiting animals. Biol Lett 11(6):20150204.

56.

Pikovsky

Rosenblum

Kurths

(2001) Synchronization: A Universal Concept in Nonlinear Sciences. Cambridge: Cambridge University Press.

57.

Pittendrigh

Daan

(1976) A functional analysis of circadian pacemakers in nocturnal rodents: IV. Entrainment: pacemaker as clock. J Comp Physiol 106(3):291–331.

58.

Ramkisoensing

van Engeldorp Gastelaars

Michel

Deboer

Rohling

Meijer

(2014) Enhanced phase resetting in the synchronized suprachiasmatic nucleus network. J Biol Rhythms 29(1):4–15.

59.

Roberts

Leise

Noguchi

Galschiodt

Houl

Welsh

Holmes

(2015) Light evokes rapid circadian network oscillator desynchrony followed by gradual phase retuning of synchrony. Curr Biol 25(7):858–867.

60.

Roenneberg

Chua

Bernardo

Mendoza

(2008) Modelling biological rhythms. Curr Biol 18(17):826–835.

61.

Roenneberg

Wirz-Justice

Merrow

(2003) Life between clocks: daily temporal patterns of human chronotypes. J Biol Rhythms 18(1):80–90.

62.

Rougemont

Naef

(2007) Dynamical signatures of cellular fluctuations and oscillator stability in peripheral circadian clocks. Mol Syst Biol 3(1).

63.

Saini

Liani

Curie

Gos

Kreppel

Emmenegger

Bonacina

Wolf

Poget

Franken

Schibler

(2013) Real-time recording of circadian liver gene expression in freely moving mice reveals the phase-setting behavior of hepatocyte clocks. Genes Dev 27(13):1526–1536.

64.

Schmal

Leloup

Gonze

(2014) Modeling and simulating the Arabidopsis thaliana circadian clock using XPP-AUTO. In: Staiger

(ed.) Plant Circadian Networks. Methods in Molecular Biology no. 1158. New York: Springer. pp. 337–358.

65.

Schmal

Myung

Herzel

Bordyugov

(2015) A theoretical study on seasonality. Frontiers in Neurology 6:94.

66.

Schmal

Myung

Herzel

Bordyugov

(2017) Moran’s I quantifies spatio-temporal pattern formation in neural imaging data. Bioinformatics 33(19):3072–3079.

67.

Sonnenschein

Schimansky-Geier

(2013) Approximate solution to the stochastic Kuramoto model. Phys Rev E 88(5):052111.

68.

St John

Doyle

(2015) Quantifying stochastic noise in cultured circadian reporter cells. PLoS Comput Biol 11(11).

69.

Strogatz

(2000) From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators. Physica D 143(1-4):1–20.

70.

Strogatz

(2003) SYNC: The Emerging Science of Sponteneous Order. New York: Hyperion.

71.

Strogatz

Mirollo

(1991) Stability of incoherence in a population of coupled oscillators. J Stat Phys 63(3–4):613–635.

72.

Takahashi

Hirata

Aihara

Mas

(2015) A hierarchical multi-oscillator network orchestrates the Arabidopsis circadian system. Cell 163(1):148–159.

73.

Taylor

Wang

Granados-Fuentes

Herzog

(2017) Resynchronization dynamics reveal that the ventral entrains the dorsal suprachiasmatic nucleus. J Biol Rhythms 32(1):35–47.

74.

Thain

Hall

Millar

(2000) Functional independence of circadian clocks that regulate plant gene expression. Curr Biol 10(16):951–956.

75.

Tokuda

Ono

Ananthasubramaniam

Honma

Herzel

(2015) Coupling controls the synchrony of clock cells in development and knockouts. Biophys J 109(10):2159–2170.

76.

Ullner

Buceta

Díez-Noguera

García-Ojalvo

(2009) Noise-induced coherence in multicellular circadian clocks. Biophys J 96(9):3573–3581.

77.

Webb

Angelo

Huettner

Herzog

(2009) Intrinsic, nondeterministic circadian rhythm generation in identified mammalian neurons. Proc Natl Acad Sci U S A 106(38):16493–16498.

78.

Webb

Taylor

Thoroughman

Doyle

Herzog

(2012) Weakly circadian cells improve resynchrony. PLoS Comput Biol 8(11):e1002787.

79.

Welsh

Logothetis

Meister

Reppert

(1995) Individual neurons dissociated from rat suprachiasmatic nucleus express independently phased circadian firing rhythms. Neuron 14(4):697–706.

80.

Welsh

Takahashi

Kay

(2010) Suprachiasmatic nucleus: cell autonomy and network properties. Annu Rev Physiol 72(1):551–577.

81.

Westermark

Welsh

Okamura

Herzel

(2009) Quantification of circadian rhythms in single cells. PLoS Comput Biol 5(11):e1000580.

82.

Wiener

(1961) Cybernetics: Or Control and Communication in the Animal and the Machine, 2nd ed. Cambridge, MA: MIT Press.

83.

Winfree

(1990) The Geometry of Biological Time. New York: Springer-Verlag.

84.

Woelders

Beersma

DGM

Gordijn

MCM

Hut

Wams

(2017) Daily light exposure patterns reveal phase and period of the human circadian clock. J Biol Rhythms 32(3):274–286.

85.

Wright

Gronfier

Duffy

Czeisler

(2005) Intrinsic period and light intensity determine the phase relationship between melatonin and sleep in humans. J Biol Rhythms 20(2):168–177.

86.

Yamaguchi

Isejima

Matsuo

Okura

Yagita

Kobayashi

Okamura

(2003) Synchronization of cellular clocks in the suprachiasmatic nucleus. Science 302(5649):1408–1412.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

1.33 MB