# Neuronal avalanches in temporal lobe epilepsy as a noninvasive diagnostic tool investigating large scale brain dynamics

Marie-Constance Corsi, Emahnuel Troisi Lopez, Pierpaolo Sorrentino, Alberto

Danieli, Simone Cuozzo, Paolo Bonanni, Gian Marco Duma

► To cite this version:

Marie-Constance Corsi, Emahnuel Troisi Lopez, Pierpaolo Sorrentino, Alberto Danieli, Simone Cuozzo, et al. Neuronal avalanches in temporal lobe epilepsy as a noninvasive diagnostic tool investigating large scale brain dynamics. Scientific Reports, 2024, 14 (1), pp.14039. 10.1038/s41598-024-64870-3. hal-04618262

HAL Id: hal-04618262

https://inria.hal.science/hal-04618262

Submitted on 20 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Neuronal avalanches in temporal lobe epilepsy as a noninvasive diagnostic tool investigating large scale brain dynamics

Marie-Constance Corsi, Emahnuel Troisi Lopez, Pierpaolo Sorrentino, Simone Cuozzo, Alberto Danieli, Paolo Bonanni, Gian Marco Duma

<figure><img src="image_1.png" /></figure>

## To cite this version:

Marie-Constance Corsi, Emahnuel Troisi Lopez, Pierpaolo Sorrentino, Simone Cuozzo, Alberto Danieli, et al.. Neuronal avalanches in temporal lobe epilepsy as a noninvasive diagnostic tool investigating large scale brain dynamics. Scientific Reports, 2024, 14 (1), pp.14039. 10.1038/s41598-024-64870-3. hal-04618262

HAL Id: hal-04618262

https://inria.hal.science/hal-04618262v1

Submitted on 1 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Neuronal avalanches in temporal lobe epilepsy as a diagnostic tool: a noninvasive
investigation of intrinsic resting state dynamics

Marie-Constance Corsi 1*†, Emahnuel Troisi Lopez 2†, Pierpaolo Sorrentino 3,4*, Alberto
Danieli 5, Simone Cuozzo 5, Palo Bonanni 5, Gian Marco Duma 5‡

Authors email

Marie-Constance Corsi: marie-constance.corsi@inria.fr;

Pierpaolo Sorrentino: pierpaolo.SORRENTINO@univ-amu.fr

Gian Marco Duma: gianmarco.duma@lanostrafamiglia.it

1 Sorbonne Université, Institut du Cerveau – Paris Brain Institute -ICM, CNRS, Inria, Inserm, AP-HP, Hopital de la Pitié Salpêtrière, F-75013, Paris, France

2 Institute of Applied Sciences and Intelligent Systems of National Research Council, Pozzuoli,
Italy

3 Institut de Neurosciences des Systèmes, Aix-Marseille Université, 13005 Marseille, France

4 University of Sassari, Department of Biomedical Sciences, Viale San Pietro, 07100, Sassari, Italy

5 IRCCS E. Medea Scientific Institute, Epilepsy Unit, Via Costa Alta 37, 31015, Conegliano, Treviso, Italy

† Co-First authors

* Co-Corresponding authors

‡ Last author
<figure><img src="image_3.png" /></figure>

OPEN

# Neuronal avalanches in temporal lobe epilepsy as a noninvasive diagnostic tool investigating large scale brain dynamics

Marie-Constance Corsi¹,⁶✉, Emahnuel Troisi Lopez²,⁶, Pierpaolo Sorrentino³,⁴✉, Simone Cuozzo⁵, Alberto Danieli⁵, Paolo Bonanni⁵ & Gian Marco Duma⁵

The epilepsy diagnosis still represents a complex process, with misdiagnosis reaching 40%. We aimed at building an automatable workflow, helping the clinicians in the diagnosis of temporal lobe epilepsy (TLE). We hypothesized that neuronal avalanches (NA) represent a feature better encapsulating the rich brain dynamics compared to classically used functional connectivity measures (Imaginary Coherence; ImCoh). We analyzed large-scale activation bursts (NA) from source estimation of resting-state electroencephalography. Using a support vector machine, we reached a classification accuracy of TLE versus controls of 0.86 ± 0.08 (SD) and an area under the curve of 0.93 ± 0.07. The use of NA features increase by around 16% the accuracy of diagnosis prediction compared to ImCoh. Classification accuracy increased with larger signal duration, reaching a plateau at 5 min of recording. To summarize, NA represents an interpretable feature for an automated epilepsy identification, being related with intrinsic neuronal timescales of pathology-relevant regions.

Epilepsy is one of the most diffused neurological conditions, affecting more than 50 million people worldwide, representing a relevant burden to health services¹. Even if consensus statements and international guidelines are present, the diagnosis of epilepsy still represents a complex scenario, with a rate of misdiagnosis, which varies widely across different centers, reaching 40% in specific cases²,³. In this context, medical doctors may need support in the diagnostic process⁴. Recent years have witnessed the introduction of artificial intelligence (AI) methodologies as a supporting tool for clinicians both in the diagnostic and treatment processes. In the field of epilepsy, vast efforts were conducted in order to translate the large amount of multimodal data that are nowadays collected during the clinical practice, towards the improvement, and potentially the automation, of the diagnostic process⁵,⁶. Across the multiple neuroimaging modalities, electroencephalography (EEG) represents the election tool in epilepsy diagnosis⁷. Previous studies have identified alteration of EEG power spectrum as a potential biomarker helping the diagnosis⁸⁻¹⁰. Importantly, the majority of applications of AI on the EEG signal have focused on the automatic recognition of epileptiform/seizure activity or seizure forecasting¹¹⁻¹³. However, the simple but fundamental diagnostic question, namely the identification of epileptic patients, is far from being efficiently addressed in an automated fashion. Very few studies exploited AI trying to provide an automatic tool aiding the clinician in the epilepsy diagnosis¹⁴⁻¹⁸. In fact, these investigations displayed a large variability in the accuracy and sensitivity which can be related to the EEG features used to train the model, opening the question of which is the feature of the EEG signal better suited for training an AI model for diagnosis. Until now, most of the studies attempting to perform AI-based diagnosis exploited functional connectivity in the sensor space, derived by the covariation between signals and/or spectral power from electrodes activity, as well as the frequency-specific phase locking (i.e. synchronization)¹⁹. However, such techniques rely on the assumption of stationarity of the EEG signals, which may represent an important limitation, and might affect the stability of the estimates. In the present work, we aimed at identifying a better suited EEG signal feature to train an AI-model-based classification

¹Sorbonne Université, Institut du Cerveau – Paris Brain Institute -ICM, CNRS, Inria, Inserm, AP-HP, Hopital de la Pitié Salpêtrière, 75013 Paris, France. ²Institute of Applied Sciences and Intelligent Systems of National Research Council, Pozzuoli, Italy. ³Institut de Neurosciences des Systèmes, Aix-Marseille Université, 13005 Marseille, France. ⁴Department of Biomedical Sciences, University of Sassari, Viale San Pietro, 07100 Sassari, Italy. ⁵Epilepsy Unit, IRCCS E. Medea Scientific Institute, Via Costa Alta 37, 31015 Conegliano, Treviso, Italy. ⁶These authors contributed equally: Marie-Constance Corsi and Emahnuel Troisi Lopez. ✉email: marie-constance.corsi@inria.fr; pierpaolo.SORRENTINO@univ-amu.fr
Abstract

The epilepsy diagnosis still represents a complex process, with misdiagnosis reaching 40%.
We aimed at building an automatable workflow, helping the clinicians in the diagnosis of
temporal lobe epilepsy (TLE). We hypothesized that neuronal avalanches (NA) represent a
feature better encapsulating the rich brain dynamics compared to classically used functional
connectivity measures (Imaginary Coherence; ImCoh). We analyzed large-scale activation
bursts (NA) from source estimation of resting-state electroencephalography. Using a support
vector machine, we reached a classification accuracy of TLE vs. controls of 0.86 ± 0.08 (SD)
and an area under the curve of 0.93 ± 0.07. The use of NA features increase by around 16% the
accuracy of diagnosis prediction compared to ImCoh. Classification accuracy increased with
larger signal duration, reaching a plateau at 5 minutes of recording. To summarize, NA
represents an interpretable feature for an automated epilepsy identification, being related with
intrinsic neuronal timescales of pathology-relevant regions.

Introduction

Epilepsy is one of the most diffused neurological conditions, affecting more than 50 million people worldwide, representing a relevant burden to health services ¹. Even if consensus statements and international guidelines are present, the diagnosis of epilepsy still represents a complex scenario, with a rate of misdiagnosis, which varies widely across different centers, reaching 40% in specific cases ²,³. In this context, medical doctors may need support in the diagnostic process ⁴. Recent years have witnessed the introduction of artificial intelligence (AI) methodologies as a supporting tool for clinicians both in the diagnostic and treatment processes. In the field of epilepsy, vast efforts were conducted in order to translate the large amount of multimodal data that are nowadays collected during the clinical practice, towards the improvement, and potentially the automation, of the diagnostic process ⁵,⁶. Across the multiple neuroimaging modalities, electroencephalography (EEG) represents the election tool in epilepsy diagnosis ⁷. Previous studies have identified alteration of EEG power spectrum as a potential biomarker helping the diagnosis ⁸⁻¹⁰. Importantly, the majority of applications of AI on the EEG signal have focused on the automatic recognition of epileptiform/seizure activity or seizure forecasting ¹¹⁻¹³. However, the simple but fundamental diagnostic question, namely the identification of epileptic patients, is far from being efficiently addressed in an automated fashion. Very few studies exploited AI trying to provide an automatic tool aiding the clinician in the epilepsy diagnosis ¹⁴⁻¹⁸. In fact, these investigations displayed a large variability in the
(patients with temporal lobe epilepsy vs. healthy controls). We hypothesized that the aperiodic activation bursts (neuronal avalanches) of source-derived EEG activities may represent a better measure for an AI-model, since they are not based on the assumptions of stationarity. In particular, neuronal avalanches (NA) belong to the framework of criticality, which lends itself nicely to the investigation of the microscopic dynamical behavior of neuronal assemblies and its (non-linear) relation to large-scale alterations such as, in the case of epilepsy, seizure initiation, propagation and termination<sup>20,21</sup>. Recent findings suggest that NA spreading on the large scale represents a more stable and reliable feature for the individualized investigation and task classification of brain dynamics<sup>22,23</sup>. Relevantly, our previous work has highlighted the altered spreading of NA and its relation to brain morphology in temporal lobe epilepsy (TLE), using the source-derived high-density EEG (hdEEG) activity at rest<sup>24,25</sup>. In the light of this, here we hypothesized that the spreading of the aperiodic activation bursts at the whole-brain level may represent a relevant feature towards the automatic classification of patients with TLE vs. controls. To this end, we recorded resting-state high-density EEG (hdEEG) from patients with temporal lobe epilepsy and a control group. The model of TLE is of particular interest as it represents a rather well-defined group of electroclinical conditions with lower clinical heterogeneity compared to other epilepsy forms. Using source-reconstructed EEG signals, we individuated the presence of NA and used them to calculate the probability of consecutive activations between brain regions (i.e., the topographical spreading of the avalanches) (see Supplementary Fig. 1). This information was stored within an adjacency matrix named avalanche transition matrix (ATM)<sup>26</sup> which was then used as a feature for a support vector machine (SVM) classifier. We hypothesized that ATMs would have performed better, compared to the power-based connectivity metric generally used in literature, differentiating with a larger accuracy the two conditions (epilepsy vs. non-epilepsy). Relevantly, we computed the ATM from a resting state activity free from epileptiform activity, in order to investigate if the basal functional organization of the system, without the presence of characteristic graphic-elements (i.e., seizure and/or interictal epileptiform discharges), may contribute to a reliable classification of the patients. We then looked into the most relevant feature pattern, namely the most informative elements used by the model to differentiate the two groups. We expected that the dynamics of those areas linked with the pathophysiology of TLE, would have informed the model more than other brain regions. This endows our results with interpretability and allows the identification of clusters of regional interactions that might be relevant to the pathology. Finally, we investigated the classification performance with respect to the duration of the brain signals. Hence, we have performed the classification based on the growing amount of data used to compute ATMs, ranging from 5 to 300 s. To summarize, the purpose of the present work is to exploit the spreading of the aperiodic activities at the whole-brain level as a feature to inform an AI-based model aimed at the automatic identification of epileptic patients, in order to provide the clinicians with a practical and non-invasive tool supporting the diagnostic process.

# Results

## Model prediction performance

