Rotating waves during human sleep spindles organize global patterns of activity that repeat precisely through the night

1Memories are stored in distributed networks across the cortex. In the two-stage model of memory consolidation, memories are integrated in the hippocampus and then linked in the neocortex for long-term storage, where information represented in visual, auditory, somatosensory, or cognitive regions must be bound into a coherent whole. It is well established that sleep oscillations actively contribute to this process: during stage 2 sleep spindles, the thalamus generates a rhythmic activity pattern that becomes widespread through large-scale thalamocortical loops, and spindles are critical to sleep-dependent memory consolidation. Long-range connections in cortex result primarily from excitatory pyramidal cells, but precisely how sleep oscillations aid strengthening of these excitatory connections between distributed cortical networks through spike-time dependent plasticity (STDP) remains unclear, particularly in the presence of long axonal conduction delays. Here, we identify a global activity pattern repeatedly observed during sleep spindle oscillations in human neocortex that could serve this role.

2We study intracranial electrocorticogram (ECoG) recordings of five clinical patients in stage 2 sleep and apply recently developed computational methods to classify spatiotemporal dynamics at the level of individual oscillation cycles. ECoG arrays were implanted in subjects undergoing evaluation for resective surgery of epileptogenic cortex (Figure 1A, left). Over several days of recording, these subjects exhibit long periods without major epileptic events. During that time, subjects express a relatively normal sleep architecture, with well-defined sleep oscillations. Stage 2 sleep epochs were then manually identified by an expert rater, and sleep spindles recorded on the ECoG were isolated using automated techniques. These spindles appeared physiologically normal and well-isolated from background noise (Figure 1A, right and Figure 1—figure supplement 1). Our algorithmic approach classifies spatiotemporal patterns as expanding waves, defined as a significant linear increase in phase offset with distance from a point source (Figure 1—figure supplement 2; see Materials and methods – Spatiotemporal dynamics), or rotating waves, defined as a significant increase in phase offset with rotation about a wave center (Figure 1—figure supplement 3). In 41,860 spindle oscillation cycles tested across subjects, a large proportion (50.8%) was classified as rotating waves, along with a smaller subset (15.6%) as expanding. After inspecting these results, we observed further that the rotating waves exhibited a clear bias towards travel in the temporal parietal frontal (TPF) direction (69.5%, p<10−10, one-tail binomial test against equal occurrence, 14,796 TPF cycles, 21,272 total; for each individual subject p<10−3, see Figure 1—source data 1 for individual wave totals) (Figure 1B and Video 1). Propagation speed distributions peaked between 2–5 m/s (Figure 1C), varying within a narrow range from the 20th to the 80th percentiles (3–9 m/s for the full distribution; 4–10, 3–8, 3–10, 2–3, and 4–13 m/s for individual subjects, respectively), within the range of conduction speeds for the short and long white matter association fibers. Further, this rotating TPF organization occurred consistently across subjects and implantation hemispheres (Figure 1D; see also Figure 1—figure supplements 5–7).

3Spike-time dependent plasticity is a well-studied mechanism for regulating synaptic strengths that depends on the relative timing of presynaptic inputs and postsynaptic spikes, but for establishing large-scale neural assemblies during sleep oscillations through synaptic plasticity, axonal conduction delays pose a specific problem. For example, cortical white matter association fibers have conduction delays up to 50 milliseconds across the cortex (Figure 2A, left).

4It is well established that spindles cause pyramidal cells and interneurons in cortex to fire preferentially at the peak of the surface-positive (depth-negative) LFP oscillation, both in intracellular and extracellular recordings. If cortical spindles were perfectly synchronized, spikes emitted during one cycle of the spindle oscillation would arrive at their post-synaptic targets with this temporal delay, leading to a pairing within the window for persistent long-term depression (LTD) that would progressively weaken long-range connections (Figure 2A, right). If, however, spindles are self-organized into large-scale wave-like activity patterns, with phase speeds matching those of the underlying fiber networks and stereotyped, precisely repeating trajectories (Figure 2B, left), then EPSPs caused by spikes traveling along pyramidal axons to distant regions in the cortex would align with the local burst of population activity (Figure 2B, right), creating the conditions necessary for synaptic strengthening to occur.

