pii: sp- 00581-15

http://dx.doi.org/10.5665/sleep.6100

BASIC SCIENCE

Spatiotemporal Organization and Cross-Frequency Coupling of Sleep Spindles in Primate Cerebral Cortex Saori Takeuchi, PhD1,5; Rie Murai, MMedSc1; Hideki Shimazu, MD, PhD3; Yoshikazu Isomura, MD, PhD4; Tatsuya Mima, MD, PhD2; Toru Tsujimoto, MD, PhD1,5 Supportive Center for Brain Research, National Institute for Physiological Sciences, Okazaki, Japan; 2Human Brain Research Center, Kyoto University Graduate School of Medicine, Shogoin, Sakyo-ku, Kyoto, Japan; 3McGovern Institute for Brain Research, Massachusetts Institute of Technology, Cambridge, MA; 4Tamagawa University Brain Science Institute, Tokyo, Japan; 5The Graduate University for Advanced Studies (SOKENDAI), Shonan Village, Hayama, Kanagawa, Japan 1

Study Objectives: The sleep spindle has been implicated in thalamic sensory gating, cortical development, and memory consolidation. These multiple functions may depend on specific spatiotemporal emergence and interactions with other spindles and other forms of brain activity. Therefore, we measured sleep spindle cortical distribution, regional heterogeneity, synchronization, and phase relationships with other electroencephalographic components in freely moving primates. Methods: Transcortical field potentials were recorded from Japanese monkeys via telemetry and were analyzed using the Hilbert-Huang transform. Results: Spindle (12–20 Hz) current sources were identified over a wide region of the frontoparietal cortex. Most spindles occurred independently in their own frequency, but some appeared concordant between cortical areas with frequency interdependence, particularly in nearby regions and bilaterally symmetrical regions. Spindles in the dorsolateral prefrontal cortex appeared around the surface-positive and depth-negative phase of transcortically recorded slow oscillations (< 1 Hz), whereas centroparietal spindles emerged around the opposite phase. The slow-oscillation phase reversed between the prefrontal and central regions. Gamma activities increased before spindle onset. Several regional heterogeneities in properties of human spindles were replicated in the monkeys, including frequency, density, and inter-cortical time lags, although their topographic patterns were different from those of humans. The phaseamplitude coupling between spindle and gamma activity was also replicated. Conclusions: Spindles in widespread cortical regions are possibly driven by independent rhythm generators, but are temporally associated to spindles in other regions and to slow and gamma oscillations by corticocortical and thalamocortical pathways. Keywords: gamma oscillations, slow wave, cortical field potential, monkey, non-rapid eye movement sleep Citation: Takeuchi S, Murai R, Shimazu H, Isomura Y, Mima T, Tsujimoto T. Spatiotemporal organization and cross-frequency coupling of sleep spindles in primate cerebral cortex. SLEEP 2016;39(9):1719–1735. Significance The sleep spindle, an electroencephalographic feature of non-rapid eye movement stage sleep, has been implicated in memory consolidation and other forms of information processing. These functions likely depend on specific spatiotemporal patterns of brain activity, but these issues are difficult to address in human subjects. Consequently, we recorded sleep spindles directly from the monkey cortex using implanted electrodes during natural sleep via telemetry. We determined that spindles in widespread cortical regions are likely driven by independent rhythm generators, but are temporally associated to spindles in other regions and to slow and gamma oscillations. We suggest that spindle-spindle interactions are mediated by local corticocortical, commissural, and thalamocortical pathways and that these spindle-slow wave-gamma interactions are involved in cortical information processing.

INTRODUCTION The mammalian brain processes information not only during wakefulness but also during sleep. For example, sleep consolidates memory and enhances motor skill learning.1,2 Although the detailed neural mechanisms of such information processing are still unclear, they may be reflected by characteristic sleep-related electroencephalographic (EEG) patterns, such as sleep spindles, k-complexes, and slow waves. Sleep spindles are brief (0.5–2 s) 11–16 Hz rhythmic waves on the cortical EEG.3,4 They often occur with kcomplexes during non-rapid eye movement (NREM) sleep stage 2 (N2) and slow oscillations during slow wave sleep (SWS).5 In humans, qualitatively similar sleep spindles show regional heterogeneity in frequency,6–10 with higherfrequency (faster) spindles measured from centroparietal cortex and slower spindles from frontal cortex. EEG studies have shown that the faster centroparietal spindles (13–15 Hz) often preceded the slower frontal spindles (9–12 Hz) by about 0.2–0.5 s.11,12 Spindles are also temporally related to slow oscillations below 1 Hz.13 In humans, centroparietal fast spindles ride on the positive phase of slow oscillations, whereas slower SLEEP, Vol. 39, No. 9, 2016

spindles in the anterior middle region ride on the negative phase.12 It is reported that slow oscillations are traveling waves that occur more frequently in the frontal cortex and propagate from medial anterior to posterior14 and that they occur locally.15 Spindles are also temporally associated with higher frequency gamma oscillations. Gamma oscillations (30–120 Hz) are measured from neocortex and hippocampus during both sleep and wakefulness.16–18 In humans, cortical gamma oscillations occur during SWS, although the power is higher in rapid eye movement (REM) sleep.19,20 Phase-amplitude coupling of gamma activity with the spindle rhythm has been reported in humans.21,22 Recent studies have suggested that cross-frequency coupling between neural oscillations may serve to transfer information in large-scale brain networks.23 However, few studies have examined the coupling of spindles with slow and gamma oscillations. Several lines of evidence indicate that spindles may be related to the consolidation of declarative and procedural memories. In humans, spindle density increases after learning,24–26 especially during SWS,27 and a similar phenomenon has been reported in rodents.28 Furthermore, it is reported that auditory closed-loop in-phase stimulation

1719

Sleep Spindles in Primate Cerebral Cortex—Takeuchi et al.

Figure 1—Arrangement of recording electrodes. (A) Symbols indicate the position of electrodes. Symbols placed adjacent to sulci stand for recording sites in the banks of sulci. (B) Schematic examples of surface (S) and depth (D) electrodes. (C) The number of electrodes is shown for each area. Abbreviations: A, area; SMA, supplementary motor area; M1, primary motor area; S1, primary somatosensory area.

of slow oscillations during SWS enhances slow oscillation rhythms, phase-coupled spindles, and consolidation of declarative memory.29 In the two-stage model of memory consolidation, two separate memory stores are assumed, temporary and long-term.30 Both REM and NREM sleep have been proposed to participate in memory transfer between these two stores.1,31–33 According to the sequential hypothesis of memory consolidation, recently encoded memory representations are reactivated and stored in both the temporary and cortical long-term stores during SWS for integration into long-term memory, while REM sleep stabilizes transformed memories. In agreement with this hypothesis, temporal correlation between cortical activities (spindles and slow waves) and hippocampal gamma oscillations (ripples) was reported during SWS in rodents.34 During SWS in monkeys, activation of wide cortical regions associated with hippocampal ripples was reported using functional magnetic resonance imaging.35 These findings suggest that the neuronal activity underlying sleep spindles is part of a cortical information processing system during sleep that also involves slow (< 1 Hz) and gamma oscillations. In this study, we recorded field potentials during natural sleep directly from the monkey cortex using implanted electrodes across the dorsolateral and medial frontal cortex, as well as the parietal cortex, and investigated the physiological properties of spindles, including the rate, frequency, and timing of occurrence, inter-cortical interaction, and relation with slow and gamma oscillations. SLEEP, Vol. 39, No. 9, 2016

METHODS Subjects The experiment was conducted on 3 adult (6–8 kg) Japanese monkeys (Macaca fuscata), one male (SA) and 2 females (SB and SC) in accordance with the National Research Council/ Institute for Laboratory Animal Resources Guide for the Care and Use of Laboratory Animals and as approved by the institutional ethics committee of the National Institute of Natural Sciences. The present data were recorded as a part of a larger study protocol for other purposes. Surgical Operation Electrodes were implanted into cortical areas 9, 46, 8, 6, 4, 1, 5, 7, 32, and 24 (Figure 1) under ketamine (5 mg/kg administered intramuscularly) and pentobarbital sodium anesthesia. The initial dose (25 mg/kg) was administered intravenously, followed by additional injections as needed. Electrodes were made of silver wire (diameter 0.2 mm) coated with Teflon. They were cut obliquely at the tip (45 degrees) and the bare cross-section surfaces served as contact areas. Electrodes were arranged in pairs, with one at the surface (S) and the other (D) 2.5–3.5 mm below the cortical surface (Figure 1B), and fixed to the skull by acrylic resin and screws. Electrooculogram (EOG) electrodes were buried in the superolateral part of the orbital bone, and 2 coupled reference electrodes were set in the marrow of the parietal bones behind both ears. Electromyogram (EMG) electrodes were made of two

1720

Sleep Spindles in Primate Cerebral Cortex—Takeuchi et al.

