Source-level EEG and graph theory reveal widespread functional network alterations in focal epilepsy

hypersynchronous


Introduction
The human brain is a complex system which relies on coordinated and flexible network activity involving its constituent regions to regulate and sustain physiological processes, reaching from movement (King et al. 2018) and autonomic functions (Fan et al. 2012) to advanced cognition (Shine et al. 2019).Disruption of such networks is increasingly considered an important mediator of physiological dysfunction (Uhlhaas and Singer 2006).Thus, an aberrant network model may account for symptoms occurring in the absence of observable structural brain pathology.Characteristic in this respect is cognitive impairment, such as memory issues and executive dysfunctions (Holmes 2015).
The paragon of brain network pathology is epilepsy (Kramer and Cash 2012), which encompasses syndromes that are heterogeneous in aetiology and symptomatology, but share the clinical hallmark of acute transient disruption of brain function caused by hypersynchronous neuronal activity.These physiological processes often manifest as epileptic seizures, but may cause additional widespread cognitive dysfunction also present in interictal periods (van Diessen et al. 2013).In primarily focal epilepsies, the epileptic activity putatively originates in a limited pathological area (i.e., the epileptogenic zone) and may rapidly propagate to other cortical regions.Importantly, from a network perspective, local aberrancy may impact the organisation of networks at a larger scale, for example mediated by redistribution of network hubs (Pedersen et al. 2015), potentially causing subtle symptoms which are not directly associated with the classical functions of the epileptogenic brain area (Tailby et al. 2018).Identification and monitoring of these functional network changes in the brain may provide important clinical biomarkers for treatment of the wider range of symptoms associated with epileptic disorders (Haneef and Chiang 2014;Hallett et al. 2020).
In recent years, the methodological frameworks of functional brain connectivity (Bastos and Schoffelen 2015) and graph analysis (Behrens and Sporns 2012) have been increasingly used to characterise networks in the human brain.Functional connectivity (FC) is concerned with the statistical dependence between spatially separated signals obtained from electro-and magnetoencephalography (EEG/MEG) or functional magnetic resonance imaging (fMRI), and may be estimated on the basis of features related to the signals, such as amplitude, frequency and phase (van Diessen et al. 2015;Rossini et al. 2019).
From a clinical perspective, evidence suggests that network characteristics are capable of detecting differences between pathological conditions and healthy brains (van Straaten and Stam 2013;Olde Dubbelink et al. 2014;Stam 2014;Engels et al. 2015).In epilepsy, interictal network abnormalities have been consistently reported across imaging modalities, both in terms of alterations compared to healthy peers (Quraan et al. 2013;Niso et al. 2015;Vy ´tvarova et al., 2017) and as a biomarker for cognitive status (Vlooswijk et al. 2011;Rodríguez-Cruces et al. 2020).However, the clinical use of FC and graph analysis is arguably still in a rudimentary stage (Douw et al. 2019), as methodology and findings diverge substantially across studies (Vlooswijk et al. 2011;Quraan et al. 2013;Pedersen et al. 2015), and there remain unresolved issues concerning analysis parameters, such as network threshold scheme (van den Heuvel et al. 2010;van Wijk et al. 2010;Garrison et al. 2015) and underlying connectivity measures (Horstmann et al. 2010;Niso et al. 2015).
Nevertheless, the most prominent of these issues with regard to clinical epileptology concerns the choice of functional imaging modality.Currently, investigations targeting connectivity and network analysis in epilepsy are dominated by fMRI (Vlooswijk et al. 2011;McCormick et al. 2013;Vaughan et al. 2016;Vy ´tvarova et al., 2017).However, the rapidly evolving dynamics of brain networks might be better captured by the millisecond timescale of MEG (Elshahabi et al. 2015;Niso et al. 2015) and EEG (Quraan et al. 2013;Vecchio et al. 2015).Importantly, the electrophysiological approach allows investigations into oscillatory activity in specific frequency bands (Rossini et al. 2019), and by extension, their associated cognitive functions (for a review, see Lopes da Silva 2013).Despite its arguably superior signal-to-noise ratio, MEG remains scarcely available in clinical settings worldwide, and the position of EEG in epilepsy diagnosis is pervasive.However, graph theory studies of focal epilepsy via EEG-based functional networks are few, and sample sizes are generally limited (Horstmann et al. 2010;Quraan et al. 2013;Vecchio et al. 2015).Thus, to achieve a feasible means of utilising network measures in applied neurology, EEG-based network methodology warrants development.
Important in this regard is to mitigate the impact of noise and volume conduction on scalp-recorded EEG (Brunner et al. 2016), which constitute an ubiquitous, non-specific perturbation of connectivity estimates.Methods of source reconstruction have been developed (Schoffelen and Gross 2009;Vorwerk et al. 2014;Hassan and Wendling 2018), which have been shown to reduce these effects in EEG (Besserve et al. 2011).Provided further rigorous methodological assessment and development, EEG/MEG source-space connectivity is envisioned a future central role in both diagnostics and treatment of epilepsy (van Mierlo et al. 2019).
The aim of the present study was to examine brain network organisation in a cohort of focal epilepsy patients with good seizure control and high quality of life, and age-matched healthy control subjects.Specifically, we targeted the role of global network efficiency and local hub distribution in relation to epilepsy.The investigated networks were constructed with bivariate estimates of EEG-based source-space phase-locking value (PLV), under the hypothesis of phase synchronisation; i.e., functionally connected brain regions generate signals whose phases evolve together (Rosenblum et al. 1996).To assay the potential effect of network density on our findings, we evaluated the global and local graph metrics under a range of predefined density parameters.Provided the well-functioning epilepsy patients in the study, our hypothesis was that only subtle interictal network alterations would be observed compared to non-epileptic subjects.Importantly, it was not within the scope of this study to investigate clinical factors relating to the individual patient, but rather to assess group-level network alterations in a heterogenous focal epilepsy cohort.This exploratory effort aims to contribute to the following issues: (1) corroborate previous reports of alterations in global functional network organisation in focal epilepsy; (2) aid the use and interpretation of macro level local graph metrics in terms of the ''network hub" principle, in focal epilepsy; and (3) evaluate the impact of network density threshold parameters for the consistency of group differences in a potential clinical graph analysis application.

