Investigating brain connectivity in epilepsy has attracted considerable attention (Kramer and Cash, 2012) since multiple different networks are involved in this neurological disorder (González Otárula et al., 2019; Kramer et al., 2010; van Diessen et al., 2013b). Large-scale epileptogenic networks consist of cortical and subcortical areas that are involved in seizure generation and propagation (Bartolomei et al., 2017). These networks can either represent an increase or decrease in brain synchrony (Khambhati et al., 2017; Kramer et al., 2010; Wendling et al., 2003; Yaffe et al., 2015) or elucidate regions traversed by ictal propagation which is referred to as traveling waves in studies with microscopic recordings (Martinet et al., 2017; Muller et al., 2018; Proix et al., 2018; Schevon et al., 2012; Smith et al., 2016; Weiss et al., 2013). This underlying association between brain areas in epilepsy can describe the seizure spread and termination processes (Kramer and Cash, 2012), explain seizure semiology (Chauvel and McGonigal, 2014), or assist in identifying the seizure onset zone (SOZ) (Burns et al., 2014). Traditionally, seizures have been considered to be characterized by a state of hypersynchrony. Yet recent work (Kramer and Cash, 2012) describes an early stage of desynchronization (Kramer et al., 2010; Schindler et al., 2010; Wendling et al., 2003) followed by synchronization amid seizure termination (Schindler et al., 2007; Schindler et al., 2008). Constructing connectivity maps with various recording techniques and different computational measures over a wide range of frequencies has given rise to different controversial perspectives on how seizures should be characterized (Jiruska et al., 2013). Few studies have examined the functional connectivity between SOZ and other areas of the brain during ictal periods. Electrocorticographic (ECoG) data suggested that for some patients the SOZ is isolated from the rest of the network early in the seizures allowing for SOZ detection in this way (Burns et al., 2014). It has also been suggested that ictal periods can be delineated by a steady series of states (Burns et al., 2014), although whether this is true in all patients remains controversial. Others have shown a decreased synchrony between SOZ and normal brain regions (Warren et al., 2010). It remains unclear how the degree of desynchronization is correlated with physiological parameters such as age and duration of epilepsy (van Diessen et al., 2013a).

Several studies have used multiunit recordings to explore ictal propagation networks (Martinet et al., 2017; Schevon et al., 2019; Schevon et al., 2012; Smith et al., 2016; Weiss et al., 2013). At the microscopic spatial scale, the hypersynchronous ictal core with high neural firing can be distinguished from the penumbra with relatively small and sparse firings (Schevon et al., 2012). In the early part of the seizure, the slow-moving ictal wavefront involves the core and surrounding areas. After recruitment, the low-frequency traveling waves rapidly propagate to other cortical regions (Smith et al., 2016) in a two-dimensional spatial scale (Martinet et al., 2017). However, the mechanism by which ictal discharges spread in the brain volume is still undetermined. The intense firing of neurons in the ictal core is characterized by high-frequency oscillations (HFOs) (Jefferys et al., 2012) in local field potentials (LFPs) (Weiss et al., 2013). The term HFOs has been attributed to brain activity, with multiple possible physiological and pathological neural mechanisms, between 80–500 Hz (Jacobs et al., 2012) or 30–600 Hz (Engel and da Silva, 2012). This includes high-gamma neural activities (80–200 Hz) (Ray and Maunsell, 2011) and broad-band high frequency (Arnulfo et al., 2020).

These oscillations have been analyzed for their value in SOZ localization (Höller et al., 2015; Liu et al., 2018) during ictal (Weiss et al., 2013) and interictal (Fedele et al., 2017; Gliske et al., 2018) periods. Nevertheless, manual HFO (ripple) detection is time-consuming (Gliske et al., 2018) and automatic approaches produce a large number of false positives (Bénar et al., 2010). Moreover, some interictal analyses of slow-wave sleep questioned the utility (Roehri et al., 2018) and accuracy (Jacobs et al., 2018) of HFOs to serve as a biomarker for epileptogenic tissue identification. These observations can be further supported by a resting-state SEEG study which has delineated the long-range high-frequency synchronization of physiological HFOs among the non-epileptogenic regions (Arnulfo et al., 2020). While HFOs have been employed in analyzing functional (Schindler et al., 2010) and propagation (González Otárula et al., 2019) networks, the spatiotemporal dynamics of ictal high-frequency synchronization (HFS) at macroscopic scales remain largely unknown.

In this paper, we analyzed SEEG recordings of cortical and subcortical regions. In SEEG, the EZ not only takes into account the earliest ictal EEG change, it emphasizes an anatomo-electro-clinical analysis (Kahane et al., 2006). This concept incorporates both the anatomic region that initiates the epileptic discharge as well as the ‘primary organization’ (Talairach and Bancaud, 1966) that leads to the manifestation of the clinical seizure itself (Wyllie et al., 2015). The gold standard method of confirming EZ localization is based on whether seizure freedom has been achieved by resection or ablation. The actual ground truth for the EZ location is unknown since in many cases the resection volumes may extend well beyond the EZ. In recent work, the EZ has been considered as part of a network (Jehi, 2018). These epileptogenic networks in focal epilepsy have been invoked in explaining the underlying pathogenesis of epilepsy, seizure initiation, ictal propagation, and disease progression as well as various associated comorbidities (Nair et al., 2004). This perspective is utilized in our work to analyze seizures in the context of a distributed network of interacting regions that include the EZ. Because of the importance of HFOs in epilepsy and SOZ localization, we construct synchronization networks in the 80–200 Hz range, to be in line with similar studies (Höller et al., 2015; Schindler et al., 2010). Fast rhythmic bursting neurons, which have the highest tendency to initiate seizures, are largely responsible for generating ultra-fast oscillations or ripples (80–200 Hz) (Timofeev and Steriade, 2004). We hypothesize that during seizures, the EZ has an abnormal and unique pattern of connectivity with other brain areas.

Graph analysis provides a mathematical framework for the quantification of brain connectivity (Bullmore and Sporns, 2009). Brain networks may be represented as a graph, ${\displaystyle G=\left(V,E\right)}$, in which nodes ($V$) characterize anatomical regions or electrodes, and edges ($E$) reflect structural or functional connections among them. Traditionally, brain connectivity in each time sample (layer) of dynamic networks has been evaluated independently, via single graph analysis. However, multilayer analysis allows us to model the entire data with a single super-graph, in which individual graphs for each time-sample are linked. Assuming there are T single graphs (time-samples) with N nodes in each, the super-graph has NT nodes. The multilayer structure has considerable methodological advantages over single-layer analysis (Betzel and Bassett, 2017). First, the interlayer coupling between neighboring time points in this model allows us to incorporate the continuity in neural dynamics. Second, by tuning the coupling parameter, processes with different timescales can be distinguished. Third, the extracted measures on these networks are less susceptible to noise in the data or spurious connectivity. Additionally, there exist several neurological rationales for employing a multilayer approach in seizure analysis. First, the concept of dynamic network reconfiguration has been studied in brain networks (Bassett et al., 2011; Braun et al., 2015). Previous research has shown state transitions during seizures, either through brain connectivity analysis (Burns et al., 2014) or microelectrode recordings (Smith et al., 2016). Multilayer networks have the capability to delineate these transitions and identify network reconfiguration (Mucha et al., 2010). Second, electrophysiological signals are highly non-stationary during ictal periods. As a result, traditional analysis of time-varying networks based on isolated graphs would be affected by instantaneous fluctuations rather than the underlying spatiotemporal networks. Third, seizure propagation is one of the key elements of ictal activity. Recent work has indicated the importance of multilayer modeling of complex systems when encountering spreading processes (De Domenico et al., 2016).

Consequently, we modeled spatiotemporal high-frequency connectivity using multilayer networks. We explore the question of whether the EZ can be identified by unsupervised clustering of nodes (representing SEEG contacts) in the feature space of these multilayer networks. Our connectivity-based EZ identification results show reasonable consistency with a previous approach (described as the fingerprint of the EZ) (Grinenko et al., 2018; Li et al., 2020) which uses three ictal features for EZ localization, namely: low-voltage fast activity (LFD), preictal spiking, and suppression of lower frequencies.

Understanding the connectivity dynamics of EZ and surrounding areas with the rest of the brain can shed light on underlying processes including seizure propagation and termination. We evaluate these interactions in different frequency bands. To use consistent terminology, here the term ‘high-frequency synchrony’ is mostly utilized to describe brain networks at higher frequencies (80–200 Hz). On the other hand, the word ‘connectivity’ has been assigned to both propagation networks in low frequency (3–50 Hz) and synchronization networks in high frequency (80–200 Hz). Our findings illustrate early high-frequency desynchronization and late increase in brain connectivity during seizures. Although, seizures usually initiate with a loss of synchronization in a small area, their termination process demands widespread brain connectivity across multiple scales. During seizure cessation, the brain experiences a critical transition with a hysteresis behavior (Kramer et al., 2012), indicating the future state depends on the current one. Accordingly, we hypothesize that post-termination high-frequency synchrony is correlated with pre-termination connectivity in low frequency.

Our results investigate large-scale ictal brain connectivity as it relates to EZ localization and seizure generation, propagation, and termination. We showed that the EZ exhibits differential dynamics compared to other brain regions during seizures. This observation is independent of seizure type, etiology, presence or absence of MRI lesion, and occurred in both temporal and extra-temporal lobe epilepsy. High-frequency oscillations have been widely studied as a potential biomarker in epilepsy (Zijlmans et al., 2012), including both early- and delayed-onset ictal HFOs (Weiss et al., 2013). However, there is still some debate on their spatial specificity, timing, and appropriate detection algorithms. As an alternative marker, we computed high-frequency synchrony (HFS) over the entire seizure period. HFS was first modeled using a multilayer network and the potential EZ was identified using a novel measure of centrality. Our approach to identifying the EZ is complementary to other approaches avoiding explicit assumptions regarding seizure patterns and dynamics.

The multilayer network approach used here for EZ prediction and identification of state transitions allows us to explore seizure dynamics in a graph-theoretic context. Traditionally in brain network studies, researchers explore graph measures like clustering and centrality, in single-layer (non-dynamic) networks (Ridley et al., 2015; van Diessen et al., 2014). However, the multilayer framework with an adjustable coupling parameter can reveal processes with different timescales and facilitate defining of new measures (Betzel and Bassett, 2017). By adjusting the coupling between layers, we can overcome the current challenges of selecting an appropriate time interval between different samples of a temporal network. Here, we assumed a relatively short time interval between layers (500 ms) to have a high temporal resolution. By increasing the coupling parameter between adjacent layers, we can investigate processes with slower timescales.

We used an unsupervised clustering algorithm to find a set of nodes among all channels, i.e., the target cluster, as the predicted EZ. The presented method was based on hierarchical clustering and a cost function that combines the separation of clusters and the compactness of the target group. The optimized coupling parameter varies substantially between patients (Figure 2—figure supplement 1), verifying the necessity of analyzing multiple timescales. In some cases, this automatic unsupervised clustering may fail to find all nodes (contacts) with features that would identify them as belonging to the EZ using other methods, such as the fingerprint method that finds a distinctive combination of pre-ictal spiking, low-frequency suppression, and low-voltage fast activity (LFD). This is reflected in the differences in contacts identified as EZ between the two approaches, which is reported in Figure 2—source data 1. However, the complementary dynamic-synchrony-based approach described here also finds plausible EZ contacts in cases where the fingerprint does not, as mentioned in the result section. Finally, based on the predicted EZ and information regarding the resection areas, the SEEG electrodes were categorized into three groups; EZ, Resected non-EZ (RnEZ), and non-resected (nR). We employed this categorization to compute the dynamics of regional ictal connectivity in low and high-frequency bands.

There is substantial evidence of multiscale alterations in structural and functional brain networks in epilepsy (Bernhardt et al., 2016; Englot et al., 2016; Tavakol et al., 2019). A meta-analysis of a dozen interictal studies with variant imaging methodologies showed increased clustering and path length (van Diessen et al., 2014). Resting-state fMRI research has revealed decreased inter-hemispheric functional connectivity in medial and lateral temporal regions among patients with TLE (Maccotta et al., 2013; Pittau et al., 2012). Ictal SEEG recordings of focal cortical dysplasia (FCD) type II have shown that nodes within the lesion have higher values of out-degree and betweenness centrality in the gamma range (30–80 Hz)(Varotto et al., 2012). Using cortico-cortical evoked potentials (CCEPs), we recently presented the differences in effective connectivity between FCD types I and II (Shahabi et al., 2021). These observations have created a new field of research known as connectivity-based biomarkers in epilepsy (Larivière et al., 2021), which is mostly performed by automatic approaches and machine learning algorithms. Applications include estimating neurocognitive performance (Mazrooyisebdani et al., 2020), lesion detection in FCD type II (Jin et al., 2018), lateralization of seizure focus (Barron et al., 2015), and ictal onset zone identification (Van Eyndhoven et al., 2019). Although these measures have been extracted using non-invasive resting-state imaging techniques, in many cases they are limited to specific types of epilepsy, such as FCD type II or lesional patients. On the other hand, we demonstrated that during seizures nodes belonging to the EZ share a unique pattern of centrality, which was verified in the majority of our patients with different etiologies. Our results are consistent with iEEG (Warren et al., 2010) and ECoG (Burns et al., 2014) studies that demonstrated the isolation of SOZ from the rest of the network. More importantly, we investigated the relation between EZ-nR connectivity and patient demographics. We observed a negative correlation between these parameters, in which the EZ and non-resected regions become more desynchronized by aging and the duration of epilepsy (Figure 5). This finding can be explained by abnormal neuroplasticity, where the continuous recruitment of epileptogenic networks intensifies their anomaly (Tavakol et al., 2019). Related work has shown a negative correlation between TLE duration and overall interictal PLI at low frequencies (0.5–48 Hz) (van Dellen et al., 2009). A resting-state fMRI study indicated that the connectivity between two contralateral regions in the mesial temporal lobe is negatively correlated with the duration of epilepsy (Liao et al., 2010). Our results along with the studies cited suggest that network abnormalities can portray the disease severity (Tavakol et al., 2019). To further verify these findings, additional studies should be performed on a larger dataset while considering other influencing factors, including sex, handedness, and lateralization of the SOZ.

The origin and frequency range of HFOs in ictal and interictal states are still disputable (González Otárula et al., 2019; Jacobs et al., 2012; Korzeniewska et al., 2014; Tamilia et al., 2018; Weiss et al., 2013). Recent research has emphasized that narrow-band physiological HFOs, not pathological HFOs or broad-band multi-unit activity (MUA), are responsible for long-range high-frequency synchronization in interictal recordings (Arnulfo et al., 2020). Also, it has been suggested that widespread HFO synchronization should be a characteristic of healthy brain activity (Arnulfo et al., 2020). In our findings, the ictal HFS between the EZ and non-resected areas was the most distinctive pattern among others and showed a significant desynchronization in the early and middle parts of the seizure (Figure 4). In this period, the exclusive presence of pathological HFOs inside the EZ resulted in decreased EZ-nR synchrony while the nR-nR connectivity was left intact. Towards seizure termination, this unique activity faded (Smith et al., 2016) and elevated the similarity of signals in the EZ and non-resected areas, which is consistent with a previous study in the same frequency range (Schindler et al., 2010). Once the seizure terminated and the brain resumed a healthy activity, we observed a widespread HFS irrespective of pathology, a possible outcome of physiological HFOs.

