SGM: a novel time-frequency algorithm based on unsupervised learning improves high-frequency oscillation detection in epilepsy

Objective. We propose a novel automated method called the S-Transform Gaussian Mixture detection algorithm (SGM) to detect high-frequency oscillations (HFO) combining the strengths of different families of previously published detectors. Approach. This algorithm does not depend on parameter tuning on a subject (or database) basis, uses time-frequency characteristics, and relies on non-supervised classification to determine if the events standing out from the baseline activity are HFO or not. SGM consists of three steps: the first stage computes the signal baseline using the entropy of the autocorrelation; the second uses the S-Transform to obtain several time-frequency features (area, entropy, and time and frequency widths); and in the third stage Gaussian mixture models cluster time-frequency features to decide if events correspond to HFO-like activity. To validate the SGM algorithm we tested its performance in simulated and real environments. Main results. We assessed the algorithm on a publicly available simulated stereoelectroencephalographic (SEEG) database with varying signal-to-noise ratios (SNR), obtaining very good results for medium and high SNR signals. We further tested the SGM algorithm on real signals from patients with focal epilepsy, in which HFO detection was performed visually by experts, yielding a high agreement between experts and SGM. Significance. The SGM algorithm displayed proper performance in simulated and real environments and therefore can be used for non-supervised detection of HFO. This non-supervised algorithm does not require previous labelling by experts or parameter adjustment depending on the subject or database considered. SGM is not a computationally intensive algorithm, making it suitable to detect and characterize HFO in long-term SEEG recordings.


