Deep learning using electroencephalogram (EEG) data for diagnosing and predicting SSRI response in major depressive disorder

Response definition

Across all the sites, treatment response was defined as a > 50% decline in depressive symptoms as assessed with the rater-based questionnaire, either the Montgomery-Åsberg Depression Rating Scale (MADRS)¹⁶, or the Hamilton Depression Rating Scale (HDRS)¹⁷. Duration until assessment of response was defined within the studies of the corresponding datasets and varied between 4-8 weeks. Written informed consent to participate in the study was obtained from participants at all sites. Details are provided in Table 1. Ethical approval for this analysis was obtained from the faculty of medicine, University Zurich.

CANBIND dataset

As part of the Canadian Biomarker Integration Network in Depression (CAN-BIND-1, details provided elsewhere¹⁸), eyes closed (EC) resting state EEG data (8 minutes) from N = 309 patients with MDD and N = 146 HC recorded between 2013 and 2017 were accessed through a controlled data release (10.60955/mnwq-sq07). Participants were labelled with unique subject IDs, and all identifiable information was removed from the data files. All patients were treated with the SSRI escitalopram (10-20 mg/d) for 8-weeks; symptom severity was assessed at baseline and after 8 weeks using the MADRS. Data from N = 175 patients and 54 HC were used for EEG analysis. From initially 309 recorded datasets from patients at the CANBIND study, 212 sufficient EEG datasets from baseline were available for this analysis, with 37 datasets missing MADRS scores after 8 weeks, resulting in a total amount of 175 patient datasets to be included in the analysis for this study. This secondary data analysis study was approved by the Royal Ottawa Health Care Group Research Ethics Board.

Ottawa dataset

As part of a larger clinical trial¹⁹, 3-min of EC resting state EEG data from N = 51 unmedicated patients with MDD and N = 43 HC were recorded between 2008 and 2011. Patients were randomized (double-blind) to one of three antidepressant regimens: escitalopram + bupropion, escitalopram + placebo or bupropion + placebo. Only patients (N = 19) randomized to receive the SSRI escitalopram were included in the current analysis. Baseline symptom severity and outcome after 12 weeks were assessed using the MADRS. This study was approved by the Royal Ottawa Health Care Group Research Ethics Board.

Leipzig dataset fMRI

As part of the LIFE-Adult project Leipzig²⁰, Germany EC resting state EEG data from N = 31 patients with MDD and N = 32 HC were recorded between 2012 and 2014. Patients were treated with different medications, in total, N = 6 received SSRIs were used for the predictive analysis. Baseline symptom severity and outcome after 4 weeks were assessed using the HRRS. This study was approved by the University Leipzig Ethics Board.

Leipzig dataset prediction

As part of the Vigilance Algorithm Leipzig (VIGALL) project²¹, Germany, EC resting state EEG data from N = 22 patients with MDD and 17 HC were recorded between 2013 and 2014. Patients were treated with different medications, in total, N = 15 received SSRIs and were used for the predictive analysis. Baseline symptom severity and outcome after 4-8 weeks were assessed using the Hamilton Depression Rating Scale (HDRS). This study was approved by the University Leipzig Ethics Board.

Praha dataset 250 Hz

As part of the II-D-AD-QEEG project, Czech Republic, EC resting state EEG data from 67 patients with MDD were recorded between 2005 and 2011. Patients were treated with different medications, in total N = 21 received SSRIs and were used for the predictive analysis. Baseline symptom severity and outcome after 4-6 weeks were assessed using MADRS. This study was approved by the Ethics Board of Praha, Charles University.

Praha dataset 1000 Hz

As part of ther research projects carried out at NIMH Klecany, Czech Republic, resting state EC EEG data from 93 patients with MDD were recorded between 2015 and 2019. Patients were treated with different medications, in total N = 36 received SSRIs and were used for the predictive analysis. Baseline symptom severity and outcome after 4-6 weeks were assessed using the MADRS. This study was approved by the Ethics Board of Praha, Charles University.