Figure 2—Automated detection of spindles. (A) Segments of the raw surface (S) – depth (D) waveforms recorded at A8a, S1 (two sites), A6, and A24 in Monkey SA during stage N4. Pink and green circles indicate START and END of spindles automatically detected. Orange lines are band-pass filtered (0.1–1 Hz) waveforms. (B) The first step of spindle detection. Raw data (black) were band-pass filtered (light blue) and the envelope was extracted (blue). Spindles were detected (red envelope) at the mean + 3SD threshold (red broken line) of the envelope. START (pink circle) and END (green circle) were determined by the mean + SD line (green broken line) (see text). (C) The second step of detection. The phase difference between S and D within a spindle was measured at every peak of D (vertical blue lines). (D) The phase difference between S and D is displayed for the spindle in C by a circular histogram. (E) The third step of detection. The spectral profile of the spindle was examined on a time–frequency representation. The interval between START and END (pink and green vertical lines) was separated into 4 sections (a-d) and the marginal spectra were calculated for each section. (F) The marginal spectra. Note that all spectra have a peak within the spindle band (12–20 Hz: the interval between horizontal lines).

Teflon-insulated 7-strand stainless steel wires (individual wire diameter of 0.05 mm, total external diameter including the coating of 0.23 mm), with the coating peeled off 5 mm from the tip, and placed subcutaneously on the posterior nuchal muscles 20 mm apart. Recording Monkeys could move freely in their cages and slept from approximately 22:00 to 06:00. Room lights were turned off from 21:00 to 08:00. Starting 2 weeks after surgery, cortical field potentials, EOG, and EMG were recorded using a telemetry system (time constant, 1 s; third-order Butterworth low-pass filter, 200 Hz; sampling rate, 500 Hz; 16 channels). Nocturnal behavior was monitored by infrared video cameras. Each video frame was time stamped to telemetric data for offline analysis. SLEEP, Vol. 39, No. 9, 2016

Sleep Stage Classification Data were segmented into 20-s epochs and visually classified into wakefulness, REM, or NREM stages N1, N2, N3, or N4 as for human sleep classification.5 Automatic Detection of Sleep Spindles The current standard for sleep studies, The AASM Manual for the Scoring of Sleep and Associated Events,4 states the frequency range of the human sleep spindle to be 11–16 Hz (most commonly 12–14 Hz). However, the monkey counterpart for this value may not be the same. Preliminary visual inspection revealed waxing and waning sinusoidal 12–20 Hz oscillations of about 1-s duration, predominantly in NREM sleep stages and frequently preceded by k-complexes (Figure 2A and Figure S1, S2, and S3 in the supplemental material). Accordingly, we set the frequency range for sleep spindle detection to 12–20 Hz in

1721

Sleep Spindles in Primate Cerebral Cortex—Takeuchi et al.

this study. Candidate sleep spindles were detected (Figure 2B) according to a previous human study11: (1) S-D field potentials were band-pass filtered (12–20 Hz), (2) the envelope was calculated using the Hilbert transform, (3) spindles were detected when the envelope exceeded mean envelope amplitude + 3SD, and (4) the START and END of spindles were defined as the time points where the envelope crossed mean + SD. To confirm that the source currents of the spindles were generated at the recording site and not from distant cortical regions, the phase of the S signal relative to the D signal was calculated at every D signal peak. Spindles were accepted only when the mean phase angle from START to END was in antiphase (180° ± 90°) (Figure 2C, 2D). To avoid false detection caused by wide frequency band events, all candidate spindles were examined on a time-frequency distribution calculated using the Hilbert-Huang transform.36,37 Time-dependent Hilbert amplitude spectra were mapped on a plane with the frequency axis divided into 200 logarithmic bins from 0.1 to 125 Hz (Figure 2E). A Gaussian filter with a full width at half maximum (FWHM) of 10 bins × 20 ms was then applied. Based on the time-frequency distribution, marginal spectra were calculated for every fourth segment of the time region between START and END (Figure 2F). Spindles were accepted only when all 4 spectra reached local maxima within the range 12–20 Hz. Spindles that passed all the aforementioned criteria were subjected to further analysis. Frequency of Sleep Spindles The frequency of individual spindles was taken as the maximum measured within the range 12–20 Hz in marginal spectra from START to END. Temporal Order of Multiple Sleep Spindles To examine the temporal relationships among spindles measured from different cortical regions, cross-correlograms of START pulses were calculated for all recording site combinations. To estimate the probability density function of spindle occurrence, the cross-correlogram was smoothed using a Gaussian kernel with FWHM of 0.5 s. Significance of correlation was evaluated by comparison to the mean and SD of the same data after random shuffling. The time lag between peaks in the cross-correlogram yields information on how spindles at one recording site tend to precede or follow those at another site. To estimate the timing sequence of spindles over all recording sites using the telemetric data, which could cover only a subset of recording sites simultaneously, we assumed that the time lags are constant for each site pair in a given monkey. By repeatedly measuring the time lags for all available recording site pairs across days, simultaneous linear equations were constructed. We collected sufficient numbers so that the equations became an overdetermined system. The best estimate of the relative time lags between all recording site pairs was determined by a leastsquares method.

evaluation of spindle concordance between recording sites, the concordance index Q was defined as Q=

where NA and NB are the total number of spindles observed at recording sites A and B, and NAB is the number of spindles at A temporally correlated with spindles at B. NAB is calculated as the area exceeding the mean line under the cross-correlogram peak within the time interval in which the cross-correlogram exceeds mean + 2SD. The value of Q ranges between 0 and 1. A Q = 1 indicates complete correlation where all spindles in A and B take place in a correlated manner, while Q = 0 indicates that none of spindles are correlated. Q is a measure to assess the ratio of associated spindles to total spindles, taking into account the time lags between spindles. Frequency of Concurrent Spindles To examine whether spindle frequency is affected by spindles generated at other regions, the frequency distribution of concurrent spindles was compared to that of non-concurrent spindles, and the null hypothesis of no effect was tested by two-way analysis of variance (ANOVA) using recording site and concurrence as main factors. In addition, linear correlation of frequency between concurrent spindles was assessed using Pearson product-moment correlation coefficient, and the null hypothesis of no linear correlation was tested. In both cases, results were assessed at a significance level of 5%. Two spindles were regarded as concurrent when their START-END intervals overlapped. Cross-Frequency Coupling between Spindles and Slow Oscillations To examine the coupling of spindles to slow potential shifts, S-D potentials aligned at each spindle START were averaged. In addition, to analyze phase-amplitude coupling between spindles and slow oscillations, we measured the phase of slow oscillations (0.1–1 Hz) at each spindle START using only time regions in which slow oscillations phase was reversed between S and D. The instantaneous phase of slow oscillations was calculated using the Hilbert transform. Phase Relationship between Slow Oscillations The phase relationships among slow oscillations (0.1–1 Hz) at various recording sites were also examined. S-D, S, and D potentials were separately band-pass filtered (0.1–1 Hz) using Fourier transform and instantaneous phases were calculated using the Hilbert transform. To minimize the effects of current spread, the phase difference between two S-D potentials was calculated exclusively for data segments in which slow oscillation phases were reversed (180° ± 90°) between S and D.

Ratio of Concordant Spindles between Separate Regions Statistical significance in the cross-correlogram is not a quantitative measure of spindle concordance rate. For quantitative SLEEP, Vol. 39, No. 9, 2016

NAB × 2 NA + NB

1722

Statistical Analysis of Phase Angle For statistical analysis of phase locking, the null hypothesis of uniformity for the phase angle was tested using the Rayleigh statistic38: Sleep Spindles in Primate Cerebral Cortex—Takeuchi et al.

|∑r| R= ∑|r| where r is a unit vector of phase angle. The value of R ranges from 0 to 1, and R = 1 indicates complete phase locking. Cross-Frequency Coupling between Spindles and Gamma Oscillations To examine amplitude-amplitude coupling between gamma oscillations and spindles, the time-frequency distribution was averaged separately at the time point of spindle START, at the maximal amplitude of the spindle envelope (MAX), and at the time point of spindle END. The result was normalized to z-scores using a surrogate time-frequency distribution calculated from randomly sampled non-spindle regions of the same data. The threshold for detection of gamma increase was set to 95% confidence interval after Bonferroni correction for multiple comparisons. The start and end of gamma activity were defined as the earliest and latest time points of the time-frequency region in which the gamma activity (40–120 Hz) exceeded the threshold. To examine whether the recorded gamma band activity was generated at the recording site, the phase of S signal relative to the D signal was measured at every D signal peak in the time region from gamma start to spindle START. The time region after spindle START was excluded to avoid possible interference by spindle harmonics. The significance of result was assessed using the Rayleigh statistic R. To obtain an overview of coupling between spindle phase and gamma amplitude, time-frequency distribution was averaged separately at the time point of the first negative spindle peak after START, at the time point of the nearest spindle peak to MAX, and at the time point of the last spindle peak before END. To assess the coupling throughout the spindle periods, phase-dependent Hilbert amplitude spectra were mapped on a plane with the phase axis divided into 100 bins from 0° to 360°. Computational Procedures To calculate band-pass filtered signals, we used a frequency domain filter based on Fourier transform. In the first step, time series data were segmented into epochs of 32.768 s (16,384 points, giving a resolution of 0.031 Hz) and transformed into frequency domain using discrete Fourier transform. Secondly, a Tukey window (a rectangular window with the first and last 10% of the width equal to parts of a cosine function) was applied to data for band-pass operation in frequency domain. Finally, data were returned back to time domain using inverse discrete Fourier transform. The phase shift in this filtering was zero. Only the central 20 s of filtered signal was used for subsequent analysis. Continuous data were processed with sliding the data window in 20-s steps. The central 20 s corresponded to the 20-s epoch used in the visual sleep stage classification. Time-frequency distribution was calculated using the Hilbert-Huang transform. In the first step of this method, the oscillatory signal is decomposed into a set of intrinsic mode functions (IMF) using empirical mode decomposition invented by Huang.36 An IMF is a monocomponent signal with a zero mean. In the second step, the instantaneous phase and SLEEP, Vol. 39, No. 9, 2016