5Next, we wanted to understand whether these population activity patterns repeat with the temporal precision required for strengthening of large-scale assemblies. We defined the correlation magnitude over phase values on the electrode array between individual oscillation cycles to be a pairwise similarity index (see Materials and methods), in order to detect similar spatiotemporal patterns across oscillation cycles (Video 2). By calculating this metric over all cycle pairs in different wave classes (all cycles, expanding, rotational), we can directly compare the temporal precision mediated by each type. The cumulative distribution function (CDF) of similarity indices among identified rotational waves is highly shifted to the right (Figure 3A, black) compared to the CDF for all cycles (5 subjects, 42/54 sleep epochs, 77.8% significant, one-tailed two-sample Kolmogorov-Smirnov test,, Bonferroni correction), indicating higher intra-class similarity between these cycles than for other wave types. Note that this is not simply a consequence of the rotational phase pattern itself, as expanding waves emanating from a consistent point source could certainly exhibit higher intra-class similarity than rotational waves with a varying center. Further, the median similarity index consistently increases in individual subjects when rotational waves of progressively increasing strength are considered (Figure 3B). This indicates the observed rotating waves strongly modulate temporal precision in repeated patterns of population activity. To be specific, by utilizing the average temporal frequency for these spindle oscillations (13.5 Hz), we can estimate that in two cycles whose similarity index falls into the highest bin in Figure 3C (0.9–1.0), 50% of electrodes will experience an alignment of the spatiotemporal activity pattern within a 5 millisecond temporal window. Recent experiments have shown a tight temporal link between field potentials and synaptic currents (both EPSPs and IPSPs; ). By detecting these precisely recurring activity patterns in ECoG recordings, we can infer that distributed networks composed of local excitatory and inhibitory groups, whose firing is modulated by thalamocortical fibers during the sleep spindle, are repeatedly activated with a millisecond accuracy that is well within the temporal precision required for STDP.

6If this millisecond precision in fact mediates formation and maintenance of corticocortical assemblies, we would then expect spiking associated with these synaptic currents to drive increased reverberation throughout the night, as excitatory connections between local groups of pyramidal cells and interneurons are strengthened and in turn promote more replay of the expressed activity pattern. High gamma-band power (HGP, 80–120 Hz), a reliable electrophysiological correlate of spiking activity (;;; similar in nature to the 'broadband power shift' described in ), consistently increases around spindles (Figure 4A). Further, HGP is modulated by spindle phase (Figure 4B), increasing towards the surface-positive (depth-negative) peak, consistent with previous animal and human recordings. Finally, by studying repeats of rotating waves over 2.5 hr of continuous sleep recording in Subject 1, we observe a preliminary indication of increased reverberation consistent with strengthening of distributed excitatory networks: similarity in the next identified rotating wave is highly predictive of the number of strong reverberations throughout the night (black dots, Figure 4C). Randomizing the relationship between the next rotating wave and the rest of the sleep recording eliminates this effect (shuffling control, Figure 4C), and such increased reverberation is not observed for expanding waves under similar conditions (Figure 4—figure supplement 1). Similar observations are consistent across subjects (Figure 4—figure supplement 2). These results support the hypothesis that precisely repeating rotating waves may enable strengthening of large-scale corticocortical assemblies throughout the night.

7Early animal sleep spindle studies, using up to 8 electrodes in a linear array, in addition to preliminary EEG evidence in the human, proposed that spindles involve global synchronization of cortical circuits, raising the possibility that this sleep oscillation places neocortex into a specialized state for consolidation of long-term memories. In recent years, several studies have reported a mixture of 'local' and 'global' spindles using amplitude-duration thresholding approaches. By carefully studying the phase information in the spindle frequency band recorded on large-scale ECoG arrays, we have uncovered that a substantial number of spindle oscillation cycles are organized into global, hemisphere-spanning patterns of rotating and expanding waves (Figure 1—figure supplement 7). These patterns most likely represent the characteristic spatiotemporal organization of the 'global' spindles observed in Andrillon et al. (≥40% involvement, cf. their Figure 5C), with more localized patterns left unclassified by our detection approach (Figure 1—figure supplements 9 and 12). These global patterns are likely established through widespread thalamocortical loops, placing the cortex into a state of large-scale coherence, shaped into rotating and expanding waves through corticocortical white matter connections with axonal conduction speeds consistent with the observed propagation speeds (Figure 1C). Future computational modeling work will address in detail the role of thalamocortical, corticocortical, and corticothalamic connections in generating the spatiotemporal activity patterns reported here.

