The choroid plexus is an important circadian clock component

We screened Bmal1-ELuc expression in brain loci by imaging enhanced bioluminescence from the Emerald Luciferase¹⁶ in acute cultures of sagittal sections. The circadian clock molecule, BMAL1, is a positive regulator of the core-clock elements, PER and CRY, which form a transcriptional feedback loop at the circadian timescale. Therefore, BMAL1 can be considered the substrate for the molecular circadian clock and its expression indicates the presence of a clock. Bioluminescence reporter-based circadian clock screening has been performed previously^17,18, and the clock locations found in these studies agree with findings from systematic quantification of RNA expression¹⁹. Several loci were identified for putative clock hosts, including the mediobasal hypothalamus²⁰, areas surrounding the 3V and the brain stem²¹ (Fig. 1a; see Supplementary Fig. 1 for full scans). Putative clock loci were then dissected and cultured as isolated explants from animals expressing PER2::LUC, a bioluminescence reporter based on a fusion protein of PER2 and firefly luciferase¹⁸. A strong bioluminescence signal was detected in the cerebellum²², but continued circadian rhythmicity of clock gene expression could not be confirmed in this study due to the difficulty of maintaining cerebellar tissue in culture. Weak circadian clock activities were observed in subsets of cells in the hippocampus and the cortex, but found insufficient for consistent characterization. Strong, persistent circadian expression of PER2::LUC was found mostly in the CVOs, including the CP, the pituitary gland, the median eminence (ME), the pineal gland, and the area postrema (AP), all of which are in direct contact with the CSF (Fig. 1b, c). We ranked the strength of the circadian clock by resolving the circadian oscillatory component and comparing it with the noise component in the time series (see Methods and Supplementary Fig. 2). This indicates that the strongest oscillation occurs in the CP, followed by the pituitary gland and the SCN (Fig. 1d). The robust circadian rhythmicity of the CP stands out among all explants sampled and in all CP structures sampled from the third ventricle (3V), the lateral ventricle (LV), and the fourth ventricle (4V), despite their heterogeneous molecular identities²³. Notably, the robustness of the CP oscillation exceeds that of the master circadian clock in the SCN.

The 4V CP is particularly valuable as a model system because despite its large size, it is possible to culture and observe the whole structure. The CP can be dissected in its intact form and monitored through a cooled-CCD camera. By tracking circadian expression of PER2::LUC and Bmal1-ELuc at the single-cell level, we found that cellular circadian clocks in the CP maintain much better phase synchrony than in the SCN (Fig. 2a, b). Phases within the SCN have been shown to encode daylength and other second-order characteristics of day–night cycles²⁴. Therefore, strong complete synchronization is not necessary in the neuronal network of the SCN. However, in the non-neuronal tissue network of the CP, component cellular clocks are highly synchronized. In the networks of astrocytes, endothelial cells, or endocrine cells, gap junctions mediate local coupling and create population-level dynamics, such as calcium waves²⁵. Expression of gap junctions has been shown by high-throughput analysis in the CP of the mouse²⁶ and the rat²⁷. Hence, a likely mediator of this synchronization is gap junction coupling, and we find strong expression of connexin 43 (Cx43) in the CP. Gap junctions connect neighboring cells, indicated by the hexagonal expression around cells (Fig. 2c). The CP epithelial cells, driven developmentally by FOXJ1 transcription factor, are known to express Cx43²⁸.

Mathematically, gap junction connections can be described by nearest neighbor coupling between single-cell circadian oscillators. Cells on the boundary end-up getting less coupling input, while cells in the central area of the tissue receive concentrated coupling input. Coupling can simultaneously affect the period and the amplitude of circadian oscillation²⁹ due to nonlinear interaction of transcriptional circuits³⁰. A straightforward way to model the interdependence between the frequency (inverse of period) of oscillation and its amplitude is by generic oscillator models^31–33. One of the simplest conceptual models for circadian oscillators is the modified Poincaré oscillator^34,35. In the simplest case, the interdependence between the oscillation frequency and amplitude is assumed to be linear (see Methods). A radially outward deviation of the oscillator from its steady-state orbit causes a proportional change in period (or increase in angular velocity). This dependence relation has been referred to as “twist” or “shear” in the literature, and by convention, a negative twist parameter (ε) can describe an oscillator in which an amplitude increase is accompanied by a decreasing period³⁶ (Fig. 2d).