Preprocessing pipeline and standardization

Export to edf-format was done using Brain Vision Analyzer 2.2 (Gilching, Germany). All closed-eyes resting state EEG files were visually screened using DeepPSY software (version 1.02, Zollikerberg, Switzerland). Additional calculations of EEG-features, filtering and artefact removal were done using MNE version 1.2.0. The overlap of EEG channels between all studies was identified (10 channels: ‘F7’, ‘F4’, ‘P3’, ‘O1’, ‘F3’, ‘C4’, ‘F8’, ‘O2’, ‘P4’, ‘C3’) and only these ten channels were kept for all datasets and analyses. A detailed description of the electrode positions can be found in the supplement (Supplementary 1.1). Further, the lowest sampling rate was chosen (250 Hz) to standardize across all the EEG datasets and a low-pass filter (100 Hz) was applied before downsampling to avoid aliasing of frequency components. All EEG files then were filtered between 0.5 (high-pass cut-off for low frequency noise) and 45 Hz (low-pass cut-off to exclude muscle artefacts usually starting at 30 to 60Hz²²) and segmented into 2 second segments after generating an average reference. We used a 0.5–45 Hz band-pass filter as a compromise to reduce low-frequency noise while retaining potential gamma-range EEG activity ( > 30 Hz), acknowledging in the limitations that residual muscle artifacts may persist in this range. (Additional analysis on a 0.5-30 Hz band pass filter can be found in the Supplementary 1.11). This resulted in a different number of segments for the patients and control subjects of the different datasets, ranging from 36 segments per subject to 450 segments. Additionally, we applied an eye-artefact correction method²³ based on electrooculogram EOG) data or reconstructed EOG channels (F7-F8) and a bad segment rejection method based on the min-max amplitude criterion ( ± 100 μV threshold) for deep learning and machine learning analysis. Further, we conducted a comparison of the amount of eye-related events/minute for all datasets (details in Supplementary 1.12). Details regarding the different datasets is outlined in Table 1. To achieve a higher generalization of the results for future research, a more detailed description of the recording settings from the different sites can be found in the supplement (Supplementary 1.7). Due to the strongly unbalanced number of EEGs across sites, enforcing uniform representation of all datasets in each batch would require extensive oversampling of smaller datasets, thereby skewing site weights and increasing the risk of overfitting rather than improving interpretability; thus, site-specific accuracies are provided in the supplement (Supplementary 1.8) instead.

Deep learning models

Using a Keras Backend for TensorFlow (version 2.4.1) with Python 3.8 on a RTC 3090 GPU, the model used in previous work^24,25 was implemented for both differentiating (i.e., discrimination of HC and people with MDD) and predictive tasks (i.e., responder vs. non-responder). The Convolutional Neural Net (CNN) included six convolutional layers with 200 filters, employing kernel sizes from (2, 2) to (1, 2), followed by max pooling and dropout layers, with a dropout rate of 0.3 after the first five convolutional blocks and 0.5 after the sixth. Hyperparameters were adjusted using a random search approach. A hyperparameter grid, including the number of filters, dense layer size, dropout rate, and learning rate, was defined (see Supplementary 1.6). Using RandomizedSearchCV, we trained multiple configurations on the training data, evaluated them using cross-validation, and selected the best-performing model. The final model was trained on the optimized hyperparameters and evaluated on the test set.

The network was configured to receive inputs of shape (10, 500, 1) and outputs a binary classification. To minimize inter-dataset variability, all EEG segments were amplitude-normalized by subtracting the mean and dividing by the standard deviation of the training set prior to model training, with the same scaling applied to the test data. Training was conducted over 900 epochs with an Adamax optimizer set to a learning rate of 5 × 10⁻⁵, and early stopping was implemented with a patience of 70 on validation loss to prevent overfitting with a total possible maximum duration of 700 epochs. Weight change was adjusted to fit the unbalanced datasets, resulting in a weight of 2:3 for HC versus MDD and 17:20 for Responders versus Non-Responders. A listing of all convolutional network layers can be found in Fig. 1.

