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

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.

(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 20,21 .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 22,23 .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 24,25 .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 welldefined 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) 26 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 wholebrain 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.

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-betalow-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 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 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 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 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 wellknown metabolic or connectivity alterations, such as the entorhinal cortex, the superior and inferior temporal gyri, the cingulate cortex and the prefrontal dorsolateral cortex [28][29][30][31] (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).www.nature.com/scientificreports/

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 26,32 .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 24 .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 [28][29][30][31] .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 33 .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 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 [34][35][36] ).Furthermore, these results have a potential direct application in clinical practice.In fact, the rate of epilepsy misdiagnosis varies between 20 and 40% 4 and remains an issue also in specialized centers 2 .Some populations carry a particularly high risk of misdiagnosis, for example, patients with psychogenic non-epileptic seizures 37 .Misidentification of EEG abnormalities reaches over 50% in the PNES population 38 .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 9,10,17,18 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 39 .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.

EEG pre-processing
Signal preprocessing was performed via EEGLAB 14.1.2b 40.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 24,25 .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 www.nature.com/scientificreports/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 41 from Brainstorm 42 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) 43 .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 44 .Finally, we used the weighted minimum norm imaging 45 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 46 .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 47 .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 26,48 .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 49 .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 50 .More specifically: where σ 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 avalanchespecific 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 elementwise 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 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 (1) www.nature.com/scientificreports/demonstrated to be robust to volume conduction, improving the identification of true interactions in EEGderived functional connectivity 53,54 .We computed ImCoh as it follows: 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 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, betalow-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 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 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/ mccor si/ Neuro nalAv alanc hes_ Tempo ralLo beEpi lepsy_ EEG.We performed classification training on the model with ATMs and ImCoh features (see Fig. 4 for the analysis pipeline representation).

Figure 1 .
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 crossvalidation 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).

Figure 2 .
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.

Figure 3 .
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. https://doi.org/10.1038/s41598-024-64870-3www.nature.com/scientificreports/

Figure 4 .
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.

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.

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.