1. Neuroscience
Download icon

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

  1. Lyle Muller
  2. Giovanni Piantoni
  3. Dominik Koller
  4. Sydney S Cash
  5. Eric Halgren
  6. Terrence J Sejnowski Is a corresponding author
  1. Salk Institute for Biological Studies, United States
  2. Harvard Medical School, United States
  3. University of California, San Diego, United States
Short Report
Cited
1
Views
2,138
Comments
0
Cite as: eLife 2016;5:e17267 doi: 10.7554/eLife.17267

Abstract

During sleep, the thalamus generates a characteristic pattern of transient, 11-15 Hz sleep spindle oscillations, which synchronize the cortex through large-scale thalamocortical loops. Spindles have been increasingly demonstrated to be critical for sleep-dependent consolidation of memory, but the specific neural mechanism for this process remains unclear. We show here that cortical spindles are spatiotemporally organized into circular wave-like patterns, organizing neuronal activity over tens of milliseconds, within the timescale for storing memories in large-scale networks across the cortex via spike-time dependent plasticity. These circular patterns repeat over hours of sleep with millisecond temporal precision, allowing reinforcement of the activity patterns through hundreds of reverberations. These results provide a novel mechanistic account for how global sleep oscillations and synaptic plasticity could strengthen networks distributed across the cortex to store coherent and integrated memories.

https://doi.org/10.7554/eLife.17267.001

eLife digest

When you wake up in the morning after a good night's sleep you feel refreshed. You can also think more clearly because your memory has been re-organized, a process called memory consolidation. The problem that the brain has to solve during sleep is how to integrate memories of experiences that happened during the day with old memories, without losing the older memories.

Scientists know that waves of electrical activity, referred to as spindles, help to consolidate and integrate memories during sleep. Spindles are active in the cerebral cortex, the part of your brain used for thinking, in the time between dream sleep and deep sleep. Yet it is not known exactly how these bursting patterns of electrical activity help to strengthen memories.

Now, Muller et al. explored how the spindles could strengthen and connect parts of memories stored in distant parts of the brain. First, a computer algorithm analyzed electrical recordings of brain activity taken while five patients with epilepsy slept. The patients were being monitored to help with their seizures, and the recordings showed that spindles do not occur at the same time throughout the cortex as previously thought. Instead, the spindle is a wave that begins in portion of the cortex near the ear, spirals through the cortex toward the top of back of the head and then on to the forehead area before circling back.

These repeated circular waves of electrical activity strengthen connections between brain cells in distant parts of the brain. For example, these waves may help strengthen connections between the cells of the cortex that separately store memories of the sound, sight and feel of an event during the day, whether that’s being bitten by a dog or talking with a friend. Next, Muller et al. plan to develop computer models of the spindles and verify whether their models make accurate predictions by studying spindles in sleeping mice and rats.

https://doi.org/10.7554/eLife.17267.002

Introduction

Memories are stored in distributed networks across the cortex. In the two-stage model of memory consolidation (McClelland et al., 1995; Rasch and Born, 2007), 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 (Wheeler et al., 2000; Horner et al., 2015). 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 (Contreras et al., 1996), and spindles are critical to sleep-dependent memory consolidation (Gais et al., 2002; Mednick et al., 2013). Long-range connections in cortex result primarily from excitatory pyramidal cells (Sholl, 1956; Schüz et al., 2002), 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 (Lubenov and Siapas, 2008). Here, we identify a global activity pattern repeatedly observed during sleep spindle oscillations in human neocortex that could serve this role.

We study intracranial electrocorticogram (ECoG) recordings of five clinical patients in stage 2 sleep and apply recently developed computational methods (Muller et al., 2014) 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 (Hagler et al., 2016). 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 (Girard et al., 2001) and long (Schüz et al., 2002; Swadlow and Waxman, 2012) white matter association fibers. Further, this rotating TPF organization occurred consistently across subjects and implantation hemispheres (Figure 1D; see also Figure 1—figure supplements 57).

Figure 1 with 12 supplements see all
Rotating waves during spindles.

(A) Electrode placement for subject 1 (left), with a stereotypical spindling epoch observed on the array (right). The right panel depicts the average over channels (black) together with the individual channels (gray). (B) When visualized on the cortex, individual spindle cycles are often organized as rotating waves traveling from temporal (+0 ms, top) to parietal (+20 ms, middle) to frontal (+40 ms, bottom) lobes. (C) Phase speed distributions across subjects. Plotted is the kernel smoothing density estimate for individual subjects (gray dotted lines) and for the full distribution (black line). (D) The field of propagation directions, aligned on the putative rotation center and averaged across oscillation cycles and across subjects, shows a consistent flow in the temporal parietal frontal (TPF) direction. The center point is marked in red.

https://doi.org/10.7554/eLife.17267.003
Video 1
Rotating waves over five spindle oscillation cycles.

Normalized activity for bandpass filtered timeseries is plotted in falsecolor at electrode positions on the cortical surface of Subject 1. The cortical electrode marked with a red dot (bottom) corresponds to the black timecourse in the inset (top). The other ECoG channels are plotted in gray. The time period visualized corresponds to approximately 300 milliseconds, or five cycles of the spindle oscillation. Note that no spatial smoothing is applied in these data.

https://doi.org/10.7554/eLife.17267.017

Spike-time dependent plasticity is a well-studied mechanism for regulating synaptic strengths that depends on the relative timing of presynaptic inputs and postsynaptic spikes (Markram et al., 1997; Bi and Poo, 1998), but for establishing large-scale neural assemblies during sleep oscillations through synaptic plasticity, axonal conduction delays pose a specific problem (Lubenov and Siapas, 2008). For example, cortical white matter association fibers have conduction delays up to 50 milliseconds across the cortex (Figure 2A, left) (Girard et al., 2001; Schüz et al., 2002; Swadlow and Waxman, 2012).

Schematic of spindles and axonal delays.

(A) Spikes emitted from region A will arrive at B with a temporal delay of 20 milliseconds (left). If spindle oscillations were perfectly synchronized across the cortex, EPSPs from region A would occur after the spikes in region B, within the window for long-term depression (right). (B) In contrast, if spindles are spatiotemporally organized with stereotyped trajectories (left), then EPSPs from region A would align with population spiking in region B, allowing for synaptic strengthening to occur.

https://doi.org/10.7554/eLife.17267.018

It 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 (Contreras and Steriade, 1995, 1996; Kandel and Buzsáki, 1997) and extracellular (Peyrache et al., 2011) 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.

