2.1 Epileptic EEG database
The CHB-MIT Scalp EEG Database (URL: https://physionet.org/content/chbmit/1.0.0/) [13] was selected for this study because it was annotated by clinical investigators, broadly adopted, and tested in a large number of studies [14–17]. It contained 24 cases of EEG data from 23 epileptic patients (the first and the 21st cases were from the same patient recorded one and a half years apart). The patients ranged from 3 to 22 years old. Most of the data contained 18 or 23 channels and a few 24–26 channels. The data recording followed the international 10–20 system for EEG electrode positions. The sampling rate was 256 Hz and the digitization resolution was 16 bits. The number of seizure events for the 24 cases ranged from 3 to 40, and the shortest seizure events lasted 6 s. The total seizure time was 11,611 s, and the total non-seizure time 680,705s. Table 1 summarized the basic statistics of the database. We analyzed the channels common to all the cases for the power spectral analysis and then selected the one channel based on its power spectral characteristics most relevant to the seizure events to test our approach to seizure detection. These channels include FP1-F7, F7-T7, T7-P7, P7-O1, FP1-F3, F3-C3, C3-P3, P3-O1, FP2-F4, F4-C4, C4-P4, P4-O2, FP2-F8, F8-T8, T8-P8, P8-O2, FZ-CZ, CZ-PZ, T7-FT9, FT9-FT10, FT10-T8, for a total of 21 channels. Please note that channels T7-P7 and P7-T7 are of the same pair of electrodes but with different polarities. Since they give the same power spectral density (PSD) function, from which the features are derived for this study, we only need to select one
of them.
[Table 1 goes about here]
2.2 EEG preprocessing
The EEG data were filtered using a zero-phase high-pass filter (-3dB cutoff frequency at 0.5 Hz) and a low-pass filter (-3dB cutoff frequency at 45 Hz). Then the data was partitioned into ictal and inter-ictal segments according to the expert annotations provided with the database. As literature suggested that inter-ictal data may be further divided into pre-ictal and post-ictal periods although their specific durations of them were not specified [18]. In addition, the inter-ictal period was much longer than the ictal period, leading to an imbalanced data problem for machine learning. We thus followed previous studies [19, 20] to define pre-ictal and post-ictal periods respectively as the 10 minutes before and after the ictal state and used the rest parts as the inter-ictal period. The data during the inter-ictal and ictal periods were then segmented into 5-second epochs and their direct current component removed. This choice of epoch length was considered to be a compromise between the ability to reliably capture the EEG features and the signal stationarity assumptions [21].
2.3 Power spectral density estimation and parameterization
The PSD of each epoch was estimated using Welch’s methods [22], with a window size of 2 seconds (512 data points) with 50% overlapping. Without zero padding, the algorithm gave a frequency resolution of 0.5 Hz.
A striking property of EEG signals revealed in previous studies was that their PSDs followed a 1/f-like distribution [23]. The recent development of PSD analysis highlighted the value of decomposing the PSD into aperiodic and periodic components, among which the aperiodic components refer to the overall 1/f profile of the PSD, and the periodic components are the bumps superimposed on the aperiodic components. Therefore, the periodic components are the difference between the total PSD and the aperiodic components. The aperiodic and periodic components are respectively modeled as a Lorentzian function and multi-Gaussians [12]. The Lorentzian function for the aperiodic components is formulated as:
$$L=b-\text{l}\text{o}\text{g}\left(k+{F}^{\chi }\right)$$
where \(b\) is the broadband offset, \(\chi\) the exponent, \(F\) as the vector of input frequencies, and \(k\) the ‘knee’ parameter that controls the bend in the aperiodic component. The parameter k determines two modes of fitting: (1) in the fixed mode, k is set to 0; (2) in the knee mode, k is to be identified with the data. While the fixed mode intends to fit the aperiodic component using a straight line in the log-log space, the knee mode is recommended when the data extends a broad frequency range. The aperiodic components can thus be fully characterized by only two parameters, offset b and exponent χ in the fixed mode, with an additional parameter k in the knee mode.
The multi-Gaussians for the periodic components is formulated as the combination of N Gaussian functions:
$$G=\sum _{i=0}^{N} {G}_{i}$$
Each Gaussian \({G}_{i}\) was a curve fit to a periodic component in the power spectrum:
$${G}_{i}={a}_{i}\times \text{e}\text{x}\text{p}\left(\frac{-(F-{{c}_{i})}^{2}}{2{{w}_{i}}^{2}}\right)$$
where \({a}_{i}\) was the power of the peak, in log10(power) values; \({c}_{i}\) and \({w}_{i}\) were the central frequency and the standard deviation of the Gaussian \({G}_{i}\); and \(F\) was the vector of input frequencies.
The total PSD function was thus finally parameterized as the sum of the aperiodic components \(L\) and the periodic components \(G\):
As most of the EEG energy resides within 0.5 ~ 30 Hz and the epileptic discharge waveforms have relatively low frequencies [11], we chose to parameterize its PSD in the frequency range of 0.5 ~ 30 Hz. Figure 1 depicts the procedure and typical results for the analysis. Figure 1A shows a one-hour EEG recording of channel 18 from patient case 11. We first segmented the inter-ictal and ictal EEG into epochs of 5 seconds. The typical inter-ictal and ictal epochs are shown respectively in the left and right panels in Fig. 1B. Their corresponding PSDs estimated with Welch’s method are plotted in black solid lines respectively in the left and right panels in Fig. 1C, within which the full models of the PSDs (solid red lines), the aperiodic components (blue dashed lines), and the detected peaks of the periodic components (green triangles) are also presented. Figure 1D and E show respectively the models for the aperiodic and periodic components. The key parameters of the models are summarized at the end of Fig. 1.
The above procedure was applied to all the 21 channels (discussed in section 2.1) for each case and the parameters were analyzed for their relationship to the seizure events. Some of the parameters were then selected and tested for their utility for seizure detection.
[Fig. 1 goes about here]
2.4 Statistical analysis
We performed statistical analysis on the PSD parameters about their relevance to the ictal and inter-ictal states channel by channel and then case by case. The distribution normality of the parameters was tested first. If the parameter followed a normal distribution, its descriptive statistics mean ± standard deviation (std) were computed, or else the median and the interquartile range (IQR) were computed. For parameters meeting the normality criteria, an independent-samples t-test was performed; for parameters not following normal distribution, the Kolmogorov-Smirnov test was used. P-value less than 0.05 for a two-sided test was considered statistically significant.
2.5 The selection of EEG channel
For each patient, we selected one channel of EEG to represent his/her inter-ictal and ictal states based on the P-value ranking for the four key parameters (offset, exponent, and the first and second highest peaks in PSD). The choice of the above four was based on the systematic statistical analysis of all relevant candidate parameters as the output of the PSD parameterization model (see Results for details). Taking case 11 as an example, we systematically examined the parameters offset and exponent for the aperiodic components, and parameters PW, CF, and BW for each periodic component and screened them down to six parameters (Tables 2 ~ 4) for their relevance to inter-ictal and ictal states. After comparing the results in Tables 3 and 4, only the parameters in Table 4 were maintained.
We first ranked the channels by the P-values for each of the four parameters in the ascending order - the channel with the lowest P-value was ranked first (Table 5A). According to the current results, only the channels ranking in the top 8 were necessary to be considered for all cases in this study. We then sorted out the channels with all four parameters included in the above tables and scored the channels by their rankings for each parameter. The total scores (Rtotal) were calculated and ranked in ascending order (RCH) (Table 5B). The channel ranked first was finally selected as the representative EEG data for that patient case. If there was no channel with four parameters ranking in the top 8, the channels with three parameters in the top 8 were ranked as above (the second half in Table 5B).
[Table 2 goes about here]
2.6 Seizure detection
To evaluate the utility of the parameters obtained from PSD parameterization, the four parameters extracted from the EEG signal from the channel selected using the above procedure were used for automatic seizure detection with support vector machine (SVM). In addition, to demonstrate the potential clinical application of our methods, continuous recording EEG were taken as input to the classification model. That is, except the ictal data, the rest of the data were regarded as inter-ictal.
Due to the severe imbalanced data problem shown in Table 1, the ratio between the seizure- and the non-seizure-time for most patients was less than 4%, synthetic minority oversampling technique (SMOTE) [24] was used for data augmentation to alleviate the bias toward the majority class from the SVM classifier. The SMOTE algorithm was implemented through an iterative procedure. For each minority sample, its k-nearest neighbors were firstly identified, among which N neighbors were randomly chosen. A synthetic sample was generated by adding the minority sample with a fraction
of its difference to each chosen neighbor. The range of fractions is between 0 and 1. In this study, k was set to 5, and N was set adaptively to the integer ratio between the number of the majority and minority samples. If N was greater than k, then N = k.
The support vector machine was employed to classify the ictal and inter-ictal states. The classifier options were set with polynomial kernel function, two polynomial orders, and standardization. Five-fold cross-validation was used to train and validate the classifier.