instantaneous amplitude of each IMF are calculated using the Hilbert transform. Instantaneous Hilbert spectrum is calculated as a time derivative of the phase. Instantaneous spectrum of the original signal is defined as superposition of IMF spectra and can be calculated at every sampling point with high temporal resolution. The high temporal resolution of spectra is an advantage of the Hilbert-Huang transform37 over other methods, such as short-time Fourier transform, continuous wavelet transform, and filtering techniques, in which time-frequency uncertainty is theoretically inevitable.39 Recording Sites To confirm record sites after recording, the monkeys were deeply anesthetized with pentobarbital sodium (50 mg/kg administered intravenously) and perfused through the heart with 10% formaldehyde neutral buffer solution. The electrode tracks were identified using postmortem magnetic resonance imaging and plotted on a standard brain map. RESULTS Recording Conditions and Sleep-Wake Cycles Cortical EEG, EOG, and EMG were recorded from all 3 monkeys for a total of 95 days (Monkey SA, 25 days; SB, 19 days; SC, 51 days). Cortical EEG was recorded from 88 recording sites (SA, 30; SB, 32; SC, 26) (Figure 1A and 1C). Total sleep duration was 763.4 h (SA, 200.3 h; SB, 149.8 h; SC, 413.3 h), and the total durations of individual sleep stages were 133.8 h for N1 (SA, 31.6 h; SB, 28.7 h; SC, 73.5 h), 161.7 h for N2 (SA, 51.6 h; SB, 39.4 h; SC, 70.2 h), 171.4 h for N3 (SA, 37.6 h; SB, 30.2 h; SC, 103.6 h), 207.8 h for N4 (SA, 50.1 h; SB, 35.7 h; SC, 122.0 h), and 89.1 h for REM (SA, 29.4 h; SB, 15.8 h; SC, 44.0 h). The EEG, EMG, and EOG patterns of each monkey exhibited characteristics in each sleep stage similar to those observed in human polysomnography, including spindles, kcomplexes, slow waves, REM-associated rapid eye movement, and the expected periods of EMG activation and quiescence (for example, Figure S1, S2, and S3 in the supplemental material). The number of NREM-REM cycles was (mean ± SD) 7.9 ± 2.1/night (SA, 9.1 ± 2.2; SB, 7.6 ± 1.7; SC, 7.4 ± 2.0). General Statistics of Spindles A total of 183,424 sleep spindles were recorded during N2–N4 (N2, 48,607; N3, 60,718; N4, 74,099). Mean spindle density (number per minute) per electrode was 1.9 ± 1.3/min with no significant difference between sleep stages (N2, 1.9 ± 1.0/min; N3, 2.0 ± 1.3/min; N4, 1.8 ± 1.5/min; one-way repeated measures ANOVA, F = 0.47, df = 2, P > 0.6). Spindles were identified at 83 of 88 recording sites (SA, 28 of 30; SB, 30 of 32; SC, 25 of 26), including areas 46d, 46v, 9, 8a, 8b, 6, 5, 7, 32, and 24, as well as primary motor (M1), primary somatosensory (S1), and supplementary motor areas (SMA) (Figure 3C). Mean spindle duration was 0.83 ± 0.26 s (N2, 0.83 ± 0.27 s; N3, 0.84 ± 0.27 s; N4, 0.81 ± 0.25 s). Spindle Amplitude Amplitude of spindle was measured as the amplitude of the envelope of band-passed spindle (12–20 Hz). Peak spindle

1723

Sleep Spindles in Primate Cerebral Cortex—Takeuchi et al.

Figure 3—Frequency and density of spindles. (A) The mean frequency spectra of spindles (red line) and non-spindle periods (blue line) at A8a of Monkey SA. Logarithmic scale is used for frequency axis. (B) The frequency distributions of spindles at A8a, S1, and A24 of Monkey SA. The positions of recording sites are shown in C with corresponding asterisk, dagger, and double-dagger. (C) The mean frequency and density of spindles at the indicated recording sites. Symbol color indicates frequency (blue for higher and red for lower) and size stands for density. Open symbols mark sites with few spindles (< 0.1/ min). (D) Plots of spindle density (red) and frequency (blue) for each recording site. The dots and bars represent mean and standard error of the mean. (E) Correlation between spindle density and frequency.

amplitude was 39 ± 26 µV (mean ± SD) during N2, 33 ± 23 µV during N3, and 31 ± 23 µV during N4. Although the amplitude of individual spindles was largely overlapped, the difference was statistically significant between N2 and N3, and between N3 and N4 (one-way repeated measures ANOVA, F = 822, df = 2, P < 10−300; Scheffe test, P < 10−16). Average spindle amplitude was 23 ± 14 µV during N2, 20 ± 13 µV during N3, and 19 ± 13 µV during N4. The difference was statistically significant between N2 and N3, and between N3 and N4 (one-way repeated measures ANOVA, F = 838, df = 2, P < 10−300; Scheffe test, P < 10−16). Spindle Frequency Spindle frequency was examined on marginal spectra obtained from time–frequency distribution. To confirm spindle detection, the mean spectrum during the START-END interval was compared to that of non-spindle periods during N2-N4 (Figure 3A). The START-END spectrum (red) showed a peak in the spindle band (12–20 Hz), while the non-spindle period spectrum (blue) showed no clear peak, confirming successful spindle detection. The frequency of individual spindles was measured and the distribution was plotted for areas 8a (A8a), 1 (S1), and 24 SLEEP, Vol. 39, No. 9, 2016

(A24) for Monkey SA (Figure 3B). The frequency of spindles was significantly different between cortical areas (one-way ANOVA, F = 17.1, df = 12, P < 3 × 10−16). Mean frequencies for each channel in all three monkeys are displayed on a brain map (Figure 3C and Figure S4 in the supplemental material) to illustrate the topographic frequency distribution. The fastest spindles (> 15 Hz, blue) were distributed primarily over the dorsolateral prefrontal cortex (DLPFC) (areas 46 and 8a); spindle frequencies in areas 46d and 46v were significantly higher than in all other areas, except area 8a (Scheffe test, P < 0.05). Spindles of intermediate frequency (14–15 Hz, green-orange) were seen mainly in the primary motor and parietal cortices (M1/ PC) and in the medial prefrontal cortex (mPFC); the slowest spindles (< 14 Hz, red) appeared primarily in the anterior cingulate cortex (ACC) (areas 32 and 24). There were also significant topographic differences in spindle density (one-way ANOVA, F = 11.6, df = 12, P < 3 × 10−12, Figure 3D) and a significant positive linear correlation between the frequency and the density (r = 0.864, P < 0.0002; Figure 3E). When compared between cortical regions, the density was highest in the DLPFC and lowest in the ACC (one-way ANOVA, F = 22.7, df = 3, P < 2e−10; Scheffe test P < 0.05) and was intermediate in the mPFC and M1/PC.

1724

Sleep Spindles in Primate Cerebral Cortex—Takeuchi et al.

In the histogram of A24 (Figure 3B), the left tail seems to be cut, suggesting that spindles of the lowest frequencies may be precluded by the criterion of frequency range (12–20 Hz). However, we abandoned the attempt to use a wider frequency range (for example, 10–20 Hz), because it introduced falsepositive spindle detection. Spindle Density Changes throughout the Night and during the NREM-REM Cycle In some recording sites, spindles density was not constant throughout the night. To visualize the changes in density, the cumulative number of spindles was plotted against the cumulative N2–4 time throughout the night (Figure S5A in the supplemental material). If the spindle density was constant, there should be a linear correlation. The resulting plots were often concave up or down (Figure S5A top). To combine data across recording days and to calculate the mean cumulative function, the duration of N2–4 in each recording day was normalized to 1, and the total number of spindles in each recording site was also normalized to 1 (Figure S5A below, blue lines). The 95% confidence intervals of the cumulative function were estimated using bootstrapping (Figure S5A below, red lines).40 The normalized median time of < 0.5 indicates that the spindles were generated more in the early night than in the late night, whereas that of > 0.5 indicates the opposite tendency. The null hypothesis that the median time equals to 0.5 was statistically tested by bootstrapping after Bonferroni correction for multiple comparison (P < 0.05). Among 83 spindle sources, a total of 22 sites (SA, 8; SB, 0; SC, 14) were classified to be the normalized median time of < 0.5, and 21 (SA, 2; SB, 18; SC, 1) were classified to be of > 0.5 (Figure S5B). Although it appears difficult to find a general trend from the result, the site with the largest median time was in S1 in all three monkeys. To examine spindle density changes during the NREMREM cycle, the same analysis was carried out on a sleep-cycle basis (Figure S6A in the supplemental material). A total of 19 sites (SA, 8; SB, 6; SC, 5) were of < 0.5 normalized median time and 5 sites (SA, 1; SB, 4; SC, 0) were of > 0.5. There seems no clear topographical pattern in the result (Figure S6B).