Introduction
Epilepsy is a neurological disorder characterized by abnormal synchronous electrical activity that can cause recurrent seizures. The first-line therapy for epilepsy consists of antiepileptic drugs, but up to 30% of patients remain refractory to pharmacological treatment [1]. Some of these patients having focal seizures may undergo surgery for resection or disconnection of the epileptogenic area, that is, the minimal area of the cortex that must be resected or disconnected to achieve freedom of seizures. Data from interictal and ictal EEG, as well as from structural and functional imaging, are used to infer the location and extent of the epileptogenic area. During interictal periods, the peculiar activity of the epileptic brain can be measured using invasive and noninvasive electrophysiological signal recordings. Recent studies analyzed high-frequency oscillations (HFO) as promising biomarkers of epileptic tissue [2,3]. HFO can be measured in invasive stereoelectroencephalography (SEEG) and recent studies have assessed its detectability in non-invasive electroencephalography (EEG) and magnetoencephalography (MEG) [4][5][6][7][8][9][10].
HFO are described as spontaneous events in the frequency range between 80 and 1000 Hz that stand out from baseline activity and last at least four cycles. HFO are commonly classified into ripples (80-250 Hz) and fast ripples (higher than 250 Hz) [11]. As HFO also appear during interictal periods there is no need to record ictal activity, thereby reducing discomfort and risks for the patient [4]. Some studies have reported that the removal of tissue associated with pathological interictal HFO is a predictor for good surgical outcome, better than tissue related to the seizure onset zone, which can only be determined using ictal EEG [3,9,12,13]. Furthermore, there is increasing interest in evaluating if changes in the appearance rate, morphology, frequency, or other HFO characteristics could be associated with seizure generation. Some of these effects have already been observed in animal models [14]. The relationship between HFO and epilepsy is still an open field of study and therefore a reliable identification of these electrophysiological biomarkers is of key importance.
The detection of HFO is either carried out by visual detection performed by experts or using automatic detectors. As HFO are short-time and lowamplitude events, visual detection is very timeconsuming: it can take up to 10 h to label HFO on a 10-minute and 10-channel recording [11]. In addition, visual detection may be influenced by potential reviewer bias [15,16]. Recently, the improvement of recording techniques has allowed the acquisition of long-term EEG, and current SEEG recordings comprise weeks of multi-channel signals, usually more than 100 channels, rendering visual detection a very difficult and tedious task. A recent prospective multicentric clinical trial evaluated the use of interictal HFO in epilepsy surgery for prediction of postsurgical seizure outcome [15]. This study concluded that given the current alternatives and resources for HFO detection, the removal of the area producing HFO was not associated with successful surgery outcome on an individual patient basis. One of the major issues influencing these negative results, as stated by the authors, was that the methodology for the detection of HFO needed improvement.
Consequently, being manual labelling of HFO not feasible in large-scale studies, the development of reliable automatic procedures is crucial for practical clinical use. In addition, electrophysiological recording may produce some channels with continuous high frequency activity [17] that may hinder detector performance. During the last 15 years, several detectors have been proposed and improved [18][19][20][21][22] but none of them has achieved a large consensus in the scientific community. One of the main drawbacks of some detectors is that they rely on fixed or variable thresholds, obtained using standard deviations or percentiles on the same data under study, to detect events that stand out from baseline activity [18,20,21]. Other more sophisticated detectors use timefrequency features to decide if an event is an HFO, achieving better results. However, they use fixed parameters, which may not work in different scenarios, to tell HFO apart from non-HFO [21,23]. In contrast, several machine-learning techniques have been proposed in the last years to detect HFO, but most of them rely on supervised methods [8,[24][25][26][27], and therefore need previous labelling of the data to perform model training.
In this study, we propose a novel automated method for detecting HFO in SEEG that combines the strengths of different families of previously published detectors. This algorithm, named S-Transform Gaussian Mixture detection (SGM), does not require fine tuning of the parameters used for detection. After a selection of events of interest (EOI), time-frequency characteristics are obtained and non-supervised classification determines if the events that stand out from baseline activity are HFO or not. The algorithm was programmed using open-source software, Python and MNE package [28].

Simulated and real signals
Freely available simulated data, validated an explained in detail in [23], were created by using real recordings of non-REM slow wave sleep of drug resistant epileptic patients undergoing pre-surgical examination with SEEG. The sampling rate was 2048 Hz with an antialiasing filter set at one third of the sampling frequency (688 Hz). Using a graphical user interface, experts selected spikes, ripples, fast-ripples, and background activity visually, using a bipolar montage. Simulated signals were generated by integrating randomly these events into a background signal. The background activity was obtained using several pieces of baseline selected by experts. These pieces were used for estimating the coefficients of an autoregressive model. Background signals were then generated by filtering Gaussian white noise using the averaged autoregressive coefficients. Events were randomly inserted (3 events min −1 ) to the background activity so that two events can appear simultaneously (e.g. a spike co-occurring with a ripple; or a ripple co-occurring with a fast ripple). HFO were individually inserted into the background signals scaled to fit the selected signal-to-noise ratio (SNR) in their relative frequency band. Simulated signals consisted of 30 channels for each value of SNR (0, 5, 10 and 15 dB).
Even though simulated signals provide a controlled testing environment, it is important to assess automatic algorithms in real conditions, as signals may have different artifacts not considered in simulations, and channels without interictal activity because they fall outside the seizure onset zone. In this study, non-rem sleep recordings of eleven patients (one woman) diagnosed with intractable epilepsy undergoing presurgical examination with intracranial electrodes were selected. Patient ages ranged from 8.3 to 51.3 years (mean 27.8, median 25.5) at the time of the study with refractory epilepsy because of focal cortical dysplasia (5 FCD type 1, 6 type 2) in frontal (4 patients), temporal (4 patients), temporo-occipital (2 patients) or parietal (1 patient) lobe. They were on 1 (2 patients), 2 (6 patients) or 3 (3 patients) antiepileptic drugs, the more frequent being levetiracetam (6 patients), carbamazepine (4 patients) and oxcarbazepine (3 patients). In all cases antiepileptic drugs were transiently stopped or significantly reduced during the study in order to register seizures.
All of them presented interictal activity coming from only one focal generator. SEEG signals were acquired at the University Hospital Freiburg using a sampling rate of 1024 Hz and 35 channels/subject. The general recommendation when inspecting signals for presence of HFO is to use a sampling frequency at least four times higher than the frequency of the oscillation of interest [29] and, for this reason, the maximum HFO frequency was set to 250 Hz. These oscillations were detected visually by two different experts. For each channel, the first three minutes of data were independently analyzed, simultaneously inspecting the original and the band-pass filtered signal (80 to 250 Hz), as well as the timefrequency distribution of the events considered relevant. Well-defined HFO were identified following the definition of basic HFO characteristics [11]. Only events selected by the two experts were selected as HFO. The number of HFOs visually detected for each subject ranged between 13 171 and 7252.

Unsupervised algorithm for the automated detection of HFO
The S-Transform Gaussian Mixture Detector (SGM) consists of three steps. First, the algorithm detects EOI, that is, potential HFO in the high-frequency range. Subsequently, the S-Transform is computed for each EOI and meaningful features are extracted from the time-frequency distribution. Finally, a Gaussian mixture model (GMM) based on these features is fitted to cluster EOI into two groups: the HFO and the non-HFO group .

First stage: detection of events of interest.
For each channel, EOI were selected as those events above the baseline activity which there were putative HFO. A reliable threshold must be defined to detect these events. Most detectors in the literature have relied on measures extracted from the distribution of the data showing a critical dependence of the threshold on some events like the number of spikes, HFO, or artifacts [18][19][20][21]. To overcome this effect, the detection of EOI in our algorithm is based on an energy function, which is more robust to outliers and more stable to SNR changes [23,31]. In detail, we followed a similar procedure to the estimation used by researchers of the Montreal Neurological Institute [31]. The steps for this first stage were as follow, for each channel we: (a) Band-pass filtered the signal (80 to 250 Hz) using a finite impulse response filter of 80th order [21]. (b) Divided the filtered signal into 50-ms segments and calculated the wavelet entropy of the autocorrelation for each segment. Computed baseline as the 95 percentile of the entropy value distribution [31]. (c) Used Hilbert transform to compute the envelope of the band-pass signal. (d) Detected events when the envelope exceeded the 99.99% of the baseline value [31]. The duration of the event was then set as the interval between the upward and downward crossings through half the detection threshold. (e) If the duration of the event was at least 6 ms, the event was qualified as EOI. (f) Merged EOI within 10 ms of each other. (g) Rejected EOI without a minimum of six peaks higher than half the detection threshold.

Second stage: S-Transform and feature extraction.
The S-Transform [32], also known as Stockwell transform, is a generalization of the short-time Fourier transform, with a window that is function of frequency and no cross-term issues. The variable window provides a good compromise between time and frequency resolutions and reduces the computational cost. The S-Transform has been previously used to distinguish HFO among other similar activity [4,21]. To do so, HFO were considered shortlived events with an isolated spectral peak at a distinct frequency [21]. Following this idea, different features were defined considering their potential usefulness for measuring frequency and time short-lived events. The steps for this second stage were as follow, for each EOI we: (a) Selected signal segments considering 500 ms before and after the EOI (EOI ± 500 ms) (b) S-Transformed each signal segment (st). (f) Selected the region of interest (ROI). Among all the connected components in the binary image, the ROI corresponded to the region including the (max_t, max_f ) coordinates. (g) Obtained the following features for the selected ROI: (1) Area, that is, the number elements (pixels) in the ROI. (2) Entropy. Set all non-ROI connected components and the background to zero (black) and the ROI to one (white) and compute the Shannon entropy of the resulting image.

Third stage: Gaussian mixture model clustering.
A Gaussian mixture model (GMM) is a probabilistic model that assumes that a series of data points are generated by a mixture of a finite number of Gaussian distributions. GMMs do not require a priori knowledge of subpopulation membership of each data point, and for this reason they constitute a form of unsupervised learning. Each Gaussian distribution is determined by its covariance and mean. The covariance matrix determines the directions and length of the contours of the distribution, all of which are ellipsoids. GMMs perform a soft clustering, that is, the GMM algorithm assigns a probability of belonging to each cluster for each data observation. We used a full covariance matrix method assuming that all the components might adopt any position and shape. Furthermore, we set the number of clusters or Gaussians to two: one for the HFO group and another for the non-HFO group. Soft clustering allowed being more or less restrictive when assigning a data point to each group.
The following feature configurations were tested when using GMM-based clustering: Area. A measure of isolation, that is, if the event is isolated in time and frequency, as should HFO be; or if it covers larger regions, such as is the case of artifacts and spikes.
Entropy. A measure of irregularity of the binary time-frequency estimation. Although similar to area, it can measure some nonlinearities that might not be captured by area.
FW. Any HFO should have a narrow frequency width, while spikes and artifacts should have wider ones. FW was analyzed to see if this is enough to cluster HFO and non-HFO.
FW & TW. Although these features are highly linked with area, they can be useful to rule out cases with narrow frequency width but long duration in time. These cases could yield low area values but cannot be considered HFO. This configuration required 2-dimensional clustering, while all the previous used only one feature (1-dimensional clustering).

Performance evaluation
The performance of the current algorithm was evaluated using simulated and real data. The simulated signals and labelled events used in this work are freely available [23].
As mentioned before, for the assessment of real signals only those events selected as HFO by the two experts were considered as gold standard. For both real and simulated data analyses, a time window of 100 ms centered on each inserted HFO was defined as confidence interval (CI). All CI containing detections were considered true positives (TP), whereas CI without detections were false negatives (FN) and detections falling outside CI were labelled as false positives (FP). Percentage values of recall (R), precision (P) and F1-score (F1) were computed as: Recall measures the algorithm's ability to detect true HFO (minimizing the number of FN). Precision is related to the ability to reject events that might stand out from the baseline but are not HFO (minimizing the number of FP). On the other hand, the F1score is the harmonic average of precision and recall and can be interpreted as a measure of accuracy.
For simulated signals under different SNR scenarios, the SGM algorithm was tested using the feature configurations previously mentioned. These configurations were also used for real signals to assess if the results obtained for simulated signals could be extrapolated to a real scenario, in particular to situations where channels with few or no HFO are present, and other types of artifacts might appear.

Analysis of the stages of the detection algorithm
3.1.1. First stage: detection of events of interest. Figure 1 shows an example of the original signal, the band-pass filtered signal, the corresponding S-Transform, and the binarized S-Transform for six different EOI detected by the first stage of the algorithm in real signals. The binary images show the selected ROI in red. Some HFO (a ripple, a fast ripple, and a spike + ripple) were visible in the non-filtered signal and were clearly distinguished from the baseline in the filtered signal. The time-frequency representation for these HFO showed a distinguishable peak clearly disconnected from low frequencies. This peak was selected in the binarized S-Transform as the ROI. Non-HFO events (a spike, a low-amplitude artifact, and a high-amplitude artifact) did not show visible oscillations in the original signal but they appeared as events that emerged from the baseline when filtered. The time-frequency representation did not show any peak standing out or separate from low frequencies.
The binarized S-Transform showed a much more extensive ROI for these events, covering a large part of the high-frequency spectrum. The tested features (area, entropy, FW, and TW) were low for HFO, while non-HFO events showed much higher values.
Simulated signals allowed the evaluation of the performance of the SGM algorithm in each one of the stages. Figure 2 shows the percentage of detected events for each analyzed SNR and for each of the provided channels. For high and medium SNR values (10 and 15 dB), the first stage of our detector exhibited high sensitivity values, around 90% and 95% respectively. As expected, this first-stage sensitivity decreased with SNR. For 0 dB, EOI had the same amplitude as the baseline and were not distinguishable from noise, infringing the definition of HFO. As expected, the sensitivity dropped considerably in such a low SNR scenario.

Second stage: feature extraction.
All EOI from the first stage underwent feature calculation (area, entropy, FW, and TW). Figure 3 shows an example of their distribution for real data. Green values belong to first-stage EOI that were labelled by the experts as HFO, and red values show those not labelled as HFO. All the selected characteristics for HFO showed lower values than for the non-HFO events. This was also the case for simulated signals in all assessed SNR values. Standard deviation values were also higher for non-HFO features, a result that was expected given that HFO are events isolated from the low-frequency spectrum and bounded in time and frequency. Non-HFO, on the contrary,

Third stage: Gaussian mixture model clustering.
A GMM with two groups using the full-covariance method was applied to each subject in the case of real signals, and to each SNR scenario in the case of simulated signals. By doing so, the singularities of each dataset could be considered. Without any prior assumption, each dataset could contain different types of HFO and artifacts that would lead to different distributions (as seen in figure 3). An event was classified as HFO when its posterior probability was higher for the cluster with lower feature value (mean value of the fitted Gaussian) and as non-HFO if its probability was higher for the belonged to the group with higher feature. A GMM can be used as a softclustering method tuned by assigning data points into each group using the posterior probabilities. However, to cluster between HFO and non-HFO groups, we used the highest probability of membership.

Performance evaluation 3.2.1. Simulated signals.
In terms of precision, the SGM algorithm showed better performance as SNR increased (see figure 4). In the low SNR scenarios, the algorithm could be accepting EOI that are not actual HFO. The noisy timefrequency distribution might present isolated peaks in high frequencies that could be interpreted as HFO, which would produce an increase in FP values and therefore a decrease in precision.
On the other hand, recall values differed between stages because FN values could arise from events discarded by the first stage, or from events passing the first stage but incorrectly labelled as non-HFO in the second and third stages. Table 1 shows the recall values for stages 2 and 3 of the SGM algorithm. Precision values did not differ between stages because precision is the capability of minimizing FP, in our case, non-HFO deemed HFO, which can only be done in the second and third stages.
For low SNR (0 dB), where the amplitude of HFO was similar to the background activity, the low performance was mainly caused by the first stage of the algorithm. The second stage was able to discern HFO from non-HFO in more than 50% of the cases, with almost a 70% in the case of FW & TW configuration. For high SNR (10 and 15 dB), the restriction imposed by the first stage was not so critical because values of recall were similar to those of the second and third stages. For high SNR (15 dB) the best feature configuration for recall values was area (p < 0.01), but entropy and FW & TW also showed high recall values FW was significantly lower than the three other features 15 dB (p < 0.01).
Global performance was evaluated in terms of F1score. For medium SNR (5 dB), all feature configurations showed performance (see figure 4), whereas for high SNR area provided better values than the rest of configurations (p < 0.01). Figure 5 shows the average performance of the SGM algorithm on 11 real subjects. Precision decreased from values higher than 94% (high SNR simulated data) to values ranging from 75%-80% depending on the feature configuration. This may seem contradictory, but this result could be expected as the number of FP increased due to events detected by SGM not labelled as HFO by the experts. These events were probably not as visually noticeable as the ones actually selected by the experts. To illustrate this effect,   we computed the ratio between the average envelope amplitude of TP and FP and the baseline they were significantly higher (p < 0.01) for TP (2.16 ± 0.82) than FP (1.59 ± 0.36).

Real signals.
On the other hand, figure 5 shows that recall values were very high. Area and entropy provided the best results (95% and 96%, respectively). With respect to F1-score, entropy and area showed the highest values (87% and 85%, respectively).
In addition to strong global performance, it is also important to ensure correlation between the results obtained from visual and automated analysis [33]. To measure this effect, Spearman rank correlation was computed between the automatic detection (SGM algorithm) and the manual selection (visual inspection by experts) for all subjects. Only entropy was considered, as it was the best performing feature (p < 0.01). All correlation values were significant (p < 0.001), with a correlation coefficient of 0.84 ± 0.07 (ranging from 0.73 to 0.94).

Main contribution of the study
This article presents a new method, namely the S-transform Gaussian Mixture detection algorithm (SGM), for the automated detection of high frequency oscillations (HFO), a promising biomarker for epilepsy.
Some detectors in the literature rely on thresholds that need tuning prior to application, and depend highly on the database or signal under study [21,31]. The SGM algorithm does not need any a priori definition, calculation or parameter tuning, allowing the detection in signals recorded from different subjects in different settings, with different levels of noise and artifacts.
By computing the baseline activity robustly against outliers, SGM is preferable than the methods using thresholds extracted directly from data such as standard deviations or percentiles. This procedure improves HFO detection in highly noisy channels, as well as in signals with a high number of epileptic events [23].
Time-frequency characteristics were obtained using the S-Transform. This variant of short-time Fourier transform provides good resolutions in time and frequency at a low computational cost, and it has proven its effectiveness in the detection of HFO [4,21]. In simulated signals with high SNR (10 and 15 dB), the area feature performed significantly better than the other three features tested (F1 -score around 95%). In real signals, entropy showed a better performance (F1-score around 90%). However, in both simulated and real data, entropy, area and FW & TW showed similar strong performance, with FW alone performing significantly worse. These three features (area, entropy and FW & TW) directly or indirectly measure the relationship between the time and frequency span of detected events. Other detectors in the literature have classified events into HFO and non-HFO measuring exclusively the frequency width [21,34], but our results suggest that the combination of time and frequency spans is an important HFO characteristic that can improve automatic HFO detection significantly.
Finally, a Gaussian mixture model (GMM) was used to cluster the data into HFO and non-HFO. Features area and entropy presented two clear peaks in their distributions (see figure 3), justifying the use of two clusters. GMM do not require previous modelling or labelling of the data, as supervised detectors do need [33]. Using such a non-supervised method makes no a priori assumptions on the characteristics, and clusters can be obtained considering only the actual data of each individual. Advanced computational intelligence methods can adapt better to different specific data distributions, and that is why they may be more suitable to deal with individual differences or the variability among databases.
The validation an automated HFO detector is not trivial. Simulated signals provide a controlled environment where true HFO are known a priori. However, simulated environments may not consider all the possible artifacts arising in real electrophysiological data. Furthermore, real HFO may vary in morphology and rate along subjects [15], making it necessary to test the automatic detection algorithm in real data.

Validation with simulated signals
We tested the SGM algorithm on a publicly available database consisting of simulated signals with varying SNR [23]. We expected low recall for low SNR (0 dB) because of the definition of HFO, which states that they should be events clearly distinguished from noise. Unsurprisingly, the percentage of detected events increased with SNR. The first stage of the SGM algorithm was based on the wavelet entropy of the autocorrelation, a procedure like the one used by the MNI detector [29]. In this way, the background activity estimation is unbiased and robust to outliers, and provides better recall results.
Decreases in SGM performance in low SNR scenarios were mainly caused the first stage, where EOI were identified, causing subsequent clustering errors (see table 1 and figure 2). In the second and third stages, recall values reached values up to 70% for 0 dB SNR, suggesting that even in very noisy environments the proposed clustering is capable of minimizing the number of FN. Lower precision values of SGME in the low SNR scenario could be explained by isolated noisy events in the time-frequency distribution that were incorrectly interpreted as HFO, increasing the number of FP. A possible solution to improve the SGM algorithm in low SNR scenarios could be the whitening method proposed in [35] before extracting the time-frequency characteristics. Even so, if the detection of HFO in low SNR scenarios is to be improved, the definition of HFO must be reviewed because the it requires an event that stands out from the baseline in the filtered signal above 80 Hz [36].

Validation with real signals
We tested the SGM algorithm on real signals from patients with focal epilepsy in which HFO were detected by experts. While in simulated signals there is complete certainty of which events correspond to HFO and which do not, in real signals HFO were labelled by experts who had to agree. Only those HFO labelled by both experts were considered as so, and presumably, the most explicit HFO were labelled and the ambiguous were discarded. While SGM obtained high precision values for simulated data (higher than 94% in all cases), precision dropped (75%-80% depending on the feature) for real signals (see figure 5). As experts were not probably labelling the ambiguous events, the number of FP increased. We confirmed our hypothesis that most ambiguous events (FP) presented lower amplitude with respect to their baseline than the most explicit events (TP). This effectively suggests that those FP were not selected by the experts because they did not stand out from the baseline as clearly as the ones that were indeed labelled.
On the other hand, recall was high for real signals (90%-96% depending on the feature). The SGM algorithm was highly effective in detecting all the events labelled by experts, minimizing the FN rate. Area and entropy provided the best results (95% and 96%, respectively). By observing F1-scores, and taking into account the results obtained for simulated signals, we conclude that entropy (87%) and area (85%) could be the best features fed to the GMM to cluster HFO and non-HFO.
Finally, we evaluated the agreement between the visual detection by experts and the SGM automated detection. We observed higher count values for automated detections in some channels, but the general trend of the detections for automated and visual inspections stayed the same. Correlation coefficients were high for all the subjects, suggesting that there was a strong correlation between the HFO selected by experts and the HFO detected by the SGM algorithm.

Comparison with other detectors
The SGM algorithm proved its effectivity in the detection of HFO both in simulated and real environments. The publicly available simulated signals were previously tested in [23]. The performance of SGM algorithm showed much better performance than the STE [18], SLL [19], HIL [20], and MNI [31] detectors, lower performance than the Delphos detector for low SNR, and a slight increase in performance with respect to the Delphos detector for high SNR conditions.
The Delphos detector distinguishes oscillations and spikes in a pre-whitened time-frequency representation by analyzing the time width and frequency spread of peaks above a threshold. Detections are classified as HFO if their frequency spread is similar to wavelet's and their time width is greater than the one of a Dirac impulse.
Delphos shows better performance in low SNR scenarios probably due to the pre-whitening of the time-frequency decomposition. On the other hand, SGM presented a slight improvement in high SNR scenarios that could be explained because of the use of non-supervised clustering (GMM in our case), given that these methods present a more flexible alternative than fixed models to classify these kinds of events. Unsupervised methods allow the model itself to discover differentiated patterns in the data and therefore it is not necessary to define a priori which model will be used (as happens with Delphos). On the one hand, patients with epilepsy have considerable variations in their electrical signals, and on the other, the measurement with different techniques can be a source of variability in the signal [37,38]. The use of unsupervised methods to adapt to the variations of each subject and technique constitutes a promising approach to this problem [39].
We were not able to test the Delphos algorithm on real signals because the algorithm was not publicly available. The real data used in this study was previously tested with an automatic detector based in artificial neural networks and radial basis functions (RBF) [27]. Recall and Precision values for RBF detector were 50% and 73.6% in the best-case scenarios. In all cases, values were significantly lower than SGM detector. Correlation values of SGM (ranging from 0.73 to 0.94) were significantly higher than RBF detector (ranging from 0.71 to 0.85).
Furthermore, the comparison with the visual detection made by experts was remarkable and higher than the comparisons documented for other detectors [21,31], even using advanced computational intelligence methods [33]. SGM overcomes some of the most common problems of the current HFO detectors, as it was designed following a three-stage approach. The stages were developed considering previous detectors and combining the most reliable and effective methodologies in the literature.
The first step of the algorithm provides a robust estimation of the background activity to delineate events of interest. Detectors based on thresholds derived from standard deviation [18,20,21] or percentiles [19] are very sensitive to SNR changes and to the number of HFO present in the data [23]. MNI, Delphos and SGM methods compute a baseline histogram in order to set the threshold. While Delphos applies a normalization in the time-frequency spectrum, MNI and SGM find baseline segments based on wavelet entropy. Using this entropy, a non-biased estimation of the baseline boundaries is obtained, preventing a higher weight of the high-SNR events. Furthermore, the baseline estimation and delineation of events of interest in the time domain can be more computationally efficient than the estimation based on time-frequency maps.
The second step of the SGM algorithm extracts features from a time-frequency estimation computed with the S-Transform. Time-frequency characteristics have been used by several HFO detectors [20,21,34,35,[40][41][42], as time-frequency maps allow HFO to be distinguished from other events that co-exist in similar frequency ranges such as artifacts or ringing filtering effects produced by epileptiform spikes [21]. The two main strategies used in the literature for timefrequency estimation are the Short-Time Fourier transform (STFT) and the wavelet transform (WT). In general, the multiresolution nature of the WT produces more accurate results when detecting epileptic events. However, STFT is far more computationally efficient [43] and therefore more suited for real-time analysis or long-term multichannel EEG recordings. The S-Transform used in this work preserves the simplicity of the STFT but includes the multiresolution strategy of the WT [44] Finally, the third step of SGM uses a computational intelligence technique to detect HFO. The main advantage of these detectors is that they do not require or rely on fixed parameters to differentiate groups that share the feature distributions of HFO from those feature distributions that correspond to other events present in the high-frequency spectrum. During the last few years, several detectors based on computational intelligence technology have been developed [8,[24][25][26][27]34], proving their superiority over classic detection methods. A major drawback of these detectors is that few of them provide a non-supervised solution [34]. Supervised methods require that the possible outputs of the algorithm are already known and a training phase using correctly labeled cases. In the case of HFO detection, manual visual labelling of multichannel signals by experts is a tedious and highly time-consuming task. This affects the results considerably, as pointed out by recent study in which the authors observed that interrater reliability of visually evaluated HFO was poor, making supervised methods prone to bias [45]. SGM uses the distribution of features clustered with a GMM to define both groups, HFO and non-HFO, without any need of a priori labelling.

Limitations and further work
GMM is a soft-clustering method. To cluster between two groups, the SGM algorithm used the maximum probability of belonging to a group. However, depending on the application one may be interested in maximizing the precision or the recall of the algorithm. The probability of belonging to one cluster or another could be adjusted to accomplish this purpose. For example, in the case of real signals, we may prefer to detect only the most predominant events to minimize the number of FP. This could be achieved by making the boundaries HFO more restrictive, accepting only events with a very high probability of belonging to that group, discarding all other events.
Furthermore, the SGM algorithm was modularly designed, consisting of three stages, each one independent from the others. The modular concept of the algorithm allows the improvement of each stage keeping the others intact. For example, to improve the performance in low SNR scenarios, the whitening time-frequency method proposed in [35] could be included for the detection of EOIs in the first stage. Other features (based on time, frequency or, for example cross-channel information) could be added to the GMM to improve the classification of HFO in the future. In addition, the GMM could consider the peak frequency of each HFO and use it to automatically separate high-gamma events, ripples, and fast ripples, always adapting to the particular data of each individual and not using predefined and tuned thresholds.

Conclusion
High frequency oscillations are promising biomarkers for epilepsy. HFO can be used to help delineating the epileptogenic zone or to analyze the changes on the dynamics of the epileptic brain. During the last years, the improvement in computational capacity and data storage has allowed the acquisition of longterm EEG signals, increasing the number of channels and sampling frequencies. In the case of invasive EEG recordings, signals can include more than 100 channels and last for weeks. In this context, the manual labelling of HFO events cannot be considered possible due to the large amount of data.
On the other hand, many automatic detectors are based on fixed thresholds that may depend on the characteristics of the input signals, which may vary between patients or between databases. Therefore, it is necessary to adopt automatic detectors that do not depend on fixed parameters but have the ability to adapt to the particular data of each dataset or individual.
In this study, we propose a new algorithm called S Transform Gaussian Mixture (SGM) based on time-frequency characteristics obtained from the S-Transform and a subsequent clustering with Gaussian mixture models. Other algorithms based on advanced computational intelligence methods require a previous model of the data (supervised learning), and they need labelled databases to train properly the model before application.
The SGM algorithm is non-supervised, does not require previous labelling by experts, and is not based on any threshold or parameter that needs a priori tuning. We tested its performance on simulated signals with varying SNR and real data, obtaining good precision and recall values in both scenarios. SGM has a good computational performance and does not require long computation times. Consequently, it can be used to characterize HFO in long-term invasive EEG recordings. The algorithm was developed using freely available software (Python along with MNE toolbox). We plan to make this algorithm publicly available in the near future so that it can be used freely and tested on other databases. Finally, its modular design makes it possible to improve any stage of the algorithm without changing the others, allowing its future testing and improvement.