Interictal Gamma Event Connectivity Differentiates the Seizure Network and Outcome in Patients after Temporal Lobe Epilepsy Surgery

Abstract Studies of interictal EEG functional connectivity in the epileptic brain seek to identify abnormal interactions between brain regions involved in generating seizures, which clinically often is defined by the seizure onset zone (SOZ). However, there is evidence for abnormal connectivity outside the SOZ (NSOZ), and removal of the SOZ does not always result in seizure control, suggesting, in some cases, that the extent of abnormal connectivity indicates a larger seizure network than the SOZ. To better understand the potential differences in interictal functional connectivity in relation to the seizure network and outcome, we computed event connectivity in the theta (4–8 Hz, ThEC), low-gamma (30–55 Hz, LGEC), and high-gamma (65–95 Hz, HGEC) bands from interictal depth EEG recorded in surgical patients with medication-resistant seizures suspected to begin in the temporal lobe. Analysis finds stronger LGEC and HGEC in SOZ than NSOZ of seizure-free (SF) patients (p = 1.10e-9, 0.0217), but no difference in not seizure-free (NSF) patients. There were stronger LGEC and HGEC between mesial and lateral temporal SOZ of SF than NSF patients (p = 0.00114, 0.00205), and stronger LGEC and ThEC in NSOZ of NSF than SF patients (p = 0.0089, 0.0111). These results show that event connectivity is sensitive to differences in the interactions between regions in SOZ and NSOZ and SF and NSF patients. Patients with differential strengths in event connectivity could represent a well-localized seizure network, whereas an absence of differences could indicate a larger seizure network than the one localized by the SOZ and higher likelihood for seizure recurrence.


Introduction
Multimodal techniques and new signal processing approaches, such as functional connectivity analysis, have advanced the concept of epilepsy as a brain network disorder (Spencer, 2002;Bartolomei et al., 2008aBartolomei et al., ,b, 2017Amiri et al., 2020;Amorim-Leite et al., 2020;Gupta et al., 2020) and suggestions to not only find epileptogenic tissue but the network generating the seizures (Spencer, 2002;Prasad et al., 2003;Boling et al., 2009;Jehi, 2015). Motivation for identifying the seizure network is readily found in cases of medication-resistant epilepsy, where in current practice, removal of the seizure onset zone (SOZ) does not always control seizures (Prasad et al., 2003;Boling et al., 2009). Presently, however, the extent of structural anomalies and functional disturbances that characterize the seizure network, how these disturbances generate seizures, and which critical portions of the network need to be removed to abolish seizures, is unknown.
Studies of the seizure network using interictal EEG functional connectivity suggest brain regions in the SOZ are more strongly connected than regions not part of the SOZ (NSOZ) and possibly disconnected from the NSOZ (Warren et al., 2010;Bettus et al., 2011). Also, more connectivity alterations in NSOZ correlate with a larger epileptogenic network (Lagarde et al., 2018). Undoubtedly multiple, complex mechanisms contribute to differences in the strength of connectivity, and we believe the basis for this involves the synchrony of excitatory and inhibitory activity that could be greater in regions generating seizures than those not (Jiang et al., 2022). If this hypothesis is correct, we reasoned gamma-band connectivity might detect differences in synchrony since gamma involves coordinated synaptic activity of excitatory and inhibitory cells (Bartos et al., 2007;Buzsáki and Wang, 2012;Chen et al., 2017;Gu et al., 2021), is associated with excitatory-inhibitory balance (Gao et al., 2017), and power positively correlates with neuronal spiking rate (Mukamel et al., 2005;Manning et al., 2009). In addition, we computed theta-band connectivity because others had found differences in the theta power between mesial temporal and extratemporal lobe regions involved in generating seizures (Bettus et al., 2008).
There are a number of approaches to measure functional connectivity, including correlation (Adey et al., 1961;Alonso et al., 1996;Pfurtscheller and Andrew, 1999), phase-based methods (Lachaux et al., 1999;Mormann et al., 2000;Reijneveld et al., 2007), and information theory (Afshani et al., 2019;Ursino et al., 2020). Among these approaches is event connectivity, which combines aspects of correlation and information theory (Kheiri et al., 2013). Although not used extensively and, to our knowledge, not in patient studies of epilepsy, gamma event connectivity in rats produces stable values within behavioral states, correlates with neuronal discharges, and is sensitive to changes in excitatory and inhibitory synaptic activity (Bragin et al., 2014). Based on these data, event connectivity appears well suited for our purposes to measure functional connectivity and indirectly the synchrony of inhibitory and excitatory activity associated with the seizure network.
In the current study, we computed the strength of event connectivity in theta (4-8 Hz), low-gamma (30-55 Hz), and high-gamma (65-95 Hz) bands from interictal EEG recorded between pairs of contacts on intracerebral electrodes implanted in patients who had resective surgery or received an electrical stimulation device to control their seizures. We predicted stronger synchrony and thus event connectivity between brain regions in the SOZ than NSOZ, and larger differences in strength of connectivity between SOZ and NSOZ in seizure-free than not seizure patients, which we suspected could be due increased synchrony in some regions of NSOZ that are involved in generating seizures in not seizure-free patients.

