The laminar profile of sleep spindles in humans.

Sleep spindles are functionally important NREM sleep EEG oscillations which are generated in thalamocortical, corticothalamic and possibly cortico-cortical circuits. Previous hypotheses suggested that slow and fast spindles or spindles with various spatial extent may be generated in different circuits with various cortical laminar innervation patterns. We used NREM sleep EEG data recorded from four human epileptic patients undergoing presurgical electrophysiological monitoring with subdural electrocorticographic grids (ECoG) and implanted laminar microelectrodes penetrating the cortex (IME). The position of IMEs within cortical layers was confirmed using postsurgical histological reconstructions. Many spindles detected on the IME occurred only in one layer and were absent from the ECoG, but with increasing amplitude simultaneous detection in other layers and on the ECoG became more likely. ECoG spindles were in contrast usually accompanied by IME spindles. Neither IME nor ECoG spindle cortical profiles were strongly associated with sleep spindle frequency or globality. Multiple-unit and single-unit activity during spindles, however, was heterogeneous across spindle types, but also across layers and patients. Our results indicate that extremely local spindles may occur in any cortical layer, but co-occurrence at other locations becomes likelier with increasing amplitude and the relatively large spindles detected on ECoG channels have a stereotypical laminar profile. We found no compelling evidence that different spindle types are associated with different laminar profiles, suggesting that they are generated in cortical and thalamic circuits with similar cortical innervation patterns. Local neuronal activity is a stronger candidate mechanism for driving functional differences between spindles subtypes.

I ntroduction

Sleep spindles are EEG oscillations arising from the interaction of thalamocortical, corticothalamic and reticular thalamic networks (Steriade, 2003; Lüthi, 2013). Sleep spindles are characterized by neural firing patterns which are highly conductive for long-term synaptic changes (Lüthi, 2013). Sleep spindles have been implicated in learning (Gais et al., 2002; Gais and Born, 2004; Clemens et al., 2005), play a role in general cognitive ability (Ujma, 2018), but they are also important clinical markers in both neurological and psychiatric con-ditions (Bódizs et al., 2012; Ferrarelli, 2015; Manoach et al., 2015; Gorgoni et al., 2016; Berencsi et al., 2017) .

Based on the existence of ‘core’ and ‘matrix’ thalamocortical con-nections (Jones, 1998, 2001) – which project to cortical layer IV and I-II, respectively – it has been hypothesized that sleep spindles can arise through either or both of these networks (Piantoni et al., 2016), and spindles from different networks may serve different functions. If different thalamocortical projections terminate in different cortical layers, the thalamic source of an electroencephalographic event can be inferred from sinks and sources in the current source density of signals recorded from electrodes which penetrate the cerebral cortex (Freeman and Nicholson, 1975; Ulbert et al., 2001). Studies in human epileptic patients have successfully revealed the laminar profile of the K-complex (Cash et al., 2009) and the slow wave (Csercsa et al., 2010), both of which were revealed to originate mainly from superficial layers. Results from animal laminar recordings (Spencer and Brookhart, 1961; Kandel and Buzsaki, 1997) reported that various thalamocortical net-works are involved in sleep spindle generation. Recently, a human study (Hagler et al., 2018) found that sleep spindles recorded within the cortex occur with variable topographies.

Sleep spindles differ in at least two characteristics which may be associated with a different network of origin, and consequently a different laminar profile. First, sleep spindles can be divided into function-ally different slow and fast subtypes (Andrillon et al., 2011; Fogel and Smith, 2011), Second, there is substantial variability in the spatial ex-tent (globality) of spindles (Andrillon et al., 2011; Nir et al., 2011; Piantoni et al., 2017).

Some evidence is available about divergent thalamocortical genera-tors of slow and fast spindles or spindle with a different degree of cortical involvement. First, pharmacological manipulations affect slow and fast spindles differently. The Ca2+ channel antagonist flunarazine selectively reduces fast spindle activity, while the voltage-gated Na+ channel antagonist carbamazepine reduces fast spindle activity and enhances both slow spindle and ∼0.75 Hz slow wave activity (Ayoub et al., 2013). Notably, Ca2+ channels are implicated in the thalamic generation of spindles (Astori and Luthi, 2013; Lüthi, 2013), especially in the reticular nucleus (TNR): thus, the resilience of slow spindles to the blocking of these channels suggests that they can be generated independently of the thalamus (Timofeev and Chauvette, 2013). Second, an in silico simulation (Bonjean et al., 2012) demonstrated that more widespread tha-lamocortical connections to superficial layers contribute to the greater spatial extent and cortical synchrony of spindles. It has also been hypothesized that widespread (global) spindles have greater reliance on thalamocortical and possibly cortico-cortical networks terminating in superficial layers (“matrix” network), while fast spindles and local spindles preferentially rely on the more spatially focused “core” network (Piantoni et al., 2016). A recent empirical study (Hagler et al., 2018) in-deed found that spindle-frequency oscillations recorded from the depth of the cortex occur with a heterogeneous topography, but did not tie this to spindle subtypes.

Based on the available evidence, it can be hypothesized that slow spindles and widespread (global) spindles have greater reliance on net-works terminating in superficial layers (“matrix” network), while fast spindles and local spindles preferentially rely on the more spatially focused “core” network (Piantoni et al., 2016). In line with this theory, relatively more superficial activations would be expected in case of slow and global spindles, while the activation should be higher in deeper cortical layers during fast and local spindles.

M aterials and methods

Patients and data selection

Sleep electrophysiological data from 4 patients (1 female, 3 males, age 12–34 years old, detailed data in Table 1) with drug-resistant epilepsy undergoing presurgical electrophysiological monitoring was used. All interventions were approved by the Hungarian Medical Sci-entific Council and the ethical committee of the National Institute of Clinical Neuroscience. Clinical procedures were not biased for scientific purposes. All patients gave informed consent in line with the Declaration of Helsinki. Data from these patients was previously used for research (Cash et al., 2009; Csercsa et al., 2010; Hagler et al., 2018).

Table 1

