Spindler: a framework for parametric analysis and detection of spindles in EEG with application to sleep spindles

Spindler takes a single channel of raw EEG data as input. The results reported in this paper use channel Cz. Spindler bandpass filters the signal in the spindle frequency range (11 to 17 Hz for this paper) using EEGLAB's pop_eegfiltnew prior to signal decomposition. Motivated by spindle geometry, Spindler approximates the signal, $s\left(t \right),$ as a weighted sum of Gabor atoms [28]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle s\left(t \right)\approx \underset{n=1}{\overset{N}{\mathop \sum}}\,{{a}_{n}}g\left({{f}_{n}},{{\sigma}_{n}},{{\tau}_{n}},{{\varphi}_{n}} \right).\nonumber \end{align}$$

The Gabor atoms, $g$, are characterized by four parameters:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle g\left(\,f,\sigma ,\tau ,\varphi \right)=~\frac{1}{\sigma \sqrt{2\pi}}{{e}^{-2{{\left(\frac{t-\tau}{2\sigma} \right)}^{2}}}}\cos \left(2\pi f\left(t-\tau \right)+\varphi \right).\nonumber \end{align}$$

Spindler uses an MP dictionary consisting only of atoms with the parameters specified in table 1. N, the total number of atoms used to represent the signal, is N_s times the length of the signal in seconds. f is an allowed atom oscillation frequency in Hz (in the range 11 to 17 Hz in steps of 0.5 Hz for sleep spindles), σ is half of the width of the atom's Gaussian envelope (0.0625, 0.125, or 0.25), and φ is the phase of the sinusoidal oscillation relative to the peak of the Gaussian (assumed to be 0). Temporal shifts, τ, corresponding to placement of atoms centered at every time point are allowed.

Spindler uses an efficient greedy strategy [29, 30] to project candidate atoms on the signal and to successively remove the contribution of the atom whose projection has the largest absolute value. The algorithm stops after selecting a specified maximum number of atoms/second (N_s), not when the representation reaches a specified reconstruction error. This greedy strategy assures that a representation with a larger number of atoms is a superset of representations with smaller numbers of atoms. Thus, the MP representation only has to be computed once to generate the signal on a parameter grid of N_s. An efficient representation merely has to retain the values of (a, f, σ, τ, φ) for each atom added to the representation between the current value of N_s on the parameter grid and the next. The reconstruction of the signal for each selected N_s can be quickly recomputed. By default, Spindler uses thirty-nine evenly spaced N_s values in [0, 0.4]. (The spindle rate diverges for N_s  >  0.3 for all of the examples we have examined.)

It should be noted that a more traditional implementation of MP uses an unrestricted dictionary of Gabor atoms over the full range of frequencies and then applies post-processing to isolate the desired signal. The difficulty with this approach is that when a signal has significant artifactual content, it is difficult to use post-processing to eliminate the effects of these artifacts. By first band-pass filtering the signal and then restricting the representation to only atoms that are likely to represent spindles, Spindler can handle very noisy signals.

To detect the presence of a spindle in a reconstructed signal at (N_s, T_b), Spindler first takes the absolute value of the reconstructed signal and then divides by the 95th percentile of these values to reduce the effect of signal outliers. Spindler sets all scaled values above the 95th percentile to 1. This scaling restricts T_b to the interval [0, 1], allowing a common grid of parameter values to be used for all datasets. Spindler forms a mask marking data points in the reconstructed signal with absolute value that exceeds the specified threshold, T_b. Spindler then applies the minimum spindle separation criteria to the mask to join adjacent pieces and the minimum and maximum length criteria to eliminate spurious pieces. Spindles are detected as regions of consecutive ones in the mask, leading to a list of spindle events for each (N_s, T_b) on the parameter grid.