Statistics and reproducibility: training, validating and testing

The EEG files suitable for further analysis from all datasets (CANBIND; Prague I/II, Leipzig I/II and Canada) were separated into a training set, a validation set and a test set. “Files” here refers to the complete recording of one subject while “segments” refer to a 2-second segment of on subject. The test set consisted only of EEG files from patients or controls that had not been used for training or validation to make sure the network did not learn individual features. This is a very important step of the testing procedure since it guarantees the generalizability of the results toward unseen patients and controls. The test set was randomly chosen with 40 EEG files set aside for classification test set MDD versus HC (balanced 20 files each group) and 20 files set aside for classification test set Responder versus Non-Responder (balanced 10 files each group). The files were randomly distributed across the recording sites.

Within the remaining corpus of EEG files, 5% of the segmented files (all files were split into 2 s segments) were used for validation using a random shuffle and splitting routine. To prevent shuffling segments from one subject into the training and testing sets, all segments received identifiers according to their subjects of origin. An early stopping of leaning epochs was implemented following a failed decrease of the validation loss of the optimizer function for consecutive 70 epochs. The maximum of epochs for training was set to 900, but never was used since validation loss criteria for early stopping were fulfilled earlier in all trials. It is important to note, due to an overlap of terminology in EEG research and Deep Learning research, that, in this paper, an epoch is not an EEG segment, but a run of the convolutional model through all training data. A short period of timeseries data from the whole EEG is called a segment instead (i.e., 2 s segment).

It is also important to note that the main outcome was set to the correct labelling of subjects, not single segments. That means, a subject in the test set was counted as correctly classified if 50% or more of single segments were classified according to the appropriate label of the subject (i.e., patient or control; responder or non-responder). Thus, there were different accuracies for single segment classification and for subject-wise classification. Although this approach could over- or underestimate the accuracy of the model at the segment-level, it seems to be the more meaningful approach since it is clinically relevant to identify subjects (e.g. responders), and not an EEG-segment.

Grad cam and important EEG topographies and frequencies

The method employs Gradient-weighted Class Activation Mapping (Grad-CAM) to visualize areas of interest in the 2D EEG arrays that influence the predictions of the CNN. By averaging the activation maps from convolutional layers, it generates heatmaps that highlight the significant regions for the model’s decisions. The code snippet iterates over a set of images, computing Grad-CAM heatmaps for each, then overlays these averaged heatmaps onto the original images for visualization. Finally, it computes and displays an overall averaged heatmap across a subset of images to identify common patterns of activation. This approach is instrumental in interpreting the model’s behavior by revealing which features contribute most to its predictions. We binarized normalized Grad-CAM maps at τ = 0.3 to generate seeds, following common practice in WSSS (Weakly Supervised Semantic Segmentation)/CAM pipelines, and prior imaging work that reports/visualizes Grad-CAM at 0.3²⁶. We then calculate a Fourier spectrum for every marked sequence from every two-second segment. Thus, a power spectrum is obtained for every channel that highlights the EEG frequencies with the largest discriminative power. Further, the visualisation of the weighted Grad-CAM outputs allows for the comparison of the diagnostic and predictive tasks by means of topography.

Conventional machine learning algorithms