8Spindles have recently been specifically and causally implicated in the sleep-dependent consolidation of long-term memories. While some memories integrate content from single sensory modalities, requiring consolidation in only single cortical regions (such as motor cortex, ), many memories integrate multimodal sensory and cognitive information, and require 'global' integration of distributed networks across the cortex. In this work, we have identified a novel mechanism by which this process could occur: the stereotyped activity patterns reported here may enable STDP to establish large-scale neuronal assemblies at scales where axonal conduction delays are long relative to the oscillation cycle, and repeat many times throughout sleep with millisecond accuracy. While the schema illustrated in Figure 2 is a highly simplified view of the microscale interactions between long-range excitatory projections and local networks during spindle oscillations, computational and theoretical studies have previously obtained a detailed understanding of STDP dynamics with neurons receiving sequenced, bursting, and oscillating inputs. This theoretical understanding of the interplay between STDP and population activity can allow in future work a precise account of how microscale synaptic interactions are shaped by global oscillation patterns, and how variability in these patterns (e.g. variation in wave speed, Figure 1C) will affect this mechanism. Taken together, these results provide insight into how distributed information stored across cortical regions may be bound into a coherent, integrated, but specific memory through spike-time dependent synaptic plasticity.

9Patients with longstanding pharmacologically resistant complex seizures gave fully informed consent according to NIH guidelines as monitored by the local Institutional Review Board (Massachusetts General Hospital). Electrocorticogram (ECoG) recordings during natural sleep were made over the course of clinical monitoring for spontaneous seizures. Electrode placement was determined solely by clinical criteria, with electrode grids usually spanning the Sylvian fissure and multiple lobes of the cerebral cortex (frontal, parietal, and temporal). Patients were informed that participation in the research would not alter their clinical treatment in any way, and that they may withdraw their consent at any time without jeopardizing clinical care.

10ECoG contacts (Ad-Tech Medical Instrument Corp., Racine, WI) were 3 mm platinum-iridium (90% platinum) discs arranged in a two-dimensional grid (8 rows and 8 columns, Subjects 1, 3, and 4; 8 rows and 12 columns, Subject 2; 8 rows and 6 columns, Subject 5) implanted semi-chronically on the pial surface in an effort to localize the seizure origin. Within the grid, electrodes were spaced 10 mm apart. In some patients, linear ECoG arrays provided additional spatial coverage; application of our multichannel detection approach, however, focused on the two-dimensional electrode grid. One strip of electrodes positioned over the pial surface and facing the skull served as the reference during the recordings; results were additionally verified using an average reference. Note that due to reference artifacts, an average reference was employed for the recordings in Subject 4. We note as well that the temporal extent of the waves, over tens of milliseconds, makes electrophysiological artifacts such as volume conduction an unlikely explanation for the observations reported here. Recordings were performed with clinical EEG monitoring equipment (XLTEK, Natus Medical Inc., Pleasanton, CA) and sampled at 500 or 512 Hz.

11Post-implantation electrode localization utilized coregistration of preoperative magnetic resonance imaging (MRI) with postoperative computed tomography (CT), as described by Dykstra et al.. Cortical surfaces were computed with FreeSurfer. To account for the misalignment between the MRI and CT due to the craniotomy, the locations of the grid electrodes were projected onto the cortical surface. Geodesic electrode distances, which take into account the folded geometry of the cortical surface and were used in some calculations (e.g. estimation of spatial correlation values), were estimated using a shortest paths approach on the cortical surface mesh.

