Application of an EEG-based deep learning model to discriminate children with epileptic spasms from normal controls

Standard protocol approvals

This study was approved by the institutional review board at the University of California, Los Angeles. Given the retrospective design and deidentification of EEG data, the requirement for informed consent was waived.

Study design and subjects

This is a retrospective cohort study, utilizing the same subjects and EEG samples that were used in our prior studies evaluating candidate biomarkers. Our methods for subject selection are described in detail therein.^10,11 Briefly, our cohort consisted of 50 patients with epileptic spasms (cases) and 50 normal controls. All cases underwent overnight video-EEG which demonstrated epileptic spasms. The normal controls all underwent overnight video-EEG for evaluation of suspected seizures, but all normal controls were found to have events that were not seizures (e.g., gastoesophageal reflux, normal movements, etc.) and were deemed neurologically normal on a clinical basis.

Data acquisition

Our methods for EEG data collection have been described previously.^10,11 In brief, all clinical data were obtained from the electronic medical record. For each study, four EEG clips were abstracted (two in wakefulness and two in sleep, approximately 10 each). As such, a total of 400 EEG samples were used in the analyses described here. Importantly, the specific epochs selected for clipping were dictated by a randomization algorithm such that the abstractor selected EEG data beginning at a specific randomly selected time.

EEG preprocessing and feature extraction

The overall methodologic framework is illustrated in Figure 1. EEG data from electrodes labeled Fp1, Fp2, F3, F4, C3, C4, P3, P4, O1, O2, F7, F8, T3, T4, T5, T6, FZ, CZ, and PZ were used. After applying a bandpass [1 Hz, 50 Hz] filter, EEG samples were segmented into windows of 15-seconds width. We applied the Fast Fourier Transform (FFT) to each channel to create a 2D channel-frequency image for each window. The input of the deep learning model is the 2D channel-frequency image of size 19 × 15000 (channels x frequency bins), with normalization to a standard gray scale image with pixel values ranging from 0 to 1.

• Download figure
• Open in new tab

Figure 1. Methodologic Framework

The workflow consists of three parts. In part (i), each EEG recording is segmented into consecutive 15- second windows, which window then converted into a normalized FFT map. In part (ii), the FFT maps are used in 5-fold model training. In part (iii), the model inference procedure is illustrated, whereby non-overlapping 15-second segments are systematically sampled from each test recording and entered into the trained neural network. Subsequently, the results for each 15-second segment are consolidated to generate a classification outcome for the patient.

Deep learning model architecture

To develop our model, we utilized the convolutional neural network (CNN), a state-of-the-art deep learning architecture that has demonstrated superior performance in the analysis of imaging¹² and EEG data.⁹ To mitigate common gradient vanishing and facilitate superior model convergence, we used a 34- layer CNN with residual connections block. The output of the model is the probability of each EEG recording having been derived from a case.

Model training and evaluation

