Using network analysis to localize the epileptogenic zone from invasive EEG recordings in intractable focal epilepsy

Treatment of medically intractable focal epilepsy (MIFE) by surgical resection of the epileptogenic zone (EZ) is often effective provided the EZ can be reliably identified. Even with the use of invasive recordings, the clinical differentiation between the EZ and normal brain areas can be quite challenging, mainly in patients without MRI detectable lesions. Consequently, despite relatively large brain regions being removed, surgical success rates barely reach 60–65%. Such variable and unfavorable outcomes associated with high morbidity rates are often caused by imprecise and/or inaccurate EZ localization. We developed a localization algorithm that uses network-based data analytics to process invasive EEG recordings. This network algorithm analyzes the centrality signatures of every contact electrode within the recording network and characterizes contacts into susceptible EZ based on the centrality trends over time. The algorithm was tested in a retrospective study that included 42 patients from four epilepsy centers. Our algorithm had higher agreement with EZ regions identified by clinicians for patients with successful surgical outcomes and less agreement for patients with failed outcomes. These findings suggest that network analytics and a network systems perspective of epilepsy may be useful in assisting clinicians in more accurately localizing the EZ.


INTRODUCTION
Epilepsy is one of the most common brain disorders, characterized by chronically recurrent seizures resulting from excessive electrical discharges from groups of neurons (Brodie et al., 1997). Epilepsy affects over 50 million people worldwide, and over 30% of all individuals with epilepsy have intractable seizures, which cannot completely be controlled by medical therapy (Berg, 2009;Berg & Kelly, 2006;Kwan & Brodie, 2000). That is, seizures continue to occur despite treatment with a maximally tolerated dose of at least two antiepilepsy drugs (AEDs). The direct cost of assessing and treating patients with medically intractable focal epilepsy (MIFE) Focal epilepsy: A subset of epilepsy in which seizures originate from a focal area of the brain. ranges from $3 to $4 billion annually ($16 billion in direct and indirect costs) in the United States (Murray, Halpern, & Leppik, 1996). Eighty percent of these costs are incurred by patients whose seizures are not adequately controlled by AEDs (Begley et al., 2000). The burden of MIFE, however, is much greater than heavy financial costs. MIFE is a debilitating illness where individuals lose their independence, causing profound behavioral, psychological, social, financial, and legal issues (Ferro & Speechley, 2009;F. Gilliam et al., 1999;F. G. Gilliam, 2005;Hermann et al., 2006;Schuele & Lüders, 2008). Cognitive performance may be impaired by MIFE as well as by side effects of AED therapy (Ferro & Speechley, 2009;F. Gilliam et al., 1999;F. G. Gilliam, 2005;Hermann et al., 2006;Schuele & Lüders, 2008).
Despite the heavy sequelae from MIFE, there is a potentially curative procedure, surgical resection of the epileptogenic zone (EZ), which can be defined as the minimal area of brain tissue Epileptogenic zone: The minimal amount of brain region that must be removed to stop seizures from occurring. responsible for generating the recurrent seizure activity (Lüders, Najm, Nair, Widdess-Walsh, & Bingman, 2006). However, to be effective, this procedure depends on correct anatomical identification of the EZ, which is often poorly defined. A comprehensive presurgical evaluation is necessary to better delineate the EZ as well as to identify the risk of neurologic morbidity such as motor, visual, or speech impairment. Various noninvasive methods are currently applied in the attempt of defining the EZ, the eloquent cortical and subcortical areas and, consequently, the optimal resective surgical strategy. Noninvasive techniques include scalp EEG and video-EEG monitoring, neuropsychological tests, speech-language studies, and brain imaging (MRI, PET, ictal SPECT). Of these methods, the highest predictor of surgical success is identification of a single visible MRI lesion Jeha et al., 2007Jeha et al., , 2006McIntosh et al., 2004;See et al., 2013;Urbach et al., 2004).
Localization and surgical success in seizure control are even more challenging in patients with nonlesional MRI. When the noninvasive methods of localization fail to identify the EZ, an invasive monitoring evaluation may be indicated, involving the implantation of subdural grid electrodes (SDE) through open craniotomies or stereoelectroencephalography (SEEG; Nair, Burgess, McIntyre, & Lüders, 2008;Önal et al., 2003;Widdess-Walsh et al., 2007). The process of identifying the EZ then involves visually inspecting tens to hundreds of invasive EEG signals without much assistance from computational tools. Epileptologists currently study the onset of seizure events that occur over several days. Early presence of beta-band activity (beta buzz) or bursts of high-frequency oscillations (HFOs) in the 100-300 Hz range, which typically occur High-frequency oscillations: Electrophysiological phenomena that involve oscillatory activity at very high frequencies (such as greater than 200 Hz). milliseconds before the clinical onset of seizures are localizing of the seizure onset (Fisher, 2012). Channels where seizure onset features first appear are commonly defined as the seizure onset zone (SOZ), the current best estimate of the unknown EZ. This is based on the assumption that the epileptic cortex generates epileptiform activity, which then entrains other regions into a clinical seizure (Fisher, 2012). Electrodecremental responses (loss of rhythmic activity) are also often observed. In general, epileptologists look at a variety of signatures to make their decision (Fisher, 2012). Despite all of these possible EEG signatures, determination of the EZ may remain unclear for nonlesional patients (Gonzalez-Martinez et al., 2013;Jung, Pacia, & Devinsky, 1999;Niedermeyer & Silva, 2004;Wieser, 1998). See Figure 1 for a schematic of a current clinical process of localizing the EZ.
Network analysis of intracranial EEG data has been heavily used to study brain activity (Bassett & Bullmore, 2006;Braun, Muldoon, & Bassett, 2015;Bullmore & Bassett, 2011;Deuker et al., 2009). Networked-based analysis assumes that signals from different EEG channels are samples of activity from brain regions that are structurally and/or functionally connected and therefore dependent (Kerr et al., 2011;Santaniello et al., 2011;Yaffe et al., 2012). Several important prior studies have looked at network dynamics in epileptic cortex during seizure events. Some works investigate correlation structure over seizure events and note changes in network coherence over events without relating metrics back to clinically annotated EZ (Kramer et al., 2010;Schindler, Leung, Elger, & Lehnertz, 2007). Other studies apply network methods, computing interelectrode coherence, and relate these measures back to clinically annotated EZ or resection regions, but on data collected from a relatively small set of patients (Khambhati, Davis, Lucas, Litt, & Bassett, 2016;Korzeniewska et al., 2014;Schevon et al., 2007;Sinha et al., 2017). Studies that incorporate computational modeling to explain mechanisms of seizures and the EZ include Khambhati et al. (2016) and Sinha et al. (2017) .
Here, we show a novel network-based algorithm that takes advantage of a certain type of signal evolution (ranked eigenvector centrality) and utilizes preictal, ictal, and postictal data Eigenvector centrality: is a graph theoretic measure of how influential a node is within a network.
for tissue suspected to be within the EZ. Our study combines data from four centers and analyzes a total of 113 seizures from 42 patients. We compute network-based statistics and Figure 1. Clinical process for implantation of SDE and seizure onset localization. Clinicians expose the brain through a craniotomy, then implant electrodes on the cortical surface of the brain, monitor patient electrocorticography (ECoG) for days/weeks, and then attempt to localize the EZ visually. Clinical teams look at recorded data on computers and annotate signals from certain electrodes and time periods. relate the eigenvector centrality (EVC) patterns back to clinically annotated EZ in patients with both successful and failed outcomes. We recently demonstrated that intracranial EEG (iEEG) is rich in network information beyond the typical signatures clinicians use to identify the EZ (Burns et al., 2014;Kerr et al., 2011;Santaniello et al., 2011;Yaffe et al., 2012). In particular, we modeled the epileptic brain as a dynamic networked system where EEG signals are correlated both temporally and spatially. We constructed a set of network-based statistics whose temporal evolution distinguishes the epileptic nodes from the nonepileptic nodes within specific epileptic networks, thus defining an electrophysiological signature of the EZ (Kerr et al., 2011;Yaffe et al., 2012). The electrophysiological signature of the EZ has a characteristic arch shape when visualized in a two-dimensional principal component (2D PC) space Principal component analysis (PCA): A dimensionality reduction technique commonly used in data analysis. described below. The arch shape is significant because it indicates that the electrodes have lower centrality before a seizure, become highly central during a seizure, and then become less central after seizure offset. This suggests that the EZ is a brain region that becomes highly centralized when seizures occur, recruiting many other brain regions to participate in epileptic activity. We used these time series network-based statistics and the identified EZ arch signature to develop an algorithm that takes as inputs iEEG data and the patient's brain image after electrode implantation and outputs the likelihood of an electrode being in the EZ.
We hypothesized that a network-based algorithm will show higher degrees of agreement with the clinically labeled EZ for successful surgical outcomes and lower degrees of agreement with the labeled EZ for failed surgical outcomes. Our hypothesis is based on our expectations that a network-based algorithm will perform favorably because epilepsy is a network disease of the brain and simply looking at biomarkers of individual electrodes ignores this fact. To test our hypothesis, we evaluated our algorithm in a blind, retrospective study on 42 patients that had undergone invasive monitoring and in most cases were followed by surgery. EEG data on one to three seizures were analyzed by our algorithm without knowledge of the seizure outcomes. Clinically identified EZ nodes were then compared with the most central nodes as defined by our algorithm. We found that the algorithm agreed more with clinical annotations for patients with successful surgical outcomes and less for patients with failed surgical outcomes. Since HFO is considered a gold standard for localization of high-frequency power, we wanted to compare our results with such a method. We also applied qHFO algorithm presented in qHFO: A quality high-frequency oscillation (i.e., highly likely to be an actual HFO) detected by an algorithm described in the paper. Gliske et al. (2016) to all patients whose EEG recordings met the requirements of the qHFO algorithm. We found that there were many patient datasets that could not be easily applied to the qHFO algorithm because of limitations on data available and sampling rates of equipment. However, on the datasets that could be compared with our network algorithm, there was a higher degree of agreement (DOA) with clinicians using a network algorithm versus only the qHFO algorithm.
Localization of the EZ is currently a time-consuming process since clinicians and technicians visually inspect fairly large datasets. In today's data science era, it is important to develop and test computational tools to assist in localization of the EZ. An assistive computational tool would not only likely reduce extraoperative monitoring time in the EMU, thereby cutting medical costs and decreasing complications associated with invasive monitoring, but could also improve seizure freedom rates, especially in the more difficult to localize patients (i.e., nonlesional MRI patients). In addition, the underlying network-based algorithm that performs EZ detection favorably will further our understanding of the organization and dynamics of brain networks in epilepsy disease. Our results suggest that epilepsy changes how the different nodes in the brain are connected, and that diseased nodes are more likely to be highly central in the neuronal network and have a high centrality signature.

