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

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

Introduction

For patients with focal drug resistant epilepsy (DRE), neurosurgery is the best available treatment and achieves seizure freedom in up to 60% of cases1. The neurosurgical outcome depends upon precise and accurate delineation of the epileptogenic zone (EZ), the brain area that is indispensable for the generation of seizures1,2. Yet, the EZ is a theoretical concept since its elucidation before resection is uncertain3. Even after a successful resection, the EZ cannot be precisely defined since the resected area can also be larger than its actual extent4. 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 nonconcordant5. The seizure onset zone (SOZ), the brain area where seizures initiate, currently serves as the best approximator of the EZ1. This area is defined by analyzing ictal intervals of intracranial electroencephalography (iEEG) data1. Yet, removal of the SOZ does not always predict surgical outcome6,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 resources8. An alternate interictal biomarker of the EZ is therefore of paramount importance9.

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 ripples10,11. These abnormalities currently serve as established interictal biomarkers of epilepsy12. However, they lack specificity in delineating the EZ since they can include areas larger than the EZ, which should be preserved during surgery13. Moreover, IEDs and ripples are identified through visual inspection of iEEG data or with the aid of automated software14. Yet, both approaches are prone to errors and biases15. 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 patterns16,17,18,19, or propagate across contacts20,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 EZ28. Furthermore, IEDs are absent in ~10% of patients with DRE29. 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 activity30.

Numerous artificial intelligence (AI) tools have recently been proposed to facilitate the presurgical evaluation of DRE patients31,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 features30,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 data18,35,36,37. These methods can identify relevant epileptiform patterns and inherent “structures” in the data without explicitly training the model37. 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.

Results

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 cm3 (range: 11.5–36.9) in good- and 20.3 cm3 (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.

Table 1 Patient demographic information categorized by post-surgical outcome
Full size table

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).

Fig. 1: Overall processing pipeline.
Machine learning on interictal intracranial EEG predicts surgical outcome in drug resistant epilepsy

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 ripples21,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 time18. Since IEDs and ripples occur coherently in multiple regions16,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 time18. 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)39. Since IEDs and ripples occur transiently, alternating with background activity10,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 score40, a measure that quantifies the similarity and overlap between networks. A Dice score (ge)0.8 is considered to have almost perfect agreement41. 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).

Fig. 2: Median network power in good- and poor-outcome patients.
figure 2

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 25th and 75th 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 (Fnet), overlap with resection (Ores), and distance from resection (Dres) (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 Fnet (in (theta), (beta), (gamma), and sb; p < 0.05, d ≥ 0.38), higher Ores (in all bands; p < 0.01, d ≥ 0.40), and lower Dres (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 Fnet (in all bands except (delta); p < 0.05, d ≥ 0.38), higher Ores (in all bands except (delta) and (alpha); p < 0.05, d ≥ 0.26), and lower Dres (in all bands except (delta) and (alpha); p < 0.01, d ≥ 0.32) (Fig. 3). Finally, the entire network had higher Ores (in (delta), (theta), and (alpha); p < 0.01, d ≥ 0.27) and lower Dres (in (theta) and (alpha); p < 0.01, d ≥ 0.06) compared to the background (Fig. 3). There were no differences in Fnet between the entire and background networks in all bands (Fig. 3).

Fig. 3: Network properties in good-outcome patients.
figure 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 25th and 75th 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).

Fig. 4: Network performance metrics for predicting resection.
figure 4

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 25th and 75th 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).

Fig. 5: Prediction of resection and seizure onset zone (SOZ).
figure 5

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 25th and 75th 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).

Fig. 6: Consistency of interictal epileptogenic networks (IENs) across variable interictal epileptic discharge (IED) rates and implantation types.
figure 6

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 25th and 75th 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 Fnet, lower Ores, and lower Dres to resection in good- compared to poor-outcome patients. At the population level, the IEN had higher Fnet 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 Ores 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 Dres 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).

Fig. 7: Interictal Epileptogenic network (IEN) properties.
figure 7

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 25th and 75th 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 Fnet and Ores were higher in most patients with good (compared to poor) outcome (Fig. 7b). Similarly, Dres was lower in most patients with good outcome. We report Fnet, Ores, and Dres 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 Fnet(t{h}_{F}) in (theta) band (Fig. 7b); contrarily, only 4 out of 16 (25%) poor-outcome patients had Fnet(t{h}_{F}). We also found that 85.1% of good-outcome patients (23 out of 27) had Ores(t{h}_{O}) with resection in (theta) band (Fig. 7b); contrarily, only 31.2% (5 out of 16) of poor-outcome patients had Ores ≥ (t{h}_{O}) with resection. Finally, we found that 92.6% of good-outcome patients (25 out of 27) had Dres ≤ (t{h}_{D}) in (theta) band; contrarily, only 37.5% of poor-outcome patients (6 out of 16) had Dres ≤ (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)42 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).