Basic clinical parameters and ECoG spindle counts by patient. The name of the drugs used in anti-epileptic medication regimens are abbreviated; CBZ: carbamazepine, CLON: clonazepam, LEV: levetiracetam, LTG: lamotrigine, OXC: oxcarbazepine, PHT: phentoin.

	Patient 1	Patient 2	Patient 3	Patient 4
Age	34	15	12	19
Sex	Female	Male	Male	Male
MRI finding	Normal	Right frontal dysgenesis	Right parietal dysgenesis	Right fronto-centro-opercular dysgenesis
Seizure type	Focal, tonic postural	Focal, tonic postural	Focal, sensory-motor hemiconvulsive	Focal hypermotor
Pharmacotherapy	OXC, LEV, CLON	PHT, LTG	CBZ, OXC	CBZ, LTG, PHT
IME location	Right superior frontal gyrus	Right frontal	Right gyrus precentralis	Right gyrus precentralis
Slow spindle frequency range (Hz)	11.45–12.61	10.43–11.54	12.26–13.5	11.28–12.22
Fast spindle frequency range (Hz)	13.3–14.67	12.49–13.38	13.64–15.28	13.49–14.44
Recording duration	16 min	2.66 h	4.43 h	2.22 h
Nslow, global	73	335	34	404
Nslow, local	19	256	97	167
Nf ast, global	76	460	28	324
Nfast, local	19	356	84	241

Patients underwent electrophysiological recordings using implanted laminar microelectrodes (IME, 24 electrodes) and subdural grid and strip electrodes, from which only grids were analyzed (ECoG, 20-64 electrodes) (Figure 1). The selection criteria for the patients in the current study were the following: 1) simultaneous recordings of electrocorticography (ECoG) data from subdural grids and at least one implanted mi-croelectrode (IME) 2) existing post-surgical histological reconstruction of the neural tissue surrounding the IME (electrode track). The selection criteria for EEG data were the following: 1) seizure-free data with clearly identifiable NREM sleep, defined by the presence of sleep spindles and slow waves on ECoG and (if available) scalp EEG channels 2) adequate signal quality of both ECoG and IME data, indicated by the absence of continuous, broad-frequency artifacts. Furthermore, the absence of high-frequency artifacts (verified by visual inspection) and the presence of visible single-unit peaks were the prerequisite for considering the data from a patient for multiple unit (MUA) and single-unit activity (SUA) analysis, respectively. For the simultaneous analysis of ECoG EEG and IME MUA data in Patient 1-2, only ECoG spindle data concomitant to high-quality IME signal was used.

Fig. 1.

Panel A: ECoG and IME implantation procedure. Panel B: schematic representation of the localization of IME contacts with postoperative histology. The histological reconstruction of the tissue surrounding the IME in Patient 1 is shown. Panel C: postoperative MR reconstructions of the location of ECoG and IME electrodes from all patients.

Electrophysiology and electrode positioning

IMEs had a diameter of 350 μm, inserted perpendicularly to the cortical surface, penetrating the cortex to a depth of 3.5 mm. 40 μm platinum/iridium electrodes were built into the IME, spaced evenly at 150 μm (Ulbert et al., 2001). A silicone sheet at the top of the mi-croelectrode prevented the array from shifting below the pial surface. ECoG hardware filters were set to 0.1–200 Hz (ECoG). ECoG data was recorded with a contralateral mastoid reference. IME data was recorded with a bipolar reference, with each contact referred to the one inferior to it (local field potential gradient, LFPg). IME data was separated into the EEG range (0.1–300 Hz) and the single/multiple unit range (300–5000 Hz) at the level of the amplifier and recorded as two separate data files. ECoG was recorded with a sampling frequency/precision of either 2000 Hz/16 bit (Patient 1) or 1024 Hz/16 bit (Patients 2, 3 and 4), while IME data was recorded with a sampling frequency/precision of 20 kHz/12 bit. ECoG and microelectrode signals were recorded separately and synchronized post-hoc based on an automatic synchronization signal applied to the deepest IME contact (which was henceforth excluded from further analyses) and a randomly selected ECoG channel. This signal appeared at regular intervals (usually 10 s), its amplitude exceeded any physiological activity on the channel by several orders of magnitude and it encoded the precise date and time. A semi-automatic algorithm was used to detect these synchronization signals and find corresponding pairs in ECoG and microelectrode recordings with data point precision.

LFPg was equivalent to the data originally recorded from the IME, corresponding to the voltage difference between neighboring IME channels at a distance of 150 μm (μV/150 μm, generally referred to as μV for simplicity). Current source density (CSD) was calculated as the negative of the second spatial derivative of the field potentials (Csercsa et al., 2010). However, since LFPg already constituted the first spatial derivative, CSD was calculated by calculating the negative of the first spatial derivative of the LFPg. Due to differences in the physical properties of the tissue surrounding the electrode environment iwe expected a different signal voltage in each patient. Therefore, we generally z-transformed signal amplitude within patient for statistical analyses and indicate when it was otherwise. Unstandardized signal voltage, when reported, is not directly comparable across patients.

Electrode positions were confirmed by intraoperative navigation, the comparison of pre- and postoperative MR scans, as well as the comparison of photographs taken during the initial surgery and the ones taken during resective surgery. The brain tissue containing the IME was re-moved during surgery, cut to 2–5 mm blocks and chemically fixated (Ulbert et al., 2004; Csercsa et al., 2010). The laminar topography of the cortical tissue blocks was reconstructed from these samples in all patients and it was used to determine the location of IME electrodes within cortical layers, taking into account the shrinkage of brain tissue during preparation (Wittner et al., 2006; Csercsa et al., 2010).

The cortical layer from which an IME signal originated was deter-mined based on the histological reconstruction of the electrode track. Postoperative examinations confirmed that the IME was implanted out-side the seizure onset zone in Patients 1–3, but within the seizure onset zone in Patient 4.

Sleep scoring and artifact rejection