Participants
Twenty-two patients diagnosed with uni-or bilateral focal epilepsy (FE; fourteen females, age 55.1 ± 4.3 years) were recruited from neurological outpatient clinics, in connection with routine follow-up visits.All patients were characterised by high levels of daily functioning.In addition, sixteen age-matched healthy control (HC) subjects (twelve females, age 55.9 ± 6.9 years) participated.The HC participants were recruited from the FE patients' social networks, providing controls with similar socio-economic background to the patients.All FE subjects had used the same anti-seizure medication (ASM) for at least six months prior to their participation in the study.Neither FE nor HC subjects had any history of epilepsy surgery, psychiatric disorders, developmental disorders, nor any other debilitating diseases.The patients' clinical data, including ASM therapy, aetiology, EEG/MRI pathology and focus localisation, are listed patient-wise in Table 1.Of the eight patients with MRI-verified pathology, two had mesial temporal sclerosis, whereas the remaining six had other local structural pathology.None had indications of progressive neurological disease nor tumour.All patients presented with past or present focal seizures, of whom twelve also had focal seizures with secondary generalisation.Neither FE nor HC participants were compensated for their study participation.Ethical approval for the study was granted by the Regional Ethics Committee of South-Eastern Norway.All participants provided informed written consent in accordance with the Declaration of Helsinki.

EEG acquisition and preprocessing
EEG was recorded with a BioSemi 128-channel system (sampling rate of 2048 Hz) during a resting-state period.The participant was comfortably seated in a chair, resting with his or her eyes closed, but awake, for 4 minutes.To minimise between-subject variation regarding task comprehension, instructions were given in written form on a screen, as recommended by van Diessen and colleagues (2015).
High-precision information on the spatial locations of the EEG electrodes was acquired using an IO Structure Sensor (Occipital, Inc.) scanner device for iPad (Apple, Inc.).From the 3D model, the electrodes were spatially identified by an operator using MATLAB (The MathWorks, Inc.) and FieldTrip functions (ver. 2019-01-16;Oostenveld et al. 2011).For enhanced accuracy in the manual alignment of electrode positions and head model (see below), the head shape mesh was also extracted.
The EEG data were preprocessed using EEGLAB functions (ver.2019.1;Delorme and Makeig 2004).The signals were downsampled to 512 Hz, and re-referenced to an average reference, obtained by iteratively removing noisy signals (amplitude SD larger than 25 mV) from the reference signal composite.Segments containing high-amplitude signals (amplitude larger than 150 mV in more than 25% of the signals) were removed, and consistently noisy signals were removed with the PREP Pipeline toolbox (Bigdely-Shamlo et al. 2015).Line noise (50 Hz and harmonics) was suppressed with Zapline (de Cheveigné 2020).To correct artefacts related to ocular and muscular activity, the signals were decomposed with the Second Order Blind Identification (SOBI; Belouchrani et al. 1993) algorithm, andcategorised with ICLabel (Pion-Tonachini et al. 2019).The signals were band-pass filtered between 1 and 45 Hz (EEGLAB default settings).An additional 10 seconds of mirrored data was temporarily applied at both onset and offset to avoid edge artifacts.Finally, the data were segmented into non-overlapping epochs of 4 seconds, and visually evaluated.Noise-and artefact-free epochs were included for source reconstruction.