For comparison of the differentiation between MDD and HC as well as for the differentiation between Responders and Non-Responders to SSRI treatment, several conventional EEG-parameters that have been found to be important for diagnostic or predictive purposes in MDD^11,27,28 were computed for all subjects after artefact—and eye movement removal. This included EEG-power of the Delta (0.5-4 Hz), Theta (4-8 Hz), Alpha (8-12 Hz), Beta (12-20 Hz) and Gamma (20-45 Hz) frequency bins at all 10 electrodes (‘F7’, ‘F4’, ‘P3’, ‘O1’, ‘F3’, ‘C4’, ‘F8’, ‘O2’, ‘P4’, ‘C3’ used for the combined analysis, Alpha peak frequency and averaged lagged linear connectivity measures between frontal, parietal, central and occipital lobes (total 6 variables), resulting in a total of 47 variables. The EEG-variable- matrices were then used for a machine learning (ML) approach using ten different established ML algorithms: In the past, EEG-based classification utilized various machine learning algorithms, each with distinct strengths. Linear models like Logistic Regression and Naïve Bayes offer interpretable decision-making but may struggle with complex patterns. Tree-based methods such as Decision Trees, Random Forest, Gradient Boosting, LightGBM, and XGBoost enhance accuracy through ensembling, with Random Forest and XGBoost being particularly robust for EEG feature differentiation. Distance and margin-based classifiers, including K-Nearest Neighbors (KNN) and Support Vector Machines (SVM), are commonly used for tasks like motor imagery classification and seizure detection. Boosting techniques like AdaBoost and Gradient Boosting iteratively improve weak classifiers but may be sensitive to noise. Thus, we implemented these ten algorithms with a “leave-one-out” cross-validation approach to compare the results to the deep-learning approach.

Usage of other deep learning architectures (EEGNet)

To compare the presented results with the performance with other state of the art deep learning networks, we used the EEGnet design to compute a model based on the same input data for the comparison between Responders and Non-Responders. The structure of the model was derived from the original work²⁹. EEGNet is a lightweight deep convolutional neural network (CNN) designed specifically for EEG signal classification tasks such as brain-computer interfaces (BCIs), motor imagery, and seizure detection. It employs depthwise and separable convolutions to efficiently learn both spatial and temporal EEG features, reducing computational complexity while maintaining high performance. The network typically consists of three main layers: (1) a temporal convolutional layer (2) depthwise convolution to capture spatial patterns across EEG channels, and (3) a separable convolution. Details can be found in the supplement (Supplemenary 1.4).

MDD versus HC

Classification of MDD subjects versus HC subjects suitable for further analysis from all datasets (CANBIND; Prague I/II, Leipzig I/II and Canada) showed a train accuracy of 73.16% and a validation accuracy of 71.13% (see Fig. 2, top right). Test accuracy for all pooled segments of unseen subjects from the test set was 64.84% (SD 2.0%) in total (not subject-wise) and the percentage of correctly classified unseen subjects based on the number of correct classified segments of a subject was 67.5% (average for 5-fold validation with best fold 70%, SD 1.8%) with a standard deviation of 30.4% across subjects. The confusion matrix yielded a sensitivity of 68.4% and a specificity of 66.7% (see Fig. 3, left).

When using an eye-artefact reduction method²³ and exclusion of artefact segment using a min-max approach, results were worse than without extensive preprocessing (correctly classified unseen subjects 65%, SD 28.4%). No differences in eye-related movements were found between groups (Supplementary 1.12). A slightly decreased accuracy was also found when applying a 30 Hz low-pass filter instead of 45 Hz low pass filter (see Supplementary 1.11) with 65%.

MDD SSRI prediction

Classification of responders versus non-responders of SSRI treatment showed an average of 71.8% accuracy for single segments for the training set after 732 training epochs. Validation set accuracy was 66.4%. The test set in these trials reached an average of 64% (SD 4.2%) for pooled segments and a mean value of 79% (5-fold average, best fold 85%, SD 2.5%) for correct subject-wise classification. The confusion matrix yields a sensitivity of 75% and a specificity of 87.5% with a standard deviation of 24.7% across subjects (see Fig. 3, right).

Again, when using extensive preprocessing procedures, the results were found inferior with correctly classified unseen subjects 70%, SD 21.0%). Regarding the fact that different times until response were used in the different datasets, we decided to do a sensitivity analysis. Thus, the accuracies for correct labelling were compared for responders and non-responders in dependence of the timepoint of the MADRS or HDRS assessment. There was no significant difference between sites with assessments from weeks eight to twelve (Canada and CANBIND) and other sites with assessments between weeks four to eight (t-statistic: −0.84, p > .41). A decreased accuracy of 55% was also found when applying a 30 Hz low-pass filter instead of 45 Hz low pass filter (see Supplementary 1.11). No differences in eye-related movements were found between responders and non-responders (Supplementary 1.12).

