Synchronization and Propagation of Global Sleep Spindles

Sleep spindles occur thousands of times during normal sleep and can be easily detected by visual inspection of EEG signals. These characteristics make spindles one of the most studied EEG structures in mammalian sleep. In this work we considered global spindles, which are spindles that are observed simultaneously in all EEG channels. We propose a methodology that investigates both the signal envelope and phase/frequency of each global spindle. By analysing the global spindle phase we showed that 90% of spindles synchronize with an average latency time of 0.1 s. We also measured the frequency modulation (chirp) of global spindles and found that global spindle chirp and synchronization are not correlated. By investigating the signal envelopes and implementing a homogeneous and isotropic propagation model, we could estimate both the signal origin and velocity in global spindles. Our results indicate that this simple and non-invasive approach could determine with reasonable precision the spindle origin, and allowed us to estimate a signal speed of 0.12 m/s. Finally, we consider whether synchronization might be useful as a non-invasive diagnostic tool.


Introduction
Synchronization is a robust and widespread phenomenon in complex systems that usually occurs in oscillatory interacting units. There are many examples of synchronization in various cells and organisms, including cardiac pacemaker cells, rhythmically flashing fireflies, chorusing crickets and neurons [1]. Since Crick and Koch proposed that conscience phenomena could be related to synchronization of different cerebral regions that impose a temporary global unity to the brain [2], several other studies have investigated the role that synchronization plays in different aspects of neuronal physiology [3][4][5][6][7][8][9].
Electroencephalography (EEG) measures electrical oscillations of the brain through electrodes placed on the scalp [10]. This technique can be seen as a non-invasive window into the internal rhythms of the brain. EEG signals are rich in detail, and by nature are non-stationary [11][12][13]. Information gained from EEG during sleep may be even richer than measurements taken during wakefulness, and EEG represents an important method to assess normal and altered physiology because sleep EEG transients may be affected by some diseases. Of all EEG elements, sleep spindles (SS) are certainly the most studied short transient [14,15]. SS were first documented in the 1930s [16] and are now known to be associated with memory processing and learning [17][18][19].
SS are defined as wave packets between 11 Hz and 16 Hz that have a duration of less than 2 s and primarily occur during non-REM (NREM) sleep stage II, although they may also occur during other sleep stages [20,21]. SS can be represented as an oscillatory function (wavefunction) with an average frequency modulated by an envelope [14]. Although this envelope is not exactly Gaussian for true SS, this model can nevertheless be used to describe the time limits for their occurrence [14,22]. By determining the ratio between this wave and its envelope, we can perform a normalization to retain only the oscillatory component of the signal. As such, the signal envelope can be used for wave propagation analysis (spread across electrodes) and the oscillatory component can be applied to synchrony analysis. SS can be easily detected and characterized, which allows an independent investigation of envelopes and oscillatory behaviors. Thus, SS are an excellent model for evaluating brain synchronization. The widespread availability and non-invasiveness of EEG could allow it to be used to explore the relevance of synchronization issues to different mental diseases, provided that powerful and automated accompanying diagnostic tools can be developed.
During the course of one night's sleep, spindles appear thousands of times and they likely have a thalamic origin. Moreover, scalp-based measurement of SS can reflect the activity of the reticulo-thalamic-cortex communication network [23][24][25]. SS have rich behaviors with amplitude and frequency modulation (short time chirp) that shows high non-stationarity [22,26,27]. Indeed, there are two types of SS, slow and fast, which appear to have different physiological characteristics [14]. In a study using Ca and Na channel-blocking drugs, different pharmacological responses for fast and slow spindles have suggested that different neuronal populations may be implicated in the production of those two spindle types [28]. It is worthwhile to note that a systematic increase in negative SS chirp factor can be associated with a shift towards lower frequencies if the FFT power spectrum is analyzed. As with any event detected on the surface of the scalp, SS represent the oscillations of a large set of neurons, and are prone to spread from a central point of occurrence to other locations across the brain surface [29][30][31][32][33]. Depth recordings have shown that low voltage spindles may occur in different cortical areas in a local, asynchronous manner, whereas higher voltage spindles may behave as diffuse phenomena [34]. As such, SS are complex events that may be simultaneously measured at different scalp positions and, at the same time, behave like a propagating wave. In recent years SS propagation and frequency modulation (chirp) have been studied in greater detail [18,19,22,27]. Given the diffuse nature of SS, it is relevant to determine whether the different brain regions that participate in global spindles present some level of phase synchronization and how this synchronization relates to dynamic variables associated with SS.
The aim of this work was to extract unambiguous, high voltage, diffuse scalp EEG spindles from human subjects and to investigate spindle dynamics in terms of phase synchronization, chirp and signal propagation, while taking into consideration slow and fast spindle types. A set of independent methodologies showing reciprocally consistent results was used in this work. We believe that this set of procedures (rather than only one methodology) can be useful because of the lack of stationarity inherent to SS segments. We will measure both envelope and oscillatory characteristics to investigate how SS propagate through the brain and discuss how these characteristics are correlated.

