Machine learning on interictal intracranial EEG predicts surgical outcome in drug resistant epilepsy

For patients with focal drug resistant epilepsy (DRE), neurosurgery is the best available treatment and achieves seizure freedom in up to 60% of cases¹. The neurosurgical outcome depends upon precise and accurate delineation of the epileptogenic zone (EZ), the brain area that is indispensable for the generation of seizures^1,2. Yet, the EZ is a theoretical concept since its elucidation before resection is uncertain³. Even after a successful resection, the EZ cannot be precisely defined since the resected area can also be larger than its actual extent⁴. Therefore, there is currently no method of delineating the EZ directly; instead, this area is defined indirectly through several diagnostic tests whose results are often inconclusive or nonconcordant⁵. The seizure onset zone (SOZ), the brain area where seizures initiate, currently serves as the best approximator of the EZ¹. This area is defined by analyzing ictal intervals of intracranial electroencephalography (iEEG) data¹. Yet, removal of the SOZ does not always predict surgical outcome^6,7. Moreover, seizures may take several days (or even weeks) to occur, extending the duration of invasive recording at the expense of considerable human and financial resources⁸. An alternate interictal biomarker of the EZ is therefore of paramount importance⁹.

Several electrophysiological abnormalities are observed during interictal intervals, where background brain activity is altered by transient interictal epileptiform discharges (IEDs), such as spikes and sharp waves, and ripples^10,11. These abnormalities currently serve as established interictal biomarkers of epilepsy¹². However, they lack specificity in delineating the EZ since they can include areas larger than the EZ, which should be preserved during surgery¹³. Moreover, IEDs and ripples are identified through visual inspection of iEEG data or with the aid of automated software¹⁴. Yet, both approaches are prone to errors and biases¹⁵. Previous studies have shown that both IEDs and ripples do not occur only as isolated events in iEEG; in contrast, they often exhibit coherent or co-activated patterns^16,17,18,19, or propagate across contacts^20,21,22,23 forming brain networks (i.e., interconnected regions characterized by covariations and correlations among iEEG time-series)^{16,24,25,26,27}. Moreover, they occur transiently, alternating with background activity. Thus, the characterization of the brain network through the analysis of iEEG data segments, which contain both epileptiform as well as background activity, may be problematic for the accurate delineation of the EZ²⁸. Furthermore, IEDs are absent in ~10% of patients with DRE²⁹. Hence, the epileptologists cannot visually pinpoint all activity related to the disease in all patients. Thus, there is an ultimate need for a method that can capture both prominent and subtle epileptiform signatures and automatically distinguish the interictal epileptogenic network (IEN) from background activity³⁰.

Numerous artificial intelligence (AI) tools have recently been proposed to facilitate the presurgical evaluation of DRE patients^31,32. Supervised machine learning (ML) models have been used to extract epileptiform patterns from multimodal neuroimaging data and automatically delineate the EZ from interictal iEEG data using morphological, spectral, and temporal waveform features^30,33,34. Yet, automated delineation of the EZ through supervised ML models is challenging since it requires a priori information regarding the exact EZ location to train the model. Unsupervised ML techniques are attractive alternatives aiming to identify dominant networks in the iEEG data^18,35,36,37. These methods can identify relevant epileptiform patterns and inherent “structures” in the data without explicitly training the model³⁷. Thus far, these methods have only been applied to small cohorts and their potential use for outcome prediction remains unexplored.

Here, we introduce a novel AI-based framework that identifies the EZ with high precision and predicts surgical outcome in DRE patients. The framework is patient-specific and requires little to no input from clinicians; thus, it offers a significant advancement in the personalized surgical planning for DRE patients undergoing resective neurosurgery. The framework transforms interictal iEEG data into brain networks and their corresponding temporal maps in different frequency bands through dynamic signal decomposition (DSD) and unsupervised ML. The extracted networks characterize highly correlated and/or quasi-synchronized coherent areas that are repeatedly active across time in continuous interictal iEEG data. The temporal maps indicate when these networks are active in time. Our method classifies the derived networks into epileptogenic and background depending on the frequency of their activity in the temporal maps. We hypothesize that the identified IEN delineates the EZ with high precision; resection (or ablation) of this brain area yields good surgical outcome. To test our hypothesis, we compared the spatial distribution of the IEN with the resection and clinically defined SOZ in a relatively large cohort of children and young adults with DRE. We evaluated the ability of these networks to identify the resection and SOZ and to predict outcome. More specifically, we explored the concordance of each network’s active time-windows in temporal maps with IEDs and ripples, compared each network’s power inside vs. outside the resection and SOZ, and used support vector machine to predict outcome based on the IEN’s focality and proximity to resection. We further assessed the applicability of our framework for data with frequent and sparse IEDs as well as for different implantation types. This evaluation assesses the role of the IEN in delineating the EZ and predicting outcome in patients with DRE.

Patient cohort

We retrospectively analyzed interictal iEEG data from 43 children and young adults (39.5% female) with DRE having good (27 patients, Engel I) and poor (16 patients, Engel (>)I) outcome following epilepsy surgery [median follow-up: 4 years (2–6)]. Patient demographics are summarized in Table 1. The patients had a median age at surgery of 13 years (range: 9.3–16.8) with median seizure onset age of four years (range: 1–7.5). The median time between diagnosis and surgery was eight years (range: 3.6–11.8). Our dataset included patients with three implantation types: (i) subdural electrodes [electrocorticography (ECoG)] (15 patients; 10 with good outcome); (ii) stereotactic EEG (sEEG) implantation (10 patients; 6 with good outcome); and (iii) both subdural and depth electrodes (18 patients; 11 with good outcome). The median number of implanted electrodes was 99 (range: 88–119). The median percentage of resected electrodes was 18% (range: 8.8–31.1) and the percentage of electrodes determined to be part of the SOZ was 9% (range: 4.6–23.7). The median resected volume determined by coregistering the preoperative and postoperative MRIs was 23.1 cm³ (range: 11.5–36.9) in good- and 20.3 cm³ (range: 8.3–43.9) in poor-outcome patients. Eleven patients (26%) had non-lesional MRI, 26 patients (60%) had a malformation of cortical development, and six patients (14%) had acquired brain injury. Eight patients (19%) had temporal epilepsy. The median IEDs rate was 40 IEDs per minute (range: 13–58.5). No significant differences in these characteristics were found between good- and poor-outcome patients.

Construction of brain networks and temporal maps

The proposed framework transforms multimodal brain imaging data (preoperative and postoperative MRIs, CT-scans, and interictal iEEG) into interpretable spatiotemporal maps comprised of active brain networks and their corresponding activity across time in the form of temporal maps using DSD and unsupervised ML. By categorizing the networks as epileptogenic and background, we perform a comparative analysis to validate the ability of the IEN to delineate the EZ and predict outcome (Fig. 1).

a Coregistration of MRI and computed tomography with intracranial implantations identified intracranial EEG (iEEG) electrode coordinates. Pre- and post- operative MRIs defined the resected volume. The resection and clinically defined seizure onset zone (SOZ) were used as gold standards for the epileptogenic zone (EZ) prediction. b Five minutes of n-channel iEEG were filtered in two frequency bands: spike band (sb) [1-80 Hz] and ripple band (rb) [80–250 Hz] and dissected into d ms time-windows with 95% overlap. Dynamic mode decomposition (DMD) with (h) time delay embeddings extracted ({r}_{1}) (for sb) and ({r}_{2}) (for rb) oscillatory components and spectra per channel across time-windows. c Feature matrices were constructed by averaging the DMD spectra across seven physiologically relevant frequency bands [delta ((delta) = 1–4 Hz), theta ((theta) = 4–8 Hz), alpha ((alpha) = 8–12 Hz), beta ((beta) = 12–3 Hz), gamma ((gamma) = 30–80 Hz), spike band (({sb}) = 1–80 Hz), and ripple band (({rb}) = 80–250 Hz)], and scaled between 0 and 1. The entire network was computed by averaging the feature matrices across time-windows in each frequency band. Non-negative matrix factorization (NNMF) then extracted in each frequency band two brain networks and a temporal map. d Networks were categorized as epileptogenic (IEN) (red-colored) and background (blue-colored) whose focality, overlap, and distance from resection were then computed. e Temporal maps showed active IEN (red-colored) segments and background network (blue-colored) segments across time which were concordant with interictal epileptiform discharge (IED) (pink-colored) and ripple annotations (purple-colored).