for multiple comparisons). Among 210 available recording site pairs, 170 pairs showed significant cross-correlation. The area in the elevated cross-correlation corresponds to the number of spindles that occurred concordantly. The concordance index (Q) is defined so that it indicates the ratio of concordantly occurring spindles to the total spindles. The concordance index (Q) distribution (Figure 4B) showed that 75% of pairs had Q values below 0.1, with highest value of 0.48. Higher Q values (Q > 0.1) were seen in several bilaterally symmetrical regions (bilateral areas 9–9 and S1–S1) and between nearby regions in the same hemisphere (e.g., area pairs 46–9, 46–46, 9–9, 9–8, 8–8, 9–32, 9–SMA, 8–SMA, SMA–SMA, 5–S1, and S1–S1; Figure 4C, 4D and Figure S7 in the supplemental material). Thus, although most spindles occurred independently, there were significant correlations between certain regions. Temporal Order of Spindles in Different Regions Concordant spindles tended to occur in a certain temporal order. To examine the temporal relationship between spindles in different regions, we measured time lags of the peaks in crosscorrelograms for all available recording site pairs (e.g., −0.12 s for 8a–S1, 0.64 s for 24–S1; Figure 4A) and mapped them on the brain image as relative values to the spindles in right S1 (Figure 5A and Figure S8 in the supplemental material). Oneway ANOVA revealed a significant time lag between cortical areas (F = 14.6, df = 12, P < 4 × 10−14). In Figure 5B, the areal pairwise mean time lags are shown in the order from earliest to latest (earliest, A8a; latest, A24). In general, it shows that spindles in the DLPFC and M1/PC precede those in the mPFC and that those in the ACC are the latest.