12During the monitoring period for spontaneous seizures, the subjects slept in the clinical environment and expressed relatively normal sleep patterns. ECoG recordings that did not have a seizure in the preceding or following 12 hr were scored visually by an expert rater following the standard sleep stage classification. For each patient, we obtained from 15 to 101.5 min of NREM stage 2 sleep, when the spindles are most prevalent. Individual sleep spindles were then detected during stage 2 sleep using one of several complementary methods, either based on amplitude-duration thresholding or a similar wavelet-based approach with additional verification steps. The number of spindles detected was in agreement with previous reports of spindle density.

13The results were additionally verified using a novel approach quantifying the signal-to-noise ratio (SNR) of power in the bandpass (9–18 Hz, 8th-order Butterworth filter) versus the bandstop (1–100 Hz bandpass, with 9–18 Hz bandstop) signal. In this approach, the SNR metric is calculated on short (500 ms) sliding windows in each channel. When the SNR metric reaches 0 dB, the signal and noise power are at parity, corresponding to a sharp, narrowband epoch in the recording. Picking a constant SNR threshold (5 dB) corresponds roughly to the constant false alarm rate (CFAR) technique in radar. This approach yields a conservative but approximately amplitude-invariant method for detecting arbitrary narrowband epochs in multichannel data.

14Following spindle detection, we made a verification analysis by calculating the average power spectral density (PSD) over isolated spindle epochs in each subject. Data were initially filtered to remove line noise artifacts, and PSDs were then calculated in 1 s intervals during the spindle and matched non-spindle epochs. PSDs for individual channel and spindle epochs were concatenated into a large array and averaged in each case. A clear peak in the 11–15 Hz frequency band for the spindle epochs can be seen, while no peak is observed in the matched non-spindle epochs (Figure 1—figure supplement 1A). Divisive normalization is calculated by dividing the power at each frequency in the spindle epochs by the power in the matched non-spindle epochs, and expressing the result in dB (Figure 1—figure supplement 1A, inset). If the divisive normalization over an epoch of stage 2 sleep reached 5 dB, then the spindles were taken to be well-isolated and possessing the spectral characteristics necessary for an accurate phase representation, and were then included in further analysis. Based on this calculation, 54 individual epochs of stage 2 sleep, varying from 30 s to 35 min in duration, were selected in five clinical subjects.

15Temporal filtering of stage 2 sleep recordings was carried out with an 8th-order digital Butterworth bandpass filter (9–18 Hz), forward-reverse in time to prevent phase distortion (see MATLAB function filtfilt). All results were checked with multiple cutoff frequencies to ensure against parametric sensitivity.

16To assess spatial correlation during spindle oscillations as a function of distance in the cortex, we adapted standard methods with a Monte Carlo implementation more suited for sampling correlations on two-dimensional electrode arrays. To calculate this metric, one electrode is first selected at random, and a second is then selected from the set of electrodes within a binned distance from the first. The temporal correlation between these electrode pairs is then computed in the bandpass timeseries between the start and end points of the spindle. This process is repeated for a given number of iterations at each distance bin, and the average spatial correlation is computed as the mean of the correlation values for the individual epochs (Figure 1—figure supplement 1B, black). The spatial correlation values were computed for non-spindle epochs matched to the temporal extent of the individual tested spindle epochs (Figure 1—figure supplement 1B, red). The average spatial correlation values are elevated during spindle oscillations with respect to the matched non-spindle periods of stage 2 sleep, indicating that a coherent, large-scale increase in global activity occurs during spindles, in agreement with previous studies.

17To study spatiotemporal dynamics in these neural recordings, we adapted our previously introduced method for detecting arbitrarily shaped traveling waves to multisite ECoG arrays (Figure 1—figure supplements 2 and 3). This approach allowed us to characterize and classify the spatiotemporal dynamics during thousands of episodes of spindling activity in many hours of sleep recordings. The method proceeds in three steps: (1) analytic signal representation for characterization of instantaneous signal characteristics at each electrode, (2) center localization at each individual oscillation cycle, and (3) quantification of the spatiotemporal pattern in each oscillation cycle as a function of distance from (or rotation about) the isolated center point. In the following, we describe in detail the method for isolating expanding and rotating waves in multichannel data.