Full size image

We initially processed the preoperative MRI, postoperative MRI, and electrode implantation CT to define the exact location of iEEG electrodes and resection in relation to brain anatomy (Fig. 1a). Next, we filtered the multichannel interictal 5-minute iEEG segment in spike (sb) (1–80 Hz) and ripple (rb) (80–250 Hz) bands (Fig. 1b). These two bands contain most interictal activity and correspond to the frequency components of IEDs and ripples^21,38. We then dissected the data in each band into time-windows and processed them using dynamic mode decomposition (DMD) (Fig. 1b). The DMD has been previously used to extract coherent brain regions by decomposing neural recordings in both space and time¹⁸. Since IEDs and ripples occur coherently in multiple regions^16,17,18, DMD is a suitable approach to characterize their dynamics. We used DMD to extract power features from interictal iEEG that characterize coherent channels oscillating at specific frequency components in the form of DMD power spectra (Fig. 1b). These components refer to common frequencies present across multichannel signals in a time-window. We further processed the DMD power spectra to construct feature data matrices containing average DMD powers in physiologically relevant frequency bands (Fig. 1c) [delta ((delta) = 1–4 Hz), theta ((theta) = 4–8 Hz), alpha ((alpha) = 8–12 Hz), beta ((beta) = 12–3 Hz), gamma ((gamma) = 30–80 Hz), spike band (({sb}) = 1–80 Hz), and ripple band (({rb}) = 80–250 Hz)] (see “Methods”).

To characterize the dynamics of interictal iEEG data and extract recurrent active networks across time-windows, we processed the feature data matrices in each band separately. We first derived the entire network by averaging DMD spectral power of channels across time-windows in each band. We then extracted brain networks and temporal maps (Fig. 1c). Brain networks represent coherent iEEG electrodes which are repeatedly active across time¹⁸. The active electrodes in the network are displayed in Fig. 1d as spheres whose radii are proportional to the activation strength. Temporal maps characterize when these networks are active across time; they are displayed in Fig. 1e as color-coded maps indicating which network is active during a specific time-window. We define ‘active time-windows’ as time periods in the temporal map where a network is active, and ‘inactive time-windows’ as periods where a network is not active. To extract the networks and their corresponding temporal map, we used an unsupervised ML approach, namely non-negative matrix factorization (NNMF)³⁹. Since IEDs and ripples occur transiently, alternating with background activity^10,11, we presumed that there are two networks which are active. We then categorized the less frequently active network as IEN and the more active as background (see “Methods”). Finally, we evaluated the ability of the networks (i.e., entire, interictal epileptogenic, and background) to delineate the EZ and predict outcome. The details of the framework are provided in “Methods”.

We first investigated whether time-windows corresponding to the IEN were concordant with the activity of IEDs and ripples (annotated by EEG experts) across various frequency bands (see “Methods”). Our analysis revealed that time-windows corresponding to the IEN on the temporal map had superior performance in detecting IEDs and ripples compared to the background network (p values < 0.0001, Wilcoxon signed-rank test) (Supplementary Note 1). Figure 1e depicts the temporal concordance of the active segments of the IEN with the annotated IEDs and ripples for different frequency bands in a portion of iEEG data from patient #1. Notably, time-windows of iEEG data with frank epileptiform activity (i.e., IEDs and ripples) were temporally concordant with those corresponding to the IEN.

To verify spatial stationarity of the identified networks over time, we assessed their consistency when networks were derived from different iEEG segments using the Dice score⁴⁰, a measure that quantifies the similarity and overlap between networks. A Dice score (ge)0.8 is considered to have almost perfect agreement⁴¹. We extracted the IEN and background networks from different 1-minute-long segments and computed their Dice scores with the networks identified from the remaining 4-minute-long segments (see “Methods”). Both the IEN and background network were consistent in terms of spatial distribution in all bands when different segments were analyzed with a median Dice score (ge)0.8 (Supplementary Figure 1).

We finally performed a randomization test to evaluate whether the IENs were meaningful and not simply due to random fluctuations (or noise) in the data. Therefore, we randomized the feature data matrices, extracted random networks using our framework, and compared them with the original IEN and background network using the Dice score (see “Methods”). The resulting Dice scores were consistently low ((le)0.3), indicating minimal overlap between the random networks and the IEN and background network (Supplementary Figure 1). Therefore, our framework generates robust IENs that are unlikely due to random data fluctuations.

Higher power of IEN inside resection and SOZ

To assess the relationship between the network power distribution and EZ, we estimated and compared the power of the entire, interictal epileptogenic, and background networks inside vs. outside the resection and SOZ in good- and poor-outcome patients, separately. In good-outcome patients, we observed increased power with moderate effect size for the entire network inside (compared to outside) resection for (delta)(p = 0.0004, d = 0.42) and (theta) (p = 0.0005, d = 0.34) bands (Fig. 2a). For these patients, we also observed higher power for the IEN inside resection for (delta) (p = 0.0011, d = 0.34), (theta) (p < 0.0001, d = 0.55), (alpha) (p = 0.0012, d = 0.33), (beta) (p = 0.0001, d = 0.47), and ({sb}) (p = 0.0066, d = 0.40) (Fig. 2a). We found no differences between inside and outside resection for the background network in good-outcome patients for any band (Fig. 2a). Poor-outcome patients did not show differences in power for the entire, interictal epileptogenic, and background networks across all bands, except for the entire network in (delta) band (Fig. 2a).

a Median spectral power inside vs. outside resection for the entire, interictal epileptogenic (IEN), and background networks in different frequency bands in good- (green-colored) and poor- (orange-colored) outcome patients. b Median dynamic mode decomposition spectral power inside vs. outside the SOZ for the entire, interictal epileptogenic (IEN), and background networks in different frequency bands in good- (green-colored) and poor- (orange-colored) outcome patients. The notch of the box plots represents the median values, lower and upper edges represent the 25^th and 75^th percentiles. Whiskers extend to the minimum and maximum values after omitting the outliers. The colored circles represent the median power of the networks of different patients and the black lines connect paired values. Significant differences are indicated by asterisks: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001 (Wilcoxon signed-rank test, significance levels were corrected for multiple comparisons using Bonferroni correction). The number displayed below the asterisks represents the effect size. Frequency bands: delta ((delta) = 1–4 Hz), theta ((theta) = 4–8 Hz), alpha ((alpha) = 8–12 Hz), beta ((beta) = 12–3 Hz), gamma ((gamma) = 30–80 Hz), spike band (({sb}) = 1–80 Hz), and ripple band (({rb}) = 80–250 Hz).

Full size image

Good-outcome patients showed higher power for the entire network inside (vs. outside) the SOZ in (delta) (p = 0.0002, d = 0.42), (theta) (p = 0.0006, d = 0.39), and (alpha) (p = 0.0036, d = 0.32) bands (Fig. 2b). Similarly, the IEN showed higher power inside the SOZ for (delta) (p = 0.0025, d = 0.30), (theta) (p < 0.0001, d = 0.55), (alpha) (p = 0.0002, d = 0.45), (beta) (p < 0.0001, d = 0.53), (gamma) (p = 0.0042, d = 0.37), and ({sb}) (p = 0.0008, d = 0.53) (Fig. 2b). We found no differences between inside and outside the SOZ for the background network in good-outcome patients for any band (Fig. 2b). Notably, poor-outcome patients did not show differences in power for the entire, interictal epileptogenic, and background networks across all bands (Fig. 2b).

Network properties