The implemented SVM model classified the participants as TLE or controls based on the global system dynamics, namely the ATMs, computed on the broadband signal (3-40 Hz), with a mean accuracy of (0.86 ± 0.08). We observed that the prediction accuracy of the model trained on the ATMs was larger as compared to the one trained with ImCoh (0.71 ± 0.12). These results were also confirmed for the signals filtered in the narrowband (3-14 Hz) (ATMs accuracy: 0.87 ± 0.08; ImCoh accuracy: 0.82 ± 0.11). The better classification performance for ATMs versus ImCoh was also confirmed by an increased area under the curve (AUC), both when considering the broadband and the narrowband filtered signal. The improvement of classification related to the ATMs is also confirmed when considering multiple frequency bands (i.e., alpha-beta, beta-low-gamma and alpha-beta-low-gamma; see Supplementary Table 1). In general, the results display a better classification across all the performance indices in the model trained with the ATMs as compared to the ImCoh (see Table 1).

We note that, as suggested by recent findings, providing exclusively the mean across splits of the classification accuracies might be misleading for the interpretation of the model performance<sup>27</sup>. Hence, here we display the full distribution of the classification accuracies across the multiple splits of the cross-validation. Notably, the distribution of the prediction accuracies with the ATMs is narrower, and skewed towards larger accuracy values (see Fig. 1A). These result are confirmed across all the frequency bands (see Supplementary Fig. 2). Additionally, in Fig. 1B, we also provide the mean receiver operating characteristic (ROC) curve and the corresponding AUC, providing further evidence of the increased performance of the ATMs compared to the ImCoh. In fact, in the

<table><caption>**Table 1.** Classification performance. The present table displays the indices of the classification performance when considering the model trained with avalanche transition matrix (ATM) or imaginary coherence (ImCoh) obtained from the narrowband (3-14 Hz) or broadband (3-40 Hz) signal. Higher performance values are in bold.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Broadband (3-40 Hz) mean ± SD</th><th colspan="2">Narrowband (3-14 Hz) mean ± SD</th></tr><tr><th>ImCoh + SVM</th><th>ATM + SVM</th><th>ImCoh + SVM</th><th>ATM + SVM</th></tr></thead><tbody><tr><td>Accuracy</td><td>0.71 ± 0.12</td><td><strong>0.86 ± 0.08</strong></td><td>0.82 ± 0.11</td><td><strong>0.87 ± 0.08</strong></td></tr><tr><td>ROC AUC</td><td>0.80 ± 0.13</td><td><strong>0.93 ± 0.07</strong></td><td>0.86 ± 0.12</td><td><strong>0.92 ± 0.07</strong></td></tr><tr><td>F1-score</td><td>0.63 ± 0.15</td><td><strong>0.82 ± 0.13</strong></td><td>0.78 ± 0.14</td><td><strong>0.83 ± 0.13</strong></td></tr><tr><td>Precision</td><td>0.68 ± 0.19</td><td><strong>0.87 ± 0.15</strong></td><td>0.79 ± 0.18</td><td><strong>0.86 ± 0.17</strong></td></tr><tr><td>Sensitivity</td><td>0.69 ± 0.22</td><td><strong>0.82 ± 0.18</strong></td><td>0.81 ± 0.18</td><td><strong>0.84 ± 0.18</strong></td></tr><tr><td>Specificity</td><td>0.75 ± 0.17</td><td><strong>0.90 ± 0.11</strong></td><td>0.83 ± 0.16</td><td><strong>0.90 ± 0.11</strong></td></tr></tbody></table>
accuracy and sensitivity which can be related to the EEG features used to train the model,
opening the question of which is the feature of the EEG signal better suited for training an AI
model for diagnosis. Until now, most of the studies attempting to perform AI-based diagnosis
exploited functional connectivity in the sensor space, derived by the covariation between
signals and/or spectral power from electrodes activity, as well as the frequency-specific phase
locking (i.e. synchronization) <sup>19</sup>. However, such techniques rely on the assumption of
stationarity of the EEG signals, which may represent an important limitation, and might affect
the stability of the estimates. In the present work, we aimed at identifying a better suited EEG
signal feature to train an AI-model-based classification (patients with temporal lobe epilepsy
vs. healthy controls). We hypothesized that the aperiodic activation bursts (neuronal
avalanches) of source-derived EEG activities may represent a better measure for an AI-model,
since they are not based on the assumptions of stationarity. In particular, neuronal avalanches
(NA) belong to the framework of criticality, which lends itself nicely to the investigation of the
microscopic dynamical behavior of neuronal assemblies and its (non-linear) relation to large-
scale alterations such as, in the case of epilepsy, seizure initiation, propagation and termination
<sup>20,21</sup>. Recent findings suggest that NA spreading on the large scale represents a more stable and
reliable feature for the individualized investigation and task classification of brain dynamics
<sup>22,23</sup>. Relevantly, our previous work has highlighted the altered spreading of NA and its relation
to brain morphology in temporal lobe epilepsy (TLE), using the source-derived high-density
EEG (hdEEG) activity at rest <sup>24,25</sup>. In the light of this, here we hypothesized that the spreading
of the aperiodic activation bursts at the whole-brain level may represent a relevant feature
towards the automatic classification of patients with TLE vs. controls. To this end, we recorded
resting-state high-density EEG (hdEEG) from patients with temporal lobe epilepsy and a
control group. The model of TLE is of particular interest as it represents a rather well-defined
group of electroclinical conditions with lower clinical heterogeneity compared to other epilepsy
forms. Using source-reconstructed EEG signals, we individuated the presence of NA and used
them to calculate the probability of consecutive activations between brain regions (i.e., the
topographical spreading of the avalanches) (see Supplementary Figure 1). This information
was stored within an adjacency matrix named avalanche transition matrix (ATM) <sup>26</sup> which was
then used as a feature for a support vector machine (SVM) classifier. We hypothesized that
ATMs would have performed better, compared to the power-based connectivity metric
generally used in literature, differentiating with a larger accuracy the two conditions (epilepsy
vs. non-epilepsy). Relevantly, we computed the ATM from a resting state activity free from
epileptiform activity, in order to investigate if the basal functional organization of the system,
Classification Performance

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Accuracy and receiver operating characteristic curves. Panel (A) of the present picture displays the full distribution of the classification accuracy of the support vector machine (SVM) model across the cross-validation splits, when trained with avalanche transition matrix (ATM, in salmon) and imaginary coherence (ImCoh, in purple). On the left the accuracy classification using narrow-band filtered signal (3–14 Hz) and on the right the broadband signal (3–40 Hz). Panel (B) shows the mean receiver operating characteristic (ROC) curves across cross-validation splits and the corresponding area under the curve (AUC), both for ATM (salmon line) and ImCoh (purple line).</figcaption></figure>

supplementary material, we also provide the ROC curves distribution in cross-validation splits for each of the frequency band (see Supplementary Fig. 3). However, in the manuscript, we display for simplicity the mean ROC across splits. To summarize, both broadband and narrow-band filtered signals display similar behavior, namely a better performance when using the aperiodic part of the signal as opposed to using a metric that assumes stationarity and is purely power-based. From here on we will limit ourselves to reporting the results in the broadband signal. Nonetheless, the results in the remaining frequency bands are reported in the Supplementary Materials.

Feature selected by the model

A relevant topic in automatic learning is related to interpretability, that is the ability to trace back the specific information that is being used to classify items. In our case, our input data was represented by a node-by-node matrix, where each edge indicates the power of the cross-covariance (ImCoh) or the transition probability (ATMs) between two nodes. We observed a difference in the features used by the model to classify between
without the presence of characteristic graphic-elements (i.e., seizure and/or interictal
epileptiform discharges), may contribute to a reliable classification of the patients. We then
looked into the most relevant feature pattern, namely the most informative elements used by
the model to differentiate the two groups. We expected that the dynamics of those areas linked
with the pathophysiology of TLE, would have informed the model more than other brain
regions. This endows our results with interpretability and allows the identification of clusters
of regional interactions that might be relevant to the pathology. Finally, we investigated the
classification performance with respect to the duration of the brain signals. Hence, we have
performed the classification based on the growing amount of data used to compute ATMs,
ranging from 5 to 300 sec. To summarize, the purpose of the present work is to exploit the
spreading of the aperiodic activities at the whole-brain level as a feature to inform an AI-based
model aimed at the automatic identification of epileptic patients, in order to provide the
clinicians with a practical and non-invasive tool supporting the diagnostic process.

**Results**

Model prediction performance

The implemented SVM model classified the participants as TLE or controls based on the global system dynamics, namely the ATMs, computed on the broadband signal (3-40Hz), with a mean accuracy of (0.86 ± 0.08). We observed that the prediction accuracy of the model trained on the ATMs was larger as compared to the one trained with ImCoh (0.71 ± 0.12). These results were also confirmed for the signals filtered in the narrowband (3-14 Hz) (ATMs accuracy: 0.87 ± 0.08; ImCoh accuracy: 0.82 ± 0.11). The better classification performance for ATMs vs. ImCoh was also confirmed by an increased area under the curve (AUC), both when considering the broadband and the narrowband filtered signal. The improvement of classification related to the ATMs is also confirmed when considering multiple frequency bands (i.e., alpha-beta, beta-low-gamma and alpha-beta-low-gamma; see Supplementary Table 1). In general, the results display a better classification across all the performance indices in the model trained with the ATMs as compared to the ImCoh (see Table.1).
TLE patients and healthy controls between the functional measures we used. When considering the ImCoh, we observed that the distribution of the relevant features is narrow, which means that most of the edges hold similar importance for the classification. Conversely, we observed a larger distribution of the importance of the features in the ATMs-based classification, which demonstrates that a specific subset of edges carries the information that is used to differentiate the two groups (see Fig. 2A). This difference is also visualized in the chord diagrams in Fig. 2B. One can observe similar importance values across most of the connections when considering ImCoh while, in the case of the ATMs, specific edges emerge from the background. The edge's relevance appeared to be clustered into specific nodes, which were different between ATM and ImCoh. On the one hand, when considering the model trained with the ATMs, the features with the largest importance were clustered around brain regions of fundamental relevance in TLE, involved both in seizure initiation and/or propagation, or characterized by well-known metabolic or connectivity alterations, such as the entorhinal cortex, the superior and inferior temporal gyri, the cingulate cortex and the prefrontal dorsolateral cortex<sup>28-31</sup> (see Fig. 2C). On the other hand, the model with ImCoh highlighted as more important the edge's clustering on pre- and postcentral gyri and rostral middle frontal regions. Similar results are also present when considering narrow-band filtered signal (see Supplementary Fig. 4).

Model Interpretability

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Feature importance for model interpretability. The present figure represents the relevance of the information used by the model to perform classification, namely the feature importance, both using avalanche transition matrix (ATM) and imaginary coherence (ImCoh). Panel (A) is a histogram of the probability of the feature importance value in ATM and ImCoh. The histogram shows a narrow distribution for ImCoh and a broader distribution for ATM, suggesting that in ATM certain edges drive the majority of information necessary for differentiating the two groups. Panel (B) is the edge representation in a chord plot, showing the importance of each edge in the classification. Finally, panel (C) is the mean importance value of each edge of a specific brain region. This representation highlights which regions mainly impact the classification.</figcaption></figure>
<table><thead><tr><th rowspan="2"></th><th colspan="2">Broadband (3-40Hz)<br>mean ± SD</th><th colspan="2">Narrowband (3-14Hz)<br>mean ± SD</th></tr><tr><th>ImCoh+SVM</th><th>ATM+SVM</th><th>ImCoh+SVM</th><th>ATM+SVM</th></tr></thead><tbody><tr><td>Accuracy</td><td>0.71 ± 0.12</td><td><strong>0.86 ± 0.08</strong></td><td>0.82 ± 0.11</td><td><strong>0.87 ± 0.08</strong></td></tr><tr><td>ROC AUC</td><td>0.80 ± 0.13</td><td><strong>0.93 ± 0.07</strong></td><td>0.86 ± 0.12</td><td><strong>0.92 ± 0.07</strong></td></tr><tr><td>F1-score</td><td>0.63 ± 0.15</td><td><strong>0.82 ± 0.13</strong></td><td>0.78 ± 0.14</td><td><strong>0.83 ± 0.13</strong></td></tr><tr><td>Precision</td><td>0.68 ± 0.19</td><td><strong>0.87 ± 0.15</strong></td><td>0.79 ± 0.18</td><td><strong>0.86 ± 0.17</strong></td></tr><tr><td>Sensitivity</td><td>0.69 ± 0.22</td><td><strong>0.82 ± 0.18</strong></td><td>0.81 ± 0.18</td><td><strong>0.84 ± 0.18</strong></td></tr><tr><td>Specificity</td><td>0.75 ± 0.17</td><td><strong>0.90 ± 0.11</strong></td><td>0.83 ± 0.16</td><td><strong>0.90 ± 0.11</strong></td></tr></tbody></table>

**Table 1. Classification performance.** The present table displays the indices of the classification performance when considering the model trained with avalanche transition matrix (ATM) or imaginary coherence (ImCoh) obtained from the narrowband (3-14Hz) or broadband (3-40Hz) signal.