Given a list of spindle events for each point on the (N_s, T_b) grid, Spindler forms three property surfaces: spindle rate, fraction of time spindling, and spindle length. The geometry of these surfaces, as shown in figure 1 for two representative sleep datasets, is similar across datasets and is used as the basis for parameter selection. N_s is plotted on the horizontal axis and each thin line corresponds to a fixed value of T_b from table 1 ranging from 0 (dark blue) to 1 (dark red) with the extreme values T_b  =  0 and T_b  =  1 emphasized using dashed and dotted teal lines, respectively. The thick solid teal line corresponds to the centered value of the respective parameter, which is computed by averaging the respective parameter values at T_b  =  1 and T_b  =  0 for each value of N_s.

Zoom In Zoom Out Reset image size

Spindler heuristics attempt to select algorithm parameters so that small changes in parameter values do not result in large changes in spindle property values. Observe that for small values of N_s, all of the spindle rate cross-sections overlap (first row of figure 1) indicating that the spindle rate is independent of T_b in this range and changes approximately linearly with N_s. Geometrically this is because the greedy MP representation adds atoms sequentially, so the representation captures only the most energetic (relatively long) spindles regardless of threshold when there are a limited number of atoms to represent the signal. However, as N_s increases, the spindle rate becomes quite sensitive to the threshold.

Spindle rate saturation occurs when all of the spindles that could be detected for a given threshold have been detected and adding more atoms to the representation does not result in more spindles. We observe a downturn and then levelling off of the spindle rate curves (top row of figure 1) as this happens. Spindle rate saturation first occurs for the T_b  =  1 (dotted teal curve) and eventually for other T_b thresholds. We restrict N_s to be less than the first saturation point, which we detect as the lowest value of N_s for which the spindle rate curve for T_b  =  1 has zero slope. This condition implies that adding additional atoms to the representation produces no more spindles. We also require that N_s be large enough that the standard deviation of the spindle rate across thresholds is greater than 0. When it is zero, the spindle rate is completely independent of threshold, indicating that there are not enough atoms to represent even the most prominent spindles. The resulting eligible range of N_s is marked by the thick gray horizontal line on the graphs in figure 1.

To obtain a stable representation with respect to threshold, Spindler calculates the central spindle fraction curve by averaging the curves for T_b  =  0 and T_b  =  1. Spindler then finds $T_{b}^{*}$, the threshold whose spindle fraction curve remains closest to the central fraction curve in the eligible range of N_s. Finally, Spindler selects the value of $N_{s}^{*}$ as the value that minimizes the spindle length for $T_{b}^{*}$ in the eligible range of N_s. As we shall see in the results, no single point in the (N_s, T_b) plane optimizes performance for all methods of matching detected spindles to ground truth. The overall strategy of the Spindler method is not to try for the optimal value in the sense of agreement with expert ratings, but rather to pick a point on the parameter grid that is central so that small changes in parameter values do not result in large changes in property values. The range of eligible points is specifically selected to 'stay away from the edges'.

The Count method [19] tabulates the number of spindles detected in segments of a specified length (30 s for this paper). A spindle that overlaps multiple segments only counts in the first interval of overlap. The number of true positives (TP) in each segment is the minimum of the number of detected spindles and the number of true spindles in that interval. If there are more true spindles than detected spindles, the difference is the number of false negatives (FN) in the segment. Conversely, if there are more detected spindles than true spindles in the segment, the difference is the number of false positives (FP) in the segment. TN is calculated using the default method. The overall confusion matrix is computed as the sum of the confusion matrices computed in each segment.

The Hit method detects a TP if there is any overlap between a true spindle and any detected spindle. Hit detects a FP if a given detected spindle does not overlap any true spindle and a FN if a true spindle does not overlap any detected spindle. Spindler finds the complementary intervals (start and ends) for the true spindles and the detected spindles and counts overlaps in the complementary intervals as true negatives. However, we caution users against using performance measures based on TN, particularly for the Hit method. Hit matching tends to estimate performance more optimistically than other methods for non-TN metrics and less optimistically for TN metrics. One detected spindle can overlap multiple true spindles without penalty (and vice versa). Further, a long period with no spindles counts as a single TN, rather than weighted by the interval length. Another difficulty with the Hit method is that labelling the entire dataset as a single spindle (TP) gives precision and recall values of 1. Hit performance is relatively insensitive to parameter selection but should be used with caution if there are too many very long events.