Hypnograms were visually scored for ECoG data on a 20 s basis based on standard criteria (Iber et al., 2007). Since the standard scoring cri-teria are generally only applicable to scalp EEG channels with a full polysomnography setup (including EOG and EMG), we restricted our hypnograms to the identification of NREM sleep (regardless of stage) and the separation of it from other sleep states and wakefulness, based on the presence of slow waves and spindles. Artifacts were excluded from ECoG data on a 4 s basis using visual inspection. Only artifact-free data from NREM sleep was considered for further analysis.

Sleep spindle detection

ECoG data was analyzed using the Individual Adjustment Method (IAM) (Bόdizs et al., 2009; Ujma et al., 2015) in order to identify sleep spindles. The IAM essentially defines sleep spindles as waveforms with a specific frequency and sufficient amplitude that contributes to the spectral peaks of the NREM sleep EEG. In the first step, the high-resolution (0.0625 Hz) amplitude spectrum (AS) of the visually scored, artifact-free NREM sleep EEG was calculated for each ECoG channel free of interictal discharges (IIDs, “spikes”) based on 4 s segments (Hanning-tapered, zero padded to 16 s in order to ensure sufficient resolution in the fre-quency domain). The second-order derivatives of the averaged AS were calculated in order to identify spectral peaks. The resulting second-order derivatives were averaged across all IID-free ECoG channels. (This was a deviation from the original IAM methodology which averages second-order derivatives across frontal and centro-parietal scalp channels separately. This change was necessary because ECoG arrays did not always contain both frontal and centro-parietal channels.) A spectral peak was defined as the interval within which second-order AS derivatives were below zero. In line with the assumptions of the IAM about the spectral characteristics of sleep spindles, two clear spectral peaks were visible for each patient, one for slow spindles and another for fast spindles. The areas between these spectral peak boundaries across the frequency domain were defined as the slow and fast sleep spindle frequency range of the patient, respectively. In the second step, ECoG data was filtered to the patient’s individual slow or fast sleep spindle frequency range, and sleep spindles were identified as events in which the amplitude exceeded an electrode-specific threshold defined as the average of the AS values at the spectral peak boundaries, multiplied by the number of high-resolution frequency bins within the frequency range. If this dynamically defined threshold was not exceeded for at least 0.5 s, no sleep spindle was detected. A detailed illustration of the IAM is available on Fig. 1 in a previous paper (Ujma et al., 2015). Sample detections and average spindle waveforms are shown on Supplementary figure S1.

Sleep spindles were detected in this manner for IID-free ECoG channels and IME channels with no permanent artifacts. An ECoG spindle was considered as “local” if it was detected in less than or equal to 33% of the IID-free ECoG channels, and “global” if it occurred in more than 33% of IID-free ECoG channels. This cutoff point was chosen to ensure a similar number of local and global spindles (Table 1), but analyses confirmed the similarity of the laminar profile of all spindles (see Results), rendering the precise definition of the local-global cutoff point irrelevant.

For the analysis of the laminar profile of ECoG spindles, ECoG data from the channel spatially closest to the IME was filtered to the individual slow and fast spindle frequency, respectively, as it was determined in the first step of the IAM (two-way FIR filtering using the eegfilt() function in MATLAB [The MathWorks, Inc., Natick, Massachusetts] EEGLab [Swartz Centerl for Computational Neuroscience, San Diego, California]). IME data was triggered to the maximum positive deflection of the spindle cycle closest to the point where the envelope (defined as the modulus of the Hilbert transform of the filtered signal) of the filtered signal of this spatially closest ECoG electrode was maximal (spindle peak). IME LFPg data segments from −500 ms to 500 ms relative to the identified ECoG spindle peaks were selected for further analysis using Neuroscan (Compumedics USA, Charlotte, NC). All selected data segments were visually inspected again in order to identify artifacts which were not present in the ECoG data. Artefactual data segments were excluded. All data segments were filtered to a generic sleep spindle range (10–16 Hz) and smoothed across layers using a Hamming window. CSD was calculated by data segment after selection, filtering and smoothing. In order to quantify laminar differences in field potential and current strength, respectively, the root mean square (RMS) of the LFPg and CSD signal of each channel was calculated for each 1-second IME data segment (referred to as LFPg/CSD magnitude). Channel-wise magnitudes were averaged across channels in the same cortical layer, yielding an average magnitude of LFPg and CSD signals for each spindle and for each cortical layer.

In co-occurrence analyses we considered two spindle events from two electrophysiological sources s1 and s2 to co-occur if a sleep spindle detected in s2 started before a concomitant sleep spindle in s1 terminated, s1 always indicating the channel where the earlier spindle was detected. For cortical layers, only spindles co-occurring on at least 1/3 of the channels within the cortical layer using the above criteria were considered in order to control for the higher prior probability of a sleep spindle occurring at a given time (and thus co-occurring with others) as a function of the number of channels within the same layer. As for the co-occurrence between IME and ECoG, a sleep spindle was considered to co-occur if it occurred simultaneously using the above criteria in any cortical layer and an ECoG electrode.

Multiple unit and single-unit activity

Multiple unit activity (MUA) was estimated by calculating the filtered (300-3000 Hz, zero phase shift, 48 dB/octave) and rectified IME LFPg sampled at 20 kHz, and passing it through a low-pass filter (20 Hz, 12 dB/octave). MUA was triggered to ECoG spindle detections in a similar manner as LFPg data from the EEG range and averaged across detections by spindle type (slow/fast, local/global). MUA was only computed from the two patients (Patient 1 and 2) with adequate high sampling rate signal quality.