Subjects and clinical recordings
All 43 subjects (26 females, 17 males) for this retrospective study were patients with medication-resistant focal seizures suspected to begin in the temporal lobe and candidates for epilepsy surgery, but required intracranial depth electrode EEG (iEEG) studies to localize the brain area of seizure onset. All patients were bilaterally implanted with seven-to nine-contact clinical depth electrodes (Ad-Tech Medical Instrument) oriented perpendicular to the lateral surface of the temporal bone and positioned to sample amygdala, entorhinal cortex, hippocampus, and parahippocampal gyrus, as well as extratemporal areas such as orbitofrontal cortex, anterior cingulate gyrus, supplementary motor areas, or parietal cortex (Table 1). Patients were recorded for 7-14 days in the epilepsy monitoring unit until a sufficient number of the patient's habitual spontaneous seizures were captured. Depth EEG recordings were reviewed by the attending neurologist who identified the electrode contacts where seizures first appeared, which were labeled as the seizure onset zone (SOZ). All remaining contacts were considered outside or non-SOZ (NSOZ). Informed consent was obtained from each patient before the implantation of depth electrodes and participating in this research, which was approved by the Medical Institutional Review Board 3 (10-001452).

Depth electrode recordings and localization
Interictal depth EEG recordings were reviewed to remove signals containing electrical noise and the remaining signals notched filtered at 60 Hz. For each patient postoperative CT scans were registered to preoperative MRI to identify electrode contacts within gray matter, and those contacts fully in white matter were excluded from the analysis as they would induce spurious connections that are a result of volume conduction. These steps yielded a total of 2055 electrode contacts without electrical noise with an average of 49 6 16 contacts per patient. For each patient a 10-to 15-min interictal depth EEG recordings were selected using the following criteria: (1) .24 h after electrode implantation, (2) before tapering of anti-seizure medications, (3) at least 6 h before the first recorded seizure, and (4) period of quiet wakefulness with eyes open or closed. Only seizures, as an epileptiform activity, were avoided. All other interictal discharges, including epileptic spikes could have appeared in the selected data portions. Fifteen patients were recorded with a sampling frequency of 200 Hz, and 28 patients had recordings of 2 kHz. All recordings were resampled to 1 kHz using the MATLAB anti-aliasing resample function before connectivity measure was calculated. To verify sampling rate did not affect connectivity, especially with high gamma (65-95 Hz), we compared (1) the ratio of low (30-55 Hz) to high gamma power, and (2) the ratio of number of events of low to high gamma detected using the MATLAB function findpeaks between patient data with different sample rates. For each patient we calculated ratios from a randomly selected 30-s window on five channels and repeated the procedure 10 times, which generated 1400ð8 Â 5 Â 10) datapoints for the patients sampled at 2 kHz and 750ð15 Â 5 Â 10Þ datapoints for the patients sampled at 200 Hz. Results show a significant, but small effect, of sampling rate on the ratios of low to high g power (Wilcoxon test, p = 7:595 Â 10 À35 , Cohen's d = 0.2001) and number of events (4:6258 Â 10 À28 , Cohen's d = 0.105; Extended Data Fig. 1-1A,B), suggesting the anti-aliasing filter had only a small effect on connectivity. Also, oversampling to 2 kHz produced more events needed in the perievent histogram (see below, Connectivity metrics) but did not affect the spectral frequency components of the signal (see Extended Data Fig. 1-1C,D).