Subjects and Sleep Studies
The eight subjects included in the present study were selected from a database that was previously designed for studies on spindle characteristics in obstructive sleep apnea (OSA) [27,35]. Consecutive patients aged 34-60 with clinically suspected OSA [20] were prospectively enrolled for polysomnography (PSG) at a university hospital-based sleep clinic between April 2007 and July 2009 (for additional details see [35]). Subjects provided informed written consent and the study was approved by the local ethics committee (Comitê de Ética do Grupo de Pesquisa e Pós-Graduação do Hospital das Clínicas de Porto Alegre, GPPG HPCA approval number: 100248). Six of these patients were considered to be non-OSA subjects (global Apnea-Hypopnea Index, AHI below 5). Subject 7 was a 54 year-old woman with body mass index 40 and mild lung emphysema. Subject 8 was a 43 year-old male with moderate OSA (AIH 15.9). Patients were taking fluoxetine (1), amitriptyline (1) and several non-psychotropic drugs for non-neurological co-morbidities.

Sleep Spindle Detection and Selection
Spindle detection was carried out with a matching pursuit (MP) program obtained from http:// eeg.pl [36]. MP was used exclusively to identify SS candidates and was applied only in NREM sleep stages II and III. All subsequent analyses were performed on the original time series (without subsampling). MP has been previously described in detail [37,38] and was suitable for sleep spindle representation [22,36,[39][40][41]. After subsampling to 128Hz, each whole-night EEG series was segmented into juxtaposed bins of 2048 digital points and subjected to MP decomposition with a dictionary size of 10 5 atoms, stopping at 96 iterations. Each atom obtained with MP has a central point in time and frequency, and time and frequency full widths at half maximum (FWHM) corresponding to ±σ on a gaussian curve. FWHM duration was used as one parameter for atom selection. Atoms with FWHM duration between 0.5s and 2s and central frequency between 11Hz and 16Hz, hereafter called spindles, were collected in the procedure. Notably, individual MP atom fulfilling detection criteria is not conceptually equivalent to a visual sleep spindle, and the procedure is robust and reliable at a statistical level [41]. The number of atoms that obey these criteria is usually on the average of 2,000/channel/night. Each extracted segment had a 2 s duration, with the highest amplitude SS positioned on the center.
To ensure inclusion of high voltage diffuse spindles while avoiding inclusion of overlapping spindles, only the 100 highest-amplitude, isolated MP atoms were selected for each subject, regardless of electrode position. The electrode containing this atom was defined as the leading electrode.
Segments with artifacts on any channel were discarded. Some EEG segments that still contained multiple superimposed frequency peaks on any given channel were a posteriori discarded (see below in Synchrony Across Channels subsection). In that case, the next highest amplitude atom was chosen in order to reach 100 SS for each subject. Therefore, the total segment number obtained was 760 × 8 channels (subject number 2 contributed only 60 SS due to noise).

Signal Envelope
Initially we used EEG signals corresponding to SS segments,S k orig , from each evaluated channel k (F3, F4, C3, C4, P3, P4, O1, O2) and submittted them to a filter as shown in Fig 1a and 1b. Signals from each channel were then normalized using Eq 1: where the average of SS intervals was taken. The envelope for S k was obtained by numerically estimating |S k (t)| local maxima and generating an interpolated curve, M k (t). The envelope center was defined as max(M k ) and SS duration was defined as the envelope Full Width Half Maximum (FWHM). We also defined a phase, θ j (t), to characterize the oscillatory part of the signal (Eq 2):