We used the right singular vectors in mlEVC decomposition to explore the network topology during seizures (Figure 3). This feature demonstrated the reconfiguration of brain networks and was later used for clustering these dynamics into brain states. The profound alterations of these features revealed the extensive changes in network topology, even when the nR-nR connectivity was relatively stable. Several studies have also displayed the presence of different brain states in the ictal period (Burns et al., 2014; Khambhati et al., 2015). However, our findings are distinctive in that they present dissimilar brain states among different seizures of the same individual, suggesting variant structures for seizure generation and propagation. It has been shown that TLE seizures can be divided into several categories and in an overwhelming majority of patients they belong to more than one category (Bartolomei et al., 2010). Consequently, in each individual, seizures can be generated by various and separate epileptogenic structures (Bartolomei et al., 2017; Bartolomei et al., 2010). Another study of TLE has identified patients with different seizure semiologies (Loesch et al., 2015), a possible outcome of variations in propagation networks. Recently, the Fingerprint features of the EZ were employed for clustering the seizures of each patient. In the presence of these variations in epileptogenic networks, it is essential to record multiple seizures (Struck et al., 2015) and prolonged interictal (Gliske et al., 2018) data. Also, seizure dissimilarities have been assessed by brain network measures in a wide frequency range (Schroeder et al., 2020). Here, we presented a new approach to group seizures based on the features of high-frequency synchronization networks. Additional work can be done to associate these findings with the patient’s semiology.

While we mainly focused on the dynamics of high-frequency synchrony, the low-frequency (3–50 Hz) propagation networks were also computed, to have a better understanding of large-scale brain interactions during seizures. Here, the PLI was used to address delays between ictal discharges in different channels. In early-seizure, the EZ-nR PLI was significantly smaller than connectivity among non-resected areas. This finding can be related to early-onset suppression of low frequencies in the EZ, as demonstrated in recent studies (Aupy et al., 2019; Grinenko et al., 2018). The low-frequency suppression can be described by rapid inhibitory synaptic currents which resist the seizure spread, causing inhibitory restraint (Schevon et al., 2012; Trevelyan, 2009). By mid-seizure, the ictal wavefront has recruited more areas and ictal discharges are traversing farther regions (Smith et al., 2016). Our results indicated an extensive increase in brain connectivity approaching seizure termination (Figure 6). This observation is consistent with multiscale recordings, reporting an increase in spatial and temporal correlation before seizures cease (Kramer et al., 2012; Martinet et al., 2017), as well as macroscopic cortical networks (Jiruska et al., 2013; Kramer et al., 2010; Schindler et al., 2007). The underlying mechanisms of seizure generation, evolution, and termination remain a matter of debate (Martinet et al., 2017). In addition to studies based on microscopic recordings and modeling, macro-electrode studies and a large-scale explanation of brain dynamics are necessary to develop a full understanding of this issue. Our findings suggest that at seizure-onset, the EZ loses synchronization with the rest of the network and the brain enters a desynchronized state. This condition continues through mid-seizure until the ictal wavefront has recruited many areas, including the core and penumbra (Smith et al., 2016). Subsequently, brain connectivity increases temporally and spatially, indicating an upcoming critical transition (Kramer et al., 2012). It has been shown theoretically and experimentally that brain networks require a higher coupling strength to restore synchronization than to lose synchronization (Kim et al., 2018). Our results support this idea since considerably stronger connectivity was needed to terminate seizures in comparison to that in the pre-ictal period.

It is important to remember that the recording techniques, signal processing approaches, and frequencies of interest are crucial in any study of brain connectivity. These parameters can be the reasons for current controversies (Jiruska et al., 2013). Recent research has investigated the effect of incomplete electrode sampling on network metrics (Conrad et al., 2020), in which EVC showed the highest reliability in comparison to other network metrics when randomly removing network nodes. However, they also showed that confidence intervals in network measures can vary significantly between patients. For this reason, care should be taken when interpreting and comparing different studies.

Employing the bipolar montage and a robust measure of synchrony, helped us to eliminate the effect of volume conduction and muscle artifacts (see Methods). In contrast, classic measures like coherence strikingly increase spurious connectivity. Additionally, our findings do not simply reflect the distance among the electrodes. We normalized the synchronization matrices with respect to the pre-ictal period based on an element-wise approach, previously suggested in Burns et al., 2014. This process reduces the chance of distance alone affecting the value of connectivity.

The results and codes generated during the current study are available in the following repository, https://data.mendeley.com/datasets/t8bvh5m8bp/1.

1. 1. Shahabi H
   2. Leahy R
   3. Nair D

   (2022) Mendeley Data

   Brain Connectivity During Seizures.

   https://doi.org/10.17632/t8bvh5m8bp

1. 1. Arnulfo G
   2. Wang SH
   3. Myrov V
   4. Toselli B
   5. Hirvonen J
   6. Fato MM
   7. Nobili L
   8. Cardinale F
   9. Rubino A
   10. Zhigalov A
   11. Palva S
   12. Palva JM

   (2020) Long-range phase synchronization of high-frequency oscillations in human cortex

   Nature Communications 11:5363.

   https://doi.org/10.1038/s41467-020-18975-8

   • PubMed
   • Google Scholar
2. 1. Aupy J
   2. Wendling F
   3. Taylor K
   4. Bulacio J
   5. Gonzalez-Martinez J
   6. Chauvel P

   (2019) Cortico-Striatal synchronization in human focal seizures

   Brain 142:1282–1295.

   https://doi.org/10.1093/brain/awz062

   • PubMed
   • Google Scholar
3. 1. Barron DS
   2. Fox PT
   3. Pardoe H
   4. Lancaster J
   5. Price LR
   6. Blackmon K
   7. Berry K
   8. Cavazos JE
   9. Kuzniecky R
   10. Devinsky O
   11. Thesen T

   (2015) Thalamic functional connectivity predicts seizure laterality in individual TLE patients: Application of a biomarker development strategy

   NeuroImage. Clinical 7:273–280.

   https://doi.org/10.1016/j.nicl.2014.08.002

   • PubMed
   • Google Scholar
4. 1. Bartolomei F
   2. Cosandier-Rimele D
   3. McGonigal A
   4. Aubert S
   5. Régis J
   6. Gavaret M
   7. Wendling F
   8. Chauvel P

   (2010) From mesial temporal lobe to temporoperisylvian seizures: A quantified study of temporal lobe seizure networks

   Epilepsia 51:2147–2158.

   https://doi.org/10.1111/j.1528-1167.2010.02690.x

   • PubMed
   • Google Scholar
5. 1. Bartolomei F
   2. Lagarde S
   3. Wendling F
   4. McGonigal A
   5. Jirsa V
   6. Guye M
   7. Bénar C

   (2017) Defining epileptogenic networks: contribution of SEEG and signal analysis

   Epilepsia 58:1131–1147.

   https://doi.org/10.1111/epi.13791

   • PubMed
   • Google Scholar
6. 1. Bassett DS
   2. Wymbs NF
   3. Porter MA
   4. Mucha PJ
   5. Carlson JM
   6. Grafton ST

   (2011) Dynamic reconfiguration of human brain networks during learning

   PNAS 108:7641–7646.

   https://doi.org/10.1073/pnas.1018985108

   • PubMed
   • Google Scholar
7. 1. Bénar CG
   2. Chauvière L
   3. Bartolomei F
   4. Wendling F

   (2010) Pitfalls of high-pass filtering for detecting epileptic oscillations: A technical note on “false” ripples

   Clinical Neurophysiology 121:301–310.

   https://doi.org/10.1016/j.clinph.2009.10.019

   • PubMed
   • Google Scholar
8. 1. Bernhardt BC
   2. Bernasconi A
   3. Liu M
   4. Hong SJ
   5. Caldairou B
   6. Goubran M
   7. Guiot MC
   8. Hall J
   9. Bernasconi N

   (2016) The spectrum of structural and functional imaging abnormalities in temporal lobe epilepsy

   Annals of Neurology 80:142–153.

   https://doi.org/10.1002/ana.24691

   • PubMed
   • Google Scholar
9. 1. Betzel RF
   2. Bassett DS

   (2017) Multi-scale brain networks

   NeuroImage 160:73–83.

   https://doi.org/10.1016/j.neuroimage.2016.11.006

   • PubMed
   • Google Scholar
10. 1. Bonacich P

    (1972) Factoring and weighting approaches to status scores and clique identification

    The Journal of Mathematical Sociology 2:113–120.

    https://doi.org/10.1080/0022250X.1972.9989806

    • Google Scholar
11. 1. Braun U
    2. Schäfer A
    3. Walter H
    4. Erk S
    5. Romanczuk-Seiferth N
    6. Haddad L
    7. Schweiger JI
    8. Grimm O
    9. Heinz A
    10. Tost H
    11. Meyer-Lindenberg A
    12. Bassett DS

    (2015) Dynamic reconfiguration of frontal brain networks during executive cognition in humans

    PNAS 112:11678–11683.

    https://doi.org/10.1073/pnas.1422487112

    • PubMed
    • Google Scholar
12. 1. Bullmore E
    2. Sporns O

    (2009) Complex brain networks: graph theoretical analysis of structural and functional systems

    Nature Reviews Neuroscience 10:186–198.

    https://doi.org/10.1038/nrn2575

    • PubMed
    • Google Scholar
13. 1. Burns SP
    2. Santaniello S
    3. Yaffe RB
    4. Jouny CC
    5. Crone NE
    6. Bergey GK
    7. Anderson WS
    8. Sarma SV

    (2014) Network dynamics of the brain and influence of the epileptic seizure onset zone

    PNAS 111:E5321–E5330.

    https://doi.org/10.1073/pnas.1401752111

    • PubMed
    • Google Scholar
14. 1. Chauvel P
    2. McGonigal A

    (2014) Emergence of semiology in epileptic seizures

    Epilepsy & Behavior 38:94–103.

    https://doi.org/10.1016/j.yebeh.2013.12.003

    • PubMed
    • Google Scholar
15. 1. Conrad EC
    2. Bernabei JM
    3. Kini LG
    4. Shah P
    5. Mikhail F
    6. Kheder A
    7. Shinohara RT
    8. Davis KA
    9. Bassett DS
    10. Litt B

    (2020) The sensitivity of network statistics to incomplete electrode sampling on intracranial eeg

    Network Neuroscience 4:484–506.

    https://doi.org/10.1162/netn_a_00131

    • PubMed
    • Google Scholar
16. 1. De Domenico M
    2. Granell C
    3. Porter MA
    4. Arenas A

    (2016) The physics of spreading processes in multilayer networks

    Nature Physics 12:901–906.

    https://doi.org/10.1038/nphys3865

    • Google Scholar
17. 1. Engel J
    2. da Silva FL

    (2012) High-frequency oscillations - where we are and where we need to go

    Progress in Neurobiology 98:316–318.

    https://doi.org/10.1016/j.pneurobio.2012.02.001

    • PubMed
    • Google Scholar
18. 1. Englot DJ
    2. Konrad PE
    3. Morgan VL

    (2016) Regional and global connectivity disturbances in focal epilepsy, related neurocognitive sequelae, and potential mechanistic underpinnings

    Epilepsia 57:1546–1557.

    https://doi.org/10.1111/epi.13510

    • PubMed
    • Google Scholar
19. 1. Fedele T
    2. Burnos S
    3. Boran E
    4. Krayenbühl N
    5. Hilfiker P
    6. Grunwald T
    7. Sarnthein J

    (2017) Resection of high frequency oscillations predicts seizure outcome in the individual patient

    Scientific Reports 7:13836.

    https://doi.org/10.1038/s41598-017-13064-1

    • PubMed
    • Google Scholar
20. 1. Gliske SV
    2. Irwin ZT
    3. Chestek C
    4. Hegeman GL
    5. Brinkmann B
    6. Sagher O
    7. Garton HJL
    8. Worrell GA
    9. Stacey WC

    (2018) Variability in the location of high frequency oscillations during prolonged intracranial EEG recordings

    Nature Communications 9:2155.

    https://doi.org/10.1038/s41467-018-04549-2

    • PubMed
    • Google Scholar
21. 1. González Otárula KA
    2. von Ellenrieder N
    3. Cuello-Oderiz C
    4. Dubeau F
    5. Gotman J

    (2019) High-Frequency oscillation networks and surgical outcome in adult focal epilepsy

    Annals of Neurology 85:485–494.

    https://doi.org/10.1002/ana.25442

    • PubMed
    • Google Scholar
22. 1. Gonzalez-Martinez J
    2. Mullin J
    3. Vadera S
    4. Bulacio J
    5. Hughes G
    6. Jones S
    7. Enatsu R
    8. Najm I

    (2014) Stereotactic placement of depth electrodes in medically intractable epilepsy

    Journal of Neurosurgery 120:639–644.

    https://doi.org/10.3171/2013.11.JNS13635

    • PubMed
    • Google Scholar
23. 1. Grinenko O
    2. Li J
    3. Mosher JC
    4. Wang IZ
    5. Bulacio JC
    6. Gonzalez-Martinez J
    7. Nair D
    8. Najm I
    9. Leahy RM
    10. Chauvel P

    (2018) A fingerprint of the epileptogenic zone in human epilepsies

    Brain 141:117–131.

    https://doi.org/10.1093/brain/awx306

    • PubMed
    • Google Scholar
24. 1. Gutkin BS
    2. Laing CR
    3. Colby C
    4. Chow CC
    5. Ermentrout GB

    (2001) Turning on and off with excitation: The role of spike-timing

    Journal of Computational Neuroscience 11:121–134.

    https://doi.org/10.1023/A:1012837415096

    • Google Scholar
25. 1. Halkidi M
    2. Batistakis Y
    3. Vazirgiannis M

    (2001) On clustering validation techniques

    Journal of Intelligent Information Systems 17:107–145.

    https://doi.org/10.1023/A:1012801612483

    • Google Scholar
26. 1. Höller Y
    2. Kutil R
    3. Klaffenböck L
    4. Thomschewski A
    5. Höller PM
    6. Bathke AC
    7. Jacobs J
    8. Taylor AC
    9. Nardone R
    10. Trinka E

    (2015) High-frequency oscillations in epilepsy and surgical outcome: a meta-analysis

    Frontiers in Human Neuroscience 9:574.

    https://doi.org/10.3389/fnhum.2015.00574

    • PubMed
    • Google Scholar