At the whole-tissue level, we observed spatial patterning of phases, with the leading phase appearing in the central area of the tissue (Fig. 2e, upper; Moran’s I ~ 0.5 for initial 3–4 days, p< 0.05, Supplementary Fig. 4; see also Supplementary Figs. 6–8 F and G)³⁷. We can explain the phase advance of the central area of the CP by the accumulation of coupling effects from the boundary (Fig. 2e, lower left). Due to resonance effects, amplitudes expand more in the central part of the chain. In a system where each oscillator has a negative twist amplitude–period relationship, the central part exhibits the shortest period. It has been previously demonstrated that short periods lead to early phases, while long periods lead to late phases upon synchronization to a mean field or an external signal^33,38. Thus, in the fully synchronized state, the short-period central part phase advances in comparison to lateral parts of the 4V CP. A chain of nearest neighbor-coupled Poincaré oscillators with negative twist (ε< 0) recaptured the observed phase relations (Fig. 2e, lower right; Supplementary Fig. 5) (see Methods for simulation details).

The negative twist represents an inverse relationship between the period and the amplitude of an oscillator. When extended to a system of oscillators, the inverse relation strengthens as a result of an amplitude expansion as the accumulated strength of coupling grows. In the cultured 4V CP, we weakened the coupling strength by increasing concentrations of the gap junction blocker, meclofenamic acid (MFA). MFA is a water-soluble drug that potently and reversibly blocks gap junctions³⁹. We took advantage of the bilaterally symmetric anatomy of 4V CP to overcome tissue-to-tissue variability. When separated into left and right sides in two dishes, both sides independently maintained perfectly identical circadian rhythms (Fig. 3a, 0 µM; n = 9 pairs; p = 0.55 for periods, paired t-test). We, therefore, compared PER2::LUC oscillations between the pairs. As a control, one side of the CP was cultured in a control medium, and the other side of the CP in a MFA-containing medium (Fig. 3a, top inset). This pairwise comparison was the experimental paradigm that we applied throughout (see Methods). As predicted by the twist model (Fig. 2d, e lower, and Supplementary Fig. 5), MFA dose dependently increased the period and decreased the amplitude (50–300 µM; n = 4–13 pairs; all data points except for outliers are indicated in Fig. 3b–d). A grid-based analysis of two-dimensional (2D) time-lapse recordings of PER2::LUC oscillations before application of MFA, during application of MFA, and after washing MFA from the medium, reveals that application of the gap junction blocker reversibly decreased global phase coherence and increased the standard deviation of the single-cell period distribution, indicating a reversible effect of the gap junction blocker (Supplementary Fig. 9). If removal of local coupling by gap junction inhibition shifts the amplitude and period along an inverse relationship, physical dissociation of the tissue is expected to cause a similar shift, albeit to a lesser degree. We cut the test side of the CP in 6–30 equal pieces, with an increment of 6. The amplitude decreased and period increased with the number of pieces, but the trend saturated after 18 pieces. The period continued to increase after enzyme dissociation, although the amplitude decrease could not be compared due to the cell loss during the dissociation process (for all physical and enzyme dissociation, n = 5–9 pairs; Fig. 3e, f). An analogous increase of the ensemble period can be observed in a grid-based single-cell analysis, comparing a surgically dissected 4V CP with a control (Supplementary Fig. 10). These characteristics are recaptured by the one-dimensional (1D) chain of coupled oscillators with negative twist (simulations with Gaussian random samples of periods and initial phases, n = 4; Fig. 3g–i). A 2D model of nearest neighbor coupled Poincaré oscillators, using a realistic geometry (see Methods for simulation details), mimics the experimentally observed phase patterns, as well as the relationship between the coupling strength, the ensemble amplitude, and the ensemble period as indicated by the above described results (Supplementary Fig. 11). It should be noted that spatial gradients of single-cell oscillator properties or coupling strength can lead also to the experimentally observed spatial patterns as discussed in the previous paragraph (Supplementary Fig. 12). Importantly, we were able to reproduce the experimentally observed dependency between coupling, amplitude, period, and phase, solely under the assumption of a negative oscillator twist (compare Fig. 3d, f with Fig. 3g–i, Supplementary Figs. 11 and 12).