Next, 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, α=0.01, 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; Haider et al., 2016). 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.

Video 2
Rotating waves with high spatiotemporal similarity.

Two rotating waves with high phase similarity on the ECoG array, separated by 5.62 min of stage 2 sleep. Bandpass filtered timeseries are normalized to their maximum within the interval and plotted in falsecolor (bottom panels). Activity for each channel is plotted as a function of time (top panels), with an indication of temporal progression (red dotted line).

https://doi.org/10.7554/eLife.17267.019
Phase pattern analysis.

(A) Rotating waves exhibit higher intra-class similarity. Cumulative distribution functions (CDFs) for shuffled data (purple), expanding waves (red), all cycles (blue), and rotating waves (black) are given for an example 15 min epoch of stage 2 sleep (subject 5). (B) Spindle cycles exhibiting stronger rotating patterns also express greater intra-class similarity. Gray lines indicate the median similarity index (ordinate) for the population of oscillation cycles expressing rotational waves above a threshold strength (abscissa), averaged over individual sleep epochs. Red dots and error bars indicate the median and median absolute deviation for the full distribution, respectively. (C) Spindle cycles exhibiting high similarity index are temporally precise. The distribution of phase difference at each electrode across spindle cycles is given as a function of the similarity index (indicated by colors, inset). By utilizing the mean spindle oscillation frequency (13.5 Hz), the midspread (interquartile range) of each distribution is given in units of time (inset).

https://doi.org/10.7554/eLife.17267.020

If 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 (Ray et al., 2008; Ray and Maunsell, 2011; Ray, 2015; similar in nature to the 'broadband power shift' described in Manning et al., 2009), 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 (Peyrache et al., 2011) and human (Andrillon et al., 2011) 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.

Figure 4 with 3 supplements see all
Spiking activity and increased reverberation.

(A) High gamma-band power (HGP) consistently increases around spindle onset. Plotted are the normalized amplitude envelopes for spindles (black) and HGP (red), averaged over 186 spindles in Subject 1. (B) HGP is modulated by spindle phase. Plotted is the mean high gamma-band power at each phase of the spindle oscillation (20 bins), for varying amplitudes of the spindle oscillation (see colorbar), each normalized by the mean HGP in matched non-spindle epochs. (C) Strength of the first repeat predicts the number of strong reverberations. The number of similar rotational patterns (above similarity index 0.7) following a spindle oscillation cycle is given as a function of the next cycle’s similarity index (black dots, mean + SEM) over 2.5 hr of sleep in Subject 1. Error bars are obscured by markers. Light blue lines indicate results from a shuffling permutation test (10 iterations).

https://doi.org/10.7554/eLife.17267.021

Early animal sleep spindle studies, using up to 8 electrodes in a linear array (Andersen et al., 1967; Kim et al., 1995; Contreras et al., 1996, 1997), in addition to preliminary EEG evidence in the human (Achermann and Borbély, 1998), 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 (Nir et al., 2011; Andrillon et al., 2011). 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. (Andrillon et al., 2011) (≥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 (Contreras et al., 1996), 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.

Spindles have recently been specifically and causally implicated in the sleep-dependent consolidation of long-term memories (Mednick et al., 2013; Hennies et al., 2016). While some memories integrate content from single sensory modalities, requiring consolidation in only single cortical regions (such as motor cortex, Khazipov et al., 2004), many memories integrate multimodal sensory and cognitive information (Gibson and Maunsell, 1997), 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 (Fries, 2005; Lubenov and Siapas, 2008), 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 (Rao and Sejnowski, 2001, 2003), bursting (Song et al., 2000), and oscillating inputs (Muller et al., 2011; Luz and Shamir, 2016). 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.

Materials and methods

Subjects

Patients 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.

Electrodes

ECoG 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.

Electrode localization

Post-implantation electrode localization utilized coregistration of preoperative magnetic resonance imaging (MRI) with postoperative computed tomography (CT), as described by Dykstra et al. (Dykstra et al., 2012). Cortical surfaces were computed with FreeSurfer (Dale et al., 1999; Fischl et al., 1999). 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 (Dykstra et al., 2012). 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.

Sleep spindle detection

During 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 (Silber et al., 2007). 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 (Gais et al., 2002; Warby et al., 2014) or a similar wavelet-based approach with additional verification steps (Hagler et al., 2016). The number of spindles detected was in agreement with previous reports of spindle density (Gais et al., 2002; Warby et al., 2014).

The 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.

Power spectral density analysis

Following 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.

Temporal filtering

Temporal 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.

Spatial correlation analysis

To assess spatial correlation during spindle oscillations as a function of distance in the cortex, we adapted standard methods (Destexhe et al., 1999) 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 di 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 Nk at each distance bin di, 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 (Destexhe et al., 1999).

Spatiotemporal dynamics

To study spatiotemporal dynamics in these neural recordings, we adapted our previously introduced method for detecting arbitrarily shaped traveling waves (Muller et al., 2014) 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.

To 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 vx,y,t is a real-valued, narrowband timeseries at a point (x,y),x[1,Nc],y[1,Nr], where Nc and Nr denote the number of rows and columns, and t[1,Nt] is the sample number, then its analytic signal representation is

(1) Vx,y,t=vx,y,t+iv^x,y,t

where i is the complex unit and f^ denotes the Hilbert transform of a signal f. The instantaneous phase of vx,y,t is then the argument at each point in this complex sequence

(2) ϕx,y,t=Arg(Vx,y,t),

and 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 (Marple, 1999). We evaluated phase values at a set of time points 𝒯={t1,t2,,tK} in each spindle near the positive oscillation peaks. The phase fields were then smoothed using a robust approach (Garcia, 2010) 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.

These 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

(3) gx,y,tj-ϕx,y,tj

with tj𝒯. 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 (Feldman, 2011; Muller et al., 2014).

Detection of expanding waves

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

(4) dx,y,tj=gx,y,tj,

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

(5) 𝒮(x,y;tj)=argmaxx,ydx,y,tj,

where argmaxf(a,b){a,b|p,q:f(p,q)f(a,b)}. 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 ρϕ,δ (Jammalakadaka and Sengupta, 2001Berens, 2009) between signal phase ϕ and radial distance δ from the source point in the original, unsmoothed phase field

(6) ρϕ,δ=rcδ2+rsδ22rcδrsδrcs1rcs2,

where rcδ represents the Pearson correlation between the cosine of the circular variable ϕ and the linear variable δ, rsδ between the sine of ϕ and the variable δ, and rcs 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).

