Resected Brain Tissue, Seizure Onset Zone and Quantitative EEG Measures: Towards Prediction of Post-Surgical Seizure Control

Background Epilepsy surgery is a potentially curative treatment option for pharmacoresistent patients. If non-invasive methods alone do not allow to delineate the epileptogenic brain areas the surgical candidates undergo long-term monitoring with intracranial EEG. Visual EEG analysis is then used to identify the seizure onset zone for targeted resection as a standard procedure. Methods Despite of its great potential to assess the epileptogenicty of brain tissue, quantitative EEG analysis has not yet found its way into routine clinical practice. To demonstrate that quantitative EEG may yield clinically highly relevant information we retrospectively investigated how post-operative seizure control is associated with four selected EEG measures evaluated in the resected brain tissue and the seizure onset zone. Importantly, the exact spatial location of the intracranial electrodes was determined by coregistration of pre-operative MRI and post-implantation CT and coregistration with post-resection MRI was used to delineate the extent of tissue resection. Using data-driven thresholding, quantitative EEG results were separated into normally contributing and salient channels. Results In patients with favorable post-surgical seizure control a significantly larger fraction of salient channels in three of the four quantitative EEG measures was resected than in patients with unfavorable outcome in terms of seizure control (median over the whole peri-ictal recordings). The same statistics revealed no association with post-operative seizure control when EEG channels contributing to the seizure onset zone were studied. Conclusions We conclude that quantitative EEG measures provide clinically relevant and objective markers of target tissue, which may be used to optimize epilepsy surgery. The finding that differentiation between favorable and unfavorable outcome was better for the fraction of salient values in the resected brain tissue than in the seizure onset zone is consistent with growing evidence that spatially extended networks might be more relevant for seizure generation, evolution and termination than a single highly localized brain region (i.e. a “focus”) where seizures start.

Using inter-ictal iEEG recordings of five epilepsy patients, Andrzejak et al. [31] have recently shown that EEG signals derived from the epileptogenic cortex are less random, more nonlinear-dependent and more stationary than those recorded from non-focal brain regions. Previously, similar findings have been demonstrated for the hemispheric, but not sub-lobar level [32][33][34]. However, despite being precondition for SOZ localization, a difference in average qEEG measures between epileptogenic and non-epileptogenic (or focal/non-focal or ipsilateral/contralateral to the SOZ) cortex alone is not sufficient for precise delineation of candidate tissue for surgical resection. For practical diagnostic applications the overlap of distributions of qEEG measures between epileptogenic and non-epileptogenic brain tissue is crucial. In the present study we set out to investigate the extent to which the prevalence of saliency of certain iEEG channels with respect to four qEEG measures at a given time was associated with resection of the corresponding brain tissue and with surgery outcome in terms of seizure control. The selected qEEG measures have previously been used in the Bern group for iEEG analysis and are representative of four different main categories of signal analysis without making claims for completeness: the absolute signal slope S, the number N of forbidden ordinal patterns, a surrogate corrected cross-correlation matrix C and a surrogate corrected mutual information matrix M.
The absolute value S of the first temporal derivative ("slope") of the signals represents a simple linear and univariate measure to quantify epileptiform EEG signals [35] and is-except for a scaling factor-identical to the "line length" feature [36]. The absolute slope S has for example been used to define objectively the onset and termination of seizures [37,38] as well as the spatial extent of secondary seizure generalization [39,40]. The so-called number (or fraction) of forbidden ordinal patterns N [41], another univariate but nonlinear qEEG measure, quantifies the degree of signal determinism. For qEEG analysis it was used in [42][43][44].
qEEG measures do not only allow to quantify features of individual signals (univariate measures), but also help to assess directed or non-directed signal interrelations and dependences (bivariate measures) and network properties (multivariate measures). As a non-directed, linear, bivariate qEEG measure we here use the "cross-correlation strength" matrix C [45,46] as a corrected version of Pearson's correlation matrix. This correction accounts for random effects due to the relative power in low frequencies and the limited duration of the moving window used for time-resolved signal analysis. A nonlinear, mutual information based analog M of the matrix C, where linear univariate properties as well as the effects of linear signal interrelations are compensated for, was introduced by Rummel et al. [46]. Here we apply the matrices C and M to a larger set of iEEG recordings for the first time.
The aim of our study is to test whether locally salient ("focal") EEG features as detected by the four qEEG measures X = {S,N,C,M} at any time step can be used to identify brain parts that should be included within the RBT to achieve seizure freedom. To this end, we compare the dynamics of the qEEG measures X immediately before, during and directly after seizures in patients with focal epilepsy syndromes. Specifically, focal EEG features that stand out from the background are defined in a dynamic and data-driven manner from the distribution of X across iEEG channels at each time step. To link the qEEG results with clinical findings we assess not only the sub-group of patients with favorable post-surgical outcome but as a contrast also unfavorable outcomes. An important aspect of our study is that by using image coregistration procedures, the anatomical positions of iEEG contacts are precisely localized and the corresponding iEEG signals are categorized into the zones Z = {RBT,SOZ,OVL,NON}. Here, OVL is defined as the overlap of RBT and SOZ. These channels also contribute to the RBT and the SOZ individually. The zone NON comprises all channels that neither recorded from the RBT nor from the SOZ.
We find that the four qEEG measures X reveal different properties in the different zones Z. For S, N and M salient values are overrepresented in the RBT of patients with post-surgical seizure control. In contrast, in those patients who have not been rendered seizure-free post-surgically, salient values are often localized outside the RBT. We conclude that part of the investigated qEEG measures provide clinically relevant and objective markers of target tissue for optimizing epilepsy surgery.