27. 1. Jacobs J
    2. Staba R
    3. Asano E
    4. Otsubo H
    5. Wu JY
    6. Zijlmans M
    7. Mohamed I
    8. Kahane P
    9. Dubeau F
    10. Navarro V
    11. Gotman J

    (2012) High-Frequency oscillations (hfos) in clinical epilepsy

    Progress in Neurobiology 98:302–315.

    https://doi.org/10.1016/j.pneurobio.2012.03.001

    • PubMed
    • Google Scholar
28. 1. Jacobs J
    2. Wu JY
    3. Perucca P
    4. Zelmann R
    5. Mader M
    6. Dubeau F
    7. Mathern GW
    8. Schulze-Bonhage A
    9. Gotman J

    (2018) Removing high-frequency oscillations: A prospective multicenter study on seizure outcome

    Neurology 91:e1040–e1052.

    https://doi.org/10.1212/WNL.0000000000006158

    • PubMed
    • Google Scholar
29. 1. Jefferys JGR
    2. Menendez de la Prida L
    3. Wendling F
    4. Bragin A
    5. Avoli M
    6. Timofeev I
    7. Lopes da Silva FH

    (2012) Mechanisms of physiological and epileptic HFO generation

    Progress in Neurobiology 98:250–264.

    https://doi.org/10.1016/j.pneurobio.2012.02.005

    • PubMed
    • Google Scholar
30. 1. Jehi L

    (2018) The epileptogenic zone: Concept and definition

    Epilepsy Currents 18:12–16.

    https://doi.org/10.5698/1535-7597.18.1.12

    • PubMed
    • Google Scholar
31. 1. Jin B
    2. Krishnan B
    3. Adler S
    4. Wagstyl K
    5. Hu W
    6. Jones S
    7. Najm I
    8. Alexopoulos A
    9. Zhang K
    10. Zhang J
    11. Ding M
    12. Wang S
    13. Pediatric Imaging, Neurocognition, and Genetics Study
    14. Wang ZI

    (2018) Automated detection of focal cortical dysplasia type II with surface-based magnetic resonance imaging postprocessing and machine learning

    Epilepsia 59:982–992.

    https://doi.org/10.1111/epi.14064

    • PubMed
    • Google Scholar
32. 1. Jiruska P
    2. de Curtis M
    3. Jefferys JGR
    4. Schevon CA
    5. Schiff SJ
    6. Schindler K

    (2013) Synchronization and desynchronization in epilepsy: Controversies and hypotheses

    The Journal of Physiology 591:787–797.

    https://doi.org/10.1113/jphysiol.2012.239590

    • PubMed
    • Google Scholar
33. 1. Kahane P
    2. Landré E
    3. Minotti L
    4. Francione S
    5. Ryvlin P

    (2006)

    The Bancaud and Talairach view on the epileptogenic zone: a working hypothesis

    Epileptic Disorders 8:S16–S26.

    • PubMed
    • Google Scholar
34. 1. Khambhati AN
    2. Davis KA
    3. Oommen BS
    4. Chen SH
    5. Lucas TH
    6. Litt B
    7. Bassett DS

    (2015) Dynamic network drivers of seizure generation, propagation and termination in human neocortical epilepsy

    PLOS Computational Biology 11:e1004608.

    https://doi.org/10.1371/journal.pcbi.1004608

    • PubMed
    • Google Scholar
35. 1. Khambhati AN
    2. Bassett DS
    3. Oommen BS
    4. Chen SH
    5. Lucas TH
    6. Davis KA
    7. Litt B

    (2017) Recurring functional interactions predict network architecture of interictal and ictal states in neocortical epilepsy

    ENeuro 4:ENEURO.0091-16.2017.

    https://doi.org/10.1523/ENEURO.0091-16.2017

    • PubMed
    • Google Scholar
36. 1. Kim H
    2. Moon JY
    3. Mashour GA
    4. Lee UC

    (2018) Mechanisms of hysteresis in human brain networks during transitions of consciousness and unconsciousness: Theoretical principles and empirical evidence

    PLOS Computational Biology 14:e1006424.

    https://doi.org/10.1371/journal.pcbi.1006424

    • PubMed
    • Google Scholar
37. 1. Kiselev VY
    2. Kirschner K
    3. Schaub MT
    4. Andrews T
    5. Yiu A
    6. Chandra T
    7. Natarajan KN
    8. Reik W
    9. Barahona M
    10. Green AR
    11. Hemberg M

    (2017) Sc3: consensus clustering of single-cell RNA-seq data

    Nature Methods 14:483–486.

    https://doi.org/10.1038/nmeth.4236

    • PubMed
    • Google Scholar
38. 1. Kivela M
    2. Arenas A
    3. Barthelemy M
    4. Gleeson JP
    5. Moreno Y
    6. Porter MA

    (2014) Multilayer networks

    Journal of Complex Networks 2:203–271.

    https://doi.org/10.1093/comnet/cnu016

    • Google Scholar
39. 1. Korzeniewska A
    2. Cervenka MC
    3. Jouny CC
    4. Perilla JR
    5. Harezlak J
    6. Bergey GK
    7. Franaszczuk PJ
    8. Crone NE

    (2014) Ictal propagation of high frequency activity is recapitulated in interictal recordings: Effective connectivity of epileptogenic networks recorded with intracranial EEG

    NeuroImage 101:96–113.

    https://doi.org/10.1016/j.neuroimage.2014.06.078

    • PubMed
    • Google Scholar
40. 1. Kramer M.A
    2. Eden UT
    3. Kolaczyk ED
    4. Zepeda R
    5. Eskandar EN
    6. Cash SS

    (2010) Coalescence and fragmentation of cortical networks during focal seizures

    The Journal of Neuroscience 30:10076–10085.

    https://doi.org/10.1523/JNEUROSCI.6309-09.2010

    • PubMed
    • Google Scholar
41. 1. Kramer M.A
    2. Cash SS

    (2012) Epilepsy as a disorder of cortical network organization

    The Neuroscientist 18:360–372.

    https://doi.org/10.1177/1073858411422754

    • PubMed
    • Google Scholar
42. 1. Kramer MA
    2. Truccolo W
    3. Eden UT
    4. Lepage KQ
    5. Hochberg LR
    6. Eskandar EN
    7. Madsen JR
    8. Lee JW
    9. Maheshwari A
    10. Halgren E
    11. Chu CJ
    12. Cash SS

    (2012) Human seizures self-terminate across spatial scales via a critical transition

    PNAS 109:21116–21121.

    https://doi.org/10.1073/pnas.1210047110

    • PubMed
    • Google Scholar
43. 1. Larivière S
    2. Bernasconi A
    3. Bernasconi N
    4. Bernhardt BC

    (2021) Connectome biomarkers of drug-resistant epilepsy

    Epilepsia 62:6–24.

    https://doi.org/10.1111/epi.16753

    • PubMed
    • Google Scholar
44. Conference

    1. Li T
    2. Ding C

    (2008) Weighted Consensus Clustering

    Proceedings of the 2008 SIAM International Conference on Data Mining. pp. 798–809.

    https://doi.org/10.1137/1.9781611972788.72

    • Google Scholar
45. 1. Li J
    2. Grinenko O
    3. Mosher JC
    4. Gonzalez-Martinez J
    5. Leahy RM
    6. Chauvel P

    (2020) Learning to define an electrical biomarker of the epileptogenic zone

    Human Brain Mapping 41:429–441.

    https://doi.org/10.1002/hbm.24813

    • PubMed
    • Google Scholar
46. 1. Liao W
    2. Zhang Z
    3. Pan Z
    4. Mantini D
    5. Ding J
    6. Duan X
    7. Luo C
    8. Lu G
    9. Chen H

    (2010) Altered functional connectivity and small-world in mesial temporal lobe epilepsy

    PLOS ONE 5:e8525.

    https://doi.org/10.1371/journal.pone.0008525

    • PubMed
    • Google Scholar
47. 1. Liu S
    2. Gurses C
    3. Sha Z
    4. Quach MM
    5. Sencer A
    6. Bebek N
    7. Curry DJ
    8. Prabhu S
    9. Tummala S
    10. Henry TR
    11. Ince NF

    (2018) Stereotyped high-frequency oscillations discriminate seizure onset zones and critical functional cortex in focal epilepsy

    Brain 141:713–730.

    https://doi.org/10.1093/brain/awx374

    • PubMed
    • Google Scholar
48. 1. Loesch AM
    2. Feddersen B
    3. Tezer FI
    4. Hartl E
    5. Rémi J
    6. Vollmar C
    7. Noachtar S

    (2015) Seizure semiology identifies patients with bilateral temporal lobe epilepsy

    Epilepsy Research 109:197–202.

    https://doi.org/10.1016/j.eplepsyres.2014.11.002

    • PubMed
    • Google Scholar
49. 1. Maccotta L
    2. He BJ
    3. Snyder AZ
    4. Eisenman LN
    5. Benzinger TL
    6. Ances BM
    7. Corbetta M
    8. Hogan RE

    (2013) Impaired and facilitated functional networks in temporal lobe epilepsy

    NeuroImage. Clinical 2:862–872.

    https://doi.org/10.1016/j.nicl.2013.06.011

    • PubMed
    • Google Scholar
50. 1. Martinet LE
    2. Fiddyment G
    3. Madsen JR
    4. Eskandar EN
    5. Truccolo W
    6. Eden UT
    7. Cash SS
    8. Kramer MA

    (2017) Human seizures couple across spatial scales through travelling wave dynamics

    Nature Communications 8:14896.

    https://doi.org/10.1038/ncomms14896

    • PubMed
    • Google Scholar
51. 1. Mazrooyisebdani M
    2. Nair VA
    3. Garcia-Ramos C
    4. Mohanty R
    5. Meyerand E
    6. Hermann B
    7. Prabhakaran V
    8. Ahmed R

    (2020) Graph theory analysis of functional connectivity combined with machine learning approaches demonstrates widespread network differences and predicts clinical variables in temporal lobe epilepsy

    Brain Connectivity 10:39–50.

    https://doi.org/10.1089/brain.2019.0702

    • PubMed
    • Google Scholar
52. 1. Modur PN
    2. Zhang S
    3. Vitaz TW

    (2011) Ictal high-frequency oscillations in neocortical epilepsy: Implications for seizure localization and surgical resection

    Epilepsia 52:1792–1801.

    https://doi.org/10.1111/j.1528-1167.2011.03165.x

    • PubMed
    • Google Scholar
53. 1. Mucha PJ
    2. Richardson T
    3. Macon K
    4. Porter MA
    5. Onnela JP

    (2010) Community structure in time-dependent, multiscale, and multiplex networks

    Science 328:876–878.

    https://doi.org/10.1126/science.1184819

    • PubMed
    • Google Scholar
54. 1. Muller L
    2. Chavane F
    3. Reynolds J
    4. Sejnowski TJ

    (2018) Cortical travelling waves: Mechanisms and computational principles

    Nature Reviews. Neuroscience 19:255–268.

    https://doi.org/10.1038/nrn.2018.20

    • PubMed
    • Google Scholar
55. 1. Nair DR
    2. Mohamed A
    3. Burgess R
    4. Lüders H

    (2004)

    A critical review of the different conceptual hypotheses framing human focal epilepsy

    Epileptic Disorders 6:77–83.

    • PubMed
    • Google Scholar
56. 1. Nissen IA
    2. van Klink NEC
    3. Zijlmans M
    4. Stam CJ
    5. Hillebrand A

    (2016) Brain areas with epileptic high frequency oscillations are functionally isolated in MEG virtual electrode networks

    Clinical Neurophysiology 127:2581–2591.

    https://doi.org/10.1016/j.clinph.2016.04.013

    • PubMed
    • Google Scholar
57. Preprint

    1. Pascual-Marqui RD

    (2007) Coherence and phase synchronization: Generalization to pairs of multivariate time series, and removal of zero-lag contributions

    arXiv.

    https://doi.org/10.48550/arXiv.0706.1776

    • Google Scholar
58. 1. Pittau F
    2. Grova C
    3. Moeller F
    4. Dubeau F
    5. Gotman J

    (2012) Patterns of altered functional connectivity in mesial temporal lobe epilepsy

    Epilepsia 53:1013–1023.

    https://doi.org/10.1111/j.1528-1167.2012.03464.x

    • PubMed
    • Google Scholar
59. 1. Prichard D
    2. Theiler J

    (1994) Generating surrogate data for time series with several simultaneously measured variables

    Physical Review Letters 73:951–954.

    https://doi.org/10.1103/PhysRevLett.73.951

    • PubMed
    • Google Scholar
60. 1. Proix T
    2. Jirsa VK
    3. Bartolomei F
    4. Guye M
    5. Truccolo W

    (2018) Predicting the spatiotemporal diversity of seizure propagation and termination in human focal epilepsy

    Nature Communications 9:1088.

    https://doi.org/10.1038/s41467-018-02973-y

    • PubMed
    • Google Scholar
61. 1. Ray S
    2. Maunsell JHR

    (2011) Different origins of gamma rhythm and high-gamma activity in macaque visual cortex

    PLOS Biology 9:e1000610.

    https://doi.org/10.1371/journal.pbio.1000610

    • PubMed
    • Google Scholar
62. 1. Ridley BGY
    2. Rousseau C
    3. Wirsich J
    4. Le Troter A
    5. Soulier E
    6. Confort-Gouny S
    7. Bartolomei F
    8. Ranjeva JP
    9. Achard S
    10. Guye M

    (2015) Nodal approach reveals differential impact of lateralized focal epilepsies on hub reorganization

    NeuroImage 118:39–48.

    https://doi.org/10.1016/j.neuroimage.2015.05.096

    • PubMed
    • Google Scholar
63. 1. Roehri N
    2. Pizzo F
    3. Lagarde S
    4. Lambert I
    5. Nica A
    6. McGonigal A
    7. Giusiano B
    8. Bartolomei F
    9. Bénar CG

    (2018) High-frequency oscillations are not better biomarkers of epileptogenic tissues than spikes

    Annals of Neurology 83:84–97.

    https://doi.org/10.1002/ana.25124

    • PubMed
    • Google Scholar
64. Software

    1. Rstudio Team

    (2018) RStudio: integrated development environment for R

    RStudio.

    http://www.rstudio.com/
65. 1. Rubinov M
    2. Sporns O.

    (2010) Complex network measures of brain connectivity: Uses and interpretations

    Neuroimage 52:1059–1069.

    https://doi.org/10.1016/j.neuroimage.2009.10.003

    • Google Scholar
66. 1. Schevon CA
    2. Weiss SA
    3. McKhann G
    4. Goodman RR
    5. Yuste R
    6. Emerson RG
    7. Trevelyan AJ

    (2012) Evidence of an inhibitory restraint of seizure activity in humans

    Nature Communications 3:1011–1060.

    https://doi.org/10.1038/ncomms2056

    • PubMed
    • Google Scholar