Detection of rotating waves

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

(7) cx,y,tj=×gx,y,tj,

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

(8) 𝒞(x,y;tj)=argmaxx,y||cx,y,tj||.

This 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 θ (Fisher, 1993; Berens, 2009) in the original, unsmoothed phase field

(9) ρϕ,θ=-xysin(ϕxy-ϕ¯)sin(θxy-θ¯)xysin2(ϕxy-ϕ¯)sin2(θxy-θ¯),

where overbar indicates circular mean

(10) ϕ¯=Arg[xyeiϕxy].

Similar 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).

Shuffling controls

To 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.

A 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 di[1,dm] from the putative rotation center, where dm indicates the maximum distance on the electrode array from that point. We then shuffled channels at distances di for all i, 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.

Simulated data controls

Using simulated expanding waves of the form

(11) f(t,δ)=Aei(ωt-κδ)+ση(t),

and simulated rotating waves

(12) f(t,θ)=Aei(ωt-γθ)+ση(t),

where A is the oscillation amplitude, ω is the oscillation angular frequency, κ is the wavenumber, γ is the polar wavenumber, and η(t) 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.

In 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).

In 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 (Peter-Derex et al., 2012). 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).

Average vector field controls

The 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 (Fisher, 1993; Berens, 2009) 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.

Phase map correlation analysis

To 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 αx,y and βx,y, the circular correlation is defined as above (Fisher, 1993; Berens, 2009)

(13) ρα,β=xysin(αxy-α¯)sin(βxy-β¯)xysin2(αxy-α¯)sin2(βxy-β¯).

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

Cij=[1ρ1,2ρ1,3ρ1,4ρ1,Mρ2,11ρ2,3ρ2,4ρ2,Mρ3,1ρ3,21ρ3,4ρ3,Mρ4,1ρ4,2ρ4,31ρ4,MρM,1ρM,2ρM,3ρM,41]

where all elements in the diagonal and lower triangle (Cijij) 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.

Code availability

A 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

References

  1. 1
  2. 2
  3. 3
  4. 4
    CircStat: a Matlab toolbox for circular statistics
    1. P Berens
    (2009)
    Journal of Statistical Software, 31, 10.18637/jss.v031.i10.
  5. 5
    Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type
    1. G Bi
    2. M Poo
    (1998)
    Journal of Neuroscience 77:551–555.
  6. 6
    Cellular basis of EEG slow rhythms: a study of dynamic corticothalamic relationships
    1. D Contreras
    2. M Steriade
    (1995)
    Journal of Neuroscience 15:604–622.
  7. 7
  8. 8
  9. 9
    Spatiotemporal patterns of spindle oscillations in cortex and thalamus
    1. D Contreras
    2. A Destexhe
    3. TJ Sejnowski
    4. M Steriade
    (1997)
    The Journal of Neuroscience 17:1179–1196.
  10. 10
  11. 11
    Spatiotemporal analysis of local field potentials and unit discharges in cat cerebral cortex during natural wake and sleep states
    1. A Destexhe
    2. D Contreras
    3. M Steriade
    (1999)
    The Journal of Neuroscience 19:4595–4608.
  12. 12
  13. 13
    Hilbert Transform Applications in Mechanical Vibration (1st ed)
    1. M Feldman
    (2011)
    West Sussex: Wiley.
  14. 14
  15. 15
    Statistical Analysis of Circular Data (1st ed)
    1. N Fisher
    (1993)
    Cambridge: Cambridge University Press.
  16. 16
  17. 17
    Learning-dependent increases in sleep spindle density
    1. S Gais
    2. M Mölle
    3. K Helms
    4. J Born
    (2002)
    Journal of Neuroscience 22:6830–6834.
  18. 18
  19. 19
    Sensory modality specificity of neural activity related to memory in visual cortex
    1. JR Gibson
    2. JH Maunsell
    (1997)
    Journal of Neurophysiology 78:1263–1275.
  20. 20
    Feedforward and feedback connections between areas V1 and V2 of the monkey have similar rapid conduction velocities
    1. P Girard
    2. JM Hupé
    3. J Bullier
    (2001)
    Journal of Neurophysiology 85:1328–1331.
  21. 21
    Heterogeneous Origins of Human Sleep Spindles in Different Cortical Layers
    1. D Hagler
    2. S Cash
    3. E Halgren
    (2016)
    Heterogeneous Origins of Human Sleep Spindles in Different Cortical Layers.
  22. 22
  23. 23
  24. 24
  25. 25
    Topics in Circular Statistics (1st ed)
    1. S Jammalakadaka
    2. A Sengupta
    (2001)
    Singapore: World Scientific.
  26. 26
    Cellular-synaptic generation of sleep spindles, spike-and-wave discharges, and evoked thalamocortical responses in the neocortex of the rat
    1. A Kandel
    2. G Buzsáki
    (1997)
    Journal of Neuroscience 17:6783–6797.
  27. 27
  28. 28
    Spindle waves are propagating synchronized oscillations in the ferret LGNd in vitro
    1. U Kim
    2. T Bal
    3. DA McCormick
    (1995)
    Journal of Neurophysiology 74:1301–1323.
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
    Self-organizing neural systems based on predictive learning
    1. RPN Rao
    2. TJ Sejnowski
    (2003)
    Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 361:1149–1175.
    https://doi.org/10.1098/rsta.2003.1190
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
    Cortical Areas: Unity and Diversity
    1. A Schüz
    2. V Braitenberg
    3. R Miller
    (2002)
    377–385, The human cortical white matter: quantitative aspects of cortico-cortical long-range connectivity, Cortical Areas: Unity and Diversity, London, UK, Taylor Francis.
  48. 48
    The Organization of the Cerebral Cortex
    1. DA Sholl
    (1956)
    John Wiley.
  49. 49
    The visual scoring of sleep in adults
    1. MH Silber
    2. S Ancoli-Israel
    3. MH Bonnet
    4. S Chokroverty
    5. MM Grigg-Damberger
    6. M Hirshkowitz
    7. S Kapen
    8. SA Keenan
    9. MH Kryger
    10. T Penzel
    11. MR Pressman
    12. C Iber
    (2007)
    Journal of Clinical Sleep Medicine 3:121–131.
  50. 50
  51. 51
  52. 52
  53. 53

