Machine Learning for Guiding Neurosurgery in Epilepsy

Justin Dauwels School of Electrical and Electronic Engineering Nanyang Technological University (NTU), Singapore [email protected] Emad Eskandar, Andy Cole, Dan Hoch, Rodrigo Zepeda, Sydney S. Cash Massachusetts General Hospital and Harvard Medical School Cambridge, MA, USA

Abstract For approximately 30% of epilepsy patients, seizures are poorly controlled with medications alone. For those patients, surgery may be an option: the aim is to remove the brain area(s) where the seizures originate. The key to success is to be able to reliably and accurately localize the seizure onset zone. In order to delineate that zone, one heavily relies on brain signals recorded during seizures (ictal recordings), by electrodes that are semi-chronically implanted in the patient’s brain. Since seizures occur only occasionally, patients often have to stay 5 days or more in the hospital until sufficient seizures have been recorded. Obviously, this procedure is costly, uncomfortable, and not without risk for infection and other side effects. In this paper, decision algorithms are developed that use periods of intracranial non-seizure (interictal) EEG to localize epileptogenic networks. The proposed algorithms combine spectral and multivariate statistics in a decision-theoretic framework to automatically delineate the seizure onset area. In the case of depth recordings, we apply standard binary classification algorithms, including linear and quadratic discriminative analysis. For the surface recordings, we derive novel decision algorithms, based upon graphical models; the latter encode common patterns among patients and also patient specific abnormalities. By conducting inference in the graphical model, we derive collaborative filtering algorithms. The outcomes from our algorithms for both depth and surface recordings are in good agreement with the determination of the seizure focus by clinicians from ictal EEG, which is the current gold standard. In the long term, our approach may lead to shorter hospitalization of intractable-epilepsy patients, since it does not rely on ictal EEG.

1

Introduction

Approximately 50 million people worldwide (2.5 million in the United States alone) have epilepsy. More than 50 percent of those suffer from localization-related epilepsy. Unfortunately, 30% of these patients continue to have seizures despite maximal medical therapy (see, e.g., [1]). Furthermore, many patients suffer from considerable side effects of the medications. On the other hand, regional surgical resection may provide seizure reduction or even cure [2]. However, it is of crucial importance to reliably localize the epileptic brain area(s). At present, one relies mostly on (scalp or intracranial) EEG that contains seizure activity (“ictal EEG” or “seizure” EEG) to determine the seizure onset area; since seizures usually do not occur frequently, recordings must last a long time (from several days to several weeks) until sufficient seizures have occurred (typically between 3 and 5). 1

In this paper, we investigate whether non-seizure (a.k.a. interictal) intracranial EEG can be used to localize epileptic brain tissue. In particular, we develop novel statistical decision algorithms that detect, combine, and leverage various abnormal interictal electrophysiological characteristics to delineate epileptic brain tissue (see Fig. 2(a)). In order to delineate the seizure onset zone from non-seizure EEG, we will exploit two phenomena: slowing and locally enhanced EEG synchrony (“hypersynchrony”). Signals recorded from damaged cortex often seem to be “slower ”, i.e., contain more power at low frequencies. Second, a unifying principle emerging from decades of intense research is that seizures are a property of abnormally firing neurons that, entrained by an imbalance of excitation and inhibition, discharge synchronously in a critical ensemble [3]. Moreover, several studies have suggested that multivariate analysis of seizure-free (rest) EEG, in which the relationship between different channels of activity are compared, may help to delineate epileptogenic cortex (e.g., [4, 5]; see also [6] for a recent review). No study so far, however, has combined several electrophysiological signatures to localize the seizure onset zone from interictal EEG. Virtually all studies focus on one particular signature, and report negative results, showing that one signature does not suffice to delineate seizure onset reliably from interictal EEG (e.g., [4, 5]). This paper demonstrates that a combined approach, founded on algorithmic decision making procedures, substantially improves performance. In future work, we will integrate more signatures into our method (e.g., high-frequency bursts and interictal spikes) to further improve its reliability. In the long term, one may therefore no longer need to rely on seizure EEG, but instead use short non-seizure EEG recordings to determine the seizure onset area; this would drastically reduce the hospitalization time for intractable-epilepsy patients. Our methodology would be useful for focal epilepsies regardless of underlying etiology. This includes focal epilepsy secondary to cortical dysplasia, tuberous sclerosis, a stroke, tumor, vascular malformation, or trauma. This paper is structured as follows. In Section 2 we describe our EEG data and the preprocessing we carried out. In Section 3 we investigate the phenomena of slowing and hypersynchrony, and in Section 4 we show that those two electrophysiological signatures correlate with the seizure onset zone. In Section 5 we briefly describe our algorithms that merge univariate and multivariate measures to infer the seizure focus, and discuss their performance. At the end of the paper we offer some conclusions.