Single-unit activity (SUA) was detected in high sampling rate IME recordings (20 kHz) in case of two patients (Patient 1 and 2), where individual cell activities could be separated reliably from MUA/field ac-tivity and background noise. SUA was detected in epoched raw IME data sampled at 20 kHz, initially extracting two-second epochs (±1000 ms around all ECoG spindle peaks). We epoched and analyzed data from slow, fast, local and global spindles separately. After DC offset removal, we detected SUA based with an amplitude threshold adjusted manually according to the magnitude of background noise of each channels. Multiple individual neurons were identified as the generators of SUA on each channel based on clustering by action potential morphology and amplitude in a 0.4 msec timeframe. We applied a principal component analysis based on waveform characteristics to refine clustering and reduce false detections. All SUA detections were performed in Spike2 (version 7 software (CED Limited, UK). SUA occurrence relative to LFPg phases was calculated based on the instantaneous phase angle of the Hilbert transform of the LFPg signal. All SUA detections within ±500 msec from ECoG spindles of the same spindle type were pooled from all cells on the same channel and from all IME channels in the same layer, while those outside this time range relative ECoG spindle peaks were discarded.

Experimental design and statistical analysis

We performed statistical analyses using STATISTICA 10.0 (StatSoft, Tulsa, OK) and IBM SPSS Statistics 20 (IBM Corporation, Armonk, New York). The statistical analyses of MUA and SUA were performed using Oriana (v2.02c, Kovach Computing Services) and the CircStat toolbox (Berens, 2009) for MATLAB. We used α 0.05 for all analyses. Error bars on all plots illustrate the 95% confidence interval of the mean. The sample size for statistical tests was always based on the number of spindles, not the number of patients, and it is indicated for each analysis.

We used parametric statistical tests in STATISTICA (Pearson’s point-moment correlation coefficients and one-way ANOVA, as indicated in the text) for the calculations illustrated on Fig. 3-5 and in Table 2. We used a X2 test in SPSS to assess the homogeneity of slow and fast spindle counts across layers. We used a repeated-measures ANOVA to assess the effects of Layer, Spindle type and Globality on signal magnitude across layers (Fig. 6 and 7), treating Layer as a 5-level within-subject (more precisely, within-spindle) variable. Due to large between-patient differences in signal amplitude and the resulting violation of the assumption of homoscedasticity, we repeated this test in each patient. For z-transformed magnitudes where amplitude differences were normalized, we only ran one test with Patient as a covariate. We ran a PCA in STATISTICA to determine the multicollinearity of LFPg and CSD magnitude across layers (Fig. 8). It must be noted that our application PCA to IME signal magnitude during ECoG spindles is not directly comparable to a previously used method (Hagler et al., 2018) applying it to the sigma-range IME signal itself. Our method specifically investigated the similarity of the IME signal during ECoG spindles, while the other method estimated the similarity of spindle-frequency IME activity. Circular statistics concerning the phase coupling of LFPg and SUA/MUA (Fig. 9 and 10) were performed in MATLAB using the CircStat toolbox (Harrison-Kanji test) and Oriana (v2.02c, Kovach Computing Services) (Rayleigh and Mardia tests).

Fig. 3.

Panel A: sleep spindle amplitude on the IME. For each layer, sleep spindle amplitude is given as the mean value of all channels in the layer. Panel B: sleep spindle density on the IME. For each layer, sleep spindle density is given as the mean value of all channels in the layer. Panel C: Histogram of the spatial extent of sleep spindles on the IME, given as the number of spindles (axis y) simultaneously detected in a given number of cortical layers (axis x). Only spindles occurring on at least 1/3 of the channels in the same cortical layers were considered.The four histograms are not significantly different (2-sample Kolmogorov test, pm^ = 0.69 from all possible comparisons). Panel D: histogram of the spatial extent of sleep spindles on the ECoG, given as the number of spindles (axis y) co-occurring on a given number of ECoG channels (axis x). The first 20 bins (this number is limited by the lower number of ECoG electrodes in Patient 2) of the four histograms are not significantly different (2-sample Kolmogorov test, pmin = 0.27 from all possible comparisons). Panel E: the association between ECoG sleep spindle amplitude (axis y) and extent defined as the number of ECoG electrodes with simultaneous detections (axis x). All correlations are statistically significant after correcting for multiple comparisons using Bonferroni corrections. Data is shown for the ECoG channel closest to the IME. N = 88, N = 94 (patient 1) and N = 78, N = 110 (patient 2) for slow and fast spindles, respectively. Correlation coefficients are not significantly different (Fisher’s r-to-z method, pmin = 0.23 from all possible comparisons). Panel F: the association between IME sleep spindle amplitude (axis y) and extent defined as the number of IME electrodes with simultaneous detections (axis x). All correlations are statistically significant after correcting for multiple comparisons using Bonferroni corrections. Data is shown from a representative channel with high spindle density. N = 74, N = 94 (patient 1) and N = 107, N = 108 (patient 2) for slow and fast spindles, respectively. In Patient 2, slow spindle amplitude is nominally more strongly correlated with extent than fast spindle amplitude (p = 0.044), but this effect does not survive Bonferroni correction. No other correlation differences are significant (Fisher’s r-to-z method, pmin = 0.13 from all other possible comparisons).

Fig. 5.

The extent (number of ECoG channels a spindle was detected on) of ECoG spindles co-occurring with IME spindles detected from each cortical layer in Patient 1 (left) and Patient 2 (right). Error bars denote 95% confidence intervals. Slow and fast spindles were pooled for this analysis.

Table 2

The correlation between spindle amplitude and spatial extent. The table shows Pearson’s correlation coefficients between the amplitude of IME sleep spindles and the spatial extent (expressed in the num-ber of ECoG electrodes the spindle was simultaneously detected on) of concomitant ECoG spindles. Results are presented by the cortical layer of IME spindle detection, patient and spindle type.

Patient 1
		Layer I	Layer II	Layer III	Layer IV	Layer V
All spindles	Pearson’s r	0.194	0.319	0.345	0.276	0.276
	N	110	126	140	127	135
	p	0.004	<0.001	<0.001	<0.001	< 0.001
Slow spindles	Pearson’s r	0.195	0.266	0.343	0.261	0.297
	N	45	53	63	59	62
	p	0.020	0.006	0.001	0.005	0.002
Fast spindles	Pearson’s r	0.261	0.403	0.389	0.373	0.334
	N	65	73	77	68	73
	p	0.004	0.000	0.000	0.000	0.000
Patient 2
		Layer I	Layer II	Layer III	Layer IV	Layer V
All spindles	Pearson’s r	0.419	0.467	0.457	0.427	0.476
	N	75	96	112	95	139
	p	<0.001	<0.001	<0.001	<0.001	< 0.001
Slow spindles	Pearson’s r	0.348	0.415	0.311	0.430	0.479
	N	42	54	57	54	65
	p	0.002	<0.001	0.002	<0.001	< 0.001
Fast spindles	Pearson’s r	0.358	0.321	0.316	0.277	0.107
	N	33	42	55	41	74
	p	0.004	0.004	0.002	0.008	0.036

Fig. 6.

Panel A: The laminar profile of slow, fast, local and global spindles, including LFPg, CSD and MUA, in representative Patient 2. MUA amplitude refers to the rectified 300–3000 Hz activity filtered through a 20 Hz lowpass filter. Color maps illustrate the average amplitude fluctuations ±100 msec before and after ECoG spindle peaks. Layer centroids are marked with Roman numerals. Black lines show average traces across all channels. Note the inverted color coding of CSD values. Panel B: Grand average LFPg and CSD magnitude across layers for slow and fast spindles from all patients. Values are expressed as z-scores across layers. Error bars represent 95% confidence intervals.

Fig. 7.

Normalized LFPg (Panel A and C) and CSD (Panel B and D) laminar profiles of sleep spindles. Panel A and B illustrate the mean LFPg and CSD magnitude of each patient in each layer, respectively. Panel C shows LFPg magnitude in each layer for each patient by spindle type (local, global, slow and fast), while Panel D shows the same statistics for CSD magnitudes. Error bars represent 95% confidence intervals on all plots. Note the similarity of the individual laminar profile across spindle types.

Fig. 8.

Scatterplot of LFPg (Panel A) and CSD (Panel B) magnitudes in each cortical layer during ECoG spindles. The three spatial dimensions are used to illustrate magnitudes in layer I-III, while marker sizes and colors illustrate magnitudes in layer IV and V/VI, respectively. Marker sizes and marker colors are transformed to nonlinear scales to optimize visibility. Supplementary video 1 (LFPg) and Supplementary video S2 (CSD) contain a more illustrative, rotating version of Panel A and B, respectively.

Fig. 9.

Panel A: the preferred phases of SUA during various spindle types, superimposed to a stereotypical spindle-frequency oscillation. Only those spindle types are shown for which the LFPg-SUA coupling (based on an at least nominally significant Rayleigh’s Z-test for circular uniformity) was significant. Spindle type is coded by marker colors and layer of occurrence is coded by marker shapes. Panel B: an illustration of the occurrence of SUA during positive-negative LFPg phase transitions using data from fast global spindles and a layer III channel in Patient 1. The rose plot illustrates the preferred firing phase slightly after 90° (see also the first data marker on panel A). The line-histogram shows the frequency of occurrence of SUA in temporal bins superimposed to the mean spindle LFPg on the same channel.

Fig. 10.

Representative LFPg and MUA averages, triggered to ECoG spindle events (Panel A and C) from channel 10 (Layer III) in Patient 1 and channel 6 (Layer II) in Patient 2. Note the fluctuations of the difference of instantaneous phase differences. Panels B and D illustrate the similarity of the LFPg and MUA signals as well as their phase difference by layer and spindle type in Patient 1 and 2, respectively. On the compass plots the length of the arrows corresponds to the modulus of the LFPg-MUA correlation coefficient, while the orientation of the arrows corresponds to the circular mean of the difference of the instantaneous phase of the two signals.

Data and code availability statement

Data and custom MATLAB scripts are available upon request from the Corresponding Author, Péter P. Ujma (peteru88@gmail.com).

Results

Spindle occurrence on the IME

Transient, waxing and waning oscillations with morphological characteristics resembling EEG and ECoG sleep spindles were only visible in the IME recordings of 2 out of our 4 patients (Patient 1 and 2), with an identical frequency as ECoG spindles. However, even in patients without visible IME spindles, once IME recordings were triggered to the peak of the spindles detected on the nearest ECoG electrode, a spindle-frequency oscillation emerged in the averaged signal, indicating that a sub-threshold spindle-frequency signal was always concomitant to ECoG spindles (Fig. 2). Due to the absence of spectral peaks necessary for IAM spindle detection, we only used data from Patient 1 and 2 to detect sleep spindles directly on the IME. However, due to the presence of a spindle-frequency average signal even in case of Patient 3 and 4, we also used IME data during ECoG spindles from these patients.

Fig. 2.

Sleep spindles are not always visible on IME LFPg recordings and do not always create characteristic spectral peaks, but they always emerge when ECoG spindle detection-triggered IME LFPg data is averaged. In Patient 1, both slow and fast spindles visible on ECoG recordings are also visible in IME derivations (Panel A), and spectral analysis reveals that the spectral peaks identified from ECoG are also present in IME data (Panel C, showing power spectral density from each IME channel after channel-wise z-transformation). Blue and red vertical lines represent the slow and fast sleep spindle frequency ranges, respectively, as determined from ECoG data). Averaged IME LFPg triggered to ECoG spindle detections also reveals the presence of spindles on the IME during ECoG spindles (Panel D). In Patient 4, however, sleep spindles are not visible in the raw IME data during ECoG spindles (Panel B), do not leave characteristic spectral peaks (Panel E), but appear when averaged IME LFPg is triggered to ECoG spindle detections (Panel F).