The inverse relationship between period and amplitude sets an operating limit for the CP circadian clock. In other words, we can predict that changing the coupling strength will alter the period of circadian gene expression oscillation in the whole tissue by simultaneously altering the amplitude. If this is true, the developmental distribution of period and amplitude in the ensemble of control samples would fall onto the inverse relationship in the same operational range (period of 23–26 h and amplitude of 20–150 kcpm or thousand photon counts per minute, Fig. 3d, f). The 4V CP is especially suitable for this kind of experiment since the entire tissue can be extracted in intact form and its size does not vary greatly compared to the changes in brain size throughout development (i.e., the CP is a relatively large structure in embryonic and neonatal brains). We pooled data from 4V CP and SCN explants from mice at various developmental stages (postnatal days P1–P500). When we superpose PER2::LUC oscillations from ensembles of CP tissues from the same zeitgeber time, we discover that in P0, the oscillations damp quickly and the period is long, often longer than 24 h. Oscillations are not synchronized to the external light–dark cycles, as is evident in the dispersed peaks in the time series (Fig. 4a, top). The damping becomes less severe and the peaks become more coherent at P7 and the trend continues progressively as the age of the animal increases (Fig. 4a). Improvements in robustness (inverse of damping) and coherence, which often indicate stronger coupling, are associated with an ontogenetic decrease of period (Fig. 4b). The critical transition of the circadian period seems to occur at around P10 (characteristic time of 9.47 days as found by exponential fit). However, no such transition could be found in the SCN samples (Fig. 4c). This supports the idea that coupling in the CP strengthens with age and because of the twist in each of the component oscillators, the period-amplitude pair shifts from the lower right to upper left as an animal ages. The period-amplitude scatterplot of all control data is consistent with this hypothesis (R = −0.43 ± 0.06, Pearson’s correlation coefficient, p = 1.65 × 10^–11, bootstrap z-test; Fig. 4d, left). The inverse relationship is not as steep in SCN samples (R = −0.28 ± 0.14, p = 0.023, bootstrap z-test; Fig. 4d, right) or as significant as the CP when the sample size is controlled (CP: R = −0.43 ± 0.13, p = 5.89 × 10^–4, after size-controlled resampling). The intact SCN is anatomically less distinct and more difficult to sample than the CP and the amplitude comparison between the SCN samples may be less reliable than in the CP samples. Nonetheless, the clear twist-mediated developmental decrease of period and increase of robustness in the CP circadian clock cannot be found in the SCN, which attains robust PER2::LUC oscillation at P7⁴⁰. This suggests potential differences in the coupling mechanism between the two systems.

The possibility that the CP circadian clock is flexible within an operational range leads us to speculate on its functional impact on whole-brain circadian rhythms. The SCN was long thought to be the major generator of behavioral and physiological changes in circadian rhythms of the body, until it was found that autonomous circadian clocks exist in multiple tissues¹⁷ and cell types, including fibroblasts⁴¹. Although it is now thought that the SCN is a master coordinator of these clocks, we know that most of the neural connectomes contain feedback connections. We have therefore evaluated a possible interaction between the CP and the SCN in vitro. Tissue-to-tissue interaction in the culture system has been demonstrated to be feasible by placing one tissue directly above the other⁴². We once again used one lateral side of the 4V CP and the SCN as the control, and the other side as the test group. The effect of circadian clock activity in one tissue on the other tissue can be evaluated in a pairwise manner by coculturing non-bioluminescent (LUC (−)) tissue on the bioluminescent (LUC (+)) effector tissue and comparing the bioluminescence oscillation with the contralateral partner without coculture (Fig. 5a, upper inset). The SCN under coculture with LUC (−) CP showed accelerated PER2::LUC oscillation (two representative pairs are shown in Fig. 5a, lower; complete statistics in Fig. 5b, left, ^***p < 0.001, paired t-test, n = 4 pairs). On the other hand, no change in oscillation period was observed in the LUC (+) CP cocultured with LUC (−) SCN (Fig. 5b, right, n = 4 pairs). This suggests that there is an internal entrainment stream from the CP to the SCN. Interestingly, the average period of PER2::LUC oscillation in the SCN, normally longer than 24 h, became “normalized” to the behavioral circadian period of 23.7 ± 0.1 h⁴³ when cocultured with the CP (Fig. 6c). The behavioral period under constant darkness is a consistent phenotype in the C57BL/6J background strain (23.83 ± 0.04 h in ref. ⁴⁴; PER2::LUC/+23.61 ± 0.09 h in ref. ¹⁸).