Patients
We included 16 patients of the epilepsy surgery program of the Inselspital Bern in this study (11 female, 5 male; median age 31.0y, IQR 15.3y, range 19-59y). Demographic data on the selected patients is compiled in Table 1. Six patients fell into Engel class I (free of disabling seizures), five into Engel class II (rare disabling seizures) and five into Engel class IV (no worthwhile improvement). Inclusion criteria were that post-surgical follow up was available for at least one year (median 3.0y, IQR 1.3y) and high resolution pre-and post-operative T1-weighted (T1w) MRI had been registered. To avoid introducing a bias towards patients with the largest number of seizures during intracranial EEG monitoring (median 4.5, IQR 5.5, range 2-14), only two seizures of each patient were included in the study (five patients from all Engel classes had only two seizures during long-term monitoring with iEEG, see Table 1). If a patient had several seizure types (pat. I-5, IV-3 and IV-4), two seizures of each type were included. For each seizure type we selected the first two seizures occurring during iEEG monitoring for qEEG analysis.
This study was approved by the Internal Review Board of the Inselspital (approval No. 159399, dated 26 th of November, 2013). All patients gave written informed consent that imaging and EEG data may be used for research purposes. The decision on the necessity for intracranial EEG diagnostics, the chosen electrode implantation scheme and the decision on surgical epilepsy therapy was made entirely on clinical grounds. These decisions were taken prior to and independently from the present retrospective study.

EEG acquisition
EEG signals were recorded intracranially by strip, grid and depth electrodes (all manufactured by AD-TECH, Wisconsin, USA), using a NicoletOneTM recording system with a C64 amplifier (VIASYS Healthcare Inc., Madison, Wisconsin, USA). An extracranial electrode, localized between 10-20 positions Fz and Cz, was used as reference for signal recording. EEG recordings were either sampled at 512 or 1024 Hz, depending on whether they were recorded with less or more than 64 contacts. The latter were down-sampled to 512 Hz prior to further analysis, and EEG signals were re-referenced against the median of all the channels free of permanent artifacts as judged by visual inspection. In addition to anti-aliasing filtering needed for proper sampling and down-sampling, EEG signals were digitally band-pass filtered between 0.5 and 150 Hz using a fourth-order Butterworth filter prior to analysis. Forward and backward filtering was applied to minimize phase distortions. that were continuously corrupted by artifacts of any kind were identified and excluded from further analysis. The number of remaining iEEG channels is denoted as n (median 60.0, IQR 25.8, range 32-100, see Table 1). Second, the time points of seizure onset and termination were determined for all included seizures. From this information the following epochs were deduced: Three consecutive epochs of one minute duration immediately before seizure onset, three consecutive epochs of one minute duration immediately after seizure termination and six equally sized epochs during seizure evolution (three termed "early ictal" and three "late ictal"). Third, the iEEG channels were identified that recorded from the visually defined SOZ. Finally, the number of seizure types per patient was determined (different SOZ and/or different sequence of seizure states). Seizure onsets were visually identified as the time of earliest EEG change associated with seizures following a systematic approach as previously described by Litt et al. [47]. We first identified unequivocal seizure activity and then moved backward in time to the moment when the first sustained change of EEG relative to the background pattern occurred. Repetitive high amplitude spike-and-wave signals or a flattening followed by low voltage fast activity [48] were the most frequently observed changes. Unequivocal seizure onset was defined as any epileptiform signal pattern, which became clearly identifiable without knowing that a seizure followed. Seizure termination on the other hand was identified as the time point when clearly identifiable epileptiform signals vanished, which typically manifested as a sudden drop of amplitude and frequency.