2

EEG Data and Preprocessing

We investigated data from 11 patients with intractable localized epilepsy who had undergone intracranial investigation using intracortical electrodes (depth; see Fig. 1(a)) and/or subdural surface electrodes (grid; see Fig. 1(d)). In particular, we analyzed data from 5 patients with only depth electrodes and 6 patients with both grid and depth electrodes. (In the near future, we will investigate data from patients with non-localized epilepsy as well, which will serve as control population.) The depth electrodes contain multiple channels (6 to 8); the grids typically contain 8 × 8 electrodes. In each case a segment of data 1 hour long and at least 24 hours separated from seizure activity was examined. (In later studies, we will investigate whether smaller separation from seizure activity affects the predictability). The data was band-pass filtered between 1 and 200Hz, and a notch filter was applied to remove the 60Hz power signal components. Before computing the univariate and multivariate measures, each EEG signal was normalized (mean subtracted, divided by standard deviation). No additional preprocessing was carried out, and no EEG channels have been selected based on any pre-existing knowledge, except that clearly dysfunctional data channels were discarded. So far, we have applied one univariate measure (relative power), and several multivariate measures, including Pearson correlation coefficient [7], magnitude coherence [7], phase synchrony [8], and omega complexity [9]. Note that Pearson correlation coefficient [7], magnitude coherence [7], phase synchrony [8] are bivariate measures, whereas omega complexity [9] can be applied to more than two signals simultaneously. For depth electrodes, the univariate measures were applied to all channels of a given depth electrode, and then averaged over all channels, resulting in one average value per depth 2

electrode array (see Fig. 1(b)). Pairwise synchrony measures were applied to all pairs of channels and then averaged, while omega complexity was calculated from all channels simultaneously. This resulted in a single synchrony value per depth electrode array (see Fig. 1(c)). The measures varied considerably across patients necessitating a normalization procedure. We divide each measure by the computed average of all depth electrodes of a given patient. For grid electrodes, the univariate measures were applied to the signal recorded by each electrode; this leads to univariate maps across the grid (see Fig. 1(e)). For the pairwise synchrony measures, we calculated local synchrony by averaging the pairwise value between the electrode of interest and all of its neighbors. We applied omega complexity simultaneously to the electrode of interest and all its nearest neighbors. Because electrodes at the edge of the grid have fewer neighbors, we normalized the synchrony measures by the number of adjacent electrodes; we divided the measures by the average computed over all electrodes with the same neighborhood size: the 4 corners (neighborhood size 3), the other channels at the grid boundaries (size 5), and the inner channels (size 8). From these calculations we obtained local synchrony maps (one for each synchrony measure; normalized and unnormalized), which display the local synchrony of the grid electrodes (see Fig. 1(f)).

3

Stable patterns of slowing and hypersynchrony occur in most patients

