Kurtosis-Based Detection of Intracranial High-Frequency Oscillations for the Identification of the Seizure Onset Zone

Pathological High-Frequency Oscillations (HFOs) have been recently proposed as potential biomarker of the seizure onset zone (SOZ) and have shown superior accuracy to interictal epileptiform discharges in delineating its anatomical boundaries. Characterization of HFOs is still in its infancy and this is reflected in the heterogeneity of analysis and reporting methods across studies and in clinical practice. The clinical approach to HFOs identification and quantification usually still relies on visual inspection of EEG data. In this study, we developed a pipeline for the detection and analysis of HFOs. This includes preliminary selection of the most informative channels exploiting statistical properties of the pre-ictal and ictal intracranial EEG (iEEG) time series based on spectral kurtosis, followed by wavelet-based characterization of the time–frequency properties of the signal. We performed a preliminary validation analyzing EEG data in the ripple frequency band (80–250 Hz) from six patients with drug-resistant epilepsy who underwent pre-surgical evaluation with stereo-EEG (SEEG) followed by surgical resection of pathologic brain areas, who had at least two-year positive post-surgical outcome. In this series, kurtosisdriven selection and wavelet-based detection of HFOs had average sensitivity of 81.94% and average specificity of 96.03% in identifying the HFO area which overlapped with the SOZ as defined by clinical presurgical workup. Furthermore, the kurtosis-based channel selection resulted in an average reduction in computational time of 66.60%.