Quantitative EEG analysis
Independent from visual EEG analysis for clinical purposes we analyzed peri-ictal epochs of the selected seizures with four qEEG measures. Quantitative signal analysis was done with C and MATLAB (MathWorks, Natick, MA., USA) programs developed by the authors.
Univariate linear qEEG measure: absolute signal slope S. The absolute value of the first temporal derivative of the signals ("slopes") was used to detect epileptiform activity in single iEEG channels in an objective way, as previously described in [37][38][39]. In brief, the slope of the iEEG signals was estimated by the differences where i and t denote the channel and the temporal sampling points and Δt is the sampling interval, respectively. The absolute value S of the signal slope defined in Eq 1 provides an appropriate characterization of epileptiform EEG, since it is large for both slow, high amplitude signals as well as for fast, low amplitude signals [35]. The absolute signal slopes S were averaged over T = 1024 sample points (2 seconds) and normalized by dividing by their channel-wise standard deviation during a pre-ictal reference period of 60 seconds duration starting three minutes before visual seizure onset. Univariate nonlinear qEEG measure: number of forbidden ordinal patterns N. We quantified signal determinism (as opposed to stochasticity) by the number of forbidden ordinal patterns N. To compute this quantity we followed an approach proposed by Amigó et al. [41] and recently applied to EEG data by [42][43][44]. A univariate iEEG signal x i (t) was mapped to a finite number of ordinal patterns by first choosing two parameters, a pattern order d>1 and a time delay τ!1. Then d observations x i (t) spaced at τ sampling intervals of length Δt were determined to generate an embedding vector of dimension d: The elements of the vector x ! i;t were mapped uniquely onto the permutation π = (π 0 , π 1 , . . ., π d-1 ) of (0, 1,. . ., d−1) that fulfilled x i ðt À p 0 tDtÞ x i ðt À p 1 tDtÞ x i ðt À p dÀ1 tDtÞ ð 3Þ Equal values were ordered according to the time of their appearance in x ! i;t . As an example, the 5-dimensional embedding vector x ! i;t = (1.26, 6.38, 0.63, 1.26, 4.92) was mapped onto the ordinal pattern π = (2, 4, 1, 0, 3). In d dimensions we have d! different permutations π in total and the sample size of ordinal patterns generated with maximal overlap from a time series of length T is T−(d−1)τ. Useful inequalities for the parameters d, τ and T in relation to the power spectrum of the signals have been derived in [44]. In the present paper we chose d = 5, τ = 1 and T = 1024 sample points (2 seconds) for the calculation of N. Normalization was performed in the same way as for the absolute signal slope S using the same pre-ictal reference epoch.
Multivariate linear qEEG measure: surrogate corrected cross-correlation matrix C. Pearson's equal-time (zero-lag) cross-correlation coefficient ρ is a linear measure for the dependence of two time series (or data sets in general). For infinitely long and independent time series of arbitrary power spectrum, ρ equals zero. However, time series of finite length-as typically result when using a moving window approach for time-resolved EEG analysis-may give rise to rather large random values of ρ, even when the time series are completely independent. In particular, such spurious correlations can occur when low frequencies (with regard to the short duration of the observed time series) dominate the power spectrum.
To estimate the strength of this random correlation we generated a set of n surr independent realizations of univariate iterated amplitude adjusted Fourier transform (IAAFT) surrogate time series [49]. In mathematical terms, the univariate IAAFT surrogates are used to estimate the width of the zero-centered distribution of Pearson coefficients ρ under the null hypothesis that (i) the time series are two independent stationary linear stochastic auto-correlated Gaussian processes. (ii) The measurement functions by which the signals were derived from the dynamics are invertible but potentially nonlinear. (iii) The auto-correlations, means, and variances of the underlying Gaussian processes are such that the measurements result in the autocorrelations and amplitude distributions of the observed time series.
For positive ρ a formula that compensates for random correlations can heuristically be written as [45,46]: Here, ρ is the Pearson coefficient of the original data and hρ surr i is the median of the values obtained from the set of surrogate time series. s is a significance factor that assumes the value one if the null hypothesis that hρ surr i has larger absolute value or is equal to ρ can be rejected and zero otherwise. Details on the generalization of Eq 4 to negative ρ can be found in [45,46]. Surrogate based baseline correction strategies similar to Eq 4 were also followed in [31][32][33][34][50][51][52][53][54].
In the present study we followed a moving window approach to analyze the non-random correlation pattern of the entire peri-ictal iEEG recordings in time-resolved manner. To this end windows of 4096 sample points length (i.e. 8 seconds) were shifted over the recording with 512 sample points displacement (i.e. 1 second). On each time step the "genuine cross-correlation strength matrix" C [45,46] was constructed element-wise from Pearson's correlation matrix with coefficients ρ. For each window and EEG channel n surr = 10 univariate IAAFT surrogates of T surr = 4096 sampling points length were generated independently. In generalization of the heuristic formula Eq 4 the correlation coefficients ρ and ρ surr were calculated from partially overlapping subsegments of T = 1024 sampling points length (i.e. 2 seconds). Within the windows these subsegments were distributed with minimal overlap to generate ensembles of size n ens = 10 for the original time series and n ens n surr = 100 for the surrogates. A non-parametric Mann-Whitney-Wilcoxon U-test was performed to assess the significance of different medians of ρ and ρ surr and to determine the significance factor s. To account for multiple comparisons (the n-dimensional correlation matrix has n Ã (n-1)/2 different channel combinations) Bonferroni corrections were applied independently on each time step. As a quantifier of the contributions of individual channels to the correlation pattern we calculated the node strength known from network analysis [55], i.e. the sum over the absolute values of all matrix elements connecting any iEEG channel with the other channels, and normalized it by dividing by the maximally possible value n-1.
Multivariate nonlinear qEEG measure: surrogate corrected mutual information matrix M. Mutual information μ is a model-free information theory based measure for interrelation and is in particular not restricted to Gaussianity or linear dependence. It quantifies the deviation of the observed joint distribution of signal amplitudes x i (t) and x j (t) from the product of the marginal distributions, which would imply statistical independence. In the present study mutual information between channel pairs was estimated from multivariate EEG time series x i (t) using the k-nearest-neighbor algorithm [56] with parameter k = 3 as implemented in the publicly available MILCA package (http://www.ucl.ac.uk/ion/departments/sobell/Research/ RLemon/MILCA/MILCA). Especially for short (T 1000) and noisy time series this algorithm has been shown to be superior to other estimators [57]. Using a transformation suggested by Joe [58] μ can be normalized to the interval [0,1]. For positively correlated linear stochastic and Gaussian distributed data this transformation in addition warrants identity of μ with ρ. If the restrictions are relaxed, the Pearson coefficient represents a lower bound for normalized mutual information [59].
To correct for the influence of the linear ρ on the more general μ a very similar baseline correction as the one described in detail above was followed [46]. It differed only by using multivariate rather than univariate IAAFT surrogates [49] to sample the distribution of μ and by replacing ρ by μ and C by M in the heuristic formula Eq 4. Multivariate IAAFT surrogates resample the observed time series to test the following null hypothesis. (i) the dynamics is a stationary multivariate linear stochastic auto-and cross-correlated Gaussian process. (ii) The measurement functions by which the signals were derived from the dynamics are invertible but potentially nonlinear. (iii) The auto-and cross-correlations, means and variances of the underlying Gaussian process are such that the measurements result in the auto-and cross-correlations and amplitude distributions of the observed time series. For practical implementation of multivariate IAAFT surrogates we used the same segment lengths and sample size parameters as for the univariate IAAFT surrogates. Like for the linear matrix C we used Bonferroni correction to account for multiple comparisons and the normalized node strength as channel-wise quantifier of nonlinear interrelation.

MR and CT image acquisition
All MRI scans were performed on a 3T Siemens Magnetom Trio (Erlangen, Germany) equipped with a 12-channel radio frequency head coil. High-resolution T1w MR images were obtained with a 3D Modified Driven Equilibrium Fourier Transform (MDEFT) [60] sequence. The optimized acquisition parameters included: 256 × 224 x 176 matrix points with a noncubic field of view of 256 mm × 224 mm x 176 mm, yielding a nominal isotropic resolution of 1 mm 3 (i.e. 1mm × 1mm × 1mm), repetition time TR = 7.92 ms, echo time TE = 2.48 ms, flip angle = 16°, inversion with symmetric timing (inversion time 910 ms), fat saturation, 12 minutes total acquisition time.
Within 24 hours after electrode implantation, every patient underwent a CT scan to control electrode localization and rule out subdural or intracranial hematoma. Following the safety protocols at our institution post-implantation MRI were not acquired. A follow up MRI using the same sequence and parameters as before electrode implantation was acquired three to four months after resective surgery.
Image processing and coregistration. First, T1w pre-and post-operative MR images as well as CT images were coregistered (without reslicing) to MNI-templates from each modality using SPM8 (http://www.fil.ion.ucl.ac.uk/spm/software/spm8/). The MNI-CT template was taken from the Clinical Toolbox [61] (http://www.nitrc.org/projects/clinicaltbx/). In a second step, CT and post-T1w images were coregistered to the pre-T1w image maximizing normalized mutual information [62], and resliced to 1 mm 3 isotropic resolution. The pre-T1w image was then segmented (in native space) into tissue probability maps of grey matter (GM), white matter (WM) and cerebrospinal fluid (CSF) using VBM8 (http://dbm.neuro.uni-jena.de/vbm/). The segmentation routine contains an additional step to calculate the partial volume estimates (PVE) of the main GM, WM and CSF as well as mixed GM-WM and GM-CSF tissue classes, following methods described in [63]. A PVE image was thus generated where each voxel was labeled according to the tissue class it belongs to. This labeled PVE image was thresholded at a value of 1 to generate a binary mask of the brain and surrounding CSF (since parts of the iEEG electrodes are located in the subdural space). We then applied this mask to the coregistered CT images and thresholded the resulting image at a threshold of 2,500 Hounsfield units. This procedure segregated the intracranial electrodes into an "electrode only"-image, essentially using the CT-artifacts generated by the implanted electrodes as a proxy for their location. Finally, the RBT was identified in the post-T1w images and manually segmented using MRIcron (http:// www.mccauslandcenter.sc.edu/mricro/mricron/), resulting in a binary RBT mask.
Three-dimensional renderings of the pre-T1w, electrode image and RBT mask were used in MRIcron to match each contact to the corresponding iEEG channel, and to determine their spatial relationship to the RBT. Channels recording from the RBT are indicated in red in all figures. Clinical reports of the intracranial recordings and previous visual analysis of the iEEG were used to identify iEEG channels belonging to the SOZ (blue in all figures). Channels that corresponded to the overlap (OVL) of RBT and SOZ are colored in magenta in all figures. Finally, 3D visualizations of the RBT, electrode image, and critical iEEG channels were generated using BrainNetViewer (http://www.nitrc.org/projects/bnv/) [64].

Statistics
To analyze focal EEG features as quantified by the four qEEG measures X = {S,N,C,M} we separated channels with saliently high values from those yielding normal and small values on the same time step. The separation, i.e. the "saliency", was obtained by determining the first and third quartile (Q 1 and Q 3 ) as well as the inter-quartile range IQR = Q 3 -Q 1 of the distribution of X separately on each time step. All channels with values larger than Q 3 +w Ã IQR with w = 1.5 were defined as salient channels. The advantage of this dynamic and data-driven approach is that minimal assumptions have to be made about the shape of the underlying distributions, which may even change from time step to time step. Such an approach is well known from the definition of "outliers" at the upper end of box-and-whisker-plots and is applicable also for very skew distributions. In addition, w is the only parameter, for which we use the standard choice for outlier definition. Importantly, no assumption is made about the prevalence of salient values indicated by measure X = {S,N,C,M}. As they need to be larger than Q 3 , their number satisfies 0 n X n/4 in a sample of size n and crucially depends on the shape of the distribution of X.
The fraction of focal salient EEG channels indicated by measure X found inside one of the channel zones of interest Z = {RBT,SOZ,OVL,NON} with size 0 n Z n was calculated as the relative size of the channel intersection F X Z ¼ n Z\X =n X . The values of F X Z are dependent on the number of salient channels n X , which typically varies as a function of time in the peri-ictal evolution, as well as on the spatial extent n Z of the zones Z, which varies from patient to patient (RBT) or may even vary from seizure to seizure (SOZ, OVL and NON). As both n X and n Z may independently become zero (and in consequence also n Z\X = 0), a simple normalization of F X Z to the relative size n Z /n of zone Z may be ill defined. Thus, we complementarily assessed the significance of the fraction of salient values within zone Z using score statistics. The probability to find exactly ν salient values inside Z and n X − ν outside Z (no replacements allowed) is given by the probability density function of the hypergeometric distribution [65]: Values of ν are restricted to the interval between ν min = max(0,n Z + n X -n) and ν max = min(n Z , n X ), with a positive lower limit in conditions with more salient channels of measure X than channels outside zone Z (n X > nn Z ). The cumulative probability of finding m or more salient values inside zone Z yields the significance of a hypergeometric test for over-representation: We report our results as the negative logarithms of Eq 6 with large values indicating high significance of the observation: L X Z ¼ Àlog 10 ðP n X ;n Z Þ. All statistics were performed on significance level α = 0.05. As this is an explorative study in a small group of patients, no correction for multiple comparisons (four measures X, three zones Z, three Engel classes) was performed for group comparison. For comparison of categorial data (gender, syndrome etc.) the χ 2 -statistics was used. As the χ 2 -test is strictly applicable only if the smallest expected bin filling is five or larger, we resampled the group allocation randomly 10,000 times to estimate the p-values. To assess differences in the distribution of ordinal quantifiers (age, qEEG measures etc.) over several outcome classes the non-parametric Kruskal-Wallis test was applied. Conditional on the test result we performed post-hoc pair-wise comparisons with a Mann-Whitney-Wilcoxon U-test to check which classes indeed differed.

Demographics and clinical information
Patient information is summarized in Tables 1 and 2. At a significance level of α = 0.05, none of the quantifiers differed between the three outcome classes. Trends were only observed towards different fraction of iEEG channels in the overlap OVL between of RBT and SOZ (p = 0.076) and in their union (p = 0.073).

Peri-ictal qEEG analysis
We illustrate the temporal evolution of the qEEG measures using as an example the surrogate corrected mutual information matrix M for the first seizure of patients I-2 and IV-1 in Figs 1 to 4. Results for all four measures and both included seizures of these two patients are compiled in S1-S4 Figs. These figures clearly demonstrate that the peri-ictal dynamics of most qEEG measures is highly stereotypical, i.e. it is very similar for different seizures of the same type in an individual patient. Patient I-2 remained seizure free after surgery during the follow-up period of three years. Fig 1 shows the anatomical relationship of the RBT and intracranial contacts (panels a and b), as well as the peri-ictal mean contribution of each electrode to the qEEG measure M (panel c). All n = 64 iEEG channels were free of permanent artifacts and therefore were all included in the analysis. The overlap OVL between the SOZ (11 channels) and the RBT (13 channels) was 9 channels (set union n-NON = 15 channels, Jaccard index J = OVL/(n-NON) = 0.6). The periictal evolution of the normalized node strength of the surrogate corrected mutual information matrix M during the first seizure recorded is displayed in panel t of Fig 2. Ahead of the seizure, channels DEL01 to DEL05 of the left mid-temporal depth electrode revealed the largest node strength, followed by the left temporo-polar channel TPL01 and the right temporo-polar channel TPR02. The channels DEL01 to DEL05 and TPL01 were all in the overlap of the SOZ and the RBT, whereas channel TPR02 was recorded from brain tissue that was not resected. During the first 2/3 of the seizure the average normalized node strength of M decreased, before a reincrease occurred during the last 1/3 of the seizure. Interestingly, the node strength pattern before seizure termination was much more uniformly distributed than the one observed at earlier times. After seizure termination the strongest contributions to the normalized node strength of M were due to different channels than before the seizure: TPL01 to TPL04 (all in the overlap of SOZ and RBT) and the temporo-basal channel TBL01, which was in the RBT but not in the SOZ. Bar plots of temporal means over the pre-ictal, early ictal (first half of the seizure), late ictal (second half of the seizure) and post-ictal phase are shown in panels p, q, r and s, respectively.
The temporal evolution of saliency in the normalized node strength of the mutual information matrix M is shown in binarized form in panel o (white: salient channel, black: non-salient channel) and bar plots of temporal means over the four phases are given in panels k, l, m and n. In the pre-ictal as well as in the early and late ictal phases the saliency pattern was very stable over time, with largest frequency of salient values in the overlap OVL of RBT and SOZ (channels DEL01 to DEL06 and TPL01 to TPL04). Ahead of seizure termination the amount of salient channels decreased and concentrated on channels recording from the right hemisphere (contralateral to the visually defined SOZ) which were consistently neither in the SOZ nor in the RBT. After seizure termination, the main prevalence of salient values was in channels TPL01 to TPL04 and TBL01. All of them were in the RBT and most of them in the SOZ. Nevertheless, the post-ictal strength pattern was substantially different from the pre-ictal and early ictal pattern. Peri-ictal evolution of quantifiers derived from the normalized node strength of the surrogate corrected mutual information matrix M of iEEG signals. The seizure starts at time point zero. Visually determined seizure onset and termination are indicated by vertical lines in the panels e and j. iEEG channels recording from RBT are indicated in red, the SOZ in blue and the overlap of both in magenta on the very left. Channels belonging to none of these zones are indicated in black. The color-scale figure in panel t shows the temporal evolution of each iEEG channel's normalized node strength of M. The vertical bar plots in panels p, q, r and s display the mean channel contribution in each of the following phases: pre-ictal, early ictal, late ictal and post-ictal. All these bar plots are scaled for optimal display. The bars are light shaded, whereas the standard error of the mean contribution is displayed in full color. The temporal evolution of salient (white) and normal channels (black) is shown in panel o. The vertical bar plots in panels k, l, m and n show the channel-wise mean prevalence of salient values in the four peri-ictal phases. For better comparison, here, all bar plots are shown in the same range 0 to 1.1. Panel e shows the temporal evolution of the fraction of salient values falling into the RBT (red) and the SOZ (blue). Panels a, b, c and d show the means over all four zones Z = {RBT,SOZ,OVL,NON} during the four peri-ictal phases. The 95% confidence interval for the mean is displayed as whiskers. The negative logarithm of the probability for randomly finding the observed or a larger amount of salient channels in the RBT (red), the SOZ (blue) and the overlap OVL of both (magenta) is shown in panel j as a function of time and in panels f, g, h and i as mean over the four peri-ictal phases. The 95% confidence interval for the mean is displayed as whiskers and the horizontal line indicates an approximation to the significance threshold -log 10 (0.05) = 1.30 for temporal means. The fraction of salient channels of the node strength of M falling into each of the zones RBT and SOZ varied between 0.5 and 1 pre-ictally (panel e). None of the fractions in the RBT, the SOZ or in OVL was distinguished, whereas hardly any salient channel fell into NON (panel a). The configuration changed substantially during the seizure: The fraction of salient values F M Z in the three zones RBT, SOZ and OVL decreased during seizure evolution and F M SOZ became larger than F M RBT and F M OVL . In contrast, the fraction of salient values F M NON in the channels belonging neither to the RBT nor to the SOZ increased (panels b and c). After seizure termination the fraction of salient values in the zones RBT, SOZ and OVL was unstable, but on average remained smaller than in the pre-ictal and early ictal phase. Now, F M Z had a similar value in the zones NON and RBT and was larger than in the zones SOZ and OVL (panel d).
The temporal evolution of the log probability L M Z is displayed in panel j. Ahead of the seizure this value varied for between 2 and 6 for zones Z = {RBT,SOZ}, equivalent to p-values in the range between 0.01 and 10 −6 and thus indicating a non-random association. Also on average L M Z was larger than -log 10 (0.05) = 1.30 (panel f), a value that can heuristically be seen as the significance threshold for temporal averages. None of the zones Z = {RBT,SOZ,OVL} had a clearly larger value than any other. During seizure L M Z increased up to the value 9 for the SOZ (overlap significance P = 10 −9 ). After the seizure L M Z dropped and varied around the heuristic significance threshold. The largest post-ictal values were obtained for L M RBT . In contrast to patient I-2 surgery outcome was unfavorable in patient IV-1 (Fig 3). Out of 62 implanted iEEG contacts, n = 59 channels recorded signals free of permanent artifacts. As the visually defined SOZ overlapped with the functional language area, it was decided to restrict palliative resection to the vicinity of the SOZ to minimize the risk of post-operative neurological deficits. In consequence, the overlap OVL between the visually defined SOZ (4 channels) and the RBT (2 channels) was empty. Retrospective qEEG analysis showed that in patient IV-1 large normalized node strength of the surrogate corrected mutual information matrix M were much more confined in space and time (panel t of Fig 4) than for patient I-2. Before and after the seizure salient values were mainly located in the left temporo-latero-posterio-basal iEEG channels TLPBL1 to TLPBL6, some of which were in the RBT and some of which in the SOZ (panels k and n). During the seizure, salient values were also visible on channels of the left temporo-latero-polar (TLPL) and left temporo-basal-lateral (TBLL) strip electrodes (panels l and m), which neither corresponded to the SOZ nor to the RBT. Throughout the peri-ictal

Summary statistics for qEEG measures
In Figs 5 and 6 we illustrate the peri-ictal evolution of our qEEG quantifiers F X RBT and L X RBT for all 38 seizures of all 16 patients. These radar plots enable an intuitive clockwise visualization of the temporal evolution before (upper right quadrant), during (lower half) and after seizures (upper left quadrant). Analog representations for the SOZ, its overlap OVL with the RBT and channels contributing to none of these zones NON are given in S5-S10 Figs.
The class-wise mean fraction F X RBT of salient channels that were resected during epilepsy surgery is shown in Fig 5 together with the value f expected from the size of the RBT under uniform distribution (fully drawn circles). The class-averaged peri-ictal mean <F> (proportional to the square root of the area included in the polygons) was largest for X = M (last column), followed by S (first column), N (second column) and C (third column). Outcome class II showed the largest resected fraction of salient values for all qEEG measures (second row), followed by class I (first row) and Engel class IV (third row). Note that in patients of Engel class II also the absolute and relative size of the RBT was largest though the difference lacked significance ( Table 2). For the unfavorable outcome class IV the average resected fraction of salient values was hardly ever larger than 25% for any qEEG measure in the peri-ictal evolution. In contrast, for the favorable outcome classes I and II we found time periods where the class-wise mean F X RBT was close to 50% despite the fact that the RBT comprised less than 25% of the channels on average. In patients who became completely seizure free (Engel class I) roughly half of the salient channels of the surrogate corrected mutual information matrix M found during the pre- . The data is arranged clockwise with the upper right quarter corresponding to three consecutive epochs of one minute duration immediately before seizure onset. The lower half corresponds to the scaled seizure time between seizure onset and termination, and the upper left quarter corresponds to three consecutive one-minute epochs immediately after seizure termination. The red polygons illustrate the temporal profile of the mean F X RBT and the parametrically estimated 95% confidence interval of the mean is displayed by the thickness of the polygon outline. Broken circular lines correspond to fractions 0.25, 0.5, 0.75 and 1.0 from inside to outside. The continuous circular line corresponds to the fraction f of salient channels expected from the size of the RBT. The mean <F> as estimated from the circle area is given in the lower right corner of each plot. Note that the means given here can deviate from the medians used in the text. Analog figures for the visually defined seizure onset zone (SOZ), its overlap OVL with the RBT and channels contributing to none of these zones NON are given in S5, S6 and S7 Figs. For the fraction of salient values in the SOZ and in the overlap OVL between SOZ and RBT we found similar results (S5 and S6 Figs) with the difference that large values of F X SOZ were more peaked in the early ictal phase, especially for X = S. Channels that neither recorded from the RBT nor form the SOZ had a consistently smaller fraction of salient values than expected from the size of zone NON only for qEEG measure M in the full peri-ictal epoch and for S and N in the early ictal phase (S7 Fig). With exception of the pre-ictal and early ictal phases of M and the early ictal phase of S in Engel classes I and II the class-averaged peri-ictal mean of L X RBT (polygon areas in Fig 6) hardly ever during the peri-ictal evolution reached the heuristic border of significance (L>1.3 or P<0.05, fully drawn lines). In terms of qEEG measures it decreased in exactly the same sequence as F X RBT , i.e. M>S>N>C. The number of forbidden patterns N and the surrogate corrected cross-correlation matrix C (second and third column, respectively) had much smaller polygon areas, indicating on average smaller significance of the overlap between salient channels of these two measures and the RBT. In terms of outcome classes, the polygon areas decreased from Engel class II (second row) over I (first row) to IV (third row).
The association between salient and SOZ channels became significant in the early ictal phase for qEEG measure S and in the pre-ictal phase for M (S8 Fig). To lesser extent a similar observation was made for the overlap OVL between RBT and SOZ (S9 Fig). In contrast, no association was found between salient channels and those neither in the RBT nor in the SOZ (S10 Fig). For all Engel classes the salient channels of the absolute signal slope S were stronger associated with the SOZ in the early ictal phase than with the RBT (panels a, e and i of S8 Fig  and Fig 6).
We estimated the peri-ictal means of the fraction of salient values F X Z of measure X in zone Z and the corresponding log probability L X Z seizure-wise from the area included in the radar plots. p-values of a Kruskal-Wallis test for class-wise different rank sum of F X Z (Table 3) were significant only for the RBT and qEEG measures M, S and N. For all these measures post-hoc testing for pair-wise differences showed that the medians in Engel classes I and II were significantly larger than in class IV, whereas there was no significant difference between classes I and II (S: 25.8%, 34.6%, 19.9%; N: 27.0%, 31.5%, 11.1%; M: 29.1%, 32.7%, 13.9% for classes I, II and IV, respectively). In contrast, for the SOZ and the overlap OVL significance was never reached. Table 4 shows the test results for class differences of L X Z . Here significant class differences were found for measures M and N in zone NON. For measure M the Engel class II had significantly larger median log probability than classes I and IV (0.096, 0.174, 0.061) and Engel classes I and II had significantly larger median log probability in measure N than class IV (0.110, 0.092, 0.064). Note that none of the log probabilities leading to class differences came even close to the heuristic significance threshold L>1.3 for means.  Statistics and result representation is as in Table 3. doi:10.1371/journal.pone.0141023.t004

Summary and Discussion
We have analyzed the spatial relation between iEEG channels with focally salient values in four qEEG measures and the resected brain tissue (RBT) as well as the seizure onset zone (SOZ). The exact localization of the RBT and of the electrodes recording the iEEG signals was confirmed radiologically by coregistration of pre-and post-surgical MRI and CT imaging. In contrast to this quantitative approach, the SOZ as well as seizure onset and termination times were defined by the current clinical gold standard, i.e. visual EEG reading by experienced epileptologists. The qEEG measures were selected to be representative for four different classes of signal analysis methods: linear univariate (S), nonlinear univariate (N), linear multivariate (C) and nonlinear multivariate (M). The methodology was retrospectively applied to 38 peri-ictal iEEG epochs from 16 epilepsy patients (three patients had two seizure types) with favorable (6 patients in Engel class I and 5 patients in class II) and unfavorable (5 patients in Engel class IV) post-surgical seizure control.

Main findings
Our main findings are the following. First, regardless of the used qEEG measure, the median fraction of salient channels in the RBT was larger in the favorable outcome classes Engel I and II than in the unfavorable outcome class IV, where on average less than 20% of the channels showing salient values were resected (Fig 5). For X = {M,S,N} the fraction of salient values in the RBT was significantly larger in the favorable than in the unfavorable outcome classes, whereas in the SOZ or in the overlap OVL we did not find a class difference. Second, in all outcome classes the average fraction of salient values in the RBT was larger for the absolute signal slope S and the normalized node strength of the surrogate corrected mutual information matrix M than for the number of forbidden ordinal patterns N and the node strength of the surrogate corrected cross-correlation matrix C (Fig 5). Third, the association between iEEG channels defined as salient and three out of four of the studied zones Z = {RBT,SOZ,OVL} intermittently reached beyond than chance level only for the absolute signal slope S and the node strength of matrix M and only for the favorable outcome classes I and II (Fig 6). Association between salient iEEG channels and the zone NON was never significant. Despite larger median values in favorable than in unfavorable outcome groups, L X RBT and F X RBT was often larger in Engel class II than I. The 95% confidence intervals of the means were also wider in class II. This rather counter-intuitive observation may be explained by two reasons. First, our patients with outcome class II could form a less homogeneous group than those of classes I and IV. Second and more likely, the fraction of iEEG channels in the overlap between SOZ and RBT was largest in class II (median fraction of 3.4%, 3.9% and 0.8% in Engel classes I, II and IV, respectively). This might have introduced a bias towards larger L X Z and F X Z in this subgroup. Note that despite different channel fraction in the overlap we did not observe a significant class difference in the Jaccard index for overlap between RBT and SOZ ( Table 2). The reason is that also the union of channels recording from the RBT and the SOZ represented a larger channel fraction in outcome class II than in IV. We interpret this finding in the sense that post-surgical seizure control tended to be better for patients where a larger fraction of iEEG channels recorded from the hypothetical EZ, i.e. for those patients with better preimplantation hypotheses about the localization of the epileptogenic brain areas.
Apart from these trends we did not find any statistical differences between the clinical parameters of the three outcome groups. We interpret this as corroborating evidence that the above summarized results are not an artifact of patient selection but underline the value of quantitative analysis of peri-ictal iEEG recordings added to classical expert visual inspection.

Discussion
For our study we selected qEEG measures that represent four classes of signal analysis algorithms. It is plausible that alternative ways of assessing similar signal properties would yield similar results. The absolute signal slope S [35] is a simple univariate and linear, but of course not the only qEEG measure that has been designed to detect epileptiform signals by its high frequency content. As alternatives the "joint sign periodogram event characterization transform algorithm" (JSPECT) [66], the "epilepticity index" [16] or high frequency oscillations [18] are conceivable. Alternatives to the univariate but nonlinear number of forbidden ordinal patterns N [41] for quantification of signal determinism could be the amplitude-based "nonlinear prediction error" [67,68] or its rank-based extension, the "nonlinear prediction score" [69,54]. Related are attempts to determine approximations to the correlation integral from EEG time series [70,71]. Also for quantification of linear or nonlinear signal interrelation there are many alternatives to the chosen surrogate corrected measures C and M of [46], see [9,10] for reviews. Without completeness we mention linear cross-coherence [72], the nonlinear correlation coefficient [73,74], phase synchronization [75], Granger causality [76] and cross-predictabilities [77].
Applying qEEG measures and identifying the RBT radiologically we centered our study around quantifiable information. We tried to avoid any potential ambiguity or inter-rater variability as might occur during visual EEG or MRI interpretation. The proposed method might be especially important in the subgroup of patients, where the current clinical gold standards, i.e. EEG "reading" by trained experts and visually defining the SOZ remain ambiguous. The safety of high-field MRI with implanted electrodes is still a matter of debate, thus co-registration of pre-and post-surgical MRI with post-implantation CT is currently the method of choice to define the precise position of intracranial electrodes and the RBT.
Our finding that the differentiation between favorable and unfavorable outcome was possible for F X RBT but not for F X SOZ is consistent with growing evidence that widespread networks might be more relevant for seizure generation, evolution and termination than a single brain region (i.e. the "focus") where seizures actually start from [78,22]. In a recent study [79] it was found that despite larger average extent of SOZ resection in pediatric patients with favorable outcome (2/3 resected as compared to 1/3 resected in unfavorable outcomes), complete resection of the SOZ was required only in one out of eight cases to achieve seizure freedom.
The dynamic and data-driven definition of salient qEEG values has the advantage that these prominent values can be analyzed whenever they occur. Specifically, this implies that no assumption has to be made when focal saliency might be most prevalent (e.g. in the early ictal phase, where "early" has to be specified somehow). This concept has the additional advantage that inter-ictal or sub-clinical events that might escape the observer's attention can be detected objectively and can be used for analysis.
Reporting F X Z and L X Z is better suited to situations without a "ground truth" than the use of classification accuracies like sensitivity and specificity, positive and negative predictive values, or likelihood and odds ratios [80]. In epilepsy surgery there is not yet a method to simulate different surgical procedures in order to evaluate if they are effective. Moreover, in successful cases the actual RBT is most probably larger than minimally needed (i.e. the epileptogenic zone EZ), which negatively biases most accuracy quantifiers. The fraction of salient values, however, has the property to saturate at 100% even if the RBT is larger than the EZ.

Limitations and Outlook
Our study has several limitations. First, the number of patients included into the study is limited (38 seizures from 16 patients), retrospective in nature, derived from a single center epilepsy surgery program and heterogeneous in terms of etiology and electrode implantation scheme. Although we were able to show that our results are robust against a number of confounding demographic and disease-related variables, it is unclear at present how our results will generalize to larger or prospective cohorts.
Second, the chosen scheme for detection of salient values is suitable for skewed distributions, but it is still a parametric concept. A different choice of the whisker parameter w produces a different number of salient values. We used the standard choice w = 1.5 and did not investigate the influence of different values of w on our results. As a further technical aspect, the IAAFT surrogate generation is computationally expensive, especially in the multivariate case where phase randomization is performed under the constraint of conserved phase relationships between all iEEG channels. For larger studies or prospective application of the surrogate corrected mutual information matrix M it has to be investigated whether faster algorithms, e.g. constrained phase randomization without amplitude adaption or the iterative procedure [81], are sufficient for our purposes.
In this study we applied the Engel seizure outcome classification, being aware of its limitations due to the dependency on the patient's and his/her family's perception and long-term reports of seizure rate, which must be considered with caution [82].
Finally, one inherent difficulty in this type of study-also mentioned above-is that even the RBT, as delineated by MRI, may be considered as an imperfect benchmark, since we cannot exclude that any other surgical approach might have led to a similar post-surgical outcome in a given patient, i.e. in a situation where the EZ is unique. In our opinion it will be hard to address or even resolve this issue satisfactorily. From a more pragmatic point of view, the RBT can easily and objectively be assessed with current imaging methods.
The clinical interpretation of some qEEG measures and results may be complex despite of their objective mathematical definition. Whereas the absolute signal slope S has a straightforward interpretation in terms of high frequency or large amplitude signals, the interpretation of the number of forbidden ordinal patterns N or of the surrogate corrected measures C and M is less obvious. Bringing qEEG closer to clinical understanding will be a task for future investigations.
At present we cannot prove that resection of brain tissue generating salient values in certain qEEG measures automatically leads to post-surgical seizure control. Rather, we believe that our findings suggest that special attention should be given to saliency-generating iEEG channels during pre-surgical evaluation and surgery planning. Primary candidates are salient channels of the absolute signal slope S in the early ictal phase or of the normalized node strength of the surrogate corrected mutual information matrix M in immediate pre-ictal epochs (Figs 5 and 6). Occasionally, it might be possible to adapt the intended resection target to increase the fraction of salient values in the anticipated RBT accordingly. The proposed methods may be further investigated in epilepsy surgery to estimate the expected outcome by simulation of the anticipated RBT. In situations where different surgical approaches are under debate, the alternatives could be modeled and could contribute to the decision making on the best strategy to render the patient seizure free. (TIFF) S1 Movie. 3D animation of the spatial layout of the intracranial electrodes implanted into the left hemisphere of patient I-2. Each sphere represents the position of one intracranial electrode. The sphere color corresponds to the analyzed zone, i.e. red, resected brain volume (RBT); blue, visually determined seizure onset zone (SOZ); magenta, overlap (OVL); black, neither of the above (NON). The sphere radius is proportional to the peri-ictal channel-wise mean of each qEEG measure during the first seizure (see main text for details). The layout is as in S1-S4 Figs: upper row, absolute signal slope S (left), number of forbidden ordinal patterns N (right); lower row, surrogate corrected cross-correlation matrix C (left), surrogate corrected mutual information matrix M (right). (AVI) S2 Movie. Same as S1 Movie but for the first seizure of patient IV-1.

Supporting Information
(AVI)