EEG source reconstruction and functional connectivity
From each participant, a T1-weighted MRI image was acquired using either one of two scanners: a 3 T Philips Achieva or a 1.5 T Siemens Magnetom Aera.This anatomical image was segmented using the SPM12 Unified Segmentation algorithm (Ashburner and Friston 2005) and the NY Head tissue probability map (Huang et al. 2016) to create a three-layer boundary element conduction model (BEM).The surfaces were generated using the iso2mesh toolbox (Qianqian Fang and Boas 2009) for MATLAB.
The source model consisted of a homogeneous regular grid of dipoles defined in MNI space, with a uniform separation of 10 mm, and labeled according to the Automated Anatomical Labeling (AAL; Tzourio-Mazoyer et al. 2002) atlas.The final grid consisted of 1210 dipoles belonging to one of the 80 cortical areas in the AAL atlas.This grid was linearly transformed to subject-space using the MRI, obtaining an individual source model.We combined this source model, the BEM conduction model and the individual electrode positions using OpenMEEG (Gramfort et al. 2010), obtaining an individual lead field.For nine subjects (five patients and four controls; all female) we were unable to obtain the individual head shape and electrode positions, and we used a set of electrode positions based on the average of the rest of participants.As the inverse method, we used a spatial filter based on linearly constrained, minimum variance beamformers (Van Veen et al. 1997).The data used to build this spatial filter was band-pass filtered between 2 and 45 Hz by means of a FIR filter applied using 2 additional seconds of real data at each side as temporary padding.
The source-space FC was evaluated under the hypothesis of phase synchronisation (Rosenblum et al. 1996) using PLV (Lachaux et al. 1999;Bruña et al. 2018), which was estimated between pairs of 26 composite regions based on the AAL atlas (listed in Table 2; for spatial layout, see Fig. 1).This was done by calculating the root-mean-square of the PLV for all the links connecting one source position in one area with one source position in the other.Calculations were done separately for five frequency bands: theta (4-8 Hz), alpha (8-12 Hz), low-beta (12-20 Hz), high-beta (20-30 Hz) and gamma (30-45 Hz).

Network analysis
The network analysis was conducted with functions from the Brain Connectivity Toolbox (BCT; ver.2019-03-03; Rubinov and Sporns 2010) and in-house MATLAB code (available upon request).Each frequency band was analysed separately.Analyses were conducted on a range of predefined density values (Lopez-Sanz et al. 2017;Sion et al. 2020).Graph metrics were calculated on both the whole-network (global) and the node (local) levels.In terms of network analysis, the FC matrices were weighted (PLV) and undirected.
All matrices were constructed with the constraint of being fully connected across all thresholds (i.e., no node was disconnected from the main network component).Full connectedness is a formal requirement in the definition of several global graph metrics.To implement this, we first computed the minimum spanning tree (Stam et al. 2014;Tewarie et al. 2015) of the inverse FC matrix (formally, the maximum spanning tree of the FC matrix), and then added weights to the minimum spanning tree backbone incrementally by descending weight order until the density threshold was reached.The predefined density threshold (DT) range extended from 25 to 75%, with increments of 5%.
For each FC matrix, one hundred re-wired null models with preserved weight, degree and strength distributions were generated from the dense matrix (BCT: null_model_und_sign; Rubinov and Sporns 2011).These null models were processed identically to the empirical network, and the mean global network metrics calculated from them were used to normalise the corresponding metric for the empirical network.
The following graph metrics were calculated: The clustering coefficient, a quantification of the density of connections between a node's neighbours.The characteristic path length, which describes the network's tendency towards global integration and efficiency.The ratio between the normalised values of a network's clustering coefficient and characteristic path length is often referred to as the small world index, a measure of to what degree the network architecture is balanced between local segregation and global integration (Watts and Strogatz 1998;Humphries and Gurney 2008).The node strength, defined as the sum of all edges connecting a node to other nodes.Finally, the eigenvector centrality, which is a measure of a node's relative importance to the network; it takes into account the number and strength of the node's edges, and whether these edges connect the node to other central nodes.For comprehensive accounts of graph metrics and their mathematical definitions, the reader is referred to Rubinov and Sporns (2010) and Newman (2008).

Power spectra analysis
PLV has been observed to be reliable in test-retest scenarios (Colclough et al. 2016;Garcés et al. 2016).However, it has some caveats, the most important being its sensitivity to volume conduction and source leakage.Here, the former was mitigated by calculating PLV in source-space.However, it is impossible to remove the effects of source leakage.Individual differences in leakage may lead to differences in the magnitude of FC (Schoffelen and Gross 2009).The varying leakage is usually generated by differences in power in the involved regions (Muthukumaraswamy and Singh 2011).To account for this possibility, we conducted comparisons of frequency band-specific relative power spectra between groups.In the regions displaying consistent local graph metric differences, we complemented the results with post hoc comparisons of power spectra.