Connectivity metrics
Connectivity measures used in previous studies are diverse. Generally, functional connectivity methods can be divided into three main categories: (1) amplitude-based measures such as different variants of the well-known amplitude correlation/coherence in time/frequency domains (Adey et al., 1961;Alonso et al., 1996;Pfurtscheller and Andrew, 1999); (2) phase-based measures where phase locking value (Lachaux et al., 1999), mean phase coherence (Mormann et al., 2000), and phase lag index (Reijneveld et al., 2007) are most frequently used; and finally; (3) connectivity originating from information theory like mutual information (Afshani et al., 2019) and transfer entropy (Ursino et al., 2020). Connectivity methods based on information theory can capture the nonlinear interactions between pairs of brain regions without a prior assumption of a predefined linear or nonlinear model that the oscillatory phase/amplitude coupling methods are usually bound to. To exploit the benefits of both correlation and information theory, we chose to use a stable connectivity measure called gamma event coupling, initially described by Bragin et al. (2014). This method is very similar to transfer entropy and mutual information where all of them use Shannon entropy to assess the strength of connectivity of a joint probability distribution but differ in the way the distribution in constructed from the available data.
The method was adapted with different windows sizes to accommodate connectivity for theta (Theta, 4-8 Hz), low gamma (LGEC,, and high gamma (HGEC, 65-95 Hz; Fig. 1A) event connectivity. Event connectivity was estimated based on the temporal relation between individual cycles or events of theta, low gamma, and high Figure 1. EEG analysis pipeline. A, Unfiltered intracerebral EEG signals are bandpass filtered to extract spectral frequencies in theta (u ), low gammaðLg ), or high gamma bands (Hg Þ. Functional coupling between a pair of channels (ch x , ch y ) is illustrated in second row. A frequency-dependent time interval L (30 s for Hg , 60 s for Lg , and 5 min for u ) is chosen, and from the signals on ch x and ch y , the local event's amplitude maxima e i (i = 1 ... n) in L are detected (represented in red traces). For each e i in ch x , the lead or lag in relation to events in ch y within time interval [-T, T] is quantified in a peri-event histogram (bottom left). The distribution of the histogram is evaluated using Shannon entropy, and a low entropy value is an indication of a peak in the histogram, which represents the strength of functional coupling for every pair of channels in the connectivity matrix (bottom right). Patients with 2-kHz sampling rate (N = 15) and those with 200-Hz sampling rate (N = 28) were both used in this study, refer to Extended Data Figure 1-1 for detailed justification. B, Spikes are detected from unfiltered interictal data using an automatic detector based on signal whitening. The gray boxes show the detected spikes on different channels. For every pair of contact coupling strength is computed as a rate of the sum of spikes on each channel divided by the recording duration in minutes. C, Statistical model includes EEG recordings to generate functional connectivity matrix (black box) and the spikes matrix (red box), patients information and test results to assess seizure onset zone (SOZ), surgery outcome, and other measures (e.g., seizure frequency), and CT scans co-registered to MRI scans to localize electrode contacts, group contacts with respect to brain region (green box), and calculate the distance between each pair of contacts to generate distance matrices (red box). Ipsilateral and contralateral grouping was ignored (see Extended Data Fig. 1-2). gamma recorded on every pair of electrodes contacts. Note that the band 55-65 Hz was omitted to reduce the chances of spurious connections resulting from 60-Hz powerline noise contamination. First, data were either down-sampled or up-sampled to 1 kHz then bandpass filtered (FIR, order) into the respective spectral frequency bands (see Fig. 1A). For each frequency band local amplitude maxima were detected using the "findpeaks" algorithm in the MATLAB toolbox where we used a "Threshold" parameter of value 0.1 (see Fig. 1A, red signal). To measure connectivity between contacts or "channels" (ch x and ch y ), peri-event histograms were used to quantify the lead or lag between each local maxima on ch y and ch x (e i where i = 1...n, where n is total number of events in an interval of length L) within time interval [-T, T]. The values of L and T were adjusted as a function of the targeted frequency bands. According to the Nyquist rate, the highest observable frequency of events should be half the sampling rate (1000 Hz/ 2 = 500 Hz). A time resolution of 2 ms (1/500 Hz) can be used to distinguish two cases. As a result, we chose a bin size of 2 ms. A frequency dependent time window T was chosen as 1/fmin where fmin is the minimum frequency at which a related event might occur. In case of low gamma (Lg ), the selected frequency band has a minimum of fmin = 30 Hz, thus T is 1/30 Hz%34 ms. A statistically valid histogram contains at least 30 data points per bin; therefore, 1020 (34 bins Â 30 events) events need to be collected with minimum file duration of 24 s [%1020 events/42.5 Hz, where 42.5 Hz = (fmin 1 fmax)/2]. Based on these calculations, a window length of L = 30, 60, and 300 s was used for Hg ; Lg , and u , respectively. This resulted in a 3D matrix of size M Â N c Â N c , where N c is the number of channels and M is the number of matrices corresponding to different windows. An average over all M windows was then calculated for each patient.
When a peri-event histogram had a large peak, the two channels from the histogram were considered functionally related. Shannon entropy (S) was used to determine the peak's power in the histogram, which is defined as: where N is the total number of bins and p i is the probability of an event belonging to the ith bin. A lower S signifies stronger connectivity and a higher S represents weaker connectivity and a quasi-uniform distribution of events. Hence the maximum value of S would occur when all events have the same likelihood of occurrence (p i ¼ 1=N), thus S max is defined as: The Shannon entropy value of each pair of channels (i and j) were then normalized by subtracting it from its maximum S max and then dividing by it by S max as follows: The obtained value h ij represent the connectivity index (strength) between two channels (i and j), it has a minimum value of "0" that means fully disconnected and maximum value of "1" representing a fully connected pair. Connectivity was organized into a "connectivity matrix" where the ith row and jth column of the matrix correspond to the connectivity strength between channels i and j: Note that the connectivity matrix is a symmetrical matrix, i.e., h ij ¼ h ji .

Spikes coupling rate
In each patient's recording, interictal spikes (IIS) were detected using an automatic algorithm based on whitening of the power spectrum (Roehri et al., 2016). The output of the algorithm was visually inspected to ensure correct detection of spikes (Fig. 1B). For a quantitative validation, we calculated the percentage of channels with top 5% spike rates from the automated spike detection and compared these channels to those labeled by the neurologist as channels with interictal discharges. Results are summarized in Table 2. Like functional coupling in the theta and gamma frequencies, the strength of IIS coupling between two channels ch i and ch j was computed as the sum of the total number of spikes on each channel divided by the total duration in seconds. The spikes rates were organized into matrices such that the coupling rate r ij found at the row i and column j represented the spikes coupling rate between channels i and j.

Euclidean distance connectivity
After electrodes contacts were localized. First, the anatomic image is co-registered with the CT image to mask nonbrain signals. The masked CT image is then processed (thresholded, eroded, Gaussian filtered, multiplied) to highlight electrode locations. This highlighted CT is then transformed to MNI space and loaded into iElectrodes toolbox (Blenkmann et al., 2017), where electrodes contacts were localized, labeled, and indexed. iElectrodes toolbox is a comprehensive open-source toolbox for depth and subdural grid electrode localization. The x, y, and z coordinates for each contact in gray matter was extracted according to the MNI system of coordinates whose origin is situated anterior commissure and has an RAS orientation. The unit of measurement was the millimeter. The Euclidean distances were then arranged into a distance matrix (Fig. 1C), where the distance d ij found on row iand column j represented the distance between channels i and j.

Exponential model
To assess the change in connectivity strength in relation to distance we used an exponential decay model of the form: where s represents the strength of the connectivity measure, A is the hypothetical strength at distance zero, d is the Euclidean distance between channels, and t is the constant representing the rate at which the strength decays. As value of t increases the connectivity strength weakens and reaches zero faster. For each patient, the model was fitted to the connectivity strength for each frequency band as a function of distance.

Grouping of contacts and networks
Electrode contacts were grouped into mesial (M), lateral (L), and extratemporal lobe (E), which largely were in frontal lobe and rarely in parietal or occipital lobes. A pilot analysis showed there was no difference between ipsilateral NSOZ and contralateral NSOZ (P NSOZÀcontraÀinpsi ¼ 0:086; h 2 ¼ 0:000125; see Extended Data Fig. 1-2), and for this reason, the ipsilateral and contralateral NSOZs were combined. Brain network connectivity in relation to the SOZ was labeled as "inside" when both contacts were part of the SOZ, "outside" when both contacts were part of the NSOZ, or "between" when one contact was part of the SOZ and the other part of the NSOZ. A similar approach was adopted for network connectivity in relation to brain regions. Since all channels were labeled M, L, or E the six  Percentage  1  RAH1-3, RA1-3, REC1-3, RPHG1-3, LA1-3  RAH1, RAH 2, REC 1  100%  4  RA1-3, REC1-3, RAH1-3, RPHG1-3, LA1-3, LEC1-3, LAH1-3, LPHG1-3 RA1, RA2, LAH2, LAH3, LAH1  100%  5  RAH1-3, RA1-3, RPH1-3, LAH1-3, REC1-3  RPH3, RAH3, RA2, ROF3  75%  7 RA5-6, REC5-6, RAH1-6, RPHG3-6, RSTGP2-3, RMC1-4, LAH1-2 RA2, RAH1, RAH2, RPHG6, RPSTG2, Twelve patients showed 100% correspondence in the top 5% of channels with the highest spike rates between the automated spike detector and those manually identified channels containing interictal discharges. Twenty patients showed at least 50%, 11 of which with .70% correspondence. Only four patients ,40% possible regional networks were M-M, M-L, M-E, L-L, L-E, and E-E. Initially, the mean connectivity strength for each brain region and SOZ network (i.e., SOZ, NSOZ, SOZ-NSOZ) was computed to evaluate connectivity between seizurefree and not seizure-free patients. However, we were not able assess the interactions between the brain region and SOZ. For example, if we consider one contact of a given pair, it might be in the mesial temporal region (M) and the other in the Lateral temporal region (L), i.e., part of the M-L network. At the same time, both electrodes might be in the SOZ and thus the connectivity is part of the SOZ network. This does not hold for all electrodes in the M-L network, i.e., not all contacts in the M-L network are necessary in the SOZ. Thus, calculating an average value for SOZ connectivity means ignoring the regions networks or vice-versa, and to consider the interaction between anatomic regions and zones, individual nonaveraged connectivity values were considered.

Statistical analysis
To examine the effects of SOZ, brain region, and seizure outcome on HGEC while controlling for IIS rates and interelectrode distance, a linear mixed model was used with (1) HGEC as dependent measure; (2) SOZ, brain region, and seizure outcome as independent variables (fixed slopes); and (3) IIS and distance as covariates. The intercepts arising from different fits for each subject was set to be random. Dependent variables that could not be transformed into normal distribution were analyzed with nonparametric Wilcoxon test. The magnitude of difference was calculated as the difference between the estimated marginal means of groups. Cohen's d (Cohen, 1988) was used to compute effect size for Wilcoxon test. Bonferroni was used to correct for multicomparisons. Pearson correlation was used to assess the linear relationship between interelectrode distance and event connectivity measures. All statistical tests were performed using SPSS software (IBM SPSS Statistics for Windows; IBM Corp, 2020) except for the nonparametric tests, which were conducted using the Statistics Toolbox of MATLAB software (The MathWorks, MATLAB, version 2020a).

Code accessibility
Gamma event connectivity code is freely available at: https://github.com/MohamadShamas/GEC.git. To help in replication of the results, we provide a small dataset of five patients (one patient with 200-Hz sampling frequency and four patients with 2-kHz sampling frequency) on the same link. Instructions on how to use the code are listed in the readme.md file.

Patient cohort
Forty-three patients (n = 26 females; mean age of 44.3 6 10.6 years) with predominately temporal lobe seizures, surgical treatment, and seizure outcome were included in the study (Table 3). Results from depth electrode recording showed seizures began in unilateral or bilateral temporal lobe structures of 30 and 8 patients, respectively, temporal and ipsilateral frontal or parietal lobe in three patients, and bilateral temporal and frontal lobe in two patients. Eighteen were seizure free with an Engel score of IA or IB, and 25 continued to have clinical seizures after resective or RNS surgery (Engel Class IC to IVC; Extended Data Table 3-1) with average follow-up of 3.25 (62.05) years. The proportion of females to males and median age at surgery was similar between seizure outcome groups, median age in seizure-free group was 52 years and in not-seizurefree group was 39 years old (nonparametric Wilcoxon test, P age = 0.0784). There was no significant difference in frequency of seizures or auras between the seizurefree and not seizure-free groups (Wilcoxon, P Seizure Freq = 0.434 and P auras Freq = 0.832, respectively) or in the duration of epilepsy (median duration 26 vs 13 years; Wilcoxon, P duration = 0.0883).
Connectivity in relation to SOZ, brain region, and seizure outcome For each patient we constructed connectivity matrices computed from a 10-to 15-min period of interictal depth EEG (see Materials and Methods). Inspection of the matrices, like the example of HGEC illustrated in Figure  2A, revealed stronger connectivity in seizure-free than not seizure-free patients (see Extended Data Fig. 2-1). Arranging the electrodes in relation to the SOZ ( Fig. 2A, top row) or brain region ( Fig. 2A, bottom row) also indicated differences in connectivity in many, but not all, patients (Fig. 2B). To verify these observations, we used linear mixed model analysis to evaluate the effect seizure outcome as well as SOZ and brain region on Theta, LGEC, and HGEC. We included the rate of interictal spikes and intercontact distance as covariates in the model since each of these could affect connectivity (Ren et al., 2015;Lagarde et al., 2018).
Results from the linear mixed model found seizure outcome did not have a significant effect on connectivity nor did SOZ, except on HGEC, and brain region did have a significant effect on HGEC, LGEC, and Theta (see Table 4). No significant differences were obtained when comparing different zones and seizure outcomes for all three frequency bands (see Extended Data Fig. 3-1A,B). Delving into the model's results (i.e., interactions), results show stronger HGEC and LGEC in the SOZ than NSOZ of seizure-free patients but no difference in not seizure-free patients (see Fig. 3). Furthermore, seizurefree patients had stronger HGEC and LGEC in the SOZ between mesial and lateral temporal lobe than not seizure patients (see Fig. 4A,B). By contrast, seizure-free patients had weaker LGEC in the extratemporal NSOZ than not seizure patients (see Fig. 4B). Also, there was weaker Theta in SOZ than NSOZ of seizure-free and not seizure-free patients. Lastly, seizure-free patients had weaker Theta in lateral temporal lobe NSOZ than not seizure-free patients (see Fig. 4C). These results derive from a seizure-free group that included patients without and with aura (i.e., Engel IA and IB outcomes). When the same analysis was performed with a seizure-free cohort consisting of Engel IA only (n = 8 patients) there was no difference in connectivity between seizure-free and not seizure patients.

Interictal spikes and connectivity
Previous studies found interictal spikes could affect connectivity, especially in the gamma band (Ren et al., 2015;Lagarde et al., 2018), and for this reason we included the rate of spikes in the model. The current analysis found a significant, albeit small, effect of interictal spikes on the strength of HGEC, LGEC, and ThEC (see Table 4). Consistent with the small estimated coefficients, overall analysis found a higher rate of spikes was weakly to moderately correlated with HGEC, LGEC, and Theta (r = 0.19, 0.29, and 0.27, respectively; see Fig. 5A for example of HGEC). At the level of the individual patient, 4 of the 43 patients had a strong correlation between interictal spikes and HGEC (r . 0.5; Fig. 5C, top scatterplot), but in all others, it was moderate (r , 0.5) or weak (r , 0.25; Fig.  5C, bottom scatterplot). The modest correlation between connectivity and spikes was unexpected and could be because of a limited sample of spikes in short duration recordings. However, there was not a significant correlation between individual r values of HGEC and spike rates (r = 0.315, p = 0.0576), suggesting a limited sampling of spikes alone can explain the modest correlation. Similar results were found with LGEC and Theta.

Interelectrode distance and connectivity
The distance between electrode contacts could affect connectivity strength, i.e., shorter distances correspond with stronger connectivity (Lagarde et al., 2018). The statistical model found interelectrode distance had a significant large effect on the strength of HGEC, LGEC, and ThEC (see Table 4). Overall, shorter distances between contacts correlated with stronger HGEC, LGEC, and Theta ( Fig. 4B, r = 0.45, 0.69, and 0.68 for HG, LG, and theta, respectively). In 42 out of 43 patients, there was a strong correlation between interelectrode distance and strength of connectivity (see example HGEC in Fig. 5D, top scatterplot), which is consistent with large estimation coefficients, and only one patient with a weak correlation (Fig. 5D, bottom scatterplot).

Connectivity in relation severity and duration of epilepsy
Difference in history or severity of epilepsy could affect connectivity; thus, we performed correlation analysis between strength of connectivity and measures of epilepsy severity and burden. Analysis found no correlation between strength of connectivity and duration of epilepsy (P durÀgec = 0.15, 0.028, 0.58), seizure frequency (P seizureFreqÀgec = 0.59, 0.99, 0.47), age of epilepsy onset (P onsetÀgec = 0.81, 0.84, 0.25), or burden of disease, i.e., seizure frequency/year X duration of epilepsy in years (P burdenÀgec = 0.71, 0.37, 0.72; Fig. 7). Also, there were no differences in the strength of connectivity between patients who received a resection and those who received an RNS, or between patients   with MRI lesion and those without a lesion (see Extended Data Fig. 7-1).

Discussion
The main findings in this study are (1) stronger HGEC and LGEC in SOZ than NSOZ of seizure-free patients; (2) stronger HGEC and LGEC between mesial and lateral temporal SOZ in seizure-free than not seizure-free patients; and (3) stronger LGEC and ThEC in extratemporal and lateral temporal NSOZ of not seizure-free than seizure-free patients. These results were unrelated to interictal spikes, clinical features of epilepsy, or MRI abnormality but were affected by interelectrode distance, which was adjusted for in the analysis. These relative differences in interictal event connectivity could indicate abnormal synchrony within and beyond the SOZ that contributes to seizure recurrence.

Differential event connectivity with respect to SOZ and NSOZ
Studies of functional connectivity in epilepsy commonly use linear or nonlinear correlation to assess the dependency between bandpass filtered EEG signals recorded from pairs of scalp or intracranial electrodes. Several Figure 2. High gamma event coupling in the SOZ and different brain regions. A, Examples of connectivity matrices of high gamma event coupling (HGEC) for patient 13 who was seizure-free (SF) and patient 19 who was not seizure-free (not-SF). The matrices are organized with respect to SOZ. If both electrode contacts are in SOZ, then the connectivity value is part of the SOZ, if both contacts are outside SOZ, then it is part of the NSOZ (complement), otherwise it is between the SOZ and NSOZ. The lower row illustrates HGEC organized by brain region (M: mesial temporal, L: lateral temporal, E: extratemporal). B, Violin plot and box plot (inside) shows the distribution, median and interquartile range of HGEC values for patients 13 and 19 with respect to SOZ (top rows) f and brain regions (bottom rows). In most cases, HGEC is stronger in patient 13 than patient 19. Check Extended Data Figure 2-1 for GEC matrices of all patients. studies found stronger interictal functional connectivity in the mesial temporal or extratemporal lobe SOZ than NSOZ (Bettus et al., 2008;Bartolomei et al., 2013;Lagarde et al., 2018). Stronger connectivity was found in conventional EEG frequency bands, including gamma, which is consistent with evidence of increased gamma power in the SOZ (Worrell et al., 2004;Medvedev et al., 2011;Cimbalnik et al., 2018;Zweiphenning et al., 2019). In the current study, we computed a form of connectivity using perievent time histograms to quantify the correlation between local maxima of individual events recorded from pairs of depth electrode contacts; a method previously used to assess event connectivity in rats (Kheiri et al., 2013;Bragin et al., 2014). With this approach we, too,    found stronger LGEC and HGEC in the SOZ than NSOZ, chiefly between the mesial and lateral temporal SOZ in seizure-free patients. Furthermore, we found stronger ThEC in NSOZ than SOZ, especially in lateral temporal lobe of not seizure-free than seizure-free patients, which could be related to the reduced theta power in mesial temporal than extratemporal lobe SOZ (Bettus et al., 2008). Differences in event connectivity associated with lateral temporal lobe found in our analysis are consistent with this region's involvement in some forms of temporal lobe epilepsy, especially those where the SOZ includes entorhinal cortex and MRI is normal or contains a lesion other than hippocampal sclerosis (Bartolomei et al., 2010), which characterizes many of the patients in the current study. To better understand the implications of these results to the seizure network and seizure outcome, it would be helpful to first explain what we believe event connectivity represents, which we discuss in the following paragraph.
What could event connectivity represent? Most brain rhythms like theta-band and gamma-band activity involve inhibition that can coordinate regular fluctuations in neuronal excitability, which generates coherent extracellular current flows measured in the EEG (Buzsáki and Watson, 2012). Gamma oscillations, for example, involve coordinated activity between inhibitory and excitatory cells (Buzsáki and Wang, 2012), but if there is inhibitory dysfunction, then there is greater excitatory asynchrony and increased gamma-band fluctuations (Yizhar et al., 2011;Cho et al., 2015). In the current study, it is likely LGEC and HGEC chiefly represent spontaneous gammaband fluctuations in multiunit activity, which was shown in rats (Bragin et al., 2014) and suggested to occur in humans (Burke et al., 2015). Regarding theta, which can be recorded in human mesial temporal lobe and neocortex (Kahana et al., 2001), it is possible that ThEC could correspond with Figure 5. Correlation between HGEC and interictal spike rate and electrode distance. A, B, Scatter plots illustrating HGEC in relation to spike rate (A) and electrode contact distance (B). Values are represented as normalize z scores. C, Specific examples of high (top) and low (bottom) correlation between spike rate and HGEC. The vertical bar to the right shows the percentage correlation coefficients for all patients. High r . 0.5 (shaded green), medium 0.25 , r , 0.5 (red), and low correlation r , 0.25 (blue). In most patients, the correlation between spike rate and HGEC was low. D, Same as panel C but correlation with electrode distance. In most patients, there was a high correlation between electrode distance and HGEC, i.e., as electrode distance decreases, HGEC increases. All correlations shown had a p , 0.0001. coordinated inhibitory and excitatory activity like LGEC and HGEC, but it involves a larger volume of tissue and/or greater spatial distribution of sources. Thus, we propose, in the epileptic brain, the strength of event connectivity corresponds with synchrony of inhibitory and excitatory activity such that relatively stronger event connectivity is associated with stronger synchrony and weaker event connectivity is associated with weaker or asynchronous inhibitory and excitatory activity.

Event connectivity and seizure recurrence
With this understanding of event connectivity, we interpret our results as follows. We assume that in seizure-free patients, brain regions corresponding with SOZ and NSOZ were completely identified, but in not seizure-free patients, the brain area responsible for generating seizures was incompletely identified and includes regions labeled SOZ and some in NSOZ (Fig. 8A). Prior work found increased excitability and synchrony in the SOZ (Staba et al., 2002;Schevon et al., 2007), and if this were because of deficits in inhibition, then it might be greater in not seizure-free than seizure-free patients to explain the recurrence of seizures. If this is correct and in the context of our current results, we should find stronger event connectivity in SOZ than NSOZ in seizure-free patients, which we do, and little difference between SOZ and NSOZ in not seizure-free patients, which also is consistent with our results. Furthermore, we should find weaker event connectivity, especially in NSOZ, of not seizure-free than seizure-free patients, but our results found stronger LGEC and ThEC in the NSOZ of not seizure-free patients. An alternative possibility is that rather than deficits in inhibition, there is a compensatory increase in the synchrony of inhibitory activity that is proportional to excitatory activity during interictal episodes (Fig. 8B). This explanation is more compatible with our results, particularly the stronger event connectivity in the NSOZ of not seizure-free patients and could correspond to increased synchrony of inhibitory and excitatory activity from an actual or potential SOZ (Lüders et al., 2006;Jehi, 2018).
Frequency band-specific sensitivity for the SOZ Analysis found more differences in g -band than thetaband event connectivity. One reason could be that unlike theta activity in rats (Buzsáki, 2002), the mechanisms generating theta are unclear in humans. However, like rats, theta can be recorded from several subcortical and cortical areas, which we suggested could correspond with large or distributed neuronal sources. It is possible that some of our recording contacts recorded neuronal activity from a common theta source that overlapped with SOZ and NSOZ making it less sensitive to detect differences between SOZ and NSOZ than LGEC and HGEC. Also, we computed event connectivity from low-gamma and highgamma bands like in previous rat studies ( Kheiri et al., 2013;Bragin et al., 2014) and as is often done in studies on gamma (Buzsáki and Watson, 2012). Although we Figure 6. Comparing event connectivity in three frequency bands. A, Violin plots show median Euclidean distance between pairs of contacts in relation to SOZ (abbreviation same as Fig. 2) and seizure outcome for all patients (abbreviation and shading same as Fig. 4). B, Violin and box plots of correlation coefficient between electrode distance and strength of functional coupling in high gammaðHg ), low gammaðLg Þ, and theta ðu Þ frequency bands in all zones and regions. Note that each patient has one correlation value, i.e., the violin plots are for 43 points each. C, The decay constant (t ) of the exponential decay model [EC = A*exp(-t *d)] relating the variation of event coupling strength (EC) of different frequency bands (Hg , Lg , u ) with the distance (d) between channels is illustrated in form of violin plots each representing 43 patients. See Extended Data Figure 6-1, which illustrates the difference between slow and fast decays. D, Coupling strength for Hg , Lg , and u are compared; p , 0.05 denoted by asterisks (*).
found similar results with LGEC and HGEC, they were not identical, and we plan future studies to investigate this further. Guiding this future work will be evidence suggesting that low gamma activity could involve inhibitory-inhibitory interactions and high gamma activity more dependent on inhibitoryexcitatory interactions (Kay, 2003). The potential differences in the contribution of local (inhibitory) and projection cells (excitatory) between gamma and theta might be related to the differences we found in the correlation between strength of event connectivity and distance. The strength of LGEC declined more rapidly with longer distances than HGEC or ThEC, which could be explained by greater contributions of local inhibitory cells in the former and more involvement of projecting excitatory cells with the latter.  . Relating connectivity to neuronal circuits mechanisms. A, A schematic illustrating brain regions involved in generating seizures (red dots) or those not involved (green dots). Clinically-defined seizure onset zone (SOZ; shaded orange) and not seizure onset zone (NSOZ; shaded blue). In an ideal seizure outcome, i.e., seizure free, all regions involved in generating a seizure are in the SOZ. The synchrony between brain regions is illustrated as connections (black lines), and a greater number of lines indicates greater synchrony. In not seizure-free patients, the SOZ is incompletely identified and a portion of the NSOZ contains regions involved in generating seizures. B, Prediction of the differences in the event connectivity when the strength of connectivity corresponds with increased synchronous inhibitory activity (blue dots and black lines) that is proportional to increased synchronous excitatory activity (red triangles and green lines). An assumption is greater synchrony associated with brain regions involved in generating seizures, which leads to the following predictions: (1) in seizure-free (SF) patients, stronger connectivity in SOZ than NSOZ; (2) in not seizure-free (NSF) patients, little or no difference in connectivity between SOZ and NSOZ; and (3) stronger connectivity in the NSOZ of NSF than SF patients. Results from our analysis are presented as three squares for each frequency band (theta = u , low gamma = LG, high gamma = HG), which are colored red and white for actual results that are consistent or inconsistent, respectively, with the aforementioned predictions.