IME spindles had the highest amplitude in the superficial layers (Fig. 3, Panel A), but their density was relatively uniform across layers (Fig. 3, Panel B). Both ECoG and IME spindles often only occurred in a few layers and more global spindles were less common (Fig. 3 C).

Spindles co-occur between sources as a function of amplitude and proximity

A large proportion of sleep spindles co-occurred between cortical layers (typically 50–80%). ECoG spindles were usually (>90%) also detected in at least one cortical layer. It was more common for ECoG spindles to be also detectable on the IME than vice versa (Fig. 4). In other words, many but not all sleep spindles were co-registered by multiple electrophysiological channels, and while there were many spindles which were only detected on the IME or within a specific cortical layer, if a spindle was detected on an ECoG source it was usually also detected on the IME. Detailed co-occurrence patterns are reported in the form of Venn diagrams on Supplementary figure S2.

Fig. 4.

Panel A: the proportion of ECoG spindles with a concomitant IME spindle in at least one cortical layer. The 6 × 8 and 4 × 5 patterns are topographically accurate representations of the implanted subdural electrode grid with a white dot indicating the location of the IME. Electrodes with unusable or pathological data are colored pure white. Panel B: the number and proportion of ECoG spindles co-occurring between specific cortical layers mapped by the IME and the nearest ECoG channel in Patient 1 and Patient 2. Only sleep spindles simultaneously detected on at least 1/3 of the IME channels within the same cortical layer were considered. The diagonals indicate the number of spindles occurring in each cortical layer or on the nearest ECoG channel. The upper half of the matrix shows the proportion of spindles in the source indicated on axis y co-occurring with spindles in the source indicated on axis x, while the lower half shows the proportion of spindles in the source indicated on axis x co-occurring with spindles in the source indicated on axis y. An example of this interpretation is provided on the figure. Panel C, left: the negative relationship between source distance (in layers) and the probabilities of co-occurrence depicted in Panel B (Pearson’s r(80)=−0.31, p = 0.005). All spindle types from both patients were pooled for this analysis, within-patient and within-type correlations (N = 40) range between −0.28 and −0.37. Panel C, middle: The probability of ECoG spindles co-occurring on the IME within a cortical layer (ECoG on IME) is higher than the probability of IME spindles within a cortical layer co-occurring on the ECoG (IME on ECoG) (One-way ANOVA F(1,38)=22.65, p < 0.001). The squares and error bars show means and 95% confidence intervals, individual data points (co-occurrence probabilities between ECoG and each layer from Panel B) are shown. Data from all layers, patients and spindle types were pooled for this analysis. Panel C, right: individual IME spindles that co-occur on ECoG have higher amplitude than those that do not (One-way ANOVA F(1,4133)=166.6,p < 0.001). Data from all layers, patients and spindle types were pooled for this analysis, amplitudes are z-transformed by patient, layer and spindle type. Differences marked by an asterisk are statistically significant. The squares and error bars show means and 95% confidence intervals, individual data points except 30 spindles with standardized amplitude > 4 are shown. This pattern was present in all layers (F(1,640–1193)=13.889–48.191, p < 0.001 in all cases).