Overall functional connectivity
Differences in the overall level of FC can spuriously generate differences in network metrics (van den Heuvel et al. 2017).In this context, overall FC is defined as the mean PLV in the individual FC matrix, discarding intra-regional estimates, and calculated separately for each frequency band.We explored possible differences in overall FC prior to conducting the analyses based on network metrics.These did not differ between groups.

Statistical tests
Taking into account the current study's exploratory approach with relatively small samples, group differences in graph metrics on both global and local levels were tested with nonparametric permutation tests (Maris and Oostenveld 2007).For each comparison, 5000 permutations were carried out.The obtained p values were corrected for multiple comparisons across frequency bands (i.e., five comparisons) with the Benjamini-Yekutieli false discovery rate with assumed positive test correlation (FDR+) procedure (Benjamini and Yekutieli 2001;Genovese et al. 2002).FDR adjusted p values (i.e., q values) below 0.1 (Niso et al. 2015) were considered statistically significant.Uncorrected p values are reported in figures.Group differences were quantified as the corrected standardised mean difference, with the effect size estimate termed Hedges' g (Lakens 2013).Directionality of effects is consistently given such that positive values of g indicate larger values for the FE group than for the HC group.All comparisons were conducted independently for each DT level.For local graph metrics, findings were considered consistent, and thus of interest, if significant frequency bandspecific group differences were observed for a particular node in three or more consecutive DT increments.Mean q and g values were calculated across the significant DT levels of these node differences.
Group differences in overall FC and power spectra were analysed in an identical permutation procedure.In contrast to the exploratory analyses, the resulting p values from the overall FC tests were not corrected for multiple comparisons.This also applied to the relative power spectra comparison which was defined as a post hoc procedure for complementing the frequency band/node combinations which yielded significant group differences in local graph metrics.Furthermore, in additional post hoc procedures, we examined the possible effects of age, sex and education (total number of years) on global and local graph metrics where a significant effect of group was evident.These analyses of covariates were implemented as ANOVAs.Similarly, for a subgroup of the patients characterised by a clearly defined ictal hemisphere (n = 18), the potential effects of focus lateralisation, age of first seizure, and epilepsy duration on global and local graph metrics were analysed.For the local metrics, the mean measure across the significant DT levels was defined as the dependent variable.All statistical analyses were conducted with MATLAB and SPSS (IBM Corp.; ver.27).

Global graph metrics
On the global level, group differences in normalised clustering coefficient and characteristic path length, and the small world index, were analysed.In line with our hypothesis, we observed subtle, yet consistent differences between patients and controls.
The networks of both patient and control groups displayed evidence of small world architecture (index larger than 1), which was represented in all frequency bands and across the majority of DT levels (Fig. 2, top row).However, in the upper areas of the DT range, the small world index fell below 1, suggesting that higherdensity networks failed to display this global organisation.Upon comparison between the groups, the patients consistently showed higher indices, although at some DT levels, the direction of the difference was reversed (Fig. 2, bottom row).The group effects were most pronounced for the DT range between 60 and 70%, achieving statistical significance (Fig. 2, middle row) in the theta (DT = 65%: q = 0.074, g = 0.730), alpha (DT = 65%: q = 0.096, g = 0.629), highbeta (DT = 70%: q = 0.089, g = 0.706) and gamma bands (DT = 65/70%: q = 0.074/0.089,g = 0.766/0.714).
In addition, we analysed the normalised clustering coefficient (Fig. 3) and the normalised characteristic path length (Fig. 4) separately.Both patient and control groups displayed normalised val- ues above 1, indicating that the observed network metrics were higher than the corresponding metrics in the randomly re-wired networks, in line with our expectations.The patients displayed consistently higher clustering coefficient and lower characteristic path length than the controls, but the effects did not reach statistical significance after FDR adjustment.Yet, these group differences showed evidence of stability across all frequency bands and DT levels, with the differences being the smallest in the low-beta band and most pronounced in the gamma band (clustering coefficient) and high-beta band (characteristic path length).
For all assessed global metrics, the normalised values shifted closer to random with increased network density.Also, the results demonstrated a tendency towards increased effect sizes in the upper end of the DT range.
In the post hoc ANOVA procedures examining the possible effect of covariates on the global graph metrics where a significant effect of group was evident, we did not observe any significant effect of neither demographic variables nor epilepsy factors.The group comparisons were not dependent on age (p > 0.300), sex (p > 0.185) nor number of years of education (p > 0.325).In the subsample of patients defined by clear epilepsy focus lateralisation, neither age of epilepsy onset (p > 0.208), epilepsy duration (p > 0.198) nor ictal hemisphere (p > 0.233) mediated the global graph metrics.