The Hit method has some similarities to the method used by Nonclercq et al [12]. These authors decomposed the signal into 0.5 s overlapping windows spaced 125 ms apart and classified each window as a spindle or a non-spindle based on frequency-amplitude properties of the windows. They categorized windows based on whether a particular window had matching labels for true or detected spindles. The Nonclercq algorithm categorizes windows with no label as true negatives. Because the Nonclercq algorithm relies on knowing the distribution of frequency-amplitude values for true positives, it is not suited for a general implementation.

The Intersect method uses the union/intersection rule similar to the algorithm proposed by Warby et al [23] for sleep spindle evaluation. The Intersect method calculates the ratio between the intersection and the union of the time intervals representing two compared events (true and detected). If this ratio is greater than the intersect tolerance (0.2 s by default), the corresponding spindles are potential matches. The Intersect method chooses the event pair with the greatest eligible ratio as a match, categorizes the match as true positive, and removes the event pair from further consideration. Intersect categorizes unmatched true events as false negatives and unmatched detected events as false positives. Warby et al do not provide pseudo code for calculating true negatives. Our implementation uses the default method for calculating true negatives for the Intersect method.

The Onset method determines whether the starting times of true and detected spindles agree to within a specified onset tolerance. A true positive occurs when a detected spindle starts within onset tolerance of the start of a true spindle. A false positive occurs when a detected spindle does not start within onset tolerance of any true spindle. A false negative occurs when a true spindle does not have a detected spindle within onset tolerance of it. Our implementation of Onset uses the default method for true negatives. As with the Intersect method, Onset removes any matched events from further consideration. This method is frequently used for the assessment of sleep spindles with an onset tolerance of 0.5 s [14, 24]. We have found that the performance results are extremely sensitive to the onset tolerance and that 0.5 s gives results that appear to be far too optimistic. Spindler uses 0.3 s as the onset tolerance by default.

The Time method, introduced by Lawhern et al [31] to evaluate alpha spindles, uses the actual length of time spent in each state rather than treating matches as events. To account for the difficulty in determining the exact starting time of a spindle, Time allows the extension of the events in either direction by a timing tolerance to get a better match. By default, Time uses 0.2 s as the timing tolerance. Time categorizes the periods during which a detected spindle does not overlap any true spindle as false positive time, and periods during which a true spindle does not overlap any detected spindle as false negative time. Any periods not labelled as TN, FP, or FN accumulate as Null agreement or true negative times. The Time method does not distinguish between single overlaps and multiple overlaps, but rather calculates true positives based on the amount of time that true spindles and detected spindles overlap.

The continuous wavelet probabilistic estimate methods proposed by Tsanas and Clifford [14] decompose the signal using a continuous Morlet wavelet transform based on 131 Morlet scales and pseudo-frequencies in the range 5.4 to 40.5 Hz. Twenty-four (24) of the basis wavelets fall in the sleep spindle frequency range of 11–17 Hz. At each time sample, the algorithm estimates the probability of a spindle based on how many of the top ten wavelet coefficients at that instant have Morlet scales in the spindle frequency range. The algorithm smooths the instantaneous probability estimate using 10 nearest neighbors and combines adjacent points into spindles by applying several rules governing spindle probabilities and durations. An example of a rule is: if the two intervals have average spindle probability  >0.7, are each at least 0.1 s long and less than 0.3 s apart, the intervals are merged into a single spindle. Tsanas and Clifford have provided two versions of their algorithm (denoted by CWT-a7 and CWT-a8) using different methods of smoothing the instantaneous probability function. The implementations take the frequency range as a parameter, but hardcode all thresholds, spindle duration criteria and spindle separation criteria. These thresholds and the Morlet wavelet parameters heavily rely on a signal-sampling rate of 100 Hz. The Spindler toolbox wrapper for the implementation resamples the signal to 100 Hz prior to execution. On return from the call, the Spindler toolbox wrapper combines adjacent events, matches to ground truth, and calculates performance measures.

