Automatic Identification of Interictal Epileptiform Discharges in Secondary Generalized Epilepsy.

Ictal epileptiform discharges (EDs) are characteristic signal patterns of scalp electroencephalogram (EEG) or intracranial EEG (iEEG) recorded from patients with epilepsy, which assist with the diagnosis and characterization of various types of epilepsy. The EEG signal, however, is often recorded from patients with epilepsy for a long period of time, and thus detection and identification of EDs have been a burden on medical doctors. This paper proposes a new method for automatic identification of two types of EDs, repeated sharp-waves (sharps), and runs of sharp-and-slow-waves (SSWs), which helps to pinpoint epileptogenic foci in secondary generalized epilepsy such as Lennox-Gastaut syndrome (LGS). In the experiments with iEEG data acquired from a patient with LGS, our proposed method detected EDs with an accuracy of 93.76% and classified three different signal patterns with a mean classification accuracy of 87.69%, which was significantly higher than that of a conventional wavelet-based method. Our study shows that it is possible to successfully detect and discriminate sharps and SSWs from background EEG activity using our proposed method.


Introduction
Epileptiform discharges (EDs) are signal patterns frequently observed in interictal electroencephalogram (EEG) of patients with epilepsy [1]. EDs are generally utilized by expert epileptologists to assist in the diagnosis of epilepsy [2]. However, manually labeling various types of EDs is time-consuming and labor-intensive because EEG data are generally recorded from many electrodes during continuous monitoring of patients.
For this reason, the automated detection of EDs has drawn researchers' attention. Automatic detection algorithms introduced in the literature mostly focus on either sharpwave (sharp, sometimes referred to as spike) or sharpand-slow-wave (SSW) discharges. Epileptiform sharps are detected using a variety of algorithms such as deterministic finite automata (DFA) [3], geometrical features and artificial neural networks [4], cross-correlation [5], and dynamic time warping [6]. A recent review reported that these algorithms can detect epileptiform sharps with approximately 90% accuracy [7]. Algorithms for the automatic detection of SSWs have also been introduced in the literature. Although the number of studies to design algorithms for SSW detection is relatively less than that for sharps detection, a variety of approaches have been used for the automatic detection of SSWs, such as Fourier transforms [8], wavelet transforms [9][10][11], Gotman and Gloor's spike detection method [12], a rulebased system [13], and a detection method using a number of geometrical features [14]. It was also shown that sharp-andwave discharges can be distinguished from sleep spindles and background activities using wavelet-based methods [9,[15][16][17][18].
There have been few studies to detect different types of EDs together without identifying their individual types. D' Attellis et al. detected EDs including spikes and spike-wave complexes using a wavelet transform [19]. They also detected three types of EDs (spike, slow-wave, and spike-and-slowwave) simultaneously by using wavelet transform coupled with a multilayer perceptron model [20]. Bergstrom et al. proposed a method utilizing total wavelet decomposition and signal variation to identify four types of EEG segments: normal, seizure, spike, and other abnormal types [21]. They successfully distinguished spike and seizure segments from background activity on EEG; however, the accuracy for identifying other abnormal types was relatively low (approximately 44% were classified as a normal segment).
In this study, we propose a new method based on a novel wave detection algorithm for the detection and identification of different types of EDs. The proposed method was evaluated with intracranial EEG (iEEG) data recorded from a patient with Lennox-Gastaut syndrome (LGS), a type of secondary generalized epilepsy in which the patient suffers from multiple seizure types, therefore generating a variety of EDs on iEEG [22].
LGS is known to be treatment resistant, but recent studies have shown that surgical treatment is possible for some patients with LGS that have focal epileptogenic zones [23]. However, because of their generalized ictal iEEG discharges, surgical resection areas are determined mostly based upon interictal iEEG patterns and functional neuroimaging results [23].
Although the automatic identification of interictal EDs is of great importance for surgical planning in patients with LGS, no studies have attempted to develop algorithms for detecting and identifying interictal EDs. Among the various interictal EDs, we have focused on repeated sharps and runs of SSWs, as these forms are hard to detect with existing algorithms. These characteristic EEG patterns have great potential to aid in surgical planning of patients with secondary generalized epilepsy. The performance of our proposed algorithm was evaluated and compared with that of conventional frequency domain approaches.

Materials and Method
2.1. Materials. We used iEEG data recorded from a pediatric patient with LGS (17-year-old male) at Severance Children's Hospital. The data was recorded at a sampling rate of 1600 Hz with a band-pass filter at 1-500 Hz. A 60 Hz notch filter was also applied before the analog-to-digital conversion. As shown in Figure 1, 128 grid electrodes were placed over the cerebral cortex of the right hemisphere to cover the prefrontal, frontal, and temporal lobes.
A medical doctor manually labeled the types of EDs and a biomedical scientist who was not involved in developing the algorithm in this study adjusted the ranges. We extracted fivesecond epochs including repeated sharps and runs of SSWs from the background activity. We identified 73 segments of repeated sharps, 1,752 segments of SSWs, and 91,701 normal segments without any specific EDs ( Figure 2).

Wavelet Features.
Wavelet transform is one of the most popular approaches for the analysis of EDs [9-11, 19, 20]. In this study, we adopted techniques introduced byÜbeyli et al. which were simple to use while providing promising results [10]. The features used were mean, maximum, minimum, and standard deviation of the wavelet coefficients in each scale. A wavelet coefficient ( ) of data at scale and position is defined as where represents a Morlet mother wavelet. The number of scales was empirically set to 10 in this study.

Detection of Sharps and SWs.
We have developed an algorithm for detecting sharps and SWs. Although there is a conventional method for the detection of SSWs [14], it assumed that a sharp always accompanies a SW. Though this assumption is useful in distinguishing SSWs from normal EEG data, it is not appropriate for classifying repeated sharps and runs of SSWs because both types of EDs commonly include sharps.
We designed an algorithm that distinguishes sharps and SWs based on the sharpness of the waveform pattern (Figure 3). To determine the range of a wave that can potentially be either sharp or a slow-wave, a line is drawn from a local maximum to a moving point (MP) near the maximum. Then, as shown in Figure 3(a), two areas can be calculated for each MP, which are (1) area under the signal and above the line (denoted by AAL), and (2) area above the signal and under the line (denoted by AUL). The AAL and AUL from a local maximum to an MP are calculated with the following equations: where { | } denotes a set in which elements are generated by a function while a condition is satisfied, ∑{⋅} denotes the summation of all the elements in a set, and denote the data indices of the local maximum and MP, respectively, is the th data point of the signal, and is a normalization factor. The purpose of AUL is to find candidates for two  boundary points (leftmost and rightmost points) of a wave.
To search for the candidates of the current wave under consideration, the MP is gradually moved outward (to the left or to the right) from the local maximum. A point is determined to be a candidate of the wave boundary if the AUL at the point is smaller than the preset threshold ( AUL ) and the AUL at the next point is greater than or equal to the threshold. After finding all candidates for a local maximum, all possible pairs of preceding and succeeding candidate points (one from the candidate points preceding the local maximum and the other from the candidate points succeeding the local maximum) are tested to identify the wave types (e.g., in Figure 3(b), two pairs A-B and A-C are tested). A wave between two candidate points is classified as a sharp if AALs of both points are smaller than the threshold ( AAL ), and their amplitudes are greater than or equal to the threshold ( amp ). The pattern is classified as a SW if the AAL of a candidate point is bigger than the threshold ( AAL ) and the amplitude is greater than or equal to the threshold ( amp ). In addition, any sharp or SW with length (duration) outside expected range of an ED is discarded since combining two candidates may indicate a waveform with longer duration than a sharp. Waveforms with unbalanced ascending and descending length are discarded as well. For example, a wave is discarded if bal amp < bal amp , where bal amp = small /( small + big ), denotes amplitude of ascending or descending slope, small denotes the smaller amplitude, and big denotes the bigger amplitude. We call bal amp the "balance in amplitude." The "balance in time," which is used for further processing, can be calculated in a similar manner: bal time = short /( short + long ), where short denotes the shorter time length and long denotes the bigger length. After testing for all possible pairs, all the overlapping ranges found for a local maximum were merged into a single range, which was then used as a feature for the automated classification of EDs.
In later experiments, we used different amplitude thresholds for the sharps and SWs because the amplitudes of sharps and SWs generally differ from each other. Moreover, the moving range of MP was restricted between [ , ], where both and are the distances from the local maximum under consideration.
The wave detection algorithm can be summarized as follows: (1) Find all the local maxima of the signal, which are candidates of a wave peak.
(2) For each local maximum point, the following steps are conducted: (A) For each direction (left or right), (i) set an MP at away from the local maximum. Let the number of MPs be denoted as , (ii) draw a line from the local maximum to the MP, (iii) calculate AUL, (iv) set MP −1 as a candidate of the wave-end if AUL −1 < AUL and AUL ≥ AUL , (v) move the MP a step outward, (vi) repeat steps (ii) to (v) until the distance of MP from the maximum reaches .
(B) Check for all possible pairs of candidates if they satisfy conditions of sharps or SWs and merge all the (overlapping) ranges of the satisfied pairs. (Figure 4). The wave detection algorithm in Section 2.3 was designed to detect waves whose deflections are downward. Since inverting signals would be inconvenient to satisfy this condition, we introduce a simple automatic procedure for a signal to be set upright (the deflection of waves is downward). This procedure is conducted for each data segment independently. First, the sharps are detected from both the original and inverted signals, and then we determined which signal is upright by evaluating the "dominance" of the detected sharps. The dominance of a signal is calculated using the following equation:

Procedure of the Proposed Method
where is the number of detected sharps and ascend and descend denote the lengths (duration) of ascending and descending blocks of the th sharp. The signal with the larger value is deemed to be the upright signal. Please note that the channel information of the segment is hidden to the wave detection algorithm and the sharps in a segment are detected when determining the dominance. Afterwards, the SWs were detected from the upright signal.
The parameters for the detection of sharps and SWs were derived from a sample dataset, which was randomly selected from all the data segments. Ten segments were selected from each type of ED and the ranges of EDs were labeled manually. Geometrical features such as AUL and AAL were calculated for each of the EDs, and the parameters were determined based on the mean and standard deviation to include most EDs in the sample dataset, as follows: where the subscript type represents either sharp or SW, ER type means the expected range of an ED, AUL/AAL type , amp type , bal type , and type denote AUL/AAL, amplitude, balance, and wave-width of an ED, respectively. The operator ⋅ denotes mean and denotes standard deviation. Please note adding or subtracting standard deviation depends on the types of thresholds. The standard deviation is added when waveforms with geometrical feature values bigger than a threshold are to be discarded and is subtracted if waveforms with the values smaller than the threshold are to be discarded. The derived parameters are shown in Table 1. The next step for identifying EDs is calculating geometrical and relational features from the source signal and   [14]. Some features were used as they were in the literature, while others were slightly modified. For example, an original feature "ratio of duration of a wave's increasing part to the amplitude" was changed to "ratio of duration of a wave to the amplitude." The new features, 2 and 14, were introduced to check the presence of suspected EDs and 26 and 27 were introduced to describe detailed characteristics of SSWs because sharps and SWs are often accompanied by each other.
Lastly, the support vector machine (SVM) was used for the classification of EDs. A Gaussian radial basis function was empirically chosen for a kernel. To classify three different types of segments using SVM, a two-step classification approach was employed (denoted by three-class classification in Results). This approach was used to distinguish between normal signals and EDs and subsequently between repeated sharps and SSWs. We also tested two-class classifications for each pair of segment types (denoted by pairwise classification in Results). We used a statistics and machine learning toolbox implemented in Matlab (ver. 10.1, MathWorks, USA). The parameters used for SVM are listed in Table 3.
The proposed method was evaluated using 10-fold crossvalidation. The segmented iEEG data were shuffled and divided into ten groups so that each group has the same number of EDs and normal segments. SVM was trained using data from nine groups and then the remaining data were tested. The training and testing were repeated ten times so that all the data were included. To avoid any failures in training caused by severely unbalanced data, only 20,000 Table 2: A full list of features used for the classification of ED segments: ALL denotes all data (source signal) in a segment; NS denotes data in the nondetected (as a sharp or SW) block only within the segment. Some features are paired (e.g., numbers 9 and 21), which can be calculated for sharps and SWs independently.   For validation of the proposed method, the wavelet transform introduced in Section 2.2 was also used for the classification with SVM. After the 10-fold validation was finished, the accuracy of each method was averaged over all validation iteration. The accuracy of a method was calculated as follows for each iteration: where denotes the number of target classes, Sens denotes the sensitivity of a target class, Spec denotes the specificity of the class, Sel denotes the selectivity [24] of the class, TP denotes the number of correctly classified segments for the class, and FN denotes the number of target segments that were not classified correctly. We used the mean of the sensitivities for all the classes instead of the general definition of the accuracy because the numbers of segments were severely unbalanced in our data. Statistical significance was tested using Wilcoxon's signed rank test. Table 4 lists the classification accuracies of the methods used for the distinction between segments with EDs and normal patterns. As shown in the table, the proposed method could stably detect the EDs from normal segments. The average accuracy of the proposed method was 93.76%, which is 2.51% higher (Bonferroni corrected = 0.006) than the conventional wavelet-based method (average accuracy = 91.25%). The accuracy was significantly decreased (91.96%) when the wavelet features were utilized together with the proposed method (Bonferroni corrected = 0.0117). One of the main reasons for this decrease could be that the number of features was increased dramatically to 67 features, whereas a decrease could have been overcome by reducing dimensions or selecting optimal feature sets [25]. The selectivity was increased when the wavelet features were utilized together with the proposed method. This is because that the sensitivity and specificity were unbalanced and the number of normal segments was relatively bigger than ED segments. Table 5 lists the accuracy of classification between repeated sharps and SSW segments. The accuracies were generally lower compared to the previous binary classification between EDs and normal segments. In particular, the sensitivities of the wavelet-based method were severely inconsistent over iteration, reporting a number of extremely low sensitivities. The mean sensitivity for the sharp segments barely reached 57.14%. Although the mean sensitivity and selectivity for the SSW segments exceeded 90%, it was caused by the low sensitivities of the sharps identification.

Binary (Pairwise) Classification.
In contrast to the wavelet-based method, our proposed method kept the sensitivities balanced in both classes. The mean sensitivities for the repeated sharp and SSW segments over iteration were 79.76% and 73.99%. Table 6 lists the accuracies of the classification of three different types of iEEG segments: normal, repeated sharps, and SSW. We used the same parameters as the previous pairwise tests. The overall accuracy of our proposed method was 87.69%, which was significantly greater than that of the wavelet-based method (81.97/86.96%, = 0.0176). The mean accuracy of the combined method (utilizing proposed and wavelet features together) was slightly lower (0.74% ) than that of the proposed method, although the difference was statistically insignificant. The sensitivities of each class trended similar to those of the previous pairwise Table 5: Classification accuracy between repeated sharps and SSW segments: the notations for the methods are the same as those used in Table 4. tests. Since the classifier for the three classes was composed of binary classifiers, the mean sensitivities for the normal class were similar to those in Table 4. The slight difference (less than 0.1% ) originated from the random selection of normal segments in each iteration number. The low accuracy of the wavelet-based method was mainly caused by a failure in identifying sharps. The sensitivity of the waveletbased method for identifying repeated sharps was 32.86 ± 16.14%, which was dramatically lower than the previous binary classification between repeated sharps and SSWs (57.14% ± 28.10%), implying that the detection of sharps was less successful than that of SSWs. In contrast, our proposed method demonstrated well balanced sensitivities in detecting repeated sharps and SSWs, as in the previous binary classification results (Table 5). Table 7 lists the distribution of outputs for each type of iEEG segments. Repeated sharp segments were misclassified as normal segments more frequently than the SSW segments when the wavelet-based method was used. In contrast, the confusion rates of repeated sharps and SSWs into normal segments were relatively similar to each other, when the proposed method was utilized solely or together with the wavelet features.

Three-Class Classification.
To evaluate the relative importance of the features used, we evaluated accuracies of each feature group ( Figure 5). The best feature group in this test was the number of detected waveforms (2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14), which resulted in a classification accuracy of 78.9%. The second and third feature groups (6-18 and 11-23) also showed high classification accuracies similar to the accuracy of the first feature group. The second-best group was the "mean ratio of amplitude to wave-width of detected waveforms" and the third-best group was the "difference between the maximum and minimum amplitude of undetected blocks," which showed 78.6% and 78.5% classification accuracies, respectively.
A limitation of this work would be that the proposed method was tested only for a single patient's data, which makes it difficult to claim that the proposed method is better than a wavelet-based method. Instead, we would like to emphasize that we developed a new method to classify two different types of EDs and background activities. Since the  proposed method achieved fairly high classification results with the parameters determined from 20 randomly selected segments only, we expect that the same approach would work for other patients' data if the shapes of EDs are stable.

Conclusion
In this paper, we have proposed a novel method to distinguish two types of EDs, repeated sharps and runs of SSWs, from normal iEEG background activity in a patient with LGS. Our proposed method was developed to detect sharps and SWs independently and to classify the categories based on the geometrical and relational features calculated from the detected sharps and SWs. The results showed that we can differentiate EDs with fairly high accuracy (87.69 ± 4.11%), which is greater than a conventional frequency-analysisbased approach (accuracy of 81.97 ± 2.24%). Our method can be used as an auxiliary tool for the surgical planning of intractable epilepsy such as LGS.

Appendix
This appendix shows the detailed procedure to find local maxima from a given discrete signal.

8
Computational and Mathematical Methods in Medicine  Let be a time-series signal and a value of the th data point of .
(1) Set a checkpoint at 2 and a variable bUP to be false, where bUP is a Boolean variable to represent if the preceding signal from the checkpoint satisfies a specific condition. (2) Test whether the checkpoint is a local maximum with the following steps. Let the data index of the check point be . (3) Move the checkpoint forward until it reaches the end of the signal.