Local graph metrics
Analyses on the local level revealed differences (significant effects in three consecutive DT increments) in several nodes across multiple frequency bands and graph metrics (see Fig. 5).Here, we focus on effects in two nodes: The left central region (LCR; node 1 in Fig. 1) and the left parietal lateral surface (LPL; node 13).These nodes displayed specific alterations in the alpha band on character-istic path length, node strength and eigenvector centrality (Fig. 6), all three characteristics implied in the evaluation of hubness (see discussion).The effects reported in the following are considered statistically significant, having survived FDR adjustment.Table 3 contains a list of all nodal group differences evident in three or more consecutive DT increments.
Similarly to the global metrics, the demographic variables (age, p > 0.398; sex, p > 0.086; education, p = 0.204) and epilepsy factors (age of epilepsy onset, p > 0.201; epilepsy duration, p = 0.303; ictal hemisphere, p > 0.362) did not show any significant effect on the local graph metrics on which a significant effect of group was evident.

Spectral power
In a post hoc approach, we analysed potential group differences in region-wise relative power at the nodes which revealed stable group effects in local graph metrics (Fig. 6, rightmost column).This was done on the basis of delineating potential confounding effects of power on graph metrics (see discussion on source leakage).The   False detection rate (FDR) adjusted q values from the mass permutation testing of group differences in local graph metrics.Significant (q < 0.1) difference in a node is flagged with a white asterisk.below reported p values are uncorrected with regard to multiple comparisons.
For the patient group, increased alpha power was observed in both LCR (p = 0.022, g = 0.761) and LPL (p = 0.045, g = 0.678) nodes.The LPL node also showed a significant decrease for the patient group in high-beta power (p = 0.002, g = -0.984).In addition, in the node representing the right occipital and medial inferior surfaces, the patient group displayed an increase in relative power (p = 0.041, g = 0.667).
The full analysis of all regions across all frequency bands is presented in Fig. 7.

Discussion
To our knowledge, the current study is the first to present a comprehensive analysis of both global and local graph metrics in adult focal epilepsy via FC estimated from source-space scalp EEG.In recent years, analysis of functional brain networks has shown clinical potential in the diagnostics and treatment of epilepsy.Here, we demonstrate that changes to functional network topology on both the global and local levels are evident even in well-functioning patients, suggesting that altered network organisation is a stable trait of focal epilepsy.In our results, epilepsy patients showed consistently higher small world indices, and the investigation of individual nodes revealed evidence of altered hubness of the left central region and the left parietal lateral surface.Interestingly, in our results, these network alterations did not appear mediated by neither demographic variables nor clinical epilepsy factors related to ictal lateralisation, age of epilepsy onset or epilepsy duration.

Global metrics: Increased small worldness in focal epilepsy
Network-wide changes in epilepsy are consistently reported, yet their clinical significance remains inadequately understood.It is generally accepted that epileptic symptoms arise from abnormal neuronal synchronisation and excitability (Engel et al. 2013).Although these phenomena can be observed experimentally through FC, the understanding of their underlying mechanisms remains scarce (van Diessen et al. 2013).Graph theoretical characteristics of topological network organisation could remedy this shortcoming (Bernhardt et al. 2015).However, the links between graph metrics and the pathophysiological mechanisms of epilepsy have yet to be firmly established.
Here, we observed an increase in the small world index of the theta, alpha, high-beta and gamma bands in patients relative to controls.These observations partly corroborate previous reports, where focal epilepsy patients have displayed increased small world index in the theta and alpha bands (Quraan et al. 2013).Interestingly, the shift towards increased small worldness contrasts earlier findings suggesting that epileptic activity is associated with increased regularisation of brain networks, manifested as a dual increase in clustering coefficient and characteristic path length (Horstmann et al. 2010;Vecchio et al. 2015).The latter tendency was also reported in an across-modality meta-analysis (van Diessen et al. 2014), and in an fMRI study by Pedersen and colleagues (2015), who suggested that epileptic network activity may arise in accordance with the fault tolerant network (Johnson 1984) hypothesis.In this scenario, a network with the potential of evolving dysfunctional nodes will compensate by altering its topology towards increased regularity.This shift in favour of network regularisation, however, was not evident in our results, which in contrast demonstrated a decrease in characteristic path length.The reason for this discrepancy is not clear; however, both the referenced EEG studies differed methodologically from ours, in where one (Vecchio et al. 2015) focused exclusively on fronto-temporal areas in dense weighted network matrices, whereas the other (Horstmann et al. 2010) employed sensorspace connectivity estimates, thus more affected by volume conduction confounds.Conversely, it is possible that heterogeneity with regards to individual clinical factors in the current patient sample gave rise to the current discrepancy.Notably, in this vein, another study examined the temporal lobe specifically with intracranial EEG, and found that small worldness in focal epilepsy could depend on the time since initial diagnosis (van Dellen et al. 2009), an association that was not examined here.
However, the relevance of the small world index remains an issue of debate (Papo et al. 2016).As a summary metric, it captures the network topology's trade-off between local segregation and global integration, i.e., to what degree the brain retains its modular functional organisation with a high number of short communication paths, yet facilitates signal transmission between modules via long-range paths to achieve integrated functioning (Bertolero Table 3 Consistent (across three consecutive density threshold levels) group differences in local graph metrics.The metrics are characteristic path length (CPL), clustering coefficient (CC), node strength (NS) and eigenvector centrality (EvC).The values in the column ''Mean metric g" are calculated across the significant density threshold (DT) levels (positive values indicate larger values for the patient group than for the control group).The ''Relative power g (p)" column provides the frequency band/region group comparison statistics (uncorrected p values; significant differences at p < 0.05 level in bold).