We have observed that all 11 patients show focal slowing, i.e., increased relative power in the 1–5Hz frequency range at certain subsets of channels, and therefore cortical regions. In some patients, this effect was weak (relative increase of 5%) in others, very significant (relative increase of 80–90%). Examples of such patterns are shown in Fig. 1(b) (depth) and Fig. 1(e) (grid); both figures show the relative power in the 1–5Hz frequency range. Electrode 1 in Fig. 1(b) (right hemisphere), indicated by a black square, shows 12% increase in relative power; the electrode at position (5,8) in Fig. 1(e), also indicated by a black square, exhibits 60% increase in relative power. In 9 of the 11 patients, stable synchrony patterns occur with certain subsets of channels, and therefore cortical regions, showing significantly enhanced local synchrony (“hypersynchrony”) compared to other areas. Different synchrony measures typically lead to slightly different patterns. Examples of such patterns are shown in Fig. 1(c) (depth) and Fig. 1(f) (grid), in which local synchrony was calculated using the magnitude coherence method. Interestingly, all patients show either focal slowing or enhanced local synchrony; most patients show both phenomena simultaneously.

4

Slowing and hypersynchrony correlates with seizure onset zones

In most of the 11 cases examined to date, focal slowing and/or locally enhanced synchrony (“hypersynchrony”) occurred close to or within the seizure onset zone (determined by the clinical team from ictal EEG), as illustrated in Fig. 1. In the depth case of Fig. 1(a) to Fig. 1(c), slowing occurs at electrode R1. The electrodes R1 and R2 are more synchronous than the other ones (“hypersynchronous”). The clinical team had determined that the seizures originate at the right anterior temporal lobe, in particular, electrodes R1 and R2; that area indeed shows slowing and hypersynchrony. The seizure onset zone in the grid case of Fig. 1(d) to Fig. 1(f) is located at position (5,6), marked in blue in Fig. 1(d) (“actual”). Slowing and hypersynchrony occurs in the close neighborhood: slowing occurs at positions (5,6) and (5,8) (see Fig. 1(e)), and hypersynchrony can be observed at electrodes (5,7) and (5,8) (see Fig. 1(f)).

5

Decision algorithms to delineate seizure onset zone

In this section, we elaborate on our decision algorithms for localizing the seizure onset from interictal depth and surface recordings. Since those two types of recordings are quite distinct, we have developed separate algorithms for each type of recording. An important difference between depth and grid recordings is the spacing between electrodes: The spacing between depth electrodes is fairly large (several centimeters), and each depth electrode can 3

(a)

(d)

(b)

(e)

(c)

(f)

Figure 1: Examples of univariate and multivariate decision making. a) 5 depth electrodes in right hemisphere; actual and predicted onset area (by combining relative power and magnitude coherence) coincide. b) Relative power in 1–5Hz band for 10 depth electrodes (5 in left and 5 in right hemisphere); black box indicates decision based only on relative power. c) Magnitude coherence for the 10 depth electrodes, where increased synchrony can be observed in electrodes R1, R2; black box indicates decision only based on magnitude coherence. d) 8x8 subdural grid; the actual and predicted onset areas (by combining relative power and magnitude coherence) are in close vicinity. e) Relative power in 1–5Hz band for each grid electrode; black box indicates decision based on relative power. f) Magnitude coherence; black box indicates decision based on magnitude coherence. be treated separately in the decision making; in contrast, the distance between neighboring grid electrodes is much smaller (5–7mm), and consequently, the measurements from nearby grid electrodes are often correlated; this spatial correlation needs to be taken into account in the decision making. We first describe our algorithms for the depth recordings, next we elaborate on surface recordings. 5.1

Depth Recordings

We aim to determine which depth electrodes are located inside the seizure focus; that problem may be viewed as binary classification [10]. We utilize a supervised learning paradigm in which the clinical determinations of seizure focus (gold standard) are used to train binary classifiers with the (normalized) univariate and multivariate measures as input features. Specifically, we conduct linear and quadratic discriminant analysis (LDA and QDA) where we tested all possible combinations of the (normalized) univariate and multivariate measures as features [10]. In all our experiments, we compute the classification rate, specificity, and sensitivity through leave-one-out crossvalidation [10]. 5.2