Decision letter

  1. Frances K Skinner
    Reviewing Editor; University Health Network, Canada

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Repeating circular waves enable strengthening of large-scale neural assemblies during sleep spindles in human cortex" for consideration by eLife. Your article has been favorably evaluated by Timothy Behrens (Senior Editor) and three reviewers, one of whom is a member of our Board of Reviewing Editors.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

This work is an analysis of spindle oscillations recorded from human subdural grids placed for the purpose of localizing epileptic seizures. The authors examine the spatio-temporal structures of spindles during sleep in human, and their possible role in memory consolidation. The techniques used are developed from their earlier work where waves were found in monkey visual cortex. Used here in considering both amplitude and phase information, the authors show that spindles exhibit a sophisticated spatio-temporal structure, with the majority of spindles traveling from temporal to parietal to frontal (TPF) cortex. The authors suggest that such traveling waves could serve to facilitate spike-timing dependent plasticity between anatomically remote areas, which would otherwise be difficult due to long axonal delays. The analysis is rigorous and the application of methods to detect propagating waves to spindles is novel. Spindles are of high interest for memory consolidation, and this dataset is highly valuable.

The revision requirements are three-fold.

1) An expansion and clarification of methodological aspects is needed in various places. Specific points are as follows.

a) Rotating spindle waves: The presented data and analyses on spindle oscillations are clear. They rely on a large dataset of spindle oscillations and the provided movies give a good sense of the observations made. However, clarifications are needed. First, the authors should incorporate a definition of what they consider as rotating and expanding waves, as well as the specific criterion used to distinguish these. Definitions should appear in the main text, and the supplementary material should include the specific criteria used for the classification of spatio-temporal patterns of activity. In the absence of this, it remains quite difficult to validate the method.

The methods present clearly the identification of the center point of the oscillation, important to perform a polar parameterization, thus essential to extract rotating or expanding waves. This step is well described and the method relies on the computation of gradients and curl. It would be important to know:

i) If the authors pre-processed the data and if smoothing occurred to compute these differential elements. In particular, filtering is indicated in Figure 1—figure supplement 8. Was a similar filtering used in the identification of waves? If this was not the case, how did the authors correct high frequency signals that may be amplified in the computation of differential quantities?

ii) Once the polar coordinate system has been derived and the phase map computed, how did the classification take place and what were the error rates accepted for this classification? The statement about the proportion of spindle oscillations classified as rotating waves of ~50% remains a little bit hard to appreciate since there are only 2 classes (plus "complex" waves that are those not satisfying their criteria of rotating or expanding wave). Moreover, the analysis of Figure 1—figure supplement 7 provides a more contrasted view of this result: rotational waves are not necessarily the leading form of oscillation, except for subject 1 (for other subjects, spindles were generally classified as complex).

iii) The authors indicate that among the rotating waves, they identified a clear bias towards TPF direction of rotation. Was this observation supported by a statistical test?

iv) The waves observed indeed reflect a peak between 3-5 m/s, but also a very high variance. From the individual curves, it seems that we see the superposition of two speed distributions centered at distinct values. Could the authors comment on the variations of the speed value provided?

b) The analysis of intra-class similarity is nicely developed. However, for the level of synchronization, it was unclear how the millisecond precision was computed. Particularly, it was unclear what was meant by "in two cycles with the highest similarity". Could the authors elaborate in the Methods or in the text on this choice?

c) A rigorous demonstration that spindles are waves would be a highly novel finding. However, given the limited size of the dataset, and the relationship of spindles to the underlying anatomy being unclear, additional analyses and explanations are needed.

i) Discrepancy with "local spindles", different frequencies of controparietal/frontal spindles finding, and relationship to underlying anatomy. Much previous work concludes that human spindles are "local" and that spindle frequencies vary between areas (i.e. Andrillon et al. 2011, Peter-Derex et al. 2012, Nir et al. 2014). However, here the statement is made that this is in fact not so and the explanation is provided that this might be due to using a threshold for detection. It is unclear how this analysis supports this statement. What would the wave-detection approach show if these spindles would in fact be local? All that is done here is to fit a field of vectors in 2D space and then to see if a curl can be fit. No goodness of fit is assessed (only that it is significant). But presumably this would also work for fields where only subsets of the vectors are fit, with the others random. Also, spindles have different frequencies depending on anatomical origin as well as in time relative to spindle onset. Could these frequency changes "generate" the appearance of waves (or vice-versa)?

ii) Structure of waves. It was unclear whether the origin (center) of waves remained constant, changed wave-by-wave or changed slowly. Figures are shown "re-centered", so it isn't apparent where the actual center was. Also, how sensitive is this method to edge effects, i.e., is it equally sensitive to detect centers everywhere on the array or only in the center? Overall, the location of the centers and the underlying anatomy should be clarified. The authors suggest a TPF direction of spread, so does this indicate center preferentially in temporal lobe?

2) A toning down of interpretation and claims and situating the work more appropriately. That is, STDP and memory consolidation aspects aren't actually measured. Specifically:

a) The authors should consider adjusting their title somewhat in toning down the claims (i.e., "… enable strengthening of large-scale neural assemblies…") to more appropriately reflect the main aspects of the paper.

b) Mechanistic interpretation is not supported. The Abstract talks extensively about synapses and STDP, neither of which is measured here. This should be clarified (e.g., interesting suggestions, but cannot be assessed using this data). Same for "spiking activity" – all that is measured is γ band power. Statements such as "groups of spikes" are not supported by this data. Others have measured spindles and spikes simultaneously in humans (i.e. Andrillon et al., as cited), showing a complex relationship.

c) Another suggestion is that the authors expand on their Figure 2 description/illustration. For example, could the authors say something about spindle mechanisms (RT cells and mutual inhibition etc.), and how such spikings illustrated in Figure 2 translate to EPSPs (IPSPs?). Haider et al. is cited to justify LFP and synaptic current tight couplings, but this includes IPSPs. Related to this is the consideration of different STDP rules for different cell types (e.g., see Abbott and Nelson, Nature Neuroscience 2000).

In essence, since much has been published regarding STDP, spindling mechanisms, memory consolidation etc., it would be good for the authors to put these various pieces together in their suggestions for the general reader to appreciate more fully, and so that interpretations/limitations/assumptions are clear in what is being said. As it stands, it comes across as overly simplistic regarding the timing relationships and LTP/LTD (Figure 2) and traveling waves with STDP aspects and so on.