Introduction
Epilepsy is a complex and heterogeneous neurological disorder which affects approximately 50 million people worldwide and 2.4 million people are diagnosed with epilepsy every year (World Health Organization, http://www.who.int/mediacentre/factsheets/fs999/en/).In selected drug-resistant patients, a surgical option can be offered after accurate identification of the Epileptogenic Zone (EZ), i.e. "the minimum amount of cortex that must be resected (inactivated or completely disconnected) to produce seizure freedom". 1In the absence of a single diagnostic technique capable of identifying the entire EZ directly, in clinical practice, the Seizure Onset Zone (SOZ), i.e. the area of the cortex from which seizures originate, is used as a surrogate of the EZ.However, a definitive marker (functional or structural) that can exactly delineate the SOZ and therefore optimize pre-surgical evaluation and reduce pre-and post-operative morbidity is still lacking. 24][5][6] HFOs are shortlasting oscillations of bioelectrical activity seen in intracranial EEG (iEEG) at frequencies higher than 80 Hz.They have been sub-classified in to ripples (80-250 Hz) and fast ripples (250-500 Hz), can be transient (burst-like) or continuous (steady-state), occur within or outside interictal epileptiform discharges (spikes), phase-locked to a stimulus or eventrelated but not phase-locked. 7The study of this paroxysmal oscillatory behavior could lead to a better understanding of mechanisms of seizure generation as well as potentially aid the development of better medical treatments, 8 hence the importance of developing algorithms for fast quantitative detection of HFOs.
In clinical practice, the recognition of HFOs still mostly relies on visual inspection by trained experts, a highly time-consuming and rater-dependent approach.Strategies for automatic identification of HFOs in iEEG have been recently proposed.The first quantitative method for HFO analysis and detection was introduced by Staba and colleagues 9 and is referred to in the literature as the Staba detector, or the root mean square (RMS) detector or the Short Time Energy detector; it consists in the estimation of the energy of the band-pass filtered signal, expressed as RMS, and in the identification of the events lasting more than 6 ms and with RMS values greater than 5 standard deviations (SDs) above the RMS mean of the whole time series.If such events contain more than 6 peaks with an amplitude greater than 3 SD above the mean value of the rectified bandpass signal, they are classified as HFOs.The authors reported sensitivity of 84%; however, if used unsupervised, it can be prone to false detection or to classifying paroxysmal epileptiform transients as HFOs.
To address this limitation, a three-stage unsupervised classification pipeline for the identification of HFOs has been proposed. 10In the first stage, the Staba detector is used to identify candidate HFOs; in the second, all the events with significant spectral similarity with the surrounding iEEG activity are discarded from candidacy; in the third stage, features are extracted from candidate events and an unsupervised clustering algorithm is implemented to determine the different classes (clusters) of transient oscillations.The method is valuable because it is fully automated and allows the identification of four distinct clusters in the 100-500 Hz, thus supporting the evidence of the existence of mixed-frequency events between the ripple and fast-ripple classes.A subsequent study 11 used a different approach based on features extracted from visually-marked HFOs from three patients to train a radial-basis artificial neural network and evaluated the performance of this method on further eight patients with sensitivity and specificity of 49.1% and 36.3%, respectively.The main advantage of this method was in allowing automatic identification of brain areas with high HFO rate, an information of high relevance for the identification of the SOZ.A two-stage method for automatic detection of fast ripples in the 250-600 Hz band has also been proposed, 12 based on the initial global detection of events of interest defined by energy increase in the defined band, followed by estimation of the energy ratio between high frequency (HF, 250-500 Hz) and low frequency (LF) bands for differentiating fast ripples from other classes of events such as spikes or artifacts.Sensitivity and specificity of the detector were 93% and 95%, respectively.An alternative method to estimate spectral features of HFOs based on the Hilbert transformation has been proposed to identify events of interest; subsequently, the Stockwell time-frequency transformation is applied to obtain instantaneous power spectral density measures that separate HFOs from other classes of events present in the iEEG; this method offered a specificity >90% in SOZ identification. 13A sensitivity of 78.5% and a specificity of 88.5% were found for a further unsupervised detector of HFOs in long-term iEEG 14 ; in this study, the authors combine the Staba detector to select the candidate events and an automatic method to reduce the number of false HFOs and to improve the specificity of the recognition for prolonged recordings in clinical settings.A further staged approach to automatic identification of the SOZ consisted in an amplitude-based initial detection of events of interest, in the extraction of time and frequency domain features in the iEEG and in a clustering method to isolate "true" HFOs from other activities. 15This latter study reported a sensitivity of 81% and a specificity of 96%.Finally, a recent study 16 evaluated a fully automatic HFO detector based on baseline detection, HFO validation by Stockwell transform 13 and artifact detection, in a large set of intraoperative electrocorticograms.The method had a positive predictive value of clinical outcome of 100% and a negative predictive value of 62% for fast ripples.
Current detectors still present inherent weaknesses, such as the lack of robustness in detecting HFOs and the lack of an effective gold standard when automatically detected events are compared against human eye, as reported in Roheri et al. 17 Moreover, a limitation of some of the pipelines used in the selection of the channels, on which the analysis is performed, is related to potential bias associated with the use of a priori knowledge of channels involved in seizure onset, with the inherent risk of neglecting relevant activity in remaining channels.Furthermore, the HFOs detection methods described above use sophisticated algorithms such as sequentially applied clustering and time-frequency transforms, which can be computationally demanding with long analysis time especially in case of stereo-EEG (SEEG) recordings, acquired from more than 100 contacts at high sampling rate for several days.
In this study, to address the issues raised above and to facilitate data reduction and identification of candidate channels for HFO analysis, we describe a method for the detection of HFOs in iEEG that uses spectral kurtosis as criterion for restricting the search of HFOs to a subset of relevant candidate channels with specific time-frequency properties.Kurtosis is a statistical measure of the "tailedness" or "peakedness" of the probability distribution of a real-valued random variable and a powerful tool in outlier detection, particularly when the number of outliers is unknown. 18Higher kurtosis values suggest that most of the variance is due to infrequent extreme deviations as opposed to frequent small deviations.This measure has been used in the source analysis of magnetoencephalographic (MEG) signal in patients with epilepsy. 19,20High kurtosis values in band-pass EEG signal are better predictors of seizure freedom than the SOZ and high-rate/amplitude and duration of interictal paroxysms. 21We hypothesize that iEEG contacts with persistent bursts of HF activity will also be characterized by high kurtosis values in a restricted range of frequencies and that preliminary selection of only "highly-kurtotic channels" for timefrequency analysis and HFO detection will result in drastic reduction in computational time while retaining high sensitivity.For an initial validation, we analyzed, in the ripples frequency band (80-250 Hz), the iEEG data acquired during SEEG procedures of six patients with drug-resistant epilepsy and at least two-year post-operative seizure freedom or residual presence of rare disabling seizures and evaluated the concordance of our findings with the clinical interpretation of ictal SEEG data and the three-dimensional (3D) location of resected brain structures.

Patients selection, electrode type and implantation sites
The patients whose data were used in the study had focal drug-resistant epilepsy and underwent diagnostic evaluation with SEEG methodology followed by resective surgery.Five of the six patients were admitted at the Claudio Munari Epilepsy Surgery Centre, Niguarda Hospital, Milan, Italy (NIG) and one at the Birmingham Children's Hospital, Birmingham, UK.Surgical implantation strategy was decided according to the method originally described by Talairach and Bancaud 22 and later refined by Munari et al. 23 The procedure has been recently integrated with advanced computer-aided imaging and surgical techniques, [24][25][26]  scans and co-registered with the pre-operative MRI.This study was a retrospective/secondary analysis of anonymized data obtained in the context of standard clinical practice and as such was authorized by the R&D department of the respective institutions.
In Table 1, patients' characteristics, pathology, implantation sites and outcome post-surgery defined according to Engel's classification 27 are reported.Anatomical details of the implantation plan of each patient are available in Appendix A.

SEEG recordings
For NIG patients, intracerebral SEEG was recorded from intracranial multichannel electrodes (DIXI Medical, 5-18 contacts; 2 mm length, 0.8 mm diameter; leads 1.5 mm apart).The number of electrodes and the sites for implantation were decided according to anatomical and clinical data collected during the noninvasive phase of the evaluation and varied between 5 and 18 contacts per intracranial electrode (maximum number of 192 recording channels).Bandpass filter of 0.016-500 Hz was used.EEG signal was acquired continuously with a Neurofax EEG-1100 system (Nihon Koden, Tokyo, Japan) and sampled at 1 kHz with 16-bit resolution.Each channel was offline rereferenced with respect to its direct neighbor (bipolar derivations with a spatial resolution of 3.5 mm) to cancel-out effects of distant sources that spread equally to both adjacent sites through volume conduction.When appropriate for diagnostic purpose, the iEEG signal was referenced to the average signal from two electrodes identified from anatomical and neurophysiological data to be in the white matter.Two scalp EEG channels (Fz and Cz referenced to a mastoid electrode) and chin electromyogram (EMG) were recorded in addition to SEEG for sleep staging.The simultaneous video-iEEG recordings lasted between 5 and 10 days.
Electrode implantation in the BCH patient was performed according to the same clinical methodology as NIG patients.SEEG recordings from 128 contacts were obtained using a commercial video-EEG monitoring system (System Plus, Micromed, Italy).Data were acquired with band-pass filter of 0.016-1 KHz and sampled at 2 KHz.The remaining recording parameters were the same as those used in the NIG patients.

Data preprocessing
Each continuous SEEG recording was reviewed by an expert neurophysiologist from the other institution involved in the study (RM and SS) who was blinded to the clinical data of the patient; the time of seizure onset and contacts involved in the onset of each seizure were annotated.Such contacts defined the anatomical boundaries of the SOZ.For each patient, data from two seizure episodes was analyzed by one of the authors (LQ) blinded to the clinical information and to the results of the diagnostic procedure.An EEG segment containing 10 min before and 2 min after the electrographic onset of each seizure was extracted and was considered as the period of interest for HFOs identification.Data were resampled to 1024 Hz (linear interpolation) and then visually inspected to identify artifactual channels, which were disregarded in successive analyses.A bipolar montage between adjacent contacts was built and data were filtered in the 80-250 Hz frequency band.The original epoch was finally segmented in 2 min long sub-epochs, which were used for further analyses.Resampling and filtering were performed by means of functions in the EEGLAB suite. 28

Software framework
A full set of algorithms and routines for the detection of HFOs in iEEG and source-space MEG data, the evaluation of their spectral properties and the identification of the candidate SOZ were developed and coded in the MATLAB programming language (Matlab 2016b, the Mathworks Inc., Natick, MA).The code was integrated into the EEGLAB interface (sccn.ucsd.edu/eeglab)as an extension plugin named EPINETLAB, a multi-GUI user-friendly framework; it consists of three main blocks of functions summarized in Fig. 1: (i) time-frequency analysis of iEEG and MEG data, descriptive statistics of signal properties and kurtosis-basis thresholding; (ii) HFOs detection; (iii) performance evaluation of HFOs detection and comparison with the clinically defined SOZ.
Additional optional functions for data preprocessing and multifile automatic HFOs detection were also embedded in the framework.Int.J. Neur.Syst.Downloaded from www.worldscientific.comby ASTON UNIVERSITY on 04/04/18.For personal use only.

Channel selection
For each bipolar channel, a time-frequency transform was computed on consecutive 1 s windows and in the frequency interval of interest (80-250 Hz), by means of a complex Morlet wavelet, an established method for time-frequency analysis of iEEG data 29 and defined as follows: where f c is the wavelet central frequency and f b is the bandwidth parameter.The wavelet scaling limits were determined dividing f c by the band frequency limits (80 Hz and 250 Hz) and by the time resolution (1/sampling rate), while the scaling resolution was set to 0.05; the relationship between scales bins and frequency bins was determined by using the function scal2freq provided in the Matlab Wavelet toolbox.
For each channel and each window, we computed the scalogram, which represents the percentage of energy of each wavelet coefficient and then, for each frequency bin, the relative spectral kurtosis 30,31 was computed in the time domain.Spectral kurtosis is a statistical measure used to detect non-Gaussian components in a time series; it is considered to reflect the presence of transients and can identify their locations in the frequency domain.To investigate the distribution of kurtosis over all frequencies and channels in relation to the presence of transient events, the distribution of average values of spectral kurtosis was fitted against a set of known distributions available in the Matlab "Statistics and Machine Learning Toolbox" (e.g.normal, exponential, gamma, generalized extreme value (GEV), etc.).The distribution ranking first in terms of logarithmic likelihood was selected as representative of the kurtosis in the specific iEEG segment and its mean and variance values were computed.A threshold on the kurtosis was set as the value that exceeds the mean value of the fitted distribution by 3 SD.In the second step, channels were ranked according to the total number of windows in which spectral kurtosis peaks could be clearly identified and exceeded the set threshold in a restricted range of frequencies (we empirically set the value to 20 consecutive frequencies).This step allows events with high kurtosis values distributed over all the frequency bins to be discarded, as they are likely to be associated with transient activities spreading over all the frequency bands of interest, typical, for example, of interictal spikes or even nonphysiological events (e.g.artifacts).This procedure ensured that only windows with transient band-limited activity in the HFO range contributed to the ranking.
After the ranking, the median (Q2) and the 75th percentile (Q3) were computed from the distribution of the number of windows.Only channels with a number of windows >Q2 + (Q3 − Q2)/2 were retained for further analysis (the term (Q3 − Q2)/2 was added to take into account the positive skewness of the distribution).

HFOs detection
On the most informative channels, candidate HFO bursts were defined as events in which the power of the wavelet coefficients calculated over 3 ms-long consecutive windows exceeded the mean power in the whole 1 s window by 5 SD and for more than 20 ms. 32andidate events with an inter-event-interval of less than 10 ms were collapsed into a single event.

Artifacts
After the recognition of the candidate events, in order to exclude from the analysis of nonphysiological events which propagate to all or to a large number of recorded channels violating the focal hypothesis of HFOs, and to remove wide-band activities, which violated the hypothesis of HFOs' limited frequency content, two criteria were adopted: (1) all the candidate events that were synchronously located on more than three channels were considered as artifacts and disregarded; (2) all the candidate events with power fitting a normal distribution over all the frequencies of interest were considered as artifacts and disregarded.This choice was based on the suggestion that a "real" HFO is represented by an isolated peak in the time-frequency plot, while other events such as interictal spikes generate a more elongated peak which is extended in frequency. 33,34This procedure removes the need for visual inspection of the raw data if exclusion of epochs containing interictal spikes is considered appropriate.

Definition of the HFO area
The definition of the channels to be included in the HFO area was performed by considering HFO-rate data of the two seizure episodes.Channels were first sorted based on their HFO rate, and then subject to further analysis using the following five methods: (a) "Max.Value" method: the five bipolar channels of each episode with highest ranking in terms of HFO rate were selected as HFO area.Five is chosen as a reasonable compromise with respect to the assumption of focal epilepsy; (b) "Tukey's fence" method: bipolar channels of both episodes and with an HFO rate higher than the upper Tukey's fence (3rd quartile + 1.5 times the inter-quartile range (IQR)) of all HFO rates were selected.This method is typically used when outliers need to be identified in a dataset.Channels within the HFO area can be considered as data with anomalously high HFO rates; (c) "Fuzzy c-means clustering" method: a fuzzy c-means clustering 35 was implemented to find two natural clusters in the HFO rate dataset which the individual HFO rate belongs to with a certain degree.Channels belonging to the cluster with the highest HFO rates and with a probability of belonging to that cluster higher than 0.7 were selected from each seizure episode and then merged; (d) "K-means clustering" method: a k-means approach 35 was applied to find clusters in the HFO rate dataset.Channels belonging to the cluster with the highest HFO rates were selected from each seizure episode and then merged; (e) "Kernel Density Estimation (KDE)" method: a KDE (with Gaussian kernel) of the distribution of HFO rates was performed. 34Then, after a smoothing procedure to reduce oscillations in the distribution, its peaks and troughs were detected and the cutoff value on HFO rate was set to the occurrence of the lowest trough.The HFO rates exceeding that threshold were selected for each seizure episodes and then merged.

Comparison with the SOZ
The channels corresponding to electrode contacts in the HFO area and localized into the SOZ were defined as true positives (TP); the channels within the HFO area but not located in the SOZ were defined as false positives (FP); the channels outside the HFO area and not located in the SOZ were defined as true negatives (TN) and those located outside the HFO area and into the SOZ were defined as false negatives (FN).The sensitivity of the detection was defined as Sens = TP/(TP + FN), the specificity as Spec = TN/(FP+TN) and expressed in percentages.Confidence intervals were estimated using a binomial distribution.To summarize the performance of the five methods used to define the HFO area, average sensitivity and specificity were computed across patients as well as the corresponding Youden index. 36This index is a commonly used metric of the performance of diagnostic tests and was calculated according to the formula (Sens + Spec − 100)/100.3D rendering of implantation images was used to visually confirm the topographic relationship between electrode contacts, the SOZ and the resected brain area.The imaging pipeline was implemented referring to a recently described method 25 and subsequent modifications.Based on the pre-operative MR data and post-implantation computer tomography scans, image fusion was performed to locate each lead on the recording electrode trajectory.The spatial location of the SEEG contacts in the preoperative MRI space of each patient was obtained using the recent software implementation SEEG Assistant. 37The post-operative MRI was then coregistered to the pre-operative MRI to evaluate if the candidate HFO contact was within the resected volume.Co-registered data were displayed using 3D Slicer, an open source platform for medical image informatics, image processing, and three-dimensional visualization. 38

Kurtosis distribution and HFO detection
The kurtosis in the ripple frequency band fitted a GEV distribution. 39An example of such distribution is reported in Fig. 2   event and bottom graph shows the spectral kurtosis in the same window.In the panel 3(a) case, the window contributes to the computation of the total number of windows on which to rank channels, as the spectral kurtosis has a clear maximum (in the 80-160 Hz range), and exceeds the threshold (32) for a limited amount of frequencies.In panel 3(b), instead, an example of window which does not contribute to channel ranking is reported.A clear maximum cannot be detected in spectral kurtosis, which is uniformly higher than the threshold (32) over all the frequencies in the band of interest.Visual inspection of the original trace confirmed that this was due to a spike discharge.For illustrative purpose, ranking of channels based on the number of windows exceeding the kurtosis threshold of 32 is shown in Fig. 4. In Fig. 5, examples of HFO detection in patient 1 on channel T5-T6 are shown for a representative epoch.Both false and true HFO events are illustrated.

HFO area identification and comparison with the SOZ
In Table 2, the results of HFO area identification in each patient, expressed in terms of sensitivity and specificity of SOZ recognition, for the five channelselection methods are reported.Tukey's fence-based and k-means clustering approaches achieved the 1850001-9 Int. J. Neur.Syst.Downloaded from www.worldscientific.comby ASTON UNIVERSITY on 04/04/18.For personal use only.highest performances in terms of sensitivity and specificity, with the first one characterized by slightly higher overall performances.For patients 2, 3 and 6, Tukey's method allowed the identification of all the channels of the SOZ (sensitivity of 100%) with high specificity (96.19% for patient 2, 94.94% for patient 3 and 100% for patient 6).The chosen method allowed an accurate identification of the SOZ while keeping the FP low.For patient 5, the method achieved comparable sensitivity to the other methods (75%) with high specificity (92.22%).For patient 1, Tukey's method achieved sensitivity lower than clustering methods (50% versus 66.67%), but with a higher specificity (96.63% against 84.09% of fuzzy c-means clustering and 80.68% of k-means clustering) and for patient 4, the performances were similar for both Tukey's and clustering methods.The "Max.Value" method obtained high performances in terms of sensitivity and specificity, while clustering methods allowed to achieve the highest sensitivity in patient 1 (66.67%), to the detriment of a high number of FPs and of a lower specificity (84.09% for fuzzy clustering and 80.68% for k-means one).This is not a desirable result as a high number of FPs could lead to the resection of not-epileptogenic areas.The same was found for patient 2 where, despite a sensitivity of 100%, a higher number of FPs lowered the specificity (87.62% with fuzzy clustering and with k-means).
For patients 4 and 5, clustering methods were equivalent while for patient 6, sensitivity of 100% was obtained together with specificity of 100% for fuzzy-c and 87.50% for k-means.KDE distribution-based method showed highly variable performances: it was associated with some of the highest specificities for all the patients (92.13%, 97.17%, 100%, 100%, 100% and 100%), but lower sensitivities (50%, 66.67%, 40%, 16.67%, 50% and 100%), due to lower number of TPs, as already demonstrated by Gliske et al. 34 The advantage offered in all the subjects in terms of increased specificity was offset by a high number of FNs, making the method too conservative for the correct identification of the SOZ.
In Table 3, average sensitivity and specificity across all subjects are reported for each HFO area identification method, together with the relative Youden index.Tukey's fence method shows the highest overall performance with marginally higher Youden index values than k-means.In Table 4, the SOZ identified in the clinical pre-surgical workup is compared to the HFO area (channels in the HFO area are sorted according to the HFO rate) as detected with the Tukey's fence method.Please note 1850001-10 Int.J. Neur.Syst.Downloaded from www.worldscientific.comby ASTON UNIVERSITY on 04/04/18.For personal use only.Notes: T = true; F = false; P = positive; N = negative; Sens = sensitivity; Spec = specificity; CI = confidence interval.
that, for patient 1, channel L3-L4, and for patient 4, channels R6-R7 and R7-R8, were included by clinicians in the SOZ but had been excluded from further analysis due to intermittent artifacts in the pre-ictal and ictal periods.For each channel in the Tukey's fence HFOs area, the relative HFO rate, normalized over the total number of detected HFOs, is reported.
The relationship between the resected area containing the SOZ shown on the post-operative MRI and contact positions for one of the patients (patient 2) is shown in the top panels of Fig. 6.Electrode Z with its contact Z11 (the most inferior and posterior contact of the SOZ) can be seen against the resected area (left upper panel); the TP contacts are highlighted in yellow.In the right upper panel, in orange are highlighted the FPs (B1-B2, P12-P13, P11-P12 and P13-P14).The lower panels of Fig. 6 show the three electrodes on appropriately chosen MRI slices.
Table 4. Clinical SOZ defined in presurgical work-up and HFOs area defined according to Tukey's fence.

Comparison to the Staba detector
Performance of the newly developed detection method was compared with that of the Staba detector 9 on the channel identified as most informative by the kurtosis-based selection.The Staba detector computes the RMS of the signal of interest in 3 ms sliding windows and defines as candidate HFOs segments with RMS values at least 5 SD above the mean amplitude of the RMS signal lasting more than 6 ms.Candidate events closer than 10 ms are collapsed into a unique event.The final condition to be met is the presence in the candidate HFO of at least 6 peaks which are greater than 3 SD from the mean value of the rectified band-pass signal.Performance of HFOs detection by means of Staba detector is reported in Table 5.The selection of the channels to be inserted in the HFO area was done by means of the Tukey's fence-based method.A considerable decrease in sensitivity is obtained in all patients but patient 6 (average sensitivity 39.17%), with patients 1 and 2's SOZ not recognized at all.Average specificity of 92% was achieved.

Computational load
The selection of a subset of channels exceeding the spectral kurtosis threshold to be fed to the analysis process significantly reduced the computational load (66.6% reduction).This is a highly desirable feature in long-term recordings performed with high spatial density and sampling rate.In the absence of a priori information that would reliably allow discarding specific channels from further processing, a reduction in computational time linked to wavelet and kurtosis computation is only achievable through improvement in coding or the use of routines different from those available in the proprietary MATLAB libraries.
In Table 6, the total number of bipolar channels, the average number of channels retained after kurtosis thresholding and the average number of channels on which HFOs were effectively detected (computed on six 2 min epochs) are reported for each patient and for one analyzed seizure episode.A significant difference was found between the number of channels retained for the analysis after kurtosis thresholding and the number of channels on which HFOs were effectively identified (paired t-test, alpha = 0.05, p = 0.004).This suggests that the kurtosis-based thresholding was effective in selecting a subsample of channels from the whole dataset; this finding and the high sensitivity of the methodology suggest that this did not result in loss of useful information.As far as computational time, detection of HFOs and discarding

Discussion
In this study, we describe a novel method for the detection of HFOs in the ripple frequency band (80-250 Hz) of the iEEG and the identification of the SOZ developed as a complete plugin for the EEGLAB suite.The aim was to provide clinicians with an easy to use GUI-based tool integrated into a widely used analysis suite.The software supports the identification of the contacts located in the area most probably involved in the generation of seizures and could help planning the resection strategy.In this line, we tested the method in a naturalistic scenario, using data collected in the process of pre-surgical evaluation of patients with drug resistant epilepsy studied in two European epilepsy surgery centers.The main property of our method is that it does not rely on a priori knowledge about the localization of the seizure in the implanted electrodes, but uses a kurtosis-based metric to select the subset of channels most probably related to SOZ.Kurtosis has been previously used in clinical studies for first-pass selection of virtual electrodes in beamformer-based analysis of MEG signal, 20 but to our knowledge, it has not yet been exploited to identify HF activity in SEEG recordings.Its main advantage is to provide a metric for the selection of a subset of informative channels and to reduce computational complexity, an essential requirement when dealing with multichannel and long-term SEEG recordings at high sampling rate.
In the implementation described here, the method was applied to HFOs outside interictal transients such as "epileptic" spikes or electrode artifacts.This was made possible by combining wavelet transform for the identification of intervals in the signals with power peaks in a limited range of frequencies, followed by a rejection of events with power fitting a normal distribution or synchronously occurring on multiple channels (see Fig. 4).In the software release, we have left this as computational option, since there is still debate as to whether HFOs within interictal spikes have specific diagnostic value different from HFOs occurring in isolation. 5esults of the methodology applied to the sample used for this initial validation are promising.Beside an average reduction of the computational load of 66.60%, the method achieved an average sensitivity of 81.94% and a specificity of 96.03% when the Tukey's fence method was used to select channels in the HFO area.This method was associated with the highest Youden index.One of its possible limitations emerges when only modest differences between HFO rates in different channels are found (lack of outliers).Nevertheless, it was still the most robust across subjects when compared to the "Max.Value", which is not influenced by the distribution of data, and can dramatically fail in case of particularly small or large SOZ and to clustering methods which, despite the high sensitivity, were characterized by a low specificity, this possibly due to class displacement (channels in HFOs area total channels).KDE-based method was the one with highest specificity, to the detriment of sensitivity, which was among the lowest.This results in few channels being identified within the HFOs area.
Despite these valuable properties, some limitations do exist; first, the detection of HFOs is done on a channel-by-channel base, while the identification of the SOZ in the pre-surgical workup includes a wider range of input data and the support of different technologies.Some channels which are identified by the method as being part of the SOZ, due to high HFO rates, could be just secondarily involved in the seizure and not causally responsible for the ictal event.Furthermore, the method could fail in identifying all the channels belonging to the SOZ, in case of multifocal of regional sources of seizures associated with a wide SOZ.
The method was tested on pre-and ictal epochs since these have been described as periods of increased HFO activity by Zijlmans et al. 40 Unlike interictal spikes, ictal and interictal HFOs co-localize to the epileptogenic area and probably represent similar phenomena.While the complex issue of ictal versus interictal HFOs is beyond the aims of this study, the statistical properties of the HFO distribution suggest that the algorithm could be successfully applied to interictal epochs as well.
The method was tested on the ripple frequency band (80-250 Hz) only.This was due to the retrospective nature of the study and the need to select patients who had a sufficiently long period of postoperative seizure freedom to evaluate the goodness of the iEEG prediction on the topography of SOZ, which resulted in the NIG patients to be acquired with an older EEG system with maximum sampling rate of 1 KHz.Higher sampling rates (≥2 kHz) are recommended for the detection of HFOs in fast ripple 1850001-14 Int.J. Neur.Syst.Downloaded from www.worldscientific.comby ASTON UNIVERSITY on 04/04/18.For personal use only.
band. 14 Because of this initial choice, we could not validate the method on fast-ripple HFOs in the 250-500 Hz range and which are also worth investigating as they have been reported to be highly correlated to positive surgical outcome. 41This limitation will be overcome with the availability of a wider group of patients from other collaborating laboratories in a prospective study that is under process.In further developing the method and the EEGLAB plugin, the larger number of patients that will become available will allow us to test the suitability of the method to a wider range of etiologies and possibly perform group analysis.This will be part of the statistical functionalities which will be released in future versions of EPINETLAB.

Conclusion
A novel method to aid detection of HFOs for the identification of the SOZ in iEEG is presented here.Essential metrics of the method were validated with data recorded in the context of invasive pre-surgical evaluation without a priori selection of the channels

Fig. 1 .
Fig. 1.EPINETLAB functional scheme.EPINETLAB consists of three main blocks of functions: (1) time-frequency analysis and kurtosis-based statistical thresholding: this includes data segmentation, computation and display of time-frequency analysis with wavelet transform, kurtosis estimation and statistical thresholding for channel reduction; (2) HFOs detection: this includes identification of HFOs and display in both time and time-frequency domains; saving detection results in text files and (3) performance evaluation: this includes computation of HFO rates for each channel, display of HFO rate histograms, identification of HFO area, comparison of detected HFO area with the clinically identified SOZ and evaluation of sensitivity and specificity.Additional optional functions are released within the framework for data preprocessing (bipolar montage creation, file cutting, channels ordering) and multifile automatic HFOs detection.
for a representative patient.The values in the right tail of the distribution represent those windows with the highest kurtosis values, most probably an expression of abnormal brain activities.Figures 3(a) and 3(b) show examples of spectral kurtosis distribution in two windows of interest and the relative iEEG signal for the same representative patient of Fig. 2. In both panels, top row graph shows the raw iEEG channel, second graph shows its bandpass filtered version, third graph shows the relative time-frequency transform at the occurrence of an 1850001-8 Int.J. Neur.Syst.Downloaded from www.worldscientific.comby ASTON UNIVERSITY on 04/04/18.For personal use only.

Fig. 2 .
Fig. 2. (Color online) Distribution of kurtosis for a representative patient (1).Kurtosis was computed on 155 channels and 120 windows and averaged over all the frequencies.The fitted distribution was a generalized extreme value distribution.The threshold, indicated by a green bar, was set to 32.

Fig. 3 .
Fig. 3. Examples of windows (a) contributing and (b) not contributing to channel ranking and the relative spectral kurtosis.Top row: raw iEEG; second row: filtered iEEG (80-250 Hz); third row: wavelet transform; bottom row: spectral kurtosis distribution over the frequencies of interest.In the example, in panel 3(a), a clear peak in spectral kurtosis is detectable; also the kurtosis threshold value 32 is exceeded in a limited interval of frequencies.In panel 3(b), a clear peak in spectral kurtosis cannot be identified; also the threshold value 32 is exceeded in a broad frequency interval.

Fig. 4 .
Fig. 4. (Color online) Ranking of the channels based on the number of windows exceeding the kurtosis threshold of 32 in patient 1 for a representative epoch.Red dotted line at 8 indicates the minimum number of windows needed for each channel to be retained in the HFO detection.The channels indicated by blue arrows are in the SOZ.

Fig. 5 .
Fig. 5. Examples of HFO detected on a contact in the SOZ (T5-T6) in patient 1. Yellow box (a) indicates a candidate HFO event rejected due to the normal distribution of power estimates.Red box (b) indicates a candidate event rejected because occurring on more than three channels at the same time.Green box (b) indicates a true HFO event.

Fig. 6 .
Fig. 6. (Color online) Top: representation of the topographic relationship between contact positions and resected zone from the post-operative MRI for one representative patient (patient 2).Green points represent contacts.Electrode Z with its contact Z11 (the most inferior and posterior contact of the SOZ) can be seen against the lower part of the resected area (left panel); the TP contacts (Z11-Z14) are highlighted in yellow and pointed by yellow arrows.In the right panel, the FP (B1-B2, P11-P12, P12-P13 and P13-P14) are highlighted in orange.Bottom: topographic display of TP (Z) and FP (B, P) contacts on appropriate slices from the post-operative MRI with the same color coding as above.Co-registered data are displayed using 3D Slicer.

Table 1 .
Patients' information and implantation sites.

Table 2 .
Results of the HFO area identification with five different methodologies and comparison with the SOZ.

Table 3 .
Average sensitivity, specificity and relative Youden index of methodologies for HFO area definition.

Table 5 .
Results of the detection of HFOs with the Staba detector and the identification of the HFOs area with the Tukey's fence.

Table 6 .
Computational load of the kurtosis-based channel selection for each patient on six 2 min epochs.
RAM: 32GB).This corresponded to a reduction of average 1260 s (± 340.84 s) in the analysis of the two seizure episodes (six 2 min epochs for each episode).1850001-13Int.J. Neur.Syst.Downloaded from www.worldscientific.comby ASTON UNIVERSITY on 04/04/18.For personal use only.

Table A .
1. Anatomical localization of the entry point (lateral structure) and the target (mesial structure) of each implantation site for each patient.