METHODS: DATA COLLECTION
Patients included in this study were surgically treated for medically intractable seizures at four different centers: the Johns Hopkins Hospital (JHH), the National Institute of Health (NIH), the University of Maryland Medical Center (UMMC), and the Cleveland Clinic (CC). All patients included in this study underwent invasive presurgical monitoring with either subdural grid-and-strip arrays or stereotactic EEG depth electrodes for seizure localization or mapping of eloquent areas. Decisions regarding the need for invasive monitoring and the placement of electrode arrays were made independently of this work and solely based on clinical necessity. The research protocol was reviewed by the Johns Hopkins Institutional Review Board (IRB), the National Institute of Neurological Disorders and Stroke IRB, the University of Maryland Medical Center IRB, and the Cleveland Clinic IRB. The acquisition of data for research purposes was done with no impact on the clinical objectives of the patient stay. Digitized data were stored in an IRB-approved database compliant with Health Insurance Portability and Accountability Act (HIPAA) regulations (e.g., server hosted behind a firewall with SFTP and SSH access).
At all four centers, as part of routine clinical care, up to three board-certified epileptologists marked, by consensus, the unequivocal electrographic onset of each seizure and the period between seizure onset and termination. The seizure onset was indicated by a variety of stereotypical electrographic features, which include, but were not limited to, the onset of fast rhythmic activity, an isolated spike or spike-and-wave complex followed by rhythmic activity, or an electrodecremental response. Concurrently with the examination of the EEG recordings, changes in the patients' behavior were sought from the video segment of video-EEG recordings. For each patient, we combined surgical notes about the electrodes corresponding to resected regions and postoperative follow-up information about how the resection affected the patient's seizures. The surgery was deemed a success and the resected area determined to include the EZ if, at least six months after surgery, a patient reported no seizures or could manage their epilepsy with medications. Failure was defined as the inability to localize the EZ at all, or if the patient continued to have seizures that were not manageable with medications after the resection.
IEEG recordings were acquired through subdural grid arrays, subdural strip electrodes, or depth-electrode arrays in various combinations as determined by clinical assessment for patients with temporal, occipital, or frontal lobe seizures. Subdural grids have 20-64 contacts per array and were used in combination with subdural strips with 4-8 contacts or depth arrays, thus having 80-116 recording electrodes per patient overall. Intracranial contact locations were documented by postoperative CT coregistered with a preoperative MRI. Signals were acquired using continuous multichannel iEEG recordings collected over 5 days on average (min.: 2 days; max: 10 days). Clinical monitoring lasted 5-10 days per patient and included two to seven clinical seizures. Then clinicians clipped what they deemed clean sets of data and passed it through a secure transfer for the data analysis.
There were a total of 42 subjects analyzed retrospectively in this study: 7 from NIH, 20 from JHH, 7 from UMMC, and 8 from the Cleveland Clinic. There were 26 total successful surgeries and 16 total failed surgeries. The total number of electrodes per patient was 111.86 ± 23.89. The total number of electrodes used in analysis per patient (after removal of noisy/faulty channels, references, EKG, etc.) was 70.82 ± 24.84. The size of the clinically annotated EZ (# electrodes) was 8.05 ± 4.34. The onset age was 17.21 ± 13.48 years old, while all patients now are 34.68 ± 12.30 years old. The subject groups for each center are shown in Figure 2.

NIH Intracranial EEG Monitoring Technique: ECoG
Seven patients included in this study were surgically treated for drug-resistant seizures at the NIH National Institute of Neurological Disorders and Stroke and underwent invasive presurgical monitoring with subdural grids for seizure localization or mapping of eloquent areas. Recordings were acquired with a Nihon Kohden clinical EEG system. IEEG signals were sampled at a 1 kHz sampling rate and filtered using a 300 Hz antialiasing filter. Signals were referenced to a common contact placed subcutaneously on the scalp, on the mastoid process, or on the subdural grid. Each data file stores continuous iEEG data from all channels and is automatically generated by the acquisition system.

Johns Hopkins Hospital Intracranial EEG Monitoring Technique: ECoG
Twenty patients included in this study were surgically treated for drug-resistant seizures at the Johns Hopkins Hospital and underwent invasive presurgical monitoring with subdural grid and strip arrays for seizure localization or mapping of eloquent areas. Recordings were acquired with a Nihon Kohden clinical EEG system with a 1 kHz sampling rate and a 300 Hz antialiasing filter, and were converted to EDF format for storage and further processing. Each EDF file stores approximately 42 min of continuous ECoG data from all channels and is automatically generated by the acquisition system. Consecutive EDF files cover consecutive, nonoverlapping, time windows with less than 5s-lag in between. Digitized data were stored in an IRB-approved database compliant with HIPAA regulations.

UMMC Intracranial EEG Monitoring Technique: ECoG
Seven patients included in this study were surgically treated for drug-resistant seizures at the University Maryland School of Medicine and underwent invasive presurgical monitoring with subdural grid and strip arrays for seizure localization or mapping of eloquent areas. At the University of Maryland Medical Center (UMMC), recordings were acquired with a Natus/XLTEK system (Natus Medical Incorporated, Inc., USA) with 250-1,000 Hz sampling rate and 50-300 Hz antialiasing filter, and were converted to EDF format for storage and further processing. Each EDF file stores approximately 42 min of continuous ECoG data from all channels and is automatically generated by the acquisition system. Consecutive EDF files cover consecutive, nonoverlapping, time windows with less than 5s-lag in between. Digitized data were stored in an IRB-approved database compliant with HIPAA regulations.

Cleveland Clinic Stereotactic EEG Monitoring Technique: SEEG
Eight patients that underwent SEEG invasive monitoring from the Cleveland Clinic epilepsy center were included in this study. The choice of electrode location was based on a preimplantation patient management conference and was made independently of the present study. Criteria for patients undergoing SEEG implantation were reviewed by clinicians to determine patient eligibility for enrollment in the current study. If the patient met study criteria, research staff not involved in the surgery implantation or postsurgical care contacted the patient for potential participation in the study.
For each subject, approximately 8-13 stereotactically placed depth electrodes were implanted. The electrode contacts were 0.8 mm in diameter, 2 mm in length, and spaced 1.5 mm apart. Depth electrodes were inserted in either orthogonal or oblique orientations using a robotic surgical implantation platform (ROSA, Medtech Surgical, Inc., USA) allowing intracranial recording from lateral, intermediate, and/or deep cortical and subcortical structures in a three-dimensional arrangement (González-Martínez et al., 2016). The day prior to surgery, volumetric preoperative MRIs (T1, contrasted with Multihance 0.1 mmol/kg) were obtained and used to preoperatively plan electrode trajectories. All trajectories were evaluated for safety; any trajectory that appeared to compromise vascular structures was adjusted appropriately without affecting the sampling from areas of interest.
SEEG electrophysiological data were acquired using a conventional clinical electrophysiology acquisition system (Nihon Kohden 1200, Nihon Kohden America, USA) at a sampling rate of 1 kHz and 300 Hz antialiasing filter. Behavioral event data were simultaneously acquired during behavioral experiments along with the SEEG electrophysiology and stored for subsequent analysis. All signals were referenced to a contact affixed to the skull. Archived electrophysiological data were not filtered prior to offline analysis.
Each patient had electrode contacts characterized according to anatomical location. The anatomical locations of all contacts were identified through inspection of postoperative imaging, requiring agreement by two clinical experts. An example of postoperative imaging contributing toward determining contact location is shown in Figure 1. Coronal and sagittal views were available for every contact.

METHODS: COMPUTATIONAL STEPS
In this study, our raw dataset consisted of EEG recordings of seizures with 60 s of data before and after each seizure. Data were collected from 42 patients with at least two seizures per patient. We applied network analysis techniques and considered each electrode in the iEEG array to be a node in a network. The overall process of our algorithm is highlighted in Figure 3. We computed the cross-power spectrum matrix for each time window, then the cor-Cross-power spectrum: The Fourier transform between two pairs of iEEG time series (i.e., two different channels).
responding EVC, and then we trained a Gaussian weighting function that assigned a likelihood Gaussian weighting function: A Gaussian distributed function that applies weights to different points in the PCA space.
to each electrode for being within the EZ. After computing the heat map for the EZ predicted set of electrodes, we compared them to the clinical electrodes for both successful and failed surgical outcomes. We show results for each center separately, and also all patients grouped Computational steps for seizure onset localization: the algorithm processes raw ECoG to compute the sequence of adjacency matrix A(t). From this sequence, A(t), it computes the sequence of leading eigenvectors, v(t), as a network centrality measure, the EVC. Algorithm then converts EVC into the sequence of rank centrality r(t). From this sequence, r(t), algorithm computes a heatmap that generates predictions of the EZ. Yellow shading indicates the EVC of the first electrode evolving in time whose rank centrality, r 1 (t), is illustrated in the plot.
together. Note that we trained the Gaussian weighting function only using one center's patients, so that we could test our results across center. Clinical procedures can vary more from center to center versus the variability within center, so it is a conservative approach to train using one center and then test on all other centers to see whether our analysis holds across different clinical procedures. All Matlab (R2016b) and Python (v 2.7) code is publicly available online at Li (2018).

Preprocessing of Data
All data underwent digital filtering with a butterworth notch filter of order 4, implemented in Matlab with the f ilt f ilt function (frequency ranges of 59.5 to 60.5). In general, EEG data are known to be noisy and referencing schemes can play a significant role in downstream data analysis. We decided to apply a common average referencing scheme to the data before Common average referencing: An EEG referencing scheme in which a sample average of all the recording sites is taken and is used as the reference signal for all electrodes.
analysis (Ludwig et al., 2009). Here, we take an average signal from all recording electrodes and subtract it from the electrodes. This has been shown to produce more stable results and rejects correlated noise across many electrodes (Gliske et al., 2016). We made sure to exclude any electrodes from subsequent analysis if they were informed to have artifacts in their recording by clinicians.

Compute and Rank Nodal Centrality Over Time
Network centrality for each node was computed every second using a 2.5 s sliding window sliding every second 60 s before seizure, during seizure, and 60 s after seizure for at least two seizure events. For each window, the brain network was first represented by a connectivity matrix (Fisher, 2012), by computing all pairwise cross-power spectra between the signals in the gamma frequency band (30-90 Hz); that is, where P i , P j are the magnitudes of the Fourier transform of the time series in the window recorded from electrodes i, j, and A ij is the element of connectivity matrix and is the adjacency between nodes i and j. We chose the gamma band because the gamma frequency band has often exhibited the most modulation in power between nonseizure and seizure periods. It has been thought to be correlated to neuronal spiking and fMRI activity and thus carries information in such invasive recordings (Gotman, 1983;Worrell et al., 2004;Wu & Gotman, 1998).
The importance of each electrode to the network connectivity was measured by the strength and number of connections it makes with other electrodes, referred to as centrality. We used the eigenvector centrality (EVC) to measure the connectivity of each electrode, as EVC showed interesting repeatable patterns over seizure events in our prior study (Burns et al., 2014). The EVC of an electrode is defined as the sum of the EVCs of all other electrodes weighted by their connectivity, which measures the relative influence of a node within the network. The EVC of all electrodes is computed implicitly as (2) λ is the leading eigenvalue of the connectivity matrix A and the EVC is then the leading eigenvector of A. In simple terms, the EVC of a node in the network (electrode) is proportional to the sum of EVCs of its neighbors (nodes it is connected to). That is, a node is important if it is (a) connected to a few nodes that are themselves very important or if it is (b) connected to a very large number of not-so-important nodes. The leading eigenvectors of connectivity matrices were calculated numerically at each second during the recordings from the connectivity matrices. Finally, the EVC vector for each second was converted to a ranked vector containing values 1 to N, where a 1 was placed in the component of EVC that had the smallest centrality and an N was placed in the component of EVC that had the largest centrality.

Normalize Rank Evolution Signals
Next, we normalized the rank evolution signals (the EVC) for each electrode in the X (time) and Y (rank centrality, i.e., number of electrodes) directions. This was done so that we can compare signals from different patients that have varying numbers of electrodes and varying seizure durations across individuals and within individuals. To normalize along the X-axis, we either stretched (interpolated) or shrunk (simply downsampled at a lower sampling rate) each ranked EVC signal during a seizure epoch such that all signals were 500 data points in length. Most ranked EVC signals were under 500 s in length, so the majority of the rank centrality signals were stretched using linear interpolation (using the interp1 function in Matlab) preserving the shape of the signal during a seizure event. To normalize along the Y-axis, we scaled the rank centrality between 0 and 1 by dividing by the number of electrodes. Further, in order to compare the ranked EVC in a quantifiable manner, we normalized all the X, Y normalized signals such that the centrality signal integrated to 1. We divided the normalized rank centrality by area under the curve. This normalization converted each signal into a probability density function, where R(t) is the normalized rank signal in time after dividing by the number of electrodes andR(t) is the normalized rank signal at normalized timet.

Compute Feature Vector From Normalized Rank Signals
For each normalized signal, we extracted the deciles in time, the locations at which the signal integrates equally to 10% of the total area under the curve, that is, points in normalized time where the signal integrates to 0.1, 0.2, 0.3, and so on until the end of the signal is reached. This gives a 10-dimensional vector for each signal that serves as a feature vector.

Electrode Weight Assignment Based on Feature Vectors
Once we calculated feature vectors for each signal, we projected the features into a 2D principle component (PC) space. This was done by assuming that each feature vector is an observation, hence the analysis was performed in space x time. We performed PC analysis and plotted the features across all electrodes and patients projected onto the first and second PCs. Each electrode (data point in Figure 4A) was labeled according to whether the electrode was in the clinical annotated EZ region and whether the surgical resection was a success or a failure. We then created a weighting function over the 2D PC space, which would assign a weight to an electrode based on their location in PC space.
To generate this weighting function, we discretized it into equally sized square partitions (100 × 100 along first and second principal components). The mean normalized rank signature across all data points was computed for each partition. The signatures for the four corner partitions are shown in Figure 4A. The shapes of the mean normalized rank signatures across partitions change in a somewhat continuous manner. Moving vertically from the bottom of the PC space to the top, the rank signatures transition from a concave to a convex shape. Moving from left to right, the signature shifts horizontally: forward (to the right) if the partition is at the bottom of the PC space, and backward (to the left) if at the top of the PC space.
Our hypothesis is that the arch signature displayed in the bottom left of Figure 4A represents the signatures of the EZ because this is the region of the PC space that has the most isolated channels that come from patients with successful outcomes (green + points). In fact, the bottom portion of the PC grid shows the arch signature. Therefore, the weighting function is set to be highest in these regions and decay as a function of distance from these regions. We defined a weighting function to be the sum of four bivariate Gaussian-like functions (Equation 5, Figure 4B) as shown in Equation 5. The 2D PC space is divided into four quadrants defined by an origin. See Figure 4B (left) with origin (−100, −100).

Training Origin of Gaussian Weighting Function
In each quadrant, the bivariate Gaussian-like function was initialized with the shapes in Figure 4A. The covariance matrix in each quadrant was computed as the sample covariance from the data points in that quadrant. The origin of the four quadrants is the mean vector, which is trained. We followed a leave-one-out training procedure on the sample of 20 patients Leave-one-out: A training procedure for an algorithm that leaves one data point out during a pass-through of the training. collected at JHU. We chose JHU because it had the greatest number of patients collected within center and would still account for less then 50% of the total patients. The mean of all four quadrants is optimized for maximizing the DOA. In Figure 4B, this is shown as (−100, 100), which was found at the end. Once the optimized mean is found, then all four quadrants' Gaussian functions, w i (x, y), are linearly combined with a Heaviside step function to get the final Gaussian weighting function, w(x, y). This final Gaussian weighting function, w(x, y), is used to assign weights to all subsequent EVC of each electrode for every patient. This in turn produces the likelihood of every electrode being within the EZ set. where ; α i -exponential decay factor for ith quadrant; x − x y , and µ − μ x μ y define the position and mean vector, respectively;

Computing Degree of Agreement and Statistical Analysis
For every seizure event for every patient in NIH, UMMC, and CC, we generated a set of electrodes with their heatmap (defined by electrode weights; see Figure 3), which can be interpreted as their likelihood for being in the EZ. For each seizure recording, we then computed the degree of agreement between the computed EZ likelihoods and clinical annotations of the EZ. The likelihood was computed using the Gaussian weighting function trained as described in the previous subsection. Then, a threshold α = 0.3, 0.6, 0.9 was applied to each heatmap, and the set of electrodes whose likelihoods exceeded α were defined as the algorithm's EZ (AEZ). The AEZ was then compared with clinically annotated EZ (CEZ) using the following degree of agreement (DOA) statistic: Note thatS is the complement of the set S, and that D ∈ [−1, 1], where DOA = 1 implies perfect agreement and DOA < 0 is less agreement.
Across all patients, electrodes, and seizure events, we have a collection of DOA values. We then derive two distributions: (a) the distribution of DOA for all electrodes implanted in patients who had successful treatments, and (b) the distribution of DOA for all electrodes implanted in patients who had failed treatments. We then test whether there is a significant difference in DOA distribution between these two patient groups using the Wilcoxon rank sum test to test for statistical differences. This nonparametric test was selected, as the data are not guaranteed to meet the normality conditions for a Student's t test (Whitley & Ball, 2002). In addition, we also added an across-center analysis where we combine all the data and test whether the DOA distributions for successful versus failed outcomes are significantly different.
On top of this analysis, we also add a min-max scaling to normalize the degree of agreements within each center, so that success and failure could be compared at the same scale.

High-Frequency Oscillator: qHFO Detector
We compared our algorithm with the qHFO algorithm presented in Gliske et al. (2016), which uses a sensitive HFO detector, then redacts HFOs that were produced by artifacts. Previous work has shown that sampling rates of 1000 Hz are capable of recording HFOs, but only capture 60% of the events Gliske et al. (2016). Therefore, we only analyzed patients with sampling rates ≥ 1,000 Hz and with available interictal data. This resulted in three patients from NIH and two patients from JHU, with a total of 13 separate recorded datasets. The datasets here analyzed had an average recording of 7.1 min, 83 total electrodes analyzed, and 10 electrodes within the clinically annotated EZ set. Using the qHFO algorithm on this data required a few minor adaptions.
We used a single common average reference applied to all analyzed intracranial electrodes (as described earlier), rather than separating the referencing between depth-electrode channels and grid channels as was done in Gliske et al. (2016). The popDet artifact rejection method also could not be used, as it requires sampling rates of at least 2,000 Hz.

RESULTS
Every patient (n = 42) with at least two seizures was analyzed (total of 113 seizures), with 20 of the patients from JHU used to train the final Gaussian weighting function. The output of the process for each seizure recording is each electrode's likelihood of being in the EZ. These likelihood scores are in turn used to produce a heatmap that can be overlaid on a brain MRI to show the relative predicted EZ region for a certain patient. Figure 5 shows a few examples of heatmaps for three patients who had successful outcomes and three patients with failed Figure 5. This figure shows an example overlay of the algorithm's heatmap of likelihood on a brain scan for six patients (three successful and three failed outcomes). The red region shows our predicted onset zone and the black outlines represent where the clinicians performed a resection. The orange, yellow, green, and blue regions represent lower likelihoods for that specific electrode being within the EZ set as predicted by the algorithm. outcomes. For the three successful patients, the AEZ lies entirely within the resected regions, suggesting a high DOA between the AEZ and CEZ. For one of the failed patients, the resected region and the AEZ do not overlap, that is, DOA is low. For the other failed patient, the AEZ is a very small set, suggesting that the EZ may not be appropriately covered by the electrode implantation.
In our comparative HFO analysis, we analyzed 13 segments of data from five patients. Of the 13 files, most patients have no HFOs, even at 1,000 Hz sampling rate (see Table 1). Only three data segments had HFO detections, but one of them did not have an anomalous grouping suggestive of the EZ (30% of the total recording time from all 13 data segments). In JH3, there were HFOs, but no channels had an anomalous rate high enough to be predicted within the EZ set. In NIH pt1aw2 and pt3alsp3, both only had a single channel predicted to be in the EZ. This prediction was in concordance with clinically annotated EZ in pt1 but not in pt3.
The lower sampling rate and short time segments are not ideal for automated HFO analysis, as is apparent from these results. In our network analysis, we had a high DOA with pt1 (0.62), while a relatively lower DOA for pt3 (−0.16). It seemed that for pt3, HFO analysis completely disagreed with clinical annotations, while the network analysis found more electrodes then the clinically annotated EZ, which led to lower DOA. For pt1, the network analysis also highlighted the same electrode as being in the EZ set. This shows how HFO and network analysis can complement each other in analyzing different sections of the data. Based on our limited comparisons due to inherent data limitations, our analysis is more capable of identifying the full clinically annotated EZ than HFOs in this specific dataset.
In Figure 6, we show the DOA for datasets collected from the test datasets (the three clinical centers: UMMC, NIH, CC) for three different threshold values, α, that are placed on the likelihood distribution (electrodes with likelihood greater then threshold are placed in the EZ set). The resulting DOA after training the Gaussian weighting function for JHU are shown in Supplementary Table 5 . It also shows the same trend as seen in Figure 6. As illustrated in Figure 6, the general trend is that the DOA distributions for successes and failures separate more as α increases, and α = 0.9 appears to be an operative threshold that shows a positive DOA for successes and a negative DOA for failed outcomes. For α = 0.9, the statistics Figure 6. This figure shows degrees of agreement using the degree of agreement index between our algorithm and clinical annotations for successful and failed surgical resections. The dashed line at DOA = 0 represents neither agreement nor disagreement. The red line is the average DOA, and the blue box is the box plot of the DOA; −1 is a perfect disagreement between the algorithm and clinical set, while 1 is a perfect agreement between the algorithm and clinical set. Table 2. Degree of agreement results for α = 0.9 with average ± standard deviation from each clinical center and also the resulting p value from the Wilcoxon rank-sum test. All centers show a significant difference between success and failure cases. Note JHH is used in the training of the Gaussian weighting function. for DOA (mean and standard deviation) are given in Table 2 for each center and across all centers together. By applying a Wilcoxon rank-sum test, we also see a significant difference at significance level 0.05 for all centers at threshold level of 0.9. At each center, there is a trend of the DOA that is a function of clinical outcome of the patient. This is consistently shown across recording platform (ECoG for UMMC, NIH and SEEG for CC) and patient population. In all cases, as the threshold increases from 0.3 to 0.9, the difference of DOA between successful and failed cases increases. If there is low DOA with the algorithms EZ and the clinically annotated EZ and the patient is a failed outcome, then this may be a case of mislocalization. If, on the other hand, there is no visible EZ from the algorithm (all weights are low), then the EZ may not be in the vicinity of the electrode, suggesting a possible misplacement of electrodes.

Center
We also show in Figure 7 that there is no bias due to center (results are shown in Table 3. All centers, when normalized, show a significant difference between successes and failures. Note that min-max normalization scales all distributions between 0 and 1. Table 3. Degree of agreement results for α = 0.9 with average ± standard deviation from each clinical center after min-max scaling and also the resulting p value from the Wilcoxon rank-sum test. All centers show a significant difference between success and failure cases. The large variation is due to the varying number of electrodes implanted per patient and the varying size of the clinical EZ hypothesis. However, all centers show significant difference when compared with a Wilcoxon rank-sum test.

Center
In the case that a patient has failed outcomes, we would not expect to see a perfect disagreement DOA score of −1 because of the above reasons. There may have been no visible EZ recorded from the electrode network, or the EZ may not have been fully resected (but part of it was still clinically annotated). It is also important to note that when a patient has a successful surgical outcome, clinicians remove a large portion of the brain, which is a superset of the clinically annotated EZ. It is not certain that all clinically annotated EZ electrodes are actually part of the true underlying EZ, so we would expect some deviation from perfect agreement with the clinically annotated EZ (e.g., we should not expect to see a perfect DOA score of 1 for successful patients).

DISCUSSION
The definition of the EZ, including its anatomical and electrophysiological signatures, has been an evolving and controversial topic since the foundation of modern epilepsy surgery. The EZ, defined as the site of primary organization of the ictal discharge, refers to the cortical areas connected together through an excessive synchronization at seizure onset (Talairach & Bancaud, 1973;Wendling, Chauvel, Biraben, & Bartolomei, 2010). Fast activity (FA) at ictal onset has been clinically accepted as the main feature of the EZ since the beginning of the invasive monitoring era, particularly in the SEEG literature (Talairach & Bancaud, 1973). Since the development of subdural ECoG recordings, much attention has also been paid to the time precedence of phasic transients, especially spiking activities (Boonyapisit et al., 2003;Palmini et al., 1995). In the last 15 years, identification of high-frequency oscillations (HFO) during interictal and ictal periods in experimental models reoriented research interest toward high-gamma activities in human epilepsies as a potential EZ marker (Bragin et al., 2002;Matsumoto et al., 2013;Zijlmans et al., 2012). In parallel, DC recordings exemplified the concomitance of fast and ultra-slow frequencies (Gnatkovsky et al., 2014;Ikeda et al., 1996;Thompson et al., 2016;Wu & Gotman, 1998), which could be used as potential biomarkers of the EZ.
Although clinical definitions have been explored, a network-based operational definition of the EZ is currently not well defined in the literature. Novel computational network analyses may overcome some of the challenges associated with more conventional invasive monitoring recordings methods. In this study, we analyze how centrality signatures of electrode recordings within an epileptic network change over time and how they relate to clinical annotations from four different hospital centers. We take in ECoG and SEEG data from 60 s before and after a seizure instance for 42 patients and produce a frequency connectivity network over time using the cross-power spectra of the signal in the 30-90 Hz range. Then we computed the EVC for each electrode at a time window to obtain a normalized ranked centrality of every electrode over time. By overlaying a Gaussian weighting function that was trained only with patients from one center, we then obtain a likelihood for each electrode of being in the EZ. Then we computed a degree of agreement between our algorithm and clinically labeled EZ using the DOA index for all patients by setting an arbitrary threshold.
Some previous approaches for marking the EZ included FA, signal flattening, and slow potential shift. Fast activity frequently occurs quasisimultaneously in multiple areas so that visual discrimination can be cumbersome and can lead to subjective interpretations. A different approach, frequency localization, was used by Gnatkovsky et al. (2014). After defining frequencies of interest (FOIs) and plotting their power change over time, they localized the distribution of FOIs in different contacts of the depth electrodes. The EZ, defined as the area exhibiting frequency changes at seizure onset, could then be delineated. In a retrospective and prospective study of patients investigated using SEEG, the same method was applied to test three potential biomarkers of EZ, namely FA, signal flattening, and slow potential shift. These biomarkers colocalized with the EZ as defined by standard SEEG criteria and postresection seizure outcome (Gnatkovsky et al., 2014).
Other approaches for marking the EZ include HFO analyses. Interictal HFOs have been shown to have some value in identifying the EZ (Burnos et al., 2016;Jacobs et al., 2012;Xiang, 2008;Malinowska, Bergey, Harezlak, & Jouny, 2015;Usui et al., 2011;Van Klink et al., 2014). In our comparative analysis, we made modifications to the algorithm based on limitations in the data that were available at the clinical centers. First, in the 1,000 Hz sampled data, the number of HFOs is significantly reduced, although the detected HFOs are still useful to identify the EZ (Gliske et al., 2016). The lower sampling rate also required some modifications to the algorithm: the fast-transient artifact detector could not be used (as it requires sampling rates > 2 kHz) and the upper edge on the band-pass filter needed to be reduced from 500 to 400 Hz. Second, the limitation to interictal data restricts the identification of the full EZ: HFO results typically report a very small number of channels involved, which are typically much smaller than the eventual resected volume of tissue. Although HFO analyses show promise in analyzing electrophysiology of epileptic patients, they do not take into account the network nature of epilepsy. HFO analyses are important for analyses of interictal data, since our analysis is limited by requiring recorded seizure events. In future studies, it would be interesting to see how network algorithms and HFO algorithms can complement each other to improve EZ localization.
It is important to note that network-based analyses are not new to analyzing EEG recordings from epilepsy patients. Previous studies have shown that seizure activity is a dynamic multichannel process, and the correlation structure right around a seizure event also follows a typical evolution, similar to our ranked EVC signal (Kramer, Kolaczyk, & Kirsch, 2008;Schindler et al., 2007). In Kramer et al. (2008) and Schindler et al. (2007), they do not relate it back to EZ, but just look at network dynamics during seizure events. In Schevon et al. (2007), the authors compute interelectrode synchrony using the mean phase coherence algorithm and relate locally synchronous EEG channels back to the EZ, but analyzed only nine patients from a single center. A similar small-scale study was performed in Korzeniewska et al. (2014) with six epilepsy patients from one center. Other studies use computational models to understand the biophysical mechanisms related to epilepsy surgery (Khambhati et al., 2016;Sinha et al., 2017). In Khambhati et al. (2016), they applied a virtual resection model using data from 10 patients. In Sinha et al. (2017), the authors developed patient-specific dynamical network models of epileptogenic cortex (computational models). However, there were only 16 patients analyzed from one center.
This manuscript describes a somewhat large-scale research study that applies networkbased data analysis tools to invasive EEG data to explore possible EEG signatures of the EZ. In no way are we proposing that this algorithm be directly translated into the clinic. Rather, it now compares how pairwise correlations may improve over quantifying HFOs in each channel individually, which has been the most recently accepted approach. We present a network analysis related back to the annotated EZ, analyzing data from before and after seizures, and analyzing data from multiple centers (with 113 seizures from 42 patients). In our study, we showed that there is a general higher degree of agreement between our algorithm and clinically successful surgical resections of the EZ and a lower degree of agreement between our algorithm and clinically failed resections. By setting a simple threshold on the likelihood maps, we can obtain a similarity measure between our algorithm and clinical labels for both successful and failed surgeries. As the threshold increases, our algorithm becomes better at identifying whether successful resections had the correct EZ. We observed that the algorithm's performance degraded with respect to degree of agreement when patients were implanted sparsely with single strips across all four lobes (UMMC patients) and sometimes in both hemispheres. The clinicians place these strips with such wide coverage if there is no clear preimplantation hypothesis and if seizures are thought to be starting from multiple brain regions. Often, these patients do not have clear EZ localization and/or do not end up as candidates for surgery. We also found that if the electrographic onset of seizure is not close to the clinically annotated onset of seizure, then the degree of agreement with clinicians is reduced. The electrographic onset is the start of seizure that is seen on the EEG recordings but not manifested in any behavioral changes in the patient. The clinical onset is the time at which the patient exhibits noticeable behavioral changes due to seizure onset (e.g., muscle twitches).
Our results suggest that network data analytics may be a useful tool to assist in localization of the epileptogenic zone, especially when electrode implantation covers the EZ network densely. This is expected, since the threshold on the network's likelihood is essentially a threshold on the algorithm's confidence in an electrode being within the EZ set. Future work entails exploring different weighting functions applied over the rank centrality space and possibly merging features from HFO and network algorithms. Besides looking solely at gamma power (30-90 Hz) cross-power matrices, the work could expand to encompass more frequency bands that could contain signals of importance in EZ localization. In addition, a more comprehensive study that compares the outcomes between SEEG and ECoG could help understand limitations of the algorithm, and also be of clinical importance in using SEEG versus ECoG. In addition, if we had more patient data from other centers, then it would be interesting to see how a pooled training procedure may improve our results. This work is meant to supplement the growing evidence in the literature that epilepsy is a network phenomenon and therefore also requires network algorithms to better understand its manifestation.

SUPPORTING INFORMATION
Code is open source at https://github.com/ncsl/eztrack. Since this was a retrospective data study, there is no table of the JHU patients and their clinical operative notes because the data were not available from JHU.