Synchrony Across Channels
Detection of SS phase synchrony across different EEG channels was measured using an adaptation of standard methodologies [8]. Our methodology was based on the determination of the Kuramoto Order Parameter (KOP) [42,43], which in this study was calculated using the Eq 3: where θ j are signal phases for each (N = 8) of the channels.
To dynamically characterize the synchronization process we fitted r(t) by Eq 4: where A is the amplitude of the KOP parameter on the window interval; t 0 is the time shift for a given SS; T is the synchronization time; and m is synchronization duration. This functional form was chosen for its convenience. KOP evolution can be illustrated for a SS detected on different channels (Fig 1). In the specific case shown in (Fig 1a-1c), there is an increase in synchrony near 0.5 s that ends approximately one second later, around 1.5 s.
When the r(t) parameter presented as an irregular shape, the corresponding SS was discarded and the subsequent element was chosen. Residual multiple frequency peak SS (possibly representing superimposed SS [14]) were thereby removed from the sample. This procedure ensured inclusion of elements that had high interchannel synchronization duration (high m), which are expected to correspond to diffuse spindles.

Spindle Frequency Analysis
Fast Fourier Transform (FFT) was used in order to determine an average frequency for each 512-point SS segment (0.5 Hz resolution). SS frequency was defined for each channel as the frequency corresponding to the highest power spectrum peak in the 11-16 Hz range. Global SS frequency was the frequency mean for all eight channels. Diffuse spindles were divided into two groups according to average frequency, Slow (f 13 Hz) and Fast (f > 13 Hz). In order to further verify whether the selected SS sample predominantly comprised diffuse or local spindles, a correlation analysis of SS frequency was performed between each of two channels, and considering the original SS segments corresponding to each channel. This procedure was then repeated for randomly chosen segments from each of the two channels. There was high interchannel frequency correlation for original SS segments, which was lost when random segments were analyzed.

Chirp Measurement
Frequency modulation (also called chirp) was calculated by windowed FFT (WFT) applied to the filtered signal. For each signal segment, the 512-point signal was extended by 0.5 s both before and after the original signal. In this 3 s segment, a 2 s moving window was used to multiply the series by a Gaussian function centered on the middle of the window and divided into 13 steps. This Gaussian procedure was chosen in order to best estimate SS central frequency along the time frame. FFT was evaluated for each 2 s segment and the frequency at spectrum peak was analyzed to yield 13 frequency peaks along the time frame. These 13 frequency peaks were used to fit a linear slope that was considered as the spindle chirp value (0.25 Hz/s resolution).
The chirp value distribution for slow and fast spindles was analyzed for each channel (single channel chirp distribution), and the average chirp value for each SS across all channels (channel-averaged chirp distribution) was also determined. A Student's t test was performed to assess the null hypothesis that chirp value distribution had a zero mean.

Signal Propagation
We assumed a very simple model for spindle propagation. We considered that each SS is produced at a single point at positionr k and propagates as a spherical wave with velocity v k . Thus, using the electrode positions and their relative distances, we compared the time delay response (τ) obtained from these electrodes to find the most probable origin and propagation velocity of a given SS.
To investigate the delay τ, initially we considered the cross-correlation function (Eq 5): hðtÞ ¼ τ was then defined by Eq 6 max ðhðtÞÞ ¼ hðtÞ ð 6Þ The delay τ represents the time delay with maximal overlap between both signals. The spherical model allows us to infer time delays for any two channels i and j (Eq 7): where R represents channel positions. These quantities could be directly compared to delays measured experimentally for each SS k, t k i;j . For each physiologically viabler and v we defined: By minimizing with respect to position and velocity, we can infer positionr Ã k and velocity v Ã Channel positions and relative distances were estimated by positioning electrodes in a generic MRI regular head examination (obtained from the BIRN database [44]) using the 10-20 IS standard [45]. As no study subject had an abnormal head size or cranial asymmetry, the distances between the electrodes were assumed to be constant among the different subjects and the error measurement was estimated to be 0.01 m.
A detailed workflow of the proposed methodology for measuring all of the relevant parameters is shown in Fig 2. 3 Results