Region label
Frequency  2015).Mathematically, the index can increase either as a function of an increased clustering coefficient, a decreased characteristic path length, or both; however, it is unclear in what range the brain can reside and still retain its optimal functioning.From our data, an association between hypersynchronous neuronal activity and increased small worldness in epilepsy is likely.Hypersynchronous neuronal activity is classically understood as a ''problem of regions and neurons connected or communicating too readily" (Kramer and Cash 2012).More importantly, in focal epilepsies without secondary generalisation, such hypersynchronous activity originates in a delimited brain region with propagation restricted to other regions with high (functional) proximity to the epileptogenic zone.Thus, high FC between neighbouring regions will be likely to occur, and the proportional threshold scheme will result in the rejection of more long-range connections in the patient group relative to the controls.This notion is supported further by evidence that global efficiency (inverse characteristic path length) is higher in patients with generalised epilepsy than focal epilepsy, which again is higher than in healthy controls (Niso et al. 2015).

Local metrics: Altered network hubs
Node-level graph metrics are considerably less discussed in the literature than global metrics, and the interpretational framework for local network alterations remains limited.However, the concept of hubness (i.e., the degree to which a node constitutes a hub in the network) has gained some influence (van den Heuvel et al. 2010;Ridley et al. 2015;Lee et al. 2018).A hub is a node of relatively high importance to the network, being essential to efficient integration of information across the network.Consequently, local disruptions of hub nodes are more likely to cause severe impairment to global integrative processes and network organisation, due to their central positions (van den Heuvel and Sporns 2013), than disruptions to regular nodes.In terms of graph metrics, a hub node may be characterised by high node degree/strength, high centrality, and short characteristic path length, combined with low clustering coefficient (van den Heuvel et al. 2010).
The role of (changing) brain network hubs has yet to be discerned in the context of focal epilepsy.Provided that brain hubs are regions with a strong degree of participation in information transmission throughout the network, they must necessarily also be considered at risk to act as bottlenecks in the system (Marois and Ivanoff 2005;van den Heuvel and Sporns 2013).Pursuing this notion, increased hubness in a selection of nodes in epilepsy might imply that epileptic activity reduces the efficiency of brain networks by excessively routing information through these nodes, maladaptively generating communication bottlenecks.Importantly, in this scenario, a dysfunction in information flow would be an effect of epileptic neuronal activity, and a possible mediator for cognitive dysfunction and other non-seizure symptoms of epilepsy.In theory, such nodes could serve as potential targets for surgical intervention; however, their anatomical localisation might render this utility impractical.
Here, we demonstrate a shift towards increased hubness in the left central region and the left parietal lateral surface in the alpha band of focal epilepsy patients.This shift manifests as nodal changes characterised by shorter characteristic path length and increased node strength and eigenvector centrality.Interestingly, the parietal lateral surface, as defined here, includes the angular gyrus, a region in which local network alterations have been reported previously (Ridley et al. 2015).Furthermore, the angular gyrus is recognised as (part of) a functional hub of the default mode network, implicated in introspective and social cognitive processes (Andrews-Hanna et al. 2014), and dysfunctions of the default mode network have been demonstrated interictally in focal epilepsy (Burianová et al. 2017;Lee et al. 2018;Ofer et al. 2019).To our knowledge, the current study is the first to do so using sourcespace analysis based on EEG, thus advancing the recognition of specific dysfunctional networks and brain regions.Moreover, temporal lobe epilepsies are associated with dysfunctional theory of mind and facial emotion recognition (Bora and Meletti 2016), which are social cognition processes dependent on the default mode network.