18To estimate instantaneous signal characteristics, we employ the well-known analytic signal representation. This approach entails transforming a real-valued timeseries into a complex phasor, whose modulus (length) and argument (angle) in the complex plane represent the signal instantaneous amplitude and phase, respectively. Specifically, if is a real-valued, narrowband timeseries at a point, where and denote the number of rows and columns, and is the sample number, then its analytic signal representation is

19where is the complex unit and denotes the Hilbert transform of a signal. The instantaneous phase of is then the argument at each point in this complex sequence

20and instantaneous amplitude is the modulus. At several points in the analysis, results were confirmed with an FIR implementation of the Hilbert transform, in addition to the standard FFT-based approach. We evaluated phase values at a set of time points in each spindle near the positive oscillation peaks. The phase fields were then smoothed using a robust approach for center localization to reduce noise and interpolate values from missing electrodes; note that the smoothed values were not used in the calculations for detecting expanding and rotational wave patterns.

21These phase values are then used to capture spatiotemporal dynamics in the multichannel data. To isolate putative expanding or rotating wave centers in each oscillation cycle, we first assess the spatial gradient of phase

22with. For the spatial gradient, derivatives are taken across the two dimensions of space and are approximated by the appropriate forward and centered finite differences. As in previous work, phase derivatives were implemented as multiplications in the complex plane.

23To detect expanding waves, we assess the divergence of the phase gradient field

24and define the putative wave source to be that point which satisfies the arg max over space in each cycle

25where. This step allows us to find the source for a possible expanding wave in each cycle (step 2, Figure 1—figure supplement 2), about which the phase field is then evaluated to quantify the evidence for an expanding wave spatiotemporal organization (step 3, Figure 1—figure supplement 2). For this next step, we calculate the circular-linear correlation coefficient between signal phase and radial distance from the source point in the original, unsmoothed phase field

26where represents the Pearson correlation between the cosine of the circular variable and the linear variable, between the sine of and the variable, and between the cosine and sine of. This approach allows us to quantify the strength of the spatiotemporal pattern of activity on the array in a single number, which is then compared to the value produced by repeating the calculations many times under random shuffling of the data (blue dotted line, Figure 1—figure supplement 4A and Materials and methods – Shuffling Controls).

27Analogous to the above case, we start by assessing the curl of the phase gradient field

28and defining the putative center to be that point which satisfies the arg max over space

29This center point then defines an anchor about which we can pass into a polar coordinate system, describing the distance and rotation angle about that point (step 2, Figure 1—figure supplement 3). With the putative rotation center isolated in each oscillation cycle, we then proceed to calculate the circular-circular correlation coefficient between signal phase and rotation angle in the original, unsmoothed phase field

30where overbar indicates circular mean

31Similar to the previous case, this number quantifies the evidence for a rotational wave organization on the array, which is then compared to the value derived from a random-shuffling permutation test (red dotted line, Figure 1—figure supplement 4A and Materials and methods – Shuffling Controls). For each case, additional control analyses with simulated rotating waves embedded in noise were used to verify the robustness of our approach (Figure 1—figure supplement 11). Finally, in the case that both expanding and rotational elements are detected, the pattern is classified as rotational, because sub-patterns of rotational waves tend to be detected as expanding elements (verified in Figure 1—figure supplement 5; see also Materials and methods – Average vector field controls).

32To quantify the level of spatiotemporal phase flow expected in the data by chance, we implemented a shuffling procedure to establish a permutation-based threshold for both the expanding and rotating wave measures. To do this, we shuffled the phase values in each oscillation cycle randomly across space a number of times (100 or 1000 times in initial tests, then reduced to 25 without changing results), repeating each time the same calculation as for the un-shuffled data. The 99th percentile of the resulting distribution then determines a threshold above which the value for the correlation metric (either for expanding or rotational waves, considered separately) exceeds chance, with the spatial autocorrelation erased.