For each network (i.e., entire, interictal epileptogenic, and background) of good-outcome patients, we estimated and compared three properties: focality measuring the proximity of the identified network electrodes to each other (F_net), overlap with resection (O_res), and distance from resection (D_res) (see “Methods”). We presume that the IEN is more focal, has larger overlap, and is closer to the resection compared to the background network. The IEN showed higher F_net (in (theta), (beta), (gamma), and sb; p < 0.05, d ≥ 0.38), higher O_res (in all bands; p < 0.01, d ≥ 0.40), and lower D_res (in all bands; p < 0.05, d ≥ 0.33) compared to the background network (Fig. 3). Compared to the entire network, the IEN had higher F_net (in all bands except (delta); p < 0.05, d ≥ 0.38), higher O_res (in all bands except (delta) and (alpha); p < 0.05, d ≥ 0.26), and lower D_res (in all bands except (delta) and (alpha); p < 0.01, d ≥ 0.32) (Fig. 3). Finally, the entire network had higher O_res (in (delta), (theta), and (alpha); p < 0.01, d ≥ 0.27) and lower D_res (in (theta) and (alpha); p < 0.01, d ≥ 0.06) compared to the background (Fig. 3). There were no differences in F_net between the entire and background networks in all bands (Fig. 3).

Focality (({F}_{{net}}), normalized values), percentage of overlap with resection (({O}_{{res}}), %), and distance from resection (({D}_{{res}}), mm) of the entire, interictal epileptogenic (IEN), and background networks of good-outcome patients in different frequency bands. The notch of the box plots represents the median values, lower and upper edges represent the 25^th and 75^th percentiles. Whiskers extend to the minimum and maximum values after omitting the outliers. Significant differences are indicated by asterisks: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001 (Wilcoxon signed-rank test, significance levels were corrected for multiple comparisons using Bonferroni correction). The number displayed below or to the right of the asterisks represents the effect size. Frequency bands: delta ((delta) = 1–4 Hz), theta ((theta) = 4–8 Hz), alpha ((alpha) = 8–12 Hz), beta ((beta) = 12–3 Hz), gamma ((gamma) = 30–80 Hz), spike band (({sb}) = 1–80 Hz), and ripple band (({rb}) = 80–250 Hz).

Full size image

The IEN predicts the EZ

We then investigated whether the entire, interictal epileptogenic, and background networks could predict the resection and SOZ in good-outcome patients by comparing the network performance metrics. We first used the power of each network’s electrodes as predictor and the resection as target. In terms of EZ prediction, the IEN outperformed the entire and background networks and was able to delineate the resection with 65% accuracy (range: 54–72%) and 85.7% precision (range: 50–100%) in (theta) band (Fig. 4). Since resection often contains areas larger than the actual EZ and not all electrodes inside resection exhibit epileptogenic activity, the observed low accuracy (but high precision) of the IEN in the (theta) band is justified. We performed a similar analysis using the SOZ as a target. Although higher accuracy was achieved compared to when resection was used as a target [92% (range: 83–96%) in (theta)], the networks poorly performed in terms of precision [49% (range: 25–66%) in (theta)] (Supplementary Figure 2).

Sensitivity, specificity, positive (PPV) and negative predictive values (NPV), and accuracy of the entire (purple-colored), interictal epileptogenic (IEN) (red-colored), and background (blue-colored) networks to predict resection. The notch of the box plots represents the median values, lower and upper edges represent the 25^th and 75^th percentiles. Whiskers extend to the minimum and maximum values after omitting the outliers. Significant differences are indicated by asterisks: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001 (Wilcoxon signed-rank test, significance levels were corrected for multiple comparisons using Bonferroni correction). Frequency bands: delta ((delta) = 1–4 Hz), theta ((theta) = 4–8 Hz), alpha ((alpha) = 8–12 Hz), beta ((beta) = 12–3 Hz), gamma ((gamma) = 30–80 Hz), spike band (({sb}) = 1–80 Hz), and ripple band (({rb}) = 80–250 Hz).

Full size image

In terms of EZ prediction, the derived IEN was able to delineate the resection with an average receiver operating characteristic (ROC) area under the curve (AUC) of 0.7440 ((sigma) = ±0.11) in (theta) band (Fig. 5a). The average AUC with resection in all other bands was <0.7 (Fig. 5a). ROC curves of all good-outcome patients were above the no-discrimination line (i.e., line equivalent to chance) only in (theta) band (Fig. 5a). Similarly, the IEN was able to delineate the SOZ with an average AUC of 0.7474 ((sigma) = ±0.18) in (theta) band (Fig. 5a). The average AUC with SOZ in all other bands was <0.7 (Fig. 5a).

a Receiver operating characteristic (ROC) curves of interictal epileptogenic networks (IENs) as predictors of resection and SOZ for all good-outcome patients across different frequency bands. The ROC curve of each patient is shown in different color. The mean ROC curves are depicted in bold black. The diagonal line represents the line of no-discrimination or chance [area under the curve (AUC) equals 0.5]. The average AUC values are shown at the top of each panel for each band. b AUC of the entire, interictal epileptogenic (IEN), and background networks with resection and SOZ in different bands. The notch of the box plots represents the median values, lower and upper edges represent the 25^th and 75^th percentiles. Whiskers extend to the minimum and maximum values after omitting the outliers. Significant differences are indicated by asterisks: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001 (Wilcoxon signed-rank test, significance levels were corrected for multiple comparisons using Bonferroni correction). Frequency bands: delta ((delta) = 1–4 Hz), theta ((theta) = 4–8 Hz), alpha ((alpha) = 8–12 Hz), beta ((beta) = 12–3 Hz), gamma ((gamma) = 30–80 Hz), spike band (({sb}) = 1–80 Hz), and ripple band (({rb}) = 80–250 Hz).

Full size image

In predicting resection, the IEN showed superior performance with higher AUC compared to background (in (theta), (beta), (gamma), and sb; p < 0.01) and entire (in (theta) and (beta); p < 0.05) networks (Fig. 5b). The entire network had higher AUC compared to background (in (theta) and (beta); p < 0.01) (Fig. 5b). In predicting the SOZ, the IEN showed higher AUC when compared to background (in (theta), (alpha), (beta), (gamma), and sb; p < 0.05) and entire (only in (theta); p < 0.01) networks (Fig. 5b). Additionally, the entire network had higher AUC compared to background when predicting the SOZ (in (theta), (alpha), (beta), and sb; p < 0.05) (Fig. 5b).

Robust IENs across variable IED rates and implantation types

We investigated the robustness of the IEN across segments with frequent and sparse IEDs for the good-outcome patients. Only 16 out of the 27 good-outcome patients had segments of 1-minute duration with both frequent (57.2 IEDs per minute) and sparse IEDs (8.7 IEDs per minute) (Supplementary Table 1). For these patients, we estimated the IENs using our framework and calculated the AUC of the IENs with the resection and SOZ. We then compared the AUC values derived from the segments with frequent IEDs with those derived from the segments with sparse IEDs (see “Methods”). We found no differences between the AUC values from segments with frequent vs. sparse IEDs (p > 0.05) (Fig. 6a).

a Area under the curve (AUC) of IENs with resection and SOZ in segments with frequent and sparse interictal IEDs. b AUC of IENs with resection and SOZ in patients with different implantation types. The notch of the box plots represents median values; lower and upper edges represent the 25^th and 75^th percentiles. Whiskers extend to the minimum and maximum values after omitting the outliers. Significant differences are indicated by asterisks: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001 (Wilcoxon signed-rank test for rate consistency analysis and Wilcoxon rank-sum test for implantation consistency analysis, significance levels were corrected for multiple comparisons using Bonferroni correction). Frequency bands: delta ((delta) = 1–4 Hz), theta ((theta) = 4–8 Hz), alpha ((alpha) = 8–12 Hz), beta ((beta) = 12–3 Hz), gamma ((gamma) = 30–80 Hz), spike band (({sb}) = 1–80 Hz), and ripple band (({rb}) = 80–250 Hz). ECoG = electrocorticography; sEEG = stereotactic electroencephalography.

Full size image

We then examined whether our framework provides consistent findings among the three implantation types in our dataset: (i) subdural electrodes; (ii) sEEG implantations, and (iii) both subdural and depth electrodes. For good-outcome patients, we initially estimated the IENs and then calculated their AUC with the resection and SOZ. We then compared the AUC of IENs with resection and SOZ among the three implantation types (see “Methods”). We found no differences between the AUCs of IENs with resection among the different implantation types in any frequency band (p > 0.05) except the (alpha) band (p = 0.0023) (Fig. 6b). We also found no differences between the AUCs for predicting the SOZ (p > 0.05) (Fig. 6b).

Focality and proximity of IEN to resection