Amplitude was a strong moderator of co-occurrence probability. Larger sleep spindles were detected on more IME or ECoG channels (Fig. 3, Panel D-E), larger IME spindles were also more likely to be detected on the ECoG (Fig. 4, Panel C), and IME spindles with greater amplitude co-occurred with more global ECoG spindles (Table 2).

The laminar profile of different spindle types is similar

In case of IME spindles, we found little evidence that different spindle types preferentially occur in different cortical layers. The number of slow and fast spindles was not significantly different across layers, (Patient 1: X2 (4,N = 1911) = 1.588, p = 0.811; Patient 2: X2 (4,N = 2224) = 3.807, p = 0.433, p> 0.05 for all row-wise z-tests comparing slow and fast proportions in individual layers, see also Fig. 3, Panel B). While we could not assess the globality of IME spindles directly, we found no evidence that IME spindles in either layer co-occur with more global ECoG spindles (Patient 1: F(4,633)=0.328, p = 0.859; Patient 2: F(4,512)=1.349, p = 0.251; Fig. 5).

Comparing IME LFPg and CSD signal magnitude in different layers during ECoG spindles, we found only weak evidence that ECoG spindle types are characterized by different laminar profiles. Within the same patient, very similar magnitude patterns across layers in different spindle types were present (Fig. 6-7), suggesting topographically similar generating networks for all spindle types. While the laminar profile was not fully homogeneous across patients, it was always characterized by maximal LFPg and CSD magnitude in the superficial layers and lower magnitude in the Layer V (Fig. 7). We used a repeated measures ANOVA with Magnitude as the dependent variable, spindle type (slow or fast) and spindle globality (local or global) as between-subject factors and Layer as a within-subject factor to estimate spindle subtype differences in the laminar profile. A significant interaction between Layer and either Spindle type or Globality would indicate that the distribution of signal magnitude across layers is different in different spindle subtypes. We repeated this analysis using z-transformed (mean magnitude=0 and SD=1 for each spindle across layers) magnitudes and treating Patient as a between-subject factor, and also by patient using raw magnitude. Detailed statistics are reported in Supplementary table S1. Layer* Spindle type interactions were significant in Patient 2 and the all-patients model for LFPg and Patients 2 and 3 and the all-patients model for CSD. Layer* Globality interactions were significant in Patients 1-3 for LFPg and Patients 1-3 and the all-patients model for CSD. The effect sizes of these interactions were small (partial η2 η2CSD, all patients=0.002; maximal partial η2LFPg=0.012, maximal partial η2CSD=0.068, both in Patient 2), indicating that while the distribution of LFPg and CSD magnitude is not perfectly uniform across layers in each spindle subtype, these subtype effects are comparatively small. Following a reviewer’s suggestion, we re-run our analyses of the z-transformed, pooled spindle sample using linear mixed models (LMMs) with patient-wise random effects to better account for dependence of observations resulting from sampling multiple spindles from each patient. This analysis was implemented in the lme4 R package. This confirmed that including spindle subtype effects has a minor but sometimes significant effect (pLFPg=<0.001-1; pCSD=<0.001-1) on model fit, but the largest influence was attributable to the patient-dependent effect of layer (pLFPg < 0.001, pCSD <0.001). Details about LMMs and model fit are reported in Supplementary table S2.

The laminar profile is homogeneous across individual spindles