Spindle Frequency Changes throughout the Night and during the NREM-REM Cycle A linear correlation between the spindle frequency and the sleeping time was examined using analysis of covariance. We found no significant linear correlation throughout the night (F = 3.15, df = 1, P > 0.07) or during the NREM-REM cycle (F = 1.95, df = 1, P > 0.16). Concordant and Non-Concordant Spindles Spindles in different recording sites sometimes occurred at nearly the same time. In the example raw data in Figure 2A, spindles in different areas appeared sequentially within a short interval. To assess whether this occurred by chance, crosscorrelation of START times was calculated (Figure 4). In Figure 4A, cross-correlation is plotted for area 8a vs. S1 (left panel) and area 24 vs. S1 (right panel). In both cases, a significant cross-correlation was observed (area 8a vs. S1, P < 2 × 10 −104; area 24 vs. S1, P < 4 × 10−22 after Bonferroni correction SLEEP, Vol. 39, No. 9, 2016

1725

Frequency Interaction of Concurrent Spindles In the above results, we saw that spindles in different areas have different mean frequencies. If there is only one spindle rhythm generator in the brain, concurrent spindles should have the same frequency. If there are multiple rhythm generators, then are they independent? We then examined whether concurrent spindles exhibit mutual frequency modulation by comparing the frequency distribution of concurrent to non-concurrent spindles (Figure 6A). However, even if the overall frequency distribution of concurrent spindles does not differ from that of non-concurrent spindles, there remains a possibility that some individual concurrent pairs do exhibit frequency correlation, so we also calculated linear correlation coefficients between concurrent pairs. Thus, frequency interaction for each recording site pair was classified by the following criteria: (1) a significant frequency difference in a given region pair for non-concurrent spindles (frequency difference, FD), (2) a significant effect of concurrence on frequency by ANOVA (frequency shift, FS), and (3) a significant linear correlation in frequency between concurrent pairs (frequency correlation, FC). Thus, FD classifies the property of non-concurrent spindles, and FS and FC indicate significant interactions between concurrent spindles. All 210 site pairs were classified into six patterns (where “+” stands for significant, “−” for nonsignificant): (a) [FD−, FS−, FC−] (72 pairs), (b) [FD+, FS−, FC−] (107 pairs), (c) [FD+, FS+, FC−] (3 pairs), (d) [FD−, FS−, FC+] (11 pairs), (e) [FD+, FS−, Sleep Spindles in Primate Cerebral Cortex—Takeuchi et al.

Figure 4—Concordant and non-concordant spindles. (A) The cross-correlation of START pulses is shown for A8a vs. S1 and A24 vs. S1 in Monkey SA. The green and red lines are the mean and mean + 2SD levels of the same data after randomly shuffling. The positions of recording sites are shown in D with corresponding asterisk, dagger, and double-dagger. (B) The distribution of concordance index Q. (C) The mean Q for various recording site pairs shown separately for contralateral and ipsilateral combinations. An open square represent nonsignificant pairs and gray squares areas where no data are available. (D) The recording site pairs with significant correlations are connected by colored lines. Black lines indicate non-significant pairs. Color scale is for B, C, and D.

FC+] (11 pairs), and (f) [FD+, FS+, FC+] (6 pairs). Theoretically, two additional patterns are possible, [FD−, FS+, FC−] and [FD−, FS+, FC+], but were not found. Examples of patterns (a)–(f) are shown in Figure 6A(a)–(f), respectively, and the cumulative result is displayed on the brain map (Figure 6B and Figure S9 in the supplemental material). Patterns (a) and (b) indicate no significant interaction between frequencies of spindles while patterns (c)–(f) indicate significant interactions. Thus, significant interactions between concurrent spindles were found at 31 site pairs, including several bilaterally symmetrical regions (area pairs 9–9 and S1–S1) and several nearby region pairs in the same hemisphere (for example, 46–9, 46–46, 9–9, 9–8, 9–32, 6–24, 9–SMA, 6–SMA, 8–SMA, SMA–SMA, 5–S1, 7–S1, M1–S1, and S1–S1) (Figure 6C and Figure S9). In the patterns (b), (c), (e), and (f) with [FD+], the frequency difference between concurrent spindles in two regions was also statistically significant, indicating that spindles of different frequencies occur concurrently. In the patterns of (d), (e), and (f) with [FC+], the correlation coefficient was positive for all 28 pairs, indicating a positive linear correlation. In the patterns (c) and (f) with [FD+, FS+], the difference in mean frequency between two areas was narrower in concurrent spindles than non-concurrent spindles in 7 of the total 9 pairs. SLEEP, Vol. 39, No. 9, 2016

Cross-Frequency Coupling between Spindles and Slow Oscillations Spindles were coupled with slow oscillation phase. In Figure 2A, two spindles in area 8a ride on approximately the same phase of slow oscillations (orange line) extracted by a band-pass filter (0.1–1 Hz). To investigate whether spindles are associated with slow potential shifts (e.g., more likely to occur during slow oscillations or a specific phase of a slow oscillation), we averaged raw data at START (t = 0). Example results are shown for areas 8a, S1, and 24 in Monkey SA (Figure 7A). Slow oscillations that have a period of more than 1 s can be seen around the START. Essentially the same effects were observed at other recording sites (data not shown). This indicated that spindles were superimposed on the slow component (< 1Hz). The polarity of the slow component and the phase relationship with spindles were dependent on cortical area. The phase angles of slow oscillations (0.1–1 Hz) at spindle START were measured (e.g., Figure 7B for the same recording sites as in Figure 7A) and plotted on circular histograms (Figure 7C, 0° indicates the negative peak and 180° the positive). The mean phase angle was calculated for all cortical spindle sources and shown on the brain map (Figure 7D). Phase locking was statistically significant for 58 of 83 cortical spindle sources

1726

Sleep Spindles in Primate Cerebral Cortex—Takeuchi et al.

(Rayleigh test, P < 0.05 after Bonferroni correction for multiple comparisons). The grand mean of mean phase angle within cortical regions was as follows: the DLPFC, 174° (N = 17, R = 0.72, P < 5e−5); ACC, 306° (N = 9, R = 0.60, P < 0.04); M1/ PC, 331° (N = 35, R = 0.40, P < 0.005); and mPFC, 265° (N = 10, R = 0.48, P < 0.11). Thus, spindles in the DLPFC tended to start on the positive phase of slow oscillations, while those in the ACC and M1/PC tended to occur on the negative phase. Spindles in the mPFC tended to start on the transitional slope from the positive to the negative, although the regional variation is relatively large. Phase Relationship between Slow Oscillations Slow oscillations themselves showed inter-regional correlations. To examine the phase relationship between slow oscillations in various cortical areas, the phase angle (0.1–1 Hz) was measured relative to simultaneously recorded slow oscillations at right S1. The mean phase angle in each recording site and the statistical significance of phase locking are shown on the brain image (Figure 8). Phase locking was statistically significant at 35 of 37 available recording sites (P < 0.05 after Bonferroni correction for multiple comparisons). Phase angle was reversed between the central (M1/S1) region and the DLPFC, mPFC, and ACC regions with a reversal boundary between M1 and area 6. The grand mean of mean phase angle within cortical regions was as follows: the DLPFC, 198° (N = 8, R = 0.95, P < 8e−5); mPFC, 185° (N = 5, R = 0.84, P < 0.02); ACC, 184° (N = 3, R = 0.99, P < 0.04); premotor cortex, 183° (N = 5, R = 0.55, P < 0.24), and M1/S1, 16° (N = 11, R = 0.72, P < 0.002). In addition, the phase angle appeared to reverse between the M1/S1 region and area 7 (mean phase angle in area 7 was 178° as referenced to right S1), although the phase difference between regions was not statistically significant (N = 3, R = 0.58, P > 0.4). Slow oscillations in bilaterally symmetrical regions appeared nearly synchronous (mean phase difference between bilaterally symmetrical regions in the same monkey was 14°), although the bilateral synchrony was not statistically significant (N = 4, R = 0.57, P > 0.29). Cross-Frequency Coupling between Spindles and Gamma Oscillations Gamma activity preceded spindles. In the time–frequency distribution of a single spindle (Figure 2E), low amplitude diffuse gamma activity (40–120 Hz) is seen before START in parallel with 1–5 Hz activity reflecting the k-complex and slow oscillations below 1 Hz. To examine the properties of gamma activity around spindles, the time–frequency distribution was averaged separately at spindle START and END, and the resultant maps normalized to z-scores. Example results are shown for areas 8a and S1 in Monkey SA (Figure 9A). The upper row shows maps aligned at START and the lower row maps aligned at END. In both areas 8a and S1, gamma activity significantly increased before START and the increase ended before END. Slow oscillations (< 1 Hz) and 1–5 Hz activity that presumably corresponds to k-complexes were also observed around the spindle. The mean gamma phase difference between S and D was in anti-phase (179° both areas 8a and S1) (Figure 9B), suggesting SLEEP, Vol. 39, No. 9, 2016

Figure 5—Time lags of spindles. (A) The relative time lags of spindle START at a given recording site relative to spindle START at right S1 are shown. The bordered symbols in right S1 indicate the position of the reference channel for each monkey. (B) The mean time lags between areas. Statistically significant lags are marked with squares (Scheffe test P < 0.05).

that these gamma activities were likely generated at the recorded sites. Recording sites were identified as significant gamma generators when both the gamma increase and phase reversal were statistically significant. The distribution of significant gamma activity is plotted with colored marks on the brain maps of Figure 9C, with color denoting the onset time range from spindle START at each recording site (early, red; late, blue). Among 83 cortical spindle sources, 42 were identified as significant gamma generators. Gamma generators were found in areas 9, 46, 6, and 8 in three monkeys, areas 24, and 7 in Monkey SA, area 32 and SMA in Monkey SC, and S1 in Monkeys SA and SC.

1727

Sleep Spindles in Primate Cerebral Cortex—Takeuchi et al.

Figure 6—Frequency of concurrent spindles. (A) Frequency of spindles is compared between a pair of recording sites (red and light blue). Four histograms in a column show frequency of non-concurrent spindles (top and bottom) and concurrent spindles (second and third). Vertical blue lines indicate mean frequency. The images below the histograms are xy-plots to evaluate the linear correlations between concurrent spindle pairs. The x and y axes correspond to the frequency in the second (light blue) and third (red) histograms, respectively. A Gaussian filter with a full width at half maximum of 0.7 Hz × 0.7 Hz was applied to the plots for presentation purposes only. The behavior of frequency for non-concurrent and concurrent spindles was examined by frequency difference (FD), frequency shift (FS), and frequency correlation (FC). By repeating the analysis across all available recording site pairs, each was classified into one of the six patterns of (a)–(f) (see Results). The positions of recording sites are shown in B with corresponding letters a–f. (B) The classification into the six patterns is displayed on the brain map. (C) The ratio of pairs with interaction for various recording site pairs shown separately for contralateral and ipsilateral combinations. Open squares represent pairs with no interaction and gray squares represent areas for which no data are available.

The mean time of gamma onset relative to START was −0.30 ± 0.17 s. The mean duration of gamma increase was 0.51 ± 0.31 s. The timing of gamma end was 0.21 ± 0.31 s from spindle START, −0.14 ± 0.30 s from spindle MAX, and −0.55 ± 0.32 s from spindle END. The statistical significance of the temporal relationship between gamma activity and spindles was tested for each recording site using bootstrapping (N = 1,000, P < 0.05 after Bonferroni correction for multiple comparisons).40 The temporal precedence of the gamma onset over the spindle START was confirmed in 40 sites among the 42 significant gamma generators (Figure 9C, filled mark) and the precedence of the gamma end over the spindle END in all 42 generators. The results indicate that short bursts of gamma oscillations emerged before the spindle START and ended before the spindle END. SLEEP, Vol. 39, No. 9, 2016

Frequency of the maximal z-score of the gamma increase was measured and displayed on the brain image (Figure S10 in the supplemental material). One-way ANOVA revealed a significant frequency difference between cortical areas (F = 8.48, df = 10, P < 2e−6). The gamma frequency in S1 was significantly higher than in A46v, A9, A8b, SMA, A32, or A24 (Scheffe test P < 0.05), although this result was based only on the observation in Monkeys SA and SC. Coupling between Spindle Phase and Gamma Oscillations Gamma activity was coupled with spindle phase. To obtain an overview of the cross-frequency coupling between spindle phase and gamma amplitude, the time–frequency distribution was averaged separately at the time point of the first spindle peak after START, at the time point of the nearest spindle peak

1728

Sleep Spindles in Primate Cerebral Cortex—Takeuchi et al.

Figure 7—Correlation between slow oscillations and spindles. (A) The averaged raw data at spindle START for recording sites A8a, S1, and A24 in Monkey SA. The recording sites are shown in D with corresponding asterisk, dagger, and double-dagger. (B) The slow-wave phase at START for A8a, S1, and A24 in Monkey SA. (C) The phase coordinate system. (D) The mean slow-oscillation phase at START. The symbol size indicates the degree of phase locking (Rayleigh statistic R). Open symbols denote sites where phase locking was not statistically significant.

to MAX, and at the time point of the last spindle peak before END. Then the resultant maps were normalized to z-scores. Example result (Figure 10A and Figure S11 in the supplemental material) shows that gamma activity is modulated by spindle phase. The modulation can be seen around START, MAX, and END even after the gamma activity had faded, as well as in the sites where no gamma increase can be noted at all (Figure S11D). The result also shows that spindle frequency is modulated by its own phase. To assess the phase-locked spectral modulation throughout the spindle period, phase-frequency distribution was mapped and normalized by two methods, one using surrogate data of non-spindle periods (Figure 10B) and the other using shuffled data of spindle periods (Figure 10C). Figure 10B depicts amplitude modulation of gamma activity and intra-wave frequency modulation of spindle (white line) by spindle phase. The modulation of gamma activity is more clearly shown in Figure 10C. The grand mean across the spindle sources was calculated for both normalization methods (Figure 10D and 10E). They indicate that the modulation patterns found in the examples (Figure 10B and 10C) are a common feature of the recording sites. To assess the intra-wave frequency modulation of the spindle (white line in Figure 10B), the mean frequency and standard error were plotted against phase as the relative value to the mean spindle frequency in each site (Figure 10F). SLEEP, Vol. 39, No. 9, 2016

The intra-wave frequency modulation of spindle was statistically significant (ANOVA, F = 106, df = 99, P < 10−300), and the maximal relative frequency was 1.05 at 328°. The result indicates that the instantaneous frequency of spindle is the highest shortly before the negative peak of spindles. Because high instantaneous frequency means fast phase rotation, this corresponds to the fact that negative peaks of spindles are more tapered than the positive troughs. Based on the grand mean phase–frequency map normalized by shuffled data of spindle period (Figure 10E), the grand mean gamma band activity and the standard error were calculated for 40–80 Hz and 80–120 Hz (Figure 10G). The modulation of gamma amplitude by the spindle phase was statistically significant for both the bands (ANOVA; 40–80 Hz; F = 60.5, df = 99, P < 100−300; 80–120 Hz; F = 47.5, df = 99, P < 100−300). The maximal relative amplitude was 1.11 at 277° and 1.05 at 313° for 40–80 Hz and 80–120 Hz, respectively. The result indicates that the gamma amplitude is the highest shortly before the negative peak of spindles. DISCUSSION Summary of the Results The present study revealed topographic features of (1) spindle density (high in the DLPFC, medium in the mPFC and M1/PC, and low in the ACC, and correlated with the spindle frequency),

1729

Sleep Spindles in Primate Cerebral Cortex—Takeuchi et al.

Figure 8—Correlation between slow oscillations. The mean slowoscillation phase at each recording site is mapped on the brain image. The bordered symbols in the right S1 indicate the reference channels for each monkey. Open symbols denote areas where phase locking was not statistically significant.

(2) spindle frequency (the same topographic tendency as the density), (3) spindle timing (early in the DLPFC and M1/PC, medium in the mPFC, and late in the ACC), (4) phase reversal of slow oscillations between the DLPFC, mPFC, and ACC region and the M1/S1 region, (5) coupling between spindles and slow oscillations (spindles in the DLPFC tend to occur on the positive phase, while those in the ACC and M1/PC on the negative), and (6) frequency of gamma oscillations associated with spindles (low in the DLPFC, mPFC, and ACC region, and high in S1) (Figure S12 in the supplemental material). It also showed that (7) most spindles occurred independently across regions, but (8) there was a significant tendency of co-occurrence in nearby regions and bilaterally symmetrical regions, (9) frequency of most concurrent spindles was independent, but there was significant interdependence between spindle frequencies in nearby regions and bilaterally symmetrical regions, (10) gamma activity preceded spindles, (11) gamma amplitude was coupled with spindle phase, (12) spindle showed intrawave frequency changes and the instantaneous frequency was highest shortly before the negative peak, and (13) spindle density during N2–4 increased or declined throughout the night depending on recording sites. The results (1), (2), (3), (5), (7) and (11) replicated and extended the human results, although the topographic patterns were different from the human patterns. Results (4), (6), (8), (9), (10), (12), and (13) are novel findings in the present study. SLEEP, Vol. 39, No. 9, 2016

Similarity to Human Spindles and Functional Significance The sleep spindles observed in monkey cortex appeared qualitatively similar to sleep spindles in humans in terms of (1) mean frequency range (13–17 Hz) and a topographic frequency distribution (albeit with a pattern distinct from humans), (2) duration (~1 s), (3) density (~1–3/min), (4) REM-NREM distribution (> in N2–N4), and (5) a characteristic temporal lag between cortical areas. More importantly, we show that spindles occurring in bilaterally symmetrical regions and proximal regions of the same hemisphere were often concordant with frequency modulation. Furthermore, spindles were associated with the particular phase of slow oscillations and preceded by gamma activity and the gamma activity was coupled with spindle phase, suggesting complex interactions reflecting cortical information processing during sleep. Spindle frequency was similar to human sleep spindles, with individual spindle frequencies ranging from 12 to 20 Hz and mean frequencies between 13–17 Hz over the cortex, only slightly higher than for human spindles (11–16 Hz).4,8,9,41–43 However, in the present study, the frequency criterion in the automatic detection of spindles was empirically determined to be 12–20 Hz based on visual inspection of raw data. We should note that the range setting may partly bias the results. Mean duration (0.83 ± 0.26 s) was also comparable to human spindles (~0.5–1.2 s)9,44 with mean density (1.9/min) and cortical distribution comparable to finding from intracranial recordings in human (1.8/min in frontal lobe, 1.3/min in parietal lobe).11 Moreover, NREM-REM modulation and other sleep stage-specific patterns in the EOG, EMG, and EEG (such as k-complexes and slow waves) appeared similar to those in the human polysomnograph. The number of sleep cycles (7.9 ± 2.1) per night was close to the 4–6 cycles in humans.45–47 Thus, monkey spindles may be a useful model for assessing the functions of human sleep spindles. Localization of Cortical Current Generators of Spindles In human EEG, sleep spindles are most frequently recorded over the central head regions48 with sources in dorsolateral and medial frontal cortex and parietal cortex.10,49–52 The present study confirmed spindle current generators in a wide region of cortex, including the DLPFC, mPFC, ACC, M1/PC, and premotor cortex, using direct recording with transcortical electrode pairs (Figure 3C). There seems to be no apparent spindle-free area within the regions explored by us. Topographic Pattern of Spindle Frequencies In humans, the frequency of centroparietal spindles is higher than that of frontal spindles.6,10 However, the topographic pattern of spindle frequencies in monkey differed in that frequency was higher in the DLPFC than the M1/PC region, possibly due to greater functional development of human prefrontal cortex (Figure S12B in the supplemental material). In addition, a significant positive correlation between spindle frequency and spindle density was found (Figure 3E) as in the previous human study,10 suggesting a close relationship between the spindle rhythm generator and the triggering system. The temporal density of spindles was often unequally distributed in the early or late night depending on the site (Figure S5 in the supplemental

1730

Sleep Spindles in Primate Cerebral Cortex—Takeuchi et al.

Figure 9—Correlation between gamma oscillations and spindles. (A) The mean normalized time–frequency distribution aligned at START (top) and END (below) for A8a and S1 in Monkey SA. Black contour lines represent the significance threshold at P < 0.05 after Bonferroni correction. Note that the increase in the gamma band precedes START. The positions of recording sites are shown in C with corresponding asterisk and dagger. (B) The phase difference of gamma oscillations between S and D at A8a and S1. (C) The onset time of gamma increase as measured from START at each recording site. Colored symbols stand for significant gamma generators. Filled symbols indicate sites where precedence of gamma onset over START was statistically significant. The mean onset time was −0.30 ± 0.17 s.

material), suggesting regional heterogeneity of the triggering system. Previous studies in cats under ketamine/xylazine anesthesia suggested that the recurrent inhibitory loop between thalamic reticular nucleus (TRN) neurons and thalamocortical neurons is critical for generation of spindle rhythm.53 In turn, the topographic patterns of spindle rhythm, density, and density change may reflect the topographic pattern of thalamocortical neurons with distinct biophysical properties, such variation in membrane potential, and microcircuit properties.

0.2–0.5 s.11,12 In monkey, spindles tended to occur sequentially in a lateral to medial order (Figure 5 and Figures S8 and S12C in the supplemental material) with a time lag of approximately 0.5 s, and there was a tendency for fast (DLPFC) spindles to precede slower (mPFC and ACC) spindles. However, no significant difference in timing was noted between the DLPFC and the M1/PC, although spindles in these areas had different frequencies. Therefore, the time lag may be more strongly correlated to cortical location than to frequency.

Concordant and Non-Concordant Spindles In humans, the majority of spindles occur independently in local cortical regions, while relatively few occur concomitantly across multiple regions.11,15,54 In monkey, most spindles also occurred in isolation, but some bilaterally symmetrical and proximal recording site pairs (Figure 4C and 4D) exhibited spindles with significant temporal concordance and frequency dependence (Figure 6). Such regions likely have relatively strong corticocortical connections via commissural or short association fibers (U-fibers).55–58 In addition, crossed thalamocortical and corticothalamic projections59,60 and crossed projections from the TRN to the thalamus61,62 may contribute to bilateral synchronization.

Multiple Spindle Rhythm Generators It is still debated whether fast and slow spindles are derived from independent rhythm sources or are simply the result of a single rhythm generator with temporal variance in output.1,63,64 The present results indicated that spindles of different frequencies can occur simultaneously in different regions (e.g., [FD+] in Figure 6), suggesting the existence of multiple rhythm generators. This is in accordance with human studies using magnetoencephalography, low-resolution brain electromagnetic tomography, or intracranial electrodes.11,49,65 Individual TRN-thalamic projection circuits with distinct properties and projecting topographically to different regions of cortex may constitute these multiple rhythm generators.53,66–70 At the majority of recording sites, spindle frequency was not affected by other concurrent spindles at other sites ([FS−, FC−] in Figure 6), suggesting that each rhythm generator acts

Temporal Organization of Spindles In human studies, fast spindles in the centroparietal cortex tend to precede slower spindles in the anterior cortex by SLEEP, Vol. 39, No. 9, 2016

1731

Sleep Spindles in Primate Cerebral Cortex—Takeuchi et al.

Figure 10—Modulation of spectra by spindle phase. (A) The mean normalized time–frequency distribution aligned at the first spindle peak after START (top), at the nearest spindle peak to MAX (middle), and at the last spindle peak before END (below) for the same recording sites as in Figure 9A. The averaged S-D field potentials (black waveforms) are superimposed. Note that both the gamma amplitude and the spindle frequency are coupled with the spindle phase. (B) The mean normalized phase–frequency distribution showing the spectra in relation to the phase of 12–20 Hz band-passed signals. The phase axis was divided into 100 bins, and instantaneous spectra during spindle periods that were calculated using the Hilbert−Huang transform were averaged for each phase bin. The resulted spectra were normalized using surrogate spectra during non-spindle periods. The maximal amplitude for each phase bin was connected (white line). (C) The same data as B were normalized using shuffled data during the spindle periods. (D) The grand mean phase– frequency distribution across the recording sites normalized by the mean spectra in non-spindle periods. (E) The grand mean phase–frequency distribution normalized by the mean spectra during spindle periods showing a gamma amplitude modulation. (F) Intra-wave frequency modulation of spindle (white line in B) was normalized using the mean spindle frequency of each recording site, and the mean and standard error are plotted. (G) Coupling between gamma amplitude and spindle phase. Based on the phase–frequency distribution in E, the mean gamma band amplitude was calculated for 40–80 Hz and 80–120 Hz and displayed with the standard error.

independently. However, two types of interactions were observed for specific site pairs (often between bilaterally symmetrical regions and between nearby regions with higher concordance index Q), a shift to reduce the frequency SLEEP, Vol. 39, No. 9, 2016

1732

difference ([FS+] in Figure 6) and a positive linear correlation ([FC+] in Figure 6). In these site pairs, corticocortical connections and crossed projections between the thalamus and the cortex or TRN may cause entrainment. Sleep Spindles in Primate Cerebral Cortex—Takeuchi et al.

Correlation between Spindles and Slow Oscillations Spindles in each cortical region occurred preferentially around a particular phase of concomitant slow oscillations at the recording site (Figure 7 and Figure S12E in the supplemental material). In the DLPFC, spindles occurred around the positive phase of slow oscillations, while in the M1/PC, spindles appeared around the negative phase. As the DLPFC produces the faster spindles and the M1/PC generates the slower ones, this result is in accord with human studies showing that fast spindles ride on the positive phase of slow oscillations, while slow spindles ride on the negative phase.12,13 In addition, the present study revealed a phase relationship among slow oscillations in different regions (Figure 8 and Figure S12D), with a tendency for phase reversal between the DLPFC, mPFC, and ACC regions and the M1/S1 region. The phase posterior to the intraparietal sulcus was also reversed relative to sensorimotor cortex, although this was based on fewer recording sites. Given the lateral to medial time lag and the spatial distribution of slow oscillation phase and fast and slow spindles, faster spindles in the DLPFC and slower spindles in M1/PC and ACC tended to occur on the opposite phases of slow oscillations (Figure S12C–S12E). We may be able to interpret the relatively complex regional pattern of coupling between the slow-oscillation phase and spindle onset as a result of relatively simple topographic properties, i.e., spindle time lags and inter-cortical slow-oscillation phase relationship. K-complexes may contribute to the correlation between slow oscillations and spindles because they often precede spindles, predominantly in frontal regions (for example, Figure 2A and Figure S2 in the supplemental material), and contain slow components (< 1 Hz) as well as delta band (1–4 Hz) activity.48,71–73 Recent observations in humans have shown that slow waves during SWS can be regulated locally15 and can propagate over the cortex as a traveling wave.14 In the present study, phase shift of slow oscillations was noted between cortical regions. Because traveling waves show phase differences between cortical regions, the present study result is in accordance with the findings in humans.

spindles, and gamma activity is reported in the hippocampus in humans22 and is proposed as a possible framework for the neural information processing. The present results revealed a similar hierarchical structure in the primate neocortex. Functional Significance The present study showed that spindles in many cortical regions are driven by their own rhythm sources, but nevertheless can be temporally and spatially correlated both to other spindles and to slow and gamma oscillations. Such interactions may serve as a functional basis to transfer information in largescale brain network during memory consolidation.23 Although detailed explanation of the neural mechanisms for these correlations is beyond the scope of the present study, a possible mechanism has been proposed by Steriade.75 According to the theory, cortically generated slow oscillations facilitate the corticocortical and corticothalamic synchronization of other rhythms, including spindles and gamma oscillations. Indeed, the most prominent spindle correlations were found between bilaterally symmetrical regions and between nearby regions, suggesting that corticocortical connections and crossed projections between the thalamus and the cortex or TRN may mediate this coupling. The cross-frequency coupling suggests that the triggering system for spindles is linked to slow oscillations and gamma oscillations. CONCLUSIONS The present study revealed a constellation of correlations between spindles and slow and gamma oscillations in monkey cortex. These complex spatiotemporal associations may underlie information processing during SWS.

Correlation between Spindles and Gamma Oscillations We also observed increased gamma oscillation amplitude preceding spindles START, with a time lag of 0.30 ± 0.17 s (Figure 9). The increased gamma activity lasted about 0.5 s, and then decreased before the time point of spindle END. To the best of our knowledge, this increase in gamma band activity before the spindle has not been reported previously. The gamma amplitude was correlated with spindle phase, and was highest shortly before the negative peaks of spindles (Figure 10G). At approximately the same spindle phase, the instantaneous frequency of spindle was highest (i.e., phase rotation was fastest) (Figure 10F), suggesting the linkage between the gamma activity and spindle oscillations. In human epileptic patients, it is reported that cortical gamma oscillations during NREM sleep emerge synchronously at many cortical areas preferentially on a specific phase of the slow wave,20,74 while coupling of gamma activity with spindles has been observed in humans.21 Hierarchical coupling of slow oscillations, SLEEP, Vol. 39, No. 9, 2016

1733

REFERENCES 1. Rasch B, Born J. About sleep’s role in memory. Physiol Rev 2013;93:681–766. 2. Walker MP, Brakefield T, Seidman J, Morgan A, Hobson JA, Stickgold R. Sleep and the time course of motor skill learning. Learn Mem 2003;10:275–84. 3. Loomis AL, Harvey EN, Hobart G. Potential rhythms of the cerebral cortex during sleep. Science 1935;81:597–8. 4. Iber C, Ancoli-Israel S, Chesson A, Quan SQ; American Academy of Sleep Medicine. The AASM manual for the scoring of sleep and associated events: rules, terminology and technical specifications. Westchester, IL: American Academy of Sleep Medicine, 2007. 5. Rechtschaffen A, Kales A. A manual of standardized terminology, techniques and scoring system for sleep stages of human subjects. Washington, DC: Washington Public Health Service, US Government Printing Office, 1968. 6. Jankel W, Niedermeyer E. Sleep spindles. J Clin Neurophysiol 1985;2:1–36. 7. Jobert M, Poiseau E, Jahnig P, Schulz H, Kubicki S. Topographical analysis of sleep spindle activity. Neuropsychobiology 1992;26:210–7. 8. Werth E, Achermann P, Dijk DJ, Borbély AA. Spindle frequency activity in the sleep EEG: individual differences and topographical distribution. Electroenceph Clin Neurophysiol 1997;103:535–42. 9. Zeitlhofer J, Gruber G, Anderer P, Asenbaum S, Schimicek P, Saletu B. Topographic distribution of sleep spindles in young healthy subjects. J Sleep Res 1997;6:149–55. 10. Peter-Derex L, Comte JC, Mauguiere F, Salin PA. Density and frequency caudo-rostral gradients of sleep spindles recorded in the human cortex. Sleep 2012;35:69–79.

Sleep Spindles in Primate Cerebral Cortex—Takeuchi et al.

11. Andrillon T, Nir Y, Staba RJ, et al. Sleep spindles in humans: insights from intracranial EEG and unit recordings. J Neurosci 2011;31:17821–34. 12. Mölle M, Bergmann TO, Marshall L, Born J. Fast and slow spindles during the sleep slow oscillation: disparate coalescence and engagement in memory processing. Sleep 2011;34:1411–21. 13. Mölle M, Marshall L, Gais S, Born J. Grouping of spindle activity during slow oscillations in human non-rapid eye movement sleep. J Neurosci 2002;22:10941–7. 14. Massimini M, Huber R, Ferrarelli F, Hill S, Tononi G. The sleep slow oscillation as a traveling wave. J Neurosci 2004;24:6862–70. 15. Nir Y, Staba RJ, Andrillon T, et al. Regional slow waves and spindles in human sleep. Neuron 2011;70:153–69. 16. Buzsáki G, Wang X-J. Mechanisms of gamma oscillations. Annu Rev Neurosci 2012;35:203–25. 17. Wang X-J. Neurophysiological and computational principles of cortical rhythms in cognition. Physiol Rev 2010;90:1195–268. 18. Takeuchi S, Mima T, Murai R, Shimazu H, Isomura Y, Tsujimoto T. Gamma oscillations and their cross-frequency coupling in the primate hippocampus during sleep. Sleep 2015;38:1085–91. 19. Gross DW, Gotman J. Correlation of high-frequency oscillations with the sleep-wake cycle and cognitive activity in humans. Neuroscience 1999;94:1005–18. 20. Valderrama M, Crepon B, Botella-Soler V, et al. Human gamma oscillations during slow wave sleep. PLoS One 2012;7:e33477. 21. Ayoub A, Mölle M, Preissl H, Born J. Grouping of MEG gamma oscillations by EEG sleep spindles. NeuroImage 2012;59:1491–500. 22. Staresina BP, Bergmann TO, Bonnefond M, et al. Hierarchical nesting of slow oscillations, spindles and ripples in the human hippocampus during sleep. Nat Neurosci 2015;18:1679–86. 23. Canolty RT, Knight RT. The functional role of cross-frequency coupling. Trends Cogn Sci 2010;14:506–15. 24. Gais S, Mölle M, Helms K, Born J. Learning-dependent increases in sleep spindle density. J Neurosci 2002;22:6830–4. 25. Fogel SM, Smith CT. The function of the sleep spindle: a physiological index of intelligence and a mechanism for sleep-dependent memory consolidation. Neurosci Biobehav Rev 2011;35:1154–65. 26. Schabus M, Hodlmoser K, Gruber G, et al. Sleep spindle-related activity in the human EEG and its relation to general cognitive and learning abilities. Eur J Neurosci 2006;23:1738–46. 27. Cox R, Hofman WF, Talamini LM. Involvement of spindles in memory consolidation is slow wave sleep-specific. Learn Mem 2012;19:264–7. 28. Eschenko O, Mölle M, Born J, Sara SJ. Elevated sleep spindle density after learning or after retrieval in rats. J Neurosci 2006;26:12914–20. 29. Ngo HV, Martinetz T, Born J, Mölle M. Auditory closed-loop stimulation of the sleep slow oscillation enhances memory. Neuron 2013;78:545–53. 30. Squire LR. Memory and the hippocampus: a synthesis from findings with rats, monkeys, and humans. Psychol Rev 1992;99:195–231. 31. Buzsáki G. Two-stage model of memory trace formation: a role for “noisy” brain states. Neuroscience 1989;31:551–70. 32. Peigneux P, Laureys S, Delbeuck X, Maquet P. Sleeping brain, learning brain. The role of sleep for memory systems. Neuroreport 2001;12:A111–24. 33. Diekelmann S, Born J. The memory function of sleep. Nat Rev Neurosci 2010;11:114–26. 34. Sirota A, Csicsvari J, Buhl D, Buzsáki G. Communication between neocortex and hippocampus during sleep in rodents. Proc Natl Acad Sci U S A 2003;100:2065. 35. Logothetis NK, Eschenko O, Murayama Y, et al. Hippocampalcortical interaction during periods of subcortical silence. Nature 2012;491:547–53. 36. Huang NE, Shen Z, Long SR, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc R Soc A 1998;454:903–95.

SLEEP, Vol. 39, No. 9, 2016

37. Huang NE. Introduction to the Hilbert-Huang transfrom and its related mathematical problems. In: Huang NE, Shen SS, eds. HilbertHuang transform and its applications. New Jersey: World Scientific, 2005:1–26. 38. Lord Rayleigh. On the resultant of a large number of vibrations of the same pitch and of arbitrary phase. Phil Mag 1880;10:73–8. 39. Cohen L. Time-frequency analysis. Englewood Cliffs, N.J: Prentice Hall PTR, 1995. 40. Efron B. Bootstrap methods - another look at the jackknife. Ann Stat 1979;7:1–26. 41. Zygierewicz J, Blinowska KJ, Durka PJ, Szelenberger W, Niemcewicz S, Androsiuk W. High resolution study of sleep spindles. Clin Neurophysiol 1999;110:2136–47. 42. Nakamura M, Uchida S, Maehara T, et al. Sleep spindles in human prefrontal cortex: an electrocorticographic study. Neurosci Res 2003;45:419–27. 43. Silber M, Ancoli-Israel S, Bonnet M, et al. The visual scoring of sleep in adults. J Clin Sleep Med 2007;3:121–31. 44. Nicolas A, Petit D, Rompre S, Montplaisir J. Sleep spindle characteristics in healthy subjects of different age groups. Clin Neurophysiol 2001;112:521–7. 45. Kandel ER, Schwaltz JH, Jessell TM. Principles of neural science, 4th ed. New York: McGraw-Hill Health Professions Division, 2000. 46. Preud’homme XA, Lanquart JP, Mendlewicz J, Linkowski P. Characteristics of spontaneous sleep with varying NREMS Episodes in healthy men: implication for delta activity homeostasis. Sleep 2000;23:193–203. 47. Le Bon O, Staner L, Rivelli SK, Hoffmann G, Pelc I, Linkowski P. Correlations using the NREM-REM sleep cycle frequency support distinct regulation mechanisms for REM and NREM sleep. J Appl Physiol 2002;93:141–6. 48. McCormick L, Nielsen T, Nicolas A, Ptito M, Montplaisir J. Topographical distribution of spindles and K-complexes in normal subjects. Sleep 1997;20:939–41. 49. Anderer P, Klösch G, Gruber G, et al. Low-resolution brain electromagnetic tomography revealed simultaneously active frontal and parietal sleep spindle sources in the human cortex. Neuroscience 2001;103:581–92. 50. Schabus M, Dang-Vu TT, Albouy G, et al. Hemodynamic cerebral correlates of sleep spindles during human non-rapid eye movement sleep. Proc Natl Acad Sci U S A 2007;104:13164–9. 51. Manshanden I, De Munck JC, Simon NR, Lopes da Silva FH. Source localization of MEG sleep spindles and the relation to sources of alpha band rhythms. Clin Neurophysiol 2002;113:1937–47. 52. Ishii R, Dziewas R, Chau W, et al. Current source density distribution of sleep spindles in humans as found by synthetic aperture magnetometry. Neurosci Lett 2003;340:25–8. 53. Steriade M, Gloor P, Llinas R, Lopes da Silva F, Mesulam M. Basic mechanisms of cerebral rhythmic activities. Electroenceph Clin Neurophysiol 1990;76:481–508. 54. Frauscher B, von Ellenrieder N, Dubeau F, Gotman J. Scalp spindles are associated with widespread intracranial activity with unexpectedly low synchrony. Neuroimage 2015;105:1–12. 55. Karol EA, Pandya DN. The distribution of the corpus callosum in the rhesus monkey. Brain 1971;94:471–86. 56. Jones EG, Hendry SHC. Distribution of callosal fibers around the hand representations in monkey somatic sensory cortex. Neurosci Lett 1980;19:167–72. 57. Iwamura Y. Bilateral receptive field neurons and callosal connections in the somatosensory cortex. Philos Trans R Soc Lond B Biol Sci 2000;355:267–73. 58. Morris R, Pandya DN, Petrides M. Fiber system linking the middorsolateral frontal cortex with the retrosplenial/presubicular region in the rhesus monkey. J Comp Neurol 1999;407:183–92. 59. Molinari M, Minciacchi D, Bentivoglio M, Macchi G. Efferent fibers from the motor cortex terminate bilaterally in the thalamus of rats and cats. Exp Brain Res 1985;57:305–12.

1734

Sleep Spindles in Primate Cerebral Cortex—Takeuchi et al.

60. Preuss TM, Goldman-Rakic PS. Crossed corticothalamic and thalamocortical connections of macaque prefrontal cortex. J Comp Neurol 1987;257:269–81. 61. Rinvik E. Thalamic commissural connections in the cat. Neurosci Lett 1984;44:311–6. 62. Raos V, Bentivoglio M. Crosstalk between the two sides of the thalamus through the reticular nucleus: a retrograde and anterograde tracing study in the rat. J Comp Neurol 1993;332:145–54. 63. De Gennaro L, Ferrara M. Sleep spindles: an overview. Sleep Med Rev 2003;7:423–40. 64. Mascetti L, Foret A, Shaffii-Le Bourdiec A, et al. Spontaneous neural activity during human non-rapid eye movement sleep. Prog Brain Res 2011;193:111. 65. Dehghani N, Cash SS, Rossetti AO, Chen CC, Halgren E. Magnetoencephalography demonstrates multiple asynchronous generators during human sleep spindles. J Neurophysiol 2010;104:179–88. 66. Zikopoulos B, Barbas H. Prefrontal projections to the thalamic reticular nucleus form a unique circuit for attentional mechanisms. J Neurosci 2006;26:7348. 67. Guillery R, Feig S, Lozsadi D. Paying attention to the thalamic reticular nucleus. Trends Neurosci 1998;21:28–32. 68. Jones EG, Wise SP, Coulter JD. Differential thalamic relationships of sensory-motor and parietal cortical fields in monkeys. J Comp Neurol 1979;183:833–81. 69. Schell GR, Strick PL. The origin of thalamic inputs to the arcuate premotor and supplementary motor areas. J Neurosci 1984;4:539–60.

SLEEP, Vol. 39, No. 9, 2016

70. Rouiller EM, Welker E. A comparative analysis of the morphology of corticothalamic projections in mammals. Brain Res Bull 2000;53:727–41. 71. Colrain IM. The K-complex: a 7-decade history. Sleep 2005;28:255–73. 72. Amzica F, Steriade M. The K-complex: its slow (< 1-Hz) rhythmicity and relation to delta waves. Neurology 1997;49:952–9. 73. Happe S, Anderer P, Gruber G, Klosch G, Saletu B, Zeitlhofer J. Scalp topography of the spontaneous K-complex and of delta-waves in human sleep. Brain Topogr 2002;15:43–9. 74. Le Van Quyen M, Staba R, Bragin A, et al. Large-scale microelectrode recordings of high-frequency gamma oscillations in human cortex during sleep. J Neurosci 2010;30:7770. 75. Steriade M. Grouping of brain rhythms in corticothalamic systems. Neuroscience 2006;137:1087–106.

SUBMISSION & CORRESPONDENCE INFORMATION Submitted for publication October, 2015 Submitted in final revised form May, 2016 Accepted for publication May, 2016 Address correspondence to: Toru Tsujimoto, Department of Neurology, Nagahama City Hospital, Nagahama, 526-8580, Japan; Tel: +81-749-682300; Fax: +81-749-65-1259; Email; [email protected]

DISCLOSURE STATEMENT This was not an industry supported study. The authors have indicated no financial conflicts of interest.

1735

Sleep Spindles in Primate Cerebral Cortex—Takeuchi et al.