We note that, as suggested by recent findings, providing exclusively the mean across splits of the classification accuracies might be misleading for the interpretation of the model performance 27. Hence, here we display the full distribution of the classification accuracies across the multiple splits of the cross-validation. Notably, the distribution of the prediction accuracies with the ATMs is narrower, and skewed towards larger accuracy values (see Figure 1A). These result are confirmed across all the frequency bands (see Supplementary Figure 2). Additionally, in Figure 1B, we also provide the mean receiver operating characteristic (ROC) curve and the corresponding AUC, providing further evidence of the increased performance of the ATMs compared to the ImCoh. In fact, in the supplementary material, we also provide the ROC curves distribution in cross-validation splits for each of the frequency band (see Supplementary Figure 3). However, in the manuscript, we display for simplicity the mean ROC across splits. To summarize, both broadband and narrow-band filtered signals display similar behavior, namely a better performance when using the aperiodic part of the signal as opposed to using a metric that assumes stationarity and is purely power-based. From here on we will
## Time scale of classification performance

We then focused on classification performance and the amount of signal used to compute the ATMs. For that purpose, we cut the EEG recordings into segments of different lengths (namely 5, 15, 30, 60, 120, 180, and 300 s). Such segments refer to the length of the chunk used for the extraction and the classification of the features. To limit the effect of the location of the considered time window ω in the recording, we randomly picked 100 times ω across the registration, and we computed the median value of the classification accuracy obtained over the permutations. As expected, the average classification performance increases as a function of the length of the segments. However, unlike the average, we note that the variance (across splits) of the accuracy does not become monotonically smaller as a function of the lengths of the segments. Instead, as one can observe in Fig. 3, the variance of the classification is nearly flat for segments of 30 s and of 180 s. However, when one moves from 30 to 60 s, one can observe that the variance increases, since higher accuracies, that were never achieved with segments, can be reached in some splits. In other words, the classification performance saturates from 5 to 30 s and then starts saturating again from 30 to 180 s (see Fig. 3).

## Discussion

The present work aimed at improving the automatable classification of epileptic patients, utilizing machine learning on information of the brain dynamics derived no- invasively from hdEEG. To this purpose, we used the avalanche transition matrix (ATM) in order to measure how large-scale aperiodic perturbations spread across the brain. In particular, the ATMs are a recently developed analytical tool that captures the probability of any two regions being successively recruited by spontaneous neuronal avalanches, resulting in a functional network<sup>26,32</sup>. Recently, we provided evidence that the ATMs are useful in temporal lobe epilepsy (TLE), as they allow the identification, from resting state hdEEG, of the regions with functional alterations, which are the ones generally involved in seizure initiation and propagation<sup>24</sup>. As a first finding, the model classified subjects as patients with epilepsy or controls, with a mean accuracy of 0.86 ± 0.08 (SD), with an AUC of 0.93 ± 0.07 (SD) and specificity of 0.90 ± 0.11 (SD) when using the ATMs-derived features. We observed that the model trained with ATMs outperforms the ImCoh across all the classification metrics increasing the performance from 10 to 19% based on the metric of interest. The implications of this finding are twofold. From a theoretical perspective, the information conveyed by the large wave of aperiodic bursts is more informative and representative of the altered brain dynamics in TLE, as compared to the power-covariance-based connectivity measure (ImCoh). This can also be appreciated when observing the type of information that the model is using to perform classification. In the ATM-model, subsets of specific connections between brain nodes mainly drive the classification. Importantly, these edges are clustered into regions of fundamental relevance in TLE, which are typically involved both in seizure initiation and/or propagation, or characterized by well-known metabolic or connectivity alterations, such as the entorhinal cortex, the superior and inferior temporal gyri, the cingulate cortex, and the prefrontal dorsolateral cortex<sup>28-31</sup>. These findings suggest that the ATMs are a sensible measure, capable of capturing altered neural dynamics and characterizing the TLE patient. The better performance obtained with ATM as compared to the one with ImCoh, both in the broad- and narrow-band filtered signal, might suggest that the aperiodic components might be of relevance when focusing on the large-scale dynamics. In other words, the alterations are not limited to a single temporal scale, but rather affect multiple timescales. This is in line with emerging literature showing that multiple nested characteristic timescales characterize the brain dynamics at rest<sup>33</sup>. Along this line of evidence, the alterations induced by TLE might also affect multiple timescales. In fact, the classification accuracy saturates at specific durations of the time series, while it shows larger variations (across the particular segment taken into account for a patient) for intermediate lengths. In other words, for a length of 30 s, which particular segment has been used to perform the classification does not have any relevance, as the classification will have converged to a specific value. When one moves to segments of 60 s, the situation is different, since the particular segment that is used for classification has an effect (the performance varies from one particular split

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Time dependence of classification accuracy. The picture shows the distribution across splits of the cross-validation based on the amount of signal used to compute avalanche transition matrix (ATM). Each dot corresponds to the estimation of the classification performance per split.</figcaption></figure>
limit ourselves to reporting the results in the broadband signal. Nonetheless, the results in the remaining frequency bands are reported in the Supplementary Materials.

Classification Performance

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Accuracy and receiver operating characteristic curves. Panel (A) of the present picture displays the full distribution of the classification accuracy of the support vector machine (SVM) model across the cross-validation splits, when trained with avalanche transition matrix (ATM, in salmon) and imaginary coherence (ImCoh, in purple). On the left the accuracy classification using narrow-band filtered signal (3-14 Hz) and on the right the broadband signal (3-40 Hz). Panel (B) shows the mean receiver operating characteristic (ROC) curves across</figcaption></figure>
to the next), which results in the variance of the classification performance across the segments. However, the average performance tends to rise as a function of the duration of the data segments, and never goes below the classification performance of the previous convergence time-point (i.e., 30 s). At 120 s, the variance is similar, but the results of the classification start converging towards a new upper value, which is reached at 180 s. Indeed, at this duration, the classification performance converges again, and it will reach a specific value regardless of the specific segment. We interpret these findings as an effect of the interplay between multiple nested timescales and the selected window length. Intuitively, the window length imposes an upper bound on the frequencies that can fit. If all time-scales were present, one would expect that widening the window-length would yield a proportional improvement in the classification performance. Our results show that there are specific window lengths where all the available information has been captured, and no further improvement is possible unless longer time windows are considered. We interpret this as a sign that widening the window allows for capturing activities evolving over slower timescales, which would simply be overlooked utilizing shorter time windows. For this reason, in order to exhaustively sample the brain dynamics complexity and exploit it for an efficient classification, longer segments are required. Further studies should elucidate whether, considering longer durations yet, allows gaining information from slower timescales yet (e.g. the well-known circadian or monthly oscillations in epileptic patients<sup>34–36</sup>).

Furthermore, these results have a potential direct application in clinical practice. In fact, the rate of epilepsy misdiagnosis varies between 20 and 40%⁴ and remains an issue also in specialized centers². Some populations carry a particularly high risk of misdiagnosis, for example, patients with psychogenic non-epileptic seizures³⁷. Misidentification of EEG abnormalities reaches over 50% in the PNES population³⁸. The method proposed in the present work may represent a potential tool for helping clinicians with diagnostics, which translates into the possibility of reducing the misinterpretation of EEG abnormalities by taking into account global brain activities. In fact, the present pipeline can be fully automated, requiring therefore minimal effort to the clinicians and the technical staff, that is the recording of a minimum of 5 min of resting state EEG. Importantly, the data have been analyzed without including epileptiform activities, which increases the generalizability of this analytical framework. In fact, in clinical practice, the identification of seizures and IEDs is a time-consuming process often performed manually. Moreover, seizures and/or IEDs are not always recorded, even during multiple days of monitoring. In this scenario, the possibility of classifying the patients using the intrinsic brain organization, without the necessity of epileptiform graphic elements, represents a strength of this approach. Nonetheless, the present work does not aim at providing a tool to replace the fundamental knowledge and experience of a clinical expert. In fact, intensive work has been dedicated in literature in the definition and characterization of easy-to-use methods as AI in IEDs and epilepsy diagnosis⁹,¹⁰,¹⁷,¹⁸ On the contrary, this workflow should be considered as a potential user-friendly pipeline in support of clinicians in the diagnosis process, leveraging the time resolution and non-invasiveness of the EEG. The present pipeline might be of important help also in the low-income countries where patients have difficulty accessing specialist care³⁹. The manageability and cost-effectiveness of the EEG, together with the present workflow, may represent a suitable instrument for helping clinicians in their daily practice. In order to translate this pipeline to clinical practice, validation on multiple epilepsy and neurological populations is warranted. In fact, in this work we trained of our model only with patients with TLE, in line with the data availability of our center, limiting its generalizability. Additional studies may apply the present workflow in order to characterize different types of epilepsy, potentially identifying the backbone of altered dynamics characterizing different clinical conditions. Additionally, another future line of work will consist in considering features selection strategies to ensure that only the most relevant information will be taken into account by the classifier. To summarize, our results provide evidence that the ATMs better capture the subtle neural dynamics alteration in epilepsy, as compared to classic functional connectivity measures, and that this may be effectively fed to artificial intelligence algorithms to improve automated classification of TLE patients.

# Method

## Participants

We retrospectively enrolled 70 patients with temporal lobe epilepsy, who underwent high-density electroencephalography (hdEEG) for clinical evaluation between 2018 and 2021 at the Epilepsy and Clinical Neurophysiology Unit, IRCCS Eugenio Medea cited in Conegliano (Italy). The diagnostic workflow included clinical history and examination, neuropsychological assessment, long-term surface Video EEG (32 channels) monitoring, hdEEG resting-state recording, 1.5/3 T brain magnetic resonance imaging, and positron emission tomography (PET) as an adjunctive investigation in selected cases. The diagnosis of temporal lobe epilepsy was established according to the ILAE guidelines. We selected for this specific work a subgroup of 31 patients with left temporal lobe epilepsy (mean age = 37.18 [SD = 18.04]; 14 females). The control group was composed of 31 healthy participants with no history of neurological or psychiatric disorders (mean age = 34.92 [SD = 9.22]; 25 females). A description of patients' demographic and clinical characteristics is provided in Table 2. The study protocol was conducted according to the Declaration of Helsinki and approved by the Ethics Committees for Clinical Trials (CESC-Treviso; protocol n.1309/CE-Medea). Informed consent was obtained from all subjects.

## Resting state EEG recording

The hdEEG recordings were obtained using a 128-channel Micromed system referenced to the vertex. Data was sampled at 1,024 Hz and the impedance was kept below 5 kΩ for each sensor. For each participant, we recorded 10 min of closed-eyes resting state while comfortably sitting on a chair in a silent room.
cross-validation splits and the corresponding area under the curve (AUC), both for ATM (salmon line) and ImCoh (purple line).

Feature selected by the model

A relevant topic in automatic learning is related to interpretability, that is the ability to trace back the specific information that is being used to classify items. In our case, our input data was represented by a node-by-node matrix, where each edge indicates the power of the cross-covariance (ImCoh) or the transition probability (ATMs) between two nodes. We observed a difference in the features used by the model to classify between TLE patients and healthy controls between the functional measures we used. When considering the ImCoh, we observed that the distribution of the relevant features is narrow, which means that most of the edges hold similar importance for the classification. Conversely, we observed a larger distribution of the importance of the features in the ATMs-based classification, which demonstrates that a specific subset of edges carries the information that is used to differentiate the two groups (see Figure 2A). This difference is also visualized in the chord diagrams in Figure 2B. One can observe similar importance values across most of the connections when considering ImCoh while, in the case of the ATMs, specific edges emerge from the background. The edge's relevance appeared to be clustered into specific nodes, which were different between ATM and ImCoh. On the one hand, when considering the model trained with the ATMs, the features with the largest importance were clustered around brain regions of fundamental relevance in TLE, involved both in seizure initiation and/or propagation, or characterized by well-known metabolic or connectivity alterations, such as the entorhinal cortex, the superior and inferior temporal gyri, the cingulate cortex and the prefrontal dorsolateral cortex <sup>28-31</sup> (see Figure 2C). On the other hand, the model with ImCoh highlighted as more important the edge's clustering on pre- and postcentral gyri and rostral middle frontal regions. Similar results are also present when considering narrow-band filtered signal (see Supplementary Figure 4).
<table><caption>**Table 2.** Sample demographic and clinical information. The table describes the demographic and clinical characteristics of the patients with temporal lobe epilepsy. MRI abnormalities are reported by sublobar localization. The continuous variables are reported as mean ± standard deviation. Antiseizure Medication abbreviations: *BRV* brivaracetam, *CBZ* carbamazepine, *CLB* clobazam, *CZP* clonazepam, *ESL* eslicarbazepine, *LCM* lacosamide, *LEV* levetiracetam, *LTG* lamotrigine, *OXC* oxcarbazepine, *PB* phenobarbital, *PER* perampanel, *VPA* valproic acid, *NO-ASMs* no pharmacological treatment. Abbreviation of the identified anomalies on the magnetic resonance imaging: *FCD* focal cortical dysplasia, *HS* hippocampal sclerosis, *DNET* dysembryoplastic neuroepithelial tumors, *UKN* unknown.</caption><thead><tr><th></th><th>Mean ± Standard deviation</th></tr></thead><tbody><tr><td>Age</td><td>37.18 ± 18.04</td></tr><tr><td>Age of onset</td><td>23.75 ± 16.06</td></tr><tr><td>Duration of Epilepsy (years)</td><td>13.57 ± 16.00</td></tr><tr><td>Number of Antiseizure medications</td><td>1.59 ± .97</td></tr><tr><td><strong>Antiseizure medications</strong></td><td><strong>Number</strong></td></tr><tr><td>BRV</td><td>3</td></tr><tr><td>CBZ</td><td>4</td></tr><tr><td>CLB</td><td>2</td></tr><tr><td>CZP</td><td>1</td></tr><tr><td>ESL</td><td>5</td></tr><tr><td>LCM</td><td>6</td></tr><tr><td>LEV</td><td>4</td></tr><tr><td>LTG</td><td>4</td></tr><tr><td>OXC</td><td>1</td></tr><tr><td>PER</td><td>5</td></tr><tr><td>VPA</td><td>5</td></tr><tr><td>NO-ASMs</td><td>5</td></tr><tr><td><strong>MRI</strong></td><td></td></tr><tr><td><strong>Mesial</strong></td><td><strong>Number</strong></td></tr><tr><td>HS</td><td>4</td></tr><tr><td>DNET</td><td>1</td></tr><tr><td>UKN</td><td>2</td></tr><tr><td>Amygdala enlargement</td><td>4</td></tr><tr><td><strong>Anterior (temporal pole)</strong></td><td></td></tr><tr><td>FCD</td><td>4</td></tr><tr><td>Encephalocele</td><td>2</td></tr><tr><td>Gliosis</td><td>2</td></tr><tr><td><strong>Anterior + mesial</strong></td><td></td></tr><tr><td>FCD + HS</td><td>2</td></tr><tr><td>Negative MRI</td><td>8</td></tr></tbody></table>