The observation that the period of PER2::LUC oscillation in cultured SCN is longer than the behavioral freerun period of locomotor activities under constant darkness⁴³ remains unexplained. In freely moving mice, the PER2::LUC oscillation period in the in vivo SCN matches that of the behavioral period of the animal⁴⁵. It is therefore conceivable that the coculture condition in Fig. 5 minimally represents the normal brain environment and that the SCN circadian clock is feedback entrained through the CSF by as-yet-unidentified signaling molecules from the CP. We can experimentally test this hypothesis in vivo by selectively silencing the circadian clock in the CP. FOXJ1 is a transcription factor that induces differentiation of motile ciliated cells. In the brain, its expression is limited to the CP and ependymal cells⁴⁶. By crossing FOXJ1-Cre transgenic mice (Cre^FoxJ1) and floxed Bmal1 mice (Bmal1^flx/flx), we generated mice with targeted knockdown of Bmal1 in the CP (Cre^FoxJ1;Bmal1^flx/flx, Fig. 6a). The mean period of circadian locomotor activities under constant darkness in control Bmal1^flx/flx littermates was identical to that of the wild type (WT: 23.70 ± 0.04 h, mean ± SEM; n = 4, Fig. 6b). As predicted, the period increased progressively with the degree of CP-specific Bmal1 knockdown in heterozygous mice (Het.: 23.81 ± 0.05 h, n = 5) and increased significantly in the homozygous knockdown (Homo.: 23.92 ± 0.04 h, n = 4, p < 0.01, Fig. 6c). This increase is consistent with the coculture results of Fig. 5.

In the brain environment, it is plausible that the CP diffuses its clock signal to the SCN through ventricular circulation of the CSF. The SCN-to-CP signaling pathway is not known, but an array of indirect pathways is feasible, such as the indirect neuronal innervation or blood-borne signaling through the neuroendocrine system (Fig. 6d). To make sure that Bmal1 knockdown was specific to the CP, we measured Bmal1 transcripts at the time of light-off in the same animals tested for their circadian locomotor activities. We found that the Bmal1 transcripts were expressed similarly in the SCN across genotypes (Homo. vs WT: p = 0.07). However, in both the 4V CP and the LV CP, Bmal1 expression was close to the detection limit in Cre^FoxJ1; Bmal1^flx/flx (Homo. vs WT: p < 0.01). We note that when samples from the heterozygous knockdown animals (Cre^FoxJ1; Bmal1^flx/+) were considered together, the locomotor period was highly correlated with the Bmal1 expression level in the CP (Pearson’s correlation coefficient, R = −0.75 (4V CP), R = −0.74 (LV CP), p < 0.005, Fig. 6e).

Mice under C57BL/6J background, of both sexes and ages ranging from postnatal day 0 to 600, were maintained under equinox (12 h light, 12 h dark) conditions and at constant ambient temperature. Heterozygous bioluminescence reporter-expressing strains PER2::LUC (MGI ID: 3040876)¹⁸ and Bmal1-ELuc¹⁶ were used for explant culture experiments. For CP-specific Bmal1 deletion, homozygous Bmal1^flx/flx and heterozygous Bmal1^flx/+; FOXJ1-Cre (MGI ID: 3797103)⁴⁶ were crossed to produce Bmal1^flx/flx littermates with the Cre (test group) or without Cre recombinase (control group). All animal protocols were reviewed and approved by the RIKEN Brain Science Institute.

