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.

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\mathrm{ }\mathrm{h}$ at ∆ $t=0.25\mathrm{ }\mathrm{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 ${\mathrm{\theta }}_i(t=0)$ has been drawn from uniformly distributed random numbers between $[0,2\pi ]$ . The initial condition was chosen once and has then been fixed throughout all simulations.

Analysis of Simulated Time Series

The instantaneous phase $\mathrm{\theta }(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\mathrm{ }\mathscr{H}(s(t))$ , where $\mathscr{H}(s(t))$ denotes the Hilbert transform. Along these lines, the instantaneous phase and amplitude are defined by $\mathrm{\theta }(t):=\mathrm{atan2}\left(\Im (z(t)),\Re (z(t))\right)$ and $A(t):=\sqrt{\Re {(z(t))}^2+\Im {(z(t))}^2}$ , respectively, where $\Re (z(t))$ and $\Im (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 \mathcal{T}=[150\times 24\mathrm{ }\mathrm{h}\mathrm{,250}\times 24\mathrm{ }\mathrm{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{\mathrm{d}{\mathrm{\theta }}_i(t)}{\mathrm{d}t}={\mathrm{\omega }}_i(t)=\frac{2\pi }{{\tau }_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 $\mathcal{T}$ via $〈{\tau }_i^{\mathrm{HT}}〉=\frac{2\pi }{a_{\mathrm{\theta }}}$ , where $a_{\mathrm{\theta }}$ is the slope of a linear fit to the (unwrapped) instantaneous phase ${\mathrm{\theta }}_i(t)$ in the interval $\mathcal{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 $\mathcal{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 $\lambda =0.05{\left(\frac{24\mathrm{ }\mathrm{h}}{\mathrm{\Delta }t}\right)}^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{\mathrm{d}x}{\mathrm{d}t}=-\gamma x-\frac{2\mathrm{ }\pi }{\tau }y+{\mathrm{\xi }}_x$$

$$\frac{\mathrm{d}y}{\mathrm{d}t}=-\gamma y+\frac{2\mathrm{ }\pi }{\tau }x+{\mathrm{\xi }}_y$$

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

$$C(t)=\frac{D}{\gamma }\mathrm{ }{e}^{-\gamma t}\mathrm{ }\cos \left(\frac{2\mathrm{ }\pi }{\tau }\mathrm{ }t\right)\mathrm{ }\mathrm{,}$$

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 $\gamma$ , $\tau$ , 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/\gamma }$ ; 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{\mathrm{d}r}{\mathrm{d}t}=\gamma \mathrm{ }(A-r)+{\mathrm{\xi }}_r$$

$$\frac{\mathrm{d}\mathrm{\theta }}{\mathrm{d}t}=\frac{2\pi }{\tau }+{\mathrm{\xi }}_{\mathrm{\theta }}\mathrm{,}$$

with $〈{\mathrm{\xi }}_i(t){\mathrm{\xi }}_i(s)〉=2D_i\delta (t-s)$ for $j\in \left\{r\mathrm{,}\mathrm{\theta }\right\}$ , 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}\left(A^2+\frac{D_r}{\gamma }\mathrm{ }{e}^{-\gamma t}\right)\mathrm{ }{e}^{-D_{\mathrm{\theta }}t}\cos \left(\frac{2\mathrm{ }\pi }{\tau }\mathrm{ }t\right)\mathrm{ }$$

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 $\gamma =0.1\mathrm{ }{\mathrm{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{\mathrm{d}{x}_i}{\mathrm{d}t}={\gamma }_i\mathrm{ }\frac{x_i}{r_i}(A_i-r_i)-\frac{2\mathrm{ }\pi }{{\tau }_i}{y}_i+M_i$$

$$\frac{\mathrm{d}{y}_i}{\mathrm{d}t}={\gamma }_i\mathrm{ }\frac{y_i}{r_i}(A_i-r_i)+\frac{2\mathrm{ }\pi }{{\tau }_i}{x}_i$$

with $r_i=\sqrt{x_i^2+y_i^2}$ and $M=\frac{K}{N}{\Sigma }_{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^{\mathrm{R}}=k_{i\mathrm{,}\mathrm{ }\mathrm{in}}^{-1}{\displaystyle {\sum }_{j=1}^N}{K}_{ij}{x}_j(t)$ , where $K_{ij}$ is assumed to be $K$ if oscillator $j$ influences oscillator $i$ , and $0$ otherwise. $k_{i\mathrm{,}\mathrm{ }\mathrm{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{\mathrm{d}{r}_i}{\mathrm{d}t}={\gamma }_i\mathrm{ }{r}_i\mathrm{ }(A_i-r_i)$$

$$\frac{\mathrm{d}{\mathrm{\theta }}_i}{\mathrm{d}t}=\frac{2\pi }{{\tau }_i}\mathrm{,}$$

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 ${\mathrm{\theta }}_i(t)$ ; see Figure S1 for an illustration. The parameters ${\tau }_i$ , $A_i$ , and ${\gamma }_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{\mathrm{d}{x}_i}{\mathrm{d}t}={\gamma }_i\mathrm{ }{x}_i(A_i-r_i)-\frac{2\mathrm{ }\pi }{{\tau }_i}{y}_i+M_i$$

$$\frac{\mathrm{d}{y}_i}{\mathrm{d}t}={\gamma }_i\mathrm{ }{y}_i(A_i-r_i)+\frac{2\mathrm{ }\pi }{{\tau }_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}$ ${\Sigma }_{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 ${\gamma }_i=0.1\mathrm{ }{\mathrm{h}}^{-1}$ , as suggested by Westermark et al. (2009), Ramkisoensing et al. (2014), and Tokuda et al. (2015). The intrinsic free-running periods ${\tau }_i$ are chosen from a normal distribution with mean ${\mu }_{\tau }=24\mathrm{h}$ and a standard deviation (SD) of ${\sigma }_{\tau }=2\mathrm{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).

Figure 1.

Temporal order formation after coupling. Steady-state dynamics (i.e., after decay of transients) are plotted for a total length of 2 days and 3 different coupling strengths: (A) $K=0.04$ , (B) $K=0.07$ , and (C) $K=0.1$ . Numerical solutions are color-coded with respect to their intrinsic free-running periods ${\tau }_i$ . Blue colors indicate short intrinsic periods, and red colors indicate long intrinsic periods, leading to early and late phases in the case of synchronization, respectively.

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 ${\tau }_i$ tend to phase-lead with respect to oscillators having larger values of ${\tau }_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 $〈{\tau }_i^{\mathrm{HT}}〉$ is estimated for each oscillator by fitting a straight line to the (unwrapped) phase ${\mathrm{\theta }}_i^{\mathrm{HT}}(t)$ and defining the inverse slope of this line times $2\pi$ as $〈{\tau }_i^{\mathrm{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 ${\tau }_i$ in the uncoupled state and the corresponding average periods $〈{\tau }_i^{\mathrm{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 $〈{\tau }_i^{\mathrm{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 $(〈{\tau }_i^{\mathrm{HT}}〉\mathrm{,}\mathrm{ }{\tau }_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 $〈{\tau }_i^{\mathrm{HT}}〉$ decreases with increasing coupling; see Figure 3A.

Figure 2.

Increasing coupling leads to the formation of frequency plateaus. Average periods $〈{\tau }_i^{\mathrm{HT}}〉$ , determined by means of a Hilbert transformation as described in the Materials and Methods section, are plotted against their corresponding free-running periods ${\tau }_i$ for 3 different coupling strengths: $K=0.04$ (blue dots), $K=0.07$ (green dots), and $K=0.1$ (red dots).

Figure 3.

Distributions of oscillator properties. Histograms of the periods (A), phases (B), and amplitudes (C) of the individual oscillators for 3 different coupling strengths. Color coding represents the coupling strength: $K=0.04$ (blue), $K=0.07$ (green), and $K=0.1$ (red). Phase distributions are plotted for a time at the midpoint of the integration interval to reduce edge effects of the Hilbert transform (Pikovsky et al., 2001). All phase distributions have been centered such that the mean of the corresponding distribution lies at point zero.

As noted in the previous section, we can observe a re-organization of the oscillation phases with varying coupling strength. Although (instantaneous) oscillator phases ${\mathrm{\theta }}_i^{\mathrm{HT}}$ are equally distributed among their whole co-domain $-\pi \le {\mathrm{\theta }}_i^{\mathrm{HT}}<\pi$ 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^{\mathrm{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\lesssim 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 ${\sigma }_{\tau }=2\mathrm{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\gtrsim 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.

Figure 4.

Dependence of emergent network properties on the coupling strength K. (A) Standard deviations (SDs) of the period distribution of individual oscillators in dependence of the coupling strength K. (B) Circular variance (gray) of the phase distribution and global phase coherence $R(t)$ in dependence of K. Variances are estimated by fitting a von Mises distribution to the phase values. (C) Mean value (bold line) and SD (gray-shaded area) of the amplitude distribution in dependence of K. The dashed black line denotes the analytical approximation of the mean value as given by Equation (13).

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

$${\mu }_A(K)=A+\frac{K}{4\gamma }\left(1+\sqrt{1-{\left(\frac{2\mathrm{\Delta }\mathrm{\omega }}{K}\right)}^2}\right)$$

in the case of overcritical coupling (i.e., for all K ≥ $2\mathrm{\Delta }\mathrm{\omega }$ ); see Supplementary Text 1.2 for a detailed derivation. Here, ${\mathrm{\mu }}_A(K)$ denotes the mean of the distribution of individual oscillator amplitudes in the network for a given $K\ge 2\mathrm{\Delta }\mathrm{\omega }$ , whereas $A$ and $\gamma$ are the intrinsic amplitude and radial relaxation rate of the individual oscillators. Finally, $\mathrm{\Delta }\mathrm{\omega }$ equals $\frac{2\pi }{\tau +{\sigma }_{\tau }}-\frac{2\pi }{\tau -{\sigma }_{\tau }}$ , where $\tau$ and ${\sigma }_{\tau }$ 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 $\gamma$ as well as the spread of the free-running periods ${\sigma }_{\tau }$ ). 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 ${\sigma }_{\tau }$ 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, ${\sigma }_{\tau }$ has little effect on the amplitude expansion for high coupling strength $K\gg K_c$ , because Equation (13) converges to ${\mu }_A(K)\approx A+\frac{K}{2\gamma }$ for $K\gg 2\mathrm{\Delta }\mathrm{\omega }$ ; see also Supplementary Figure S5B.

For fixed values of $A_i=A$ and ${\sigma }_{\tau }=2\mathrm{h}$ , an increasing radial relaxation rate ${\gamma }_i=\gamma$ 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 ${\gamma }_i=1\mathrm{ }{\mathrm{h}}^{-1}$ ) exhibit almost no variations in amplitude with coupling, whreeas nonrigid oscillators (e.g., for ${\gamma }_i=0.01\mathrm{ }{\mathrm{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.

Figure 5.

Fitting oscillator models to experimental data. (A) Ensemble averages (black line) and standard deviations (gray-shaded area) of bioluminescence time series from automatically identified and tracked neurons in coronal slices of the SCN, as described in Abel et al. (2016). (B–D) Fits to time series from representative cells under all 3 conditions: before applying tetrodotoxin (TTX) (B), during the application of TTX (C), and after washing TTX from the culture medium (D). Black curves denote the original data, detrended by means of a Hodrick–Prescott filter as described in the Materials and Methods section. Blue curves denote simulations of the stochastic ordinary differential Equations (1)–(2) for the parameters, obtained by fits to the corresponding experimental time series as described in the Materials and Methods section. Initial conditions have been selected such that the stochastic dynamics resemble the experimental time series.

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 ${\mu }_{\tau }\approx 26.18\mathrm{h}$ and an empirical SD of $s_{\tau }\approx 0.29\mathrm{h}$ (Fig. 6A, red). During the application of TTX, the ensemble mean of the period drops to a value of ${\mu }_{\tau }\approx 24.91\mathrm{h}$ , accompanied by a broadening of its distribution to $s_{\tau }\approx 2.10\mathrm{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 ${\mu }_{\tau }\approx 25.46\mathrm{h}$ , whereas its empirical SD decreases to $s_{\tau }\approx 0.35\mathrm{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.

Figure 6.

Emergent properties in SCN slice data. Period (A), phase (B), and amplitude (C) distributions under 3 experimental conditions, that is, before applying tetrodotoxin (TTX) (red bars), during the application of TTX (blue bars), and after washing out TTX from the medium (green bars). Periods and amplitudes of single neuronal biolouminscence time series are determined by means of oscillator fits as described in the Materials and Methods section. Phases have been obtained at times $t_1=54\mathrm{ }\mathrm{h}$ (control), $t_2=180\mathrm{ }\mathrm{h}$ (+TTX), and $t_3=365\mathrm{ }\mathrm{h}$ (−TTX) by means of an analytic signal approach, equivalent to what has been used in the case of simulated time series.

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 $\gamma =0.1\mathrm{ }{\mathrm{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.

Figure 7.

Dependence of emergent network properties on the coupling strength $K$ in a network model based on fitted parameters. Same analysis as in Figure 4, using a network model that is more closely related to experimental data. $N=383$ mean-field coupled linearized Hopf oscillators are parametrized by amplitude and period values that are obtained from experimental time-series data fits.

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.