Surface Recordings

Similarly as for depth recordings, our objective is to determine which grid electrodes are localized inside the seizure onset zone. In principle, that problem may also be viewed as binary classification. However, our experiments have confirmed that such approach performs poorly (not shown here). Since the univariate and multivariate measures are prone to noise, they vary sometimes significantly across the grid; those fluctuations are random in nature, and do not necessarily reflect abnormal electrophysiological markers. To overcome this issue, we have developed a novel statistical decision algorithm that exploits the correlation between grid recordings to smoothen the random fluctuations across the grid. The algorithm makes use of the statistical distribution of the univariate and multivariate measures across electrode locations and patients; it computes how likely any given electrode is located inside 4

(a)

(b) OR

decision map x AND

AND

p0 (yi |xi ) xi exp(J xi xj )

normalized

unnormalized

measure 1

(c)

measure 2

(d)

Figure 2: Decision making algorithm and graphical representation of measures of accuracy. a) Diagram of the decision making procedure, which combines various measures to compute a probability map of seizure onset; the predicted seizure onset zone is determined as the area in which the seizure onset probability is above a certain threshold. b) Example in which the decision (in red) can be compared to what was determined clinically (in blue). c) Ising model encodes spatial continuity of the measures. d) Logical AND gates merge the normalized and unnormalized values of each measure; logical OR gates merge the different univariate and multivariate measures. the seizure focus, given the univariate and multivariate statistics, thereby taking correlations between neighboring electrodes into account. An electrode is considered to be part of the seizure onset zone if the posterior probability of belonging to the seizure onset (computed in the previous step) is above a certain threshold; the latter is computed adaptively, and depends on the posterior probabilities. In the following, we will outline the underlying mathematical formalism. For each surface electrode, we need to decide whether or not it is located inside the seizure focus. We associate a binary variable xi with each surface electrode Ei (with i = 1, 2, . . . , N, and N is the number of surface electrodes), where xi = 1 if the electrode Ei is located inside the seizure focus, and xi = -1 otherwise. Formally, our objective is to infer the binary sequence x = x1 , x2 , . . . , xN from the univariate and multivariate statistics. For now, let us consider one measure y, e.g., relative power. Similarly as in the depth cases, we normalize that measure, resulting in normalized values y ′ . Both the normalized values y ′ and unnormalized values y are taken into account in our decision algorithm. We will first focus on the former. We compute a histogram of the normalized values, across all electrodes and all EEG segments of a given length, denoted by fH (y ′ ); that histogram is an approximation of the (unknown) probability density function f (y ′ ) of the (normalized) measure y ′ . We compute such histogram for EEG segments of different lengths (e.g., 2, 5, 15, 20, 30, 45 min). When a measure is averaged over a short EEG segment (e.g., 2min), it may fluctuate more than when averaged over a long segment (e.g., 45min); the histograms corresponding to different segment lengths may therefore be considerably different. From the histogram fH (y ′ ), we compute the cumulative histogram FH (y ′ ); the latter is a mapping 5

that counts the cumulative number of observations in all of the bins up to the specified bin; it is an approximation of the probability distribution F (y ′ ) of y ′ . The cumulative histogram FH (y ′ ) provides useful statistics to detect abnormal values; in this context, we assume that From FH (y ′ ) we define the likelihood p(yi′ |xi ) that the electrode at hand is located inside the seizure focus: α ′ α ′ p(yi′ |xi = 1) ∝ FH (yi ) and p(yi′ |xi = −1) ∝ (1 − FH (yi )),