Fig. 8: Outcome prediction using machine learning.
figure 8

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)43 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 Fnet, Ores, and Dres in (theta) band had the highest importance for predicting outcome followed by Fnet and Ores in (beta), and Dres 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).

Fig. 9: Case studies of a good- and a poor-outcome patient.
figure 9

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 Fnet (0.74 >(t{h}_{F})), higher Ores (100% >(t{h}_{O})), and lower Dres (6 mm <(t{h}_{D})) for the IEN (compared to background; Fnet (=)0.51, Ores (=)18%, and Dres=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 Fnet <(t{h}_{F}), Ores <(t{h}_{O}), and Dres >(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.

Discussion

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 EZ16,21,25,30,44,45,46. However, these biomarkers rely on manual selection of IEDs by experts, which is time-consuming and prone to errors15. 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 errors48. 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 data49. Our framework overcomes these hurdles by identifying coherent IENs without the need to detect IEDs or ripples. It adopts an unsupervised ML approach18,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 methods17,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 outcome16,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 predictor16,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 ms11,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 controls53,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 mechanism55. 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) band53,56.

Automated IED detectors4,15, including AI-based tools57,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 activity18,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 consistency18,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 patient2,26,27,60. In particular, sEEG employs a network-based implantation strategy that differs fundamentally from subdural and depth electrodes26,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 level61. 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) band62, 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 sEEG21,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.

Methods

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 seizures64. This selection ensured the inclusion of segments with the highest IED rate and minimal motion artifacts65,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)67. We adjusted the electrode locations to compensate for possible brain shift occurring after electrocorticography implantation68.

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)67. 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)38. 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 DMD69. 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 decay70. 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 coherence18,71.

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

$${x}_{t+1}=A{x}_{t}+{w}_{t},t=0,1,2,…,m-1$$
(1)

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:

$${{||}{x}_{t+1}-A{x}_{t}{||}}_{2}$$
(2)

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

$$X=left(begin{array}{ccc}| & | & |\ {x}_{0} &;;;;; {x}_{1}ldots &;;;;;{x}_{m-1}\ | & | & |end{array}right)$$
(3)

(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):

$$X{prime} =left(begin{array}{ccc}| & | & |\ {x}_{1} &;;;;;;{x}_{2}ldots &;;;{x}_{m}\ | & | & |end{array}right)$$
(4)

Equation (1) can be written as:

$${X}^{prime} ={AX}$$
(5)

The solution of the least squares problem approximates A:

$$A={X}^{{prime} }{X}^{dagger }$$
(6)

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,

$$Xapprox phi {e}^{varOmega t}b$$
(7)

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) samples72. Hankel time delay embeddings can overcome this deficiency by augmenting the data by (h) measurements with past history73. To guarantee an optimal representation, we used (h) time delay embeddings satisfying ({hn}) > (2m)18. 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 ripples74,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 limitations18. 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 activity39.

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:

$$underbrace{{minimize}}_{W,H}{{||}Y-{WH}{||}}_{2}{subject}; {to}W,>,0,{H},>,0$$
(8)

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)39.

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 method76 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:

$$V={arg }mathop{max }limits_{q}({H}_{j})$$
(9)

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 clustering77 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 Fnet, Ores, and Dres (Fig. 1d). We defined Fnet as the normalized reciprocal of the average Euclidean distance between electrodes of the thresholded networks. We defined Ores as the percentage of electrodes of the thresholded network within 10 mm from resection. Finally, we defined Dres as the average Euclidean distance of the thresholded electrodes in a network to resection. For each band, we compared Fnet, Ores, and Dres 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 mm78.

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 Dice40 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 agreement41. The average Dice score of the five-fold tests were computed for both the IEN and background network.

As a form of surrogate testing79, 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 Fnet, Ores, and Dres 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 (Fnet, Ores, and Dres) 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) Ores ((t{h}_{O})); and (iii) Dres ((t{h}_{D})) that discriminate good- from poor-outcome patients. We then evaluated the variation of the three IEN properties (Fnet, Ores, Dres) 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 SVM42 using the three IEN properties (Fnet, Ores, and Dres) 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:

$$I=-log (P)$$
(10)

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 (Fnet, Ores, and Dres). 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 tests80. 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.).

Related Articles

Rapid brain tumor classification from sparse epigenomic data