3) Clearly describing/demonstrating that a differentiation between the two types of 'waves' is possible along with what is the current consensus in the literature – 'local spindles' and how it was ruled out with their analyses. This is unclear at present. Perhaps some of the unclassified waves are of that nature? That is statistics and classification are not always clear [see point (1) above].

Additional points to address:

i) In general, please provide more statistical support and more information on the statistical tests and the data analysis procedures used. It is stated that the code is available on bitbucket, but when looked for by one of the reviewers, access was not possible, thus preventing a better understanding of the types of tests used.

ii) How were the recordings referenced? Phase analysis is highly sensitive to this.

iii) What Figure 4C quantifies is unclear – please clarify how this is calculated. Can one show the same data for expanding waves and show this critical analysis for multiple patients?

https://doi.org/10.7554/eLife.17267.025

Author response

[…] The revision requirements are three-fold.

1) An expansion and clarification of methodological aspects is needed in various places. Specific points are as follows.

a) Rotating spindle waves: The presented data and analyses on spindle oscillations are clear. They rely on a large dataset of spindle oscillations and the provided movies give a good sense of the observations made. However, clarifications are needed. First, the authors should incorporate a definition of what they consider as rotating and expanding waves, as well as the specific criterion used to distinguish these. Definitions should appear in the main text, and the supplementary material should include the specific criteria used for the classification of spatio-temporal patterns of activity. In the absence of this, it remains quite difficult to validate the method.

This is an excellent suggestion; we have now inserted explicit definitions for rotating and expanding waves into the main text (second paragraph). We have further collected the specific classification criteria into two discrete Methods subsections (“Detection of expanding waves” and “Detection of rotating waves”), in addition to adding another figure supplement detailing the process for detecting expanding waves (Figure 1—figure supplement 2). We agree that these additions greatly improve the presentation of analysis techniques.

The methods present clearly the identification of the center point of the oscillation, important to perform a polar parameterization, thus essential to extract rotating or expanding waves. This step is well described and the method relies on the computation of gradients and curl. It would be important to know:

i) If the authors pre-processed the data and if smoothing occurred to compute these differential elements. In particular, filtering is indicated in Figure 1—figure supplement 8. Was a similar filtering used in the identification of waves? If this was not the case, how did the authors correct high frequency signals that may be amplified in the computation of differential quantities?

As in our previous work (Muller et al., 2014), the robust nonparametric smoothing technique of Garcia (2010) is employed before center localization (and, here, computation of differential quantities). It is important to note that after center localization, however, the phase correlation quantities are calculated on the unsmoothed phase fields, in order to prevent any smoothing effects from contaminating wave detection results. This was stated in the original manuscript, and in the revised submission we have striven further to clarify this important point (subsection “Spatiotemporal dynamics”, second paragraph and subsection “Detection of expanding waves”, first paragraph).

ii) Once the polar coordinate system has been derived and the phase map computed, how did the classification take place and what were the error rates accepted for this classification?

In order to validate our computational approach for the 48-96 channel ECoG arrays, we analyzed the performance of our algorithm in simulated data under varying degrees of noise. Because our method takes a statistical approach to wave detection, the algorithm detects waves with fidelity up to high noise levels, as opposed to algebraic methods that tend to break down in the presence of noise (Wong and Yip, Pattern Recognition42, 2009). While we did not include these verification tests in the original manuscript, this comment has emphasized to us the importance of doing so, and we now include these tests in the revised submission (see Figure 1—figure supplements 11 and 12). As demonstrated in Figure 1—figure supplement 11C, our approach robustly identifies simulated expanding and rotational waves while reporting no detections with only noise. We thank the Reviewers for this helpful comment.

The statement about the proportion of spindle oscillations classified as rotating waves of ~50% remains a little bit hard to appreciate since there are only 2 classes (plus "complex" waves that are those not satisfying their criteria of rotating or expanding wave). Moreover, the analysis of Figure 1—figure supplement 7 provides a more contrasted view of this result: rotational waves are not necessarily the leading form of oscillation, except for subject 1 (for other subjects, spindles were generally classified as complex).

This was unclear in the original submission, and we thank the Reviewers for making this important point. The critical ambiguity here was to designate cycles that are not classified as rotational or expanding waves as complex, rather than simply unclassified. Our initial logic was that, because we observe very few spindle oscillation cycles that are fully synchronized across the array (i.e. temporal extent less than 10 ms), we termed the unclassified oscillation cycles “complex”; however, this does not allow for the fact that some spindle cycles may in fact be “local”. We have thus updated the term “complex” in the manuscript simply to “unclassified”, and explained this more clearly in the revised text (main text, seventh paragraph and Figure 1—figure supplement 9). We thank the Reviewers for this helpful suggestion.

iii) The authors indicate that among the rotating waves, they identified a clear bias towards TPF direction of rotation. Was this observation supported by a statistical test?

We have taken this comment into account and now provide the results of an exact binomial test on the bias towards TPF direction (main text, second paragraph). The bias is significant for each individual subject (p < 10-3 in each case), as well as for the full population (p < 10-10, 14,796 TPF cycles, 21,272 total). We thank the Reviewers for this helpful suggestion.

iv) The waves observed indeed reflect a peak between 3-5 m/s, but also a very high variance. From the individual curves, it seems that we see the superposition of two speed distributions centered at distinct values. Could the authors comment on the variations of the speed value provided?

These speed distributions are similar in extent to those observed in previous literature (Rubino et al., Nature Neuroscience9, 2006, cf. their Supplementary Figure 2D), and are strongly peaked near values corresponding to physiological measurement. Looking further into the distributions, high variance is not apparent: from the 20th to the 80th percentiles, speeds in the full distribution vary over a tight range (3-9 m/s), corresponding approximately to 20-50 ms of travel time across the array. Speed distributions for individual subjects are similar (from 20-80th percentile: 4-10, 3-8, 3-10, 2-3, and 4-13 m/s, respectively). The variation in these distributions is due in part to measurement noise, with noisy but statistically significant events showing up as high-speed waves on the array. We thank the reviewer for this helpful comment, which we have now pointed out in the main text (second paragraph).

Further, true variations in wave speed can certainly result in better or worse alignment of spikes across distant cortical regions, as discussed in our Figure 2. As mentioned in our response to point 2c, we obtained in previous work a mathematical expression for the expected weight change by STDP in this case, which allows us to understand the effect of speed variability on the plasticity mechanisms discussed here. We fully agree to include discussion of this important point in the main text (last paragraph), and we thank the reviewers for this helpful comment.

