Two publicly available databases were used to evaluate 6 simplex and 3 hybrid automated sleep spindle detection algorithms. The first sleep EEG database was from the DREAMS Spindles Database of University of MONS—TCTS Laboratory and Université Libre de Bruxelles—CHU de Charleroi Sleep Laboratory (Devuyst, 2013). It consists of eight patients (4 men and 4 women aged between 31 and 53) with different pathologies (Devuyst et al., 2011). A segment of 30 min of the central EEG channel (C3-A1 or Cz-A1) from these 8 individual patients are publicly available on the DREAMS Database website where the sampling frequency is 200 Hz (6 patients), 100 Hz (1 patient), and 50 Hz (1 patient) respectively. Sleep spindles of these 30-min-long EEG signals are annotated independently by two experts. The second expert annotated only six out of eight datasets and the spindles are uniformly assigned 1-s-long duration. In order to build a common ground truth from multiple raters, there are several different approaches such as treating individual rater separately (O'Reilly and Nielsen, 2015), the total agreement (intersection operation; Devuyst et al., 2011), and partial agreement (union operation; Warby et al., 2014; Tsanas and Clifford, 2015). For forming (total/partial) agreement sets, we took the same approach proposed by Tsanas and Clifford (2015) specifically for the DREAMS database. In the case of overlapped spindles, we only kept the annotations by expert 1 for their better duration assessment. In total, the number of identified spindles was 355 (from 8 subjects) by scorer 1, 409 (from 6 subjects) by scorer 2, 160 by both, and 604 by either, respectively. The details are summarized in Table 1. Note that the authors of the DREAMS database did not specify which scoring rules the experts used for scoring spindles.

The second EEG sleep database was a subset (denoted as SS2) of the MASS database. The MASS database consists of 200 polysomnographic (PSG) recordings gathered from eight research protocols conducted between 2001 and 2013 in three different laboratories at the Center for Advanced Research in Sleep Medicine (CARSM), Montreal, Canada (O'Reilly et al., 2014). The MASS SS2 consists of 19 complete-night PSG recordings sampled at 256 Hz from young healthy subjects. For this subset, sleep spindles are scored by two experts on NREM stage 2 sleep epochs and on channel C3 with linked-ear reference. The first scorer used traditional AASM scoring rules. The second scorer used both broad-band EEG signals (0.35–35 Hz band) and sigma filtered signals (11–17 Hz band) to facilitate the identification of short duration, small amplitude or obscured (by delta waves or K-complexes) spindles. Also, no minimal spindle duration was used by the second scorer and four datasets (out of the 19) were not scored due to poor quality sleep or signal (O'Reilly and Nielsen, 2015). In total, the number of identified spindles was 11,147 (from 19 subjects) by scorer 1, 21,951 (from 15 subjects) by scorer 2, 9,514 by both, and 22,388 by either, respectively. The details are summarized in Table 2. Note that the total/partial (intersection/union) agreement in Table 2 is solely given as supplementary information to describe inter-rater agreement. The intersection and union operations were carried out based on each sample point (by-sample) and then the resultant spindle points were regrouped into individual sleep spindles (by-event) without duration checking.

The performance of a diagnostic test is generally characterized by sensitivity and specificity. However, sleep spindles are sparse events such that their lengths sum up to 8.2 (mean) ± 4.9 (standard deviation) and 29.4 ± 11.2 min per whole night sleep according to scorers 1's and 2's annotations of the MASS SS2 database respectively. The performance evaluation using specificity measurement will be high and therefore provide unrealistically positive results. To avoid this pitfall, Warby et al. (2014) propose using precision, recall, and F₁-score for the evaluation of infrequent, discrete events such as sleep spindles in the EEG signal. Let TP denote the amount of true positives, FN false negatives and FP false positives, precision, P, and recall, R, are defined as:

Taking the weighted harmonic average of precision and recall leads to the F-score,

If we assume a uniform prior (β = 1), then

Although P, R, and F₁ are commonly used evaluation metrics for the assessment of spindle detector performance, so far their potentials have not been fully explored yet. In the original work of Warby et al. (2014), spindle detectors were only tested with default operating parameters. This view is too narrow as a detector generally has several adjustable operating parameters that allow the user to optimize its performance between two competing objectives: minimizing FN and minimizing FP. For example, we can take the root mean square (RMS) algorithm of Martin et al. (2013) and test it out on the DREAMS database, the default operating parameter with a threshold of 0.95 (95 percentile) and a duration between 0.3 and 3 s. The result is simply a point (red circle) among many other possible performances (black crosses) on the objective space as shown in Figures 1A,B. O'Reilly and Nielsen (2015) took a step forward by making the threshold an adjustable parameter within a range between 0.7 and 0.995. The resultant FN-FP pairs, PR-curve and F₁-curve (blue dots in Figures 1B–D) provide a broader view to evaluate the detector's performance. In this study, we extended the idea of O'Reilly and Nielsen (2015) even further by making adjustable more operating parameters (such as the upper and lower duration criteria) and derived the performance metrics from the Pareto optimal solutions (green dots in Figure 1) in the multi-objective optimization context. Figure 1 demonstrates that results from a detector can be improved substantially in a multi-objective context by allowing the operating parameters to adapt to a training gold standard and that its Pareto optimal solutions enable us to define some useful performance metrics uniquely. The following section summarizes the proposed performance metrics more formally.

In a multi-objective context, the performance of a sleep spindle detector is described by a pair of raw numbers: FN and FP. The number of detections on the C3 channel of EEG signals that do not have an appropriate overlap rate (R_ov) with any true sleep spindle (SS) is defined as FP; similarly, the number of true sleep spindles that are not detected by automatic detectors is defined as FN. Note that FN and FP were defined in the event-by-event analysis context (Warby et al., 2014) throughout this study. Note also that a true positive event-detection (D) was scored based on an overlap rate,

We followed these scoring rules to make our results comparable to the work by Warby et al. (2014) because the first three Fourier-based detectors evaluated in this study were derived directly from their work.

In mathematical terms, let FN and FP be described as functions of an operating parameter set x: FP(x) and FN(x). For example, x includes threshold, lower duration, and upper duration for the RMS algorithm illustrated in Figure 1. Figure 1A shows a random selection x's of 100 possible operating parameter combinations with their results [FN(x), FP(x)] scattered around the objective space. The optimal solutions of such a detector involve minimizing both FN and FP rates simultaneously. Since low FN and FP rates are two conflicting objectives, this problem does not produce a single optimal solution but a set of possible solutions known as a Pareto optimal set, which results in a Pareto front (green dots in Figure 1A) in objective space. A Pareto front is essentially an objective boundary such that any solution on the front can only be outperformed by another solution in at most one of the two competing objectives. Therefore, a Pareto optimal set is also called a Pareto non-dominated set.

Formally the multi-objective optimization problem can be equivalently stated as minimizing a two-objective vector

where x is the vector of a detector's operating parameters. A solution x₁ is said to dominate x₂ if and only if

where x₂ ≠ x₁. From the aforementioned Pareto optimal solutions and Equations (1–4), we are able to derive PR-curves and F₁-curves that are Pareto optimal as illustrated in Figures 1C,D.

The solutions of Equation (6) satisfying the conditions listed in Equation (7) can be found rather efficiently by MOEA with genetic mechanisms. In a standard genetic algorithm, there are usually four steps in the evolutionary procedure: (1) randomly initializing the solution population, (2) evaluating and assigning a fitness value for each individual in the population according to its performance, (3) selecting individuals based on their fitness values to procreate, and (4) using crossover and mutation to produce next generation from the selected individuals. In this study, we used the SPEA2 algorithm to find the Pareto fronts of all of the 9 examined detectors. SPEA2 was selected for its fast convergence rate and good performance because it kept a relatively small, yet diverse population (Zitzler et al., 2001a).

A brief work flowchart of our proposed metrics is provided in Figure 2. It consists of three key modules: (1) detector, (2) subsample, and (3) SPEA2 modules. Six simplex and 3 hybrid detectors with their operating parameter descriptions are given in Sections Six Simplex Detectors and Three Hybridization Detectors respectively. Subsample and SPEA2 implementations are delineated in Sections Subsample Strategy and SPEA2 Module.

Fourier-based bandpass filtering (Figure 3, 2nd row) is the foundation of many automated detection methods for identifying the frequency of sleep spindles. The main difference among such bandpass automatic spindle detectors is to apply various methodological strategies for improving the identification of the “right amplitude” and “right duration” of sleep spindles. Alternatively, CWT and HHT for frequency discrimination are also proposed in more recent approaches. We have focused on six existing algorithms: three using Fourier filtering, one using CWT, and two using HHT. They are briefly reviewed as follows.

There are different ways in hybridizing different detectors for performance improvement. Simple approaches include intersection (total agreement) and union (partial agreement) of the results from multiple detectors. However, such a simple design through giving different weights to different detectors increase the computational cost cumulatively because it needs to collect the results from all participating detectors. Instead, in this study, we followed a double reading paradigm, in which the first detector screened all data to find potential sleep spindles as usual. However, the second detector only reviewed the sleep spindle candidates identified by the first detector. This approach was taken to balance between accuracy and computation cost.

We selected d₂, d₄, and d₆ from the RMS, CWT, and HHT categories respectively for composing 3 hybridization detectors. First, hybrid detector d₇ was composed of d₆ as the first detector and d₂ as the second detector. Its operating vector x_d7 = {x_d6, x_d2} had 17 parameters in total. Second, hybrid detector d₈ was composed of d₂ as the first detector and d₆ as the second detector. Its operating vector x_d8 = {x_d2, x_d6} had 17 parameters. Third, hybrid detector d₉ was composed of d₄ as the first detector and d₆ as the second detector. Its operating vector x_d9 = {x_d4, x_d6} had 27 parameters.

Since the total length of sleep spindles in the MASS SS2 datasets ranged from only 8 to 30 min per night, we decided to extract a subsample of 60 min from each whole night sleep in the training datasets to reduce SPEA2 computation time. Considering that our tested detectors used intrinsically different methods to estimate frequency, we did not want to create a common, single subsample for all of them. Instead, we tried to apply the optimal parameters derived from the DREAMS database by each individual detector to collect its own false positives from the training set of MASS database. These detector-dependent false positives and all true positives (defined by the chosen gold standard) were used to create small segments of signals. These small segments were dilated by 2.5 s from both ends and merged (if overlapped) to form a pool of bigger segments. From this pool, we randomly picked a total amount of 60 min long signals with true and false spindles. Tested with the maximal F_0.5-, F₁-, and F₂-score (beta = 0.5, 1, 2 in Equation 3), we found that the parameters achieving the highest F₂-score based on the union gold standard of the DREAMS database were a good choice for generating false detection subsamples from the MASS database.

SPEA2 software has three mandatory input operating parameters: (1) solution population size, (2) maximum number of generations, and (3) minimum and maximum boundary values of the optimized detector's operating parameters. In this study, the solution population sizes and maximum generations were empirically chosen such that they were at least 10 times of the number of operating parameters. On several occasions, some of these numbers were increased by 50–100 than the aforementioned recommendations to ensure convergence to Pareto front solutions. Note that population size, number of generations, and parameter boundary values were empirically chosen in the training stage of the hold-out experiment described in the following section. Their values are listed in Section Spindle Detection Performance on DREAMS Database. Finally, for the other SPEA2 parameters (such as mutation and crossover) we used the default values implemented by Popov (2005).

For reporting the performance of each automatic spindle detector, we first derive the vector of operating parameters that achieves the maximal F₁-score from its Pareto-optimal solutions solved by the SPEA2 algorithm based on a training dataset. Second, we apply the optimal operating parameter vector to a test dataset (unseen to the detector) to assess the resultant F₁-score (Figure 2). Two validation strategies were conducted to assess the performance of spindle detectors in the proposed Pareto optimal and gold standard adaptive context. First, a hold-out strategy was used to assess overfitting by our SPEA2 implementations. Second, a k-fold cross-validation strategy was used to assess spindle detectors.

In our hold-out experiment, we divided the datasets into three groups: (1) all DREAMS datasets, (2) the first half of MASS datasets (subjects 1–9), and (3) the second half of MASS datasets (subjects 10–19). Subjects 7 and 8 were excluded from the DREAMS database if the gold standard of scorer 2 or intersection was applied; data 4, 8, 15, 16 were excluded from the MASS database if scorer 2's gold standard was applied. There were two versions of EEG signals in the MASS database. We used the new version published in 2015 to evaluate all the 9 detectors described in Sections Six Simplex Detectors and Three Hybridization Detectors. Our hold-out paradigm was conducted such that DREAMS and the first half of MASS were freely explored to find the proper ranges of adjusted parameters for each detector and to choose the proper population size and the number of generations for SPEA2. At this stage, we also experimented with a few different ways to generate an appropriate subsample for MASS datasets. Although subjecting a detector to 9 whole-night sleep EEG was not an impossibility, the idea of “gold standard adaptive” would work for most clinicians only if the training time was within hours by a regular personal computer. Once the parameter ranges, population size, number of generations, subsample strategy were decided, they were then fixed and used to train all the detectors with every different gold standard to find their own Pareto optimal solutions. The second half of MASS datasets that remained unseen to the spindle detectors during the entire training stage were then used to test the final solutions only once for deciding a conservative error bound (Brun et al., 2008).

The 95% confidence intervals (CI) of the estimated F₁-scores by the training and test datasets in the hold-out experiment were derived by using probabilistic interpretation. Since the distributions of the precision and recall are Beta distributions (Goutte and Gaussier, 2005), the 95% CI of a F₁-score can be estimated by running Monte Carlo simulation of precision and recall that are calculated from the data. We assumed that the number of true negatives in the data was 10 times the number of true positives and estimate the 95% CI by running Monte Carlo simulation for 10,000 times.

Lastly, 3-fold cross-validation, where the subjects in the MASS SS2 database were randomly partitioned into 3 equal sized subgroups, was conducted to assess the performance of 9 detectors. Of the 3 subgroups, a single subgroup was retained as the validation data for testing the detectors, and the remaining 2 subgroups were used as training data. In each fold, multiple Pareto optimal parameter vectors were solved by using a subsample of the training data as shown in Figure 2 for computation efficiency. However, the parameter vector x which generated the maximum F₁-score was judged based on the complete-night training data. This maximum F₁-score x was then used to assess the FN and FP results on the complete-night validation data. The cross-validation process was repeated 3 times with each of the 3 subgroups used exactly once as the validation data. The validation results from all 3-folds were finally grouped together to produce a single estimation.

All the sleep spindle detection algorithms evaluated in this study were implemented in Matlab version R2015a (MathWorks, Natick, MA, USA) and C language. Fourier-based bandpass filtering software (d₁, d₂, and d₃) was originally developed by Warby et al. (2014). The CWT-based software (d₄) was developed by Tsanas and Clifford (2015). Parts of these Matlab codes were rewritten in C language for reducing SPEA2 computation time. The HHT-based software (d₅ and d₆) and hybrid software (d₇, d₈, and d₉) was developed by our group. All Matlab and C source codes are available in Supplementary Material. Third party software such as SPEA2 and EMD is available in Popov (2005) and Wang et al. (2014) respectively. Computation time estimates were performed on a Fujitsu Lifebook laptop with Intel Core i7-3632QM processors at 2.20 GHz, using 12 GB of RAM memory running a 64-bit Windows 10 operating system.

In the training stage of the hold-out experiment, we found that SPEA2 was able to converge to an adequate set of solutions with a population size of 100 parameter vectors that evolved 100 generations for detectors d₁-d₃ with 3–5 adjustable parameters. For detectors d₄-d₉ with more than 10 parameters, the population size was ascertained by 10 times the number of parameters and the number of generations was empirically determined by the population size plus 50–100 to ensure convergence to Pareto front solutions for the DREAMS database. In summary, the population sizes of d₁-d₉ were 100, 100, 100, 150, 120, 120, 170, 170, and 270 respectively. The maximum generations were 100, 100, 100, 250, 200, 200, 250, 250, and 320 respectively. The ranges of lower and upper parameter values were also empirically determined based the Pareto front results of the DREAMS database for all 9 detectors. The ascertained minimum and maximum (boundary) values of detectors d₁-d₅ parameters were

The minimum and maximum values of x_d6 were the same as x_d5. The hybrid detectors d₇-d₉ used the same minimum and maximum boundaries for their respective components d₂, d₄ and d₆. The ranges of other parameters are given in Supplementary Material for code transparency and result reproducibility. Figure 4 shows the Pareto front-derived PR-curves and F₁-scores of 9 sleep spindle detectors, simplex d₁–d₆ and hybrid d₇–d₉, solved by using the SPEA2 algorithm with the DREAMS database of 4 different gold standards listed in Table 1. The PR-curves (Figures 4A,C,E,G) and F₁-scores (Figures 4B,D,F,H) vary dramatically from the worst performance evaluated based on scorer 1's gold standard to the gold standard of intersection, union, and finally to the best performance of scorer 2's gold standard. Although the F₁-scores of all tested detectors vary, the hybrid detectors (especially d₉) seem to perform better regardless of which gold standard is used.

After we trained these 9 detectors with the subsamples extracted from the first half of MASS database by using the subsample strategy (Section Subsample Strategy), we reported their training performance based on the FN and FP calculated from the whole night training datasets with the Pareto optimal solutions derived from subsamples. The Pareto front-derived PR-curves estimated in the training procedure are shown as red dots in Figures 5, 6 for the gold standard of scorers 1 and 2 respectively. The PR-values of test results derived from the second half of MASS database are shown as blue dots (in Figures 5, 6) and are connected by a light blue line to the training results (red dots) that are estimated by using the same vector of operating parameters. The maximal F₁-scores acquired by training are marked by red squares in Figures 5, 6. Note that the blue squares in Figures 5, 6 are not the maximal test F₁-scores but the test results by using the same set of operating parameters of the corresponding red squares. Here we did not report the maximal test F₁-scores because in real practice we can choose only a small number of operating parameters after the software system has been optimized. The PR-values, F₁-scores, and 95% CI of F₁ at the operating points marked by squares in Figures 5, 6 are listed in Table 3 and their operating parameter values are given in Supplementary Material.

The F₁-scores for all 9 detectors by the gold standards of scorer 1 and scorer 2 are illustrated in Figure 7. From Figures 7A,B, Table 3, we observe that Fourier-based simplex detectors d₁–d₃ perform the worst with the estimated maximal F₁-scores 0.563–0.609 by scorer 1's gold standard. HHT-based simplex detectors perform better with the estimated maximal F₁-scores 0.622–0.699. Hybrid detectors d₇ and d₈ perform the best with the estimated maximal F₁-scores 0.722–0.738. CWT-based simplex detector d₄ and hybrid d₉ perform well with high training F₁-scores at 0.631 and 0.735 but low test results of 0.589 and 0.673 respectively. The results in Table 3 indicate that Fourier-based detectors improve significantly by scorer 2's gold standard [scored based on both broad-band EEG signals (0.35–35 Hz band) and sigma-band filtered signals (11–17 Hz band)]. Despite the existence of a biased gold standard in favor of Fourier filtering, hybrid d₇ and d₈ still outperform with maximal F₁-scores 0.722–0.744 compared to 0.637–0.7 by simplex Fourier-based detectors. Finally, automatic detectors are compared to the iso-curves of F₁-scores of 0.75 (good) and 0.67 (average) in Figure 8. The performance levels of “good” and “average” refer to the study by Warby et al. (2014), where experts' average F₁ was 0.75 and non-experts' consensus 0.67. Figures 8A,B illustrate that the optimized simplex detectors d₁-d₅ are indeed inferior to human performance of 0.67. However, HHT-based simplex detector d₆ perform comparably to the level of non-experts. Hybrid detectors d₇ and d₈ perform well on a level slightly below that of experts regardless which gold standard is used (Figures 8A–D).

The maximum F₁-scores of d₁–d₉ estimated by 3-fold cross-validation were 0.588, 0.604, 0.583, 0.61, 0.633, 0.679, 0.732, 0.727, and 0.7 respectively based on scorer 1's gold standard. They were 0.622, 0.677, 0.651, 0.651, 0.644, 0.669, 0.737, 0.726, and 0.68 respectively based on scorer 2's gold standard. Most of these numbers are substantially better than the baseline performance by using their originally published parameters. Both 3-fold validation, baseline performance, and their corresponding PR-values are also listed in Table 3 for comparison. These cross-validation estimates were comparable to the hold-out estimates except d₅'s F₁-score on scorer 2's gold standard. The hold-out test F₁-score was 0.608 (95% confidence intervals 0.601–0.616), which was significantly smaller than the cross-validation's F₁ estimate of 0.644.

The SPEA2 training time for d₁–d₃, d₄, and d₅–d₆ was 27–42 m, 5.3 h, and 1.7–2.2 h, respectively; for hybrid detectors d₇, d₈, and d₉ was 4.8, 4.4, and 16 h, respectively. Although the computation time (in seconds) to preprocess the C3-A1 channel of a complete-night MASS dataset was 0.9, 7, 16, 28, 73, 835, 838, 390, 390 for d₁–d₉ respectively, the execution time of spindle detectors d₁–d₉ implemented in C language was 0.06, 0.09, 0.16, 0.63, 0.20, 0.14, 0.22, 0.27, and 0.71, which is well below 1 s per complete-night dataset. Therefore, the training time was mainly proportional to the product of population size and number of generations, which approximated the square of a detector's parameter number.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.