We examined the focality and proximity of the IEN to resection. We presumed that the IEN would have higher F_net, lower O_res, and lower D_res to resection in good- compared to poor-outcome patients. At the population level, the IEN had higher F_net in good-outcome patients (compared to poor) in (theta) (p = 0.0002, d = 0.69), (alpha) (p = 0.0051, d = 0.52), and (beta) (p = 0.0040, d = 0.53) bands (Fig. 7a). The IEN also had higher O_res in good-outcome patients in (theta) (p = 0.0003, d = 0.54) and (beta) (p = 0.0013, d = 0.34) bands (Fig. 7a). Finally, the D_res of the IEN was lower in good-outcome in (theta) (p = 0.0004, d = 0.54), (beta) (p = 0.0022, d = 0.42), and (gamma) (p = 0.0037, d = 0.54) bands compared to poor-outcome patients (Fig. 7a).

a Network focality (({F}_{{net}}), normalized values), percentage of overlap with resection (({O}_{{res}}), %), and distance from resection (({D}_{{res}}), mm) of the IENs in different frequency bands in good- (green-colored) and poor- (orange-colored) outcome patients. The notch of the box plots represents the median values, lower and upper edges represent the 25^th and 75^th percentiles. Whiskers extend to the minimum and maximum values after omitting the outliers. Significant differences are indicated by asterisks: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001 (Wilcoxon rank-sum test, significance levels are corrected for multiple comparisons using Bonferroni correction). The number displayed below the asterisks represents the effect size. b Properties of the IEN across patients and outcomes in different frequency bands. For visualization purposes, the reciprocal of the distance from resection ((dagger)) is shown. All values were normalized between 0 and 1. Circle radii are proportional to the normalized values with larger circles indicating higher focality, higher overlap, and proximity to resection. The bar graphs on the sides of the panels for each property represent the number of cases that were above the corresponding thresholds ({th}) (where ({th}) represents (t{h}_{F}) = 0.48 for ({F}_{{net}}), (t{h}_{O}) = 49% for ({O}_{{res}}), (t{h}_{D}) = 18.25 mm for ({D}_{{res}}),) for good- (left) and poor-outcome (right) patients per frequency band. Frequency bands: delta ((delta) = 1–4 Hz), theta ((theta) = 4–8 Hz), alpha ((alpha) = 8–12 Hz), beta ((beta) = 12–3 Hz), gamma ((gamma) = 30–80 Hz), spike band (({sb}) = 1–80 Hz), and ripple band (({rb}) = 80–250 Hz).

Full size image

At the patient level, we report the IEN properties for all patients aiming to evaluate differences between good- and poor-outcome patients. The F_net and O_res were higher in most patients with good (compared to poor) outcome (Fig. 7b). Similarly, D_res was lower in most patients with good outcome. We report F_net, O_res, and D_res for each patient (Fig. 7b) after estimating optimal thresholds (see “Methods”) of each property for predicting outcome ((t{h}_{F}) = 0.48, (t{h}_{O}) = 49%, and (t{h}_{D}) = 18.25 mm, respectively). The thresholds for each band are reported in Supplementary Table 2. We observed that 20 out of 27 (74.1%) good-outcome patients had F_net ≥(t{h}_{F}) in (theta) band (Fig. 7b); contrarily, only 4 out of 16 (25%) poor-outcome patients had F_net ≥(t{h}_{F}). We also found that 85.1% of good-outcome patients (23 out of 27) had O_res ≥(t{h}_{O}) with resection in (theta) band (Fig. 7b); contrarily, only 31.2% (5 out of 16) of poor-outcome patients had O_res ≥ (t{h}_{O}) with resection. Finally, we found that 92.6% of good-outcome patients (25 out of 27) had D_res ≤ (t{h}_{D}) in (theta) band; contrarily, only 37.5% of poor-outcome patients (6 out of 16) had D_res ≤ (t{h}_{D}) (Fig. 7b).

IEN predicts surgical outcome

To evaluate the IEN as interictal biomarker of the EZ, we examined its ability to predict outcome using network properties as predictors. Particularly, we trained a linear support vector machine classifier (SVM)⁴² to estimate the outcome using the three IEN properties as features and outcome as target in each of the seven bands. We also trained another linear SVM classifier (denoted by SVM-all) using 21 properties (seven bands × three properties per band) simultaneously as features that incorporated information from all bands combined. Each SVM model was validated using five-fold cross-validation. SVM outputs the probability of being classified as good and poor (see “Methods”). Similarly, for each band, we trained SVMs using entire network properties (as well as SVM-all with 21 properties) as outcome predictors.

Predictive models using SVM revealed that the entire network predicted outcome in (theta) and sb (p < 0.05) (Fig. 8a). Moreover, the IEN demonstrated predictive capabilities in multiple bands, including (theta), (beta), (gamma), and sb (p < 0.05) (Fig. 8b). Overall, the IEN outperformed the entire network in predicting outcome. In particular, SVM trained using the IEN in (theta) band performed best with 93% sensitivity, 63% specificity, 81% precision [positive predictive value (PPV)], 83% negative predictive value (NPV), 81% accuracy, and 0.85 AUC (Fig. 8c). On the other hand, SVM trained using the entire network in (theta) band had 85% sensitivity, 44% specificity, 72% PPV, 64% NPV, 70% accuracy, and 0.71 AUC (Fig. 8c).

a Confusion matrices of the support vector machine (SVM) classifier using three properties (focality, overlap, and distance from resection) of the entire network in each band as predictors of outcome. “All” (checkered red-colored) is the result of training SVM-all using all properties of all bands as predictors. Numbers in parentheses represent p values (Fisher’s exact test). Significant differences are indicated in bold. b Confusion matrices of SVM using three network properties of the interictal epileptogenic network (IEN) in each band as predictors. “All” (red-colored) is the result of training SVM-all classifier. Numbers in parentheses represent p values (Fisher’s exact test). Significant differences are indicated in bold. c Bar graphs depict the performance metrics in percentage [sensitivity, specificity, positive and negative predictive values (PPV and NPV, respectively), accuracy, and five-fold cross-validation area under the curve (AUC)] of the confusion matrices of (a) and (b). In each frequency band, the first bar represents the entire network (checkered pattern), while the second bar represents the IEN (fully-colored). d Feature importance analysis using ANOVA where (F), (O), and (D) represent focality, overlap, and distance from resection, respectively. Significant features are indicated by asterisks: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001 (corrected for multiple comparisons using Bonferroni correction). Frequency bands: delta ((delta) = 1–4 Hz), theta ((theta) = 4–8 Hz), alpha ((alpha) = 8–12 Hz), beta ((beta) = 12–3 Hz), gamma ((gamma) = 30–80 Hz), spike band (({sb}) = 1–80 Hz), and ripple band (({rb}) = 80–250 Hz).

Full size image

SVM-all trained using all 21 IEN properties predicted outcome (p = 0.0007) (Fig. 8b) with 89% sensitivity, 63% specificity, 81% PPV, 77% NPV, 79% accuracy, and 0.86 AUC (Fig. 8c). SVM-all trained using the entire network properties predicted outcome (p = 0.0104) (Fig. 8a) with 78% sensitivity, 63% specificity, 78% PPV, 63% NPV, 72% accuracy, and 0.75 AUC (Fig. 8c).

Finally, we performed feature importance analysis using analysis of variance (ANOVA)⁴³ to rank the IEN properties in the seven bands in decreasing order of importance and identified the IEN properties and bands that best discriminated between good and poor outcome. We found that F_net, O_res, and D_res in (theta) band had the highest importance for predicting outcome followed by F_net and O_res in (beta), and D_res in (gamma) band (Fig. 8d).

Representative cases

To highlight the ability of the IEN to predict the EZ and outcome at the patient level, we report findings from two representative cases: a patient with good and one with poor outcome. For these patients, we depict a snapshot of annotated iEEG data, the corresponding temporal maps (Fig. 9a), and the IEN and the background network in (theta) band (Fig. 9b). We also report network properties, their predictive values for resection and SOZ, and the concordance of temporal maps with IEDs (Fig. 9c).