Although the intraoperative molecular diagnosis of the approximately 100 known brain tumor entities described to date has been a goal of neuropathology for the past decade, achieving this within a clinically relevant timeframe of under 1 h after biopsy collection remains elusive. Advances in third-generation sequencing have brought this goal closer, but established machine learning techniques rely on computationally intensive methods, making them impractical for live diagnostic workflows in clinical applications. Here we present MethyLYZR, a naive Bayesian framework enabling fully tractable, live classification of cancer epigenomes. For evaluation, we used nanopore sequencing to classify over 200 brain tumor samples, including 10 sequenced in a clinical setting next to the operating room, achieving highly accurate results within 15 min of sequencing. MethyLYZR can be run in parallel with an ongoing nanopore experiment with negligible computational overhead. Therefore, the only limiting factors for even faster time to results are DNA extraction time and the nanopore sequencer’s maximum parallel throughput. Although more evidence from prospective studies is needed, our study suggests the potential applicability of MethyLYZR for live molecular classification of nervous system malignancies using nanopore sequencing not only for the neurosurgical intraoperative use case but also for other oncologic indications and the classification of tumors from cell-free DNA in liquid biopsies.

EEG-based headset sleep wearable devices

The rise of wearable technology has led to EEG-based sleep monitoring devices that use electrodes placed on the forehead, ear, or neck. These devices offer promising applications in clinical and healthy populations by comparing sleep patterns, monitoring intervention responses, and examining the relationship between sleep and lifestyle factors. Despite their potential, challenges like validation against polysomnography, regulatory hurdles, data privacy, and usability hinder clinical adoption. This review explores these devices, their applications, and integration challenges in clinical practice.

Chemogenetics with PSAM4-GlyR decreases excitability and epileptiform activity in epileptic hippocampus

Despite the availability of new drugs on the clinics in recent years, drug-resistant epilepsy remains an unresolved challenge for healthcare, and one-third of epilepsy patients remain refractory to anti-seizure medications. Gene therapy in experimental models has emerged as effective treatment targeting specific neuronal populations in the epileptogenic focus. When combined with an external chemical activator using chemogenetics, it also becomes an “on-demand” treatment. Here, we evaluate a targeted and specific chemogenetic therapy, the PSAM/PSEM system, which holds promise as a potential candidate for clinical application in treating drug-resistant epilepsy. We show that the inert ligand uPSEM817, which selectively activates the chloride-permeable channel PSAM4-GlyR, effectively reduces the number of depolarization-induced action potentials in vitro. This effect is likely due to the shunting of depolarizing currents, as evidenced by decreased membrane resistance in these cells. In organotypic slices, uPSEM817 decreased the number of bursts and peak amplitude of events of spontaneous epileptiform activity. Although administration of uPSEM817 in vivo did not significantly alter electrographic seizures in a male mouse model of temporal lobe epilepsy, it did demonstrate a strong trend toward reducing the frequency of interictal epileptiform discharges. These findings indicate that PSAM4-GlyR-based chemogenetics holds potential as an anti-seizure strategy, although further refinement is necessary to enhance its efficacy.

Clinical practice recommendations for the diagnosis and management of X-linked hypophosphataemia

X-linked hypophosphataemia (XLH) is a rare metabolic bone disorder caused by pathogenic variants in the PHEX gene, which is predominantly expressed in osteoblasts, osteocytes and odontoblasts. XLH is characterized by increased synthesis of the bone-derived phosphaturic hormone fibroblast growth factor 23 (FGF23), which results in renal phosphate wasting with consecutive hypophosphataemia, rickets, osteomalacia, disproportionate short stature, oral manifestations, pseudofractures, craniosynostosis, enthesopathies and osteoarthritis. Patients with XLH should be provided with multidisciplinary care organized by a metabolic bone expert. Historically, these patients were treated with frequent doses of oral phosphate supplements and active vitamin D, which was of limited efficiency and associated with adverse effects. However, the management of XLH has evolved in the past few years owing to the availability of burosumab, a fully humanized monoclonal antibody that neutralizes circulating FGF23. Here, we provide updated clinical practice recommendations for the diagnosis and management of XLH to improve outcomes and quality of life in these patients.

Understanding learning through uncertainty and bias

Learning allows humans and other animals to make predictions about the environment that facilitate adaptive behavior. Casting learning as predictive inference can shed light on normative cognitive mechanisms that improve predictions under uncertainty. Drawing on normative learning models, we illustrate how learning should be adjusted to different sources of uncertainty, including perceptual uncertainty, risk, and uncertainty due to environmental changes. Such models explain many hallmarks of human learning in terms of specific statistical considerations that come into play when updating predictions under uncertainty. However, humans also display systematic learning biases that deviate from normative models, as studied in computational psychiatry. Some biases can be explained as normative inference conditioned on inaccurate prior assumptions about the environment, while others reflect approximations to Bayesian inference aimed at reducing cognitive demands. These biases offer insights into cognitive mechanisms underlying learning and how they might go awry in psychiatric illness.

Responses

Your email address will not be published. Required fields are marked *