## EEG pre-processing

Signal preprocessing was performed via EEGLAB 14.1.2b⁴⁰. The continuous EEG signal was first downsampled at 250 Hz and then bandpass-filtered (0.1–45 Hz) using a Hamming windowed sinc finite impulse response filter (filter order = 8250). The signal was visually inspected to identify interictal epileptiform discharges (IEDs) by the clinicians and then segmented into 1-s-long epochs. Epochs containing IEDs activities were removed. We purposely removed epochs containing IEDs since we wanted to use the intrinsic brain functional organization of the brain, independent from epileptiform activity, for patients vs. control classification. At least one IED was recorded in the 83.8% of patients with epilepsy, while no spikes were detected in the control group.

The epoched data underwent an automated bad-channel and artifact detection algorithm using the TBT plugin implemented in EEGLAB. This algorithm identified the channels that exceeded a differential average amplitude of 250 µV and marked those channels for rejection. Channels that were marked as bad in more than 30% of all epochs were excluded. Additionally, epochs having more than 10 bad channels were excluded. We automatically detected possible flat channels with the Trimoutlier EEGLAB plug-in within the lower bound of 1 µV. We rejected an average of 40.15 ± 41.24 (SD) epochs. The preprocessing analysis pipeline has been applied by our group in previous studies investigating neuronal avalanches in epilepsy with resting state EEG activity²⁴,²⁵. Data cleaning was performed with independent component analysis, using the Infomax algorithm as implemented in EEGLAB. The resulting 40 independent components were visually inspected and those related to eye blinks, eye movements, muscle, and cardiac artifacts were discarded. An average of 9.51 ± 4.68 (SD) components were removed. The remaining components were then projected back to the electrode space. Finally, bad channels were reconstructed with the spherical spline interpolation method. An average of 2 channels (range
Model Interpretability

<figure><img src="image_3.png" /><figcaption>**Figure 2. Feature importance for model interpretability.** The present figure represents the relevance of the information used by the model to perform classification, namely the feature importance, both using avalanche transition matrix (ATM) and imaginary coherence (ImCoh). Panel (A) is a histogram of the probability of the feature importance value in ATM and ImCoh. The histogram shows a narrow distribution for ImCoh and a broader distribution for ATM, suggesting that in ATM certain edges drive the majority of information necessary for differentiating the two groups. Panel (B) is the edge representation in a chord plot, showing the importance of each edge in the classification. Finally, panel (C) is the mean importance value of each edge of a specific brain region. This representation highlights which regions mainly impact the classification.</figcaption></figure>
0–4) was interpolated in the sample, with no preferential spatial location across electrodes positions. The data were then re-referenced to the average of all electrodes. At the end of the data preprocessing, each subject with epilepsy had at least 5 min of artifact-free signal, while healthy controls had at least 7 min of artifact-free signal. To avoid bias in the classification due to the difference in the signal length, we used 5 min of artifact free signal for each subject in the following analyses.

## Cortical source modeling

We used the individual anatomy MRI to generate individualized head models for the patients with TLE. The anatomic MRI for source imaging consisted of a T1 isotropic three-dimensional (3D) acquisition. For twelve patients and for the control group we used the MNI-ICBM152 default anatomy⁴¹ from Brainstorm⁴² since the 3D T1 MRI sequences were not available. The MRI was segmented into skin, skull, and gray matter using the Computational Anatomy Toolbox (CAT12)⁴³. The resulting individual surfaces were then imported in Brainstorm, where three individual surfaces adapted for Boundary Element Models (BEM) were reconstructed (inner skull, outer skull and head) and the cortical mesh was downsampled at 15,002 vertices. The co-registration of the EEG electrodes was performed using Brainstorm by projecting the EEG sensor positions on the head surface with respect to the fiducial points of the individual or the template MRI. We applied manual correction of the EEG cap on the individual anatomy whenever needed, prior to projecting the electrodes on the individual head surface. We then derived an EEG forward model using the 3-shell BEM model estimated using OpenMEEG method implemented in Brainstorm⁴⁴. Finally, we used the weighted minimum norm imaging⁴⁵ as the inverse model, with the Brainstorm's default parameters setting.

## Features extraction

### Avalanche transition matrix computation

Similarly to our previous work using neuronal avalanches (Duma et al. 2023), we extracted the activity of a total of 68 regions of interest (ROIs) from the Desikan-Killiany atlas⁴⁶. The ROIs time series were obtained by averaging the activity across the vertices composing each ROI. To study the dynamics of brain activity, we estimated “neuronal avalanches” from the source-reconstructed ROI time series. Firstly, the time series of each ROI was discretized by calculating the z-score over time and then setting positive and negative excursions beyond a threshold as 1, and the rest of the signal as 0⁴⁷. A neuronal avalanche begins when, in a sequence of contiguous time bins, at least one ROI is active (i.e., above threshold), and ends when all ROIs are inactive²⁶,⁴⁸. These analyses require the time series to be binned. This is done to ensure that one is capturing critical dynamics, if present. To estimate the suitable time bin length, for each subject, each neuronal avalanche, and each time bin duration, the branching parameter σ was estimated⁴⁹. In fact, systems operating at criticality typically display a branching ratio ~ 1. The branching ratio is calculated as the geometrically averaged (over all the time bins) ratio of the number of events (activations) between the subsequent time bin (descendants) and that in the current time bin (ancestors), and then averaging over all the avalanches⁵⁰. More specifically:

$$ \sigma_i = \prod_{j=1}^{N_{bin}-1} \left( \frac{n_{events}(j+1)}{n_{events}(j)} \right)^{\frac{1}{N_{bin}-1}} \quad (1) $$

$$ \sigma = \prod_{i=1}^{N_{aval}} (\sigma_i)^{\frac{1}{N_{aval}}} \quad (2) $$

where σᵢ is the branching parameter of the ith avalanche in the dataset, Nₜₜₜₜₜ is the total amount of bins in the i-th avalanche, nᵣₚₚₚₚ is the total number of events active in the j-th bin, and Nₐₐₐₐₐ is the total number of avalanche in the dataset. In our analyses, the branching ratio was 1 for bin = 1 (corresponding to bins of 4 ms). An avalanche-specific transition matrix (ATM) was calculated where element (i, j) represented the probability that region j was active at time t + δ, given that region i was active at time t, where δ ~ 4 ms. The ATMs were averaged element-wise across all avalanches per each participant. It is important to mention that the ATMs were computed over the broadband signal. However, given the relevance of the slow frequencies (i.e., theta and alpha bands) in the TLE⁵¹,⁵², we also computed ATMs on a narrowband filtered signal (theta and alpha bands; 3–14 Hz). Moreover, we performed classification in additional consecutive frequency bands (i.e., alpha-beta, beta-low-gamma and alpha-beta-low-gamma). This was performed to evaluate the potential changes of the accuracy of classification when including only specific frequency bands. For each tested frequency band, we optimized two parameters: the threshold applied to the z-scored signals (ranging from 1.2 to 3.0) and the minimal duration of the avalanches considered (ranging from 2 to 8 time points). The chosen judgment criterion was the classification accuracy rate reached by the chosen configuration. Finally, we computed the ATMs multiple times, each time averaging over different lengths of the signals (i.e., 5, 15, 30, 60, 120, 180, and 300 s). This was performed in order to investigate the effect of signal length used in the ATMs computation on the classification performance. To mitigate the influence of the position of a given temporal window ω in the registration on the classification performance, we randomly picked 100 times ω across the registration, and we considered the median value of the classification accuracy over the permutations.

### Imaginary coherence

In addition to the ATMs, we computed the imaginary coherence in order to compare the performance with a widely used metric of functional connectivity. The imaginary part of the coherence (ImCoh) has been
## Time scale of classification performance

We then focused on classification performance and the amount of signal used to compute the ATMs. For that purpose, we cut the EEG recordings into segments of different lengths (namely 5, 15, 30, 60, 120, 180, and 300 seconds). Such segments refer to the length of the chunk used for the extraction and the classification of the features. To limit the effect of the location of the considered time window $\omega$ in the recording, we randomly picked 100 times $\omega$ across the registration, and we computed the median value of the classification accuracy obtained over the permutations. As expected, the average classification performance increases as a function of the length of the segments. However, unlike the average, we note that the variance (across splits) of the accuracy does not become monotonically smaller as a function of the lengths of the segments. Instead, as one can observe in Figure 3, the variance of the classification is nearly flat for segments of 30s and of 180s. However, when one moves from 30 to 60 seconds, one can observe that the variance increases, since higher accuracies, that were never achieved with segments, can be reached in some splits. In other words, the classification performance saturates from 5 to 30s and then starts saturating again from 30s to 180s (see Figure 3).

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Time dependence of classification accuracy. The picture shows the distribution across splits of the cross-validation based on the amount of signal used to compute avalanche transition matrix (ATM). Each dot corresponds to the estimation of the classification performance per split.</figcaption></figure>
Analysis Pipeline

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Analysis pipeline. The present figure illustrates the analytical steps used for feature extraction and machine learning based classification using a support vector machine approach.</figcaption></figure>

demonstrated to be robust to volume conduction, improving the identification of true interactions in EEG-derived functional connectivity⁵³,⁵⁴. We computed ImCoh as it follows:

$$ImCohxy(f) = \frac{Im(Sxy(f))}{\sqrt{Sxx(f)Syy(f)}} \qquad (3)$$

where S<sub>xy</sub> is the cross-spectral density between region x and region y, while S<sub>xx</sub> and S<sub>yy</sub> represent the autospectral density of the two regions. The spectral densities were estimated using a multitaper method with digital prolate spheroidal sequence (DPSS)⁵⁵ windows of 10 s, an overlap of 50%, and a maximum frequency resolution of 0.1 Hz. To be consistent with the scheme used with the ATMs, we calculated the ImCoh over both the broadband and the narrowband filtered signal (3–14 Hz) as well as the additional frequency bands (i.e., alpha-beta, beta-low-gamma and alpha-beta-low-gamma).
Discussion