b) The analysis of intra-class similarity is nicely developed. However, for the level of synchronization, it was unclear how the millisecond precision was computed. Particularly, it was unclear what was meant by "in two cycles with the highest similarity". Could the authors elaborate in the Methods or in the text on this choice?

We agree that it is useful for us to clarify this important point. In the revised submission, we have updated the main text to describe this calculation more clearly (fifth paragraph).

c) A rigorous demonstration that spindles are waves would be a highly novel finding. However, given the limited size of the dataset, and the relationship of spindles to the underlying anatomy being unclear, additional analyses and explanations are needed.

We appreciate this positive comment on our work, and we agree that reporting spindles are spatiotemporally organized into global traveling waves is an important result with strong implications for spindles’ established function in sleep-dependent consolidation of long-term memory. In the revised submission, we have striven to clarify analysis procedures and to provide additional evidence to further support these results (detailed point-by-point in the responses below).

i) Discrepancy with "local spindles", different frequencies of controparietal/frontal spindles finding, and relationship to underlying anatomy. Much previous work concludes that human spindles are "local" and that spindle frequencies vary between areas (i.e. Andrillon et al. 2011, Peter-Derex et al. 2012, Nir et al. 2014). However, here the statement is made that this is in fact not so and the explanation is provided that this might be due to using a threshold for detection.

In accordance with point 1.ii above, we recognize the importance of clarifying the relationship of the observed waves to previous observations of “local spindles”, and we thank the Reviewers for pointing this out. In the revised submission, we now denote oscillation cycles unclassified by the algorithm simply as “unclassified”, rather than “complex” (Figure 1—figure supplement 9). This shift in terminology now allows for the possibility of local spindles, which constitute approximately half of the observed events (< 40% channel involvement) in previous work (Andrillon et al., 2011). We acknowledge the possibility that it is in fact the set of “global” spindles which are robustly (and almost completely) organized as rotational and expanding traveling waves across the cortical surface (main text, seventh paragraph).

In fact, this variation in spatial extent suggests an interesting functional role: while in some cases, memory content will be specialized and thus confined to a single cortical region (Khazipov et al., Nature432, 2004), in other cases memories will contain content across modalities (visual, auditory, or cognitive) that must be distributed across regions, an organization which requires the spatiotemporal patterns observed here (as illustrated in Figure 2). We have now updated our submission to include this point (main text, last paragraph), and we thank the Reviewers for their helpful input.

Finally, it is important to emphasize, as in our original submission, that our phase-based method represents a well-controlled, robust approach to detecting waves in noisy multichannel data. As such, it is not surprising that we are able to detect broad spatiotemporal patterns of activity in cortex where past studies have detected only noisy patches of spindle activity in a handful of electrodes with amplitude-duration detection (Mölle et al., J Neurosci22, 2002). Our approach is robust to false wave detections and sensitive to coherent activity patterns embedded in noise. In order to clarify these points, we have included several control analyses demonstrating the validity of our approach. We detail these specifically in the response to the next few points, and we have updated the manuscript to clarify these aspects of our method.

It is unclear how this analysis supports this statement. What would the wave-detection approach show if these spindles would in fact be local? All that is done here is to fit a field of vectors in 2D space and then to see if a curl can be fit. No goodness of fit is assessed (only that it is significant).

We understand and appreciate the reviewers’ concern on this point. We must emphasize, however, that the above characterization of our method is inaccurate. As stated in the Methods, we only use the curl operation to estimate the putative center of rotation (step 2, Figure 1—figure supplement 3), after which we calculate the circular-circular correlation between signal phase and rotation about this point (step 3, Figure 1—figure supplement 3), which “quantifies the strength of the spatiotemporal pattern of activity … in a single number” (subsection “Detection of expanding waves”). This number, ϱφ,θ, represents in itself a “goodness of fit”, whose absolute value ranges from 0 (no pattern) to 1 (fully organized rotational pattern).

In addition to this, we employed the average vector field quantification as an independent control (Figure 1—figure supplement 5). In both the expanding and rotating wave cases, this control demonstrates that we correctly detect and separate the two classes of waves considered, providing an additional demonstration of “goodness of fit”, calculated in a different but complementary manner to our primary wave detection approach.

Further, it is important to note that the expanding and rotating patterns extend out to the edges of the array in Figure 1—figure supplement 5, indicating in accordance with the phase correlation measure that these waves are truly global coherent activity patterns. In the revised manuscript, we have added a figure supplement illustrating the full distribution of propagation directions at each point in the average vector field for rotational TPF waves in all subjects (gray lines, Figure 1—figure supplement 6), in addition to the circular mean (black arrows, as in other figures), in order to clearly demonstrate that the observed patterns are indeed widespread and global.

We hope that these additional explanations fully address the reviewers’ concern on this point. In accordance with this, we have striven in our revision to further clarify our computational algorithm for wave detection and the controls validating this approach.

But presumably this would also work for fields where only subsets of the vectors are fit, with the others random.

In order to address this concern directly, we have included an additional control comparing detection of a global rotating wave versus a very local one (6 electrodes, or 9.4% of the 64 channel electrode array; Figure 1—figure supplement 12). This control demonstrates that if few electrodes participate in the phase pattern, as would be the case with a fully “local” spindle, this will not be detected. This control illustrates that finding high ϱφ,θ is strongly indicative of a distributed and global phase pattern. This result is further supported by the average vector field control calculations (Figure 1D and Figure 1—figure supplements 57, as detailed in the point above).

Further, it is important to note that this approach, together with the sampling density available in this study, emphasizes global phase patterns and will be less sensitive to smaller, intermediate patterns. This is perhaps why we see a robust rotating TPF wave organization in the vector field average over all cycles in each individual subject (Figure 1—figure supplement 7). Thus, while we can be confident that the rotating wave organization reported appears with at least the occurrence rate reported here, the true incidence rate is most likely higher (including any intermediate waves). This point will be pursued in future work, utilizing higher density electrode arrays.

We thank the reviewers for this interesting and helpful comment, and we believe that it has greatly improved our work.

Also, spindles have different frequencies depending on anatomical origin as well as in time relative to spindle onset. Could these frequency changes "generate" the appearance of waves (or vice-versa)?