Spindle frequency, duration, amplitude and chirp characteristics
The frequency distribution (Fig 3a) was bimodal, with a valley between the peaks that was close to the predicted 13 Hz threshold that divided the slow and fast SS groups. On the leading channel, the median frequency of the SS was 13.3 Hz and the interquartile frequency range (IQ) was 0.3 Hz. The median duration of the SS was 0.64 s (0.10 s IQ).
The chirp distribution for slow and fast spindles in averaged channels was also determined (-0.61 Hz/s ±0.49 and -0.17Hz/s ±0.49, respectively) and showed that slow spindles tended to have more negative chirp relative to fast spindles (Fig 3b). When single channel mean chirp values were analyzed, a distinctive antero-posterior gradient could be seen for both groups, with tendency of negative values and the highest p values (indicating the lowest probability of zero mean distribution) appearing in the electrodes that had a more anterior position (Fig 3b).
We also evaluated the correlations among the measured parameters ( Table 2). The most important correlations were chirp × m, chirp × frequency and duration × m. We found no . We can observe that the frequency distribution is a bimodal function. The synchronization time average is larger than spindle duration, this behaviour indicates that the synchronization remains significant even for small signal amplitude. The slow spindles chirp is biased towards negative values while fast spindles distribution are evenly distributed around zero. Finally the amplitudes synchronization are concentrated on values close to 1, meaning that in general global spindles show significant phase locking behaviour. We could not detect any velocity differences for slow and fast spindles.  Inter-channel scatter-plots are displayed in Figs 4, 5 and 6 for frequency, chirp and velocity variables, respectively. The array in these figures forms two triangular matrices, symmetric in relation to the main diagonal (which represents the scatter-plot for the same channel, included just for comparison purposes). We can observe that correlations decrease as we compare more distant channels. The behavior is basically the same for all the considered variables. These results show that these spindles share most of its defining characteristics, they do not only occur at the same time, but are different manifestations of a single phenomenon. This is an indication that global and local spindles have distinct physiology, since most local spindles are not synchronous and are restricted to specific brain areas. The decreasing correlation indicates that the local population of neurons plays a role by making small but measurable perturbations on the detected signal. We do not observed any pattern on these variations, indicating that they might be intrinsically random [46].

Signal propagation and locating spindle generators
Upon calculating the signal propagation velocity distribution we found that the median propagation velocity was 0.12 m/s (0.08 m/s IQ, Fig 3c). By superimposing heatmaps atop MRI tomographic images deposited in the BIRN database, the SS propagation center location can be visualized in axial and sagital views (Fig 7). In these images the SS originates in the thalamic area, which is consistent with previous studies [23,25,47]. Thus, the proposed methodology generates results that are compatible with those obtained using much more sophisticated techniques [48]. However, we found no statistically significant difference in the origins of slow and fast SS spindles (data not shown).

Spindle synchrony across time
We next determined the distributions for spindle synchrony time (T) and synchrony duration (m) (Fig 3d and 3e, respectively). The SS duration as measured by FWHM was also evaluated (Fig 3e), as were the cumulative sums for synchrony amplitude (Fig 3f). The median synchronization time was 0.11 s (0.09 s IQ) and the spindle synchrony was fast and robust over time. The majority of spindles synchronized in less than 0.2 s, which is during the first oscillation. The median synchrony duration was 1.17 s (0.21 s IQ) and the median synchrony amplitude was 0.98 (0.02 IQ). Among all spindles, 80% had synchronization amplitudes that were larger than 90% (Fig 2f).