Theoretically, it could be the case that the average laminar profile is the aggregate of many heterogeneous laminar profiles during individual ECoG sleep spindles: that is, individual spindles have various laminar profiles, but when averaged, these converge on a stereotypical pattern which is not representative for many or most individual spindles. In order to test this hypothesis, we z-transformed both LFPg and CSD magnitude (root-mean square of the amplitude of the 1-second wave-form around the ECoG spindle peak) within each layer of each patient across all detections to eliminate between-layer and between-patient voltage differences, and performed a principal component analysis on the magnitude values in each layer. Both LFPg and CSD magnitudes were strongly correlated across layers: if a spindle was characterized by low magnitude in one layer, it typically also had low magnitude in other layers (and vice versa) (Fig. 8). In each patient, only one principal component with an eigenvalue >1 emerged, explaining 62–82% of the variance in LFPg magnitude and 60–76% of the variance of CSD magnitude across layers.

We next used cluster analysis in order to investigate whether sub-groups of spindles exist based on their concomitant laminar activation, even if these subgroups do not specifically consist of slow, fast, global or local spindles. We used SPSS’s TwoStep clustering algorithm to assess how well ECoG spindles can be assigned to clusters based on the laminar distribution of LFPg or CSD magnitude. We ran this analysis separately by patient in order to avoid creating clusters based on patients and also in an attempt to replicate cluster structures across patients. We also repeatedly ran this analysis for LFPg and CSD activations, always using z-transformed values (mean magnitude = 0 and SD = 1 for each spindle across layers). Neither the number of clusters in the best solution (2 or 3) nor the laminar activation patterns of the spindles belonging to each cluster were consistent across patients (Supplementary figure S3), and the average silhouette coefficients (<0.5 in all cases) suggested a less than ideal separation of clusters.

In sum, all IME and ECoG spindle types (slow/fast, local/global) were characterized by a similar, stable pattern of signal amplitude in cortical layers, making it unlikely that neuronal networks with substantially different cortical innervation patterns contribute to different spindle types or individual spindles.

Single-unit activity is coupled to LFPg

High-quality single-unit (SUA) and multi-unit (MUA) activity could be obtained from Patient 1 and 2. Based on a statistically significant deviation from circular uniformity using Rayleigh’s test, in Patient 1 SUA was significantly coupled to LFPg during all spindle types and in all layers except slow local, fast local and slow global spindles in Layer II and during fast local spindles in Layer V. In Patient 2, both slow and fast global spindles were coupled to LFPg in Layer II and III. Slow local spin-dles in Layer II and fast global spindles in Layer VI were also significantly coupled. In case of significant coupling, SUA preferentially occurred during positive-negative LFPg transitions in all patients (Fig. 9). The phase relationship between SUA and LFPg by patient, layer and spindle type and average SUA waveforms is reported in the Extended data 9-1.

A circular factorial ANOVA using the Harrison-Kanji test revealed a significant omnibus effect of both layer and spindle type on the preferred phase of SUA-LFPg coupling in both patients (Patient 1: X2 layer (6,N = 41756)= 158.89, X2type (6, N = 41756)= 149.63; Patient 2: X2 layer (8,N = 19496)= 106.82, X2 type (6, N = 19496)= 32.95, all p< 0.001), with a significant interaction in Patient 1 (X2 (9,N = 41756)= 22.69, p = 0.031), but not in Patient 2 (X2 (12,N = 19496)= 9.83, p = 0.364). As post-hoc tests, we compared the preferred phases of SUA- LFPg coupling across layers and across spindle types using Mardia’s test of the difference of circular means. In Patient 1, slow and fast spindles within the same layer typically had similar preferred phases, which were however different for the same spindle types in different layers, and for local and global spindles even in the same layer. Local spindles prefer-entially occurred at a later phase of the spindle. In Patient 2, the most prominent effect was the later preferred phase of layer V-VI compared to superficial global spindles. Detailed statistics are provided in the Ex-tended data 9–1.

Multiple-unit activity

Multiple-unit activity (MUA) of adequate quality was similarly avail-able from Patient 1 and 2. In order to assess its coupling to LFPg, both MUA and LFPg from each IME channel was triggered to ECoG spindle events and averaged across events. The similarity of the resulting aver-age signals was assessed by computing their Pearson’s correlation coef-ficient, and their phase difference was estimated by the angular mean of the difference of the phase angle of their Hilbert transforms. LFPg and MUA during spindle events were generally highly similar and exhib-ited an approximately antiphase (90°<ΔΦ<270°) relationship (Fig. 10). However, the LFPg-MUA phase angle difference was not constant dur-ing spindle events, especially in Patient 2, leading to a downward bias of correlation coefficients and the poor representativeness of the average phase difference.

In sum, SUA and MUA was generally maximal during LFPg nega-tive phases or positive-negative transitions, but this relationship was not fully homogeneous across patients, layers and spindle types.

Discussion

This study analyzed co-detected sleep spindles from both the cortical surface (ECoG) and from individual cortical layers (IME). Our findings highlight several characteristics of sleep spindle occurrence which are important for the understanding of both the generation of sleep spindles and the functional differences between sleep spindle types.

First, we found that sleep spindles occur in all electrophysiological derivations and often co-occur between them. Co-occurrence, however, is neither ubiquitous nor random. Sleep spindles detected in cortical lay-ers are often local events (Fig. 3, Panel C), localized to a single cortical layer. They co-occur in other cortical layers about 50–80% of all cases, and in the nearest ECoG derivation about one-third of the time (Fig. 4). Co-occurrence is more probable if the layers in question are close to each other or if the IME spindle has higher amplitude. ECoG spindles, on the other hand, co-occur in at least one cortical layer about 90% of the time. Previous research has demonstrated that larger sleep spindles appear on more derivations (Andrillon et al., 2011; Nir et al., 2011) or that only about 50–75% of sleep spindles co-occur between cortical layers (Hagler et al., 2018). Our study, however, is the first to extend this finding to a combination of subdural and intracortical electrodes, showing that these co-occurrence patterns are general phenomena. The asymmetrical co-occurrence patterns between IME and ECoG is reminis-cent of the observation that spindle oscillations in electrophysiological channels with larger recording areas (EEG, especially monopolar EEG) are usually also seen in channels with smaller recording areas (MEG, especially gradiometers) but the reverse is not true (Dehghani et al., 2010; Dehghani et al., 2011): that is, a large degree of sleep spindles are extremely local events and monopolar EEG derivations only detect their largest, most synchronized instances. IME data are essentially bipo-lar EEG recordings from inside the cerebral cortex, with each channel having an extremely small recording area measurable in hundreds of mi-crometers at most (Ulbert et al., 2001). Therefore, IME-exclusive spin-dles likely represent the extreme left tail of the distribution of the spatial extent of sleep spindles, encompassing the most local instances.