The present work aimed at improving the automatable classification of epileptic patients,
utilizing machine learning on information of the brain dynamics derived no- invasively from
hdEEG. To this purpose, we used the avalanche transition matrix (ATM) in order to measure
how large-scale aperiodic perturbations spread across the brain. In particular, the ATMs are a
recently developed analytical tool that captures the probability of any two regions being
successively recruited by spontaneous neuronal avalanches, resulting in a functional network
<sup>26,32</sup>. Recently, we provided evidence that the ATMs are useful in temporal lobe epilepsy (TLE),
as they allow the identification, from resting state hdEEG, of the regions with functional
alterations, which are the ones generally involved in seizure initiation and propagation <sup>24</sup>. As a
first finding, the model classified subjects as patients with epilepsy or controls, with a mean
accuracy of 0.86 ± 0.08 (SD), with an AUC of 0.93 ± 0.07 (SD) and specificity of 0.90 ± 0.11
(SD) when using the ATMs-derived features. We observed that the model trained with ATMs
outperforms the ImCoh across all the classification metrics increasing the performance from
10 to 19% based on the metric of interest. The implications of this finding are twofold. From a
theoretical perspective, the information conveyed by the large wave of aperiodic bursts is more
informative and representative of the altered brain dynamics in TLE, as compared to the power-
covariance-based connectivity measure (ImCoh). This can also be appreciated when observing
the type of information that the model is using to perform classification. In the ATM-model,
subsets of specific connections between brain nodes mainly drive the classification.
Importantly, these edges are clustered into regions of fundamental relevance in TLE, which are
typically involved both in seizure initiation and/or propagation, or characterized by well-known
metabolic or connectivity alterations, such as the entorhinal cortex, the superior and inferior
temporal gyri, the cingulate cortex, and the prefrontal dorsolateral cortex <sup>28-31</sup>. These findings
suggest that the ATMs are a sensible measure, capable of capturing altered neural dynamics
and characterizing the TLE patient. The better performance obtained with ATM as compared
to the one with ImCoh, both in the broad- and narrow-band filtered signal, might suggest that
the aperiodic components might be of relevance when focusing on the large-scale dynamics.
In other words, the alterations are not limited to a single temporal scale, but rather affect
multiple timescales. This is in line with emerging literature showing that multiple nested
characteristic timescales characterize the brain dynamics at rest <sup>33</sup>. Along this line of evidence,
the alterations induced by TLE might also affect multiple timescales. In fact, the classification
accuracy saturates at specific durations of the time series, while it shows larger variations
(across the particular segment taken into account for a patient) for intermediate lengths. In other
## Classification pipeline

To perform the classification, we used a Support Vector Machine (SVM) classifier⁵⁶. We optimized the SVM with a fivefold cross-validated grid search to find optimal parameters, namely the kernels (linear or Radial Basis Function-RBF) and the regularization hyperparameter C during the training step. The kernel functions enable the system to work in a high-dimensional space and to allow for more complex decision functions when the data is not linearly separable⁵⁷. The C parameter adjusts the misclassification of training examples against the simplicity of the decision. The classification scores were evaluated by cross-validation consisting of a random permutation cross-validator with a training size of 80%, a test size of 20%, and 50 re-shuffling and splitting iterations. To avoid biased results, the data was divided across patients. Samples were first shuffled and then split into a pair of training and testing sets. Such an operation has been performed 50 times (corresponding to the number of splits we fixed). To ensure the reproducibility of our results, we set the randomness with a fixed pseudo random number generator. To assess the performance of the classification we reported several metrics, namely the accuracy classification score, the Area under the Receiver Operating Characteristic Curve (ROC AUC) from the prediction scores, the F1, the precision, the sensitivity, and the specificity scores. To investigate the interpretability of the classification performance, we studied the relative importance of the features, derived from the absolute value of the classification coefficients of the classification model. Here, we computed the median value over the splits of the cross-validation of the testing set. The analysis was conducted with Python. The code is publicly available at:

https://github.com/mccorsi/NeuronalAvalanches_TemporalLobeEpilepsy_EEG.

We performed classification training on the model with ATMs and ImCoh features (see Fig. 4 for the analysis pipeline representation).

## Data availability

The data that support the findings of this study are available on request to the corresponding author. The raw data are not publicly available due to privacy or ethical restrictions. All the scripts are available at the following github page: https://github.com/mccorsi/NeuronalAvalanches_TemporalLobeEpilepsy_EEG.

Received: 19 February 2024; Accepted: 13 June 2024

Published online: 18 June 2024

## References

1. Juarez-Garcia, A., Stokes, T., Shaw, B., Camosso-Stefinovic, J. & Baker, R. The costs of epilepsy misdiagnosis in England and Wales. *Seizure* 15(8), 598–605. https://doi.org/10.1016/j.seizure.2006.08.005 (2006).

2. Uldall, P., Alving, J., Hansen, L. K., Kibæk, M. & Buchholt, J. The misdiagnosis of epilepsy in children admitted to a tertiary epilepsy centre with paroxysmal events. *Arch. Dis. Child.* 91(3), 219–221. https://doi.org/10.1136/adc.2004.064477 (2006).

3. Zaidi, A., Clough, P., Cooper, P., Scheepers, B. & Fitzpatrick, A. P. Misdiagnosis of epilepsy: Many seizure-like attacks have a cardiovascular cause. *J. Am. Coll. Cardiol.* 36(1), 181–184. https://doi.org/10.1016/S0735-1097(00)00700-2 (2000).

4. Oto, M. The misdiagnosis of epilepsy: Appraising risks and managing uncertainty. *Seizure* 44, 143–146. https://doi.org/10.1016/j.seizure.2016.11.029 (2017).

5. Ilias, L., Askounis, D. & Psarras, J. Multimodal detection of epilepsy with deep neural networks. *Expert Syst. Appl.* 213, 119010. https://doi.org/10.1016/j.eswa.2022.119010 (2023).

6. Jin, L. et al. Multimodal image integration for epilepsy presurgical evaluation: A clinical workflow. *Front. Neurol.* 12, 709400. https://doi.org/10.3389/fneur.2021.709400 (2021).

7. Tatum, W. O. et al. Clinical utility of EEG in diagnosing and monitoring epilepsy in adults. *Clin. Neurophysiol.* 129(5), 1056–1082. https://doi.org/10.1016/j.clinph.2018.01.019 (2018).

8. Schmidt, H. et al. A computational biomarker of idiopathic generalized epilepsy from resting state EEG. *Epilepsia* 57(10), e200–e204. https://doi.org/10.1111/epi.13481 (2016).

9. Thangavel, P. et al. Improving automated diagnosis of epilepsy from EEGs beyond IEDs. *J. Neural. Eng.* 19(6), 066017. https://doi.org/10.1088/1741-2552/ac9c93 (2022).

10. Myers, P. et al. Diagnosing epilepsy with normal interictal EEG using dynamic network models. *medRxiv* https://doi.org/10.1101/2023.08.12.23294018 (2023).

11. Acharya, U. R., Hagiwara, Y. & Adeli, H. Automated seizure prediction. *Epilepsy Behav.* 88, 251–261. https://doi.org/10.1016/j.yebeh.2018.09.030 (2018).

12. Zhang, Y., Guo, Y., Yang, P., Chen, W. & Lo, B. Epilepsy seizure prediction on EEG using common spatial pattern and convolutional neural network. *IEEE J. Biomed. Health Inform.* 24(2), 465–474. https://doi.org/10.1109/JBHI.2019.2933046 (2020).

13. Zhou, M. et al. Epileptic seizure detection based on EEG signals and CNN. *Front. Neuroinform.* https://doi.org/10.3389/fninf.2018.00095 (2018).

14. Rijnders, B., Korkmaz, E. E. & Yildirim, F. Hybrid machine learning method for a connectivity-based epilepsy diagnosis with resting-state EEG. *Med. Biol. Eng. Comput.* 60(6), 1675–1689. https://doi.org/10.1007/s11517-022-02560-w (2022).

15. Varone, G. et al. A machine learning approach involving functional connectivity features to classify rest-EEG psychogenic non-epileptic seizures from healthy controls. *Sensors* 22(1), 129. https://doi.org/10.3390/s22010129 (2021).

16. Wang, B. et al. An integrative prediction algorithm of drug-refractory epilepsy based on combined clinical-EEG functional connectivity features. *J. Neurol.* 269(3), 1501–1514. https://doi.org/10.1007/s00415-021-10718-z (2022).

17. Jing, J. et al. Development of expert-level automated detection of epileptiform discharges during electroencephalogram interpretation. *JAMA Neurol.* 77(1), 103–108. https://doi.org/10.1001/jamaneurol.2019.3485 (2020).

18. Thomas, J. et al. Automated adult epilepsy diagnostic tool based on interictal scalp electroencephalogram characteristics: A six-center study. *Int. J. Neur. Syst.* 31(05), 2050074. https://doi.org/10.1142/S0129065720500744 (2021).

19. Saeidi, M. et al. Neural decoding of EEG signals with machine learning: A systematic review. *Brain Sci.* 11(11), 1525. https://doi.org/10.3390/brainsci11111525 (2021).

20. Liu, S., Li, F. & Wan, F. Distance to criticality undergoes critical transition before epileptic seizure attacks. *Brain Res. Bull.* 200, 110684. https://doi.org/10.1016/j.brainresbull.2023.110684 (2023).

21. Moosavi, S. A. & Truccolo, W. Criticality in probabilistic models of spreading dynamics in brain networks: Epileptic seizures. *PLOS Comput. Biol.* 19(2), e1010852. https://doi.org/10.1371/journal.pcbi.1010852 (2023).

22. Corsi, M. C. et al. Measuring neuronal avalanches to inform brain-computer interfaces. *iScience* 27(1), 108734. https://doi.org/10.1016/j.isci.2023.108734 (2024).
words, for a length of 30 s, which particular segment has been used to perform the classification
does not have any relevance, as the classification will have converged to a specific value. When
one moves to segments of 60s, the situation is different, since the particular segment that is
used for classification has an effect (the performance varies from one particular split to the
next), which results in the variance of the classification performance across the segments.
However, the average performance tends to rise as a function of the duration of the data
segments, and never goes below the classification performance of the previous convergence
time-point (i.e., 30s). At 120 seconds, the variance is similar, but the results of the classification
start converging towards a new upper value, which is reached at 180 seconds. Indeed, at this
duration, the classification performance converges again, and it will reach a specific value
regardless of the specific segment. We interpret these findings as an effect of the interplay
between multiple nested timescales and the selected window length. Intuitively, the window
length imposes an upper bound on the frequencies that can fit. If all time-scales were present,
one would expect that widening the window-length would yield a proportional improvement
in the classification performance. Our results show that there are specific window lengths
where all the available information has been captured, and no further improvement is possible
unless longer time windows are considered. We interpret this as a sign that widening the
window allows for capturing activities evolving over slower timescales, which would simply
be overlooked utilizing shorter time windows. For this reason, in order to exhaustively sample
the brain dynamics complexity and exploit it for an efficient classification, longer segments are
required. Further studies should elucidate whether, considering longer durations yet, allows
gaining information from slower timescales yet (e.g. the well-known circadian or monthly
oscillations in epileptic patients ³⁴⁻³⁶).

Furthermore, these results have a potential direct application in clinical practice. In fact, the
rate of epilepsy misdiagnosis varies between 20-40%⁴ and remains an issue also in specialized
centers². Some populations carry a particularly high risk of misdiagnosis, for example, patients
with psychogenic non-epileptic seizures³⁷. Misidentification of EEG abnormalities reaches
over 50% in the PNES population³⁸. The method proposed in the present work may represent
a potential tool for helping clinicians with diagnostics, which translates into the possibility of
reducing the misinterpretation of EEG abnormalities by taking into account global brain
activities. In fact, the present pipeline can be fully automated, requiring therefore minimal
effort to the clinicians and the technical staff, that is the recording of a minimum of 5 minutes
of resting state EEG. Importantly, the data have been analyzed without including epileptiform
activities, which increases the generalizability of this analytical framework. In fact, in clinical
23. Sorrentino, P. et al. Brain fingerprint is based on the aperiodic, scale-free, neuronal activity. *NeuroImage* **277**, 120260. https://doi.org/10.1016/j.neuroimage.2023.120260 (2023).

24. Duma, G. M. et al. Altered spreading of neuronal avalanches in temporal lobe epilepsy relates to cognitive performance: A resting-state hDEEG study. *Epilepsia* **64**(5), 1278-1288. https://doi.org/10.1111/epi.17551 (2023).

25. Duma, G. M. et al. Altered spread of waves of activities at large scale is influenced by cortical thickness organization in temporal lobe epilepsy: A MRI - high-density EEG study. *Brain Commun.* **2023**, 348. https://doi.org/10.1093/braincomms/fcad348 (2023).

26. Sorrentino, P. et al. Whole-brain propagation delays in multiple sclerosis, a combined tractography-magnetoencephalography study. *J. Neurosci.* **42**(47), 8807–8816. https://doi.org/10.1523/JNEUROSCI.0938-22.2022 (2022).

27. Shafiezadeh, S. et al. Methodological issues in evaluating machine learning models for EEG seizure prediction: Good cross-validation accuracy does not guarantee generalization to new patients. *Appl. Sci.* **13**(7), 4262. https://doi.org/10.3390/app13074262 (2023).

28. Bernasconi, N. et al. Entorhinal cortex in temporal lobe epilepsy: A quantitative MRI study. *Neurology* **52**(9), 1870–1870. https://doi.org/10.1212/WNL.52.9.1870 (1999).

29. Kubota, Y. et al. In vivo human hippocampal cingulate connectivity: A corticocortical evoked potentials (CCEPs) study. *Clin. Neurophysiol.* **124**(8), 1547–1556. https://doi.org/10.1016/j.clinph.2013.01.024 (2013).