67. 1. Schevon CA
    2. Tobochnik S
    3. Eissa T
    4. Merricks E
    5. Gill B
    6. Parrish RR
    7. Bateman LM
    8. McKhann GM
    9. Emerson RG
    10. Trevelyan AJ

    (2019) Multiscale recordings reveal the dynamic spatial structure of human seizures

    Neurobiology of Disease 127:303–311.

    https://doi.org/10.1016/j.nbd.2019.03.015

    • PubMed
    • Google Scholar
68. 1. Schindler K
    2. Leung H
    3. Elger CE
    4. Lehnertz K

    (2007) Assessing seizure dynamics by analysing the correlation structure of multichannel intracranial EEG

    Brain 130:65–77.

    https://doi.org/10.1093/brain/awl304

    • PubMed
    • Google Scholar
69. 1. Schindler KA
    2. Bialonski S
    3. Horstmann MT
    4. Elger CE
    5. Lehnertz K

    (2008) Evolving functional network properties and synchronizability during human epileptic seizures

    Chaos 18:033119.

    https://doi.org/10.1063/1.2966112

    • PubMed
    • Google Scholar
70. 1. Schindler K
    2. Amor F
    3. Gast H
    4. Müller M
    5. Stibal A
    6. Mariani L
    7. Rummel C

    (2010) Peri-ictal correlation dynamics of high-frequency (80-200 hz) intracranial EEG

    Epilepsy Research 89:72–81.

    https://doi.org/10.1016/j.eplepsyres.2009.11.006

    • PubMed
    • Google Scholar
71. 1. Schroeder M
    2. Diehl B
    3. Chowdhury FA
    4. John S
    5. De TJ
    6. Trevelyan AJ
    7. Forsyth R
    8. Taylor PN
    9. Wang Y
    10. Schroeder GM
    11. Diehl B
    12. Chowdhury FA
    13. Duncan JS
    14. De TJ

    (2020) Erratum: Neuroscience biophysics and computational biology

    PNAS 117:e1922084117.

    https://doi.org/10.1073/pnas.1922084117

    • Google Scholar
72. 1. Shahabi H
    2. Taylor K
    3. Hirfanoglu T
    4. Koneru S
    5. Bingaman W
    6. Kobayashi K
    7. Kobayashi M
    8. Joshi A
    9. Leahy RM
    10. Mosher JC
    11. Bulacio J
    12. Nair D

    (2021) Effective connectivity differs between focal cortical dysplasia types I and II

    Epilepsia 62:2753–2765.

    https://doi.org/10.1111/epi.17064

    • PubMed
    • Google Scholar
73. 1. Smith EH
    2. Liou J
    3. Davis TS
    4. Merricks EM
    5. Kellis SS
    6. Weiss SA
    7. Greger B
    8. House PA
    9. McKhann GM
    10. Goodman RR
    11. Emerson RG
    12. Bateman LM
    13. Trevelyan AJ
    14. Schevon CA

    (2016) The ictal wavefront is the spatiotemporal source of discharges during spontaneous human seizures

    Nature Communications 7:11098.

    https://doi.org/10.1038/ncomms11098

    • PubMed
    • Google Scholar
74. 1. Solá L
    2. Romance M
    3. Criado R
    4. Flores J
    5. García del Amo A
    6. Boccaletti S

    (2013) Eigenvector centrality of nodes in multiplex networks

    Chaos 23:033131.

    https://doi.org/10.1063/1.4818544

    • PubMed
    • Google Scholar
75. 1. Struck AF
    2. Cole AJ
    3. Cash SS
    4. Westover MB

    (2015) The number of seizures needed in the emu

    Epilepsia 56:1753–1759.

    https://doi.org/10.1111/epi.13090

    • PubMed
    • Google Scholar
76. 1. Tadel F
    2. Baillet S
    3. Mosher JC
    4. Pantazis D
    5. Leahy RM

    (2011) Brainstorm: A user-friendly application for MEG/EEG analysis

    Computational Intelligence and Neuroscience 2011:879716.

    https://doi.org/10.1155/2011/879716

    • PubMed
    • Google Scholar
77. 1. Talairach J
    2. Bancaud J

    (1966) Lesion, irritative zone and epileptogenic focus

    Confinia Neurologica 27:91–94.

    https://doi.org/10.1159/000103937

    • PubMed
    • Google Scholar
78. 1. Tamilia E
    2. Park EH
    3. Percivati S
    4. Bolton J
    5. Taffoni F
    6. Peters JM
    7. Grant PE
    8. Pearl PL
    9. Madsen JR
    10. Papadelis C

    (2018) Surgical resection of ripple onset predicts outcome in pediatric epilepsy

    Annals of Neurology 84:331–346.

    https://doi.org/10.1002/ana.25295

    • PubMed
    • Google Scholar
79. 1. Tavakol S
    2. Royer J
    3. Lowe AJ
    4. Bonilha L
    5. Tracy JI
    6. Jackson GD
    7. Duncan JS
    8. Bernasconi A
    9. Bernasconi N
    10. Bernhardt BC

    (2019) Neuroimaging and connectomics of drug-resistant epilepsy at multiple scales: From focal lesions to macroscale networks

    Epilepsia 60:593–604.

    https://doi.org/10.1111/epi.14688

    • PubMed
    • Google Scholar
80. 1. Timofeev I
    2. Steriade M

    (2004) Neocortical seizures: Initiation, development and cessation

    Neuroscience 123:299–336.

    https://doi.org/10.1016/j.neuroscience.2003.08.051

    • PubMed
    • Google Scholar
81. 1. Trevelyan AJ

    (2009) The direct relationship between inhibitory currents and local field potentials

    The Journal of Neuroscience 29:15299–15307.

    https://doi.org/10.1523/JNEUROSCI.2019-09.2009

    • PubMed
    • Google Scholar
82. 1. Truccolo W
    2. Donoghue JA
    3. Hochberg LR
    4. Eskandar EN
    5. Madsen JR
    6. Anderson WS
    7. Brown EN
    8. Halgren E
    9. Cash SS

    (2011) Single-Neuron dynamics in human focal epilepsy

    Nature Neuroscience 14:635–641.

    https://doi.org/10.1038/nn.2782

    • PubMed
    • Google Scholar
83. 1. van Dellen E
    2. Douw L
    3. Baayen JC
    4. Heimans JJ
    5. Ponten SC
    6. Vandertop WP
    7. Velis DN
    8. Stam CJ
    9. Reijneveld JC

    (2009) Long-Term effects of temporal lobe epilepsy on local neural networks: A graph theoretical analysis of corticography recordings

    PLOS ONE 4:e8081.

    https://doi.org/10.1371/journal.pone.0008081

    • PubMed
    • Google Scholar
84. 1. van Diessen E
    2. Diederen SJH
    3. Braun KPJ
    4. Jansen FE
    5. Stam CJ

    (2013a) Functional and structural brain networks in epilepsy: What have we learned?

    Epilepsia 54:1855–1865.

    https://doi.org/10.1111/epi.12350

    • PubMed
    • Google Scholar
85. 1. van Diessen E
    2. Hanemaaijer JI
    3. Otte WM
    4. Zelmann R
    5. Jacobs J
    6. Jansen FE
    7. Dubeau F
    8. Stam CJ
    9. Gotman J
    10. Zijlmans M

    (2013b) Are high frequency oscillations associated with altered network topology in partial epilepsy?

    NeuroImage 82:564–573.

    https://doi.org/10.1016/j.neuroimage.2013.06.031

    • PubMed
    • Google Scholar
86. 1. van Diessen E
    2. Zweiphenning WJEM
    3. Jansen FE
    4. Stam CJ
    5. Braun KPJ
    6. Otte WM

    (2014) Brain network organization in focal epilepsy: A systematic review and meta-analysis

    PLOS ONE 9:e114606.

    https://doi.org/10.1371/journal.pone.0114606

    • PubMed
    • Google Scholar
87. 1. Van Eyndhoven S
    2. Hunyadi B
    3. Dupont P
    4. Van Paesschen W
    5. Van Huffel S

    (2019) Semi-automated EEG enhancement improves localization of ictal onset zone with EEG-correlated fmri

    Frontiers in Neurology 10:805.

    https://doi.org/10.3389/fneur.2019.00805

    • PubMed
    • Google Scholar
88. 1. Varotto G
    2. Tassi L
    3. Franceschetti S
    4. Spreafico R
    5. Panzica F

    (2012) Epileptogenic networks of type II focal cortical dysplasia: A stereo-EEG study

    NeuroImage 61:591–598.

    https://doi.org/10.1016/j.neuroimage.2012.03.090

    • PubMed
    • Google Scholar
89. 1. Warren CP
    2. Hu S
    3. Stead M
    4. Brinkmann BH
    5. Bower MR
    6. Worrell GA

    (2010) Synchrony in normal and focal epileptic brain: the seizure onset zone is functionally disconnected

    Journal of Neurophysiology 104:3530–3539.

    https://doi.org/10.1152/jn.00368.2010

    • PubMed
    • Google Scholar
90. 1. Weiss SA
    2. Banks GP
    3. McKhann GM
    4. Goodman RR
    5. Emerson RG
    6. Trevelyan AJ
    7. Schevon CA

    (2013) Ictal high frequency oscillations distinguish two types of seizure territories in humans

    Brain 136:3796–3808.

    https://doi.org/10.1093/brain/awt276

    • PubMed
    • Google Scholar
91. 1. Wendling F
    2. Bartolomei F
    3. Bellanger JJ
    4. Bourien J
    5. Chauvel P

    (2003) Epileptic fast intracerebral EEG activity: Evidence for spatial decorrelation at seizure onset

    Brain 126:1449–1459.

    https://doi.org/10.1093/brain/awg144

    • PubMed
    • Google Scholar
92. Book

    1. Wilcox R

    (2011) Modern statistics for the social and behavioral sciences

    CRC press.

    https://doi.org/10.1201/9781466503236

    • Google Scholar
93. Book

    1. Wyllie E
    2. Gidal BE
    3. Goodkin HP
    4. Loddenkemper T
    5. Sirven JI

    (2015)

    Wyllie’s treatment of epilepsy: Principles and practice (6th ed)

    Lippincott Williams & Wilkins.

    • Google Scholar
94. 1. Yaffe RB
    2. Borger P
    3. Megevand P
    4. Groppe DM
    5. Kramer MA
    6. Chu CJ
    7. Santaniello S
    8. Meisel C
    9. Mehta AD
    10. Sarma SV

    (2015) Physiology of functional and effective networks in epilepsy

    Clinical Neurophysiology 126:227–236.

    https://doi.org/10.1016/j.clinph.2014.09.009

    • PubMed
    • Google Scholar
95. 1. Zijlmans M
    2. Jiruska P
    3. Zelmann R
    4. Leijten FSS
    5. Jefferys JGR
    6. Gotman J

    (2012) High-Frequency oscillations as a new biomarker in epilepsy

    Annals of Neurology 71:169–178.

    https://doi.org/10.1002/ana.22548

    • PubMed
    • Google Scholar
96. 1. Zweiphenning WJEM
    2. van ’t Klooster MA
    3. van Diessen E
    4. van Klink NEC
    5. Huiskamp GJM
    6. Gebbink TA
    7. Leijten FSS
    8. Gosselaar PH
    9. Otte WM
    10. Stam CJ
    11. Braun KPJ
    12. Zijlmans GJM

    (2016) High frequency oscillations and high frequency functional network characteristics in the intraoperative electrocorticogram in epilepsy

    NeuroImage. Clinical 12:928–939.

    https://doi.org/10.1016/j.nicl.2016.09.014

    • PubMed
    • Google Scholar

Decision letter after peer review:

Thank you for submitting your article "Multilayer brain networks can identify the epileptogenic zone and seizure dynamics" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Michael Frank as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Mangor Pedersen (Reviewer #2); Christian Benar (Reviewer #3).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1. The prediction of the EZ seems to be completely within-sample. Please provide some indication of cross-validated out-of-sample prediction accuracy.

2. Please evaluate network properties with respect to appropriate null models. An important component of assessing the presence of network structure is to evaluate whether that same structure or nodal importance is not evident in null model where the underlying network structure is randomized. This is key to making sure that clustering results are not just sensitive to random fluctuations in the data. One suggestion is to use a phase-randomization of the voltage trace before computation of HFS to help better evaluate the null distribution of the mlEVC metric. A related reference is: https://pubmed.ncbi.nlm.nih.gov/22343126/

3. The authors may consider that at least one recent study has discussed the challenges and considerations of different recording and electrode sampling techniques on network metrics. And may explain discrepancies in the findings across studies: https://pubmed.ncbi.nlm.nih.gov/32537538/

4. The methods used in this study are complex, and rationale behind the series of steps used to perform the study are not provided. Given that this paper is primarily showcasing a novel method, the authors may want to consider incorporating more methodological details and rationale in-line with the Results.

5. Providing intuition behind abstract mathematical concepts would be helpful. For example (lines 139-141): Help the reader interpret the hierarchical unsupervised clustering applied to the left singular vectors of the concatenated and quantized mlEVC.

6. Line 96: Consider describing the concepts of graph/super-graph, and single/multi-layer in a general sense, first, as this is the first use of these terms in the manuscript.

7. Line 262: Typo, should "EZ-RnRZ" read "EZ-RnEZ".

8. The seizure evolution and state transition plots can be difficult to follow, due to the combination of colors and hues. Please consider a flow field or separate plots per seizure help demonstrate the overlapping and non-overlapping dynamics between seizures?

9. The authors may want to consider incorporating additional details in Table 1, pertaining to: (a) size of the network; (b) number of nodes overlapping with the resected tissue; (c) durations of the seizures; (d) seizure type.

10. One previous study motivated the application of phase and lagged-based connectivity metrics on the analytical amplitude of the broadband signal to avoid computing phase relationships within the asynchronous broadband range of the voltage signal. See: https://pubmed.ncbi.nlm.nih.gov/12631571/. Note that this is distinct from applying phase-based connectivity metrics directly to the high-frequency broadband component of the voltage signal or to high-frequency narrowbands (such as HFOs).

11. If I understand the methodical approach correctly, the authors generate an ordinal graph optimal when using fast intracranial electrophysiological data (500 or 1000Hz). It would be interesting to see a correlation plot between the current eigenvector centrality with an ordinal coupling parameter and cross-correlation of this data at multiple time-lags. In addition to tuning the coupling strength, the current work would be improved to quantify the optimal time interval between matrix 'layers'. It would be useful for the authors to consider/discuss this point in detail.

12. Related to the previous comment, it would be good to head further why the authors choose to use the Phase Lag Index for slow frequency estimation, a measure that does not operate proximately to the zero-lag.

13. A clearer rationale for using a spatiotemporal graph if warranted in this paper. The authors outline and justify the use of multilayer graphs from a methodological point of view, but it would also be good to hear why the authors choose this methodology from a clinical point of view. i.e., what can multilayer metrics tell us about epilepsy that single-layer cannot?

14. Related to the previous point, have the authors thought about comparing eigenvector centrality between multilayer and single-layer graphs, to see whether the results differ between the two?

15. Please provide more details about the filtering procedure, conducted before the Hilbert transform, including the type (and order) of filtering, how the DC offset was accounted for and whether the Hilbert transform in the given frequencies satisfy the Bedrossian's criteria for narrow-band frequency estimation (Xu and Yan, 2006).

16. It would have been nice to see in more detail examples from one or two patients with the actual resection area and the detected networks in volume. Figure 1 would benefit from an intuitive description of the method.

17. P7 It seems to me the term 'primary organization' has to be credited to Talairach and Bancaud. Also, early work on network organization were performed by Wendling and Bartolomei and could be cited (including the desynchronization at seizure onset in Wendling et al. 2003).

18. P8 it is very interesting to see that in a few patients one method fails whereas the other gives good results. This pleads for a complementarity of strategies to cope with inter-patient variability. Can the authors make hypotheses on which signal features/seizure patterns could explain the difference in results across methods and patients?

19. P11 « false positives from the EZ » this really depends on the definition of the EZ. If this is the minimum of cortex to be resected to render the patient seizure free, then the statement may hold. If this is the « primary region of organization of seizures » (definition of Bancaud, given p7), it is possible that not all these regions need to be resected. It would be useful to clarify the two definitions.

20. p17 it would be very interesting to perform a comparison of HF and LF synchrony in predicted the EZ. Moreover, as mentioned below too, part of the HF synchrony can arise from trains of sharp spikes that have energy in all frequencies, which would not be synchrony of high frequency oscillations but an emphasis on the sharp part of ictal spikes. This is not a problem per se, but more an issue of interpretation of the results.

21. Is data processing on monopolar or bipolar montage? Please clarify

22. Is it possible that part of the synchrony arises from filtered spikes? Can the authors provide a (normalized across frequencies) time-frequency representation of a representative seizure in order to assess the relative contribution of oscillations and spikes?

23. P25, eigenvalue centrality. While the rank of the matrix is likely to be T, wouldn't one expect some kind of temporal correlation across points (arising from both the high overlap and some actual temporal persistence of the networks), which would result in an elbow in the eigenvalues after some dimension? If it is not the case, and each time point is independent from the other ones (linear eigenvalue spectrum), what is the added power of multilayer method? Wouldn't a more simple method as summing the degrees across time be equivalent? Please clarify.

24. It would be interesting to visualize the eigenvalues of a representative example in order to measure this. Also, and importantly, it would be useful to compare the multilayer method to a simple mean degree across time (Wilke et al. 2010, van Mierlo et al. 2013, Courtens et al. 2016, Li et al. 2016, Balatskaya et al. 2020).

References:

Kramer, M.A., Eden, U.T., Kolaczyk, E.D., Zepeda, R., Eskandar, E.N., Cash, S.S., 2010. Journal of Neuroscience. 30, 10076-10085.

Schindler, K.A., Bialonski, S., Horstmann, M.-T., Elger, C.E., Lehnertz, K., 2008. Chaos 18, 033119.

Xu, Y and Yan, D, 2006. Proceedings of the American Mathematical Society 134, 2719-2728

Reviewer #1:

In this study, Shahabi and colleagues develop a new computational algorithm based on graph theory to study changes in brain network connections during drug-resistant, epileptic seizures. The structure and organization of these network connections have been the subject of many studies over the past two decades with the principal purpose of mapping seizures and identifying targets for resective surgery. A novel analysis of distributed brain network interactions between the local epileptogenic zone and remote, healthy brain areas outside of this zone is facilitated by invasive stereo-electroencephalographic (sEEG) recordings, which are a key asset to the present study. The authors construct time-dependent functional brain networks based on inter-areal synchrony in the high-frequency band and combine a new graph theory metric called multilayer eigenvector centrality (mlEVC) to quantify the changing pattern of connectivity amongst brain areas within and outside the epileptogenic zone. Using their innovative technique, the authors recapitulate a widely reported result that seizure progression invokes a robust alteration in network organization in which connectivity between the epileptogenic zone and healthy brain areas desynchronize. However, a key contribution of the present study is that the network nodes in the epileptogenic zone whose subsequent surgical resection led to seizure freedom could be predicted using unsupervised machine learning. These findings suggest that nodes involved in desynchronization during seizures may serve as putative surgical targets for epilepsy treatment. The authors also demonstrate that if patients are left untreated, then this abnormal desynchronization process during seizures intensifies with age and duration of epilepsy. This is a compelling scientific advancement in the field that begins to tackle the long-standing question of whether seizures beget seizures.

Generally, the study and analysis are presented as exploratory and proof-of-principle. A major strength of this study is the development of a new methodology to describe complex properties of seizures. However, I have concerns regarding possible overfitting of the data as no cross-validation nor testing was performed. Key claims could be better motivated with concrete hypotheses that are contextualized by prior work. My specific comments are detailed below:

1. Key questions posed in the study are provided and/or stated in the introduction but in many cases the importance or relevance in the broader context are not fully developed. Examples below:

a. Lines 50-52: "It has also been suggested that ictal periods can be delineated by a steady series of states(Burns et al., 2014), although whether this is true in all patients remains controversial."

b. Lines 52-55: "It remains unclear how the degree of desynchronization is correlated with physiological parameters such as age and duration of epilepsy(Van Diessen et 55 al., 2013a)"

2. Much of the introduction focuses on high frequency oscillations and pros/cons of their utility in localizing seizure foci. The authors present ictal high frequency synchronization as a phenomenon of interest, but to the general audience the differences between HFO and HFS are not concretely provided. The rationale for studying ictal HFS over HFOs is also not provided. Are these the same phenomenon? Are they distinct? Why do authors deem HFS analysis important for this study?

a. Lines 75-77: "While HFOs have been employed in analyzing functional(Schindler et al., 2010) and propagation(González Otárula et al., 2019) networks, the spatiotemporal dynamics of ictal high frequency synchronization (HFS) at macroscopic scales remain largely unknown."

b. Lines 238-240: A later point conflates HFOs and HFS as the same phenom: "Early-onset and late ictal HFOs have been considered biomarkers for seizure onset zone identification (Weiss et al., 2013), with the latter found to be a more reliable metric (Modur et al., 2011). Our EZ localization technique considers both features."

3. The interpretation of connectivity in the high-frequency broadband is unclear in the context of previous studies that demonstrate that the broadband activity is due to asynchronous, non-oscillatory neural firing (see references below). Consequently, it is unclear how one should interpret the notion of "synchronization" in the broad 80-200 Hz frequency range.

a. https://pubmed.ncbi.nlm.nih.gov/23283342/

b. https://pubmed.ncbi.nlm.nih.gov/29167419/

4. Was a null network model employed to ascertain whether the clustering procedure and the mlEVC metric identified target epileptogenic zone areas more reliably than chance?

5. What was the trade-off in the sensitivity and specificity of the EZ prediction algorithm? A receiver-operator curve analysis would help here.

6. One concern is that the analysis corresponding to the prediction of the epileptogenic zone is performed at the population level – aggregating nodes across all patients (line 160). What was the mean and standard error of the performance across patients?

7. The comparison between the present method and Fingerprinting method is of great value here. However, a clear conclusion regarding the comparison is not provided because the two algorithms provided results that only partially overlapped. How do these results compare to more conventional measures of the epileptogenic zone such as HFOs or spikes derived from the same dataset?

8. Results regarding the occurrence of different states of connectivity during the ictal period (lines 177-180) are certainly very interesting, but it is not clear how this finding advances previous studies:

a. https://pubmed.ncbi.nlm.nih.gov/20668192/

b. https://pubmed.ncbi.nlm.nih.gov/23366973/

9. A related worry regarding statistical power of the study is reflected in the stereotypy analysis (lines 192-195). To what extent is the stereotypy / lack of stereotypy in the network a function of the dimensionality of the feature space and overfitting of the data? Specifically, as the number of features / complexity of the model increases, are you more likely to find that each seizure event network is different from the others? A cross-validation and out-of-sample prediction approach would help mitigate these concerns.

10. The finding that stereotypy does not necessarily occur at high-frequencies is very interesting (lines 206-209), but the claim is not supported by any statistical testing. Furthermore, the explanation provided "the brain experiences divergent topologies, which might be the result of dissimilar EEG recordings" (line 209) is rhetorical as the topologies are based on features derived from the EEG recordings.

11. Why was the mid-seizure period used to normalize network measures for age/epilepsy duration-related analysis? How is this baseline period more advantageous than a pre-seizure baseline?

12. Were any other connection groups tested as being predictive of age or duration of epilepsy? In particular, is it possible that prolonged epilepsy might reorganize areas that are purely outside the EZ? Was there any relationship between the size of the EZ and age/duration?

13. Can the authors be sure that the emergence in low-frequency connectivity is not related to an increase in the rate or amplitude of ictal spiking? Are they independent phenomena?

14. How are the phases of the seizure (pre-termination, termination, early-seizure) defined? Are they based on the state analysis or on expert ratings of the seizures? A clear definition of these periods would be helpful.

15. Specific steps and choice of parameters in the unsupervised clustering pipeline described in the Methods are not justified. Such as

a. Usage of just the first four singular vectors

b. Limited combinations of the singular vectors

c. Determination of the rank T for the A matrix

d. How the mlEVC approach compares to the traditional EVC approach similar to that studies in the Burns et al. (2014).

Reviewer #2:

The authors combines intracranial EEG data with multilayer network modelling to delineate the seizure onset zone and spatiotemporal network activity during seizures. The major strength of this study is that the networks in this study incorporates spatial and temporal connectivity during seizures, and the relative importance of different intracranial electrodes. This is a novel and interesting approach that is likely to have an impact in the field. Below I discuss two topics (near zero-lag connectivity and previous research in the field) that is worth bearing in mind regarding the current manuscript.

- The authors generate an ordinal graph (connectivity between two adjacent time-points) to quantify spatiotemporal connectivity from intracranial electrophysiological data with a high sampling rate. An ordinal graph computes the connectivity between two adjacent time points, which may represent synchronicity close to a zero-lag in this data. It is worth bearing in mind graphs with connectivity data where time-points that are further apart, and how this may reflect neuronal connectivity at different time-scales.

- The current set of findings are interesting and they show good clinical accuracy at predicting the seizure onset zone. Spatio-temporal connectivity during seizures in this study also align with previous graph theoretical research using intracranial electrophysiological with single-layer graphs. Previous studies show distinct time-varying peri-seizure changes in the clustering coefficient and path length metrics. During seizure initiation and propagation, there is an increase of clustering coefficient and shortest path length, resembling a regularized, or isolated, network topology (Kramer et al., 2010; Schindler et al., 2008). In a network that is regularized, local clusters of neurons may isolate themselves from the rest of the network. It would be good to evaluate the current findings in the context of previous graph theoretical research using intracranial electrophysiological data in epilepsy.

Reviewer #3:

This manuscript introduces a new method for identifying the epileptogenic zone from intracerebral EEG data, based on graph measures in high frequency bands, multilayer graph analysis and clustering.

Graph measures and multivariate methods are promising tools in the electrophysiological characterization of the epileptogenic zone; the strategy proposed here falls within this timely topic. The strength of the approach is to be fully automatic, and to rely on multivariate measures that can capture the overall structure of the data better than the usual monovariate measures. The method identified electrodes within the resected area in 88% of the patients, with only few contacts detected outside of the resection.

It could be interesting to also test the method versus the clinician EZ, as the resected area may be larger this would be consistent with the fact that the 'target' area found by clustering is much smaller that the resected contacts. This would be a further probe of the sensitivity of the method. Another interesting measure would be to probe whether the proportion of contacts outside the resected area increase in non seizure-free patients. Yet another interesting test would be to see whether the multivariate method outperforms more classical measures such as graph strength. In other words, it what does the multilayer approach add to the single layer measures?

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Multilayer brain networks can identify the epileptogenic zone and seizure dynamics" for further consideration by eLife. Your revised article has been evaluated by Michael Frank (Senior Editor) and a Reviewing Editor.

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below. Please address these comments comprehensively in your response.

1. In their response and amendments to the main text, the authors claim that their method does not require cross-validation and testing over independent datasets since the approach does not involve training the model parameters across patients and does not access information regarding the resection zone. While it is true that the method is applied to each patient, separately, the methodological choices including, but not limited to, functional connectivity metrics, frequency band selection, percentile threshold selection, still effectively "train" the general algorithm and may be highly specific to the cohort studied here. This issue presents itself as a driver of inter-patient variability, as demonstrated by the conflicting findings between the mlEVC method and Fingerprint method. While this limitation detracts from the methodological advances put forth by this study, but some caution in making the claim that cross-validation/testing is not necessary is warranted. Specifically, the manuscript should include a rationale for why cross-validation was not applied in this context, coupled with an acknowledgement that such an approach in future will be important to verify the current results.

2. The relationship between the coupling parameter and network time-scale should be further elaborated, given the importance of this parameter to the authors' claim that patients may have different rates of network change. Please present data examining how sensitive the findings of the EZ are to the choice of this parameter. Please also provide examples of how different coupling parameters would yield more insight into the different time-scales of network change

3. The findings comparing a single-layer model to the multi-layer model should be explicitly quantified to concretely justify the claim in the Discussion section that single-layer models yielded poorer EZ identification than multi-layer models. How much more accurate was the multi-layer model to the single-layer model?

4. More attention should be paid in the text to the risk of contamination of high frequency synchrony by two processes: (1) Filtering of spikes; and (2) Harmonics of non-sinusoidal oscillations. These topics are addressed in the following study, which could be cited in support of the discussion.

Pitfalls of high-pass filtering for detecting epileptic oscillations: a technical note on "false" ripples. Bénar CG, Chauvière L, Bartolomei F, Wendling F. Clin Neurophysiol. 2010 Mar;121(3):301-10

5. The time frequency figures should be accompanied by a presentation of the signal time course in order to fully appreciate the presence/absence of spikes. While there indeed seems to be spikes in the preictal period (vertical lines) , they present much less energy than the high frequency activity during the seizure. This large increase is quite blurred and indistinct and it is not clear what this actually represents. The original signal time course would help understand which seizure pattern this corresponds to.A consideration of the potential effect of harmonics would also be helpful.

Essential revisions:

1. The prediction of the EZ seems to be completely within-sample. Please provide some indication of cross-validated out-of-sample prediction accuracy.

We were inspired by (Burns et al., 2014) while working on our hypothesis. They state, based on exploratory observations, that:

“We found that, in patients with positive surgical outcomes (i.e., patients who were seizure-free after surgery; Table 1), the focal areas were the least connected in the network (least important) shortly after the onset time of the seizure and we denoted this state as the “isolated focus” (IF) state. We also found that, in many of these patients, the focal areas became most connected in the network in a brain state that occurred during the middle or end of the seizures, and we defined this state as the “connected focus” (CF) state.”

They used the clinical resection information to investigate the dynamics of eigenvector centrality in the seizure onset zone. They observed specific patterns of centrality among different groups of patients. However, (Burns et al., 2014) did not introduce an algorithm that directly predicts the seizure onset zone. We believed that there is a gap in this part of literature and an automated algorithm could be developed to help predict the epileptogenic zone (EZ) from ictal network measures separately in each subject.

We therefore proposed an unsupervised framework that aims to identify the epileptogenic zone without any knowledge about the actual resection volume. Based on the observations in (Burns et al., 2014) and other studies, we hypothesized that the EZ should have a unique pattern of brain connectivity. Our algorithm for EZ identification is the same for all patients and does not involve training using resection labels. Since the algorithm was not trained across subjects, there is no need for a separate out-of-sample validation.

The steps employed for unsupervised clustering are common in the literature. For example, thresholding and quantization have been widely used in different brain network studies. For selecting a proper coupling parameter in our multilayer model, we followed the guidelines in (Betzel and Bassett, 2017).

“Without strong evidence to select one weighting scheme over another, interlayer links are usually assigned the same value, ω, that is sometimes varied over a narrow range. Ideally, there would be a principled, data-driven approach for selecting this value.”

Considering neural mechanisms with different timescales, we hypothesized that the proper interlayer parameter might be different among patients. As a result, in our data-driven approach, we first defined a performance function (Halkidi, 2001) for the unsupervised clustering algorithm linkage in MATLAB. As we describe in the method section, the performance function was defined based on our hypothesis that nodes inside the “presumptive EZ” should have a similar pattern of brain connectivity. For the final analysis, we combined the results of performance function over the full range of values for thresholding and interlayer coupling parameter, in an identical manner for each subject. These ranges were held constant among all patients.

To clarify this point in the paper, we added the following to the Results section:

“Cross-validation and testing over a separate independent set were not required for our EZ identification technique since the unsupervised approach was applied independently to each subject without training within or across subjects. Furthermore, resection labels were not employed as part of the prediction algorithm. The same performance function and ranges of coupling and thresholding values were used for all subjects.”

In addition to examining the accuracy of our method on 16 patients, we also performed several other analyses to evaluate validity of our technique, which are described later in this response, with associate modifications in the manuscript. Finally, we note that in the initial version of this manuscript, we used “classification” or “classify” interchangeably with “grouping” or “clustering”. We have clarified the text to make it clearer that all clustering was performed in an unsupervised manner.

2. Please evaluate network properties with respect to appropriate null models. An important component of assessing the presence of network structure is to evaluate whether that same structure or nodal importance is not evident in null model where the underlying network structure is randomized. This is key to making sure that clustering results are not just sensitive to random fluctuations in the data. One suggestion is to use a phase-randomization of the voltage trace before computation of HFS to help better evaluate the null distribution of the mlEVC metric. A related reference is: https://pubmed.ncbi.nlm.nih.gov/22343126/

We consulted the suggested reference (Zalesky et al., 2012) and a related paper for generating surrogate data by time series randomization using Fourier transform and phase scrambling (Prichard and Theiler, 1994). The latter article was employed to create surrogate data for several patients, using their SEEG signals during seizures. The mlEVC metric for phase-randomized signals was calculated with several coupling parameters. We compared results of original times series and surrogate data. Results are presented in new Figure 2—figure supplement 2 – see added text below. The left column portrays the mlEVC for the original data and the right column depicts the same metric when the phase scrambled data was used. We sorted nodes such that those inside the resected volume are on the top of the figures. We observe a clear distinction between the two scenarios, in which the mlEVC of the original time series displays characteristic patterns both in resected and non-resected regions that are not apparent in the scrambled data.

We added this paragraph to the results:

“To explore the validity of our approach, we constructed a null model by also computing the mlEVC from phase-randomized SEEG signals (Prichard and Theiler, 1994). In figure 2—figure supplement 2, we show typical results of the mlEVC measures calculated from the original time-series and phase-randomized data for a single subject. These result show significant differences between original and randomized data in which the characteristic patterns of brain connectivity both in resected and non-resected areas are lost when signals are phase-randomized. We performed this analysis for different patients and seizures with similar findings.”

3. The authors may consider that at least one recent study has discussed the challenges and considerations of different recording and electrode sampling techniques on network metrics. And may explain discrepancies in the findings across studies: https://pubmed.ncbi.nlm.nih.gov/32537538/

In the introduction and discussion, we note possible sources of discrepancies among the results of studies on ictal connectivity. In particular, based on the reviewers’ suggestion we added the following:

“Recent research has investigated the effect of incomplete electrode sampling on network metrics (Conrad et al., 2020), in which eigenvector centrality showed the highest reliability in comparison to other network metrics, when randomly removing network nodes. However, they also showed that confidence intervals in network measures can vary significantly between patients. For this reason, care should be taken when interpreting and comparing different studies.”

4. The methods used in this study are complex, and rationale behind the series of steps used to perform the study are not provided. Given that this paper is primarily showcasing a novel method, the authors may want to consider incorporating more methodological details and rationale in-line with the Results.

5. Providing intuition behind abstract mathematical concepts would be helpful. For example (lines 139-141): Help the reader interpret the hierarchical unsupervised clustering applied to the left singular vectors of the concatenated and quantized mlEVC.

6. Line 96: Consider describing the concepts of graph/super-graph, and single/multi-layer in a general sense, first, as this is the first use of these terms in the manuscript.

The reviewers suggested adding more description regarding the methodology in the results and introduction section of this paper. Before explaining our changes in detail, we note that the proposed methodology in this paper is dominantly focused on EZ identification. Following the eLife format, we incorporated the detailed algorithm description at the end of the paper. In Response 1 above, we elaborated how these steps were chosen based on published studies and common practice in the current literature. To improve the readability of the manuscript, we made the following modifications:

In the introduction, we introduce mathematical concepts including graphs/super-graphs/multilayer.

“Graph analysis provides a mathematical framework for quantification of brain connectivity (Bullmore, 2009). Brain networks may be represented as a graph, ${\displaystyle G=(V,E)}$,­ in which nodes (${\displaystyle V}$) characterize anatomical regions or electrodes, and edges (${\displaystyle E}$) reflect structural or functional connections among them. Traditionally, brain connectivity in each time sample (layer) of dynamic networks has been evaluated independently, via single graph analysis. However, multilayer analysis allows us to model the entire data with a single super-graph, in which individual graphs for each time-sample are linked. Assuming there are T single graphs (time-samples) with N nodes in each, the super-graph has NT nodes.”

In the Results section, we added references to the eigenvector centrality (EVC) and its formulation. EVC usage in ictal studies and its limitation is also discussed. We then describe why and how a new measure (mlEVC) has been defined based on our multilayer model:

“As discussed in the introduction, we aimed to model ictal brain connectivity by multilayer networks. Additionally, we were interested in quantifying network dynamics using centrality metrics (Rubinov and Sporns, 2010). Among those measures, eigenvector centrality, a rank-based metric which assesses the importance of each node in the network (Bonacich, 1972), has been employed in studying seizures (Burns et al., 2014). Mathematically, the leading eigenvector of the adjacency (connectivity) matrix has been assigned as the eigenvector centrality of the graph when there is one clique (component) in the matrix (Bonacich, 1972). However, in our multilayer model of time-varying brain connectivity, a single vector cannot explain the complex structure of the spatiotemporal networks. Therefore, we introduced a new measure called multilayer eigenvector centrality (mlEVC) that incorporates the top T eigenvectors of the adjacency matrix of the super-graph (see Methods). This allows us to evaluate patterns of nodal centrality and identify regions with similar connectivity characteristics to the rest of the graph. Further, because mlEVC is a function of the interlayer coupling parameter (c), we can explore neural processes at different timescales by varying c.”

In the results, a new section entitled “Algorithm for predicting EZ using mlEVC” was inserted to explain our algorithm, the rationale behind each step, and the way it was executed:

“We hypothesize that the EZ can be identified as the set of nodes in the graph that exhibit a characteristic and distinct pattern of connectivity to other areas during the seizure. To explore this question, we first quantized the mlEVC of each seizure into three levels based on percentile thresholding (d). The top d/2-portion of elements were assigned a value of ‘1’, the bottom d/2-portion a value of ‘-1’, and the remainder a value of ‘0’. For each subject, quantized measures of mlEVC for two high-frequency bands and all seizures were concatenated to a single matrix with dimension N by T_tot (twice the total number of sample points) (see Methods). We applied the singular value decomposition (SVD) to the concatenated matrix and used the left singular vectors (${\displaystyle u_i\in {\mathbb{R}}^N}$) to identify nodes (SEEG channels) with similar features. As an illustration, Figure 2c depicts ${\displaystyle u_1}$, ${\displaystyle u_2}$, and ${\displaystyle u_3}$ for patient 17.

Next, we applied an unsupervised clustering algorithm using the ${\displaystyle u_i}$ vectors to detect a target cluster that represents the EZ. Following our initial hypothesis, the target cluster should portray a dense and distinctive set of nodes in the feature space with a significant distance from nodes in the non-target group (Figure 2c). We describe the clustering algorithm in detail in the method section. Briefly, we designed a data-driven framework to cluster nodes into two groups using an agglomerative hierarchical clustering technique (function linkage in MATLAB). This process was performed for different combinations of ${\displaystyle u_i}$ as features of the clustering algorithm, and a range of values for c and d, resulting in 440 clustering runs. We weighted each run using a performance function (Halkidi, 2001) that examines the tightness of the target cluster and its separation from other nodes. Finally, using the weighted sum of performance metric for all runs, SEEG contacts were divided into two groups; target and non-target.”

In the methods, we added the following description:

“Historically, the EVC is defined as a weighted sum of all eigenvectors of the matrix, but simplified to the leading eigenvector for matrices with one clique (Bonacich, 1972). In our model of time-varying brain connectivity, this simplification does not hold because of the time-varying nature of connectivity during the seizure that is embodied in the super-adjacency matrix.”

7. Line 262: Typo, should "EZ-RnRZ" read "EZ-RnEZ".

Thanks for pointing this out. Corrected.

8. The seizure evolution and state transition plots can be difficult to follow, due to the combination of colors and hues. Please consider a flow field or separate plots per seizure help demonstrate the overlapping and non-overlapping dynamics between seizures?

We updated the state transition plots based on your recommendation. Please check figure 3d and figure 3—figure supplement 1.

9. The authors may want to consider incorporating additional details in Table 1, pertaining to: (a) size of the network; (b) number of nodes overlapping with the resected tissue; (c) durations of the seizures; (d) seizure type.

We updated Table 1 with the requested information for parts a, b, and c. Unfortunately, we do not have the seizure type in our database.

10. One previous study motivated the application of phase and lagged-based connectivity metrics on the analytical amplitude of the broadband signal to avoid computing phase relationships within the asynchronous broadband range of the voltage signal. See: https://pubmed.ncbi.nlm.nih.gov/12631571/. Note that this is distinct from applying phase-based connectivity metrics directly to the high-frequency broadband component of the voltage signal or to high-frequency narrowbands (such as HFOs).

We understand the reviewers’ questions about the terminology, relation of HFOs and high-frequency synchronization and applying phase-based connectivity metric on broad-band signals.

To clarify the terminology and justify the methods and approach used in our analysis, we first added the following description to the introduction.

“The term HFOs has been attributed to brain activity, with multiple possible physiological and pathological neural mechanisms, between 80-500 Hz (Jacobs et al., 2012) or 30-600 Hz (Engel and da Silva, 2012). This includes high-γ neural activities (80-200 Hz) (Ray and Maunsell, 2011) and broad-band high frequency (Arnulfo et al., 2020).”

Although some epilepsy studies have investigated HFOs with specific patterns (e.g. ripples and fast ripples in the time-domain), the broad definition of HFOs remains unchanged. Our analysis of brain connectivity uses the standard HFO definition and is computed using a synchrony metric.

Regarding use of phase-based measures on broad-band signals we added references to several studies that employ them (see Results section):

“Several studies have applied phase-based connectivity metrics on broad-band high-frequency oscillations. For instance, (Zweiphenning et al., 2016) computed phase lag index (PLI) in two high-frequency bands: ripple (80-250 Hz) and fast ripple (250–500 Hz). Another study, used PLI to compute brain connectivity in broadband 80-250 Hz (Nissen et al., 2016). Last, (Burns et al., 2014) used coherence to investigate ictal networks in γ (25-90 Hz) band.”

We also note, as shown in (Aydore et al., 2013) for intracranial recordings, phase-based measures provide essentially the same information as coherence in most cases.

Importantly, we first divided the 80-200 Hz into two narrower frequency ranges (80-140 Hz, and 140-200 Hz) and then computed the connectivity matrices based on the lagged-coherence metric for each frequency range separately. The mlEVC measure was computed for each frequency range separately and the results then concatenated.

11. If I understand the methodical approach correctly, the authors generate an ordinal graph optimal when using fast intracranial electrophysiological data (500 or 1000Hz). It would be interesting to see a correlation plot between the current eigenvector centrality with an ordinal coupling parameter and cross-correlation of this data at multiple time-lags. In addition to tuning the coupling strength, the current work would be improved to quantify the optimal time interval between matrix 'layers'. It would be useful for the authors to consider/discuss this point in detail.

13. A clearer rationale for using a spatiotemporal graph if warranted in this paper. The authors outline and justify the use of multilayer graphs from a methodological point of view, but it would also be good to hear why the authors choose this methodology from a clinical point of view. i.e., what can multilayer metrics tell us about epilepsy that single-layer cannot?

14. Related to the previous point, have the authors thought about comparing eigenvector centrality between multilayer and single-layer graphs, to see whether the results differ between the two?

We respond to these three comments together since they are highly relevant. In summary the reviewers requested (a) A discussion about the optimal time interval between different layers of the network, (b) Comparison of results between single-layer and multilayer networks, and (c) A rationale for employing a multilayer approach in epilepsy specifically during seizures.

The common approach to analyzing time-varying brain connectivity is using sliding windows with a constant window size and overlap between windows across time. If the interval between windows is small, instantaneous fluctuations will dominate processes with slower timescales. In the case of large intervals between windows, fast dynamics might be discarded. As a result, choosing an optimal interval, is a challenging issue. Applying the same interval for all participants will exacerbate the problem since seizure dynamics are highly variable across subjects.

In our multilayer model of brain networks, we use a 500ms time interval between layers to allow us to capture fast seizure dynamics. By increasing the coupling parameter among adjacent layers, we can incorporate processes with slower time scales. Using this approach, we only need to tune the coupling parameter to extract desired spatiotemporal processes. So, rather than finding an optimal time interval for the multi-layer network, we instead focus on finding an appropriate coupling parameter. We use a data-driven approach to evaluate the cost function for EZ identification over a range of values for the coupling parameter (See Response #1).

We added the following explanation to the discussion.

“By adjusting the coupling between layers, we can overcome the current challenges for selecting an appropriate time interval between different samples of a temporal network. Here, we assumed a relatively short time interval between layers (500ms) to have a high temporal resolution. By increasing the coupling parameter between adjacent layers, we can investigate processes with slower timescales.”

To assert the importance of a multilayer approach when compared to single-layer, we inserted this paragraph into the results.

“Furthermore, we examined the effect of multilayer modeling and adjusting the coupling parameter on EZ prediction and super-graph structure. In figure 2—figure supplement 1, we display the left singular vectors of mlEVC. For illustration purposes, we showed projections with respect to the singular vectors for the parameter set (c,d) for which the prediction of EZ (binary labeling) was closest to the assumed ground truth (resection information). We observed that the optimal coupling value is different between patients, which supports our primary rationale for choosing a multilayer framework with an adaptive interlayer coupling parameter. In our initial analysis, we ran the algorithm for the single-layer case and EZ identification results were poor.”

After discussing the methodological reasons for using a multilayer framework, we also explained the neurophysiological rationales. The details as appeared below were added to the introduction:

Additionally, there exist several neurological rationales for employing a multilayer approach in seizure analysis. First, the concept of dynamic network reconfiguration has been studied in brain networks (Bassett et al., 2011; Braun et al., 2015). Previous research has shown state transitions during seizures, either through brain connectivity analysis (Burns et al., 2014) or microelectrode recordings (Smith et al., 2016). Multilayer networks have the capability to delineate these transitions and identify network reconfiguration (Mucha et al., 2010). Second, electrophysiological signals are highly non-stationary during ictal periods. As a result, traditional analysis of time-varying networks based on isolated graphs would be affected by instantaneous fluctuations rather than the underlying spatiotemporal networks. Third, seizure propagation is one of the key elements of ictal activity. Recent work has indicated the importance of multilayer modeling of complex systems when encountering spreading processes (De Domenico et al., 2016).

12. Related to the previous comment, it would be good to head further why the authors choose to use the Phase Lag Index for slow frequency estimation, a measure that does not operate proximately to the zero-lag.

In addition to the brief justification for choosing this measure in the initial version of the manuscript, we added the following description to the result section:

“The PLI is recognized for its superiority to measures such as phase locking value (PLV). By discarding interactions with a phase difference of zero, the PLI quantifies phase coupling between two signals while excluding volume conduction as a confounding factor. Compared to other functional connectivity metrics, PLI is therefore less susceptible to common sources. Since a major part of low-frequency (<50Hz) brain activity during seizures is derived from ictal discharges, which exhibit distance-dependent delays (Smith et al., 2016), true neurological zero-lag connectivity is highly unlikely in this scenario. Because of its robustness to volume conduction, PLI has been widely used in studies of low-frequency (<50 Hz) brain connectivity in epilepsy (Nissen et al., 2016; Schevon et al., 2012; Van Diessen et al., 2013; Zweiphenning et al., 2016).”

15. Please provide more details about the filtering procedure, conducted before the Hilbert transform, including the type (and order) of filtering, how the DC offset was accounted for and whether the Hilbert transform in the given frequencies satisfy the Bedrossian's criteria for narrow-band frequency estimation (Xu and Yan, 2006).

21. Is data processing on monopolar or bipolar montage? Please clarify.

Thanks for bringing up these important issues. Regarding the signals montage and preprocessing, we added the following paragraph to the “Patients and recordings” subsection in the methods:

“We preprocessed SEEG signals before further analyses. First, we constructed bipolar channels by subtracting unipolar signals recorded from each pair of adjacent contacts. For analysis we included bipolar channels only if both contacts were inside gray matter. Next, any DC offset was removed from each bipolar channel. Thereafter, using the Brainstorm software (Tadel et al., 2011), bipolar EEG signals were bandpass filtered in the 3-200 Hz range and notch filtered at 60, 120, and 180 Hz. Bandpass filtering was performed using a Kaiser-window linear phase FIR filter of order 7,252 (for a sample rate of 1000 Hz), using the fir1 MATLAB function. We compensated for filter-induced delay by shifting the filtered sequence, resulting in a zero-delay filter. The notch filtering used an IIR filter of order 6 with a bandstop of 1 Hz in notch frequencies. We used the filtfilt function in MATLAB to compensate for group delay. The resulting preprocessed data were employed in all analyses in this paper. When analyzing low-frequency and high-frequency brain connectivity, we bandpassed the preprocessed data using the same filter type (Kaiser-window linear phase FIR filter) with different ranges and orders.”

In Author response image 1 we show the frequency response of the bandpass and notch filters and their characteristics. The low transition band (2.5-3 Hz) is also magnified for a better visualization. In the early version of this manuscript, we listed (2-200 Hz) as the range of the bandpass filter. This has been corrected to (3-200 Hz) in the revised manuscript. The data were filtered over a sufficiently long interval either-side of the seizure window that transient effects of filters was not a matter of concern. Author response image 2 demonstrates the frequency response for the notch filter.

As stated in response #10, we bandpassed preprocessed signals in two narrower frequency ranges, (80-140 Hz) and (140-200 Hz), before applying Hilbert transform, satisfying the Bedrosian identity (Xu and Yan, 2006).

Author response image 1

Download asset Open asset

Author response image 2

Download asset Open asset

16. It would have been nice to see in more detail examples from one or two patients with the actual resection area and the detected networks in volume. Figure 1 would benefit from an intuitive description of the method.

We appreciate your suggestion. Unfortunately, we did not have access to MRIs with co-registered electrodes and resected volumes while preparing this response letter.

17. P7 It seems to me the term 'primary organization' has to be credited to Talairach and Bancaud. Also, early work on network organization were performed by Wendling and Bartolomei and could be cited (including the desynchronization at seizure onset in Wendling et al. 2003).

Thanks for the suggestions. We now cite the following publications in the article.

Talairach J, Bancaud J. 1966. Lesion, “irritative” zone and epileptogenic focus. Confin Neurol 27:91–94. doi:10.1159/000103937

Wendling F, Bartolomei F, Bellanger JJ, Bourien J, Chauvel P. 2003. Epileptic fast intracerebral EEG activity: Evidence for spatial decorrelation at seizure onset. Brain 126:1449–1459. doi:10.1093/brain/awg144

18. P8 it is very interesting to see that in a few patients one method fails whereas the other gives good results. This pleads for a complementarity of strategies to cope with inter-patient variability. Can the authors make hypotheses on which signal features/seizure patterns could explain the difference in results across methods and patients?

The result section mentioned that insufficient data (short seizures) might be one reason for mlEVC's failure. Here, we added another analysis to look for signal features. For patient 8, mlEVC could not identify any nodes inside the resected area as EZ, while Fingerprint predicted two nodes as EZ (one false positive). We looked at the time-frequency spectrum and mlEVC for one seizure of this subject, which is illustrated in Author response image 3. In the upper figure, we can see the main features of the Fingerprint algorithm (suppression, fast activity, and pre-ictal spikes) in the predicted nodes. However, we do not observe any significant HFOs for those nodes (see also our response to comments #20 & #22 regarding HFOs). As a result, the mlEVC does not have a distinctive pattern in any nodes in the resected region (Author response image 3, lower). Lack of dominant HFOs in the EZ can be another factor attributes to mlEVC's failure.

On the other hand, Fingerprint can fail for different reasons. First and foremost, its features might not be present for some seizures or patients. Second, since Fingerprint uses a trained model for predicting the EZ in the new patient, the new features can be quite different from the base model, or the detection process falls short of the presumed threshold for EZ identification.

Seizure types and locations can also affect the results in both methods, but it is hard to draw conclusions with a small number of patients in each cohort. Additional studies on a larger dataset are needed to further investigate these issues, which is beyond the scope of the current paper.

Author response image 3

Download asset Open asset

19. P11 « false positives from the EZ » this really depends on the definition of the EZ. If this is the minimum of cortex to be resected to render the patient seizure free, then the statement may hold. If this is the « primary region of organization of seizures » (definition of Bancaud, given p7), it is possible that not all these regions need to be resected. It would be useful to clarify the two definitions.

We use the first definition since with this, we can definitively state that the EZ was not resected in those patients who are not seizure free. The EZ is discussed in the Introduction as follows:

“In using SEEG to identify the epileptogenic zone (EZ) we attempt to take into account the earliest ictal EEG change based on an anatomoelectroclinical analysis (Kahane et al., 2006). This concept incorporates both the anatomic region that initiates the epileptic discharge as well as the “primary organization” (Talairach and Bancaud, 1966) that leads to the manifestation of the clinical seizure itself (Wyllie et al., 2015). The gold standard method of confirming EZ localization is based on whether seizure freedom has been achieved by resection or ablation. However, the actual ground truth for EZ location, in the sense of the minimal resection volume needed to achieve seizure freedom, is unknown since in many cases the resection volumes may extend well beyond the EZ.”

20. p17 it would be very interesting to perform a comparison of HF and LF synchrony in predicted the EZ. Moreover, as mentioned below too, part of the HF synchrony can arise from trains of sharp spikes that have energy in all frequencies, which would not be synchrony of high frequency oscillations but an emphasis on the sharp part of ictal spikes. This is not a problem per se, but more an issue of interpretation of the results.

22. Is it possible that part of the synchrony arises from filtered spikes? Can the authors provide a (normalized across frequencies) time-frequency representation of a representative seizure in order to assess the relative contribution of oscillations and spikes?

We appreciate your suggestion on using LF connectivity features to predict EZ and comparing the results with HF. While that can be a valuable investigation, we chose HF features based on the literature. Mainly, (Schevon et al., 2012) showed that core and penumbral might share similar low-frequency activity patterns. However, the high-frequency oscillations are the unique characteristics of the core. Importantly, our findings in Figures 4 and 6 in the main article presented that EZ-nR connectivity in mid-seizure has changed significantly for HF compared to the pre-ictal period. At the same time, this change was not significant in LF.

To investigate the contribution of ictal spikes and oscillations in SEEG signals, we computed the time-frequency representation (normalized across frequencies) of all signals in several seizures and patients. The following discussion and related figure (Figure 2—figure supplement 3) were added to the Results section of the manuscript.

“To investigate whether apparent high-frequency networks were produced as a result of ictal spikes rather than true oscillations, we examined the time-frequency representation of seizures. Figure 2—figure supplement 3 displays time-frequency plots for two sample seizures. The upper plot shows the spectrum averaged among nodes predicted as EZ, and the lower spectrum shows the case for randomly selected nodes outside the resection region, with the number of nodes in both groups equal. We observed pre-ictal spikes in the time-frequency plot for nodes inside the EZ, a previously reported feature of the epileptogenic zone (Grinenko et al., 2018). However, the figure shows that non-spiking high-frequency oscillations were the dominant activity after seizure onset. This observation supports our hypothesis that HFOs, rather than ictal spiking, are the main contributor to observed synchronization during seizures.”

23. P25, eigenvalue centrality. While the rank of the matrix is likely to be T, wouldn't one expect some kind of temporal correlation across points (arising from both the high overlap and some actual temporal persistence of the networks), which would result in an elbow in the eigenvalues after some dimension? If it is not the case, and each time point is independent from the other ones (linear eigenvalue spectrum), what is the added power of multilayer method? Wouldn't a more simple method as summing the degrees across time be equivalent? Please clarify.

24. It would be interesting to visualize the eigenvalues of a representative example in order to measure this. Also, and importantly, it would be useful to compare the multilayer method to a simple mean degree across time (Wilke et al. 2010, van Mierlo et al. 2013, Courtens et al. 2016, Li et al. 2016, Balatskaya et al. 2020).

Previously in this response letter, we scrutinize the virtues of a multilayer framework in response to comments (#11, #13, and #14). We mentioned why a single-layer method cannot extract the neural processes with timescales slower than the time intervals between layers. Our EZ identification method showed that by increasing the coupling between layers we can explore across temporal scales.

We performed the requested analysis for the magnitude of the eigenvalues based on different coupling parameters. We inserted the following discussion and figure in the Results section:

“Figure 2—figure supplement 4 presents the changes in super-graph eigenvalues by adjusting the coupling parameter. For c≤1, there is a falloff when the number of eigenvalues (n_e) meets the number of layers (T), indicating a super-adjacency matric with effective rank T. This is an expected result since the coupling is relatively small. When c increases, the eigenvalues become larger, and their corresponding eigenvectors would comprise several neighboring layers. Although the rate of decline in magnitude of eigenvalues accelerates as c increases, falloff can be detected when n_e ≤ 2T. By increasing the coupling value, the super-adjacency matrix transforms from a block diagonal matrix with effective rank T, to a matrix with major non-diagonal blocks and an effective rank far greater than T. In the computation of mlEVC, we considered the top T eigenvectors for all coupling parameters to avoid erroneous assumptions about the rank of super-graphs.”

References

Arnulfo, G., Wang, S.H., Myrov, V., Toselli, B., Hirvonen, J., Fato, M.M., Nobili, L., Cardinale, F., Rubino, A., Zhigalov, A., Palva, S., Palva, J.M., 2020. Long-range phase synchronization of high-frequency oscillations in human cortex. Nat. Commun. 11, 1–15. doi:10.1038/s41467-020-18975-8

Aydore, S., Pantazis, D., Leahy, R.M., 2013. A note on the phase locking value and its properties. Neuroimage 74, 231–244.

Bassett, D.S., Wymbs, N.F., Porter, M.A., Mucha, P.J., Carlson, J.M., Grafton, S.T., 2011. Dynamic reconfiguration of human brain networks during learning. Proc. Natl. Acad. Sci. U. S. A. 108, 7641–7646. doi:10.1073/pnas.1018985108

Betzel, R.F., Bassett, D.S., 2017. Multi-scale brain networks. Neuroimage 160, 73–83. doi:10.1016/j.neuroimage.2016.11.006

Bonacich, P., 1972. Factoring and weighting approaches to status scores and clique identification. J. Math. Sociol. 2, 113–120. doi:10.1080/0022250X.1972.9989806

Braun, U., Schäfer, A., Walter, H., Erk, S., Romanczuk-Seiferth, N., Haddad, L., Schweiger, J.I., Grimm, O., Heinz, A., Tost, H., Meyer-Lindenberg, A., Bassett, D.S., 2015. Dynamic reconfiguration of frontal brain networks during executive cognition in humans. Proc. Natl. Acad. Sci. U. S. A. 112, 11678–83. doi:10.1073/pnas.1422487112

Bullmore, E. and O.S., 2009. Complex brain networks: graph theoretical analysis of structural and functional systems. Nat Rev Neurosci. 10, 312. doi:10.1038/nrn2618

Burns, S.P., Santaniello, S., Yaffe, R.B., Jouny, C.C., Crone, N.E., Bergey, G.K., Anderson, W.S., Sarma, S. V., 2014. Network dynamics of the brain and influence of the epileptic seizure onset zone. Proc. Natl. Acad. Sci. 111, E5321–E5330. doi:10.1073/pnas.1401752111

Conrad, E.C., Bernabei, J.M., Kini, L.G., Shah, P., Mikhail, F., Kheder, A., Shinohara, R.T., Davis, K.A., Bassett, D.S., Litt, B., 2020. The sensitivity of network statistics to incomplete electrode sampling on intracranial eeg. Netw. Neurosci. 4, 484–506. doi:10.1162/netn_a_00131

De Domenico, M., Granell, C., Porter, M.A., Arenas, A., 2016. The physics of spreading processes in multilayer networks. Nat. Phys. 12, 901–906. doi:10.1038/nphys3865

Engel, J., da Silva, F.L., 2012. High-frequency oscillations – Where we are and where we need to go. Prog. Neurobiol. 98, 316–318. doi:10.1016/j.pneurobio.2012.02.001

Grinenko, O., Li, J., Mosher, J.C., Wang, I.Z., Bulacio, J.C., Gonzalez-Martinez, J., Nair, D., Najm, I., Leahy, R.M., Chauvel, P., 2018. A fingerprint of the epileptogenic zone in human epilepsies. Brain 141, 117–131. doi:10.1093/brain/awx306

Halkidi, M., 2001. On Clustering Validation Techniques – Springer 107–145.

Jacobs, J., Staba, R., Asano, E., Otsubo, H., Wu, J.Y., Zijlmans, M., Mohamed, I., Kahane, P., Dubeau, F., Navarro, V., Gotman, J., 2012. High-frequency oscillations (HFOs) in clinical epilepsy. Prog. Neurobiol. 98, 302–315. doi:10.1016/j.pneurobio.2012.03.001

Kahane, P., Landré, E., Minotti, L., Francione, S., Ryvlin, P., 2006. The Bancaud and Talairach view on the epileptogenic zone: a working hypothesis. Epileptic Disord. 8, 16–26.

Mucha, P.J., Richardson, T., Macon, K., Porter, M.A., Onnela, J.P., 2010. Community structure in time-dependent, multiscale, and multiplex networks. Science (80-. ). 328, 876–878. doi:10.1126/science.1184819

Nissen, I.A., van Klink, N.E.C., Zijlmans, M., Stam, C.J., Hillebrand, A., 2016. Brain areas with epileptic high frequency oscillations are functionally isolated in MEG virtual electrode networks. Clin. Neurophysiol. 127, 2581–2591. doi:10.1016/j.clinph.2016.04.013

Prichard, D., Theiler, J., 1994. Generating surrogate data for time series with several simultaneously measured variables. Phys. Rev. Lett. 73, 951.

Ray, S., Maunsell, J.H.R., 2011. Different origins of γ rhythm and high-γ activity in macaque visual cortex. PLoS Biol. 9. doi:10.1371/journal.pbio.1000610

Rubinov, M., Sporns, O., 2010. Complex network measures of brain connectivity: Uses and interpretations. Neuroimage 52, 1059–1069. doi:10.1016/j.neuroimage.2009.10.003

Schevon, C.A., Weiss, S.A., McKhann, G., Goodman, R.R., Yuste, R., Emerson, R.G., Trevelyan, A.J., 2012. Evidence of an inhibitory restraint of seizure activity in humans. Nat. Commun. 3, 1011–1060. doi:10.1038/ncomms2056

Smith, E.H., Liou, J.Y., Davis, T.S., Merricks, E.M., Kellis, S.S., Weiss, S.A., Greger, B., House, P.A., McKhann, G.M., Goodman, R.R., Emerson, R.G., Bateman, L.M., Trevelyan, A.J., Schevon, C.A., 2016. The ictal wavefront is the spatiotemporal source of discharges during spontaneous human seizures. Nat. Commun. 7, 1–12. doi:10.1038/ncomms11098

Tadel, F., Baillet, S., Mosher, J.C., Pantazis, D., Leahy, R.M., 2011. Brainstorm: A user-friendly application for MEG/EEG analysis. Comput. Intell. Neurosci. 2011, 1–13. doi:10.1155/2011/879716

Talairach, J., Bancaud, J., 1966. Lesion, “irritative” zone and epileptogenic focus. Confin. Neurol. 27, 91–94. doi:10.1159/000103937

Van Diessen, E., Hanemaaijer, J.I., Otte, W.M., Zelmann, R., Jacobs, J., Jansen, F.E., Dubeau, F., Stam, C.J., Gotman, J., Zijlmans, M., 2013. Are high frequency oscillations associated with altered network topology in partial epilepsy? Neuroimage 82, 564–573. doi:10.1016/j.neuroimage.2013.06.031

Wyllie, E., Gidal, B.E., Goodkin, H.P., Loddenkemper, T., Sirven, J.I., 2015. Wyllie’s treatment of epilepsy: principles and practice, 6th ed. Lippincott Williams & Wilkins.

Xu, Y., Yan, D., 2006. The Bedrosian identity for the Hilbert transform of product functions. Proc. Am. Math. Soc. 134, 2719–2728.

Zalesky, A., Fornito, A., Bullmore, E., 2012. On the use of correlation as a measure of network connectivity. Neuroimage 60, 2096–2106. doi:10.1016/j.neuroimage.2012.02.001

Zweiphenning, W.J.E.M., van ‘t Klooster, M.A., van Diessen, E., van Klink, N.E.C., Huiskamp, G.J.M., Gebbink, T.A., Leijten, F.S.S., Gosselaar, P.H., Otte, W.M., Stam, C.J., Braun, K.P.J., Zijlmans, G.J.M., 2016. High frequency oscillations and high frequency functional network characteristics in the intraoperative electrocorticogram in epilepsy. NeuroImage Clin. 12, 928–939. doi:10.1016/j.nicl.2016.09.014

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below. Please address these comments comprehensively in your response.

1. In their response and amendments to the main text, the authors claim that their method does not require cross-validation and testing over independent datasets since the approach does not involve training the model parameters across patients and does not access information regarding the resection zone. While it is true that the method is applied to each patient, separately, the methodological choices including, but not limited to, functional connectivity metrics, frequency band selection, percentile threshold selection, still effectively "train" the general algorithm and may be highly specific to the cohort studied here. This issue presents itself as a driver of inter-patient variability, as demonstrated by the conflicting findings between the mlEVC method and Fingerprint method. While this limitation detracts from the methodological advances put forth by this study, but some caution in making the claim that cross-validation/testing is not necessary is warranted. Specifically, the manuscript should include a rationale for why cross-validation was not applied in this context, coupled with an acknowledgement that such an approach in future will be important to verify the current results.

We thank reviewers for their due diligence in examining the external validity of our algorithm.

Regarding the methodological choices, we would like to emphasize that the rationale for selecting the "connectivity metric" (lagged Coherence) was mentioned in the paper. Coherence is one of the most utilized connectivity metrics in neuroscience. It is also used in the article that we were inspired by (Burns et al., 2014) when performing this research. However, Coherence is prone to spurious connectivity due to volume conduction. Although, in comparison to scalp EEG, volume conduction is less of an issue in intracranial recordings, we were still interested in avoiding it in our study. As a result, we selected lagged Coherence in our analysis. The elementwise normalization of connectivity matrices with respect to the pre-ictal period was also followed by (Burns et al., 2014).

For "frequency band selection", the main question of our study was to investigate the brain connectivity in HFOs, as these oscillations are widely analyzed in epilepsy and are considered a biomarker of the epileptogenic zone. Since the sampling rate of the data was 500 Hz for several patients, we decided to choose the 80-200 Hz range to be consistent among patients. We then divided this band in half to have narrower frequency bands. Although we noticed that in some seizures or patients, only one band (80-140 Hz or 140-200 Hz) was informative, we did not drop any frequency band to avoid adding complexity to our algorithm. A future study can design an algorithm for tailoring the proper frequency band for each patient (Arnulfo et al., 2020).

Regarding our prediction and Fingerprint algorithm, we emphasized that the two methods use different features and can be complementary. Both algorithms have a small number of false positives, and the main differences are in electrodes predicted as true positives, which in both cases they were a subset of the electrodes in the resected area. Since both studies are performed retrospectively, the actual EZ is unknown. One probability can be that the union of the predicted regions by two algorithms is part of the actual EZ and its resection was necessary for seizure freedom among those patients.

Finally, we carefully edited our previous response and added more explanations to address the reviewers' concerns. These descriptions are written under a new section in the discussion titled as "limitations and considerations", which reads as follows:

"We should point out that cross-validation was not performed in this study for several reasons. First, resection labels were not used as part of the prediction algorithm, and the same performance function and ranges of coupling and thresholding values were applied to all subjects. Second, as shown in this manuscript and similar works, parameters such as coupling or active frequency band are patient-specific. By using cross-validation, we would have to fix these parameters based on a subset of patients, which may not be ideal as the actual range of these parameters could differ from the training dataset. While this study was primarily conducted to propose the multilayer network methodology during seizures, a more rigorous investigation is necessary to evaluate the external validity of our algorithm. Future research could use data from different clinical centers to assess the generalizability of the proposed technique and identify optimal parameters."

2. The relationship between the coupling parameter and network time-scale should be further elaborated, given the importance of this parameter to the authors' claim that patients may have different rates of network change. Please present data examining how sensitive the findings of the EZ are to the choice of this parameter. Please also provide examples of how different coupling parameters would yield more insight into the different time-scales of network change

We will address both comments at the same time. In summary, the reviewers requested that we quantify the performance of EZ identification if we used a single-layer model for all patients. Additionally, they asked whether employing different coupling values will result in a better understanding of brain dynamics.

We ran the EZ identification algorithm using both a single-layer model (with coupling set to 0) and a multilayer model with different coupling values (c = 1, … 10). We then compared these results to our proposed method, which involves analyzing a range of coupling parameters and weighting them based on clustering performance. Figure 2-source data 3 presents three performance measures for each case. These include the number of electrodes predicted as true positive, the number of electrodes predicted as false positive, and the number of patients with zero true positives. Our results demonstrate that the proposed approach yields superior EZ identification performance when compared to both the single-layer and multilayer models with fixed coupling values across all patients.

In addition to improving EZ identification, using the appropriate coupling parameter for each patient can also provide valuable insights into seizure dynamics. Author response image 4 compares the brain state dynamics during three seizures in two conditions: using c = 1 and the best-matched coupling parameter (c = 10 for this patient), as discussed in the main manuscript. With c = 1, there is no clear distinction between brain states, making it challenging to observe transitions between them. In contrast, when clustering is performed using the proper coupling parameter, the brain states become separable, allowing us to distinguish between neural state transitions. Although seizure semiology information was not available in our dataset, exploring its correlation with brain state transitions can reveal further insights into brain dynamics in future studies. Author response image 5 displays the same plot for patient 3.

Author response image 4

Download asset Open asset

Author response image 5

Download asset Open asset

3. The findings comparing a single-layer model to the multi-layer model should be explicitly quantified to concretely justify the claim in the Discussion section that single-layer models yielded poorer EZ identification than multi-layer models. How much more accurate was the multi-layer model to the single-layer model?

4. More attention should be paid in the text to the risk of contamination of high frequency synchrony by two processes: (1) Filtering of spikes; and (2) Harmonics of non-sinusoidal oscillations. These topics are addressed in the following study, which could be cited in support of the discussion.

Pitfalls of high-pass filtering for detecting epileptic oscillations: a technical note on "false" ripples. Bénar CG, Chauvière L, Bartolomei F, Wendling F. Clin Neurophysiol. 2010 Mar;121(3):301-10

5. The time frequency figures should be accompanied by a presentation of the signal time course in order to fully appreciate the presence/absence of spikes. While there indeed seems to be spikes in the preictal period (vertical lines) , they present much less energy than the high frequency activity during the seizure. This large increase is quite blurred and indistinct and it is not clear what this actually represents. The original signal time course would help understand which seizure pattern this corresponds to.A consideration of the potential effect of harmonics would also be helpful.

Author response image 6 shows the EEG signal of a single node predicted as the EZ, along with two time-frequency spectra. Two sections of this data are enlarged to better illustrate pre-ictal spikes and HFOs after seizure onset (Author response image 7).

As expected, the pre-ictal spikes display broadband activity across frequencies. In contrast, the HFOs, especially between 5-10s after onset, do not exhibit any spikes in the time domain and show a transparent "blob" in the frequency domain at 90-140 Hz. This frequency range falls within the first range we analyzed, but we did not observe significant activity in the 45-70 Hz range. These oscillations cannot be harmonics of non-sinusoidal components.

While we did not visually inspect the spectra in all seizures, our observations are consistent with the "ripple" described in the referenced article (Bénar et al., 2010).

Regarding the power comparison between HFOs and spikes, please note that we normalized these spectra for each frequency bin, as per your request. We performed this normalization to the pre-ictal period to better emphasize the changes that occur after the onset. Since there are stronger low-frequency components (such as spikes) before the onset, and more high-frequency oscillations after the seizure starts, these plots amplify the energy related to HFOs.

To address your concern, we added the following paragraph under the section of “limitations and considerations”:

“Filtering and subsequent analysis of HFOs should be performed with extra caution due to the possibility of producing false or spurious ripples. Research has shown that high-pass (> 80Hz) filtering of sharp transitions in EEG (such as artifacts and spikes) or harmonics of non-sinusoidal signals can result in such false ripples (Bénar et al., 2010). Several approaches have been suggested to investigate this phenomenon. These include visually inspecting the epileptic waveform alongside the filtered signal, using the time-frequency transform of the original signal since real and false ripples have different spectra, and employing automated spectral decomposition techniques such as Matching Pursuit (Bénar et al., 2010). In our analysis, we visually evaluated several seizures in both the time and time-frequency domains. However, it would be beneficial to use one of these mentioned algorithms in a systematic manner in future studies.”

Author response image 6

Download asset Open asset

Author response image 7

Download asset Open asset

Author details

1. Hossein Shahabi

   Signal and Image Processing Institute, University of Southern California, Los Angeles, United States

   Contribution

   Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing

   For correspondence

   shahabi.h27@gmail.com

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0003-3600-8007

2. Dileep R Nair

   Epilepsy Center, Cleveland Clinic Neurological Institute, Cleveland, United States

   Contribution

   Conceptualization, Resources, Data curation, Supervision, Funding acquisition, Investigation, Writing – original draft, Project administration, Writing – review and editing

   Competing interests

   No competing interests declared

3. Richard M Leahy

   Signal and Image Processing Institute, University of Southern California, Los Angeles, United States

   Contribution

   Conceptualization, Resources, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Methodology, Writing – original draft, Project administration, Writing – review and editing

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-7278-5471

• 901

  Page views

• 164

  Downloads

• 0

  Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.