Second, IME spindles have highly variable laminar topographies and may occur exclusively in any layer, but EEG activity recorded from the IME during ECoG spindles is highly stereotypical and shows an ampli-tude distribution which is very stable at least within the same patient (Fig. 7). Our results are in line with a model of sleep spindle generation in which sleep spindles may be generated in multiple thalamocortical (or possibly cortico-cortical) networks with various cortical innervation characteristics, but only extremely local spindles likely invisible for the EEG or ECoG (Dehghani et al., 2010; Dehghani et al., 2011) rely on a sin-gle well-delineated generator. For larger spindles, the joint involvement of multiple networks is the norm, resulting in co-detections in multiple layers and a more stereotypical laminar profile.

Third, we found little evidence that sleep spindle subtypes - slow and fast spindles, or spindles with variable cortical involvement (local and global spindles) - were characterized by a sharply different lami-nar profile, as previously hypothesized (Timofeev and Chauvette, 2013; Piantoni et al., 2016). The occurrence of slow and fast spindles in dif-ferent cortical layers was not significantly different, and IME spindles from different layers did not co-occur with ECoG spindles of a different degree of globality. The laminar profile of ECoG spindles was not per-fectly homogeneous, but the effect of spindle subtype was small (par-tial η = 0–0.06) and non-homogeneous across participants. We did not observe the previously hypothesized (Timofeev and Chauvette, 2013; Piantoni et al., 2016), sharply different cortical profiles for different spindle subtypes. These results suggest that while sleep spindles may be generated in multiple thalamocortical networks with various cortical innervation patterns (Hagler et al., 2018), these networks are not en-gaged in a substantially different manner during the generation of dif-ferent spindle subtypes. Regardless of their frequency or spatial extent, sleep spindles on average rely on similar neuronal networks, with more variation around this average during individual spindles in case of ex-tremely local IME spindles and less variation in case of the larger, more widespread ECoG spindles. The similar laminar profile of slow and fast spindles is also in line with recent findings showing their similar cou-pling to thalamic downstates (Gonzalez et al., 2018), and these findings overall suggest that any functional differences between sleep spindle subtypes are unlikely to be driven by different generating networks.

Taking advantage of the high sampling rate of IME recordings and the presence of visible single-cell discharges on several IME channels, we performed an analysis of SUA and MUA across cortical layers dur-ing ECoG spindles. Our results replicate the only available human study (Andrillon et al., 2011): we found that SUA was strongest during sleep spindle troughs and the preferred phase of single-unit discharges was during the positive-negative transition of sleep spindles. MUA was also heavily entrained by sleep spindle activity, and exhibited an antiphase correlation pattern with local field potentials across most cortical layers. Unlike in case of LFPg and CSD, we observed notable differences in the preferred SUA phase angles and the LFPg-MUA relationship across spin-dle subtypes (but also cortical layers and patients). This suggests that local neuronal events and their contribution to the local field potential are more heterogeneous across sleep spindle subtypes than laminar to-pography, and may contribute to functional differences. However, due to the lack of consistent effects across the two patients as well as the relatively low statistical power of our study this issue requires further research. It is of note that this study is the first to utilize invasive record-ings together with the IAM. The strong entrainment of MUA/SUA indi-cators of neuronal activity by IAM spindles is evidence for the validity of this detection method.

Our study has a number of shortcomings which may affect the in- terpretability of our results. First, due to the limitations of our method-ology which required data from two different highly invasive electro- physiological sources in humans, together with a synchronizing signal and preferably with post-surgical histological reconstruction, we only had access to a limited amount of data from a small number of pa-tients. The differing age, sex, medication regimen, disease progression and possibly the different pubertal status of our patients may have also affected the results. An additional limitation is that implantation loca-tion of IMEs was biased towards frontal areas, and the results may not generalize to more posterior brain regions. However, even within this heterogeneous population the superficial amplitude maximum of sleep spindles and the similarity of the LFPg and CSD laminar profile of var-ious spindle types was a consistent finding across patients. We found a similar laminar profile and ECoG-IME occurrence and co-occurrence patterns of sleep spindles in the two patients (Patient 1 and Patient 2) with visually identifiable IME spindles, and we successfully replicated previous results about the relationship between spindle amplitude and extent as well as the entrainment of neuronal activity by the spindle oscillation. We emphasize that our results are not incompatible with small differences between the generating mechanisms or laminar pro-files of sleep spindle subtypes, especially if these are limited to spindles occurring in certain cerebral areas. However, our results are not com-patible with very strong previous hypotheses (Timofeev and Chauvette 2013, Piantoni et al 2016) which assumed strikingly different laminar profiles of the different sleep spindle subtypes regardless of electrode placement or other individual characteristics. We had a limited ability to investigate SUA differences due to the limitations of recording and analyzing such activity in humans without experimental manipulations. Possible differences in neuronal activity underlying sleep spindle sub-types warrants further research. Different neuron types, identifiable by their waveform morphology (Peyrache et al., 2012) may also be syn-chronized differently to spindles, with a possible interaction with spin-dle type. Third, our study was performed in epileptic patients under the influence of antiepileptic medications and having just undergone highly invasive surgery. While this is a common limitation of all human inva-sive electrophysiological studies, it may affect the generalizability of the results to healthy participants.

In sum, our study provides evidence that 1) IME spindles occur with a variable topography and co-occur with other IME and ECoG spindles more frequently as a function of greater amplitude and sensor proximity, 2) IME spindles may occur within any cortical layer, but their layer of occurrence is not systematically different across spindle subtypes 3) the laminar profile of various ECoG spindle subtypes is also similar and 4) SUA and MUA dynamics are, however, affected by spindle types.