30. Qin, L. et al. Alterations functional connectivity in temporal lobe epilepsy and their relationships with cognitive function: A longitudinal resting-state fMRI study. *Front. Neurol.* https://doi.org/10.3389/fneur.2020.00625 (2020).

31. Vismer, M. S., Forcelli, P. A., Skopin, M. D., Gale, K. & Koubеissi, M. Z. The piriform, perirhinal, and entorhinal cortex in seizure generation. *Front. Neural Circuits* https://doi.org/10.3389/fncir.2015.00027 (2015).

32. Rucco, R. et al. Neuronal avalanches to study the coordination of large-scale brain activity: Application to Rett syndrome. *Front Psychol.* **11**, 550749. https://doi.org/10.3389/fpsyg.2020.550749 (2020).

33. Sorrentino, P. et al. Dynamical interactions reconfigure the gradient of cortical timescales. *Netw. Neurosci.* **7**(1), 73–85. https://doi.org/10.1162/netn_a_00270 (2023).

34. Karoly, P. J. et al. Circadian and circaseptan rhythms in human epilepsy: A retrospective cohort study. *Lancet Neurol.* **17**(11), 977–985. https://doi.org/10.1016/S1474-4422(18)30274-6 (2018).

35. Khan, S. et al. Circadian rhythm and epilepsy. *Lancet Neurol.* **17**(12), 1098–1108. https://doi.org/10.1016/S1474-4422(18)30335-1 (2018).

36. Baud, M. O. et al. Multi-day rhythms modulate seizure risk in epilepsy. *Nat. Commun.* **9**(1), 88. https://doi.org/10.1038/s41467-017-02577-y (2018).

37. Xu, Y. et al. Frequency of a false positive diagnosis of epilepsy: A systematic review of observational studies. *Seizure Eur. J. Epilepsy* **41**, 167–174. https://doi.org/10.1016/j.seizure.2016.08.005 (2016).

38. Reuber, M., Fernández, G., Bauer, J., Singh, D. D. & Elger, C. E. Interictal EEG abnormalities in patients with psychogenic nonepileptic seizures. *Epilepsia* **43**(9), 1013–1020. https://doi.org/10.1046/j.1528-1157.2002.52301.x (2002).

39. Dorji, T. et al. Challenges in epilepsy diagnosis and management in a low-resource setting: An experience from Bhutan. *Epilepsy Res.* **192**, 107126. https://doi.org/10.1016/j.elepsyres.2023.107126 (2023).

40. Delorme, A. & Makeig, S. EEGLAB: An open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. *J. Neurosci. Methods* **134**(1), 9–21. https://doi.org/10.1016/j.jneumeth.2003.10.009 (2004).

41. Evans, A. C., Janke, A. L., Collins, D. L. & Baillet, S. Brain templates and atlases. *NeuroImage* **62**(2), 911–922. https://doi.org/10.1016/j.neuroimage.2012.01.024 (2012).

42. Tadel, F., Baillet, S., Mosher, J. C., Pantazis, D. & Leahy, R. M. Brainstorm: a user-friendly application for MEG/EEG analysis. *Comput. Intell. Neurosci.* **2011**, 1–13. https://doi.org/10.1155/2011/879716 (2011).

43. Dahnke, R., Yotter, R. A. & Gaser, C. Cortical thickness and central surface estimation. *NeuroImage* **65**, 336–348. https://doi.org/10.1016/j.neuroimage.2012.09.050 (2013).

44. Gramfort, A., Papadopoulos, T., Olivi, E. & Clerc, M. OpenMEEG: Opensource software for quasistatic bioelectromagnetics. *BioMed. Eng. OnLine* **9**(1), 45. https://doi.org/10.1186/1475-925X-9-45 (2010).

45. Hämäläinen, M. S. & Ilmoniemi, R. J. Interpreting magnetic fields of the brain: Minimum norm estimates. *Med. Biol. Eng. Comput.* **32**(1), 35–42. https://doi.org/10.1007/BF02512476 (1994).

46. Desikan, R. S. et al. An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest. *NeuroImage* **31**(3), 968–980. https://doi.org/10.1016/j.neuroimage.2006.01.021 (2006).

47. Sorrentino, P. et al. The structural connectome constrains fast brain dynamics. *eLife* **10**, e67400. https://doi.org/10.7554/eLife.67400 (2021).

48. Shriki, O. et al. Neuronal Avalanches in the resting MEG of the human brain. *J Neurosci.* **33**(16), 7079–7090. https://doi.org/10.1523/JNEUROSCI.4286-12.2013 (2013).

49. Haldeman, C. & Beggs, J. M. Critical branching captures activity in living neural networks and maximizes the number of metastable states. *Phys Rev Lett* **94**(5), 058101. https://doi.org/10.1103/PhysRevLett.94.058101 (2005).

50. Bak, P., Tang, C. & Wiesenfeld, K. Self-organized criticality: an explanation of the 1/f noise. *Phys. Rev. Lett.* **59**(4), 381–384. https://doi.org/10.1103/PhysRevLett.59.381 (1987).

51. Scherer, M. et al. Desynchronization of temporal lobe theta-band activity during effective anterior thalamus deep brain stimulation in epilepsy. *NeuroImage* **218**, 116967. https://doi.org/10.1016/j.neuroimage.2020.116967 (2020).

52. Sip, V., Scholly, J., Guye, M., Bartolomei, F. & Jirsa, V. Evidence for spreading seizure as a cause of theta-alpha activity electrographic pattern in stereo-EEG seizure recordings. *PLoS Comput. Biol.* **17**(2), e1008731. https://doi.org/10.1371/journal.pcbi.1008731 (2021).

53. Haufe, S. & Ewald, A. A simulation framework for benchmarking EEG-based brain connectivity estimation methodologies. *Brain Topogr.* **32**(4), 625–642. https://doi.org/10.1007/s10548-016-0498-y (2019).

54. Nolte, G. et al. Identifying true brain interaction from EEG data using the imaginary part of coherency. *Clin. Neurophys.* **115**(10), 2292–2307. https://doi.org/10.1016/j.clinphys.2004.04.029 (2004).

55. Moore, I. C. & Cada, M. Prolate spheroidal wave functions, an introduction to the Slepian series and its properties. *Appl. Comput. Harmon. Anal.* **16**(3), 208–230. https://doi.org/10.1016/j.acha.2004.03.004 (2004).

56. Smola, A. J. & Schölkopf, B. A tutorial on support vector regression. *Stat. Comput.* **14**(3), 199–222. https://doi.org/10.1023/B:STCO.0000035301.495498 (2004).

57. Faouzi, J. & Colliot, O. Classic Machine Learning Methods. In *Machine Learning for Brain Disorders* (ed. Colliot, O.) 25–75 (Neuromethods Springer, 2023). https://doi.org/10.1007/978-1-0716-3195-9_2

## Acknowledgements

This work was supported by Ricerca Corrente 2024 funds for biomedical research of The Italian Health Ministry and European Union "NextGenerationEU", (Investimento 3.1.M4. C2), project IR0000011, EBRAINS-Italy of PNRR.
practice, the identification of seizures and IEDs is a time-consuming process often performed manually. Moreover, seizures and/or IEDs are not always recorded, even during multiple days of monitoring. In this scenario, the possibility of classifying the patients using the intrinsic brain organization, without the necessity of epileptiform graphic elements, represents a strength of this approach. Nonetheless, the present work does not aim at providing a tool to replace the fundamental knowledge and experience of a clinical expert. In fact, intensive work has been dedicated in literature in the definition and characterization of easy-to-use methods as AI in IEDs and epilepsy diagnosis ⁹,¹⁰,¹⁷,¹⁸ On the contrary, this workflow should be considered as a potential user-friendly pipeline in support of clinicians in the diagnosis process, leveraging the time resolution and non-invasiveness of the EEG. The present pipeline might be of important help also in the low-income countries where patients have difficulty accessing specialist care ³⁹. The manageability and cost-effectiveness of the EEG, together with the present workflow, may represent a suitable instrument for helping clinicians in their daily practice. In order to translate this pipeline to clinical practice, validation on multiple epilepsy and neurological populations is warranted. In fact, in this work we trained of our model only with patients with TLE, in line with the data availability of our center, limiting its generalizability. Additional studies may apply the present workflow in order to characterize different types of epilepsy, potentially identifying the backbone of altered dynamics characterizing different clinical conditions. Additionally, another future line of work will consist in considering features selection strategies to ensure that only the most relevant information will be taken into account by the classifier. To summarize, our results provide evidence that the ATMs better capture the subtle neural dynamics alteration in epilepsy, as compared to classic functional connectivity measures, and that this may be effectively fed to artificial intelligence algorithms to improve automated classification of TLE patients.

# METHOD

## Participants

We retrospectively enrolled 70 patients with temporal lobe epilepsy, who underwent high-density electroencephalography (hdEEG) for clinical evaluation between 2018-2021 at the Epilepsy and Clinical Neurophysiology Unit, IRCCS Eugenio Medea cited in Conegliano (Italy). The diagnostic workflow included clinical history and examination, neuropsychological assessment, long-term surface Video EEG (32 channels) monitoring, hdEEG resting-state recording, 1.5/3T brain magnetic resonance imaging, and positron

a supprimé: to provide
**Author contributions**

Conceptualization: M.C.C., P.S., G.M.D., Methodology: M.C.C., E.T.L., P.S., S.C., Investigation: M.C.C., E.T.L., P.S., S.C., G.M.D., Visualization: M.C.C., G.M.D., Supervision: M.C.C., P.S., A.D., P.B., G.M.D., Writing—original draft: M.C.C., P.S., G.M.D., Writing—review & editing: M.C.C., E.T.L., P.S., S.C., A.D., P.B., G.M.D.

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-024-64870-3.

**Correspondence** and requests for materials should be addressed to M.-C.C. or P.S.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024
emission tomography (PET) as an adjunctive investigation in selected cases. The diagnosis of temporal lobe epilepsy was established according to the ILAE guidelines. We selected for this specific work a subgroup of 31 patients with left temporal lobe epilepsy (mean age = 37.18 [SD = 18.04]; 14 females). The control group was composed of 31 healthy participants with no history of neurological or psychiatric disorders (mean age = 34.92 [SD = 9.22]; 25 females). A description of patients' demographic and clinical characteristics is provided in Table 2. The study protocol was conducted according to the Declaration of Helsinki and approved by the local ethical committee. Informed consent was obtained from all subjects.

<table><thead><tr><th></th><th>Mean ± Standard deviation</th></tr></thead><tbody><tr><td><strong>Patients with TLE</strong></td><td></td></tr><tr><td>Age</td><td>37.18 ± 18.04</td></tr><tr><td>Age of onset</td><td>23.75 ± 16.06</td></tr><tr><td>Duration of Epilepsy (years)</td><td>13.57 ± 16.00</td></tr><tr><td>Number of Antiseizure Medications</td><td>1.59 ± .97</td></tr><tr><td><strong>Antiseizure Medications</strong></td><td><strong>Number</strong></td></tr><tr><td>BRV</td><td>3</td></tr><tr><td>CBZ</td><td>4</td></tr><tr><td>CLB</td><td>2</td></tr><tr><td>CZP</td><td>1</td></tr><tr><td>ESL</td><td>5</td></tr><tr><td>LCM</td><td>6</td></tr><tr><td>LEV</td><td>4</td></tr><tr><td>LTG</td><td>4</td></tr><tr><td>OXC</td><td>1</td></tr><tr><td>PER</td><td>5</td></tr><tr><td>VPA</td><td>5</td></tr><tr><td>NO-ASMs</td><td>5</td></tr><tr><td><strong>MRI</strong></td><td></td></tr><tr><td><strong>Mesial</strong></td><td><strong>Number</strong></td></tr><tr><td>HS</td><td>4</td></tr></tbody></table>
<table><tr><td>DNET</td><td>1</td></tr><tr><td>UKN</td><td>2</td></tr><tr><td>Amygdala enlargement</td><td>4</td></tr><tr><td><i>Anterior (temporal pole)</i></td><td></td></tr><tr><td>FCD</td><td>4</td></tr><tr><td>Encephalocele</td><td>2</td></tr><tr><td>Gliosis</td><td>2</td></tr><tr><td><i>Anterior + mesial</i></td><td></td></tr><tr><td>FCD+ HS</td><td>2</td></tr><tr><td>Negative MRI</td><td>8</td></tr></table>