Adult mice of different genotypes, including control littermates, were individually housed in a light–dark controlled chamber with a ventilating fan that provided constant environmental noise. Locomotor activities were monitored with a passive infrared sensor and logged on the computer as described previously². Circadian periods were estimated from the chi-squared periodogram with a custom code written in Mathematica⁶².

Tissue culture and bioluminescence measurement protocols were identical to those used previously except for the addition of CP sampling⁴³. PER2::LUC and Bmal1-ELuc mice were anesthetized and killed about 0–2 h before lights off and under dim light, brains were isolated to ice-cold HBSS within 15 s after decapitation. The brain was rinsed in fresh ice-cold HBSS and coronally cut into two close to the lambda point. For imaging-based scanning of Bmal1-ELuc expression, sagittal sections were made at 150 µm thickness on a Microslicer (DSK). For time-lapse imaging of the SCN, the anterior portion was coronally sliced at 250 µm thickness and trimmed to sample the mid-posterior section of the SCN. The posterior portion was transferred to a new dish containing ice-cold HBSS for sampling the 4V CP.

For imaging of the CP and tissue explant luminometry, select brain loci were dissected and isolated in ice-cold HBSS with forceps and a scalpel blade under the stereomicroscope. For collecting the 4V CP, the 4V was accessed by lifting the gap between the cerebellum and the brain stem, and the whole CP attached on the side of the cerebellum was gently taken out with forceps. The AP was sampled from the nearby area on the brain stem with a scalpel blade. The pineal gland was isolated from the inner side of the skull after removing surrounding dura mater and blood vessel. The whole pituitary gland was taken from the base of the skull with forceps. The whole SCN was dissected from the ventral side of the brain with a scalpel blade directly under and slightly posterior to optic chiasm. The ME and VLPO were collected along the base of 3V. Consistency of PER2::LUC activities was confirmed by comparing against tissue explants isolated from imaged slices. The anatomical location of the SCN was independently confirmed with GAD67 expression and precision of isolation was further assessed by Avp and Vip transcript expression patterns in dorsoventral subdivisions. Culture of each explant was maintained at 37 °C on the culture membrane (Millicell-CM; Millipore) in a vacuum-sealed 35-mm dish containing 1 mL of DMEM (Sigma) with NeuroBrew-21 supplement (Miltenyi Biotec), penicillin-streptomycin (25 U/mL–25 µg/mL; Nacalai Tesque), and 100 µM beetle luciferin (Promega).

Single-cell resolution imaging was performed on a LV200 (Olympus) with the Orca C4742-80-12AG cooled-CCD camera (Hamamatsu) and a custom-built incubator-microscope system with the Orca R2 cooled-CCD camera (Hamamatsu). Binning was set to 4 × 4 and exposure time was 15 min. For whole-tissue monitoring of bioluminescence, we used a 24 dish, double photomultiplier tube (PMT) system (LM-2400; Hamamatsu) with a sampling interval of 15 min. SCN explants for whole-tissue bioluminescence were attempted in the whole SCN dissected out of the ventral side of the brain.

For gap junction inhibition, a stock solution of MFA sodium salt (Sigma-Aldrich) was prepared at ×1000 concentration in DDW and refrigerated before experiments. MFA was chosen over carbenoxolone (CBX, 400 µM) and mefloquine (MEQ, 100 µM) due to toxicity under our experimental system. The drug effect was pairwise compared in each bilateral side of the tissue (SCN and 4V CP), by culturing one side in the control medium (control) and the other side under manipulation or the drug-containing medium (test). Micro-dissection of the CP was performed on the CP mounted on the culture membrane. Enzyme dissociation of the CP was performed using TrypLE Select ×1 solution (no phenol red, Thermo Fisher Scientific) according to the manufacturer’s protocol. The trypsin replacement dissociation enzyme has been conveniently called “trypsin” or “Tryp.” in the main text and figures.