Wendt et al [10] proposed a sleep spindle detection algorithm that takes input from two EEG channels: a central channel such as C3 and an occipital channel such as O1. The algorithm detects candidate sleep spindles using local extrema of the signal envelop from the central channel after filtering in the 11 to 16 Hz frequency band. The algorithm further corrects for 'alpha effects' by removing candidate spindles if the spindle frequency is  ⩽13 Hz and the power of the band-passed central channel in the interval defining the candidate spindle is less than the power in the occipital channel. The algorithm assumes that the input signal corresponds to NREM sleep. We use the implementation provided by Warby et al [23]. The algorithm intrinsically relies on specific frequency filters and cannot be easily adapted for detection of other types of spindles. The Spindler toolbox provides a wrapper to the author's implementation of this algorithm for combining adjacent events, matching to ground truth, and calculating performance measures.

Spinky [18, 19] provides a supervised method of sleep spindle detection that decomposes a signal using the tunable Q-factor wavelet transform. Q is a hard-coded constant chosen so that the wavelets have a minimum number of cycles corresponding to the number of oscillations in the desired frequency band 11 to 17 Hz for an interval of 0.5 s. After wavelet decomposition, Spinky determines a threshold as the maximum of Sensitivity–FDR using the Count matching method on 30 s intervals. Spinky calculates minimum and maximum allowed values of this threshold based on the signal. Spindler uses a grid of values (by default spaced 10 units apart) between these maximum and minimum values for parametric analysis. The Spindler toolbox provides a wrapper for the original Spinky MATLAB implementation of spindle detection. Note that the original Spinky toolbox also provides K-complex detection and a python user interface to facilitate manual ratings, which are not included in the Spindler toolbox.

McSleep [21] and its single-channel predecessor DETOKS [20] are non-linear methods that model EEG as the sum of a sparse transient, a low frequency signal, and a sparse oscillatory signal. The methods require four regularization parameters (λ₀, λ₁, λ₂, μ), a threshold (c), and some block size and iteration size parameters. The authors indicate that λ₀ and λ₁ can be fixed based on the signal sampling frequency. They optimized the remaining parameters using the F₁ metric with Count matching. The McSleep toolbox does not provide code for parameter optimization. Our analysis uses the parameters suggested in [21] for λ₀, λ₁, and μ for the Dreams and the MASS-SS2 datasets, respectively. Spindler uses a grid for λ₂ and c to generate parameter curves for McSleep. By default, the grid for λ₂ is [20, 50] with spacing of 1, and the grid for c is [0.5, 3] with a spacing of 0.1. This paper uses 3 channels (central frontal, central, central occipital) for the analysis.

We used the 30 min excerpts from the publicly available Dreams Sleep Spindles Database [24]. The second expert's ratings only included the spindle starting times. Following Devuyst et al, we assumed a spindle length of 1 s for the second expert. The excerpts did not account for a particular sleep stage and were not pre-processed. In each case, experts annotated the central channel corresponding to position CZ.

The Montreal Sleep Archive [22, 25] Cohort 1 subset SS2 consists of whole night recordings of 19 healthy subjects sampled at 256 Hz. We used all 19 recordings, 15 of which were rated by two experts and 4 of which were only rated by Expert 1 (subjects 4, 8, 15, and 16). The whole night recordings had also had ratings of sleep stage. We extracted the longest segment of stage 2 sleep in each recording to use in the analysis.