**Table 2. Sample demographic and clinical information.** The table describes the demographic and clinical characteristics of the patients with temporal lobe epilepsy. MRI abnormalities are reported by sublobar localization. The continuous variables are reported as mean ± standard deviation. Antiseizure Medication abbreviations: BRV = brivaracetam, CBZ = carbamazepine, CLB = clobazam, CZP = clonazepam, ESL= eslicarbazepine, LCM = lacosamide, LEV = levetiracetam, LTG = lamotrigine, OXC = oxcarbazepine, PB = phenobarbital, PER = perampanel, VPA= valproic acid, NO-ASMs = no pharmacological treatment. Abbreviation of the identified anomalies on the magnetic resonance imaging: FCD= focal cortical dysplasia, HS=hippocampal sclerosis, DNET = dysembryoplastic neuroepithelial tumors, UKN = unknown.

## Resting State EEG recording

The hdEEG recordings were obtained using a 128-channel Micromed system referenced to the vertex. Data was sampled at 1,024 Hz and the impedance was kept below 5kΩ for each sensor. For each participant, we recorded 10 minutes of closed-eyes resting state while comfortably sitting on a chair in a silent room.

## EEG pre-processing

Signal preprocessing was performed via EEGLAB 14.1.2b⁴⁰. The continuous EEG signal was first downsampled at 250 Hz and then bandpass-filtered (0.1 to 45 Hz) using a Hamming windowed sinc finite impulse response filter (filter order = 8250). The signal was visually
inspected to identify interictal epileptiform discharges (IEDs) by the clinicians and then
segmented into 1-sec-long epochs. Epochs containing IEDs activities were removed. We
purposely removed epochs containing IEDs since we wanted to use the intrinsic brain
functional organization of the brain, independent from epileptiform activity, for patients vs.
control classification. At least one IED was recorded in the 83.8% of patients with epilepsy,
while no spikes were detected in the control group.

The epoched data underwent an automated bad-channel and artifact detection algorithm using
the TBT plugin implemented in EEGLAB. This algorithm identified the channels that exceeded
a differential average amplitude of 250µV and marked those channels for rejection. Channels
that were marked as bad in more than 30% of all epochs were excluded. Additionally, epochs
having more than 10 bad channels were excluded. We automatically detected possible flat
channels with the Trimoutlier EEGLAB plug-in within the lower bound of 1µV. We rejected
an average of 40.15 ± 41.24 (SD) epochs. The preprocessing analysis pipeline has been applied
by our group in previous studies investigating neuronal avalanches in epilepsy with resting
state EEG activity <sup>24,25</sup>. Data cleaning was performed with independent component analysis,
using the Infomax algorithm as implemented in EEGLAB. The resulting 40 independent
components were visually inspected and those related to eye blinks, eye movements, muscle,
and cardiac artifacts were discarded. An average of 9.51± 4.68 (SD) components were
removed. The remaining components were then projected back to the electrode space. Finally,
bad channels were reconstructed with the spherical spline interpolation method. An average of
2 channels (range 0-4) was interpolated in the sample, with no preferential spatial location
across electrodes positions. The data were then re-referenced to the average of all electrodes.
At the end of the data preprocessing, each subject with epilepsy had at least 5 minutes of
artifact-free signal, while healthy controls had at least 7 minutes of artifact-free signal. To avoid
bias in the classification due to the difference in the signal length, we used 5 minutes of artifact
free signal for each subject in the following analyses.

Cortical Source modeling

We used the individual anatomy MRI to generate individualized head models for the patients
with TLE. The anatomic MRI for source imaging consisted of a T1 isotropic three-dimensional
(3D) acquisition. For twelve patients and for the control group we used the MNI-ICBM152
default anatomy <sup>41</sup> from Brainstorm <sup>42</sup> since the 3D T1 MRI sequences were not available. The
MRI was segmented into skin, skull, and gray matter using the Computational Anatomy
Toolbox (CAT12) <sup>43</sup>. The resulting individual surfaces were then imported in Brainstorm,
where three individual surfaces adapted for Boundary Element Models (BEM) were reconstructed (inner skull, outer skull and head) and the cortical mesh was downsampled at 15,002 vertices. The co-registration of the EEG electrodes was performed using Brainstorm by projecting the EEG sensor positions on the head surface with respect to the fiducial points of the individual or the template MRI. We applied manual correction of the EEG cap on the individual anatomy whenever needed, prior to projecting the electrodes on the individual head surface. We then derived an EEG forward model using the 3-shell BEM model estimated using OpenMEEG method implemented in Brainstorm ⁴⁴. Finally, we used the weighted minimum norm imaging ⁴⁵ as the inverse model, with the Brainstorm’s default parameters setting.

## Features extraction

### Avalanche transition matrix computation

Similarly to our previous work using neuronal avalanches (Duma et al., 2023), we extracted the activity of a total of 68 regions of interest (ROIs) from the Desikan-Killiany atlas ⁴⁶. The ROIs time series were obtained by averaging the activity across the vertices composing each ROI. To study the dynamics of brain activity, we estimated “neuronal avalanches” from the source-reconstructed ROI time series. Firstly, the time series of each ROI was discretized by calculating the z-score over time and then setting positive and negative excursions beyond a threshold as 1, and the rest of the signal as 0 ⁴⁷. A neuronal avalanche begins when, in a sequence of contiguous time bins, at least one ROI is active (i.e., above threshold), and ends when all ROIs are inactive ²⁶,⁴⁸. These analyses require the time series to be binned. This is done to ensure that one is capturing critical dynamics, if present. To estimate the suitable time bin length, for each subject, each neuronal avalanche, and each time bin duration, the branching parameter $\sigma$ was estimated ⁴⁹. In fact, systems operating at criticality typically display a branching ratio ~1. The branching ratio is calculated as the geometrically averaged (over all the time bins) ratio of the number of events (activations) between the subsequent time bin (descendants) and that in the current time bin (ancestors), and then averaging over all the avalanches ⁵⁰. More specifically:

$$ \sigma_i = \prod_{j=1}^{N_{\text{bin}}-1} \left( \frac{n_{\text{events}}(j+1)}{n_{\text{events}}(j)} \right)^{\frac{1}{N_{\text{bin}}-1}} \quad (1) $$
$$ \sigma = \prod_{i=1}^{N_{aval}} (\sigma_i)^{\frac{1}{N_{aval}}} \quad (2) $$

where $\sigma_i$ is the branching parameter of the ith avalanche in the dataset, $N_{bin}$ is the total amount of bins in the i-th avalanche, $n_{events}$ (j) is the total number of events active in the j-th bin, and $N_{aval}$ is the total number of avalanche in the dataset. In our analyses, the branching ratio was 1 for bin = 1 (corresponding to bins of 4 ms). An avalanche-specific transition matrix (ATM) was calculated where element (i, j) represented the probability that region j was active at time t + $\delta$, given that region i was active at time t, where $\delta \sim 4$ ms. The ATMs were averaged element-wise across all avalanches per each participant. It is important to mention that the ATMs were computed over the broadband signal. However, given the relevance of the slow frequencies (i.e., theta and alpha bands) in the TLE 51,52, we also computed ATMs on a narrowband filtered signal (theta and alpha bands; 3-14 Hz). Moreover, we performed classification in additional consecutive frequency bands (i.e., alpha-beta, beta-low-gamma and alpha-beta-low-gamma). This was performed to evaluate the potential changes of the accuracy of classification when including only specific frequency bands. For each tested frequency band, we optimized two parameters: the threshold applied to the z-scored signals (ranging from 1.2 to 3.0) and the minimal duration of the avalanches considered (ranging from 2 to 8 time points). The chosen judgment criterion was the classification accuracy rate reached by the chosen configuration. Finally, we computed the ATMs multiple times, each time averaging over different lengths of the signals (i.e., 5, 15, 30, 60, 120, 180, and 300sec). This was performed in order to investigate the effect of signal length used in the ATMs computation on the classification performance. To mitigate the influence of the position of a given temporal window $\omega$ in the registration on the classification performance, we randomly picked 100 times $\omega$ across the registration, and we considered the median value of the classification accuracy over the permutations.

## Imaginary Coherence

In addition to the ATMs, we computed the imaginary coherence in order to compare the performance with a widely used metric of functional connectivity. The imaginary part of the coherence (ImCoh) has been demonstrated to be robust to volume conduction, improving the identification of true interactions in EEG-derived functional connectivity 53,54. We computed ImCoh as it follows:
$$ ImCoh_{xy}(f) = \frac{Im(S_{xy}(f))}{\sqrt{S_{xx}(f)S_{yy}(f)}} \quad (3) $$

where $S_{xy}$ is the cross-spectral density between region x and region y, while $S_{xx}$ and $S_{yy}$ represent the autospectral density of the two regions. The spectral densities were estimated using a multitaper method with digital prolate spheroidal sequence (DPSS) 55 windows of 10s, an overlap of 50%, and a maximum frequency resolution of 0.1Hz. To be consistent with the scheme used with the ATMs, we calculated the ImCoh over both the broadband and the narrowband filtered signal (3-14 Hz) as well as the additional frequency bands (i.e., alpha-beta, beta-low-gamma and alpha-beta-low-gamma).

## Classification pipeline

To perform the classification, we used a Support Vector Machine (SVM) classifier 56. We optimized the SVM with a 5-fold cross-validated grid search to find optimal parameters, namely the kernels (linear or Radial Basis Function-RBF) and the regularization hyperparameter C during the training step. The kernel functions enable the system to work in a high-dimensional space and to allow for more complex decision functions when the data is not linearly separable 57. The C parameter adjusts the misclassification of training examples against the simplicity of the decision. The classification scores were evaluated by cross-validation consisting of a random permutation cross-validator with a training size of 80%, a test size of 20%, and 50 re-shuffling and splitting iterations. To avoid biased results, the data was divided across patients. Samples were first shuffled and then split into a pair of training and testing sets. Such an operation has been performed 50 times (corresponding to the number of splits we fixed). To ensure the reproducibility of our results, we set the randomness with a fixed pseudo random number generator. To assess the performance of the classification we reported several metrics, namely the accuracy classification score, the Area Under the Receiver Operating Characteristic Curve (ROC AUC) from the prediction scores, the F1, the precision, the sensitivity, and the specificity scores. To investigate the interpretability of the classification performance, we studied the relative importance of the features, derived from the absolute value of the classification coefficients of the classification model. Here, we computed the median value over the splits of the cross-validation of the testing set. The analysis was conducted with Python. The code is publicly available at:

https://github.com/mccorsi/NeuronalAvalanches_TemporalLobeEpilepsy_EEG.
We performed classification training on the model with ATMs and ImCoh features (see Figure 4 for the analysis pipeline representation).

<figure><img src="image_5.png" /><figcaption>**Figure 4. Analysis pipeline.** The present figure illustrates the analytical steps used for feature extraction and machine learning based classification using a support vector machine approach.</figcaption></figure>

## Data availability statement

The data that support the findings of this study are available on request to the corresponding author. The raw data are not publicly available due to privacy or ethical restrictions. All the scripts are available at the following github page: https://github.com/mccorsi/NeuronalAvalanches_TemporalLobeEpilepsy_EEG.
## Acknowledgements

This work was supported by Ricerca Corrente 2023 funds for biomedical research of The Italian Health Ministry and European Union "NextGenerationEU", (Investimento 3.1.M4. C2), project IR0000011, EBRAINS-Italy of PNRR.

## Author contributions:

Conceptualization: MCC, PS, GMD

Methodology: MCC, ETL, PS, SC

Investigation: MCC, ETL, PS, SC, GMD

Visualization: MCC, GMD

Supervision: MC, PS, AD, PB, GMD

Writing—original draft: MCC, PS, GMD

Writing—review & editing: MCC, ETL, PS, SC, AD, PB, GMD

## Competing interests

The authors declare no competing interests.

## References

1. Juarez-Garcia A, Stokes T, Shaw B, Camosso-Stefinovic J, Baker R. The costs of epilepsy misdiagnosis in England and Wales. *Seizure*. 2006;15(8):598-605.
doi:10.1016/j.seizure.2006.08.005

2. Uldall P, Alving J, Hansen LK, Kibæk M, Buchholt J. The misdiagnosis of epilepsy in children admitted to a tertiary epilepsy centre with paroxysmal events. *Archives of Disease in Childhood*. 2006;91(3):219-221.
doi:10.1136/adc.2004.064477

3. Zaidi A, Clough P, Cooper P, Scheepers B, Fitzpatrick AP. Misdiagnosis of epilepsy: many seizure-like attacks have a cardiovascular cause. *Journal of the American College of Cardiology*. 2000;36(1):181-184.
doi:10.1016/S0735-1097(00)00700-2

4. Oto M (Meritxell). The misdiagnosis of epilepsy: Appraising risks and managing uncertainty. *Seizure*. 2017;44:143-146.
doi:10.1016/j.seizure.2016.11.029

5. Ilias L, Askounis D, Psarras J. Multimodal detection of epilepsy with deep neural networks. *Expert Systems with Applications*. 2023;213:119010.
doi:10.1016/j.eswa.2022.119010

6. Jin L, Choi JY, Bulacio J, et al. Multimodal Image Integration for Epilepsy Presurgical Evaluation: A Clinical Workflow. *Front Neurol*. 2021;12:709400.
doi:10.3389/fneur.2021.709400