Animals from the locomotor activity monitoring chamber were anesthetized and decapitated under darkness between 1 h before and after the lights off time. The brain was isolated and immersed in ice-cold HBSS within 15 s after decapitation under dim light. Basic RNA sampling and RT-qPCR methods follow our previous protocol². Whole bilateral pairs of SCN, 4V CP, and LV CP were dissected out in HBSS and each sample was transferred to 50 µl TRK lysis buffer (Omega Bio-Tek) and stored immediately at −80 °C. Total RNA was microcolumn-purified using the ENZA total RNA kit (Omega Bio-Tek) with the final elution volume of 20 µl in DEPC water. cDNA was synthesized from total RNA of 113.4 ± 18.4 ng (SCN), 318.2 ± 3.1 ng (4V CP), and 308.5 ± 5.2 ng (LV CP) (mean ± SEM) using SuperScript II RT (Thermo Fisher Scientific) with random primers. GAPDH was chosen as the internal control for each sample. RT-qPCR was performed in triplicate with SYBR Green DNA binding dye on StepOnePlus (Applied Biosystems). The following primer sequences were used (5′–3′): GAPDH forward ACGGGAAGCTCACTGGCATGG CCTT, GAPDH reverse CATGAGGTCCACCACCCTGTTGCTG; mBmal1 forward GCAGTGCCACTGACTACCAAGA, mBmal1 reverse TCCTGGACATTGCATTGCAT. Additional methods used for Supplementary Fig. 14 are described in Supplementary Methods.

The 4V CP sample was fixed in 4% PFA, permeabilized in 0.2% Triton X-100, blocked in 5% normal goat serum and stained slowly at low temperature (4 °C) as described in our previous work⁴³. As a primary antibody, anti-connexin-43 (Cx43) rabbit polyclonal antibody (Sigma C6219, Lot. 113M4756; RRID: AB 476857) was used at dilution of 1:500–1:1000. The second antibody was anti-rabbit antibody conjugated with Alexa 488 at dilution of 1:500. The stained sample was mounted in VECTASHIELD mounting medium with DAPI (Vector Labs) before examination on an FV1000 confocal microscope (Olympus).