Machine learning classification and EEGNet results

The results of the ten different machine learning results for the differentiation of patients with MDD and HC varied between 52.1% (Naïve Bayes) to 60.6% (XGBoost) and 66.1% (Random Forest). Thus, the results of the tree-based algorithms came close to the results of the Deep Learning approach but did not reach the same accuracy. Regarding the differentiation between Responders and Non-Responders, the results of the ML approaches stayed way behind the accuracy of the Deep Learning approach with some algorithms staying below random-guess accuracy (e.g. Gradient Boosting with 47.1%) to a maximum of 60.3% with Ada Boost. More details can be found in the supplement (Supplementary 1.3).

Although the computing time was much lower in comparison to the other DL approach, the results of the EEGNet (more details on the model can be found in the Supplementary 1.4 and 1.5) for differentiation between Responders and Non-Responders stayed well below the already reported results with training and validation accuracies reaching 74% each, but achieving finally only 55% of correctly classified subjects (SD 31.6%) in the test set. Among the tested architectures, the final CNN achieved the best performance due to optimized hyperparameters, appropriate regularization, and its capacity to capture the complexity of EEG patterns, rather than model depth alone.

Grad-CAM, EEG topographies and frequencies

Using the Grad-CAM technique, the layer activations for the subjects in the unseen test set were identified at the segment-wise subject level as well as the averaged global level over all subjects. This allowed us to map the topography of EEG-channel importance for the label classification (Figs. 2 and 4, bottom). The resulting maps for the HC versus MDD classification showed (Fig. 2) a clear accentuation of the electrodes at F7, F3, P3 and C4 while electrodes O1, O2, F3, F8, P4 and C3 showed less involvement in the classification process. Results of mapping these results onto the scalp-electrode positions, the relevant topography for EEG activity that divides HC subjects from patients with MDD, can be seen in Fig. 2, bottom right.

Using the same technique for the differentiation between Responders and Non-Responders, it became apparent that only electrodes at F3 and F7 showed a high weight in correct classification of patients. The mapping results can be seen in Fig. 4, bottom right.

Following this, we tried to determine the frequencies of EEG activity that contributed most to the correct classification. We therefore computed heatmap overlays for all classified EEG segments (example Fig. 5). Then, we extracted the periods with a weight >0.3 for all channels and calculated the EEG power spectra (see Fig. 6). The most prominent peak for all relevant channels can be found between 8-12 Hz, resembling the EEG alpha rhythm.

Number-needed-to-treat

Using the confusion matrix for treatment response with an accuracy of 80%, we identified that 9 patients are true positives and 7 are true negatives. When the marker is applied, responders (9 patients) receive SSRI treatment, achieving a 100% response rate. Non-responders (7 true negatives and 3 false negatives) receive SNRI treatment, which has a 50% response rate, resulting in 5 responders. In total, 14 out of 20 patients respond to treatment, yielding an overall response rate of 70%. This represents a 20% improvement compared to the baseline response rate of 50% for SSRI treatment without the marker. Table 2 depicts the increase of responders and decrease of numbers needed to treat for a rising accuracy of correct classification. Under the assumption that SSRIs and SNRIs as alternative show a 50% response rate in an unselected patient cohort³⁰, an accuracy of 80% in the prediction of SSRI response using the EEG has a number-needed-to-treat (NNT) of five. We added the receiver-operator-characteristics and the area-under-the-curve values in the supplement (Supplementary 1.10). To make the results comparable to other studies, we computed also the normalized positive predictive value (nPPV) as:

Peer review information

Communications Medicine thanks Danish M. Khan, Shiyang Lv and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. A peer review file is available.