Identification of Preseizure States in Epilepsy: A Data-Driven Approach for Multichannel EEG Recordings

The retrospective identification of preseizure states usually bases on a time-resolved characterization of dynamical aspects of multichannel neurophysiologic recordings that can be assessed with measures from linear or non-linear time series analysis. This approach renders time profiles of a characterizing measure – so-called measure profiles – for different recording sites or combinations thereof. Various downstream evaluation techniques have been proposed to single out measure profiles that carry potential information about preseizure states. These techniques, however, rely on assumptions about seizure precursor dynamics that might not be generally valid or face the statistical problem of multiple testing. Addressing these issues, we have developed a method to preselect measure profiles that carry potential information about preseizure states, and to identify brain regions associated with seizure precursor dynamics. Our data-driven method is based on the ratio S of the global to local temporal variance of measure profiles. We evaluated its suitability by retrospectively analyzing long-lasting multichannel intracranial EEG recordings from 18 patients that included 133 focal onset seizures, using a bivariate measure for the strength of interactions. In 17/18 patients, we observed S to be significantly correlated with the predictive performance of measure profiles assessed retrospectively by means of receiver-operating-characteristic statistics. Predictive performance was higher for measure profiles preselected with S than for a manual selection using information about onset and spread of seizures. Across patients, highest predictive performance was not restricted to recordings from focal areas, thus supporting the notion of an extended epileptic network in which even distant brain regions contribute to seizure generation. We expect our method to provide further insight into the complex spatial and temporal aspects of the seizure generating process.