33A possible confound resulting from this shuffling procedure is that the data intrinsically possess some spatial autocorrelation (Figure 1—figure supplement 1B), which is ignored by the so-constructed permutation test. To address this point in the context of rotational wave detection, we conceived an additional permutation test control. In this second control, we considered the set of points at a Chebyshev (i.e. King’s chessboard) distance from the putative rotation center, where indicates the maximum distance on the electrode array from that point. We then shuffled channels at distances for all, and repeated the calculation for the rotational wave detection. The resulting permuted data have a spatial correlation function identical to that in the un-shuffled data, but with the rotational structure fully destroyed. The 99th percentile cutoff determined from this second control analysis fell within 0.01 of the originally estimated value (3% difference), validating the original shuffling permutation test employed above.

34Using simulated expanding waves of the form

35and simulated rotating waves

36where is the oscillation amplitude, is the oscillation angular frequency, is the wavenumber, is the polar wavenumber, and is a Gaussian white noise term, we verified the robustness of our detection approach under noise of varying amplitudes. Note that and are defined with respect to the wave center, left unspecified for simplicity. Oscillation amplitude was set to unity, without loss of generality, and other oscillation parameters were matched to those observed during stage 2 sleep spindles. Oscillation frequency was set to the average instantaneous frequency estimated from 702 spindles in Subject 1 (13.5 Hz), and wavenumbers were adjusted to approximate the wavelengths observed in the data. Varying systematically the level of added noise, we ran the algorithms for detecting expanding and rotating waves described above for 25 trials at each point and recorded the algorithm’s detection performance in each case (mean SEM, Figure 1—figure supplement 11C). These results illustrate the approximate invariance of our computational approach to random noise.

37In another test, we systematically varied the position of the wave center on the simulated 64 electrode array, for both expanding and rotational waves (Figure 1—figure supplement 11D). Parameters were set as above, and simulations were again run over 25 trials at each point. This test probed the sensitivity of the rotational detection approach to border effects, which is expected to be negligible at the encountered noise levels in comparison to the thresholds established by the permutation controls (dotted lines, Figure 1—figure supplement 11D).

38In a third test, we verified that the spatiotemporal patterns observed here are not due to variations in spindle frequency, which are known to occur along the rostro-caudal axis. To do this, we re-ran our analysis on one stage 2 sleep session containing 179 spindles in Subject 1, generating surrogate data as follows. Each electrode evolved in time according to its mean instantaneous frequency during the spindle, but with a randomized initial phase angle. These surrogate data thus possessed the same frequency content on average as in the original data, but with their spatial organization of phase removed. In this control, both rotating and expanding wave patterns were highly decreased (3.5% and 0.8% of cycles, respectively, compared to 64% and 14% in the original data).

39The algorithmic classification of wave patterns in individual oscillation cycles involves several steps, and we wanted to make an independent check to verify these results. To do this, we adopted a re-centered averaging approach, shifting the vector field of propagation directions from the smoothed phase fields at each oscillation cycle to the putative rotation center (red dots, Figure 1—figure supplement 5), and taking the circular mean of propagation direction at each point. Performing the calculation in this way prevents regions with noise or high phase gradient magnitude from dominating the result. The obtained vector fields for rotational TPF (Figure 1—figure supplement 5A) and expanding (Figure 1—figure supplement 5B) waves illustrate the accuracy of the algorithm and the general validity of our classification approach in separating waves into expanding and rotational groups.

40To quantify the precision of repeated spatiotemporal patterns during across spindle oscillations over several minutes of data, we calculated the circular-circular correlation between phase values in individual oscillation cycles. For two phase maps and, the circular correlation is defined as above

41This correlation value defined between individual phase maps then constitutes elements of an matrix, where is the number of isolated oscillation cycles in question:

42where all elements in the diagonal and lower triangle are not considered, without loss of generality. We then construct this matrix for all cycles in an individual wave classification (all cycles, expanding, rotational), and consider the cumulative distribution function (CDF) of values in the upper triangle (Figure 3A). To construct the CDF in the permutation case (Figure 3A, purple line), we first randomly shuffled the phase maps in all spindle cycles and then proceeded with the calculation as normal.

43A MATLAB toolbox for analysis of traveling waves and complex spatiotemporal dynamics in noisy multisite data is available as an open-source release on BitBucket: http://bitbucket.org/lylemuller/wave-matlab