Abstract
Epilepsy, a neurological disorder characterized by recurrent unprovoked seizures, significantly impacts patient quality of life. Current classification methods focus primarily on clinical observations and electroencephalography (EEG) analysis, often overlooking the underlying dynamics driving seizures. This study uses surface EEG data to identify seizure transitions using a dynamical systems–based framework—the taxonomy of seizure dynamotypes—previously examined only in invasive data. We applied principal component and independent component (IC) analysis to surface EEG recordings from 1,177 seizures in 158 patients with focal epilepsy, decomposing the signals into ICs. The ICs were visually labeled for clear seizure transitions and bifurcation morphologies (BifMs), which were then examined using Bayesian multilevel modeling in the context of clinical factors. Our analysis reveals that certain onset bifurcations (saddle node on invariant circle and supercritical Hopf) are more prevalent during wakefulness compared with their reduced rate during nonrapid eye movement (NREM) sleep, particularly NREM3. We discuss the possible implications of our results in the context of modeling approaches and suggest additional avenues to continue this exploration. Furthermore, we demonstrate the feasibility of automating this classification process using machine learning, achieving high performance in identifying seizure-related ICs and classifying interspike interval changes. Our findings suggest that the noise in surface EEG may obscure certain BifMs, and we suggest technical improvements that could enhance detection accuracy. Expanding the dataset and incorporating long-term biological rhythms, such as circadian and multiday cycles, may provide a more comprehensive understanding of seizure dynamics and improve clinical decision-making.
- automated seizure classification using machine learning
- bifurcations in epilepsy
- epileptic seizure classification
- noninvasive EEG analysis
- sleep stages and seizure dynamics
- taxonomy of seizure dynamotypes (TSD)
Significance Statement
Traditional seizure classification focuses on clinical symptoms and electrophysiological signs but often overlooks the underlying seizure dynamics. The taxonomy of seizure dynamotypes introduces a novel computational approach that links electrophysiological transition signatures to these dynamics. While previously applied to invasive recordings, this study extends the taxonomy to noninvasive electroencephalography. Our analysis reveals a relationship between sleep stages and seizure dynamics. We suggest that integrating these modeling approaches with sleep and circadian dynamics models may reveal insights into seizure timing and generalization, opening new pathways for better diagnostics. Broader adoption of this classification is limited by its labor-intensive visual inspection process. Here, we demonstrate the potential of automated classification, enabling analysis to scale to larger cohorts.
Introduction
Epilepsy, characterized by recurrent unprovoked seizures, affects millions worldwide and presents a significant challenge in neurology. Seizure classification approaches rely primarily on clinical observations and electroencephalography (EEG) visual analysis. Yet, in most cases, this classification does not provide information on the underlying dynamics. A better understanding of the transitions into and out of a seizure can shed light on the mechanisms that drive seizure activity.
From a dynamical systems point of view, the brain is constantly shifting between dynamical states in response to internal fluctuations and external stimuli (Deco et al., 2011). Epilepsy is viewed as a condition that facilitates transient oscillatory ictal states (Da Lopes Silva et al., 2003), and the transition into and out of a seizure can be considered a bifurcation—a mathematical concept that describes a sudden change in the system's dynamics. This offers a new lens through which to analyze the initiation and termination of seizures, potentially leading to novel strategies for their treatment.
Through rigorous mathematical and modeling work, research in this area has mapped out the “taxonomy of seizure dynamotypes” (TSD), introducing a structured classification approach that is based on a minimal descriptive model, connecting electrophysiological morphologies to dynamical bifurcations (Saggio et al., 2020). In this paper, we will refer to these as bifurcation morphologies (BifMs). Similar models have demonstrated clinical relevance by improving the outcome of resective surgery (An et al., 2019; Makhalova et al., 2022) or predicting the effect of electrical stimulation on seizure abortion (Szuromi et al., 2023).
The TSD approach describes seizure transitions using several time scales: a fast subsystem that intermittently generates seizure activity; a slow subsystem that modulates the fast system's parameters, facilitating mainly seizure termination; and an ultraslow subsystem that can effect seizure propensity (Saggio et al., 2020). From this modeling perspective, the fast subsystem represents the seizure onset zone (SOZ); the slow subsystem represents reactive processes limiting seizure duration; and the ultraslow subsystem corresponds to changes in the brain state that occur over significantly longer timescales. The last, which is not explicitly modeled in the TSD framework, may include sleep stages (Herman et al., 2001), as well as circadian or longer physiological rhythms, which have been shown to affect seizure propensity (Karoly et al., 2021). Sleep stages, particularly wake, nonrapid eye movement (NREM), and REM, represent inherently different dynamical states with different tendencies toward epileptic events (Ng and Pavlova, 2013; Halász et al., 2019), making them a significant target for investigation.
The TSD framework has been investigated using invasive data from animal models and human subjects (Jirsa et al., 2014; Saggio et al., 2017, 2020). Previous research explored the relationship between factors such as age, gender, etiology, and localization (Saggio et al., 2020); however, this analysis was limited to a single seizure per patient and did not account for sleep stages. Analyzing multiple seizures per patient presents several statistical challenges, including accounting for repeated measures within patients and considering both patient- and seizure-level factors. Additionally, the inherent imbalance in BifMs introduces further statistical complexities.
As noted above, previous research has primarily relied on invasive intracranial recordings, which, while effective, are limited in broader application due to their invasiveness. In contrast, surface EEG is more commonly used in epilepsy diagnosis, but it is often compromised by significant noise from nonbrain sources, particularly during seizures. To mitigate this issue, we have employed signal decomposition methods such as principal component analysis (PCA; Maćkiewicz and Ratajczak, 1993) and blind source separation techniques like independent component analysis (ICA; Bell and Sejnowski, 1995) to extract seizure-related information from the EEG signal. These methods have proven effective in tasks like focal seizure localization (De Vos et al., 2007; Stern et al., 2009; Matin et al., 2019; Rabby et al., 2021). Notably, a strong temporal and spatial correlation between noninvasive independent components (ICs) and intracranial recordings from the SOZ has been reported (Barborica et al., 2021). Therefore, selecting ICs that contain seizure information may provide a reasonable approximation of the activity generated in the SOZ.
In this study, we explore a possible methodology for extracting seizure information from surface EEG by performing a visual analysis of the resulting ICs, identifying ICs of brain origin with clear seizure transitions, and classifying the BifMs. We then examine the relationship between BifMs and clinical factors at the patient and seizure levels using Bayesian multilevel modeling (Bürkner, 2017). Recognizing that the labor-intensive process of visually classifying BifMs limits scalability, we aimed to automate the identification of seizure-related ICs and the classification of different BifMs to enhance the method's applicability.
Materials and Methods
Data and selection criteria
We utilized the EPILEPSIAE database (Ihle et al., 2012), which provides meticulously labeled seizure onset and offset times for EEG recordings from 158 patients with drug-resistant focal epilepsy undergoing presurgical evaluation. This database includes 1,214 seizure annotations (Fig. 1, top right). For each seizure, a control segment of identical duration was selected, positioned at least 5 min away from any seizure activity within the same recording. Due to insufficient control data, 38 seizures were excluded, resulting in a total of 1,177 seizures that were visually analyzed (Fig. 1, top middle).
Analysis overview. Project overview: we utilized data from the EPILEPSIAE surface database (Ihle et al., 2012), which includes 1,214 seizures from 158 patients undergoing presurgical evaluation. Using the seizure annotations provided in this database, we automatically extracted and preprocessed the seizure data, adding a margin of 30% of the seizure length on each side. Thirty-seven seizures lacking sufficient data were excluded during preprocessing. The preprocessing pipeline is detailed in Materials and Methods, Visual analysis and labeling. Initially, we reduced the dimensionality of the EEG data using PCA to retain the top eight components. ICA was then applied to achieve signal components with strong statistical independence, resulting in 9,416 components. Two independent raters visually reviewed these components, with inter-rater agreement evaluated (Tables 1, 2). Only components with agreement were used for further analysis. During the visual analysis, we focused on identifying components with clear transitions into and out of a seizure agreed upon by both reviewers and classifying the transition's BifM based on the “TSD” (Saggio et al., 2020). These agreed-upon components served as the gold standard for developing two automated classifiers. We evaluated the models using a LOOCV approach (detailed in Materials and Methods, Development and validation of automated classification; Fig. 4; Tables 5, 6; Extended Data Tables 5-1, 6-1). A seizure could contain multiple components with clear transitions (Fig. 2A). We used objective measures extracted from all components to select the one that best represented the onset and offset for each seizure, primarily based on the transition timing and the probability of containing brain information (see Extended Data Fig. 2-2 for more details). Most seizures included a single bifurcation, but we recorded cases where more than one bifurcation type was initially labeled (detailed in Table 1). We modeled the prevalence of seizure BifMs in relation to various seizure factors (vigilance, seizure length, and classification) and patient factors (anatomical localization, age, gender, etiology) using a Bayesian multilevel model (for more details, see Materials and Methods, Modeling the interactions between BifMs and clinical factors). The overall BifM prevalence (Fig. 3) served as the prior probability for the model intercept initialization. We used the posterior estimations to examine the differences between modeled factors. Extended Data Figure 4-1 explains the metrics used to compare the posterior distributions. The analysis results are detailed in Table 4 and in Extended Data Tables 4-1 and 4-2.
Table 1-1
The table includes all the patient and seizure metadata for the 1,177 seizures included in the analysis. Download Table 1-1, XLSX file.
Table 1-2
The table includes the reviewer labels for all ICs and the final labels determined by inter-rater agreement. Download Table 1-2, XLSX file.
Seizure annotations included the sleep stage before a seizure and seizure classification. In this database, seizures were classified, according to the 1981 classification scheme (Bancaud et al., 1981), as simple partial (focal aware, FAS), complex partial (focal with impaired awareness, FIAS), and secondarily generalized (focal to bilateral tonic–clonic, FBTCS). Unclassified (UC) seizures remained the same (Scheffer et al., 2017). Patient metadata included age, gender, seizure localization, and etiology. Patient ages ranged from 13 to 65 years, with a mean of 35.4 years; there were 82 male and 76 female patients. The most common seizure localization was in the temporal lobe (117), followed by the frontal lobe (17). The most common etiology was hippocampal sclerosis (54), followed by cryptogenic (46) and cortical dysplasia (28). Extended Data Table 1-1 provides detailed patient and seizure metadata used in this work.
Preprocessing for the extraction of seizure information from surface EEG data
Preprocessing was conducted using MATLAB and EEGLAB (Delorme and Makeig, 2004). The code in this work was run on a 64 bit operating system, an x64-based processor, and a Windows 10 pro operating system to obtain the results reported here. The EEG signals were downsampled to 256 Hz to ensure uniformity across the database. A high-pass filter with a 1 Hz cutoff was applied to remove low-frequency noise and drift. We utilized artifact subspace reconstruction (Mullen et al., 2015) to identify and correct strong movement artifacts in the recordings. After artifact correction, the signals were re-referenced to the mean, and a low-pass filter with a 40 Hz cutoff was applied.
Under the assumption that synchronous ictal EEG activity accounts for a significant portion of signal variance, we applied PCA (Maćkiewicz and Ratajczak, 1993). PCA transforms the data into a new set of variables called principal components (PCs), which are orthogonal (uncorrelated) and ordered by the amount of variance they capture. By selecting the top PCs, we reduced the dimensionality of the data while preserving most of the variance.
EEG data often exhibit significant amplitude variations due to differences in electrode impedance, contact, and positioning, which can bias the PCA results. For example, the first PC may be mainly influenced by the electrode with the largest amplitude. This could be due to either synchronized ictal activity or technical reasons (e.g., movement of an electrode). To address this issue, we standardized the electrode data with a 30% margin on each side of the seizure annotation. Including these edges in the standardization process (mean reduction and unit scaling, also termed pre-whitening) preserves the relative amplitude increase associated with ictal activity, which might otherwise be diminished if only ictal data were used.
In the next step, we applied infomax ICA on the resulting PCs (Bell and Sejnowski, 1995). ICA is used to find a linear transformation that makes the output signals as statistically independent from each other as possible. It is commonly used in EEG analysis to separate and remove noise sources such as eye movements, muscle activity, and heartbeats (Delorme and Makeig, 2004). Previous research has demonstrated that ICA can be effectively applied to surface EEG data to extract ICs that exhibit strong temporal and spatial correlations with ictal intracranial recordings from the SOZ (Barborica et al., 2021).
The PCA and ICA steps produce two linear transformation matrices that convert the prewhitened EEG electrode data into ICs. By multiplying these matrices, we can present the electrode weights for each component, which helps identify the relative contribution of surface electrodes to each IC, as shown in Figure 2B. A preliminary analysis of 10 randomly selected seizures indicated that seizure information was concentrated within the first eight PCs. Using more than eight PCs introduced noise and resulted in less distinct ICs with clear seizure information. Therefore, we established a cutoff of eight PCs to minimize noise and reduce the workload for visual labeling, leading to a total of 9,416 components (Fig. 1, top middle).
Information used to label the ICs. A, The time series of the ICs from an example seizure. The ICs highlighted in dark blue contain a clear transition into a seizure that includes data mainly from brain sources. The PCA–ICA procedure affects the components’ scale; hence, the units are presented as procedure-defined units (PDUs). B, For the selection of the ICs, we used custom code and EEGLAB visualization functions to present each reviewer with the component time series, the topographic distribution of the electrode weights, the spectral profile, and the IClabel classifier output. The topographic maps represent the normalized weight (NW) distribution used to transform the EEG time series into the IC time series. These maps help evaluate the scalp distribution of each component, helping with the recognition of patterns typical for brain-related sources (IC1, 3, 5), eye movement (IC2, 6), or other nonbrain sources (IC4, 7, 8; as described in Pion-Tonachini et al., 2019b). Using the power spectrum, we could better distinguish ICs with typical brain frequency profiles (IC1, 3, 5) from components containing mainly muscular activity (IC4, 6, 8). These and additional signal features are input to the IClabel classifier. This classifier is built into EEGLAB and is designed to identify and remove noise components in surface EEG preprocessing (Pion-Tonachini et al., 2019a). To support our decision and to estimate the consistency of this classifier on seizure data, we presented the classification probabilities to the user. Extended Data Figure 2-2A displays the overall brain score distribution given the different labels and shows that clear transitions are generally associated with a higher brain score. However, exceptions exist, and this does not provide a clear-cut role to identify seizure data. C, Examples of the three onset BifMs: (a) SupH, which shows an amplitude increase (from zero); (b) SNIC, characterized by increased spike frequency (from zero) or decreasing ISIs from the seizure onset; and (c) SN (expected DC shift) or SubH, in which both amplitude and ISI remain constant or random; however, due to the lack of DC recording, these are indistinguishable. Examples of the three offset BifMs include (d) SNIC or SH (expected DC shift), characterized by increasing ISIs or reduction in frequency (to zero), and these are indistinguishable like the previous example; (e) SupH, with an amplitude decrease (to zero); and (f) FLC, with amplitude and ISI that are constant or random toward seizure offset. See Extended Data Figure 2-1 for additional examples. D, An increasing amplitude at the seizure onset characterizing the SupH BifM may be attributed to propagation in noninvasive EEG. Therefore, it is important to analyze the component time series in the context of additional information beyond the time series alone. For example, if we observe only IC2, it shows a SupH BifM starting at the 20 s mark. However, when examining IC3, we see a lower constant amplitude beginning at the 16 s mark and undergoing a second change at the 20 s mark. The topographic maps of these components reveal a more localized weight distribution in IC3 than in IC2. This suggests that the pattern in IC2 may result from propagation, emphasizing the importance of selecting the earliest and most localized spatial pattern (in focal seizures) as the seizure onset BifM.
Figure 2-1
Labeling examples. Each component is a linear combination of all electrodes with varying weight ranges; hence, we present the component time series with procedure-defined units (PDUs). Download Figure 2-1, TIF file.
Figure 2-2
Automated measures and detectability characteristics. Download Figure 2-2, TIF file.
Visual analysis and labeling
Before initiating the visual analysis and labeling, two independent reviewers were trained on the classification tasks. The initial challenge was identifying ICs containing brain activity related to seizures in surface EEG data. To improve proficiency in distinguishing between brain and nonbrain sources, reviewers used the ICLable labeling tutorial application (Pion-Tonachini et al., 2019b). This step was crucial for enhancing the ability to identify various physiological and nonphysiological IC sources. Figure 2, A and B, illustrates the primary information sources used for labeling.
Reviewers were also trained to identify seizure onset and offset BifMs using empirical and simulation examples from previous work (Saggio et al., 2020). This classification was based on specific trends in frequency (interspike intervals, ISIs) and amplitude changes. In general, there are six types of BifMs: supercritical Hopf (SupH), saddle node on invariant circle (SNIC), saddle homoclinic (SH), saddle node (SN), subcritical Hopf (SubH), and fold limit cycle (FLC).
The SupH BifM prescribes an amplitude scaling law where the amplitude increases from zero at the onset and decreases to zero at the offset. SNIC and SH BifMs are characterized by a spike frequency scaling law. In the SNIC BifM at the onset, the frequency increases from zero. At the offset, SNIC and SH BifMs are characterized by frequency decreasing to zero. Finally, SN and SubH BifMs at the onset and the FLC BifM at the offset have no prescribed scaling laws and may exhibit constant or arbitrary amplitude and spike frequency. Detailed examples and expected features of each BifM are provided in Figure 2C.
The absence of a direct current (DC) component in our data limited the distinction between SN and SubH BifMs at the onset, as well as between SH and SNIC BifMs at the offset. As a result, we grouped these BifMs into combined categories, following the approach used in previous studies (Saggio et al., 2020). The onset labels were categorized as SupH, SNIC, SN/SubH, and unclear, and the offset labels were defined as SupH, SNIC/SH, FLC, and unclear.
We developed an in-house graphical user interface in which reviewers were presented with the EEG time series, the eight ICs’ time series for each seizure, their spatial activation maps, and the power spectrum (see Fig. 2A,B for additional details). This interface provided a context for selecting the most relevant ICs, which contained both brain activity and clear transitions into and out of seizures; these were then labeled with the relevant BifM. Each reviewer labeled the ICs from the 1,177 seizures—a total of 9,416 ICs (Fig. 1, IC level). Figure 2C provides labeling examples and considerations, with additional examples provided in Extended Data Figure 2-1.
After the first labeling session, the reviewers recognized a learning curve in distinguishing brain-derived data from other physiological sources within seizure data. Components not unanimously agreed upon as seizure-related were independently relabeled by each reviewer, who was blinded to the prior labels. Components that were unanimously identified as either seizure-related or not, despite differences in BifM labels, were not relabeled. The top three rows in Table 1 summarize the ICs labeled by both reviewers and the number of components with agreements at the onset and offset.
Component agreement numbers: this table includes details on the numbers of ICs and seizures included throughout the analysis
We quantified the inter-rater agreement using Cohen's kappa coefficient (Cohen, 1960), evaluating it separately for ICs identified as having a clear transition and for the BifMs labeled. To estimate the chance level for kappa values given the label imbalance, we conducted two permutation analyses. The first approach involved randomly permuting the reviewers’ labels overall, while the second involved permuting labels only within each seizure, ensuring the same number of ICs labeled by each reviewer per seizure. Each analysis included 10,000 permutations. Following this agreement evaluation, only ICs marked with a similar bifurcation by both reviewers were considered to be clear bifurcations and were included in further analysis.
In some instances, more than one IC per seizure was labeled as having a BifM. Typically, the ICs for each seizure shared the same BifM; however, there were several seizures (40 onsets and 53 offsets) where multiple ICs were labeled with different BifMs. In these cases, we determined the BifM at the seizure level by evaluating objective measures such as the transition time and a general brain score. For a more detailed description of these measures and the decision flow, see Extended Data Figure 2-2A–C.
Modeling the interactions between BifMs and clinical factors
To analyze the relationship between BifMs and clinical factors at the patient and the seizure levels, we used the BRMS package (Bürkner, 2017) in R, which implements Bayesian multilevel models. These statistical models allow us to account for data hierarchies and repeated measures within a patient. Moreover, Bayesian models provide a probabilistic framework for statistical inference, incorporating prior knowledge and the observed data to estimate the parameters of interest.
Patient-level effects included age, gender, the lobe of the epileptogenic zone, and etiology. We treated the patient as a random effect for vigilance, seizure length, and classification to account for repeated measures within subjects. We included all three BifMs—SN/SubH, SNIC, and SupH—in a categorical logistic regression model for the onset. However, we encountered only three seizures with a SupH BifM for the offset. Consequently, we created a binomial model that included only SH/SNIC and FLC BifMs, excluding seizures with SupH BifM from the analysis.
In Bayesian statistics, a prior represents our initial belief about the parameters. Here, the intercept priors reflect our expectations about the BifM distribution. By initializing the model intercepts with the overall bifurcation distribution, we effectively assume a particular and similar BifM distribution across the clinical factors. Consequently, more robust evidence is required to deviate from this initial distribution during the model fitting process. This approach helps prevent overestimating effects due to random variations and makes our conclusions more robust.
The model distribution parameters were estimated using Monte Carlo–Markov chain (MCMC). MCMC is a computational algorithm used to approximate the posterior distribution when it is difficult to calculate directly. After considering the observed data, the posterior distribution represents our updated beliefs about the parameters. For each model, we used six chains with 6,000 iterations each; the first 2,000 iterations were considered a warm-up (or burn-in period) and were excluded from the final posterior estimation to ensure that the chains had reached a stable distribution.
For adequate model fitting, we evaluated the fit quality in two ways. First, we performed a posterior predictive check (PP check) by simulating the overall posterior distribution and confirming that it accurately represented the observed data. This involves generating and comparing data from the model to the actual data to assess how well the model captures the underlying patterns. Second, we assessed the
The Bayesian modeling approach allowed us to estimate posterior distributions under two clinical conditions, marginalized across all other factors, and to analyze the difference distribution to estimate the direction and effect size. For a centrality measure, we report the median and maximum a posteriori (MAP) of the difference distribution function and the credible interval (CI), which is equivalent to the frequentist confidence interval and provides a range within which the parameter value lies with a 0.95 probability (Extended Data Fig. 4-1A). To estimate the magnitude of the effect, we use the probability of direction (pd) (Extended Data Fig. 4-1B). The pd is mathematically defined as the proportion of the posterior distribution on the side of the median, effectively quantifying the certainty about the direction of an effect. It ranges between 0.5 (maximum uncertainty) and 1 (maximum certainty). A pd larger than 0.97 is considered to be a strong effect (Makowski et al., 2019).
We also used the full region of practical equivalence (ROPE) to define the null hypothesis over a range of values considered negligible or too small to be practically important. The ROPE approach allows us to determine whether an effect is not only strong but also if it is practically meaningful. This measure quantifies the percentage of the posterior distribution within the ROPE range (Extended Data Fig, 4-1C). In this study, we defined a 10% change in the minority class as the range of practical equivalence. A ROPE >0.97 indicates strong support for the effect being practically negligible, while a ROPE <0.03 indicates a meaningful effect size (Makowski et al., 2019).
Development and validation of automated classification
Feature extraction
All features were extracted utilizing a custom MATLAB code. The offset signal was flipped on the time axis to align the nonictal to ictal segments and trends with the onset data. We applied the IClabel (Pion-Tonachini et al., 2019a) and MARA (Winkler et al., 2011) classifiers on the ICs, both globally for each IC and locally around the onset and offset bifurcations. Extended Data Figure 2-2A descriptively illustrates that components labeled with a clear BifM had a higher mean brain probability estimate from these classifiers, indicating the potential utility of this information in detecting ICs containing seizure-related data.
For precise detection of the transition time, we analyzed a 10 s window around the seizure annotations (both onset and offset) and the corresponding time periods in the control datasets. Our goal was to identify an objective shift in the statistical characteristics of the signal, referred to as a change-point (Truong et al., 2020). We used a variance-based change–point detection algorithm (MATLAB's “ischange” function) and restricted the detection to a single change-point. When a single clear change-point was identified, it served as an anchor for subsequent feature extraction, significantly improving the quality of the features, especially for estimating amplitude and ISI changes. If no change-point was found, the dataset annotation time was used as the reference start point.
A change-point can sometimes occur outside a seizure or be related to nonbrain changes within a seizure; however, Extended Data Figure 2-2B shows a higher proportion of detected change-points in seizure data, particularly in ICs selected for having clear BifMs. A limitation of this approach is that if the change occurs close to the edge of the analysis window (less than a second), it may not be detected. Extending the analysis time to accommodate this could introduce nonstationarity biases, so we opted not to extend the analysis window.
Amplitude modulation was assessed using the peak amplitude and the signal's envelope root mean square during the first 2 s of the transition time. ISIs were determined as intervals between subsequently identified peaks within 5 s of the estimated transition time during the seizure. Established morphological curves (Saggio et al., 2020) were fitted to these values, with the fit quality assessed via root mean square error (RMSE).
Over a 5 s window from the transition time within the seizure, we computed the signal's mean, median, standard deviation, and skewness. For ISIs, both the mean interval and standard deviation were calculated. The power spectrum of the signal, ranging from 1 to 40 Hz, was computed. A linear fit was applied to the logarithm of the power spectrum against frequency, from which intercept, slope, confidence intervals, RMSE, and adjusted R squared were derived. The most deviant positive point from the linear fit line was set as the frequency peak, potentially correlating with the dominant seizure oscillation frequency.
Within the same time window, the spectral and approximate entropy—recognized as robust seizure detection metrics (Boonyakitanont et al., 2020)—were employed to pinpoint seizure-associated ICs. The signal's autocorrelation function was used to delineate cyclical patterns. The peaks in the autocorrelation function were identified, and the intervals between them were analyzed to compute their mean and standard deviation. The entropy of the autocorrelation function was calculated to provide insights into signal oscillation regularity. Additionally, the width of the autocorrelation function, a metric associated with critical slowing down (Maturana et al., 2020), was estimated as it could be related to the proximity to a bifurcation event at larger time scales.
An initial evaluation of the extracted features was conducted to address various considerations. In several instances, feature extraction failed due to brief seizure durations or inadequately identified peaks, often attributed to noise interference that resulted in missing values. For each classification task, we excluded features with over 5% missing values.
Automated classification
We separated our classification task into two stages. Initially, we focused on finding ICs that included a clear transition into or out of a seizure. Subsequently, we aimed to classify the BifM features within these ICs. The two key morphological features are the amplitude and ISI changes. In surface EEG, the amplitude may not solely reflect seizure-initiating neuronal properties but could also indicate the extent of the recruited neuronal population as the seizure propagates (Fig. 2D). Moreover, focusing on a binary classification task facilitates a more straightforward model evaluation and interpretation. Therefore, we first prioritized automating the classification of ISI change.
Developing a proper classifier for these tasks involved three main challenges. The first was the limited sample size, which heightened the risk of overfitting and reduced the potential for generalization. To mitigate this, we combined onset and offset samples to increase the sample size and integrated a categorical feature to account for inherent differences between the two groups. We used a random forest (RF) model, limiting the maximum tree depth to three to reduce overfitting (Fawagreh et al., 2014). In the context of the RF model, the categorical feature allows the model to learn distinct decision flows for onset and offset transitions, preserving the distinction between these two groups in the classification task. At the same time, this approach enables the classifier to leverage common patterns shared between onset and offset cases. The second challenge was the noisy nature of the seizure data, which included potential outliers. To address this, we employed the RobustScaler from scikit-learn (Pedregosa et al., 2011), providing a scaling approach that is robust to outliers. In regard to the final challenge, the class imbalance in our data, which could bias the classifier toward the majority class, we recalibrated the class weights to ensure equal sensitivity across both classes (Sun et al., 2009).
We evaluated performance using a leave-one-subject-out cross–validation (LOOCV) paradigm, which allowed us to assess performance across patients and estimate performance confidence intervals (Tougui et al., 2021). We measured performance using balanced accuracy (BA; Brodersen et al., 2010) across all patients. The area under the receiver operating characteristic curve (AUC-ROC) was computed for patients with both classes present. However, this computation was not feasible for patients with only one class, so it was not reported in these cases. The results for the test patients were aggregated, and overall performance was evaluated based on all test predictions. Finally, we performed a SHAP analysis to identify the main driving features for each classifier (Rodríguez-Pérez and Bajorath, 2020).
Code accessibility
The code described in the paper is freely available online at https://github.com/miriamguen/miriamguen-surface_eeg_seizure_analysis.git.
Extended Data 1
This item provides the latest version of the code zipped from our GitHub repository. Within this zipped folder is a README.md file that provides the instructions to run each analysis step presented in this manuscript. Download Extended Data 1, ZIP file.
Results
In this study, we explored the use of a blind source separation approach on surface EEG data to extract seizure information and to identify the BifM of transitions into and out of seizures according to the TSD. We analyzed 1,177 seizures from 158 patients with focal epilepsy who were undergoing presurgical evaluation. Utilizing multichannel ictal EEG data, we applied a PCA–ICA preprocessing approach to separate sources within the ictal signals, resulting in a set of eight ICs for each seizure. The ICs were visually analyzed and labeled for clear transitions and their corresponding BifMs by two independent reviewers. We then used these labels to examine the interaction of the BifMs with various clinical factors. Additionally, to extend the applicability of our findings to a broader population, we used these labels as a gold standard to develop a scalable automated classification pipeline.
Label validation
Two reviewers evaluated 1,177 seizures and a total of 9,416 ICs. The top part of Table 1 presents the number of ICs labeled with a clear transition and BifM by each reviewer and the number of ICs where both reviewers agreed on identifying a clear transition and the same BifM. Extended Data Table 1-2 contains the raw labels of both reviewers, which were used to create the summary table and perform the agreement analysis.
The overall agreement evaluation in both tasks is presented in Table 2, as evaluated using Cohen's kappa coefficient, which indicates strong agreement (Landis and Koch, 1977; McHugh, 2012). We conducted a permutation analysis to account for class imbalance in our data. The results demonstrate that inter-rater agreement in both tasks was significantly better than chance (p = 0.0001) using two different permutation approaches (Table 2).
Component agreement statistics: the table below shows the inter-rater agreement statistics, including empirical and permutation analyses of kappa values
A concern in surface EEG data is the low signal-to-noise ratio; in nonictal data, ICA is typically used to identify and remove noise sources. Consequently, several automated classifiers have been developed to distinguish between brain and noise sources. In this work, we used the mean value of two approaches: MARA (Winkler et al., 2011) and IClabel (Pion-Tonachini et al., 2019a), as briefly described in Materials and Methods, Feature extraction, and in Figure 2. Extended Data Figure 2-2A demonstrates that ICs labeled with a clear BifM correspond to a higher “brain” probability, indicating that the output of these classifiers can be used to identify ICs with clear seizure transition information. However, this is not specific to seizure data and can be high in brain sources that do not include clear seizure information, as shown in the control data and in the example in Extended Data Figure 1-1D.
We used a change-point detection algorithm to identify ICs with a clear objective transition. Extended Data Figure 1-1B descriptively demonstrates the proportion of segments with a detected change-point. The proportion was lowest for the control recordings, highest for the selected seizure ICs, and in between for the nonselected seizure time ICs. While this approach is not specific to ICs containing seizure transition information of a brain source, when combined with the brain score, it can be a good indicator of relevant seizure-related information.
Detectability and BifM distribution at the seizure level
Table 1 presents the number of seizures with a clear BifM identified by each reviewer and by both reviewers, the mean number of ICs with a clear BifM per seizure, and the number of seizures with more than one BifM. Only ICs where both reviewers agreed on the BifM are defined as “clear” in the following section. Extended Data Figure 2-2C illustrates the selection flow for determining the seizure-level BifM used in cases where more than one IC included a clear BifM per seizure.
In 49.5 and 40.3% of seizures, we identified onset and offset BifMs that met our criteria. Onset and offset BifMs were not identified in 18 and 22 patients, respectively. The mean number of unique BifMs was higher in patients with more seizures (Extended Data Fig. 2-2D), suggesting that the number of observed seizures limits the number of observed BifMs. Detectability was better for seizures with a posterior SOZ and poorer for those with a frontal SOZ. Seizures with stronger clinical manifestations had a higher proportion of clear bifurcations, such as in FIAS and FBTCS, than in FAS or UC seizures (Extended Data Fig. 2-2E,F). Figure 3A demonstrates the overall distribution of onset and offset BifMs, showing that their rank aligns with previous findings (Saggio et al., 2020). Table 3 summarizes the onset and offset BifM pairs identified.
BifM distribution. The distribution of BifM in our data, encompassing all seizures marked by the onset or offset bifurcations, aligns broadly with previously reported intracranial data. However, the proportions of SN/SubH and FLC morphologies were higher in surface data. These morphologies do not necessitate specific patterns, such as ISI or amplitude changes. Additionally, the higher noise levels in surface data may obscure other morphologies, potentially leading to an overestimation of these patterns.
Table 3-1
The empirical dynamotype prevalence in the data is ordered by the mathematical complexity ranking, showing that low complexity is aligned with higher prevalence. Download Table 3-1, XLSX file.
A summary of the onset and offset seizure BifM, including seizures with at least one (onset or offset) clear BifM
The prevalence of dynamotypes is hypothesized to be inversely related to their complexity. Previous work estimated this complexity using the canonical model, assuming hysteresis-loop bursting, i.e., bistability and feedback-based slow variables driving the onset and offset bifurcations Extended Data Table 3-1 summarizes the dynamotype count in our data, ordered by this ranking. Consistent with theoretical expectations, the most prevalent dynamotypes were of lower complexity. However, some exceptions were observed, which may be attributed to several factors limiting the application of the burster ranking to the data. These include the mathematical ranking derived for hysteresis-loop bursting dynamics, while other dynamics, such as slow-wave bursting, may also exist. Additionally, the FLC at the offset and SN/SubH at the onset lack specific features for identification, combined with the noise in surface EEG, which may lead to their overestimation in the visual analysis.
Relationship with clinical factors
The comparisons that met the defined criteria are presented in Table 4, with the complete analysis results provided in Extended Data Tables 4-1 and 4-2. Robust findings were observed only in the comparisons of onset BifMs. While other trends were noted in the data, a larger cohort is necessary to draw definitive statistical conclusions.
Onset comparisons in the Bayesian posterior analysis have a high pd and fall mostly outside the ROPE
Based on the defined criteria, at the patient level, seizures associated with cortical dysplasia had a higher proportion of SupH than those associated with hippocampal sclerosis. In a previous study, the SupH proportion was higher with older age (Saggio et al., 2020). Our data showed a similar trend that did not meet the defined criteria (Extended Data Table 4-1).
At the seizure level, no differences in BifM proportions were observed across the NREM stages. However, during wakefulness, there was a higher proportion of SupH BifMs compared with NREM1 and NREM2 and a higher proportion of SNIC BifMs compared with NREM3 (Table 4). Additionally, FIAS and UC seizures had a higher proportion of SupH BifMs than FBTCS and FAS seizures.
Extended Data Table 4-3 shows the relationship between vigilance and clinical classification, replicating previous findings in focal onset seizures (Herman et al., 2001). FBTCS, which occurs more frequently in NREM, had a lower proportion of SupH BifMs than FIAS and UC, possibly indicating a shared effect. In contrast, FAS seizures mostly begin during wakefulness, in which SupH BifMs are more common. FIAS seizures occur at similar rates during wakefulness, while UC seizures are less frequent. However, FAS seizures have a lower proportion of SupH BifMs compared with both FIAS and UC seizures. These observations suggest distinct connections between dynamical patterns, vigilance states, and the clinical manifestations of seizures.
Automated classification
Automated classification of ICs for clear seizure transition
We analyzed 9,416 ICs for seizure onset and offset to identify clear transitions, resulting in 18,832 examples; 1,103 had a clear onset, and 884 had a clear offset transition. In cases where feature extraction was unsuccessful (e.g., due to noise or insufficient peaks), we used empty values to ensure output continuity. We evaluated performance using the LOOCV paradigm for each classification task, estimating performance with BA and AUC-ROC. Figure 4A presents the normalized confusion matrix from the cumulative labels of the cross-validation process. Table 5 summarizes the performance results for all subjects, while Extended Data Table 5-1 details the results by patient.
General classifier analytics. This figure presents the analysis of two classification tasks. The first classifier identifies components with a clear transition among all given components, while the second classifier identifies components exhibiting ISI slowing based on BifM. Both classifiers were trained using the same feature set and algorithm hyperparameters, with class weights adjusted differently to account for class imbalance and to achieve balanced performance. Each classifier was evaluated using a LOOCV approach. The confusion matrices (panels A and C) display cumulative classification results on the test sets. To understand the classification process, feature contributions were analyzed using SHAP analysis (Rodríguez-Pérez and Bajorath, 2020), highlighting the most influential features of each classifier and their directional effects. A, The confusion matrix for the transition detection task. B, Dominant features for the transition detection task: high brain probability scores were associated with clear components, while high probability scores for undefined noise sources (“other”) indicate a lower likelihood of being seizure components. High spectral entropy, a feature previously used in seizure detection tasks (Boonyakitanont et al., 2020), supported positive decisions. High mean peak amplitude in the first 5 s of detected onset/offset was related to clear seizure transitions. C, The confusion matrix for the ISI slowing classification task. D, Top features for the ISI slowing classification task: a value indicating an offset increased the probability of ISI slowing due to the higher likelihood of this pattern occurring at offsets. Seizure length was positively correlated with ISI probability, although no significant trend was observed in multilevel analysis. Lower dominant frequency and higher standard deviation in ISIs were associated with ISI slowing. The ISI linear intercept estimate (higher CI limit) and lower fit quality (RMSE) to the ISI trend curves were positively related to slowing. This may be due to low error in cases where the fitted line was flat and the observation that slowing usually occurred in components with higher amplitudes, as suggested by the peak trend estimation relation.
Figure 4-1
The Bayesian measures used for evaluation. This figure provides an intuitive illustration of the measures presented in Table 4. We evaluated the difference between two posterior distributions of bifurcation morphologies, given the tested clinical factors compared to the null hypothesis. The reported measures used in this work are as follows, as suggested by (Makowski et al., 2019). Download Figure 4-1, TIF file.
Table 4-1
Onset bifurcation morphology distribution analysis. This table provides all the conditions compared in the current analysis. Download Table 4-1, XLSX file.
Table 4-2
Offset bifurcation morphology distribution analysis. This table provides all the conditions compared in the current analysis. Download Table 4-2, XLSX file.
Table 4-3
Seizure clinical classification crossed with sleep stages. This table provides a contingency matrix for 572 seizures analyzed for onset bifurcations (Saggio et al., 2017). In this table, we can observe that seizures undergoing generalization have a higher prevalence in deeper stages of sleep, as previously reported. Interestingly, we also see a high portion of generalization during REM. Download Table 4-3, XLSX file.
Seizure transition classification performance metrics, as measured by LOOCV
Table 5-1
The full by-patient results of the seizure transition classification results and performance metrics, as measured by leave-one-subject-out cross-validation. Download Table 5-1, XLSX file.
Figure 4B presents the SHAP analysis results, highlighting the top 10 features influencing the model's decisions and the direction of the effect for each example. Eight of these features were derived from the MARA and IClabel classifiers. Low scores, indicating that the ICs are not from a brain source, led the classifier to an “unclear” label. In contrast, high scores showed a weak association with the “clear” label, suggesting that this information is necessary but insufficient for identifying clear brain-origin seizure transitions. Spectral entropy, a well-known feature used for seizure detection (Acharya et al., 2015), strongly identified an IC as seizure-related. A higher peak amplitude also contributed to classifying the data as containing seizure information.
Automated classification of ICs for ISI increase
For this stage, we used the ICs for which both reviewers agreed on the BifM during visual labeling. This included 1,103 onset ICs and 884 offset ICs (Table 1). Of these, 52 onset and 234 offset ICs showed an increase or decrease in ISI (from and to zero). The summary of the LOOCV results is presented in Table 6, with the detailed performance per patient provided in Extended Data Table 6-1. The overall confusion matrix reflecting performance is shown in Figure 4C, and the SHAP analysis results for this task can be seen in Figure 4D. The top feature was the onset/offset indicator, as ISI slowing was more prevalent at the offset in our data. Features quantifying ISI changes had the most significant impact on the classifier's decision. The width of the autocorrelation function emerged as a strong predictor of ISI slowing (Maturana et al., 2020).
Seizure ISI slowing classification performance metrics, as measured by LOOCV
Table 6-1
Full by-patient results of the ISI slowing classification results and performance metrics, as measured by leave-one-subject-out cross-validation. Download Table 6-1, XLSX file.
Discussion
Recent studies employing a phenomenological modeling approach to develop the TSD have demonstrated its potential to provide clinically relevant insights. These computational models have been used to construct digital simulations representing the epileptogenic network in patients undergoing presurgical evaluation, improving surgical outcomes (An et al., 2019; Makhalova et al., 2022). Additionally, they have been utilized to optimize electrical stimulation for seizure abortion (Szuromi et al., 2023). Empirical findings from intracranial recordings (electrocorticography, ECoG) in humans have demonstrated that the prevalence of onset and offset bifurcations aligns with theoretical predictions (Saggio et al., 2020). Our study extended this classification to noninvasive EEG recordings, addressing several challenges.
The first challenge is the mixture of physiological and nonphysiological sources in surface EEG signals, requiring effective methods to isolate the information of interest. We addressed these challenges by employing a blind source separation approach, which decomposes the complex signals into ICs that are as statistically independent as possible (Makeig et al., 1995). We visually identified the ICs containing seizure transition information relevant to the TSD classification while accounting for the limitation of missing DC component in our recordings. This component distinguishes several onset and offset classes and has potential clinical relevance in predicting the response to seizure-aborting stimulation (Szuromi et al., 2023). This remains a limitation of our current analysis approach, which uses a high-pass filter that would eliminate a DC component if contained in the data. However, if DC information were available, this limitation could be addressed by applying a low-pass filter to the original signal to extract the DC components. Using the learned weights, the ICA transformation could then be applied to the low-frequency components, and this additional information could be incorporated into the classification.
The second challenge hinges on the essential difference between the nature of the recorded signals in the current and previous work. When recorded on the surface, the brain signal undergoes volume conduction and smearing as it traverses various mediums, such as the cerebrospinal fluid, dura mater, skull, and skin, before reaching the sensors (van den Broek et al., 1998). This process can alter the signal generated by neuronal tissue, affecting the morphology of the data. Additionally, the neuronal population required to detect activity in surface EEG is larger compared with intracranial recordings (Tao et al., 2007a,b), potentially reflecting different dynamical behaviors, such as seizure spread and propagation, that are not evident from a single electrode in invasive recordings. Conversely, the seizure onset may appear later compared with invasive recordings or may not be apparent in surface recording, depending on the depth of the ictal source (Koessler et al., 2010; Barborica et al., 2021). Thus, the dynamics captured in the surface data may reflect a larger-scale process of the epileptogenic network, and its interaction with the global brain state must be considered in the interpretation of our results.
Due to the inherent noise in the data, we applied strict labeling criteria, one of which required agreement between two independent reviewers to classify a transition as clear. This approach minimized the risk of misidentifying noise as seizure dynamics; for more detail, see Materials and Methods. Using these criteria, we identified clear transitions in 49.5% of the onsets and 40.3% of the offsets (Table 1). Clear onset and offset BifMs could not be identified in 18 and 22 patients, respectively. Detection rates were highest in occipital localizations and decreased toward frontal areas (Extended Data Fig. 2-2E). Additionally, detectability was better in more clinically pronounced seizures (Extended Data Fig. 2-2F), highlighting the influence of the involved neuronal tissue on the prominence of seizure activity in surface EEG (Tao et al., 2007b).
Our approach was based on the assumption that seizure information consists of high-amplitude synchronized activity. Consequently, seizure onsets characterized by low-voltage fast activity were potentially missed. However, we often observed onset transitions in ICs that began with low-amplitude fast activity and evolved into high-amplitude low–frequency activity. Previous work demonstrated a strong correlation between these interictal and surface IC patterns, with detectability on the surface primarily influenced by the depth of the seizure source (Barborica et al., 2021).
Overall, the BifM ranking aligned with previous findings (Saggio et al., 2020). However, the proportion of offset FLC was notably higher in our data, at 71.5% (Fig. 3) compared with the previously reported 53.6%, and onset SN/SubH was also higher, at 88% compared with 71.1%. Several factors may explain this discrepancy.
First, prior studies used intracranial recordings, which directly capture data from smaller neuronal populations at the seizure focus, reflecting finer spatial and temporal dynamics. In contrast, scalp EEG requires synchronized activity over larger brain regions to produce a detectable signal. Consequently, scalp EEG recordings reflect seizure activity at a different temporal and spatial scale, as detailed above. This inherent difference is likely a significant factor contributing to the discrepancies observed between our data and previous studies.
Second, the noise profile in surface EEG data differs, and despite preprocessing efforts, this noise could obscure features like amplitude or spike rate modulation. This may lead to the overclassification of BifMs that do not rely on distinct morphological features and underestimation of others. Third, from a clinical perspective, the proportion of BifMs (especially SupH) may be influenced by factors such as sleep stage, seizure type, or etiology (Table 4), making direct comparisons between studies challenging without accounting for these variables.
Previous work reported low numbers of SupH offset transitions, with 3% in recordings without DC components and 6% in recordings with them, consistent with TSD expectations (Saggio et al., 2020). In our data, this rate was even lower, at 0.4% (n = 3), which is possibly attributable to the distance between the seizure source and the recording electrodes. The TSD framework models seizure dynamics from a point source; however, the signal measured in surface EEG is indirect and undergoes transformations that may hinder the accurate detection of these patterns. This is further illustrated by a trend observed in seizure pathology. While SupH offsets were excluded from statistical modeling due to the small sample size, the onset analysis revealed a higher estimate of SupH BifMs in cortical dysplasia (a superficial source) compared with hippocampal sclerosis (a deeper source; Table 4)—a trend not found in intracranial data. Depending on the etiology, such considerations may also be relevant in intracranial data, particularly when comparing ECoG and depth electrodes.
At the onset transition, the observed amplitude increase in the SupH BifM in surface EEG may reflect propagation or be influenced by noise levels in the signal, rather than exclusively representing the dynamics of the neuronal population in the SOZ. In some cases, we observed amplitude increases related to propagation (Fig. 2D). While we attempted to differentiate these instances, some propagation events may have been misclassified as SupH. Typically, higher propagation is expected in seizures clinically classified as FBTCS. However, we observed lower SupH BifM rates in FBTCS compared with FIAS and UC seizures (Table 4), with rates nearly equivalent between FBTCS and FAS (Extended Data Table 4-1). These findings suggest that this morphology in surface EEG reflects dynamics beyond mere seizure propagation.
Several approaches could enhance our understanding of the surface and source signal differences. Simultaneous recordings have previously been employed to validate the correlation between source and surface activity, offering valuable insights into the spatial and temporal dynamics of seizures (Tao et al., 2007b; Koessler et al., 2010; Barborica et al., 2021). These approaches could be particularly useful for clarifying how SupH bifurcations manifest on the surface. Source imaging has shown significant utility in localizing epileptogenic foci and improving presurgical evaluations by integrating high-density surface EEG data with structural information to infer the origins of electrical activity (Michel et al., 2004; Plummer et al., 2008). This noninvasive technique could enhance the reconstruction of the SOZ time series, potentially providing a more accurate representation of the underlying neuronal dynamics.
The modeling framework underlying the development of the TSD is phenomenological, meaning it uses a mathematical model to describe patterns of neuronal activity based on the observed behavior rather than explicitly linking them to biological mechanisms. This approach simulates the transitions into and out of seizures, with a dynamical fixed point (FP) representing the resting state and a limit cycle (LC) representing the oscillatory activity during a seizure. These states capture the general patterns of brain activity during seizures, modeling their initiation and cessation without directly specifying the underlying cellular or network-level processes (Saggio et al., 2020).
The theoretical framework describes the seizure driving forces at three main time scales: fast, slow, and ultraslow, which drive the transitions into and out of a seizure. These transitions can occur through two mechanisms. First, when multiple attractors coexist, such as a stable FP and LC (bistability), the dynamics can be altered by introducing noise that shifts the system from one attractor to another. Noise has been previously demonstrated to drive transitions into seizures in neuronal networks, particularly in bistable systems, where random fluctuations can push the system from a stable resting state to an oscillatory seizure state (Deco et al., 2009). Second, changing the system's parameters can result in a sudden alteration of its dynamical map—for instance, by eliminating an FP and creating an LC—thereby causing the transition from a resting state to an oscillatory state through a bifurcation. Bifurcations have been widely proposed as a core mechanism underlying seizure onset, where small parameter changes in neuronal systems can lead to large-scale shifts in network dynamics (Da Lopes Silva et al., 2003). This bifurcation process is central to the theoretical foundation of the TSD framework (Saggio et al., 2020).
The TSD mathematical framework explicitly models the fast and slow variables. The fast variable is responsible for modeling the oscillatory activity of seizures and can transition from FP to LC dynamics. The slow variable is modeled with hysteresis-loop bursting, essentially using input from the fast subsystem to provide feedback to the fast variable, leading to the cessation of oscillations and providing an intrinsic mechanism to limit seizure duration (Izhikevich, 2000; Jirsa et al., 2014). A one-dimensional slow system can model this feedback loop, offering a simplified mathematical structure while retaining essential dynamical features. Thus, this framework forms the basis for the dynamotype ranking used in our analysis (Saggio et al., 2017).
This modeling approach was initially developed to simulate local seizure dynamics, focusing on the transitions into and out of seizures, where the slow variable's dynamics were primarily attributed to the tissue in the SOZ or its immediate surroundings (Jirsa et al., 2014). However, in seizures occurring within the broader context of the brain, interactions between the SOZ and other brain regions must be considered. These external regions can modulate the epileptogenic network, particularly in different sleep states (Moore et al., 2021). Shifting from local dynamics to a more comprehensive view that includes these interactions can deepen our understanding of seizure generation and propagation.
Models like the “virtual epileptic brain” (Jirsa et al., 2017) have been developed to simulate large-scale brain dynamics, focusing on the onset, propagation, and termination of seizures. These models use different excitability assumptions for the slow dynamics of the epileptogenic network, the propagation network, and the noninvolved zones. However, they primarily operate on short timescales and do not consider the ultraslow modulation of excitability. To better capture seizure dynamics over longer timescales, it could be beneficial to introduce an explicit ultraslow parameter that modulates the excitability in the different zones reflecting the ultraslow physiological processes. For example, this modulation could leverage established frameworks like the two-process model of sleep–wake regulation (Borbély, 2022) or utilize the empirically predicted seizure susceptibility cycles found in many patients with epilepsy (Karoly et al., 2021). This integration has the potential to shed light on the interaction between processes at multiple timescales driving seizure dynamics.
In this work, we began by applying the taxonomy to larger-scale brain signals, acknowledging the limitations discussed earlier. This allowed us to test consistency and uncover potential insights for future modeling at this scale. First, we analyzed the distribution of BifMs in the dataset and compared these with prior findings in local recordings. While the overall rank remained similar, the absolute proportions, particularly for SupH BifMs, differed, potentially due to noise, as noted earlier. In our analysis of the relationship between clinical factors, most findings were linked to changes in SupH proportions. However, given the limitations, these findings may be influenced by technical aspects rather than seizure dynamics. Conversely, changes in the ISIs are less likely to be affected by propagation artifacts, making them a more reliable target for investigation using this modality.
Our analysis reveals that the proportion of SNIC onsets was reduced during NREM3 sleep compared with the awake state (Table 4). Similar but weaker trends were noted when comparing wakefulness to NREM1 and NREM2 (Extended Data Table 4-1). In these comparisons, the SupH proportions were consistently higher during wakefulness, while the SN/SubH proportion increased during NREM sleep (Extended Data Table 4-1). Overall, SNIC and SupH onsets were more prevalent during wakefulness, aligning more closely with intracranial findings and theoretical expectations. NREM3, or slow-wave sleep (SWS), differs significantly from wakefulness in behavior and electrophysiology. In the context of seizures, research has shown that seizures are more likely to generalize during SWS (Herman et al., 2001), a finding reflected in our data (Extended Data Table 4-3).
Moreover, synchronized NREM sleep generally facilitates epileptogenic activity (Halász et al., 2019), while desynchronized REM sleep tends to be seizure-protective (Ng and Pavlova, 2013) and shows more localizable epileptogenic activity (McLeod et al., 2020). These findings demonstrate the importance of considering sleep stages in understanding seizure dynamics. In the context of the modeling approach, these observations could be viewed as the modulation of the ultraslow variable, leading to a change in excitability affecting the seizure rate and propagation or to the movement on the state space map affecting the observed BifMs.
Although our dataset includes many seizures (1,177), not all could be labeled for their BifM due to noise contamination (Table 1). Expanding this analysis to additional noninvasive datasets would increase the number of participants and the range of clinical conditions considered (e.g., generalized epilepsies). However, the visual classification process is labor-intensive and prone to inter-rater variability.
Previous studies have demonstrated significant inter-rater variability in EEG interpretation, highlighting the importance of multiple reviewers to improve diagnostic reliability (Grant et al., 2014; Zibrandtsen et al., 2019; Aanestad et al., 2024). Unlike previous studies focusing solely on time series analysis and long-term EEG recordings for rare seizures, we utilized data with prelabeled seizure onset and offset times. Within each seizure, we identified ICs containing seizure information from brain sources, relying on time series, IC topography, spectral profile, and automated source rating. We also conducted a training protocol with a large dataset of ICs and simulated data to enhance reviewer consistency. This protocol resulted in an inter-rater agreement that was significantly better than chance and is considered moderate to high (McHugh, 2012). Nonetheless, it is important to acknowledge that involving only two reviewers may limit the robustness of our findings. Future studies should consider including additional reviewers to enhance the reliability of EEG classifications.
Automating this process would enable the extension of this analysis to larger cohorts and a broader population. In this work, we used a baseline classification algorithm to examine the utility of features for identifying seizure ICs and classifying the ISI slowing (based on the limitations mentioned above). We demonstrated good cross-subject performance in both tasks. The top contributing features for identifying seizure information included brain score indicators (Winkler et al., 2011; Pion-Tonachini et al., 2019a) and entropy measures typically used for seizure detection (Boonyakitanont et al., 2020). For classifying ISI slowing, the top features were morphological characteristics and the width of autocorrelation function (Maturana et al., 2020). This classification framework could be further tested and improved by employing a more extensive dataset.
Current and previous work has focused primarily on the seizure onset and offset bifurcations. Automated detection tools, particularly change-point analysis (Truong et al., 2020) throughout a seizure, combined with automated classification, may enable the dynamics throughout the spatiotemporal evolution of seizures in empirical data to be quantified, providing additional input for modeling epileptogenic networks in patients. Moreover, by reducing the labeling workload, such automated classifiers could extend this analysis to a more extensive database, allowing seizure dynamic changes across a broader range of clinical conditions to be examined with greater statistical power, feeding back information and improving modeling approaches.
Footnotes
The authors declare no competing financial interests.
We acknowledge Matan Ben-Shahar for his statistical consultation and Carmel Salomonski's assistance in building the graphical user interface and conducting the initial analysis. Their contributions greatly enhanced the quality of our study. We are grateful for their expertise and support. This research was supported by Grants 504/17 and 794/22 from the Israel Science Foundation (https://www.isf.org.il).
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.
References
In this issue