replaced by an on-demand therapy. This includes, among others, the application of fast-acting anticonvulsant substances or electrical stimulation (Theodore and Fisher, 2004;Morrell, 2006;Stacey and Litt, 2008).
Research in the field of seizure prediction over the last 15-20 years has identified analysis techniques that appear to be capable of detecting long-lasting spatial-temporal changes on intracranially recorded electroencephalograms (iEEG), which can be regarded as seizure precursors (see Mormann et al., 2007;Schelter et al., 2008 for comprehensive overviews). Despite the availability of continuous multiday, multichannel iEEG recordings, and statistical methods for testing the significance of the predictive performance of analysis techniques (Andrzejak et al., 2003(Andrzejak et al., , 2009Winterhalder et al., 2003;Kreuz et al., 2004;Schelter et al., 2006a;Wong et al., 2007;Snyder et al., 2008),

IntroductIon
Epilepsy affects 60 million humans worldwide, which is approximately 1% of the world's population (Hauser et al., 1996). Twothirds of epilepsy patients achieve effective seizure control from anticonvulsive medication, and another 8-10% benefit from resective surgery. Thus, presently in about 25% of epilepsy patients, seizures cannot be successfully controlled by any available therapy. The sudden, apparently unforeseen occurrence of seizures, which represents one of the most disabling aspects of epilepsy (Murray, 1993;Schulze-Bonhage et al., 2010), calls for a method that is capable of predicting the occurrence of seizures. This could significantly advance therapeutic possibilities (Elger, 2001). But also for those patients who can be treated successfully nowadays, current preventive treatment strategies could be no study has been published that demonstrates the identification of seizure precursors in prospective, randomized clinical trials with accuracy sufficient for clinical application .
One of the challenges in this field is the fact that a priori it is not obvious which regions in the brain have to be studied in order to reveal predictive information (Lehnertz, 2009). Dependent on clinical requirements, iEEG is recorded for prolonged periods from several tens of electrode contacts to identify potentially epileptogenic brain tissue and to delineate it from regions that are indispensable for defined cortical functions (Rosenow and Lüders, 2001). Nevertheless, for a retrospective identification of seizure precursors promising results have been achieved particularly with bivariate measures from linear or non-linear time series analysis. With a time-resolved analysis of multichannel iEEG recordings these techniques render time profiles of a characterizing measure, the so-called measure profiles.
Various downstream evaluation techniques have been proposed to single out measure profiles that carry potential information about preseizure states. These techniques either rely on heuristic topologic criteria like a restriction to neighboring electrode contacts (Mormann et al., 2003a(Mormann et al., ,b, 2005 or use a priori clinical information about individual recording sites (Maiwald et al., 2004;Esteller et al., 2005;Schelter et al., 2006a,b;Feldwisch-Drentrup et al., 2010). Some studies used a distinct part of the available iEEG data, e.g., from the beginning of the recording, to determine those measure profiles carrying optimal predictive information Le van Quyen et al., 2005), while others examined the predictive information of measure profiles for one seizure to select those for the prediction of a subsequent seizure (Iasemidis et al., 2001;Sackellares et al., 2006). Knowledge about the time of seizure onsets is required for the latter approaches.
A different approach bases on an evaluation of the significance of predictive information of measure profiles. When investigating synchronization phenomena in the epileptic brain (Mormann et al., 2005;Schelter et al., 2007;Osterhage et al., 2008;Lehnertz et al., 2009;Mirowski et al., 2009;Kuhlmann et al., 2010), however, a time-resolved analysis of interactions between all possible pairs of recording sites results in a large number of measure profiles. A significance evaluation will inevitably lead to the statistical problem of multiple testing. For approaches based on analytic random predictors (Winterhalder et al., 2003;Schelter et al., 2006a), this leads to high sensitivities of the random predictors and therefore to reduced statistical power. For Monte-Carlo based approaches, which use randomized surrogate data to test for statistical significance (Andrzejak et al., 2003(Andrzejak et al., , 2009Kreuz et al., 2004), the high number of measure profiles may lead to non-independent surrogate realizations.
Addressing these issues, we have developed a method to preselect measure profiles that carry potential information about preseizure states and to identify brain regions associated with seizure precursor dynamics, by exploiting specific aspects that appear to be associated with seizure precursor dynamics. We demonstrate the suitability of our method through a retrospective time-resolved analysis of the strength of interactions in long-lasting, multichannel iEEG recordings.

MaterIals and Methods clInIcal data
We analyzed continuous iEEG recordings from 18 patients who suffered from pharmacoresistant focal epilepsy with neocortical and/or hippocampal origin (cf. Table 1). All patients underwent resective surgery. Recordings were performed with strip (s), depth (d), and/or grid (g) electrodes (Spencer et al., 2007). Location of the seizure onset zone: hippocampal (H), neocortical (NC). Outcome according to a modified Engel classification (Engel and Rasmussen, 1993). important property of the frequency-adaptive approach is that the instantaneous phase/frequency always relates to the predominant frequency in the spectrum (Boashash, 1992). In iEEG time series this predominant frequency is subject to changes. From an electrophysiological point of view, an extraction of phase information that automatically adapts to the main rhythm locally in time can thus be advantageous particularly when analyzing non-stationary signals such as the EEG. For windows of 30 s duration, the first and last 5 s of each window were tapered by a Hann window to reduce spectral leakage. Windows were shifted forward in time by 10 s. The resulting measure profile was smoothed using a 5-min median filter (Mormann et al., 2003b;Schelter et al., 2006a;Kuhlmann et al., 2010). For the M electrode contacts of each patient, all K = M(M − 1)/2 possible pairs of contacts were considered, resulting in a total of 378-4005 measure profiles (average: 1815.2 profiles).

Preselection method
The preselection method introduced here is based on the idea that any measure profile suitable for studies on the predictability of seizures must contain distinctive signatures in advance of seizure onsets. In contrast to unspecific fluctuations of measure profiles on a short timescale, these changes are expected to occur on a longer timescale (cf. Lehnertz et al., 2007). As quantification thereof, we estimated the ratio of the global to local temporal variance of each measure profile x k , k = 1,…, K, by calculating (1) between subsequent samples; x k denotes the arithmetic mean of x k , and ∆x k the arithmetic mean of ∆x k . The number of data points of a measure profile is denoted by N. In Eq. 1, the factor of two in the numerator imposes a specific normalization; for independently distributed random variables x k , the expectation value of S k would be one. Thus, values of S k not equal to one indicate correlations in the data. In general, measure profiles with Intracranially recorded electroencephalograms were recorded during presurgical epilepsy monitoring at the Department of Epileptology at the University Hospital of Bonn, Germany, and at the Epilepsy Center of the University Medical Center Freiburg, Germany. The retrospective evaluation of the data received prior approval by the local Ethics Committees and informed consent was obtained from each patient.
The iEEG recordings were performed by digital video EEG systems with a sampling rate of 200 Hz (Bonn: Stellate Systems, Montreal, Canada and Schwarzer, Munich, Germany) or with 256 Hz or 1024 Hz (Freiburg: Neurofile NT, IT-Med, Usingen, Germany). The total number of electrode contacts on subdural grid and strip electrodes and/or intrahippocampal depth electrodes ranged from M = 28 to M = 90 (Table 1; see Figure 1 for an exemplary electrode implantation scheme). Data were high-pass filtered at 0.5 Hz, lowpass filtered for anti-aliasing, and digitized with a 16-bit analogto-digital converter. iEEG signals were offline low-pass filtered at 85 Hz; data recorded with 1024 Hz were resampled at 256 Hz in order to achieve a comparable number of data points per analysis window. The overall recording time amounted to 133.1 days during which 133 spontaneous seizures occurred (cf. Table 1). We visually identified the time of seizure onset on the iEEG as the time of earliest clear change from the patient's baseline or normal background iEEG that eventually led to an electrographic seizure. Subclinical seizures were neglected in our analysis. In the majority of patients, antiepileptic medication was reduced during the recording period.

Measuring the strength of interactions
We estimated the mean phase coherence R (Mormann et al., 2000), a bivariate symmetric measure for the strength of interactions between two brain regions, in a time-resolved fashion (moving-window approach). The phases of the iEEG time series can be obtained either with a frequency-adaptive (via the Hilbert transform) or with a frequency-selective approach (e.g., via the wavelet transform). Both approaches were shown to give similar results for a quantitative analysis of phase synchronization from neuronal data (Le Van Quyen et al., 2001;Quian Quiroga et al., 2002;Bruns, 2004) as well as in seizureprediction studies Osterhage et al., 2007). An FiGurE 1 | Electrode implantation scheme for patient no. 4. Intrahippocampal depth electrodes (A) were implanted stereotaxically along the longitudinal axis of the hippocampal formation from an occipital approach with the amygdala as the target for the most anterior electrode. Each catheter-like, 1 mm thick silastic electrode contained 10 cylindrical contacts of a nickel chromium alloy (2.5 mm) every 4 mm. Subdural strip and grid electrodes (B,C) consisted of either four (strips) or 4 × 8 (grid) stainless steel contacts with a diameter of 2.2 mm, embedded in a silastic (intercontact spacing of 10 mm). Strip electrodes (B) were inserted through burr holes and were placed over the left anterior and posterior inferior temporal cortex. The grid (C) was inserted after craniotomy and covered a lesion in the left temporal lobe. cases where the time between two successive seizures was less than T p + 30 min, the maximum amount of data available (i.e., from seizure onset back to the end of the postictal phase of the preceding seizure) was used instead. We derived the ROC curve by varying a threshold for amplitude values of R. The area under the ROC curve (AUC) can be used to quantify the degree to which the two distributions can be distinguished. For identical distributions AUC = 0.5, which indicates that no predictive changes before seizure onset can be detected. We here considered ROC * = |AUC − 0.5|, which allows testing the hypotheses of an increased or a decreased strength of interactions during the preictal period as compared to the interictal period (cf. Mormann et al., 2005). Larger values of ROC * represent a higher predictive performance.

Testing for statistical validity
While the appreciation of the clinical relevance of an analysis technique needs to be based on patient views and requirements of intervention strategies, a decisive test is whether the predictive performance is better than random. To test whether the performance ROC k * of a measure profile k is indeed indicative of a true predictive power, we employed the concept of seizure time surrogates (STS; Andrzejak et al., 2003), a Monte-Carlo based resampling technique for the generation of randomized seizure onset times. In order to preserve the total number of seizures as well as the distribution of time intervals between consecutive seizures, a random permutation of the original seizure intervals was performed. Additionally, the duration of the initial interval was chosen randomly between 0 and 4 h. If the predictive performance ROC k * of measure profile x k with the true sequence of seizures exceeded the maximum ROC STS,k * obtained from 19 realizations of STS, the predictive performance of this measure profile was considered significant at a level of α = 5% (cf. Schreiber and Schmitz, 2000 and references therein).

Testing the relationship between predictive performance and variance ratio
Methods for a data-driven preselection of measure profiles are sensible only if preselected profiles contain more predictive information of an upcoming seizure than the remaining profiles. We therefore hypothesized to observe a close relationship between the variance ratio S and ROC * , which would suggest that measure profiles can be identified which carry relatively high predictive information -without considering any information about seizure onset times in the preselection process. For each patient, we investigated the relationship between S k and ROC k * of all K measure profiles by determining the non-parametric correlation coefficient of Kendall's tau (Kendall, 1976).

Testing the relationship between variance ratio and brain sites
We also investigated a possible relationship between the variance ratio S k and the location of the electrode contacts that the measure profile x k was derived from. For this purpose, we assigned all electrode contacts to four location categories. Thereby, we made use of knowledge concerning location and extent of the seizure onset zone, which is defined by the electrode contacts showing initial ictal activity. We considered the following four categories of electrode contacts: low values of S k exhibit a dominance of short-term fluctuations over long-term variations, whereas profiles with higher values of S k contain more long-term fluctuations, which supposedly carry predictive information (cf. Figure 2).

Estimating predictive performance
In order to test whether measure profiles carry potential information about preseizure states, we estimated their predictive performance in terms of their ability to distinguish between the interictal period and the preictal period. To this end, we followed Mormann et al. (2005), assuming the existence of a preictal period with a duration T p = 4 h. We investigated the separability of amplitude distributions of the mean phase coherence R from the interictal and preictal period using the receiver-operating-characteristic (ROC) statistics which is widely used in studies on the predictability of seizures. We adopted evaluation scheme #2 in Mormann et al. (2005) and estimated, for each measure profile x k separately, the preictal, and interictal distributions from the pooling of all preictal and interictal periods, respectively. Hereby, we expected to find an increased performance particularly for measure profiles for which seizure precursors (either preictally increased or decreased strength of interactions) occur constantly in the same pair of recording sites and on a similar level for all the seizures of a patient.
In order to exclude effects from the postictal period, which can be accompanied by alterations in the EEG, recording periods within 30 min after the electrographical onset of a seizure (Wyllie et al., 2006;So and Blume, 2010 and references therein) were discarded from the analysis. We thus considered seizures with an inter-seizure interval of more than 30 min only ( Table 1 summarizes the number of seizures per patient that entered subsequent analysis steps). In  and S k were significantly positively correlated (p < 0.05; results not shown).
Since statistical evaluation techniques such as the STS used here allow the quantification of the significance of the predictive performance, we can further substantiate our results of the correlation analysis by restricting it to measure profiles with a significant predictive performance. The corresponding results are shown in the histograms of Figure 3. Especially for patients nos. 9 and 17, high values of S were observed for measure profiles with a significant predictive performance. Albeit the multiple testing challenge, the observed correlations indicate that measure profiles with abovechance performance exhibit -on average -a higher variance ratio S than the remaining profiles.

relatIonshIp between varIance ratIo and braIn sItes
For 17 of 18 patients (all but patient no. 15), we observed differences in the distributions of S values among all categories of electrode contact combinations f/f, f/n, f/c, f/e, n/n, n/c, n/e, c/c, c/e, and e/e (p < 0.001; cf. Figure 4). Varying between individual patients, significantly larger or also significantly smaller values of S were found for the different groups of contact combinations tested, i.e., those involving focal, neighboring, contralateral, or extra-focal contacts, than for the respective remaining groups. Details are summarized in Table 2.

coMparIson to preselectIon based on clInIcal InforMatIon
For each of the six patients from the Epilepsy Center Freiburg, an increased average predictive performance could be observed for the 15 measure profiles with highest values of S in comparison to the 15 manually preselected measure profiles (cf. Table 3). Mean predictive performance (ROC * ) of the 15 measure profiles increased for these patients on average from 0.072 to 0.107. For patients nos. 14, 15, and 17, the predictive performances of the 15 measure profiles selected by considering S were significantly higher than those of the 15 manually selected measure profiles (p < 0.01, Mann-Whitney-Wilcoxon signed rank sum test).
For four patients (patient nos. 13, 16, 17, and 18), the majority of the 15 measure profiles preselected with highest values of S belonged to the e/e combination category. For patient no. 14, combinations n/e and also e/e were prominent, and for patient no. 15, the majority of preselected measure profiles belonged to combination category c/e. For only two patients (no. 17 and no. 18), measure profiles from combination category f/f were among the 15 profiles with highest variance ratio S.

dIscussIon
We introduced a method for the data-driven preselection of measure profiles in retrospective studies on the predictability of epileptic seizures. The method is based on the ratio S of the global to local temporal variance of a measure profile. Due to the high correlation between S and the predictive performance of a measure profile (as quantified by the ROC), a preselection based on S enables a reduction of the number of profiles to a small subset containing those that carry potential information about preseizure states. In contrast to previously proposed methods (Iasemidis et al., 2001;D'Alessandro et al., 2005;Le van Quyen et al., 2005), our procedure does not Neighbor (n) Electrode contacts not more than two contacts distant to those from category "f" (on average 17.2%, varying between 6.9 and 42.9%) Contralateral (c) Electrode contacts homologous to those from "f" and "n" in the contralateral hemisphere (available for 12 patients with a bilateral electrode implantation; on average 7.8%, varying between 3.4 and 18.5%) Extra-focal (e) All remaining electrode contacts (on average 67.5%, varying between 32.8 and 89.6%) Since the mean phase coherence is a symmetric measure, measure profiles were thus derived either from f/f, f/n, f/c, f/e, n/n, n/c, n/e, c/c, c/e, or e/e combination categories. The non-parametric Kruskal-Wallis test (Kruskal and Wallis, 1952) was applied to test for differences in the distributions of S values among combination categories. For those patients for whom significant differences were observed, post hoc non-parametric Mann-Whitney-Wilcoxon signed rank sum tests (Mann and Whitney, 1947) were used to analyze whether specific groups of combination categories differed in S. We investigated whether S values of measure profiles derived from focal contacts (f/f, f/n, f/c, f/e) differed from those of the remaining ones. The same was also done for S values of measure profiles derived from neighboring contacts (f/n, n/n, n/c, n/e), from extra-focal contacts (f/e, n/e, c/e, e/e), and for the 12 patients with bilaterally implanted electrodes also for the S values of measure profiles from contralateral contacts (f/c, n/c, c/c, c/e). The significance levels were adjusted in order to correct for the application of the eightfold test.

Comparison to preselection based on clinical information
For the patients who underwent presurgical evaluation at the Epilepsy Center Freiburg, we compared the proposed preselection method to an approach that is based on clinical information (Maiwald et al., 2004;Schelter et al., 2006aSchelter et al., ,b, 2007Winterhalder et al., 2006;Feldwisch-Drentrup et al., 2010). Thereby, three electrode contacts showing initial ictal activity and three contacts not involved at all or latest during spread of ictal activity were preselected by a certified epileptologist (ASB). In order to contrast this approach with the preselection based on S, we compared the predictive performance ROC * of the 15 measure profiles derived from the manually preselected electrode contacts to the ROC * values of the 15 measure profiles preselected with highest variance ratio S.

results relatIonshIp between predIctIve perforMance and varIance ratIo
In Figure 3 we show, for two exemplary patients with unilaterally (no. 9 and no. 11) and for two exemplary patients with bilaterally implanted electrodes (no. 4 and no. 17), the predictive performances ROC k * of all measure profiles x k in dependence on their variance ratios S k . Except for patient no. 11, significant positive correlations (p < 0.05) between predictive performances ROC k * and mean phase coherence as a bivariate measure for the strength of interactions. We note, though, that our preselection method can be applied to any measure profile that is derived from univariate, bivariate, or multivariate EEG analysis techniques (see Mormann et al., 2007;Schelter et al., 2008 for an overview). Nevertheless, further studies are necessary to evaluate the impact of the choice of the duration of preictal and postictal periods and of algorithmic exploit any information about the relation of the measure profiles to the occurrences of seizures for the preselection process. Hence, it may be applied prior to an analysis of predictive performance, based on the same data.
To demonstrate the suitability of the method, we retrospectively analyzed interactions in long-lasting, multichannel intracranial EEG recordings from 18 patients with focal epilepsies using the   parameters. In addition, future studies should also investigate whether precursors of very focal subclinical seizures can be identified on the EEG. Modern EEG acquisition systems enable recordings with high spatial and temporal resolution. The large numbers of recording channels pose severe challenges due to the tremendous increase in the number of measure profiles, particularly for algorithms based on bivariate and multivariate EEG analysis techniques. In such cases, the proposed preselection strategy can be used to quickly screen measure profiles in order to identify those with a presumably high predictive performance. This is especially helpful if only few seizures are available and if a rigorous statistical evaluation of predictive performance is hardly possible due to the statistical problem of multiple testing.
The close relationship between the variance ratio S and predictive performance, which could be observed in 17 out of 18 patients, indicates that our method is beneficial for the preselection of measure profiles carrying potential information about preseizure states. In comparison to an ad hoc preselection criterion based on prior knowledge concerning the location of recording sites relative to the seizure onset zone, measure profiles preselected with S had a considerably higher predictive performance. Interestingly, these profiles mostly captured interactions between non-focal brain regions rather than interactions that involved focal regions. Our investigation concerning a possible relationship between predictive performance and the location of recording sites yielded statistically significant differences in the variance ratios between different categories of combinations of sites. Yet, for the patients investigated here, no consistent patterns for optimal combinations could be observed, which may indicate a high intra-and interindividual variability of the spatial-temporal dynamics underlying the generation of focal onset seizures (cf. Mormann et al., 2003aMormann et al., ,b, 2005Winterhalder et al., 2006;Kuhlmann et al., 2010). This clearly underlines the importance to study brain regions outside the epileptic focus, which yet might be part of an epileptic network Kuhnert et al., 2010) and contribute to ictogenesis. The concept of an epileptic network, whose interactions extend over large regions of the brain, might also be of interest for computational neuroscience approaches to epilepsy (Lytton, 2008;Soltesz and Staley, 2008).
To summarize, investigating the ratio of the global to local temporal variance of a measure profile constitutes an easy and efficient way to identify profiles that carry information predictive of impending seizures and thus to identify brain regions possibly involved in seizure generation. Hence, our method can help to improve knowledge about the complex spatial and temporal aspects of the seizure generating process and to substantiate the existence of preseizure states in epileptic brain networks. This may eventually facilitate the development of seizure-prediction techniques to be used in clinical applications.