Discussion
EEG quantification can be challenging because of several factors, including a lack of stationarity, which also makes EEG useful for gathering electrophysiological information. As EEG analysis becomes more accurate and specific, signal tools are being created and their application for EEG is being carefully explored. Since small changes in EEG signals are known to reflect changes in internal brain mechanisms, such analysis tools take on added relevance from a diagnostic perspective. For example, a recent study showed that mice genetically predisposed to epilepsy display a gradient in the frequency modulation (chirp) that uniquely correlates with Synchronization and Propagation of Global Sleep Spindles the development of the epileptic pathology [26]. The source of these signals may lie in how certain neural centers initialize and maintain their oscillation pattern.
In the present study, we used WFT to calculate SS chirp. In this approach the chirp measurement is independent from the MP spindle detection procedure and as such differs from earlier methodologies [22], wherein chirp measurement and spindle detection were performed using a modified MP dictionary that included Gabor chirplet functions. Nevertheless, our present results (slow chirps between 0.61 Hz/s ±0.49 and fast chirps between 0.17 Hz/s ±0.49) are in line with our earlier findings and reinforce the chirp effect [22,27]. Slow SS had the clear tendency to chirp around 0.5 Hz/s, which is consistent with results from our previous study [22], also employing healthy subjects, but studying only a single (central) EEG channel.
In terms of synchronization data, we found that diffuse SS synchronize in 0.11 s, and this process seems to begin before the SS can be visually detected on the scalp. Our results indicate that synchronization could be starting earlier, and thus would have an effective duration that is Heatmap for the origin of slow and fast spindles in axial and sagittal views. MRI images were obtained from the BIRN database [44]. Relative intensity is represented by the color scale bar displayed below. We can observe that in general the spindle origin is thalamic and, on the sagittal view, that slow SS origins have a broader distribution. Synchronization and Propagation of Global Sleep Spindles longer than the measured SS duration. When macroscopic SS are captured on a scalp EEG, synchronization between channels appears to be already in the process of consolidation.
Although synchronization of signals in human brains has been measured previously, to the best of our knowledge this result is the first demonstration that phase synchronization can be measured for a SS. Such behavior should be a direct consequence of the large-scale organization of brain topology, wherein even distant brain regions are connected. Indeed, the existence of these long-distance connections was shown by Hangman et al. in humans [49], as well as in other animals [50].
Mizuhara et al. investigated phase synchronization during a cognitive task and found that different brain regions presented phase lock behavior that lasted for at least 1 second [4,51]. They also found relevant synchronization of both theta and beta waves, and the most relevant locked frequencies ranged between 13 and 17 Hz, which is the same frequency range where we observed SS in this study. Meanwhile, during the performance of attention-related tasks, Gross et al. observed phase synchronizations that fell within the same frequency range [52]. In both studies, the formation of neural networks was mediated through phase locking. The relevant difference between the Gross and Mizuhara studies and our current study is that the networks we detected appeared to be global rather than local. Furthermore, the earlier studies showed that synchronization was functionally correlated with attention and awareness, but our study focused on signals arising during sleep, which suggests that phase locking might play different roles in brain physiology.
The envelope procedure presents advantages when only a few channels are used and the distances between electrodes is large. This approach differs from the methodology presented by O'Reilly and Nielsen [33] in that reflected signal propagation in all available channels, and beyond adjacent channels, can be considered. There are also other advantages to using our procedure, primarily that it is easier to perform in clinics because fewer electrodes are needed (i.e., no dense electrode arrays) and the recommended 256 Hz sampling rate is sufficient. The velocity values presented here are also consistent with those reported in earlier studies [29,30,53]. The spindle dynamics has been subject of debate. Although some studies detected standard trajectories, most reports showed complex wave trajectory formation patterns that were often not reproducible. For instance, in some measurements, a wave could be seen to spiral in one moment, and in the next moment propagate in a straight line [29,30,33,53]. As such, our methodology may be considered very conservative.
Previous work by Manjarrez et al. investigated spindle propagation by considering the signal measured at each channel as a wave packet moving at a local velocity [54]. However, even by applying the simple model proposed here, the local movement that could be observed in the channels would be complex and counter-intuitive. In general, scalp EEG wave propagation is indeed very complex since the signal first emerges at some focal point and then concomitantly spreads along the neural population. At the same time, this signal may be reaching secondary channels, which complicates the determination of the region where the signal originated.
In our case, we investigated isolated spindles with global characteristics. From results in the literature, slow SS would be expected to be more prominent in frontal channels and detected at a lower amplitude by other channels. Similarly, fast SS would likely be more prominent in parietal regions and then propagate to frontal sites. Although there are many studies that have attempted to map these propagation directions, the origin of the spindles is generally accepted to be thalamic. Therefore, we mapped the wave origin by assuming a homogeneous and isotropic system in three dimensions that allowed the verification of whether the origin of SS is congruent with physiological models.
A recent study by Frauscher et al. [55] evaluated SS in epileptic patients, and correlated surface measurements with intracerebral depth measures. Their results suggest that SS have a strong cortical influence and could be considered to be a local cortical event. However, those results may not be comparable to ours, since, as pointed in Nir et al. [46], not all SS are global and here we focused only on global SS. Our work suggests that diffuse SS can be used to characterize the thalamic source in greater detail despite the presence of many local SS. The origin of local SS can be associated with the intrinsic cortical responses such as reverberation. Considering local and global spindles together may lead to different interpretations and results when the effects of wave propagation are the focus of study.
By considering diffuse SS as global phenomena that are observed at different brain regions and yet originate in the thalamus, we may understand propagation of SS using a simple model. Excitation and suppression mechanisms also are likely to play a role in propagating SS, as evidenced by results in patients with apnea or other disorders that showed a tendency for SS to have a more constant frequency (i.e., less chirp) [27]. Recently, we investigated SS in the context of obstructive sleep apnea (OSA), and found a correlation between the degree of apnea and some quantifiers of SS such as frequency and chirp [27,35].
Moreover, the wave propagation velocity of a field potential (which may be indirectly represented by scalp EEG) across the brain is certainly an important indirect measure of environment properties. Thus, precise measurements of the speed of waves propagating in a medium are commonly used to make inferences about the conditions of the medium itself, or about any obstacles and discontinuities that exist in the wave path. Although similar inferences can be made about EEG signals, neuron populations in the cortex oscillate and propagate this oscillation [33], which implies that neighboring neuronal populations help maintain these oscillations. One natural mechanism of synchrony can be measured by the relationship between chirp effects seen among different channels. Frequency modulation, as measured in previous studies [26,27,33], could be correlated to wave propagation along the cortex mass. In any case, defects in neuronal networks could easily be reflected by alterations in both chirp and propagation velocity. Therefore, separating these two effects is important for the investigation of these values as reliable pathology quantifiers.
We have only few electrodes, and using them we are directed approximately to the same place in thalamus. Studies using more electrodes would improve spatial resolution. With our resolution it is not possible to say if very close different networks are being used to generate the two different SS, but it is noteworthy that Slow SS seems to be more diffuse than Fast SS in our analysis. We believe that there are no separated global two regions, but it seems that the networks that generate both Fast and Slow spindle types lie very close.