This control was one of the first concerns raised in validating our results. In order to test this, we re-ran our analysis on one stage 2 sleep session from Subject 1 with surrogate data, where each electrode evolved in time with its mean instantaneous frequency during the spindle, with a randomized initial phase angle. In this control, both rotating and expanding wave patterns were highly decreased (3.5% and 0.8% of cycles with rotating and expanding waves, respectively, in the control simulation versus 64% and 14% in the original). This control demonstrates that frequency shifts cannot explain our results, indicating that the observed waves are due to a true phase shift across channels. To create these phase shifts, active corticocortical, thalamocortical, and corticothalamic connections are necessary. We plan to address detailed mechanisms for such network-level generation of the spatiotemporal activity patterns in future modeling work. We thank the reviewer for addressing this important point, the control for which we now discuss in the Methods (subsection “Simulated data controls”, last paragraph).

ii) Structure of waves. It was unclear whether the origin (center) of waves remained constant, changed wave-by-wave or changed slowly. Figures are shown "re-centered", so it isn't apparent where the actual center was. Overall, the location of the centers and the underlying anatomy should be clarified. The authors suggest a TPF direction of spread, so does this indicate center preferentially in temporal lobe?

We agree that this is an important point, and welcome the opportunity to expand upon it. In the revised submission, we have included an additional figure (Figure 1—figure supplement 8) to illustrate the distribution of rotation center (across all 10,944 rotating TPF waves in Subject 1) and movement between cycles (aggregated across TPF waves in all subjects). These additional panels demonstrate that the center of rotation distribution appears to be concentrated just dorsal to the temporal lobe, aligned on the Sylvian fissure (Figure 1—figure supplement 8A). The rotation center further appears to move slowly across the surface of the cortex, if at all (Figure 1—figure supplement 8B), in agreement with the impression from Video 1. We believe that this analysis, along with our response to the following query, fully addresses the reviewers’ concerns on this point.

Also, how sensitive is this method to edge effects, i.e., is it equally sensitive to detect centers everywhere on the array or only in the center?

We have tested this point directly using simulated rotating waves in surrogate data, and in the revised submission, we have included this analysis in an additional figure (Figure 1—figure supplement 11D). In short, we are able to detect noisy waves embedded across the array. Though sensitivity is decreased at the edges, we are confident that the distribution of detected rotation centers is unaffected by issues of sensitivity.

2) A toning down of interpretation and claims and situating the work more appropriately. That is, STDP and memory consolidation aspects aren't actually measured. Specifically:

a) The authors should consider adjusting their title somewhat in toning down the claims (i.e., "… enable strengthening of large-scale neural assemblies…") to more appropriately reflect the main aspects of the paper.

While we do believe our manuscript contains intriguing preliminary evidence towards strengthening of repeating patterns, and similar evidence in visual cortex of the mouse has previously been interpreted as involving plasticity-dependent mechanisms (Han et al., Neuron60, 2008), we fully understand the reviewers’ concern on this point and agree that adjusting the title can more accurately reflect the analyses presented in the paper. In line with this suggestion, we have updated the title in our revised submission to:

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

We thank the reviewers for this helpful comment.

b) Mechanistic interpretation is not supported. The Abstract talks extensively about synapses and STDP, neither of which is measured here. This should be clarified (e.g., interesting suggestions, but cannot be assessed using this data). Same for "spiking activity" – all that is measured is γ band power. Statements such as "groups of spikes" are not supported by this data. Others have measured spindles and spikes simultaneously in humans (i.e. Andrillon et al., as cited), showing a complex relationship.

We understand this point, and we regret that our initial submission did not clearly review the background knowledge from animal electrophysiology on this point. A wealth of evidence from intracellular and single-unit recordings presented over two decades shows that during spindles, cells in neocortex receive strong excitatory input from thalamus and concentrate spikes at the peak of depth negativity (Contreras and Steriade, J Neurosci15, 1995; Contreras and Steriade, J Physiol 490, 1996; Kandel and Buzsaki, J Neurosci 17, 1997; Peyrache et al., PNAS 108, 2011). In this work, we analyze high-γ power (HGP), which is a well-accepted proxy for spiking activity in the field (Ray et al., J Neurosci 28, 2008; Manning et al., J Neurosci 29, 2009; Ray and Maunsell, PLoS Biology 9, 2011; Ray, Curr Opin Neurobio 31, 2015). Our results are fully in line both with those obtained in animal recordings and previous work in human – specifically Andrillon et al. (2011), whose Figure 4D matches very closely with Figure 1E in Peyrache et al. (2011).

Note that phase 0 in Andrillon et al. (2011) corresponds to π in Peyrache et al. (2011), since Peyrache et al. took deep layer LFP as reference. We believe that the preponderance of evidence, both from previous animal electrophysiology and from human recordings, is so strong on this point that our demonstration of preferred phase for HGP activity (Figure 4B) should be viewed as mere confirmation of previous work.

At the same time, we do fully understand that neither spikes nor STDP were directly measured in this study, only proxies for each; thus, while we took care not to misrepresent the results we provide as definitive evidence for this effect, we have revised the manuscript accordingly, ensuring that all functional statements are carefully supported by the data analyses presented here and knowledge from previous studies in animal electrophysiology. Finally, we have updated the main text to further clarify our use of high-γ power, including a clarification of terminology (“HGP” versus “broadband power shift”).

c) Another suggestion is that the authors expand on their Figure 2 description/illustration. For example, could the authors say something about spindle mechanisms (RT cells and mutual inhibition etc.), and how such spikings illustrated in Figure 2 translate to EPSPs (IPSPs?). Haider et al. is cited to justify LFP and synaptic current tight couplings, but this includes IPSPs.

Related to the previous point, we have striven in the revised manuscript to relate the known mechanisms of cortical spindling activity to our results. Specifically, it is known that neocortical networks receive strong rhythmic excitatory volleys from thalamocortical cells, driving the spindle rhythm (see citations above). These strong excitatory inputs modulate firing of cortical pyramidal cells and interneurons (Peyrache et al., 2011), whose EPSPs and IPSPs, respectively, will contribute to the cortical LFP during the spindles measured here.

In order to elaborate on the spiking activity illustrated in Figure 2, we have more clearly emphasized in the revised text that the spikes considered originate from pyramidal cells, the dominant contributor to the white matter tracts (Sholl, The Organization of the Cerebral Cortex, 1956; Schüz and Braitenberg, in Cortical Areas: Unity and Diversity, 2003). It is these excitatory long-range axons that serve to link distant cell groups across the cortex, forming the basis for distributed long-term memories. It is thus their spiking activity – and resulting EPSPs – on which we focus. In accordance with this point, we have clarified the discussion of Figure 2 in the main text (first, fourth and fifth paragraphs).