Time series from the PMT (Fig. 1c) were analyzed by means of the multiresolution analysis (MRA) as described previously⁶³, using the PyWavelets software (http://pywavelets.readthedocs.io↗) for the Python Programming language. The MRA uses the discrete wavelet transform to decompose a given time series s(t) into components (also called details) at different disjoint frequency bands⁶⁴. In order to reduce edge effects, the first and last days of data were truncated after application of the MRA for subsequent analysis. The method has been verified on our own dataset from the SCN (Supplementary Fig. 2). Rhythmicity can be quantified by the ratio of the variance explained by the circadian component (D_circadian or D₆ normalized over other components, covering from 16 h to just before 32 h) and the variance explained by the noise component (D_noise or D₁ normalized over other components, covering from 0.5 h to just before 1 h). Strong circadian rhythms are associated with a large corresponding ratio. For the rhythmicity score, the components were normalized across all tissues,

where i = 1, 2, 3, …, 7 indexes each of spectral component and all tissues indicate the 12 tissues as in Fig. 1d.

Statistical significance was tested for spatial patterns in the 1D chain of Bmal1-ELuc oscillation phases, as well as for a simulated 2D grid of nearest neighbor coupled Poincaré oscillators (see Fig. 2e and Supplementary Figs. 11–13, respectively). To this end, the spatial autocorrelation index Moran’s I (MI) under the assumption of a nearest neighbor interaction was applied. We recently demonstrated that MI reliably detects dynamic signatures of spatial order in 2D arrays of circadian gene expression data³⁷. Briefly, the MI takes into account the cyclic nature of phase variables by

where N denotes the number of observables X_i, and w_ij the spatial weights, and dθ(X1,X2):=tan−1sinX1−X2cosX1−X2 assigns a distance between two cyclic variables X₁ and X₂. The mean value X¯ of a set of cyclic variables {Xi}i=1N defined by X¯:=tan−1(S¯∕C¯) with C¯:=(1∕N)∑i=1Ncos(Xi) and S¯:=(1∕N)∑i=1Nsin(Xi). We also used a similar quantity that we called circular MI, I_circ, which is essentially Eq. 2 with X_i= exp(iθ_i) and d(x, y) = x − y. Then,

where. X ¯ = 1 N ∑ j e i θ j

Supplementary Fig. 4A shows MI values as calculated using Eq. 2 for the phase patterns in Fig. 2e (upper). Statistical significance of observed MI values was assessed by a comparison with a sampling distribution under the null hypothesis of no spatial autocorrelation, computationally determined by means of a Monte Carlo approach (see Supplementary Fig. 4B). Following a brief transitory period after transferring the tissue to the culture medium, the phase organization adopts a spatial pattern that is statistically significant for a duration of roughly 3.5 days.

A modified Poincaré oscillator,

conveniently describes conceptually important oscillator properties, such as the radial relaxation rate λ_i, amplitude A_i, twist ε_i, and angular velocity ω_i = 2π/τ_i, where τ_i is the internal free-running period of an oscillator i. Eq. 3a, b describe the dynamical evolution of the radial component r_i(t) and phase component φ_i(t).

Isochrons are a set of initial conditions in the basin of attraction of an attracting limit cycle that have the same asymptotic phase as t → ∞⁶⁵. Importantly, isochrons provide information about the resetting and entrainment behavior of an oscillating system, such as the circadian clock. Isochrons of the oscillating system given by Eq. 3a, b can be obtained analytically.

For a given set of initial conditions r₀ and φ₀ at time t₀ = 0, the solutions of Eq. 3a, b at time t can be read as

and

respectively. From this it follows that any dynamics starting on the limit cycle, i.e., r₀ = A, are given by

and hence

By definition, isochrons at a certain phase φ_0,LC are given by all combinations of (r₀, φ₀) that fulfill the constraint

i.e., φ(t), as well as φ_LC(t) approach the same phase after the decay of transient dynamics as t → ∞. Supplementary Fig. 3A–C depict examples of the vector-field, isochrons, and dependency of the oscillator period on the amplitude for different values of twist ε ∈ {−0.01, 0, 0.01} h⁻¹ and fixed values of A = 1, λ = 0.02 h⁻¹, and τ = 24 h, respectively.

When the cells are coupled to others, coupling was introduced as an additive term on both coordinates in the Cartesian form of Eq. 3a, b as described in our previous work^35,66, using the definitions x_i = r_i cos(φ_i) and y_i = r_isin(φ_i). Then, Eq. 3a, b become

If the system is sufficiently close to synchrony φi≈φj for all (i,j), cos(φj−φi)≈1−0((φj−φi)2) in Eq. 9a, and we can use mean field approximation, where r_i approximately equals a neighborhood average, i.e.​, rjj∈Ni≈ri. Then, Eq. 9a becomes

where Ki= ∑j≠iKij. Then, at a stationary state (dr_i/dt = 0), ri≈Ki∕λ+A and the phase evolution is approximately

where ωi′=ωi−εKi∕λ. Experimental data for amplitudes can constrain the ratio, K_i/λ, from the amplitude increase compared with the case in which gap junctions are blocked, and also ε from the frequency change due to nonzero K_ij. However, λ and K_ij cannot be easily determined, and we chose suitable values that can replicate synchrony and MI (see above and ref. ³⁷) in the data (see also Supplementary Fig. 13 for how synchronization of the system and circular MI evolved when λ and K_ij slowly increased while K_i/λ was held fixed in a simulation).

A 1D chain of N Poincaré oscillators was constructed, each receiving coupling inputs from immediate neighbors only. Coupling was introduced as an additive term on both coordinates in the Cartesian form of Eq. 3a, b as described in our previous work^35,66, using x_i = r_i cos φ_i, y_i = r_i sin φ_i, and ri=xi2+yi2.

The coupling constant K_ij is nonzero, K_ij = K_j, only if j = i ± 1 for i = 1, 2, 3,… N, while K₀ = 0 and K_N1 = 0. Single-cell oscillatory parameters used for the simulations in Figs. 2e and 3g–i and Supplementary Fig. 5 were A_i = 1, λ_i = 0.02 h⁻¹, and ε_i = −0.01 h⁻¹. Dose dependency of coupling inhibition was evaluated by multiples of K_i with q = 0, 1, 2,…10. In Fig. 3g–i, identical simulations were run for four random sets of Gaussian period (µ= 25.5 h, σ = 1 h) and initial phase distributions (σ = 1 h).

In Supplementary Figs. 11, 12 and 13, the effect of twist in spatially extended systems of coupled Poincaré oscillators was analyzed for different geometries, coupling topologies and oscillator properties. In Supplementary Figs. 11 and 12, it was assumed that oscillators couple solely through a 2D nearest neighbor interaction given by a von Neumann neighborhood, i.e., the four most adjacent cells. The coupling terms for dx_i/dt and dy_i/dt in Eq. 12a, b were replaced with 0.5∑j=1NKijxj and 0.5∑j=1NKijyj, respectively, where K_ij = K whenever oscillators i is influenced by oscillator j through nearest neighbor interaction. In Supplementary Fig. 13, each cell coupled to neighbors within a certain distance with the same equation.

In the case of N oscillators, we have 2N degrees of freedom in choosing initial conditions x₀ and y₀ and 4N degrees of freedom in choosing single-cell oscillator properties. In order to find a set of parameter values that qualitatively reproduced dynamical characteristics from the experimental analysis, we fixed the single-cell amplitude to A₀ = 1 throughout all simulations. We restricted ourselves to a nearest neighbor coupling topology as suggested by putative gap-junctional coupling and sampled the period from a normal distribution with mean 25.5 h and a standard deviation of 1 h. A mean value of 25.5 h has been chosen since period values saturate as they approach this value (Fig. 3d, f), i.e., after application of a gap junction blocker or surgical cuts. Subsequently, we determined an appropriate set of initial conditions and manually chose values for ε, λ, and K such that the model qualitatively reproduced the data best, i.e., the experimentally observed phase pattern is well represented and the period in case of highest coupling strength K is close to the experimentally observed value of 23.5 h as indicated by Fig. 3d, f.

Twist chain oscillator simulations were performed on Mathematica (Wolfram Research, IL), the Python programming language using scipy.integrate.odeint function in SCIentific PYthon (https://scipy.org/scipylib/↗), and MATLAB 2016b (Mathwork, VA). Visualization and other data analyses, including periodogram and fast Fourier transform (FFT)-based period analysis, were performed on Mathematica using the custom-made PMTAnalysis package (http://sourceforge.net/projects/imaginganalysis↗)⁴³. All circadian periods of the bioluminescence oscillations were estimated with FFT and all amplitudes were quantified by root mean squares of raw oscillations for equal duration of 10 recording days. Detailed methods for 2D grid-based data analysis are described in the Supplementary Methods. All statistical analyses were done on Mathematica and MATLAB 2016b by using built-in analysis functions, except for the correlation coefficients (Fig. 4d) where we computed the mean and S.E.M. by 10,000 random resampling and performed the z-test. In a size-controlled procedure for the CP data, we drew the same number of samples as the SCN data (n = 58).

The simulation codes of illustrative examples in this study are available at modelDB (https://senselab.med.yale.edu/modeldb/↗) by an accession number 238338. Data analysis codes are available at https://github.com/JihwanMyung/ImagingAnalysis↗ (ImagingAnalysis.m and PMTAnalysis.m) as Mathematica packages. Further codes and data are available on request from the corresponding authors.

The simulation codes of illustrative examples in this study are available at modelDB (https://senselab.med.yale.edu/modeldb/↗) by an accession number 238338. Data analysis codes are available at https://github.com/JihwanMyung/ImagingAnalysis↗ (ImagingAnalysis.m and PMTAnalysis.m) as Mathematica packages. Further codes and data are available on request from the corresponding authors.