7. Tatum WO, Rubboli G, Kaplan PW, et al. Clinical utility of EEG in diagnosing and monitoring epilepsy in adults. *Clinical Neurophysiology*. 2018;129(5):1056-1082.
doi:10.1016/j.clinph.2018.01.019

8. Schmidt H, Woldman W, Goodfellow M, et al. A computational biomarker of idiopathic generalized epilepsy from resting state EEG. *Epilepsia*. 2016;57(10):e200-e204.
doi:10.1111/epi.13481
9. Thangavel P, Thomas J, Sinha N, et al. Improving automated diagnosis of epilepsy from EEGs beyond IEDs. *J Neural Eng.* 2022;19(6):066017. doi:10.1088/1741-2552/ac9c93

10. Myers P, Gunnarsdottir K, Li A, et al. Diagnosing Epilepsy with Normal Interictal EEG Using Dynamic Network Models. Published online September 19, 2023:2023.08.12.23294018. doi:10.1101/2023.08.12.23294018

11. Acharya UR, Hagiwara Y, Adeli H. Automated seizure prediction. *Epilepsy & Behavior*. 2018;88:251-261. doi:10.1016/j.yebeh.2018.09.030

12. Zhang Y, Guo Y, Yang P, Chen W, Lo B. Epilepsy Seizure Prediction on EEG Using Common Spatial Pattern and Convolutional Neural Network. *IEEE J Biomed Health Inform.* 2020;24(2):465-474. doi:10.1109/JBHI.2019.2933046

13. Zhou M, Tian C, Cao R, et al. Epileptic Seizure Detection Based on EEG Signals and CNN. *Frontiers in Neuroinformatics*. 2018;12. Accessed December 4, 2023. https://www.frontiersin.org/articles/10.3389/fninf.2018.00095

14. Rijnders B, Korkmaz EE, Yildirim F. Hybrid machine learning method for a connectivity-based epilepsy diagnosis with resting-state EEG. *Med Biol Eng Comput.* 2022;60(6):1675-1689. doi:10.1007/s11517-022-02560-w

15. Varone G, Boulila W, Lo Giudice M, et al. A Machine Learning Approach Involving Functional Connectivity Features to Classify Rest-EEG Psychogenic Non-Epileptic Seizures from Healthy Controls. *Sensors*. 2021;22(1):129. doi:10.3390/s22010129

16. Wang B, Han X, Yang S, et al. An integrative prediction algorithm of drug-refractory epilepsy based on combined clinical-EEG functional connectivity features. *J Neurol.* 2022;269(3):1501-1514. doi:10.1007/s00415-021-10718-z

17. Jing J, Sun H, Kim JA, et al. Development of Expert-Level Automated Detection of Epileptiform Discharges During Electroencephalogram Interpretation. *JAMA Neurology*. 2020;77(1):103-108. doi:10.1001/jamaneurol.2019.3485

18. Thomas J, Thangavel P, Peh WY, et al. Automated Adult Epilepsy Diagnostic Tool Based on Interictal Scalp Electroencephalogram Characteristics: A Six-Center Study. *Int J Neur Syst.* 2021;31(05):2050074. doi:10.1142/S0129065720500744

19. Saeidi M, Karwowski W, Farahani FV, et al. Neural Decoding of EEG Signals with Machine Learning: A Systematic Review. *Brain Sciences*. 2021;11(11):1525. doi:10.3390/brainsci11111525

20. Liu S, Li F, Wan F. Distance to criticality undergoes critical transition before epileptic seizure attacks. *Brain Research Bulletin*. 2023;200:110684. doi:10.1016/j.brainresbull.2023.110684

21. Moosavi SA, Truccolo W. Criticality in probabilistic models of spreading dynamics in brain networks: Epileptic seizures. *PLOS Computational Biology*. 2023;19(2):e1010852. doi:10.1371/journal.pcbi.1010852

22. Corsi MC, Sorrentino P, Schwartz D, et al. Measuring neuronal avalanches to inform brain-computer interfaces. *iScience*. 2024;27(1):108734. doi:10.1016/j.isci.2023.108734
23. Sorrentino P, Troisi Lopez E, Romano A, et al. Brain fingerprint is based on the aperiodic, scale-free, neuronal activity. *NeuroImage*. 2023;277:120260.
doi:10.1016/j.neuroimage.2023.120260

24. Duma GM, Danieli A, Mento G, et al. Altered spreading of neuronal avalanches in temporal lobe epilepsy relates to cognitive performance: A resting-state fDEEG study. *Epilepsia*. 2023;64(5):1278-1288. doi:10.1111/epi.17551

25. Duma GM, Pellegrino G, Rabuffo G, et al. Altered spread of waves of activities at large scale is influenced by cortical thickness organization in temporal lobe epilepsy: a MRI - high-density EEG study. *Brain Communications*. Published online December 14, 2023;fcad348. doi:10.1093/braincomms/fcad348

26. Sorrentino P, Petkoski S, Sparaco M, et al. Whole-Brain Propagation Delays in Multiple Sclerosis, a Combined Tractography-Magnetoencephalography Study. *J Neurosci*. 2022;42(47):8807-8816. doi:10.1523/JNEUROSCI.0938-22.2022

27. Shafiezadeh S, Duma GM, Mento G, et al. Methodological Issues in Evaluating Machine Learning Models for EEG Seizure Prediction: Good Cross-Validation Accuracy Does Not Guarantee Generalization to New Patients. *Applied Sciences*. 2023;13(7):4262.
doi:10.3390/app13074262

28. Bernasconi N, Bernasconi A, Andermann F, Dubeau F, Feindel W, Reutens DC. Entorhinal cortex in temporal lobe epilepsy: A quantitative MRI study. *Neurology*. 1999;52(9):1870-1870. doi:10.1212/WNL.52.9.1870

29. Kubota Y, Enatsu R, Gonzalez-Martinez J, et al. In vivo human hippocampal cingulate connectivity: A corticocortical evoked potentials (CCEPs) study. *Clinical Neurophysiology*. 2013;124(8):1547-1556. doi:10.1016/j.clinph.2013.01.024

30. Qin L, Jiang W, Zheng J, Zhou X, Zhang Z, Liu J. Alterations Functional Connectivity in Temporal Lobe Epilepsy and Their Relationships With Cognitive Function: A Longitudinal Resting-State fMRI Study. *Frontiers in Neurology*. 2020;11. Accessed September 15, 2023. https://www.frontiersin.org/articles/10.3389/fneur.2020.00625

31. Vismer MS, Forcelli PA, Skopin MD, Gale K, Koubeissi MZ. The piriform, perirhinal, and entorhinal cortex in seizure generation. *Frontiers in Neural Circuits*. 2015;9. Accessed December 4, 2023. https://www.frontiersin.org/articles/10.3389/fncir.2015.00027

32. Rucco R, Bernardo P, Lardone A, et al. Neuronal Avalanches to Study the Coordination of Large-Scale Brain Activity: Application to Rett Syndrome. *Front Psychol*. 2020;11:550749. doi:10.3389/fpsyg.2020.550749

33. Sorrentino P, Rabuffo G, Baselice F, et al. Dynamical interactions reconfigure the gradient of cortical timescales. *Network Neuroscience*. 2023;7(1):73-85.
doi:10.1162/netn_a_00270

34. Karoly PJ, Goldenholz DM, Freestone DR, et al. Circadian and circaseptan rhythms in human epilepsy: a retrospective cohort study. *The Lancet Neurology*. 2018;17(11):977-985. doi:10.1016/S1474-4422(18)30274-6

35. Khan S, Nobili L, Khatami R, et al. Circadian rhythm and epilepsy. *The Lancet Neurology*. 2018;17(12):1098-1108. doi:10.1016/S1474-4422(18)30335-1
36. Baud MO, Kleen JK, Mirro EA, et al. Multi-day rhythms modulate seizure risk in epilepsy. *Nat Commun.* 2018;9(1):88. doi:10.1038/s41467-017-02577-y

37. Xu Y, Nguyen D, Mohamed A, et al. Frequency of a false positive diagnosis of epilepsy: A systematic review of observational studies. *Seizure - European Journal of Epilepsy*. 2016;41:167-174. doi:10.1016/j.seizure.2016.08.005

38. Reuber M, Fernández G, Bauer J, Singh DD, Elger CE. Interictal EEG Abnormalities in Patients with Psychogenic Nonepileptic Seizures. *Epilepsia*. 2002;43(9):1013-1020. doi:10.1046/j.1528-1157.2002.52301.x

39. Dorji T, Yangchen, Wangmo S, et al. Challenges in epilepsy diagnosis and management in a low-resource setting: An experience from Bhutan. *Epilepsy Research*. 2023;192:107126. doi:10.1016/j.eplepsyres.2023.107126

40. Delorme A, Makeig S. EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. *Journal of Neuroscience Methods*. 2004;134(1):9-21. doi:10.1016/j.jneumeth.2003.10.009

41. Evans AC, Janke AL, Collins DL, Baillet S. Brain templates and atlases. *NeuroImage*. 2012;62(2):911-922. doi:10.1016/j.neuroimage.2012.01.024

42. Tadel F, Baillet S, Mosher JC, Pantazis D, Leahy RM. Brainstorm: A User-Friendly Application for MEG/EEG Analysis. *Computational Intelligence and Neuroscience*. 2011;2011:1-13. doi:10.1155/2011/879716

43. Dahnke R, Yotter RA, Gaser C. Cortical thickness and central surface estimation. *NeuroImage*. 2013;65:336-348. doi:10.1016/j.neuroimage.2012.09.050

44. Gramfort A, Papadopoulo T, Olivi E, Clerc M. OpenMEEG: opensource software for quasistatic bioelectromagnetics. *BioMed Eng OnLine*. 2010;9(1):45. doi:10.1186/1475-925X-9-45

45. Hämäläinen MS, Ilmoniemi RJ. Interpreting magnetic fields of the brain: minimum norm estimates. *Med Biol Eng Comput.* 1994;32(1):35-42. doi:10.1007/BF02512476

46. Desikan RS, Ségonne F, Fischl B, et al. An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest. *NeuroImage*. 2006;31(3):968-980. doi:10.1016/j.neuroimage.2006.01.021

47. Sorrentino P, Seguin C, Rucco R, et al. The structural connectome constrains fast brain dynamics. *eLife*. 2021;10:e67400. doi:10.7554/eLife.67400

48. Shriki O, Alstott J, Carver F, et al. Neuronal Avalanches in the Resting MEG of the Human Brain. *J Neurosci.* 2013;33(16):7079-7090. doi:10.1523/JNEUROSCI.4286-12.2013

49. Haldeman C, Beggs JM. Critical Branching Captures Activity in Living Neural Networks and Maximizes the Number of Metastable States. *Phys Rev Lett.* 2005;94(5):058101. doi:10.1103/PhysRevLett.94.058101

50. Bak P, Tang C, Wiesenfeld K. Self-organized criticality: An explanation of the 1/f noise. *Phys Rev Lett.* 1987;59(4):381-384. doi:10.1103/PhysRevLett.59.381
51. Scherer M, Milosevic L, Guggenberger R, et al. Desynchronization of temporal lobe theta-band activity during effective anterior thalamus deep brain stimulation in epilepsy. *NeuroImage*. 2020;218:116967. doi:10.1016/j.neuroimage.2020.116967

52. Sip V, Scholly J, Guye M, Bartolomei F, Jirsa V. Evidence for spreading seizure as a cause of theta-alpha activity electrographic pattern in stereo-EEG seizure recordings. Taylor PN, ed. *PLoS Comput Biol*. 2021;17(2):e1008731. doi:10.1371/journal.pcbi.1008731

53. Haufe S, Ewald A. A Simulation Framework for Benchmarking EEG-Based Brain Connectivity Estimation Methodologies. *Brain Topogr*. 2019;32(4):625-642. doi:10.1007/s10548-016-0498-y

54. Nolte G, Bai O, Wheaton L, Mari Z, Vorbach S, Hallett M. Identifying true brain interaction from EEG data using the imaginary part of coherency. *Clinical Neurophysiology*. 2004;115(10):2292-2307. doi:10.1016/j.clinph.2004.04.029

55. Moore IC, Cada M. Prolate spheroidal wave functions, an introduction to the Slepian series and its properties. *Applied and Computational Harmonic Analysis*. 2004;16(3):208-230. doi:10.1016/j.acha.2004.03.004

56. Smola AJ, Schölkopf B. A tutorial on support vector regression. *Statistics and Computing*. 2004;14(3):199-222. doi:10.1023/B:STCO.0000035301.49549.88

57. Faouzi J, Colliot O. Classic Machine Learning Methods. In: Colliot O, ed. *Machine Learning for Brain Disorders*. Neuromethods. Springer US; 2023:25-75. doi:10.1007/978-1-0716-3195-9_2