with exponent α 1) and p(yi′ |xi =

′

(1)

> 1. Since our objective is to infer x for given y , only the ratio of p(yi′ |xi α ′ α ′ −1) matters, which is equal to the ratio of FH (yi ) and (1 − FH (yi )).

=

In addition to the information from the univariate and multivariate statistics, we also have prior knowledge about the seizure focus. Usually the seizure focus is located in one area or a few areas, i.e., the seizure focus typically does not consist of numerous disconnected areas. Moreover, several neighboring surface electrodes are typically located in the seizure focus. As a consequence, if an electrode is part of the seizure focus, then probably its neighbors are as well. We encode this continuity constraint with an Ising model: Y p(x) ∝ exp(Jxi xj ) (2) i,j∈N (i)

where N (i) are the nearest neighbors of electrode i and J > 0 is a real number. Suppose that xi and its neighbor xj are both equal to 1, and hence the electrodes Ei and Ej are both inside the seizure focus; the exponential factor exp(Jxi xj ) then evaluates to exp J. On the other hand, if xi and its neighbor xj are equal to 1 and -1 respectively, and hence Ei and Ej are inside and outside the seizure focus respectively, the exponential factor exp(Jxi xj ) evaluates to exp −J, and the prior p(x) will be smaller in that case. In summary, the prior model p(x) imposes spatial continuity, and favors connected seizure foci. By combining the prior p(x) with the likelihood functions p(yi′ |xi ) we obtain a statistical model p(y ′ , x) that relates the decisions x to the normalized measure y ′ : Y Y Y p(y ′ , x) = p(x) p(yi′ |xi ) ∝ exp(Jxi xj ) p(yi′ |xi ). (3) i

i

i,j∈N (i)

The model p(y ′ , x) is depicted in Fig. 2(c); more specifically, the latter shows a factor graph that represents the factorization (3) of p(y ′ , x) [11]. The green circles correspond to the variables xi , the blue and red squares represent the factors exp(Jxi xj ) and p(yi′ |xi ) respectively. A circle is connected to a square if and only if the variable associated with the circle occurs in the factor associated with the square [11]. In the factorization (3), xi occurs in (two, three, or four) factors exp(Jxi xj ) and one factor p(yi′ |xi ), and consequently each green circle in Fig. 2(c) is connected to (two, three, or four) blue squares and one red square. If both the normalized and unnormalized values are abnormal at an electrode Ei , the latter may be inside the seizure focus, otherwise it is probably not. Similarly to p(y ′ , x) (3), we introduce a model p(y, x) for the unnormalized values y. We merge the models p(y, x) and p(y ′ , x) through logical AND gates: Y δ[AND(x′i , x′′i ) − xi ], p(y, y ′ , x, x′ , x′′ ) = p(y, x′′ )p(y ′ , x′ ) (4) i ′

′′

where x and x are auxiliary binary sequences, AND(z, z ′ ) stands for the logical AND of the binary variables z and z ′ , and δ[.] is the Kronecker delta. The auxiliary binary variables x′i and x′′i (with i = 1, 2, . . . , N ) are associated with the electrodes Ei ; they may be viewed as intermediate binary decision variables that encore whether the normalized and unnormalized values y ′ and y respectively are abnormal, whereas xi encodes whether both y and y ′ are abnormal simultaneously. If both the normalized and unnormalized values are abnormal at electrode Ei (and hence xi = 1), the latter may be located inside the seizure focus, otherwise it is probably not. So far, we have assumed that only one measure y is available. In practice, we compute several univariate and multivariate measures y = y (1) , y (2) , . . . , y (M) (with M the number 6

of measures), and we wish to combine them to delineate the seizure focus. We will denote the normalized measures by y′ = y ′(1) , y ′(2) , . . . , y ′(M) , and the three sets of decision variables associated with the measures y = y (1) , y (2) , . . . , y (M) by x′ = x′(1) , x′(2) , . . . , x′(M) , x′′ = x′′(1) , x′′(2) , . . . , x′′(M) , and x ˜ = x(1) , x(2) , . . . , x(M) ; the latter are combined together through logical OR gates, resulting in the binary decision sequence x that encodes whether any of the measures is abnormal. We eventually use that sequence x to delineate the seizure focus. The overall statistical model is given by: Y  Y (1) (M) p(y, y′ , x, x′ , x′′ , x ˜) = p(x) δ OR(xi , . . . , xi ) − xi p(y (m) , y ′(m) , x(m) , x′(m) , x′′(m) ), m

i

(5) where OR(z1 , z2 , . . . , zn ) stands for the logical OR of the binary variables z1 , z2 , . . . , zn , p(x) is given by (2), and the factors p(y (m) , y ′(m) , x(m) , x′(m) , x′′(m) ) are models of the form (4). A factor graph of model (5) is shown in Fig. 2(d) for M = 2 measures; it contains an Ising model for each normalized and unnormalized measure y ′(m) and y (m) respectively (cf. Fig. 2(c)). The Ising models for the normalized and unnormalized values of a measure are coupled through logical OR gates; the different measures are merged through logical AND gates. In Fig. 2(d) logical OR and AND gates are only shown for one set of variables ′(1) ′′(1) (1) ′(1) ′′(2) (2) xi , xi , xi , xi , xi , xi , and xi , so as not to overload the figure. We wish to infer the decision variables x from the unnormalized and normalized measures y and y′ respectively. The electrodes with the largest evidence p(xi = 1|y, y′ ) are considered to be part of the seizure focus, where p(xi = 1|y, y′ ) is computed from (5) by marginalization. More precisely, an electrode Ei is considered to be located inside the seizure focus (and hence x ˆi = 1) if and only if p(xi = 1|y, y′ ) > ρ pmax , where pmax = maxi p(xi = 1|y, y′ ) is the maximum evidence, and 0 < ρ < 1. In words: only those electrodes Ei whose evidence p(xi = 1|y, y′ ) is sufficiently strong (larger than the fraction ρ from the maximum evidence) are considered to be part of the seizure focus. In principle, the marginals p(xi = 1|y, y′ ) may be computed by dynamic programming (or equivalently, by sum-product message passing on a cycle-free factor graph of (5)) [11]. However, the complexity of such algorithm is prohibitive for the model at hand, since it involves sums over an exponential number of configurations. As an alternative, we compute the marginals p(xi = 1|y, y′ ) approximately by applying sum-product message passing to the cyclic graph depicted in Fig. 2(d) [11]. The algorithm may be derived straightforwardly, since all involved variables are binary. The above decision procedure depends on several parameters: the Ising coefficients J, the exponents α, and the fraction ρ. We have optimized those parameters manually; in future work, we will design procedures to infer those parameters automatically from data. 5.3 Numerical results Here we present results from our decision algorithms. Due to space constraints, we limit ourselves to results for the most discriminative pair of features, i.e., relative power and correlation coefficient. Our experiments have shown that adding more features only slightly improves those results. In the depth cases (see Fig. 3(a)), the specificity and sensitivity using relative power alone is 78% and 88% respectively; for correlation coefficient it is 29% and 80% respectively. When combined the specificity and sensitivity increase to 80% and 88% respectively, which is a minor improvement over the results based on relative power alone. Those results were obtained with QDA, which consistently performed better than LDA for the data at hand. In the grid cases, the predicted and actual seizure onset zone may be very close to each other without any overlap (see, e.g., Fig. 1(d)), in which case both specificity and sensitivity are zero. Therefore, we also compute distance scores: for each electrode in the predicted onset area, we compute the distance to the closest electrode in the actual onset area, and vice versa, as illustrated in Fig. 2(b). For example, electrode A in Fig. 2(b) is part of the predicted onset area, and its distance to the actual onset area is kAAk; likewise, electrode B is part of the actual seizure onset area, and its distance to the predicted onset area is kBBk. The distances kAAk, kBBk, etc. are Euclidean distances computed in grid coordinates; for example, a distance of 2 corresponds to twice the distance between neighboring electrodes; 7

Specificity

0.8

Power Corr Combined

0.6 0.4 0.2 0.2

0.4

0.6 0.8 Sensitivity

(a)

1

Median distance (decision) 4

Power MagCoh Combined

3 2 1 0

10 20 30 40 Segment Length [min]

(b)

Median distance (onset) Median distance (onset)

Sensitivity and Specificity (depth) 1

Median distance (decision)

the latter is usually about 5mm. We compute those distances for each electrode in the actual and predicted onset area, and summarize the statistics by two numbers: the median distance for the actual and predicted onset area. As can be seen from Fig. 3(b) and 3(c), the median distances for actual and predicted onset areas are about 2 (±10mm) for relative power and the correlation coefficient alone, and are 0.6 (±3mm) and 0.4 (±2mm) respectively when combined. 4

Power MagCoh Combined

3 2 1 0

10 20 30 40 Segment Length [min]

(c)

Figure 3: Aggregate measures of performance of the decision heuristic based on power, correlation coefficient, or a combination of the two. a) Sensitivity and specificity for depth cases. b) Median distance between predicted onset and closest actual onset for grid cases. c) Median distance between actual onset and closest predicted onset for grid cases.