Because ground truth labels (case/control) are missing for each 15s window of the EEG segment despite the availability of subject information, we adopted the weakly supervised learning approach, i.e. all windowed samples from case patients were labeled as case and control as control. We utilized a 5-fold cross-validation in which subjects were divided into five mutually exclusive groups, with 20 subjects (10 cases and 10 controls) in each group. For each fold, 80 subjects were used for model training and 20 subjects were used as a test set. Among the 80 training subjects, 10 patients (5 cases and 5 controls) were used for validation purposes to assess training performance in each training iteration. Overall, this training/evaluation approach supports the practical application of the model to novel subjects. Since the model is required to predict whether each 15-second EEG window is derived from a case or control, we used binary cross-entropy as the loss function, L=−[y·log(x)+(1−y)·log(1−x), where x is the probability of the windowed sample from the network, and the y is the label indicating case (positive, y = 1), or control (negative, y = 0). Given that EEG samples were of varying lengths, we adopted a novel stratified training approach which is widely used in the processing time-series data in the multimedia community.¹³ In each training iteration, for each patient, we randomly selected starting time stamps and constructed FFT maps from 15 second windows, such that information contributed from each patient was equal in each training iteration, with positive and negative samples balanced during the training phase. We conducted 400 training iterations in each fold. To improve generalization, we selected the model corresponding to the global minima in the validation set over training iterations, i.e., the highest balanced accuracy score.

As it is impractical to exhaustively evaluate model prediction on the FFT maps starting at every data point in the recording, we used a non-overlapping construction, such that FFT maps were generated from consecutive non-overlapping 15-second windows. To evaluate model performance we calculated accuracy, specificity, sensitivity, and F1 score. Finally, in order to classify a subject as case or control, we used the number of median-filtered successive case predictions in the non-overlapping construction, based on model confidence scores (predicted probability of each non-overlapping window). As such, a subject was predicted to be a case if any subset of windowed samples generated a prediction of case status. To protect external validity (i.e., generalization of the model to novel cohorts), we again employed a 5-fold cross-validation approach in the downstream analyses. Due to the sequential nature of the model inference (based on successive case predictions on the non-overlapping construction), we defined the earliest detection latency as the time from the beginning of each recording to the time when the model had sufficient information to classify the entire sample as having been derived from a case patient.

Occlusion sensitivity analysis and salient band discovery

We next employed occlusion sensitivity analysis, a contemporary technique used in the interpretation of black-box classifiers such as CNN, in the attempt to identify the frequency band(s) in which critical EEG features may reside. Specifically, for each input column (frequency band) of an input image x, we evaluated the model with that frequency band missing and observed how the output (prediction) changed. In particular, a column (frequency band) is deemed important when the change of the model confidence ||f(x) - f(x_without_i)|| is large. The occlusion analyses were carried out on sleep and awake EEG models separately, using data samples from the non-overlapping construction. An occluded rectangular region (occlusion map) was defined as 19 channels x 5 Hz with a stride size of 1/15 Hz. We then generated a patient-wise occlusion template using the median value across each pixel in all patient-specific occlusion maps, and finally generated an overall sleep occlusion template across subjects by using the median value of each pixel across all patient-specific occlusion templates (n = 50). We repeated the same procedure for awake EEG to generate an overall awake occlusion template. The salient frequency band which contributes the most to the case prediction of the model was then identified using adaptive thresholding (Otsu’s method) in the overall occlusion template map.

Targeted input perturbation

As illustrated in Figure 2, after salient frequency band identification by occlusion sensitivity analysis, the model was reassessed with targeted input perturbation, similar to methods previously employed by our group.^9,14 As above, this procedure was conducted separately for awake and sleep data. For each 15 second window derived from Case EEG, a band-reject filter was applied to remove data within the salient frequency band, and these data were then replaced with average control data (average normalized waveform over the non-overlapping construction of all control patients) within the same frequency band. The effect of this perturbation across all case subjects was visually evaluated with histogram analysis of case/control predictions. Statistical assessment was accomplished with a one- tailed t-test performed on the change of output probability scores from the deep neural network.

• Download figure
• Open in new tab

Figure 2. Targeted input perturbation workflow

The comprehensive workflow for implementing targeted input perturbation involves several steps. Initially, a specific 15-second segment of an EEG recording in waveform is chosen as the subject for perturbation. Frequencies other than the “biomarker” frequencies, which were determined from a previous occlusion sensitivity analysis, are retained by employing a band-reject filter. Following this, the so-called average control waveform is calculated by computing the numerical average waveform from all 15-second segment waveforms in control patients’ recordings. This average control waveform is then subjected to a band-pass filter to emulate normal control behavior within the biomarker band. The processed case waveform and average control waveform are subsequently combined to determine whether the case biomarker, as characterized by the processed case waveform, can transform the case waveform into a control. Finally, this synthesized waveform is input into the inference pipeline, and the resulting prediction from the deep neural network (DNN) is analyzed.

Statistical Methods

All statistical calculations were accomplished with Python (v3.7.3; Python Software Foundation, USA). PyTorch (version 1.6.0; Facebook’s AI Research lab) was used to develop the deep neural network. Continuous variables are described by median (interquartile range) or mean (standard deviation). Comparisons between groups were performed using chi-square, t-test, or Wilxocon rank-sum as appropriate. Significant results were considered with p < 0.05 unless stated otherwise.

Model performance

The trained model performed satisfactorily with 5-fold cross validation. Using sleep and awake EEG data, respectively, accuracy in classifying cases and controls was 0.95 (SD 0.049) and 0.91 (0.041), recall was 0.96 (0.054) and 0.84 (0.11), precision (sensitivity) was 0.94 (0.052) and 0.98 (0.040), and F1 score was 0.95 (0.049) and 0.90 (0.052). In the final model, we enacted a threshold classifier such that a full-length EEG sample was predicted to have been derived from a case (patient) if any 6 consecutive 15 second windows were classified as “case”. With EEG preprocessing, the minimum sample required for case versus control prediction (i.e., detection latency) is 120 seconds. In patient-wise 5-fold cross-validation, the median detection latency was 120 seconds (95% CI 120 – 135).

Interpretability analysis

In the occlusion sensitivity analysis, with the adaptive thresholding procedure, we found that model prediction was most affected by occlusion of data within the 9.7 – 22.0 Hz frequency band in sleep EEG, and the 1.0 – 8.6 Hz frequency band in awake EEG. After targeted input perturbation, in which these salient frequency bands were replaced with average control data, we found that patients who had been correctly classified as cases were now consistently re-classified as control subjects. The results of the targeted input perturbation are summarized in Figure 3.

• Download figure
• Open in new tab

Figure 3. Targeted Input Perturbation

The top (a – c) and bottom (d – f) rows refer to sleep and awake datasets, respectively. Panels (a) and (d) illustrate sampled FFT maps from correctly classified cases. Panels (b) and (e) illustrate FFT maps after perturbation. Panels (c) and (f) illustrate the absolute shift in case prediction with perturbation.