Related to this is the consideration of different STDP rules for different cell types (e.g., see Abbott and Nelson, Nature Neuroscience 2000).

In essence, since much has been published regarding STDP, spindling mechanisms, memory consolidation etc., it would be good for the authors to put these various pieces together in their suggestions for the general reader to appreciate more fully, and so that interpretations/limitations/assumptions are clear in what is being said. As it stands, it comes across as overly simplistic regarding the timing relationships and LTP/LTD (Figure 2) and traveling waves with STDP aspects and so on.

In context of the previous point, we are most interested in the spiking (and plasticity) behavior of neocortical pyramidal cells, making long-range projections to pyramidal cells and interneurons in cortex. In previous work, we have studied this problem theoretically and provided an exact analytical result (Muller et al., Frontiers in Computational Neuroscience5, 2011). Further work has generalized this to different plasticity rules (Luz and Shamir, PLoS Computational Biology12, 2016), such as that observed in excitatory-to-inhibitory connections (cf. their Equation 4). Thus, while we present here a simplified schema to illustrate the functional consequences of the observed spatiotemporal activity patterns, we do also understand the system from a theoretical and mathematical perspective. We thank the reviewer for pointing out this initial oversight, and we have updated the main text to mention previous theoretical work and its connection to our submission (last paragraph). Finally, we note that this theoretical understanding also allows us to understand the functional implications for the variance in distribution of wave speed (cf. point 1.iv).

3) Clearly describing/demonstrating that a differentiation between the two types of 'waves' is possible along with what is the current consensus in the literature – 'local spindles' and how it was ruled out with their analyses. This is unclear at present. Perhaps some of the unclassified waves are of that nature?

In accordance with point 1.ii above, we appreciate this point and have updated the revised manuscript to include precisely the possibility that some of the unclassified waves are local in nature (main text, seventh paragraph), allowing us to more carefully contextualize our work in relation to previous results. We have also clarified our explanation of the definitions and method for distinguishing rotating from expanding waves, in addition to emphasizing the validation of this distinction in Figure 1—figure supplement 5.

That is statistics and classification are not always clear [see point (1) above].

In accordance with point (1) above, we have carefully considered the presentation of methods and statistical tests, and in our revised submission, we have expanded our explanations and the details provided for all calculations and statistical tests.

Additional points to address:

i) In general, please provide more statistical support and more information on the statistical tests and the data analysis procedures used. It is stated that the code is available on bitbucket, but when looked for by one of the reviewers, access was not possible, thus preventing a better understanding of the types of tests used.

We welcome the opportunity to expand on the Methods employed in this work. In line with this and other reviewers’ comments, we have moved definitions for analysis terms to the main text, included more information on statistical tests, and expanded the description of our algorithmic approach.

ii) How were the recordings referenced? Phase analysis is highly sensitive to this.

In the Methods section of the original manuscript, we stated that “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”. In the revised manuscript, we have clarified the reference for the recordings and that we replicated our results with an average reference, in order to rule out any sensitivity in our analysis to reference effects. At the same time, however, we note that while our results are unaffected by such re-referencing techniques, this is not the preferred approach when analyzing phase in electrode recordings (Shirhatti et al., Neural Computation 28, 2016).

iii) What Figure 4C quantifies is unclear – please clarify how this is calculated. Can one show the same data for expanding waves and show this critical analysis for multiple patients?

We are certainly interested in this analysis and welcome the opportunity to expand on this point. In the revised manuscript, we have clarified the presentation of Figure 4C and included a figure supplement illustrating the results for expanding waves (Figure 4—figure supplement 1). In addition, we have now included an additional analysis showing this result for multiple subjects (Figure 4—figure supplement 2). We believe these additional analyses will address the point raised by the Reviewers, and thank them for raising this important point.

https://doi.org/10.7554/eLife.17267.026

Article and author information

Author details

  1. Lyle Muller

    1. Computational Neurobiology Laboratory, Salk Institute for Biological Studies, La Jolla, United States
    Contribution
    LM, Conception and design, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon 0000-0001-5165-9890
  2. Giovanni Piantoni

    1. Department of Neurology, Massachusetts General Hospital, Harvard Medical School, Boston, United States
    Contribution
    GP, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon 0000-0002-5308-926X
  3. Dominik Koller

    1. Computational Neurobiology Laboratory, Salk Institute for Biological Studies, La Jolla, United States
    Contribution
    DK, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  4. Sydney S Cash

    1. Department of Neurology, Massachusetts General Hospital, Harvard Medical School, Boston, United States
    Contribution
    SSC, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  5. Eric Halgren

    1. Department of Radiology, University of California, San Diego, San Diego, United States
    2. Department of Neurosciences, University of California, San Diego, San Diego, United States
    Contribution
    EH, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  6. Terrence J Sejnowski

    1. Computational Neurobiology Laboratory, Salk Institute for Biological Studies, La Jolla, United States
    Contribution
    TJS, Conception and design, Analysis and interpretation of data, Drafting or revising the article
    For correspondence
    1. terry@salk.edu
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon 0000-0002-0622-7391

Funding

Swartz Foundation

  • Terrence J Sejnowski

Howard Hughes Medical Institute

  • Terrence J Sejnowski

Office of Naval Research (MURI N000141310672)

  • Sydney S Cash
  • Eric Halgren
  • Terrence J Sejnowski

National Institutes of Health (R01EB009282)

  • Terrence J Sejnowski

National Institutes of Health (5T32EY20503-5)

  • Lyle Muller

Bial Foundation (BIAL 220/12)

  • Giovanni Piantoni

Office of Naval Research (N000141210299)

  • Terrence J Sejnowski

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

The authors would like to thank the clinical patients for their participation in the research, in addition to T Bartol, F Chavane, A Destexhe, and CF Stevens for helpful discussions. The authors are grateful to the Swartz Foundation, the Howard Hughes Medical Institute, the Office of Naval Research, and NIH for support.

Ethics

Human subjects: Patients 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 (see Materials and Methods - Subjects).

Reviewing Editor

  1. Frances K Skinner, Reviewing Editor, University Health Network, Canada

Publication history

  1. Received: April 26, 2016
  2. Accepted: October 19, 2016
  3. Version of Record published: November 15, 2016 (version 1)

Copyright

© 2016, Muller et al

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,138
    Page views
  • 485
    Downloads
  • 1
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Comments

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Genomics and Evolutionary Biology
    2. Microbiology and Infectious Disease
    Kevin S Bonham et al.
    Research Article