6

Conclusions

In the present work, we have presented algorithms that automatically infer the seizure onset zone from interictal recordings, by exploiting various electrophysiological signatures. With those algorithms in hand, we will in the near future organize a multi-center trial in which short-length interictal data is used to determine the seizure onset zone and a recommended resection is offered prospectively. Outcomes in terms of both seizure freedom and complications will be analyzed. We will also try to incorporate short-term intra-operative recordings into our algorithm. Demonstration that an accurate seizure onset zone prediction could be obtained from OR recordings only would eventually obviate the need for semi-chronic invasive recordings.

References [1] T. Keranen, M. Sillanpaa, et al., “Distribution of seizure types in an epileptic population,” Epilepsia 29(1):1–7. [2] J. Engel Jr., S. Wiebe, et al., “Practice parameter: Temporal lobe and localized neocortical resections for epilepsy: Report of the Quality Standards Subcommittee of the American Academy of Neurology, in Association with the American Epilepsy Society and the American Association of Neurological Surgeons,” Neurology 60(4): 538–47, 2003. [3] A. Wyler and A. Ward, Epilepsy: A Window to Brain Mechanism, New York, Raven, 1992. [4] C. Schevon, J. Cappell, et al., “Cortical abnormalities in epilepsy revealed by local EEG synchrony,” Neuroimage 35(1):140–8, 2007. [5] G. Ortega, L. M. de la Prida, et al., “Synchronization clusters of interictal activity in the lateral temporal cortex of epileptic patients: Intraoperative electrocorticographic analysis,” Epilepsia 49(2):269–280, 2008. [6] Lehnertz K., Bialonski S., et al.,“Synchronization phenomena in human epileptic brain networks,” Journal of Neuroscience Methods 183(1):42–48, 2009. [7] P. Nunez and R. Srinivasan, Electric Fields of the Brain: The Neurophysics of EEG, Oxford University Press, 2006. [8] J.-P. Lachaux, E. Rodriguez, et al., “Measuring Phase Synchrony in Brain Signals,” Human Brain Mapping 8:194–208 (1999). [9] Saito, N., T. Kuginuki, et al., “Global, regional, and local measures of complexity of multichannel EEG in acute, neurolepticnaive, first-break schizophrenics,” Society of Biological Psychiatry 43: 794-802, 1998. [10] R.O. Duda, P.E. Hart, et al., Pattern Classification (2nd Edition), Wiley-Interscience, October 2000. [11] Loeliger H.-A., Dauwels J., et al., 2007. The factor graph approach to model-based signal processing. Proceedings of the IEEE 95(6), 1295–1322.

8