Conclusion
In this study we found no correlation between wave velocity, synchronization amplitude and frequency modulation of global SS. In a passive medium we would expect constant wave speed, no phase locking and no frequency modulation in the neural substrate. As such, this approximation seems reasonable for investigating the envelope properties of SS, but is inappropriate for investigating the oscillatory behavior of SS, since both frequency and phase are influenced by the inherent complexity of the medium. This complexity is also reflected by the amplitude distribution of scalp SS, implying that some cortex regions are more strongly connected to the thalamic region than others, such that slow spindles and fast spindles might have different effects on different regions. The finding that slow spindles are more prone to chirp is an indication that the medium reaction indeed depends on the transmitted frequencies.
Our study supports the idea that a significant proportion of SS should be considered as a global brain phenomenon, since most characteristics of SS such as frequency, duration and chirp are strongly correlated. Moreover, there is a clear indication that the system as a whole presents high levels of synchrony. Since this collective behaviour is strongly dependent on the complex topology of the neural substrate, this method may prove to be a useful tool for characterizing disease phenotypes in neurodegenerative processes (e.g., Alzheimer's disease) that are expected to influence synaptic strength. Thus, this method should be further developed in future studies.
Supporting Information S1 Data. Sleep spindle characteristics data. ZIP file with all data, figures and the script used in sleep spindle analysis. (ZIP)