Network density threshold effect
A subsidiary aim of the current effort was to probe the consistency of group differences in graph metrics across a predefined range of adjacency matrix densities, a matter which remains unresolved in the field.Functional networks are per definition arbitrary, and no basis exists on which to define an ecologically valid threshold level.The motivation to impose a threshold on the FC matrix prior to computing graph metrics, comes from the assumption that weak connections are more likely to be spurious (and thus conceal key topological properties) than strong connections (Fornito et al. 2012; for a discussion of PLV, see also Aydore et al. 2013).
Mathematically, network density directly influences a range of graph metrics (for a detailed account, see van Wijk et al. 2010).The primary concerns relating to this issue are the reproducibility of findings across studies and the comparison of graph metrics across groups within the same study.The latter issue is mostly solved by the application of proportional (i.e., density-based) thresholds, where the connection weight cut-off varies individually in order to standardise the network density across subjects.Here, we corroborate previous reports that the proportional threshold scheme produces relatively consistent group differences across the threshold continuum (Quraan et al. 2013;van den Heuvel et al. 2017).However, in contrast to absolute thresholds, the use of proportional thresholds does not exclude the possibility that differences in overall FC of the dense matrices influence the resulting graph metrics (van den Heuvel et al. 2017).Our results did not reveal altered overall FC between groups, suggesting that the topological group effects were not inflated by the network threshold scheme.
We investigated group differences across a middle range of density threshold levels, thus excluding the extremes, which have been shown to produce erroneous results in fMRI-based connectivity (Garrison et al. 2015).In accordance with previous reports, we found relatively consistent directionality of effects across density levels (i.e., there were few ''sign reversals"), when a proportional threshold scheme was used (Garrison et al. 2015;van den Heuvel et al. 2017).It should be noted that in the current approach, in contrast to the referenced studies, network matrices were constrained by the requirement to remain fully connected after threshold imposition.

Clinical relevance of network analysis in epilepsy and future aims
Epilepsy is increasingly considered a disease of brain network pathology, and procedures utilising FC and graph theory have been proposed for several clinical applications relating to epilepsy diagnostics and treatment (Stefan and Lopes da Silva 2013;van Mierlo et al. 2019).Among the most investigated are potential improvements to pre-surgical evaluation, such as localisation of epileptogenic zones (Panzica et al. 2013;van Mierlo et al. 2014), and prediction of resection outcome.An intriguing example of the latter, is a recently introduced method of virtual resection of functional brain networks derived from intracranial EEG (Kini et al. 2019).The authors report high precision in predicting which patients with drug-resistant focal epilepsy will benefit from epilepsy surgery, as well as in restricting the resection zone.In terms of non-invasive methods to improve epilepsy surgery, FC from source EEG shows promise in identification of surgical targets (Staljanssens et al. 2017), and graph metrics (in this case derived from fMRI) have been successfully employed to predict the cognitive outcome for patients undergoing temporal lobe resections (Doucet et al. 2015).The association between cognitive functioning and network metrics has also been demonstrated with EEG methods, in cohorts of Alzheimer's disease and mild cognitive impairment patients (Vecchio et al. 2016), and in healthy subjects (Langer et al. 2012).Here, we observed robust network alterations in theta, alpha, high-beta and gamma bands; all oscillation frequencies associated with different cognitive functions.Both slowwave delta and theta, and gamma frequencies have been consistently implicated in memory processes (Raghavachari et al. 2006;Axmacher et al. 2008;Sauseng et al. 2009;Jacobs and Kahana 2010;Hanslmayr et al. 2019), whereas alpha rhythm is associated with control of attention (Mathewson et al. 2009); both of which are affected in focal epilepsy (Hermann et al. 2007).However, whether aberrant network topology in certain oscillatory bands is related to dysfunction in specific cognitive domains, remains unclear, and should constitute an important target for future studies.
In sum, it is reasonable to believe that EEG-based functional network analysis in the near future will play an integral part in diagnostics and treatment of epilepsy.As demonstrated here, although the patient cohort was characterised by successful ASM therapy and high quality of life, functional network alterations were still present relative to non-epileptic peers.Thus, an aim for future studies should be to investigate specific topographical alterations in relation to customary elements of epilepsy diagnostics and treatment, such as cognitive dysfunction and ASM.Provided the small sample investigated in the current study, we did not stratify the patients with regard to medication.However, some efforts have been made towards examining the potential of network analysis as an indicator for successful treatment, although clinical validation remains sparse.One study found that carbamazepine and oxcarbazepine affected the betweenness centrality of several hubs, suggesting that these medications work, at least in part, by altering (or reverting the alterations introduced by the disease) the brain's hub organisation (Haneef et al. 2015).Another study found selective differences in clustering coefficient and global efficiency in patients using topiramate compared to other medications (van Veenendaal et al. 2017).Taken together, these studies suggest that ASMs may have differential effects on functional brain networks.If this holds true, network topology might guide the choice of ASM administered to newly diagnosed epilepsies, an important goal for clinical neurology.Thus, comprehensive studies of ASMs' effect on functional network topology are highly warranted, such as the comparison between larger groups of patients receiving different ASMs in monotherapy, also including patients tapered off medication.Furthermore, targeting other patient groups using ASMs, e.g.psychiatric or neuropathic pain, could further contribute to the delineation of ASMs' effect on functional brain network organisation.

Limitations
Some limitations with the current study should be addressed.First, our sample sizes were small.Although we were able to identify significant group differences, the low number of participants did not allow precise stratification with regard to individual clinical factors, such as ASM therapy, epileptogenic zone localisation or lateralisation, seizure status and structural brain lesions.It remains possible that one or more of these factors play important mediating roles in the functional network disruptions observed in focal epilepsy.The current sample size also limits the generalisability of our statistical approach; thus, we encourage future studies to corroborate the current findings in larger samples.Second, compared to fMRI and MEG, EEG has limitations regarding spatial resolution.We have mitigated this shortcoming by employing highquality source reconstruction, and by using a reduced version of the AAL atlas with larger regions.On the other hand, the low cost and clinical availability associated with EEG justify this trade-off.
Third, electrophysiologically derived FC is limited by the effect of source leakage.If the leakage varies across groups, differences can arise in the FC matrices, and by extension, the graph metrics.While some FC measures eliminate the effects of leakage by removing zero-lag synchronisation (Ewald et al. 2012;Bruña et al. 2018), these show low test-retest reliability (Colclough et al. 2016;Garcés et al. 2016).Here, we decided to quantify source leakage, using local power as a proxy (Muthukumaraswamy and Singh 2011).In our data, local power in the alpha band was increased in the patients in both LCR and LPL nodes.This could result from higher leakage, increasing the FC with neighbouring areas; however, the associated effect sizes were smaller than those of the graph metrics, arguing that the differences in power were less relevant.In further support of this notion, the remaining nodes displaying consistent alterations did not show significant differences in power.Overall, the power differences between groups are most pronounced in the high-beta band (Fig. 7), in which we did not observe changes in node hubness.

Conclusion
Although we investigated a relatively small sample, we have presented evidence that focal epilepsy patients, despite good seizure control and high quality of life, exhibit widespread functional network alterations relative to healthy peers.Interestingly, these discrepancies were evident at the group-level, suggesting that functional network alterations are evident in focal epilepsy regardless of (or as the sum of) individual clinical factors.Importantly, the present findings add to previous accounts to provide support for the potential clinical relevance of network analysis as a framework in which secondary epilepsy symptoms, such as cognitive dysfunction, may be understood.Moreover, this relevance might extend to several clinical applications, ranging from pre-surgical localisation of epilepsy foci and post-surgery outcome prediction, to monitoring of patients over time with a view to map ASM effectiveness and side-effect profile, and cognitive health.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Node/brain region mapping.Node indices correspond to region indices listed in Table2.Left hemisphere nodes are represented by red dots, whereas right hemisphere nodes are blue.The node positions are calculated from the centroid of the regions comprising the node.Top: Axial view; bottom: medial view.

Fig. 2 .
Fig. 2. Small world index.Group comparison of the small world index across density threshold levels (FE = focal epilepsy patients; HC = healthy controls).The grey shadings indicate significant group difference at given density threshold level.Top row: Small world index, where the vertical lines represent the 95% bias-corrected and accelerated confidence interval of the mean.Middle row: Uncorrected p values and false detection rate (FDR) adjusted q values associated with the permutation tests.The horizontal dashed line represents the critical alpha threshold (q < 0.1).Bottom row: Effect size estimates.

Fig. 4 .
Fig. 4. Normalised characteristic path length.Group comparison of the network-wide normalised characteristic path length across density threshold levels (FE = focal epilepsy patients; HC = healthy controls).Top row: Normalised characteristic path length, where the vertical lines represent the 95% bias-corrected and accelerated confidence interval of the mean.Middle row: Uncorrected p values and false detection rate (FDR) adjusted q values associated with the permutation tests.The horizontal dashed line represents the critical alpha threshold (q < 0.1).Bottom row: Effect size estimates.

Fig. 5 .
Fig. 5. Mass testing of local graph metrics.False detection rate (FDR) adjusted q values from the mass permutation testing of group differences in local graph metrics.Significant (q < 0.1) difference in a node is flagged with a white asterisk.

Fig. 6 .
Fig. 6.Node hubness.Comparison of local graph metrics and relative power in the left central region (A) and left parietal lateral surface (B) nodes (FE = focal epilepsy patients; HC = healthy controls).Columns 1-3: Shaded grey areas indicate significant difference at given density threshold level (q < 0.10).Column 4: Double asterisk indicates significant group difference when both uncorrected (p < 0.05) and false detection rate (FDR) adjusted (q < 0.10).

Fig. 7 .
Fig. 7. Spectral power.Region-wise group comparison of relative power (FE = focal epilepsy patients; HC = healthy controls).A black rectangle indicates significance at q < 0.10 level after false detection rate (FDR) adjustment of p values from the permutation test.

Table 1
Clinical data of the patients.Age given in years.Sex: F = female; M = male.Duration of epilepsy is defined as years since the patient's first seizure.

Table 2
List of composite regions based on areas defined in the Automated Anatomical Labeling (AAL) atlas.The AAL area indices correspond to the original area/index pairing in Tzourio-Mazoyer et al. (2002).The region indices correspond to the map in Fig. 1.