a Five-second 10-channel sample intracranial EEG data with annotated spikes (pink-colored) and corresponding temporal maps showing activity of the background network (blue-colored) and interictal epileptogenic network (IEN) (red-colored) for a good- (patient #1, 11-year-old male) and a poor- (patient #31,16-year-old male) outcome patient. b The background (blue-colored) and IEN (red-colored) networks in (theta) band projected on MRI with the resection (green-colored) and zoomed around the resection for both good- and poor-outcome cases. The spheres were centered at the network’s electrodes and their radii were proportional to the dynamic mode decomposition (DMD) spectral power (only electrodes greater than the threshold are shown). c For both patients, IEN and background network properties [i.e., focality (({F}_{{net}})), overlap (({O}_{{res}}), %) and distance from resection (({D}_{{res}}), mm)] are displayed for the (theta) band. Performance metrics [i.e., area under the curve (AUC) with resection (AUC-RES) and seizure onset zone (AUC-SOZ), and concordance with interictal epileptiform discharges (IEDs) (AUC-IED)] for both the background network (blue-colored) and IEN (red-colored) are displayed for the (theta) band. Frequency bands: delta ((delta) = 1–4 Hz), theta ((theta) = 4–8 Hz), alpha ((alpha) = 8–12 Hz), beta ((beta) = 12–3 Hz), gamma ((gamma) = 30–80 Hz), spike band (({sb}) = 1–80 Hz), and ripple band (({rb}) = 80–250 Hz).

Full size image

For the good-outcome patient, we observed higher F_net (0.74 >(t{h}_{F})), higher O_res (100% >(t{h}_{O})), and lower D_res (6 mm <(t{h}_{D})) for the IEN (compared to background; F_net (=)0.51, O_res (=)18%, and D_res=31 mm) in (theta) band (Fig. 9c). The IEN had higher AUC with resection (0.92) and SOZ (0.93) compared to background network (0.62 and 0.43, respectively) (Fig. 9c). In terms of network’s temporal map concordance with IED annotations (Fig. 9a), the IEN had higher AUC (0.75) compared to background (0.24) (Fig. 9c).

For the poor-outcome patient, both networks in (theta) band had F_net <(t{h}_{F}), O_res <(t{h}_{O}), and D_res >(t{h}_{D}) (Fig. 9c). The IEN was able to predict resection with an AUC = 0.76. Both networks had low AUC when predicting the SOZ (0.53 and 0.63, respectively) (Fig. 9c). In terms of network’s temporal map concordance with IED annotations (Fig. 9c), both networks had low AUC (0.43 and 0.58, respectively). Ultimately, the IEN identified candidate resection areas. Notably, while these areas overlapped with resection in the good-outcome patient, they covered wider areas and were far away from resection in the poor-outcome patient.

We propose a novel patient-specific framework that extracts the IEN from interictal iEEG data using unsupervised ML, delineates the EZ, and predicts outcome in patients with DRE. The proposed framework automatically transforms the interictal iEEG data into temporal maps and brain networks and categorizes them as background and epileptogenic. The IEN presents high-power coherent activity which corresponds to IEDs and ripples in temporal maps; surgical resection of this network predicts good outcome. Its performance is better than the one of the entire network. Our framework works on both data with frequent and sparse IEDs as well as different implantation types reducing the burden of manually identifying epileptogenic iEEG time-windows and localizing the corresponding underlying epileptogenic generator.

Interictal biomarkers in the form of elevated power, coherence, source-sink, or functional connectivity have been associated with the EZ^{16,21,25,30,44,45,46}. However, these biomarkers rely on manual selection of IEDs by experts, which is time-consuming and prone to errors¹⁵. AI-based tools that delineate the EZ have been recently proposed;^33,47 yet, they are not widely used in clinical practice possibly because they have been trained and validated only in small potentially biased datasets. Moreover, they rely on information about the EZ provided by experts as a priori, which is prone to errors⁴⁸. Furthermore, in many cases, these methods are validated without properly separating testing and training sets, which can lead to overfitting; thus, they suffer from poor generalization and may not perform well on new unseen data⁴⁹. Our framework overcomes these hurdles by identifying coherent IENs without the need to detect IEDs or ripples. It adopts an unsupervised ML approach^18,35,36 which does not rely on a model that is trained using prior information. Moreover, it separates epileptiform from background activity, enabling the identification of epileptogenic regions without the need for manual data processing; thus, it enhances the localization of the EZ in comparison to previously described methods^17,50,51, which rely on aggregated computations such as averaging. This notion is validated by our findings showing the superiority of the IEN (compared to the entire) to identify critical epileptogenic areas. The entire and IEN had higher power inside resection in several bands for patients with good outcome (Fig. 2). Moreover, the IEN showed higher focality, greater overlap with resection, and shorter distance to resection compared to the entire network (Fig. 3). This suggests that the IEN is more focal and spatially overlapped with the EZ compared to the entire network. The IEN is superior to the entire network to delineate the EZ achieving higher precision (Fig. 4) and AUC (Fig. 5). Therefore, although the entire network predicts the EZ, the IEN outperforms it.

Previous studies showed the association of elevated interictal power and functional connectivity with epileptogenicity and surgical outcome^16,17,34,45. These studies suggest that seizure-generating regions in focal epilepsy patients are isolated (from surrounding regions) and that the likelihood of good outcome is increased by resecting these regions. Previously, the overlap of the identified regions with resection was used as outcome predictor^16,20,21,50. Our findings support this notion since we observed a decrease in percentage overlap of the IEN with resection in poor-outcome patients (Fig. 7). Furthermore, we found that the distance of the IEN from resection was higher in poor-outcome patients (Fig. 7). We also observed a decrease in the IEN focality in poor-outcome patients (Fig. 7). Our findings suggest that the more focal the epileptogenic activity is the higher the chances are for the patient to become seizure-free if these regions are resected. Therefore, in addition to overlap, we employed the IEN focality and its distance from resection as outcome predictors. We also compared the outcome prediction of the IENs and entire networks. Our findings suggest that, although the entire network could predict outcome, the IEN is superior in terms of sensitivity (93% vs. 88%), PPV (81% vs. 78%), NPV (83% vs. 64%), accuracy (81% vs. 72%), and AUC (86% vs. 75%) (Fig. 8). In 25 (out of the 27) good-outcome patients, there was an agreement between the EZ (identified by the IEN) and the one delineated by the clinicians (Fig. 8). In poor-outcome patients, the IEN predicted outcome correctly in 10 out of 16 cases (Fig. 8). These findings may indicate other epileptogenic regions, distant from resection, that have not been resected due to overlap of the IEN with eloquent areas [e.g., patient #31 (Fig. 9)].

IEDs are often characterized by different morphologies that occur at varying durations. Spikes often last between 20 and 70 ms, while sharp spikes have a duration of 70–200 ms^11,52. In the spectral (frequency) domain, these durations are equivalent to components oscillating in (theta) and (beta) bands. Here, we found that the IEN in (theta) (and (beta)) bands outperformed other bands when identifying the EZ and predicting outcome. The IEN in the (theta) band identified resection with a precision of 85.7%. This is probably due to the concordance of the IENs temporal activity with IEDs whose morphological characteristics were predominantly expressed in the form of (theta) waves. Additionally, ANOVA during outcome prediction showed that the IEN in (theta) band had the best discriminative power between good- and poor-outcome patients when all networks in all bands were used as outcome predictors. Our findings are in line with previous studies showing increased functional connectivity in (theta) band in patients with epilepsy when compared to healthy controls^53,54. Increased interictal synchronization in (theta) band may not necessarily be a unique characteristic of epilepsy. Instead, it could be a consequence of compensatory mechanisms or disinhibition resulting from the disease. As the brain attempts to compensate for the damage caused by the disease, it may increase the synchronized activity in (theta) band as a compensatory mechanism⁵⁵. Several studies have related (theta) activity to abnormal plasticity in epilepsy patients with lesions; in these patients, abnormal neural connections may lead to excessive synchronization in (theta) band^53,56.

Automated IED detectors^4,15, including AI-based tools^57,58, have been designed to detect IEDs in iEEG recordings by exhibiting varying performances (specificity: 68–94%, sensitivity: 47–99%, and accuracy: 68–100%)^36,59. These tools do not have the ability to assess the spatial extent of the EZ without further processing. Tools that allow the spatiotemporal mapping of interictal data have the potential to identify the spatial and temporal extent of epileptiform activity^18,35,36. In line with this notion, our framework identified the IEN whose temporal activity was concordant with IED and ripple annotations with 77% accuracy (in sb) (Supplementary Figure 1). This was achieved by analyzing interictal iEEG recordings of short duration; this is in line with previous studies showing that short time-windows are sufficient to map the spatiotemporal dynamics of iEEG data with high consistency^18,35,36,45. Additionally, our methodology robustly estimated the IENs regardless of the IEDs frequency (Fig. 6a).

Previous studies have shown that the interpretation of epileptiform activity can be influenced by the implantation type used in each patient^2,26,27,60. In particular, sEEG employs a network-based implantation strategy that differs fundamentally from subdural and depth electrodes^26,27. sEEG has the ability to capture propagating epileptiform activity beyond the EZ leading to more distributed activity across multiple brain regions. On the other hand, subdural electrodes record activity exclusively at the cortical level⁶¹. Both modalities have limited spatial resolution since they record epileptiform activity in the direct vicinity of contacts and thus are blind to other brain regions. Therefore, the actual focus of the EZ may be missed leading to surgical failure. Here, we examined whether the implantation type can affect the robustness of our methodology. Our analysis showed no differences between the AUCs of IENs with resection among the different implantation types in any frequency band except the (alpha) band (Fig. 6b) and no differences for the SOZ. The differences between the AUCs of the IENs of subdural and sEEG implantations in the (alpha) band may be explained by the fact that subdural electrodes are unable to capture dysregulations of thalamocortical interactions observed in patients with epilepsy in (alpha) band⁶², which may be captured with sEEG. In summary, our data demonstrate that our framework provides robust findings independent of the implantation type.

Our retrospective study analyzed iEEG data from a single-center cohort of patients with DRE having different pathologies. Future prospective studies from larger multicenter cohorts would allow us to examine its applicability in clinical settings and in homogeneous groups of patients. Our analysis was performed on short data segments; analysis of longer data segments may enhance the accuracy of the IEN to delineate the EZ and predict outcome. Due to the lack of interictal data segments with absent IEDs in our database, further studies should examine the applicability of our methodology for patients where IEDs are absent in intracranial EEG recordings. Moreover, iEEG does not cover the entire brain but is rather implanted in locations agreed upon during pre-implantation analysis which can be subjective. While most of the implantation correctly covered the EZ in good-outcome patients, epileptogenic regions may have been missed in poor-outcome patients. Furthermore, the resection and SOZ were subjectively identified and may either miss critical epileptogenic regions or describe areas larger than the actual EZ. Advances in whole brain iEEG implantations may overcome these drawbacks. Our framework may not discriminate the EZ from IEDs propagation. Future studies may incorporate the phase information of DMD in order to distinguish between the primary epileptogenic activity and propagating IEDs. Additionally, the phase information may provide directionality of IED propagations providing better understanding of the epileptogenic network dynamics. Future studies could consider applying the proposed framework to a large range of virtual channels, estimated through electrical source imaging performed on iEEG recordings. This approach allows the characterization of neural activity in brain regions which are not directly sampled by the sEEG^21,63. Thus, it may allow capturing the entire phenomena of IEDs propagation potentially discriminating it from the primary EZ. Finally, our framework presumes that only two networks were active. Future studies should examine the exact number of networks that were active during interictal iEEG and duration of data analysis under different scenarios. Future directions may include clinically applicable 3D surgical navigation systems that integrate the IEN enabling the precise identification of iEEG contacts to resect or ablate.

To conclude, our proposed framework identifies an IEN which delineates the EZ and predicts surgical outcome from interictal iEEG data of patients with DRE. The framework’s automated nature reduces post-acquisition preprocessing cost and time, reducing risks associated with the presurgical evaluation. Such a framework would be particularly useful to epilepsy centers that lack the multidisciplinary expertise to delineate accurately and precisely the EZ in complex cases of patients with DRE. If validated prospectively, our framework may potentially replace ictal recordings to localize the EZ.

Ethics statement

The protocol was approved by North Texas Regional IRB (2019-166; PI: C. Papadelis) that waived the need for informed consent considering the study’s retrospective nature. All methods and analyses were performed in accordance with relevant guidelines and regulations.

Study cohort

We retrospectively analyzed iEEG data recorded from 43 children and young adults with DRE who had resective neurosurgery at Boston Children’s Hospital between June 2011 and June 2018. We selected patients based on the following criteria: (i) availability of at least 5-minute interictal iEEG data; (ii) availability of post-implantation computerized tomography (CT); (iii) availability of preoperative and postoperative MRIs; (iv) information about the resection volume and the clinically defined SOZ; and (v) availability of post-surgical outcome at least one year after surgery. The outcome was determined by a pediatric epileptologist (J.B.) after follow-up visits. We use the Engel score to classify patients as good (Engel(=)I, seizure-free) or poor outcome (Engel(>)I, non-seizure-free). The patient demographic and clinical information are provided in Supplementary Table 1.

Interictal iEEG recordings

The multidisciplinary clinical team decided the locations of implanted iEEG electrodes independently of this study. Long-term iEEG data were acquired with subdural electrodes, sEEG implantations, and subdural and depth electrodes using XLTEK NeuroWorks (Natus Inc., USA). The number and type of implanted electrodes are presented in Supplementary Table 1. Subdural electrodes were 2–3 mm in diameter with a 10 mm inter-contact distance, whereas depth electrodes were composed of 6 to 16 linearly arranged contacts ~1.5–2.5 mm apart. Each contact had a diameter of 0.8 mm and a length of 2 mm. The average acquisition time using iEEG implantations was 5.6 days, 12.6 hours, and 25.3 minutes. The sampling frequency ranged between 1000 and 2048 Hz. We selected 5-minute interictal segments with frequent IEDs from non-rapid eye movement slow-wave sleep (whenever applicable) at least one hour before/after clinical seizures or half an hour before/after electrographic seizures⁶⁴. This selection ensured the inclusion of segments with the highest IED rate and minimal motion artifacts^65,66.

Localization of iEEG electrodes

MRI scans with standard magnetization-prepared rapid acquisition gradient-echo sequences were performed before and after resection using a 3 T scanner (TIM TRIO, Siemens AG). CT scans were performed after iEEG implantation. We localized the iEEG electrode coordinates by coregistering the post-implantation CT scans with the preoperative MRI using Brainstorm (Fig. 1a)⁶⁷. We adjusted the electrode locations to compensate for possible brain shift occurring after electrocorticography implantation⁶⁸.

Defining the resection and SOZ

To define resection, we coregistered the preoperative and postoperative MRIs and manually drew the resection volume boundary on consecutive slices using Brainstorm (Fig. 1a)⁶⁷. Pediatric epileptologists defined the SOZ through visual inspection of ictal data independently from this study. iEEG contacts whose dynamics changed with seizure onsets were marked as SOZ contacts. To define the EZ, we considered as gold standards the resected electrodes and the clinically defined SOZ electrodes (Fig. 1a).

Preprocessing of iEEG data

We initially preprocessed the data by applying a DC offset removal and common average referencing. We then notch-filtered the data to remove the 60 Hz power line noise and its harmonics, and bandpass filtered them in two frequency bands: ({sb}) (1–80 Hz) and ({rb}) (80–250 Hz)³⁸. Bad channels and channels with artifacts were excluded from this study.

Feature extraction using DMD

The first step of our framework extracts coherent spatial features from interictal iEEG data. Using a sliding window approach, each time-window is processed with DMD⁶⁹. DMD decomposes high-dimensional time-series data into a set of coherent structures that exhibit similar linear dynamics in time, in the form of oscillations and exponential growth and decay⁷⁰. Since the spatial coherent patterns identified by DMD are relative to inherent temporal dynamics, these structures can be considered as functionally connected networks in terms of coherence^18,71.

DMD approximates locally (in a specific time-window) a non-linear dynamical system linearly in the discrete domain as:

where ({x}_{t}in {{mathbb{R}}}^{n}) denotes the measurements from (n) channels at instant (t), (A) is an (ntimes n) matrix, (m) is the number of samples, and ({w}_{t}) represents the residual error. DMD finds in the least square sense, a low-rank eigen-decomposition of (A) by minimizing:

The matrix (Xin {{mathbb{R}}}^{ntimes m}) resulting from horizontally stacking (m) measurements taken every (Delta t) can be represented as:

(X) can represent a time-window of iEEG with (n)-channels and m samples.

Constructing (X{prime} in {{mathbb{R}}}^{ntimes m}) which contains one-time shifted version of (X):

Equation (1) can be written as:

The solution of the least squares problem approximates A:

where † is the Moore-Penrose pseudo-inverse. The eigenvectors and eigenvalues of (A) correspond to DMD modes and eigenvalues. DMD computes a lower dimensional representation of the data to find an approximation (widetilde{A}) of (A) by keeping the leading (r) eigenvalues, eigenfunctions, and modes. The output of DMD is the DMD tuple ((lambda ,psi ,phi )) defined by:

1. 1.

   The eigenvalues ({lambda }_{j}) from which the frequency of oscillations ({f}_{j}) of the eigenfunctions can be computed where (scriptstyle{f}_{j}=|frac{{imag}({omega}_{j})}{2pi}|), ({omega }_{j}=log ({lambda }_{j})/Delta t), ({imag}(.)) is the imaginary part of a complex number where (j=mathrm{1,2},ldots ,r).

2. 2.

   The eigenfunctions (psi ={e}^{varOmega t}) oscillating (with or without growth or decay) at a frequency ({f}_{j}).

3. 3.

   The spatial modes (phi) representing the weights that reconstruct the data from the (r) eigenfunctions.

Therefore,

where (b) is computed from initial conditions ({x}_{0}) such that ({x}_{0}=phi b).

DMD fails if (m,gg, n) since (n) modes are insufficient to accurately approximate the time dynamics in (m) samples⁷². Hankel time delay embeddings can overcome this deficiency by augmenting the data by (h) measurements with past history⁷³. To guarantee an optimal representation, we used (h) time delay embeddings satisfying ({hn}) > (2m)¹⁸. For more detail about DMD and Hankel time delay embeddings, refer to Supplementary Note 2. DMD applied on a sample time-window is described in Supplementary Note 3.

In our study, we first filtered the 5-minute-long interictal iEEG segments in ({sb}) and ({rb}) (Fig. 1b). The segments were dissected into (L) (({L}_{1}) in sb and ({L}_{2}) in rb) time-windows with (m) samples each using a sliding window with 95% overlap. We used 250 ms time-windows for ({sb}) and 150 ms time-windows for ({rb}) since they were sufficient to accurately capture the characteristic waveforms of IEDs and ripples^74,75. Since a 250 ms time-window spans only one cycle of (delta) activity, two cycles of (theta) activity, and three cycles of (alpha) activity, estimating the signal power reliably with a brief time-window is challenging mostly due to the windowing edge effect. Yet, unlike the discrete Fourier and continuous wavelet transforms, which are constrained by the window size for capturing low frequency components, the DMD is not subject to these limitations¹⁸. Additional analysis to validate the ability of DMD to accurately estimate the frequency components from 250 ms time-windows are provided in Supplementary Note 4.

Each time-window ({X}_{i}) ((i=mathrm{1,2},ldots ,L)) was processed using DMD with ({r}_{1}) modes in ({sb}) and ({r}_{2}) modes in ({rb}) to extract spatial mode matrices (({phi }_{i}^{{sb}}in {{mathbb{C}}}^{ntimes {r}_{1}},{phi }_{i}^{{rb}}in {{mathbb{C}}}^{ntimes {r}_{2}})) and eigenfunctions (({psi }_{i}^{{sb}}in {{mathbb{C}}}^{{r}_{1}times m},{psi }_{i}^{{rb}}in {{mathbb{C}}}^{{r}_{2}times m})) with their corresponding oscillation frequencies (({f}_{i}^{{sb}}{{mathbb{in }}{mathbb{R}}}^{{r}_{1}},{f}_{i}^{{rb}}{{mathbb{in }}{mathbb{R}}}^{{r}_{2}})). We justified and used ({r}_{1}=50) and ({r}_{2}=100) (Supplementary Note 5).

The magnitude of the spatial mode matrices ((left|{phi }_{i}^{{sb}}right|,left|{(phi }_{i}^{{rb}}right|)) in each time-window (i) at frequencies (({f}_{i}^{{sb}},{f}_{i}^{{rb}})) are computed to generate DMD spectra of each channel (Fig. 1b). We define the following physiologically relevant frequency bands: delta ((delta =) 1–4 Hz), theta ((theta =,)[4–8] Hz), alpha ((alpha =)[8–12] Hz, beta ((beta =,)[12–30] Hz), gamma (γ = ([)30–80] Hz), and broadband (({sb}=)[1–80] Hz). For ({sb}), DMD spectra powers (left|{phi }_{i}^{{sb}}right|) are distributed in the six defined bands by collecting the average power of the modes whose frequency of oscillation ({f}_{i,j}^{{sb}}) ((j=mathrm{1,2},ldots ,{r}_{1})) falls within the band boundaries. For ({rb}), the average spectrum is computed using all ({r}_{2}) components. We denote by (Y) the feature matrices that contain the DMD spectra of all the time-windows in both ({sb}) and ({rb}) (Fig. 1c).

Extraction of the entire network

The entire network was computed by averaging the feature matrix (Y) across time-windows for each patient in each band. These networks represent average DMD power of channels across the entire iEEG segment serving as a baseline to compare with the proposed approach. The resulting vectors ((in {{mathbb{R}}}^{ntimes 1})) are thresholded using the mean plus one standard deviation as threshold.

Dominant brain network extraction using NNMF

The next step of our processing pipeline involved the identification of the consistently active interictal networks that occur recurrently across time-windows. We used NNMF, a dimensionality reduction and an unsupervised ML method, to extract these dominant spatial configurations and their corresponding temporal activity³⁹.

NNMF decomposes a non-negative matrix (Y)((in {{mathbb{R}}}^{ntimes L})) into two non-negative matrices (W) ((in {{mathbb{R}}}^{ntimes k})) and (H)((in {{mathbb{R}}}^{ktimes L})) such that (Yapprox Wtimes H) where (k) is an integer (< min {n,m}) that defines the number of spatial basis functions to extract. NNMF estimates (W) and (H) by minimizing the reconstruction error and constraining both (W) and (H) to be positive:

Typically, (k) is a user-defined parameter. (W) is termed the basis matrix and contains the (k)-dimensional representation of the data, while (H) is the coefficient matrix which contains reconstruction coefficients that reconstructs the original data (Y) from the (k) basis functions in (W)³⁹.

In our study, we used NNMF to process the feature data matrix (Y) in each band separately to extract (k) dominant spatial basis functions ((in {{mathbb{R}}}^{ntimes 1})) that represent the recurrent spatial configurations across the time-windows in the form of brain networks (columns of (W)) and their corresponding temporal activity strengths (corresponding rows in (H)). Supplementary Note 6 describes NNMF on sample data. We assumed two distinct active networks (i.e., interictal epileptogenic and background) in interictal segments and set (k=2). We justified our choice using the dispersion coefficient method⁷⁶ in Supplementary Note 7.

Determining activity of brain networks across time

By assuming that one network is active in each time-window, we quantized the coefficient matrix (H) per column by finding the index (q) with maximum coefficient. In other words, for each time-window, where (qin left{mathrm{1,2}right}), we chose the index of the maximum per column:

The vector (V) forms the temporal map, showing the indices of the active networks across time-windows. The basis functions (columns of (W)) are thresholded using the mean plus one standard deviation per column, forming the two networks (Fig. 1d) (Supplementary Note 6).

Network categorization

After identifying the two networks and their corresponding temporal maps, we categorized them as epileptogenic and background. We presumed that the IEN contaminates the background activity intermittently across time. Therefore, we considered the network with the least frequent activity in the temporal map to be the epileptogenic (Fig. 1e).

Generating stable networks and temporal maps

Since NNMF factorizes the matrices starting from random initial values, for each patient and band we repeated the following three steps 30 times: NNMF, identification of networks and temporal maps, and categorization of networks as epileptogenic and background. We used k-means clustering⁷⁷ to smooth the indices of the temporal map (V). Since k-means assigns the cluster numbers randomly, we matched the indices of (V) with the ones generated with k-means (Supplementary Note 8). We then averaged the resulting 30 IENs and background networks, as well as the temporal maps.

Network properties

We estimated and compared DMD power of the entire, interictal epileptogenic, and background networks inside and outside the resection and SOZ for good- and poor-outcome patients, separately (Wilcoxon signed-rank test). We then computed three measures to characterize the networks in reference to the resection volume. For each of the three networks (i.e., entire, interictal epileptogenic, and background), we estimated the F_net, O_res, and D_res (Fig. 1d). We defined F_net as the normalized reciprocal of the average Euclidean distance between electrodes of the thresholded networks. We defined O_res as the percentage of electrodes of the thresholded network within 10 mm from resection. Finally, we defined D_res as the average Euclidean distance of the thresholded electrodes in a network to resection. For each band, we compared F_net, O_res, and D_res of the entire, interictal epileptogenic, and background networks (Wilcoxon signed-rank test). The selection of the 10 mm cut-off was based on studies that showed that the gyral width is between 11 to 21 mm⁷⁸.

Test for consistency and robustness of extracted networks

To test the robustness of our methodology, we performed five-fold cross-validation. For each patient, we dissected the 5-minute-long iEEG segment into five 1-minute-long segments with no overlap forming five folds. For each fold, we estimated the IEN and the background networks from the 1-minute segment and the remaining 4-minute segment. For each fold, we then computed the Dice⁴⁰ score between the networks (IEN and background) identified from the 1-minute segment and the ones identified using the remaining 4-minute segments. We then computed the average Dice score across folds. Refer to Supplementary Figure 3 for an illustration of the methodology. The analysis evaluates the consistency of the identified networks when different segments of the same data are analyzed. In general, a Dice score (>)0.8 is considered to have almost perfect agreement⁴¹. The average Dice score of the five-fold tests were computed for both the IEN and background network.

As a form of surrogate testing⁷⁹, for each band and patient, we first generated 100 amplitude adjusted phase shuffled surrogates of the feature matrix derived from the 5-minute iEEG. This procedure keeps the power information while randomizing the temporal relationships by shuffling the order of the values across time-windows. The randomized feature matrices are processed using NNMF to extract the IENs and the background networks, which are termed here as random networks. The methodology can help verify that the identified networks are not due to random fluctuations in the data, but rather due to patterns that repeat consistently across time. We then computed the average Dice scores of the random networks with the IEN and the background networks estimated from the 5-minute-long data segments.

Predicting the EZ

To assess the ability of the entire, interictal epileptogenic, and background networks to predict the EZ, we used the power of the networks as predictors and the resection as target. We performed the following analysis in each band separately and only for good-outcome patients since their resected tissue is assumed to contain the actual EZ. For each network, we considered electrodes whose power was above the threshold as active and the rest as inactive. Electrodes ≤10 mm from resection were considered to be inside resection. We then defined as: (i) true positives (TP), the number of active electrodes that were inside resection; (ii) false positives (FP), the number of inactive electrodes that were inside resection; (iii) false negatives (FN), the number of active electrodes that were not inside resection; and (iv) true negatives (TN), the number of inactive electrodes that were not inside resection. We then calculated the following performance metrics per patient in each band: (i) sensitivity [TP/(TP + FN)]; (ii) specificity [TN/(TN + FP)]; (iii) precision or PPV [TP/(TP + FP)]; (iv) NPV [TN/(TN + FN)]; (v) accuracy [(TP + TN)/(TP + TN + FP + FN)]; and (vi) AUC. The terms precision and PPV are used interchangeably. We also determined the optimal threshold in each band by averaging the thresholds determined by the maximum Youden index across patients. We used this threshold for optimal thresholding of networks. We also performed ROC AUC analysis using the SOZ as the target. We then compared the AUCs of the entire, interictal epileptogenic, and background networks in predicting the resection and SOZ in different bands (Wilcoxon signed-rank test).

We then estimated F_net, O_res, and D_res of the IENs of good- and poor-outcome patients in each band. We finally compared these properties between patients with different outcomes (Wilcoxon rank-sum test). Using outcome as the actual class and each of the three IEN properties (F_net, O_res, and D_res) as predictors, we performed ROC analysis to find for each band the Youden index and the corresponding threshold for each of the three properties. We then averaged the thresholds across all bands to estimate thresholds for: (i) ({F}_{{net}}) ((t{h}_{F})); (ii) O_res ((t{h}_{O})); and (iii) D_res ((t{h}_{D})) that discriminate good- from poor-outcome patients. We then evaluated the variation of the three IEN properties (F_net, O_res, D_res) at the patient level in the seven bands. For each patient, we also reported whether their IEN properties were above ((t{h}_{F},t{h}_{O})) or below ((t{h}_{D})) the thresholds.

Robustness of IEN across variable IED rates and implantation types

To validate the applicability of our framework to segments with sparse or absent IEDs, we investigated the ability of the IENs to predict resection and SOZ for both segments with frequent and sparse IEDs. Specifically, we selected for each patient two data segments of 1-min duration with interictal activity: (i) one with frequent IEDs (≥20 IEDs per minute), and (ii) one with sparse or absent IEDs (<20 IEDs per minute). We estimated the frequency of IEDs for each data segment and tabulated the findings in Supplementary Table 1. Considering only data from good-outcome patients, we analyzed data segments with both frequent and sparse IEDs. Only 16 out of the 27 good-outcome patients had segments of 1-minute duration with both frequent and sparse IEDs. For these patients, we estimated the IENs and calculated the AUC of the IENs with the resection and SOZ. We then compared the AUC values derived from the segments with frequent IEDs with those derived from the segments with sparse IEDs.

To examine whether our framework provides consistent findings among different implantation types, we studied the ability of the IENs to predict the resection and SOZ in patients with different implantation types in our cohort. For good-outcome patients, we initially estimated the IENs and then calculated their AUC with the resection and SOZ. We then compared the AUC of IENs with the resection and SOZ among the three implantation types.

Predicting surgical outcome

We developed automated outcome predictors with linear SVM⁴² using the three IEN properties (F_net, O_res, and D_res) as features and outcome as target. We set the regularization parameter of the SVM to 1. An SVM was trained in each band separately (i.e., the properties of IENs in each band as predictors and outcome as target). We also trained another linear SVM (denoted by SVM-all) using 21 IEN properties (seven bands (times) three properties per band) simultaneously as features that incorporated information from all bands. Each SVM was validated using five-fold cross-validation by dividing patients into five random splits. We defined as: (i) TP, the number of patients predicted as good outcome who had a good outcome; (ii) FN, the number of patients predicted as poor outcome who had a good outcome; (iii) FP, the number of patients predicted as good outcome who had a poor outcome; and (iv) TN, the number of patients predicted as poor outcome who had a poor outcome. We then calculated the following performance metrics: sensitivity, specificity, PPV, NPV, accuracy, AUC, and Fisher’s exact test p values of each classifier. Similarly, we evaluated the ability of the entire network to predict outcome, and its performance was compared to that of the IEN.

Finally, ANOVA was performed to rank the IEN properties in the seven bands in decreasing order of importance to predict outcome. Feature importance is defined as:

where P represents the p values of ANOVA. The analysis identified the IEN properties and bands that can best discriminate between good and poor outcome.

Statistical analysis

The Kolmogorov-Smirnov test was used to test the normality of the power and network properties (F_net, O_res, and D_res). Cliff’s d measure was used to compute effect sizes. We applied two-sided non-parametric Wilcoxon signed-rank test for all paired comparisons (median values inside vs. outside the resection and SOZ, comparing the properties of the entire, interictal epileptogenic, and background networks). We applied two-sided Wilcoxon rank-sum test for non-paired comparisons between good- and poor-outcome patients. Bonferroni correction was applied in all multiple comparison tests⁸⁰. We used one-sided Fisher’s exact test to evaluate the predictive value of the IENs to predict outcome. Chi-squared test was used to compare different measures (gender, implantation side, epilepsy localization, and MRI findings) in terms of outcome. We assumed statistically significant results if p (le)0.05 and reported measures as median and interquartile range. All analysis was performed using MATLAB 2022b (The MathWorks, Inc.).