Factors that could affect the strength of event connectivity
There are several factors to consider when interpreting the current results. First, it is important to distinguish between connectivity that derives from information theory and that based on amplitude or phase correlation/coherence methods. For example, in amplitude correlation, the power of the signal is an important factor that affects the strength of connectivity. As noted previously, event connectivity derives from the entropy of the peri-event histogram and is affected by the timing of the individual events. However, the algorithm detecting the peak of events requires an amplitude threshold, and it is possible it missed low amplitude events. Second, connectivity was computed from a 10-to 15-min interictal recording. Like previous studies, we selected a duration and time of recording to reduce potential effects of general anesthesia, spontaneous seizures, and anti-seizure medication tapering (Bartolomei et al., 2008a,b;Bettus et al., 2008;Medvedev et al., 2011;Cimbalnik et al., 2018;Lagarde et al., 2018;Klimes et al., 2019), and like these other studies, we found comparable results. Also, there is evidence that event connectivity is stable over a period of several days in freely behaving rats (Kheiri et al., 2013), and using the same methodology, we found 84.36 13.0% (n = 5 patients) similarity in the strength of connectivity between signals from first 10 min and the last 10 min of the recording. Third, an increase in neuronal spiking firing during interictal spikes can generate g activity (Alvarado- Rojas et al., 2013;Muldoon et al., 2015;Ren et al., 2015), which might overestimate the strength of connectivity on contacts with high rates of spikes. We included the rate of interictal spikes as a covariate in our linear mixed model and found spikes have a significant, but small, effect on connectivity. The latter result is consistent with other work that found little difference in connectivity values computed from EEG signals containing spikes and the same EEG signals after spikes were removed (Bettus et al., 2008). Fourth, we realize that Euclidean distance is an imprecise measure of anatomic connectivity, yet there was a significant effect of distance on events connectivity. Interelectrode distance was shorter in SOZ and NSOZ of seizure-free than not seizure-free patients, justifying the decision to include distance as a covariate in our linear mixed model. Same as for interictal spikes, connectivity values were adjusted for differences in distances and unlikely explain connectivity results with respect to seizure outcome. Lastly, measures of seizure severity, epilepsy burden, or other features of epilepsy did not correlate with the strength of connectivity after correcting for multiple comparisons, suggesting differences in connectivity with respect to SOZ and seizure outcome do not correspond with progressive aspects of epilepsy.
In conclusion, event connectivity is sensitive to differences in the synchrony of signals recorded in the SOZ and NSOZ and between surgical patients with and without seizure control. Differences in the strength of event connectivity between SOZ and NSOZ suggest a well-localized seizure network. By contrast, little or no difference in event connectivity could indicate a larger brain area generating seizures than localized to the SOZ and higher likelihood for seizure recurrence. In future work, we plan to perform unit recordings to investigate the neuronal basis of event connectivity and how changes in the strength of event connectivity correlate with neuronal excitability in brain areas where seizures begin and spread.