<%BANNER%>

Spatio-Temporal Dependency Analysis of Epileptic Intracranial Electroencephalograph


PAGE 1

SPATIO-TEMPORAL DEPENDENCY ANALYS IS OF EPILEPTIC INTRACRANIAL ELECTROENCEP HALOGRAPH By ANANT HEGDE A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2006

PAGE 2

Copyright 2006 by ANANT HEGDE

PAGE 3

This dissertation is dedicated to my famil y, teachers and friends for their enduring love, support and friendship.

PAGE 4

iv ACKNOWLEDGMENTS It has been a privilege being Dr. Jose Pr incipe’s student. His constant guidance, encouragement and technical directions helped me tremendously in approaching research challenges. The fact that he has infinite tim e for his students makes him very special in my eyes. No wonder his enthusiasm, his quest for excellence and hi s abilities to inspire and motivate have taken him a long way. My experience as a researcher under his supervision has been very rich and fulfilli ng and I would like to thank him for teaching me how to research, for giving me enough fr eedom to explore problems independently and also for his valuable advice whenever I needed it most. I would like to thank Dr. Ch ris Sackellares for being on my committee and offering me plentiful guidance in resear ch. It is always a privilege to know someone like Dr. John Harris. I thank him for his support, both mora l and technical, and also for serving in my committee. I also thank Dr. Mingzhou Ding fo r being on my committee and his valuable feedback in my research. It is seldom that someone working in CNEL would not know or have heard of Dr. Deniz Erdogmus. I consider myself lucky having been around when Deniz was here. Besides being a great friend, he offered rele ntless guidance all along my research and his contributions in this accomplishment are as pr ecious as anyone else’s. I would also like thank Dr. Deng Shiau and Dr. Yadunandana Rao for their technical feedback. I also wish to extend my acknowledgements to all the members of CNEL who have been primarily responsible for my fruitful st ay in the lab. In particular, I would like to

PAGE 5

v cherish those productive disc ussions with Mustafa Can, Rui Yan Jianwu Xu, Puskal, Sudhir, Han, Kyu-hwa and Jack. Overall, my stay in CNEL has been very rewarding and if I look back a few years from now, I am sure I’ll say ‘No doubt, those were my formative years.’ I would like to extend my gratitude to the always cheerful Ellie Goodwin for her golden words of wisdom. He r ability to get things done was truly remarkable. I would also like to acknowl edge Linda Kahila and Shannon for their extensive support and as sistance during my stay at UF. Ka rtikeya Tripati has been very inspiring all along and I am very gr ateful to his constant support. I would like to thank my family and friends for their constant love and encouragement. They have allowed me to pursu e whatever I wanted in life. Without their guidance and affection, it would have been impossible for me to be where I am today. Lastly, I would like to tha nk my life partner Madhoo for making my life so colorful and for being on my side whenever I need her. Her everlasting love has made me a better individual.

PAGE 6

vi TABLE OF CONTENTS page ACKNOWLEDGMENTS.................................................................................................iv LIST OF TABLES...........................................................................................................viii LIST OF FIGURES...........................................................................................................ix ABSTRACT......................................................................................................................x ii CHAPTER 1 SPATIO-TEMPORAL INTERACTIONS...................................................................1 1.1 Introduction.............................................................................................................1 1.2 Different Kinds of Synchronization........................................................................3 1.3 Second Order Synchronization Measures...............................................................4 1.3.1 Cross-Correlation.........................................................................................5 1.3.2 Partial Directed Coherence (PDC)...............................................................6 1.3.3 Synthetic simulation.....................................................................................8 1.4 Mutual Information...............................................................................................11 1.5 Phase Synchronization..........................................................................................12 1.5.1 Hilbert Transformation...............................................................................13 1.5.2 Empirical Mode Decomposition (EMD)....................................................15 1.5.3 Wavelet decomposition..............................................................................18 1.5.4 Indexing Phase Synchronization................................................................19 1.6 Non-Linear synchronization measures.................................................................22 1.6.1 Signal Reconstruction.................................................................................22 1.6.2 Non-Linear Synchronization Measures......................................................23 1.7 Objectives and Author’s Contribution..................................................................25 2 SELF-ORGANIZING-MAP BASED SI MILARITY INDEX MEASURE...............27 2.1 Introduction...........................................................................................................27 2.2 The Similarity Index Technique...........................................................................27 2.3 SOM-based Similarity Index (SOM-SI)...............................................................28 2.4 Simulation Results................................................................................................31 2.4.1 Rosseler-Lorenz simulation........................................................................31 2.4.2 Intracranial EEG Simulation......................................................................35 2.5 Epileptic Intracranial EEG Data Description.......................................................38

PAGE 7

vii 2.6 Testing the application of SOMSI on Epileptic Intracranial EEG......................38 2.7 Statistical comparison between SI and SOM-SI...................................................40 2.8 Testing the Robustness of SOM-SI on Multiple SOMs.......................................42 2.9 Summary...............................................................................................................46 3 SPATIO-TEMPORAL CLUSTERING MODEL......................................................49 3.1 Introduction...........................................................................................................49 3.2 Model for Spatio-temporal clustering...................................................................50 3.3 Temporal Quantification Using Hamming Distance............................................53 3.4 Simulations with Synthetic Data..........................................................................55 3.4.1 Roessler – Lorenz System..........................................................................55 3.4.2 Linear Model..............................................................................................58 3.5 Summary...............................................................................................................64 4 APPLICATION OF SPATIO-TEM PORAL MODEL TO EPILEPTIC INTRACRANIAL EEG..............................................................................................65 4.1 Introduction...........................................................................................................65 4.2 Application of SOM-SI on Ep ileptic Intracranial EEG Data...............................65 4.3 Results...................................................................................................................6 9 4.4 Statistical Validation.............................................................................................75 4.4 Summary...............................................................................................................89 5 STATISTICAL ASSESSMENT OF SPATIO-TEMPORAL DEPENDENCY CHANGES IN EPILEPTI C INTRACRANIAL EEG................................................91 5.1 Introduction...........................................................................................................91 5.2 Spatio-temporal Changes......................................................................................91 5.3 Mantel Test of Matrix Comparison......................................................................92 5.4 Application to Epileptic Intracranial EEG data....................................................94 5.5 Statistical Approach..............................................................................................96 5.6 Inter-Hemisphere Synchronization Differences...................................................98 5.7 Statistical Assessment.........................................................................................104 5.8 Experimental Results..........................................................................................105 5.8.1 Experiment #1..........................................................................................105 5.8.2 Experiment #2..........................................................................................107 5-9 Summary............................................................................................................108 6 CONCLUSIONS AND FUTURE DIRECTIONS...................................................112 6-1 Conclusions........................................................................................................112 6-2 Future Research Directions................................................................................114 LIST OF REFERENCES.................................................................................................119 BIOGRAPHICAL SKETCH...........................................................................................129

PAGE 8

viii LIST OF TABLES Table page 2-1 Coupling strength between pairs of signa ls, the normalized similarity index and the original similarity index between them..............................................................37 2-2 Quantitative comparisons between the SOM-SI profiles obtained from SOM-1 and SOM-2. LOF3 & LOF4 data was projected on each of the SOMs and then the SOM-SI measure was applied to anal yze the dependency of LOF3 on LOF4..46 2-3 Summary of the comparisons betwee n the SOM-SI profiles from SOM-1 and SOM-2. Each row represents the statis tics (mean and varian ce) of pair-wise SOM-SI analyses of the epileptic ECOG data from 6 channels (15 combinations)...........................................................................................................48 3-1 Cluster groupings for the simulation example.........................................................58 4-1 Spatio-temporal groupings as obtai ned for seizure 11 of patient P093....................74 4-2 P093, Seizure 11: Over each 30 minute (91 samples total) window, number of times the within-cluster interaction is gr eater than between-clu ster interaction, at 95% significance level.............................................................................................77 4-3 Spatio-temporal groupings as obtained for seizures 4 and 5 of patient P093..........80 4-4 Spatio-temporal groupings as obtained for seizure 6 and 7 of patient P093............80 4-5 Spatio-temporal groupings as obtai ned for seizure 1 of Patient P092.....................86 4-6 Spatio-temporal groupings as obtai ned for seizure 3 of Patient P092.....................86 4-7 Spatio-temporal groupings as obtai ned for seizure 4 of Patient P092.....................87 4-8 Spatio-temporal groupings as obtai ned for seizure 5 of Patient P092.....................87 4-9 Spatio-temporal groupings as obtai ned for seizure 6 of Patient P092.....................88 4-10 Spatio-temporal groupings as obtai ned for seizure 7 of Patient P092.....................88

PAGE 9

ix LIST OF FIGURES Figure page 1-1 Coupling diagram showing the directi on of interactions between nodes, x1, x2 and x3......................................................................................................................... 8 1-2 Matrix layout plots for PDC, describi ng the interactions in example simulation......9 1-3 Synthetic Sinusoidal signal X(t) a nd its components x1, x2 and x3. X(t) also consists of noise w~N(0, 0.1)...................................................................................17 1-4 Decomposition of sinusoids using EMD..................................................................17 1-5 Hilbert spectrum constructed from the IMF functions.............................................18 1-6 Decomposition of sinusoids using wavelets.............................................................19 1-7 Hilbert spectrum constructed from the wavelet reconstructed narrow-band signals.......................................................................................................................2 0 2-1 Phase-space trajectories of the Rosse ler-Lorenz system for various coupling strengths...................................................................................................................33 2-2 Illustrating the dependency relationshi ps between the Rosseler (driver) and the Lorenz system (response) as a function of coupling strength..................................34 2-3 Original EEG signals X and Y. for the EEG simulation example...........................34 2-4 The nonlinearly mixed (syntheti c) EEG signals V, Z, and W..................................34 2-5 Schematic diagram of the coupling function between the signals...........................35 2-6 The phase-space trajectory of the training EEG signal (projected in two dimensions) and the weights of the trained EEG-SOM (circles).............................36 2-7 Diagram of the depth and subdural elec trode montage in an epileptic brain...........39 2-8 The phase-space trajectory of the traini ng intracranial EEG signal (solid lines) superimposed on the weights (dots) of the trained 25x25 EEG-SOM grid............40 2-9 A qualitative illustration of th e accuracy of the 25x25 EEG-SOM........................41

PAGE 10

x 2-10 SOM-SI results showing the depe ndencies between two signals. X and Y correspond to channels RTD4 and RST1, respectively............................................42 2-11 Experiment setup to compare SOM-Simila rity Indices obtained from 2 separate maps.........................................................................................................................44 2-12 Time intervals from all the seizures us ed as test data (not drawn to scale).............44 2-13 Comparing interdependencies be tween channels LOF3 and LOF4.........................47 3-1 Block diagram to extract spatio-temporal groupings information in Multivariate EEG structures..........................................................................................................50 3-2 Time profiles of the Roessler and th e Lorenz components. a) Coupling strength, C = 0 and b) C = 5....................................................................................................57 3-3 Dendrogram illustration of similarity-bas ed time series clustering. a) Coupling strength, C = 0..........................................................................................................60 3-4 Model showing the generation of linearly dependent signals using configuration 1.............................................................................................................................. ..62 3-5 Plot showing the two band pass filt ers, separated by a narrow pass band...............63 3-6 Z-planes showing the polezero confi gurations for the two band-pass filters in Fig 3-5......................................................................................................................63 4-1 Maximum average driving ability of each of the six (6) channels, nearly 100 minutes before and 70 minutes af ter Seizure-1 in patient P092...............................68 4-2 Seizure 11 of patient P093: Cluste r-similarity matrices P indicating the probability that two channels share the same cluster label in a 30 minute timeinterval......................................................................................................................7 2 4-3 Dendrogram representation of the cl uster results in Seizure 11, P093....................73 4-4 Statistical validation of the clustering re sults. In each panel, thick lines are used to represent the profiles of the three clusters in a 30 minute time interval..............78 4-5 Seizures 4 and 5 of patient P093: Cl uster-similarity matrices indicating the probability that two channels share the same cluster label in a 30 minute timeinterval......................................................................................................................8 1

PAGE 11

xi 4-6 Seizures 6 and 7 of patient P093: Cl uster-similarity matrices indicating the probability that two channels share the same cluster label in a 30 minute timeinterval......................................................................................................................8 2 4-7 Dendrograms corresponding to P092, Seizure 1......................................................83 4-8 Cluster-similarity matrices indicating th e probability that tw o channels share the same cluster label in a 30 minute time-interval........................................................85 5-1 Quantifying the spatial changes al ong time using Mantel statistics.........................99 5-2 Time smoothened profiles of the Mantel statistics during in ter-ictal states...........100 5-3 Illustrating the statistical difference be tween z-score at time‘t’ and z-score at a reference time, 90, 100 and 110 minutes before a seizure. Top to Bottom: profiles of z-score statistics for se izures 1 through 3 of patient P092....................101 5-4 Illustrating the statistical difference between z-score at time‘t’ and z-score at a reference time, 90, 100 and 110 minutes before a seizure. Top to Bottom: profiles of z-score statistics for se izures 4 through 6 of patient P092....................102 5-5 Illustrating the statistical difference be tween z-score at time‘t’ and z-score at a reference time, 90, 100 and 110 minutes before a seizure. Top to Bottom: profiles of z-score statistics for seizure 7 of patient P092......................................103 5-6 Results of synchronization analysis on seizure 1 of P092, to check the differences between the focal and non-focal hemispheric interactions, at pre and post-ictal states.......................................................................................................109 5-7 Null-hypothesis rejections rate averag ed over 10 seizures for experiment 1.........110 5-8 Null-hypothesis rejections rate averag ed over 10 seizures for experiment 2.........111

PAGE 12

xii Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy SPATIO-TEMPORAL DEPENDENCY ANALYS IS OF EPILEPTIC INTRACRANIAL ELECTROENCEP HALOGRAPH By Anant Hegde May 2006 Chair: Jose C. Principe Major Department: Electrical and Computer Engineering There is evidence that the mechanisms leading to epileptic seizures can be understood by continuously tracking the ongoing spatio-temporal mappings in the brain. However, the exact spatio-temporal change s leading to epileptic seizures, although widely studied, are not yet well understood. Pr evious studies have mostly focused on individually tracking the dynamical change s along time. However, approaches to simultaneously explain a system’s temporal ch anges in its overall spatial configuration can be much more effective in effort s to characterize epileptic events. In this dissertation, we propose a simple statistical approa ch to quantify the temporal changes in spatial patterns of an intracranial EEG. Previ ously, we developed a non-linear synchronization measure, called the SOM-Similarity Index, to quantify mutual associations between various brain regions. We propose to apply the Mantel test statistics on the SOM-similarity indices to track the temporal changes of the spatial patterns.

PAGE 13

xiii Statistical comparisons between inter-ictal and pre-ictal states suggest significant changes in the spatial connectivity prior to a seizure. Another aspect of this study investigates the regional gr oupings in th e intracranial EEG spatial networks at different ictal-states We develop a model that will allow us to study the various cluster patterns in an epil eptic brain. Studies on 10 seizures from 2 patients reveal strong connections between th e left-sub-temporal and the left-temporal depth areas. In addition, strong homologous conn ectivity is found in the orbito-frontal regions. In the third aspect of this study, we investigate the differences in the synchronization levels between the focal and the non-focal hemisphe res at pre-ictal and post-ictal states. Statistical tests confirm th e existence of significant differences at preictal states followed by a strong entr ainment, 20 minute post-ictal period.

PAGE 14

1 CHAPTER 1 SPATIO-TEMPORAL INTERACTIONS 1.1 Introduction The word “interaction” in a system can possi bly be referred to as an exchange of information between its components. A system can be understood to be made up of a large number of independent components, each acting as an individual sub-system. These sub-systems generally possess an independent environment of their own. However, in a synergetic environment, they also comm unicate or share information with their counterparts in such a way that the cumula tive effect is responsible for the overall functioning of the system. Even though the notion of a system is th e same in any discipline, the specific interpretation or representation of the same differs and is usually tailored to suit our understanding. In physics or engineering, a system is usually considered to be made up of a set of well connected oscillators that are assumed to have independent rhythms specified over periods of time and expressed as time-sequences. Harmonics induced by the oscillator rhythms encode the state of th e oscillator and perhaps the system as well. As a result of coupling with other oscillat ors in an interactive environment, a signal recorded from a transducer can be expressed in a number of ways: a) as a superposition of oscillations generated individually with osci llations shared from other oscillators or b) as a complex interactive process induced by phase, frequency or amplitude variations between oscillators in the sp atial neighborhood. Thus the pro cess where oscillators share

PAGE 15

2 complex dynamics with their counterparts ca n be characterized as a state of mutualharmony or synchronization. Referring to the degree of interacti on, it is possible that it be weak or strong, but is nevert heless capable of triggering an event. While a lot of enthusiasm is placed on understanding the exact nature of synchronization, equal importance has been devoted to relate the event manifestations caused due to synchronization changes. Finall y, synchronization analysis of the time-recordings of a physical system allows us to probe the unde rlying hidden dynamics involved in creating the system. Synchronization changes in a system can be across both space and time. Particularly, in a multi-variate system, unde rstanding the interactions among its various nodes, whose behavior can be represente d along time as time-sequences, presents numerous challenges. One of the key aspect s of tightly coupled systems with spatial extent is their ability to interact both across space and time, which complicates the analysis greatly. In biological systems such as the central nervous system, this difficulty is compounded by the fact that the component s of interest have nonlinear complicated dynamics that can dictate overall changes in the system behavior [ 22 ]. In summary, it is important to realize that the spatial structures also under go functional changes that can dictate the overall behavior of a system. With the advent of sophisticated tools, multivariate time-signal analysis has become the forefront of research in discrete time signal processing [ 66-72, 85-87 ]. Timesequences recorded from real-world system s encode huge amounts of information that are both temporally and spatially extended. Espe cially in neural systems, synchronization between cortical structures or local-field poten tials recorded at the tissue level is believed

PAGE 16

3 to be one of the prominent indicators of physiological and pathological events [ 21-22 ]. In the next section therefore, we provide a brief literature overview of the efforts in this area and then in the last part of this chapter, we try to address the funda mental question related to epileptic seizures that eventually served as a motivation for us to pursue this study. 1.2 Different Kinds of Synchronization In the previous section, we introduced the notion of coupling and abstractly described it as a process of information exchange between components of a system. Exactly how information is shared is a quest ion that remains ambiguous. Previous studies on multi-variate time series analyses have resulted in development of a wide array of signal-processing tools for quantification of coupling in systems. However, the general consensus on the best approach to quan tify this phenomenon is largely uncertain. Synchronization can be conveniently cl assified as linear or non linear [ 13 ]. Linear synchronization assumes that the processes ar e Gaussian distribute d, stationary, ergodic and therefore the couplings can be expre ssed using just the fi rst and second order moments. While some of the tools quantif y the correlations by exploiting the time parameters such as amplitude or phase of th e data, the others operate in the frequency domain to achieve the same [ 11, 60-63 ]. Recently, multi-resolution techniques have facilitated significant progress in exploiting the time and fr equency parameters jointly [ 90 ]. One of the major advantages of assuming linearity is the fact that the system interactions can be quantitatively analyzed in great depth. However, assuming a linear model can be arguably very weak and in practice, it is very restrictive to assume that the functional mapping between two coupled systems is sensitive to the differences in first and second order moments. Nonlinear coupli ng can produce changes in the higher order

PAGE 17

4 moments of stochastic signals without aff ecting appreciably the lower order moments, which raises concerns on the genera l applicability of these methods. Most real world systems are evidently non-linear processes an d their functional interactions are known to be dynamical in natu re. Unlike in the linear case, it is hard to determine the exact functional mapping and th erefore without any apriori assumptions, only the degree to which the two systems interact can be quantified [ 43 ]. This is exactly how most of non linear tools opera te. These tools, in general, assume that if two systems are related, states of one system should be described by the states of another. Even though these techniques have a wi der sense of applicability, they still suffer from the lack of analytical transparency, i. e. their inability to analyze the systems in depth. Another major setback is the huge demand in the computational complexity involved. In the literature, synchronization between systems has also been distinguished as identical synchronization, phase synchronization and generali zed synchronization. As far as the tools are concerned, they can be classifi ed into two categories, 1) bi-variate and 2) multi-variate. As stated earli er, these classifications are made to suit the hypothetical definition of what synchronization is. In the next section, we provide a brief literature overview of the various signal processing t ools designed to quan tify the degree of synchronization between systems. 1.3 Second Order Synchronization Measures The literature on quantification of the spatio-temporal couplings in multidimensional dynamical systems is rich and a bundant, and for simplicity, it can be grossly divided into linear and non linea r. This section reviews some of the widely used second order linear tools in this effort. We will discuss the common goals among these measures, their approaches and the advantag es and trade-offs along the way.

PAGE 18

5 1.3.1 Cross-Correlation Cross-correlation analysis is one of the earliest and mo st relied-on linear techniques [ 6 ]. It measures the linear coupling be tween two signals and by design, is symmetric. Consider an identically distributed sta tionary stochastic process X. Assuming ergodicity, represent a single time recording of the process as the discrete time-signal x(n) n = 1,2, …. M where M is the number of total number of samples in the signal. Cross-correlation between tw o discrete time signals x(n) and y(n) can be estimated as a function of lags, as follows: 1 ,..... 1 0 1 ,..., 1 ) ( 1 ) ( ˆ1 N N y x N RN i i i xy (1.1) Due to the symmetric property, ) ( ˆ ) ( ˆ xy xyR R If y x are the mean-estimates of x(n) and y(n) respectively, the cross-covariance between x(n) and y(n) can be estimated as: y x xy xyR C ) ( ˆ ) ( ˆ (1.2) Again, it is obvious that ) ( ˆxyC is symmetric. Cross-correlation coefficient between x(n) and y(n) can now be defined as the cro ss-covariance normalized by the product of the square root of the variances of the two signals. Mathematically, y x xy xyC ˆ ˆ ) ( ˆ ) ( ˆ (1.3) where y x ˆ ˆ are the standard-deviations of x(n) and y(n) respectively. xy ˆ is normalized between -1 and +1. The bounds repr esent strongest correlation (positive and negative) while 0 represents independence.

PAGE 19

6 Similarly, define )) ( ( ) ( n x F X and )) ( ( ) ( n y F Y as the Fourier transform equivalents of x(n) and y(n) respectively. If the cross-spectrum ) ( ˆxyC and the autospectrums ) ( ˆxxC and ) ( ˆyyC are defined accordingly, normalized cross-coherence can be mathematically represented as ) ( ˆ ). ( ˆ ) ( ˆ ) ( ˆ w C w C w Cyy xx xy xy (1.4) The cross-coherence quantifies the degr ee of coupling between X and Y at each frequency. Again, the bounds are the same as in cross-correlation coefficient. One of the problems with cross-correlati on and the cross-coherence measures is that they are symmetric. Symmetry fails to quantify the directi onal information and therefore allows us to know only the averag e amount of information exchanged between X and Y. Also, as with most other technique s, pair-wise computation is a limitation with correlation/coherence analysis as well. 1.3.2 Partial Directed Coherence (PDC) Efforts to improve on the correlation t echnique resulted in development of multivariate auto-regressive modeling techniques such as directed coherence, directed transfer functions (DTF) [ 19-20 ], and partial directed coherence (PDC) [ 2-3, 91 ]. It is a frequency domain technique based on the con cept of Granger-causali ty [27] which says that an observed time series xj(n) causes another series xi(n), if knowledge of xj(n)’s past significantly improves prediction of xi(n) Along the same lines, the reverse case may or may not be true. One way to make a quantit ative assessment of th e amount of linear interaction and the directi on of interaction among multiple time-series would be by translating the concept of Granger-causality into a mathematical model such as a multi-

PAGE 20

7 variable autoregressive model (MVAR). This will enable us to check the number of predictable variations arising from using the knowledge of the past samples of another time-series to predict the current sample of the time-series of interest. The partial directed coherence (PDC) from j to i at a frequency f is given by ) ( ) ( ) ( ) ( f a f a f A fj H j ij ij (1.5) where p r fr j ij p r fr j ij ijotherwise e r a j i e r a f A1 2 1 2, ) ( ) ( 1 ) ( (1.6) ) (r aij are the multivariate auto-regressive (MAR) coefficients at lag r obtained by finding the least-squares solution of the MAR model n r p x A xp r r 1 0) ( (1.7) Here, r NN r N r N r N r r ra a a a a a A . . . . .2 1 1 12 11 and ) ( . ) ( ) ( ) (2 1r p x r p x r p x r p xN p represents the dept h of the AR model, r represents the delay and n represents the prediction error or the noise. Note that ) (fij represents the relati ve coupling strength of the interaction of a given signal source j with regard to signal i as compared to all of j’s interactions to other si gnals. It turns out that the PDC is normalized between 0 & 1 at all frequencies. When i=j, the PDC represents the influence of a signal’s past on its current state. The PDC also differentiates between direct and indirect interactions among

PAGE 21

8 multiple time series. Evidently, the MVAR mo dels are better than most bivariate measures because they use information simultaneously from all the channels, and thus are able to unambiguously distinguish between di rect and indirect causal connectivity between nodes [ 41 ]. 1.3.3 Synthetic simulation Suppose that three simultaneously observed time-series are generated as follows ) ( ) 1 ( ) 4 ( 36 0 ) ( ) ( ) 1 ( 1 0 ) 2 ( 75 0 ) ( ) ( ) 2 sin( ) 2 sin( ) 2 sin( 1 0 ) (3 1 2 3 2 2 1 2 1 3 2 1 1n w n Cx n x n x n w n x n x n x n w n f n f n f n x (1.8) where f1=0.1fs f2=0.2fs, f3=0.44fs sampling frequency fs=100Hz; Coefficient C denotes the coupling strength between x1 and x3. w1, w2 & w3 are zero-mean uncorrelated white noise processes with identic al variances equal to 0.5. One can observe from (1-5) that x2 is influenced by x1 and x3 is influenced by x2 and can possibly be influenced by x1 as well, depending on the coupling strength C. Therefore, x3 can be influenced by x1 directly or indirectly or both directly and indirectly. The coupling state diagram is plotted in fig. 1-1. Figure 1-1. Coupling diagram showing the di rection of interactions between nodes, x1, x2 and x3.

PAGE 22

9 (a) (b) Figure 1-2. Matrix layout plots for PDC, describing th e interactions in example simulation (1-8). ij indicates the influence of xj on xi at frequency f. The xaxis limits range from 0 to pi, and the y-axis limits range from 0 to 1. a) Case when C = 0.3 and b) C = 0. Corresponding PDC matrix plots for the case when C = 0.3 are shown in fig 1-2a. The plots clearly suggest a dependency of x2 and x3 on x1. Dependency of x3 on x2 is also prominent from the plot. All the other plots (excluding the diagonal pl ots) indicate some minor influences in the reverse direction, i.e, x2 to x1 etc. However, a suitable statistical threshold may clearly indi cate an unambiguous influence among the nodes. The PDC

PAGE 23

10 results seem to be in agreement, at least qualitatively, with the c oupling diagram in Fig 11. For the case when there is no direct coupling between x1 and x3, i.e, C=0, the plots are described in fig. 1-2b. In this case, PDC indicates a much lesser influence of x1 on x3, as expected from its design. Qualitatively, the PDC plots are reasonably good indicators of directional dependence between multi-variate structures. Ho wever, it is difficult to exactly mark the directions without setting an appropriate th reshold. Quantitatively, the information on the exact amount of coupling between two signals at a particular frequency f, as determined by PDC, is ambiguous. In the simulation example, we created x2 and x3 using direct and indirect influences from x1 Since x1 is the driver and also is composed of multiple sinusoids (f1, f2 & f3), it is natural to expect the frequency-spectrum of x2 and x3 to mainly consist of these sinusoidal frequencies. This also implies that the coherence between these time-structures will be strong at these pole locations. However, the PDC plots fail to provide that information, as we can s ee from fig.1-2a and 1-2b. The coherence plots are either smeared across all the frequencies or exhibit high values at frequencies other than the expected frequency locations. The MVAR approaches have been used to determine the propagation of epileptic intracranial EEG activity in tempor al lobe and mesial seizures [ 2-3, 10, 19-20 ]. However, these models strictly require that the measur ements be made from all the nodes, or the directional relationships coul d be ambiguous. In addition, there is no clear evidence of causality relationships among the cortical struct ures and as suggested in Rodriguez et al. [ 84 ], the nature of synchronization is mostly instantaneous or w ithout any detectable

PAGE 24

11 delay. The applicability of MV AR measures is constrained to the linear dependencies even though there are strong evidences point ing to non linear pattern s of interactions [ 23 ]. 1.4 Mutual Information Linear measures are typically restricted to measure statistical dependencies up to the 2nd order. If signals are Ga ussian distributed, the 2nd order statistics are sufficient to capture all the information in the data. Howeve r, in practice, most signals are either nonGaussian or quassi-Gaussian rendering the lin ear statistical measures inadequate. As an alternative, information theory measures [ 11-12, 66-68, 94-95 ] have been widely claimed as acceptable choices to exploit the higher orde r statistical dependencies between signals. A useful quantity is the mutual information approach that measures the Kullback-Liebler (K-L) divergence between the joint probabi lity and the product of the marginal probabilities [ 11 ]. Consider that the vector samples x = {x1,….,xN} from an information source in a ddimensional signal space. Suppose that these samples are drawn from the distribution f(x). Similarly, let the vector samples y = {y1,….,yN} form another information source in a d-dimensional signal space belongs to distribution g(x). By definition, mutualinformation between x and y can be computed as follows ) ( ) ( ) ( ) ( ) ( ) ( log ) ( ) ( ) ( ) ( log ) ( ) ; (, ,y x H y H x H dy y f dx x f dxdy y x f y x f dxdy y f x f y x f y x f y x Iy xxy y x xy xy y x y x xy xy (1.9) where ) ( y x fxy is the joint probability event between x and y and ) ( ), ( y f x fy xare the marginal probabilities, respectively. H(x) can be regarded as the average measure of

PAGE 25

12 uncertainty and in the informati on theory literature is also called as Shannon’s marginal entropy, in x. Similarly H(y) is Shannon’s marginal entropy in y, while their joint entropy is denoted by H(x,y) MI has been successfully used to quan tify statistical coup lings in biological applications such as newbor n cardio-respiratory systems [ 70-72 ], event related potentials (ERPs) [ 65 ] and epileptic seizures [ 66-68 ] One of the problems with this approach, however, is that it requires large data sets for probability density estimation making it computationally unviable. Effi cient algorithms such as gene ralized mutual information function (GMIF) have been recently proposed [ 70-72 ] to enable speedy computation using very small data sequences. 1.5 Phase Synchronization Conceptually, if the rhythms of one signal ar e in harmony with that of the other, the two signals are known to be phase locked [ 99 ]. Phase synchronizati on can therefore be defined as the degree to which two signals are phase locked. Most of real world signals have broad sp ectra. For example, EEG signal recordings are usually in the range of 0.1 to 1000 Hz even though they are usually band pass filtered between 0.1 and 80 Hz since a major portion of the energy is contained in that spectrum. The EEG can be classified r oughly into five (5) different frequency bands, namely the delta (0-4 Hz), theta (4-8 Hz), alpha (812 Hz), Beta (12-16 Hz ) and the Gamma (16-80 Hz) frequency bands. Freeman [ 21 ] demonstrated evidence of phase locking between EEG frequency bands across diffe rent regions of the brain, leading to certain clinical events such as evoked potentials. Similarly, it is also believed that phase synchronization across narrow frequency EEG bands, pre-se izure and at the onset of seizure [ 81-83 ] may provide useful hints of the spatio-tempor al interactions in epileptic brain.

PAGE 26

13 The procedure normally adopted to measur e phase synchronization consists of initially tracking the instantaneous phases of the two signals of interest. This is followed by computing their differences called the relative phase. Qu antifying the distribution of the relative phase will yield an index describing the degree to which the two signals are phase locked. Hilbert transform [ 33, 49, 75, 81-83 ] is a widely used signal-transformation trick to compute the instantaneous paramete rs of a time-signal. The following subsection briefly describes the steps involved. 1.5.1 Hilbert Transformation Consider a real-valued narrow-band signal) (t x concentrated around frequency fc. Define ) ( ~ t x as t t x t x1 ) ( ) ( ~ (1.10) ) ( ~ t x can be viewed as the output of the filter with impulse response t t h1 ) ( t (1.11) excited by an input signal ) ( t x. We call this filter a Hilbert transformer. Next, we construct an analytical signal ) ( ~ ) ( ) ( t x j t x t z (1.12) which contains only positive frequencies. Our objective is to extract the instantaneous parameters of the signal such as instanta neous amplitude, instan taneous frequency and the instantaneous phase. z(t) can also be expressed in Cartesian co-ordinates as ) () ( ) (t je t A t z (1.13) where A(t) is the instantaneous amplitude and (t) is the instantaneous angle of the signal x(t). (t) can be expressed as )) ( ) ( cos( ) (0t t t (1.14)

PAGE 27

14 ) ( t is the instantaneous pha se of the signal and ) (0t is the instantaneous frequency of the signal x(t). Hilbert transforms are accurate only when the signals have a narrow-band spectrum, which is unrealistic for most real-world signals. Pre-processing of the signal such as decomposing it into narrow frequency bands is needed before we apply Hilbert transformation to compute the instantaneous parameters. There are several ways of breaking a signal into its sub-ba nd components. Previous work [ 48, 76-80 ] to achieve signal decomposition either points to a bank of digital FI R filters or the sub-band coding using wavelet transforms. In the former, a multitude of conventional FIR filters are designed to cater to the desi red pass-bands using the Hamming or the Hann window. It is well known from the signal processing literature that the transient events of the signal, such as the spike activity and the sudden transi tions that quite often exist in real signals, are not adequately captured using the convent ional FIR filters. In other words, it is essential to capture the changing events in the signal locally at the time of their occurrence. If the signal is non-stationary, the requirement for locality is certainly not met with conventional FIR f ilters. The Wavelet filters [ 101 ] have shown better capability of representing a signal in both the time and frequency domain and also capture the transient details of the signal more effectiv ely than the conventional digital FIR filters. However, in both the techniques the decom position does not always lead to narrow-band signals. In other words, the decomposed signa ls do not conform to the definitions of a narrow-band signal, even t hough they exist in a very narrow band frequency range. Certain conditions need to be met to define a meaningful instantaneous frequency on a narrow-band signal. They are as follows:

PAGE 28

15 the signal needs to be symmetric with respect to the local zero mean, and the signal should have the same numb er of zero crossings and extrema. 1.5.2 Empirical Mode Decomposition (EMD) EMD is another preprocessing technique [ 34 ] in which the data are decomposed into narrow band components called the intr insic mode functions (IMF). Consider a broad-band real-valued signal s(t) Using an iterative sifting process [ 34 ], the IMF components are extracted from the signal s(t) iteratively. The first step in the sifting process involves identifying the lo cal minima and the maxima of s(t) and then using an interpolation filter to combine the extremas. The local mean of the signal m1 is computed using the local extrema information and subtracted from s(t) i.e, 11 1) ( h m t s (1.15) The new signal is checked for the above 2 conditions. If it did not satisfy them, the whole sifting process is repeated on h11 to obtain say h12. This sifting process is iteratively repeated until an IMF component, say h1, satisfying the two criteria is obtained. The whole sifting process is now repeated on the residue signal s1(t) obtained from subtracting the IMF component h1 from the original signal s(t) i.e, 1 1) ( ) ( h t s t s (1.16) Eventually, all the successive IMFs satis fy the narrow band constraints. The second condition is particularly essent ial in order to ensure that the instantaneous frequency does not have any undesired fluctuations, introduced by asymmetric wave forms. Contrary to previously discussed methods, the EMD derives IMFs or the basis f unctions directly and adaptively from the data. This characteristic ma kes it particularly useful to detect local

PAGE 29

16 changes in the signal. The basis functions can be empirically shown to be orthonormal and complete. We demonstrate the ability of the EMD to decompose a signal into its narrow-band components by presenting a simple syntheti c simulation. Consider a non-stationary stochastic process X ( t ) consisting of a combination of sinusoidal components and noise (constructed as shown in fig. 1-3). Components x1 and x3 consist of 100 Hz ( f1) and 450 Hz ( f3) sinusoid oscillations respec tively and span over the en tire duration of the signal X ( t ). x2 is a sinusoidal component of frequency 260 Hz ( f2) spanning only for some period of the entire signal’s duration. The signals were sampled at 1024 Hz ( fs) frequency. Decomposition of X (t) into its narrow band components using EMD is shown in fig 1-4. The IMF components (fig 1-4a) and their corresponding frequency spectrum (fig 14b) reveal a few things and can be listed as follows 1. All three frequencies are well separated by the EMD and also the middle panel in fig 1-4a shows how the frequency compone nt f2 exists only between the timeduration that it was created in. 2. There exists leakage of f2 and f3 into IMF1, but by looking at the power spectrum, it is easy to observe that the energy at thes e frequencies is very low compared to f1. 3. Peak to peak amplitude of the extracted IM Fs is almost half of the peak to peak amplitudes in the original signal; however, th is is just a scaling problem and can be easily resolved by scaling up each of the IMFs. 4. The Hilbert spectrum, computed using the information from the instantaneous parameters of the IMFs, reveals the joint time-frequency structure of the signal. The extraneous frequencies are either created during the first transient period or due to the presence of noise. Energy at these frequencies is below 0dB.

PAGE 30

17 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -2 0 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -2 0 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -2 0 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -5 0 5 X = x1 + x2 + x3 + w x1 x2 x3 Figure 1-3. Synthetic Sinusoida l signal X(t) and its compone nts x1, x2 and x3. X(t) also consists of noise w~N(0, 0.1) 0 0.5 1 -2 0 2 0 500 1000 0 200 400 0 0.5 1 -2 0 2 0 500 1000 0 200 400 0 0.5 1 -2 0 2 0 500 1000 0 200 400 IMF1 IMF2 IMF3 time (sec) frequency (Hz) IMF components Frequency spectrum (a) (b) Figure 1-4. Decomposition of sinusoids us ing EMD. a) IMF components obtained from successive decomposition of x(t) using the EMD b) Corresponding frequency spectrum of the IMF components. IMF1 predominantly consists of 450 Hz, IMF2 consists of 260 Hz and IMF3 consists of 100 Hz, as desired.

PAGE 31

18 1.5.3 Wavelet decomposition As described in previous sections, one of the many ways to achieve narrow-band decomposition is by using ortho-normal wavelet filters [ 101 ]. For the same synthetic signal X (t), we apply the ‘7 length Daubechies’ filter functions, at various scales, to achieve signal decomposition (fig 1-6). Figure 1-5. Hilbert spectrum constructed fr om the IMF functions. The spectrum reveals the presence of 100 Hz and 450 Hz along the entire signal duration. Presence of 260 Hz (corresponding to 0.26 on the Y-axis) component between the 200th and 800th sample (corresponding to approximately 0.2 and 0.8 second) is also clearly visible from the spect rum plot. The leakag e frequencies are mainly due to the transient oscillations created in the first few samples of the IMFs and also due to noise. The following observations indicate the superior ability of EMD compared to wavelets in a non-stationary scenario. Two spurious frequency components of reasonably high energy, namely 150 Hz and 240 Hz are created by wavelet analysis and synthesis (recons truction) filters. 5. Qualitative comparison of Hilbert spectrum plots (fig 1-7 and fig 1-5) suggest a larger smearing of the frequencies in wavelets compared to EMD.

PAGE 32

19 6. The amount of energy leakage from the 260 Hz component into the first panel (fig. 1-4b) is considerably highe r than the corresponding am ount of energy leakage in IMF (first panel in fig. 1-4b). 0 0.5 1 -5 0 5 0 500 1000 0 500 1000 0 0.5 1 -5 0 5 0 500 1000 0 500 1000 0 0.5 1 -5 0 5 0 500 1000 0 500 1000 Wavelet reconstruction Frequency spectrum W1 W2 W3 Figure 1-6. Decomposition of sinusoids us ing wavelets. Left) Wavelet reconstructed components obtained from using Daub echies 7-scaled filters. Right) Corresponding frequency spectrum of th e wavelet reconstructed signals. W1 predominantly consists of 450 Hz, but also has a major portion of 240 Hz. W2 has more power concentrated in 260 Hz and 150 Hz frequencies and IMF3 consists of 100 Hz, as desired. Presence of spurious frequencies such as 240 Hz and 150 Hz are revealed using wavelet-decomposition. 1.5.4 Indexing Phase Synchronization As we saw earlier, the Hilbert transforma tion is a useful time domain to a timedomain transformation technique to compute th e instantaneous phase of a signal. Due to the narrow band constraints, any wide-band signal needs to be broken down into its narrow band components and we saw that EM D and wavelet filters can be used as potential pre-processors to achieve the same. Once a meaningful estimate of the instantaneous phases is obtained, we need to determine the degree of phase locking by quantifying them with a suitable metric.

PAGE 33

20 Figure 1-7. Hilbert spectrum constructed fr om the wavelet recons tructed narrow-band signals. The spectrum reveals the pr esence of 100 Hz and 450 Hz along the entire signal duration. Frequencie s in the range of 200 Hz and 350 Hz (corresponding to 0.26 and 0.35 on the Y-axis) are smeared between the 200th and 800th sample (corresponding to approximately 0.2 and 0.8 second). Presence of 260 Hz is not very clear from the spectrum plot. Let ) ( tx denote the instantaneous phase of one of the narrowband components of x(t) Similarly, let ) (ty be the instantaneous phase of the corresponding narrowband component of y(t). Relative phase (t)can be computed as the difference between ) (tx and ) (ty i.e, ) ( ) ( ) (t t ty x t = 1, 2,……, N. The relative phase or the phasedifference (t) forms a distribution and a number of metrics can be used to compute an index of phase synchronizati on by using the distribution of (t). Two popular metrics exist in the literature [ 60-63, 81-83, 85, 99 ]. The first one expresses the phase synchronization index as a devi ation in the distribution of (t) from the uniform distribution diff [ 85, 99 ]. If we use the Shannon’s defi nition for entropy, comparison of the entropy H of (t) with entropy of a unifo rm phase distribution Hunif will result in, unif unifH H H (1.17)

PAGE 34

21 Note that 0 for complete lack of synchronizati on and is 1 for perfect phase-locking. Another metric commonly used synchronization index is the phase locking value (PLV), also called as mean phase coherence [ 81-83 ]. It is given by, t t je) ( (1.18) The bounds are the same as in previous index. EMD is a useful pre-proce ssing technique for decomposi ng a time-series into its narrow band components. However, it does not provide any design parameter using which a signal can be filtered into specific bands of intere st. In phase synchrony analysis where instantaneous phase is used as the synchrony metric to evaluate phase-locking between two time-series, it is very important that the narrow band components of the time-series are in sufficiently narrow (such as 2 to 4 Hz pass band range) spectra. While it is possible to apriori design a conventional FIR filter catering to a desired pass band and stop band, the EMD does not provide the luxury of choosing the desired frequency bands. Due to the fact that EMD adaptively de rives its basis functions based on the data, the IMFs of two different signals can easil y lie in different bands. Hence, another limitation associated with using EMD is that it is not possible to compare between the IMFs of two signals. Unlike the MVAR model based approaches the applicability of phase-synchrony measures is restricted to bi-variate time-series [ 10 ]. We also saw that both the synchronization indices are symmetric and are therefore incapable of providing the directionality of dependency. Also they are restricted to id entifying only the phase locking between time-signals. Synchroniza tion between two systems is a general phenomenon and can broadly be defined as the functional relati onship governing the

PAGE 35

22 oscillation of the two systems [ 88-92 ]. The function could be lin ear invertible, linear noninvertible, non linear invertible and non lin ear non-invertible. Ph ase synchronization, being only a subset of the general synchroniza tion definition, fails to address the overall functional relationship between systems. 1.6 Non-Linear synchronization measures A lot of studies in the last two decad es, to quantify nonlinear dependencies [ 88-93 ] between two or more signals have met with reasonable success. As pointed out earlier, most real world systems can be ch aracterized as non linear processes [ 52-54 ] with their samples dynamically evolving from interacting w ith past samples of the same process as well as from the samples in the spatial ne ighborhood. Functional non linear mappings are generally very hard to charac terize and the fact that the recordings are embedded in noise makes it harder. A common procedure in any non linear measur e is to construct a system map (or the Poincare map) by embedding a scalar time series into a high dimensional Euclidean space and building a trajectory by connec ting the time points in that space [ 98 ]. This is followed by computations that are generall y directed towards determining the mean distances between the mutual neighbors in the phase space of the signals. The synchronization indices indicate the degree to which the state of one system is a function of another and vice-versa. 1.6.1 Signal Reconstruction Taken’s embedding theorem states that unde r certain conditions, it is possible to construct a finite dimensional dynamical system from a scalar time series ,1N t tx where xt= g(ut). ut is the actual dynamical system, residing in dimension say ‘d’. Taken [ 98 ]

PAGE 36

23 postulated that, using time-delay embedding the exact dynamics of a system can be reconstructed, while also preserving its topol ogical characteristics. The reconstruction of the attractor is done from a finite time seri es of the observation of a single variable x(t). Thus a multi-dimensional embedding space is constructed from the time series data, and a point in it represents the state of the system at a given time. ) ) 1 ( ( ) ) 2 ( ( . ) 2 ( ) ( ) ( d t x d t x t x t x t x xn n n n n n (1.19) Note that time-delay reconstruction relies heavily on the embedding parameters such as the embedding dimension (d) and the delay time (). Nevertheless, these parameters are important since they convey a lot of inform ation about the dynamical properties of the underlying processes. The embedding dimension d can be selected by the false-nearest neighbor criteria [ 98 ] and the delay time can be estimated by the method of mutualinformation [ 98 ], or by measuring zeros of autocorrelation. 1.6.2 Non-Linear Synchronization Measures Eckmann et al. [ 17 ] proposed the method of recurren ce plots (RPs) that represents the recurrence of states in the phase-space tr ajectory of a chaotic signal. Since chaotic systems are non linear in nature, this method has been fairly successful in detecting bifurcations and non-stationar ities in time sequences. Cross -recurrence plot (CRP) is an extension of the RP idea to multi-dimensional time signals [ 24 ]. The CRP has found applicability in describing the time-depe ndency between multiple time-series recorded

PAGE 37

24 from multiple locations. However, the lack of quantitative information and the computational complexity makes it tedious for analyzing large sets of data. One of the other drawbacks of this measure is that it fa ils to indicate the direction of information flow (or influence). Variants of recurrence plots [ 59 ] are used to measure recurrence of states in the phase-space betw een two chaotic signals. However, these plots are difficult to analyze due to the lack of quantitative measures, in addition to their computational complexity. Several other non linear measures were proposed following the definition of generalized synchronization by Rulkov and co-workers [ 88-92 ] Mutual predictabilities were defined and studied extensively by Schiff et al. [ 92 ] and Le Van Quyen et al. [ 81 ]. Principally, if two systems are functionally re lated, observing the states of one system should help predict the states of the other. A statistic constr ucted on this idea will indicate the relative strength of coupling as well as the direction of co upling between any two systems. Recently, Arnhold et al. [ 1, 76 ] introduced the similarity–index technique (SI) to measure such asymmetric dependencies be tween time-sequences. Conceptually, this method relies on the assumption that if ther e is a dependency between two signals, the neighboring points in time will also be neighbor ing points in state space. In other words, if there is a functi onal dependency between two signals, then the recurrence of dynamics of one signal will mean the recurrence of dynamics in the other signal. Arnhold et al. [ 1, 76 ] propose to search for recurrence in the signals’ state-space, which poses an enormous computational burden, especially for large data sets. In this dissertation, we propose a self -organizing map (SOM) based improvement to this method

PAGE 38

25 to reduce computational complexity, while maintaining accuracy. This is achieved by mapping the embedded data from signals onto a quantized output space through a SOM specialized on these signals, and utilizing the activation of SOM neurons to infer about the influence directions between the signals in a manner similar to the original SI technique. The approach not onl y facilitates real-time comput ation but it also enables us to derive long-range synchronization patterns. 1.7 Objectives and Author’s Contribution Mechanisms involved in the synchroni zation among neural populations have received a lot of attention in recent years [ 4-5, 8-9, 14-16, 42, 60-63, 97 ] It is agreed that even in the normal physiological states, the neuron ensembles at the cortical structures always interact and undergo a certain de gree of synchronization. The Response to sensory stimulations in the form of cohesive firing of neurons is known to lead to certain transient changes. A larger question, however, that is still debated upon is whether elevated synchrony in firing in a spatial neighborhood and its subse quent propagation to other areas in the brain can l ead to pathological situations such as seizures, stroke or Alzheimer’s. Epileptic studies on adults and neo-na tes have been widely researched [ 26, 35-40, 60-62, 82 ]. Researchers from different discip lines have contributed immensely in understanding the various neurol ogical aspects leading to seizures. Synchronization and de-synchronization between the co rtical networks are believed to be one of the plausible reasons for this neurological disorder. It is soon becoming clear that epilepsy is a dynamical disease [ 56 ], i.e. the macroscopic spatio-temporal dynamical across different regions of the brain are consistent with rapi d, sometimes gradual and often very subtle nonlinear dynamical interactions. It is believed that synchr onization occurs due to both

PAGE 39

26 local and global discharges of the neurons. Fr om the epilepsy perspe ctive, quantifying the changes in spatio-temporal interactions c ould potentially lead to the development of seizure-warning systems. This quantification wo uld also help us identify the regions that actively participate during epileptic seizures. In this dissertation, we apply the SOMSI measure on epileptic seizure data to investigate i) the temporal e volution of dependencie s at different stages of a seizure, ii) temporal changes in spatial connectivity of si gnals as a seizure is approached and iii) spatial-temporal clustering of channels at different stages of seizure. In Chapter 2, we present the SOM-based SI measure to quantify spatio-temporal patterns. The SOM based SI approach to detect spatial interdependencies can be generalized to any multi-variate dynamical syst em. Initially we present this idea from a general perspective and then proceed to de monstrate the robustness brought about by the SOM-improvisation through synthetic simulatio ns. In chapter 3, we propose a spatiotemporal clustering model that effectively u tilizes similarity indices to achieve timeseries clustering. From an application standpoi nt however, the goal of this thesis is to identify spatio-temporal connectivity’s in epile ptic seizures. In chapter 4 we present the application of our model to epileptic intr acranial EEG data and analyze the clustering results on a large set of seizures. In chapte r 5, we investigate the hypothesis that the spatial patterns in intracranial EEG undergo significant changes as the brain progresses towards a seizure state. This chapter also in vestigates if the depe ndencies in the focal hemisphere are significantly different from the non-focal hemisphere, at pre-seizure and post-seizure states. Chapter 6 presents a de tailed discussion on conclusions and possible directions for future research.

PAGE 40

27 CHAPTER 2 SELF-ORGANIZING-MAP BASED SIMILARITY INDEX MEASURE 2.1 Introduction In this chapter, we will mainly discuss the Similarity Index (SI) measure proposed by Arnhold et al. [ 1 ] and the proposed SOM based impr ovement to the SI measure. Simulations with synthetic data are used to demonstrate the appli cability and robustness of the proposed SOM based SI measure. 2.2 The Similarity Index Technique Synchronization between two signals X and Y is understood, in principle, as a functional mapping between X and Y. In most neuro-biological systems especially, it is realistic to expect the functi onal mappings to be nonlinear. The attractors of functionally synchronous systems are related such that, the trajectory of one system partly describes the trajectory of another and vice-versa. In unilateral coupling, if the trajectory of Y is influenced by the trajectory of X, then Y is said to be functionally dependent on X. However, in instances where bi-directional relationships exist, the trajectory of Y can also partly influence the trajectory of X. In such a case, it is important to quantify two aspects of synchronization: direction and strength. Using the principles of generalized synchronization defined by Rulkov [ 88 ] and, Arnold et al. [ 1 ] recently formulated a synchronization measure to characterize functional relationships between no n linearly coupled systems. Assume that X and Y are two time series generated by a system, which are embe dded into two vector signals in time using delays. N(X|Y) is defined as the average dependency of X on Y and it can be written as,

PAGE 41

28 1 0) ( ) | ( ) ( 1 ) | (N n n n nX R Y X R X R N Y X N (2.1) where Rn(X) is the average Euclidean distan ce between the state-vector of Xn and the remaining state-vectors in X. Y-conditioned Euclidean distance Rn(X|Y) measures the average Euclidean distance between Xn and the vectors in X whose corresponding timepartners are the k-nearest neighbors of Yn. This measure takes va lues in [0, 1], where 0 implies no coupling and 1 implies perfect synchr onization. In other words, 1 suggests that recurrence of a state in Y implies a recurrence in X [ 1, 7, 8, 45, 78 ]. On the same principles, 0 ) | ( Y X N implies complete independence between X and Y. By design, SI can quantify nonlinear dependencies. Similarly, it is possible to quantify the average dependence of Y on X by 1 0) ( ) | ( ) ( 1 ) | (N n n n nY R X Y R Y R N X Y N (2.2) Comparing N(X|Y) and N(Y|X), we can determine which signal is more dependent on the other. The difficulty with this approach is that at every time instant, we must search for the k nearest neighbors of the current embedded signal vectors among all N sample vectors; this process requires O(N2) operations. This high complexity hinders real-time implementation and analysis. In addition, the measure depends heavily on the free parameters, namely, the number of n earest neighbors and the neighborhood size The neighborhood size needs to be adjusted every time the dynamic range of the windowed data changes. 2.3 SOM-based Similarity Index (SOM-SI) The SOM-based SI algorithm is aimed at reducing the computational complexity of the SI technique. The central id ea is to create a statistically quantized representation of the dynamical system using a SOM [ 74, 89 ]. A SOM is a neural-network in which spatial

PAGE 42

29 patterns from the input-space are mapped onto an ordered output space consisting of a set of neurons, called processing elements (PE). Thus each neuron in the SOM, based on its location on the map, compactly models diffe rent features/dynamics of the input. For best generalization, the map needs to be trained to represent all possible states of the system (or at least with as much variation as possible). As an example, if we were to measure the dependencies between EEG signals recorded from different regions of the brain, it is necessary to create a SOM that represents the dynamics of signals collected from all channels. The SOM can then be us ed as a prototype to represent any signal recorded from any spatial lo cation on the brain, assuming that the neurons of the SOM have specialized in the dynami cs from different regions. One of the salient features of th e SOM is topology preservation [ 74, 89 ]; i.e., the neighboring neurons in the feature space corres pond to neighboring states in the input data. In the application of SOM modeling to the similarity index concept, the topology preserving quality of the SOM feature of the SOM will be of added advantage, because of the fact that the neighboring neurons in the feature space will now correspond to neighboring states in the input data. Assume that X and Y are two time series genera ted by a system, which are embedded into two vector signals in time usi ng delays. Define the activation region of a neuron in the SOM as the set of all input vectors (the embe dded signal vectors) for which the neuron is the winner based on some distan ce metric (Euclidean in most cases). Let Xn be the set of time indices of input vectors xj that are in the activation region of the winner neuron corresponding to the input vector xn at time n. Similarly define the set Yn. Then the procedure to estimate the directed SOM-SI between X and Y is as follows

PAGE 43

30 7. Train a SOM using embedded vectors from both X and Y as the input. 8. At time n, findx nW, the winner neuron for vector xn, and find y nW, the winner neuron for vector yn. 9. To find Rn(X), compute the average Euclidean distance between x nW and all the other winner neurons in the SOM. Similarly, compute Rn(Y). 10. Determine the sets Xn and Yn for x nWand y nW, respectively. 11. Determine the nearest neurons y j nW, corresponding to vectorsjy, wherenX j Determine the nearest neurons x j nW, corresponding to vectorsjy, wherenY j 12. Calculate|| || ) / 1 ( ) | (, 1 x j n x n q j nW W q Y X R where q is the number of elements in Xn. Calculate|| || ) / 1 ( ) | (, 1 y j n y n q j nW W q X Y R where q is the number of elements of Yn. 13. Compute the ratios, ) ( / ) | ( ) ( ) | ( X R Y X R X R Y X Nn n n n (2.3) ) ( / ) | ( ) ( ) | ( Y R X Y R Y R X Y Nn n n n (2.4) 14. Find interdependencies ) | ( Y X N and ) | ( X Y N as the average of ) | ( Y X Nn and ) | ( X Y Nn over all n. 15. Compute the SOM-SI as the difference, ) | ( ) | ( Y X N X Y N Positive values of indicate that influence of X on Y is more than the influence of Y on X, while negative values indicate th e opposite. Higher magnitude of indicates a stronger coupling of the signals [ 29 ]. The computational savings of the SOM appr oach is an immediate consequence of the quantization of the input (signal) vector space. The nearest ne ighbor search involves O(NM) operations as opposed to O(N2) in the original SI, where M is the number of PEs.

PAGE 44

31 Traditionally M<
PAGE 45

32 Lorenz system (for different C values) ar e shown in Fig. 2-1. Each SOM is an 8x8 rectangular grid, and is tr ained on a set of 4000 samples using a Gaussian neighborhood function for about 400 iterations. The neighbor hood radius (standard deviation of the Gaussian neighborhood function) is exponentially annealed star ting from an initial value of 4 with a time constant of 100. The step size is also annealed exponentially from 0.08 using the same time constant. Using the SOM-SI approach, the normalized indices are calculated for coupling strengths of C = 0, 0.5, 1, 2, 3 and 5. The robustne ss to varying SNR is demonstrated in the test data [ 31 ], as shown in Fig 2-2. Firstly, we observe that the asymmetric difference in the interdependencies between the Rossele r and the Lorenz system increases as the coupling (C) is increased (irrespective of th e noise level). Despite the fact that the training was done with 30dB SNR level, the ab solute values of th e dependencies do not change much for SNR up to around 20dB, i ndicating noise-robustness. However, a drastic decrease in the absolute dependen cy values for 10dB and 0dB SNR situations suggest that the structural relationship between the tw o systems is being largely overwhelmed by the stochastic noise. In summar y, we see that the idea of SOM modeling of system’s dynamics not only reduces com putational complexity but also preserves neighborhood mappings during noi sy perturbations, thus providing enhanced noise robustness.

PAGE 46

33 (a) (b) (c) (d) (e) (f) Figure 2-1. Phase-space traj ectories of the Rosseler-Lorenz system for various coupling strengths a) Rosseler b) Lorenz (C=0) c) Lorenz (C=0.5) d) Lorenz (C=2) e) Lorenz (C=3) f) Lorenz (C=5). The SO M weights (circles) for each signal are superimposed on the trajectory.

PAGE 47

34 0 2 4 6 0 0.2 0.4 0.6 0 2 4 6 0 0.2 0.4 0.6 Coupling Strength (C) 0 2 4 6 0 0.2 0.4 0.6 0 2 4 6 0 0.2 0.4 0.6 N( Roess | Lor ) N( Lor | Roess) SNR = 30dB SNR = 20dB SNR = 10dB SNR = 0dB Figure 2-2. Illustrating the dependency relationships betw een the Rosseler (driver) and the Lorenz system (response) as a function of coupling strength. Figure 2-3. Original EEG signals X and Y. for the EEG simulation example Figure 2-4. The nonlinearly mixed (synt hetic) EEG signals V, Z, and W.

PAGE 48

35 czw Cvz cvw cxv cxz X V Z W Y R Figure 2-5. Schematic diagram of the coupling function between the signals. 2.4.2 Intracranial EEG Simulation In this example [ 30 ], intracranial EEG signals recorded from epileptic patients are considered. Clinically, it is useful to know or find out the direction of information flow through intracranial EEG signals. This analysis may help locate the epileptic foci of the seizures as well as providing means of predicting them. In this signal, since the intracranial EEG measurements are generated by a closed system (the brain), we assu me that the dynamical statistic s of the signals observed at different channels can be modeled using a single SOM, unlike the previous synthetic example, where a separate SOM was used fo r each signal. The SOM, however, must be trained using data that represents all po ssible psycho-physiologica l states that the intracranial EEG signals might exhibit. In the case of an epilepsy pa tient, these include pre-ictal, itcal and post-ictal states, in addition to the inter-ictal state. In this example, intracrani al EEG signals collected from two different patients at different locations (labeled X and Y) are used. A synthetic nonlinear functional relationship with infl uence direction from X to V, W, and Z is created according to (2.6). Care is taken in choosing these functions to ma ke sure that the synt hetic intracranial EEG

PAGE 49

36 mixtures in (2.6) exhibit the characteristics of real EEG signals. Th e signals are shown in fig. 2-3. This is achieved by ve rifying that the time structure -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6 EEG SOM x(n)x(n+) Figure 2-6. The phase-space trajectory of the training EEG signal (projected in two dimensions) and the weights of the trained EEG-SOM (circles). and the power spectra of these signals are consistent with that of an EEG signal. ) ( ) 2 ( ) ( 05 0 ) ( ) 2 ( ) 2 ( ) ( ) 3 ( ) 2 ( ) ( ) ( n v c n z c n r n w n x c n v c n z n x n x c n y n vvw zw xz vz xv (2.6) Here, x(n) and y(n) denote the original time sequences and v(n), w(n), and z(n) denote the synthetic signals driven by the tw o original signals. In addition, r(n) is a zero-mean unitvariance Gaussian noise term. The synthetic in tracranial EEG signals are shown in Fig. 24 and the flow diagram representing the relations hips in (2.6) is depicted in Fig. 2-5. A 10x10 rectangular SOM (referred to as EEG-SOM) is trained using 3000 samples of embedded intracranial EEG data (with an embedding dimension of 10 and embedding delay of 30ms). The phase-space tr ajectory of the training data and the weights of the trained SOM are shown in Fig. 2-6. The normal EEG state is represented by the smaller amplitude activity (the dominant portion of th e training data), whereas the

PAGE 50

37 larger amplitudes correspond to the spiky, sh arp, and slow wave activity formed during the ictal state of the brain, or to artifacts formed due to muscle movements, etc. After training the SOM, the normalized similarity index between the original signals X and Y is evaluated to verify whether these intracrania l EEG signals are indeed independent or not. The dependencies between X, V, W, and Z are also evaluated using the SOM to calculate the normalized similarity index. The results are summarized in Table 2-1, where both the coupling strength and the estimated si milarity index between pairs of signals are presented. Table 2-1. Coupling strength between pairs of signals, the normalized similarity index and the original similarity index between them. Coupling Strength ) | ( ) | ( ) ( X Y S Y X S Y X S cxv = 1 -0.1668 V X -0.1112 V X cxz = 2 -0.0756 Z X -0.0612 Z X cvz = -0.8 0.0901 V Z 0.0772 V Z cvw = 1 -0.213 W V -0.336 W V czw = 0.3 -0.1225 W Z -0.1044 W Z The results obtained from the SOM-based SI measure and the original SI measure (in Table 2-1) is in perfect agreement. We conclude that X influences V and Z, V influences Z and W, and W influences Z. Comp aring these with the flow diagram in Fig. 2-5, it is seen that all directional couplings are consistent with the true construction except for the relationship between V and Z. Po ssibly, this discrepancy is due to some cancellations between the couplings from X a nd from V. Also, we can see that V is exclusively constructed from the X and the Y components and does not have any

PAGE 51

38 independent oscillations of its own, unlike W. These results indicate that the similarity index approach might not produce results that are consistent with what one would expect from the equations (if these are known) wh en the coupling diagram has closed loops. 2.5 Epileptic Intracranial EEG Data Description Intracranial EEG signals were record ed from hippocampus, sub-temporal and frontal-cortex structures of epileptic patients having a hist ory of complex-partial, partialsecondary and sub-clinical se izures of temporal-lobe fo cus, using bilaterally and surgically implanted electrodes (Fig 2-7) Using amplifiers with an input range ofmv 6 0 the recorded signals were converted to narrow-band using an anti-aliasing filter with a cut-off range between 0.1Hz and 70Hz. Using an analog-to-digital converter with 10-bit quantization precision, the narrowband signals were sampled/digitized at 200 samples/sec. Measurements involved recording intracranial EEGs from multiple sensors (28 to 32, with common reference channels ) and the recordings spanned over 6 continuous days. A total of 55 seizures, of te mporal lobe onset, were recorded from 5 patients, in the range of 6 to 18 seizures for each patient. 2.6 Testing the application of SOM-SI on Epileptic Intracranial EEG In this section, we demonstrate the utility of SOM-SI in epileptic intracranial EEG analysis and statistically ve rify accuracy of the results with the original SI. A 25x25 EEG-SOM is trained using 3000 input vector s constructed by embedding (dimension, m = 10, delay, = 30ms) intracranial EEG signals colle cted from various regions such as temporal, sub-temporal, and orbit frontal, of an epilepsy patient. The EEG-SOM needs to represent all possible EEG-dynamics, so the tr aining data must include samples from the inter-ictal, ictal, and the pre-ictal states of the patient. Fig. 2-8 shows the phase-space trajectory of the data and the PEs of the EEG-SOM in two-dimensions. The normal

PAGE 52

39 EEG state is represented by the smaller amplit ude activity (the domi nant portion of the training data), whereas the larger amplitudes correspond to the spiky, sharp and slow wave activity, Figure 2-7. Diagram of the depth and subdural electrode montage in an epileptic brain. Electrode strips are placed over the left orbitofrontal (LOF), right orbitofrontal (ROF), left subtemporal (LST), righ t subtemporal cortex (RST). Depth electrodes are placed on the left tempor al depth (LTD) and right temporal depth (RTD), to record hippocampus EEG activity. mostly formed during the ictal state. We note th at the distribution of the neurons is sparse in the higher amplitude region because of the density matching property of the SOM. To ensure generalization of the SOM, a test set of intracranial EEG signals were quantized by the trained SOM. As seen fr om Fig. 2-9, the quantization results successfully approximate the dynamics of the te st data set (projected in one-dimension). The correlation coefficient between the two si gnals was found to be 90.1%. For the most

PAGE 53

40 part, the correlation coefficient was betw een 80% and 95%. Note that the amplitude errors are higher in the larger amplitude re gions corresponding to sp ike and slow waves. This is expected because of the sparse distribution of the neurons in the higher amplitude regions. These errors can be compensated by using a larger SOM grid (> 25x25), but since the dynamics of the data are more important for the neighborhood information in the SI measure and computational complexity will be an issue, we chose not to increase the SOM grid size. 2.7 Statistical comparison between SI and SOM-SI Next, we quantify the accuracy of the SOM-SI measure relative to the original SI measure by comparing their results. SI values were calculated on nearly 39 minutes of data corresponding to a pair of signals obtained from the right temporal (RTD4) and the right sub temporal depth (RST 1) electrodes, corresponding to patient P093. The entire -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6 X1X2 Figure 2-8. The phase-space trajectory of th e training intracranial EEG signal (solid lines) superimposed on the weights ( dots) of the trained 25x25 EEG-SOM grid. interval of 39 minutes data was segmen ted into 230 non-overlapping windows of 10 seconds each. Fig.2-10 shows the interdependenc y values of both the measures. It is easy to see that the results from both the measures are in agreement to a large extent. There are also subtle differences, which need to be quantified using statistical tests.

PAGE 54

41 200 300 400 500 600 700 800 900 1000 -4 -3 -2 -1 0 1 2 3 4 time (samples) Figure 2-9. A qualitative illustration of the accuracy of the 25x25 EEG-SOM. Sample test signal (solid line) overlapped on the SOM-reconstructed output (dash and dot line). In this case, the co rrelation coefficient = 90.1%. The comparison will be two-fold: (i) identi fy if the number of windows in which the predicted directions of influence differ is si gnificant or not, (ii) given time instances where both measures agree on the direction, check if significant differences exist in predicted strengths of influence. Assuming that SOM-SI and SI measure values come from normal populations, we use the two-sided pa ired t-test to invest igate the extent of disagreement between the two methods. The te st was performed at a significance level of =0.05, over a size of 138 randomly selected samples out of the 230 available samples. Null hypothesis: H0:0) ( SI SI somd Alternate hypothesis: H1:0) ( SI SI somd Paired t-test is chosen, because the obse rvation in window 1 of original SI is obtained under similar conditions as the window 1 of SOM-based SI, and hence, the data may be said to occur in pairs. In this case, texp was found as -0.9441, whereas tcrit=t(0.05), 137=1.960. Since texp
PAGE 55

42 -30 -25 -20 -15 -10 -5 0 5 10 -0.2 0 0.2 0.4 0.6 Original Similarity Index Time ( minutes ) -----------------> -30 -25 -20 -15 -10 -5 0 5 10 -0.2 0 0.2 0.4 0.6 SOM-based Similarity Index Time ( minutes ) -------------->Interdependence N(X|Y) N(Y|X) Interdependence Figure 2-10. SOM-SI results showing the de pendencies between two signals. X and Y correspond to channels RTD4 and RST1, respectively. The time instant ‘0’ corresponds to the seizure onset (top). Results produced by the original SI (bottom). 2.8 Testing the Robustness of SOM-SI on Multiple SOMs The previous simulations successfully demonstrated the accuracy of the SOMbased measure by statistically comparing its re sults with results from the original SI measure. For application on seizures esp ecially, a 25 x 25 SOM grid was trained to embed all the dynamical states of the EEG attractor. SOM being one of the most important elements of this improvised measure, one of the pre-requisites of this approach is to ensure the following: a) that for data modeling purposes, the training set captures the essence of variance found in the dynamics of the ictal states from all the channels, for a given patient and b) the similarity indi ces computed using the SOM’s processing elements (PEs) are independent of the SOM and the corresponding tr aining data set. In other words, pair-wise similarity indices computed on two separate SOMs should be significantly close to each other, if not equal.

PAGE 56

43 We performed a simple test to check the latter point before proceeding with extensive data analysis. From the multivariate intracranial EEG data samples of the patient P093, two separate training sets we re constructed. One of the training sets consisted of portions of data sampled from the inter-ictal, ictal, pre and post-ictal states of seizures 1 & 2. The other training set consis ted of data portions picked around seizure 4 & 5. Using the same normalization procedures on both the sets and with the same set of training parameters as before, two separa te SOMs (called as SOM-1 and SOM-2 for convenience) were trained. Following that, the SOM-similarity indices were obtained from pair-wise analysis of interdependen ce among channels chosen from the ROF and LOF regions of the brain, as illustrated in Fig 2-11. The test data from the ROF and the LOF regions were picked from intervals surrounding seizures 6 & 7, seizures 4 & 5 a nd seizure 11, respectively. Fig 2-12 shows the exact time-intervals of the test data. The similarity index profiles {N1 (X|Y)t and {N2 (X|Y)t obtained from computing the SOM-SI on la rge intervals (say time t = 1,…,T) of seizure data are quantitatively compared us ing the classical correlation coefficient and error-percentage as the comparison metrics. Th e error-percentage is computed as follows T t t t tY X N Y X N Y X N e1 1 2 1) | ( ) | ( ) | ( 100 (2.7) where N(X|Y) is the normalized interdepende ncy of X on Y. Note that the notations X are Y are used to denote the two cha nnels of interest. Normalized error e quantifies the percentage difference between the interdep endency values from SOM-2 and SOM-1, keeping interdependency value from SOM-1 as the reference. From the error population, the fraction of the absolute er ror values less than 20% and the fraction less than 10% are

PAGE 57

44 computed to determine the degree of dependence of the SOM-SI measure on the data used to train a SOM. Figure 2-11. Experiment setup to compare SOM-Similarity Indices obtained from 2 separate maps. Figure 2-12. Time intervals from all the seizur es used as test data (not drawn to scale) For illustration, the results from analyzi ng the interdependency of LOF3 on LOF4 on various seizures are shown in fig 2-12. The histograms correspond to the error ensembles obtained from analyzing over l ong seizure intervals. Qualitatively, the superimposed traces in fig 2-13 indicate the extent of agreement or disagreement between the SOM-SI profiles. Table 2-2 compiles a summary of the agreement between the SOM-

PAGE 58

45 SI profiles for about 13 hours of intracranial EEG data from 5 seizures. A large fraction of errors less than 20%, supported by a high correlation coefficient between the two SOM-SI profiles suggests that there was very little disparity between the SOM-SI profiles from SOM-1 and SOM-2. Besides, the high pe rcentages also suggest the intracranial EEG data dynamics do not vary drastically from one seizure to another, and therefore the two SOM models produced the almost identical SI results. This finding will consequently support (or reinforce) our orig inal belief that a well-trained SOM and a well picked training data set is sufficient to carry out i ndependency analysis on all the seizures of a patient. Overall, pair-wise analyses of the interdependency among 6 channels (6 2C = 15 combinations) on 5 seizures of P093 was pe rformed on SOM-1 and SOM-2. The average correlation coefficient and the error results between the SOM-SI pr ofiles are shown in Table 2-3. Results from table 2-3 indicate that around 80% of the times, the differences between the SOM-SI results are less than 20%. This is not surprising considering that the differences are measured in percentages (2.7 ) and therefore even small discrepancies in the case of small dependency values can appear magnified. In addition, we also speculate that the discrepancies could be the outcome of the two SOMs being trained in an identical fashion instead of being fine-t uned to obtain the lowest rec onstruction error in each. In

PAGE 59

46 Table 2-2. Quantitative comparisons between the SOM-SI profiles obtained from SOM-1 and SOM-2. LOF3 & LOF4 data was projected on each of the SOMs and then the SOM-SI measure was applied to analyze the dependency of LOF3 on LOF4. P093, Interdependency N(LOF3| LOF4) Correlation Coefficient (%) Fraction of error less than 20% Fraction of error less than 10% Seizure 6 & 7 95.74 0.8504 0.5597 Seizure 4 & 5 98.45 0.9234 0.7543 Seizure 11 91.59 0.6452 0.3614 general, if the SOMs can be designed to obtain the lowest reconstruction error, by iteratively choosing the best se ts of parameters, a slight improvement in the performances can be easily achieved. But as it stands, a s light discrepancy can ne vertheless be always expected although it may have very little imp act in the overall scheme of analysis. 2.9 Summary The similarity index measure determines directional dependencies between two signals using the basic assumption that two related signals will have similar recurrences of the embedded state vector. This method has high computational complexity in terms of the number of samples, since a search for nearest neighbors must be performed in the phase-space of the signal. In th is dissertation, we proposed a SOM-based approach as an improvement over the original SI measure, for detecting functiona l dependencies among multivariate structures. This approach reduces the computational complexity drastically by exploiting the accurate quantization prope rties of the SOM in representing the dynamics of the signal in the phase sp ace. Another advantage of the SOM-based

PAGE 60

47 approach is that the difficulties that the orig inal similarity index approach encounters in handling nonstationary data (such as the nece ssity to tweak paramete rs) are eliminated by training the SOM using samples from various regimes of the nonstationary system. On 0 50 100 150 200 250 0 0.1 0.2 0.3 0.4 time (minutes)SOM-Similarity Indices)P093, Seizure 67 -100 -50 0 50 100 0 10 20 30 40 50 60 70 error (%)Patient P093, Seizure 6 & 7 0 50 100 150 200 250 300 -0.2 0 0.2 0.4 0.6 0.8 time (minutes)SOM-Similarity Indices Patient P093, Seizure 4 & 5 -50 0 50 0 5 10 15 20 25 30 35 40 Error (%) 0 50 100 150 200 250 -0.1 0 0.1 0.2 0.3 0.4 time (minutes)SOM-Similarity Indices Patient P093, Seizure 11 -100 -50 0 50 100 0 20 40 60 80 100 Error ( % ) Figure 2-13. Comparing interdependencies between channels LOF3 and LOF4. Left : SOM-similarity profiles from th e output of SOM-1 and SOM-2 are superimposed. Right: Histogram of the errors in %. Top: Seizure 4 & 5 Middle: Seizure 6 & 7. Bottom: Seizure 11.

PAGE 61

48 Table 2-3. Summary of the comparisons between the SOM-SI pr ofiles from SOM-1 and SOM-2. Each row represents the statis tics (mean and varian ce) of pair-wise SOM-SI analyses of the epileptic intr acranial EEG data from 6 channels (15 combinations). P093 Correlation Coefficient (%) Fraction of error less than 20% Fraction of error less than 10% Seizure 6 & 7 94.32 2.85 0.79 0.1 0.54 0.12 Seizure 4 & 5 97.46 1.08 0.91 0.06 0.73 0.12 Seizure 11 93.24 2.06 0.71 0.08 0.41 0.07 the other hand, the SOM-based approach might suffer from inaccuracy if the quantization is severe. Therefore, the size of the SOM could be decided by a trade-off between representation accuracy and co mputational complexity. Through simulations, we also demonstrated the noise-robustne ss features of the measure.

PAGE 62

49 CHAPTER 3 SPATIO-TEMPORAL CLUSTERING MODEL 3.1 Introduction In Chapter 2, we proposed the SOM-SI measure [ 29-31 ] and demonstrated SOM’s resourcefulness as a model infrastructure fo r computational purposes. Often, the timesequences in a multi-variate system share sim ilar information that is reflected in their interactive or synchronization abilities. By definition the word similar could mean that the information shared among a set of channe ls is stronger than the information they share with other channels. Such spatial sim ilarities could possibly be momentary up to a few seconds or could even stre tch to several minutes or hours As we postulated earlier, spatio-temporal similarity changes could be one of the driving factor s to trigger certain events in biological systems. From a seizure point of view, we suspect that analyzing the temporal changes in channel similarities could reveal some interesting aspects about the epileptic brain. Similarity-based time-series clustering [ 25, 44 ] is a well researched topic in the area of dynamical graph theory. It is an extr emely useful approach to characterize the spatial groupings in time-sequences. Simila r time sequences are typically grouped based on their mutual interactions. In this study, us ing the SOM-SI as a computational tool to derive the distance/similar ity/proximity matrix, we pr opose a clustering model to dynamically analyze the spatio-temporal gr oupings in muti-variate time sequences. Simple toy simulations are presented to vali date the robustness of the model. Since our

PAGE 63

50 goal is specific to intracranial EEG, we use this model in chapter 4 to track the seizurerelated temporal changes in groupings. 3.2 Model for Spatio-temporal clustering In this section, we propose a clustering approach to ex tract information on spatiotemporal distribution of multivariate time-meas urements. A 3-fold approach consisting of spatial-discretization of the data using spectral-clustering technique [ 57, 64, 100 ], temporal quantification using hamming dist ance, followed by application of another clustering technique, is presen ted in Fig 3-1. The rational will become apparent during the explanation. Figure 3-1. Block diagram to extract sp atio-temporal groupings information in Multivariate EEG structures Spectral clustering is one of the many clustering methods that use subspace decomposition on data-derived affinity matrix to achieve data-clustering. Using kernel methods, the data samples are projected onto a higher dimensional space where the discriminant analysis is much easier. Projec ting the data onto a feature space results in tightly formed clusters such that the between cluster en tropy is maximized and within cluster entropy is minimized. Spectral Clustering is inspired by the normalized cut theory in computer science where the distance between the nodes/data samples is interpreted as

PAGE 64

51 an affinity matrix on which subspace deco mposition yield membership labels of the nodes/data samples. In our study, we use the standard spectral clustering algorithm by Ng et al. [ 64 ] to spatially cluster the similarity-indices obtai ned by the SOM-SI technique. The algorithm is outlined below, 1. Assuming N to be the number of nodes S = {s1, s2, …. sN} and k as a priori specified cluster size, form a similarity matrix A RNxN using a distance metric. 2. Define D to be the diagonal matrix whose (i,i)th element is the sum of A’s ith row and construct the matrix L = D-1/2A D-1/2. 3. Find x1, x2, ….xk the k largest eigen-vectors of L and form the matrix E by stacking [x1, x2, ….xk] RNxK the eigen-vectors in columns. 4. Form the matrix Y from X by renormalizing each of X’s rows to have unit length. 5. Treating each row in Y as a point in Rk, cluster them into k clusters using k-means or any other distortion minimization algorithm. 6. Finally, assign original node (si) to cluster j if and only if row i of the matrix Y was assigned to cluster j. Pair-wise evaluation of SOM-SI measur e on all the possible combinations ((CN 2, where N is assumed to be the number of channels) of a portion of a multi-variate time series leads to k = 2*(CN 2) similarity indices in the bounds of [0, 1]. k is multiplied by 2 because of the asymmetric nature of the SOMSI measure. If we imagine the time-series as various inter-connected nodes in a multi-di mensional graph, the SOM-SI similarity indices represent the affinity or rather, th e weights of the connection, between those nodes. Therefore, we can translate th em into a square matrix of size N x N, where N is the # of channels. Since the wei ghting is normalized between 0 and 1, the diagonal elements representing the affinity of a cha nnel with itself, are coded as 1.

PAGE 65

52 However, to be able to perform spectraldecomposition on an affinity matrix, Ng’s algorithm [ 64 ] requires that the affinity matrix be s quare and symmetric in nature. This is because the eigen decomposition yields orthogonal column vectors (also called eigenvectors) only if the projection matrix is square-symmetric. The asymmetric matrix can be transformed to a symmetric matrix by adding it to its tran spose and dividing each entry by 2. Mathematically, if 1 represents the asymmetric affinity matrix, the transformation can be represented as, 2 ) (1 1Tk (3.1) The transformed affinity matrix now repr esents the average information exchanged between all pairs of channels. This implies th at we do not have the luxury of using the asymmetric nature of the dependencies to create a membership grouping among channels. In epileptic study however, we will see that the averaging of the dependency information will not affect channel dependencies because of the earlier observation that there is no major difference in the driving and re ceiving information of the channels. On the transformed affinity matrix eigen decomposition can be done, followed by K-means clustering. Consequently, we have a set of labeled clusters (number of clusters is set, based on significant eigen va lues) representing the membership of the channels. If the above procedure is repeated over consecutive time windows (overlapped or non-overlapped), channel groupings obtained on each time-window can be arranged as in a matrix form as in (3.2).

PAGE 66

53 2 1 . . 2 1 3 . . . . . . . . . . 2 3 . . 2 2 1 1 3 . . 2 2 3spect (3.2) To characterize the average clustering of the channels over a l onger period of time, we propose another, albeit simp le, hierarchical clustering ap proach that uses hamming distance to derive the proximity matrix. 3.3 Temporal Quantification Using Hamming Distance We showed in the previous section that the multivariate time-series can be grouped by using similarity-based cluste ring techniques such as spectr al clustering. The spectrally clustered labels specify the groups of channe ls exhibiting high degr ee of within-cluster similarities and low between-cluster similarities Often in applications such as epileptic EEG analyses, it is important to identify ch annel groupings over a l onger period of time. The epileptic seizures, in particular, are charac terized by dynamic states (inter-ictal, ictal, pre-ictal & post-ictal) that are known to possess both loca l and global spatio-temporal clusters. Channels associate and de-associat e in time; however, depending on the psycho physiological state of the brain, certain groups of channels might have a higher likelihood of sharing same channel labels thus forg ing a long-term association. In epileptic intracranial EEG, identifying such state-depende nt clusters may provide us with useful insights on the evolution of brain patterns during seizure states. State dependent connections can be qua ntified by clustering rows of the spect matrix that are similar with each othe r over a longer time interval, say T. Mutual Information is one of the widely used metrics to quantitatively evaluate the similarity between two discrete processes. However, in our study, using mutual information to

PAGE 67

54 assess the degree of similarity between two channels/rows can have serious drawbacks. For any two channels, two pairs at different time-instances, having same elements (say 3 & 1) need not bear any semblance. This is because the cluster la bels are arbitrarily labeled at each time-instance and do not ha ve a fixed representation across time. Technically, it implies that the probability termj i pv j u i ), (, where i, j are the channels/rows in the spect matrix and u, v are the channel labels/memberships, can be skewed unless u = v. In this context, a very simple statistic to determine the degree of similarity would be to compute the relative frequency of any two channels sharing the same labels/groupings. In other words, in a time window of length T, we check the average number of times when the two channels of interest, share the same cluster label. In an algebraic context, the above operat ion is equivalent to computing pair-wise hamming distance in a time window T. Of course, similarity would be obtained by subtracting the hamming di stance from 1. That is, If ham ijd is the hamming-distance between channels ‘i’ & ‘j’, similarity in probabilistic terms is obtained as, ham ij sim ijd p 1 (3.3) Thus, computing the pair-wise similarity for all i & j combinations will result in a P matrix of size N x N (N is the number of channels). For convenience, we will call the matrix P the cluster-similarity matrix in all our future references. Hierarchical clustering on the cluster-similarity matrix P will yield information on the cluster groupings over a time T. In the context of intracranial EEG data, clustering

PAGE 68

55 will thus enable us to know the groups of ch annels that have similar behavioral structure in the brain, over a wider time interval. 3.4 Simulations with Synthetic Data This section presents simple simulations to demonstrate the validity of hierarchical clustering approach to cluster a multivariate time-series data based on their mutual interactions. 3.4.1 Roessler – Lorenz System Recall the Rosseler and the Lorenz system example from section 2.3.1. The idea is to first generate non linearly coupled dynami cal signals in each system and identify mutual interactions using the SOM-SI measur e. A hierarchical clustering algorithm such as spectral clustering is then applied to group the time-series into pre-determined number of clusters. For the zero-coupling case (C = 0), 5 000 data-samples from the non-transient portion of each of the 6 components are extrac ted as shown in Fig 32a. Observe that the time profiles of the first and the second co mponents of each system have a strong phasesimilarity between them. The dynamics of each system evolve as a result of non linear interactions within their components. Ideally, only two SOMs suffice, each trai ned to represent the Rosseler and the Lorenz system individually. However, if non linear coupling is in troduced between the two systems, the attractor dynamics change depending on the coupling strength. This situation will force us to ha ve individually SOMs to re present each dimension of a system. Therefore, without any loss in gene rality, we train individual SOMs for each component of the Roessler and the Lorenz sy stem. Specific to this simulation, each SOM was a 2-d grid of size 8 x 8 processing elements.

PAGE 69

56 Affinity between the components is obt ained by pair-wise computing the SOMsimilarity indices. Further, groupings among th e 6 time-profiles are obtained by applying spectral clustering algorithm. Spect ral clustering, for a cluster-size n = 2, clear distinguishes between Rosseler and the Lore nz components. For visual purposes, we present the dendrogram (as in Fig 3-3a) ge nerated in MATLAB 7.0. As expected, the Rosseler and the Lorenz components form separate individual groupings. The fusion levels also support the fact th at the coupling strength C is ve ry low (in this case = 0), by clearly suggesting weak betw een-component interaction. The next case consists of establishi ng a coupling (C = 5) from the X2, the 2nd component of the Roessler to Y2, the 2nd component of the Lorenz system. Fig 3-2b shows the time-traces of the Roessler and th e Lorenz components. It is clear that Y2, (and therefore the Y1, the1st component as well) has an altered time-profile. In addition, since Y2 and Y3 of the Lorenz system also interact in a non linear fashion, we observe a clear change in the time-dynamics of Y3. The procedure of applying SOM-SI to compute the affinity matrix followed by spectral clustering is repeated. Dendrogram is shown in Fig. 3-3b, to graphically illustrate the cluster or ganization. Firstly, we can observe that the topology of this dendrogram (for C = 5) is drastically different from the topology of dendrogram for C = 0. We no longer see tw o clearly separated Rosseler and Lorenz groupings. The fusion levels among the Lorenz components is still pretty much at the same level as they were for C = 0. However, the fusion levels between the combined Lorenz components and X1, X2 combined is reduced cons iderably, reflecting strong coupling introduced earlier. The topology-diffe rence between the two dendrograms in fig.

PAGE 70

57 3-3b suggests that due to the non linear coupl ing introduced in (2.5 ) through the coupling constant C, the dynamics of X2 & X1 are more closer to the dynamics to Y than to X3. (a) (b) Figure 3-2. Time profiles of the Roessler and the Lorenz components. a) Coupling strength, C = 0 and b) C = 5

PAGE 71

58 Table 3-1 shows the cluster group ings, obtained as a function of n for various coupling constants. While it is hard to obtain any distinction among the groupings, with n = 4, the difference is clear for n = 2. Spectral Clustering, consistent with the model construction and the observati on in dendrogram, groups the X1, X2, Y1, Y2 & Y3, components separately from X3. Table 3-1. Cluster groupings for the simulation example. C = 0 C = 5 # of clusters, n = 2 (X1, X2, X3), (Y1, Y2 ,Y3) (X1, X2, Y1, Y2 ,Y3), (X3) # of clusters, n = 3 (X1, X2, X3), (Y1, Y2), (Y3) (X1, X2), (X3), (Y1, Y2, Y3) # of clusters, n = 4 (X1, X2), (X3), (Y1, Y2) (Y3) (X1, X2), (X3), (Y1, Y2) (Y3) 3.4.2 Linear Model Consider a set of colored noise data obtained as the band-pass filtered output of an additive white-Gaussian stochastic process of zero mean and unit-variance (as shown in Fig. 3-2). Each of the outputs can be denoted by si(t), i = 1,2,…..,8 & t = 1,2,…,T. We construe xi(t) as the observation measurements, sample d from eight (8) different channels of a spatio-temporal system, over an interval T. Observe that xi(t) is constructed by adding Gaussian white noise (of variance between 0.1 and 3) to si(t). ) ( ) ( ) ( t n t s t xi i i (3.4)

PAGE 72

59 Design of noise ni(t) is such that a variance higher than 3 will destroy the structure of xi(t) to the extent that the signals become completely un-correlated with each other. Even though the two band-pass filters denoted by BPF1 and BPF2 respectively (as shown in fig. 3-3 and fig. 3-4) share different pass-bands, th ere is a significant amount of overlap as seen by the narrow separation of th eir pass-bands. It should also be noted from fig. 3-4 that some of the output measurements xi(t) can share the same band-pass filters. However the additive noise component ni(t) can introduce some stochastic differences in their linear structures. Measurements xi(t) that share the same band-pass filter can loosely be considered analogous to time-structures having reasonably str ong linear-interactions. This is because information exchanged be tween two time-series can sometimes cause their phase dynamics to be similar. We compute pair-wise cross-correlation indices (1) among all the channels to quantify the interactions or rather, the exact amount of information exchanged among the spatio-temporal measurements xi(t). Due to the fact that the signals are stochastic, linear and do not possess any complex or chaotic phase dynamics, we felt that the crosscorrelation computation would suffice to indicate the amount of their pair-wise interactions. The cross-correlation indices as denoted by 1 8 8 1 2 8 1 1 8 1 8 2 1 2 2 1 1 2 1 8 1 1 2 1 1 1 1 1. . . . . . . . R R R R R R R R R R denote the affinity matrix for the time seri es corresponding to say window 1. Repeating this process over a number of windows W will result in the spatio-temporal affinity

PAGE 73

60 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Lor_1 Lor_2 Lor_3 Ros_1 Ros_2 Ros_3 Fusion Level Clustering on Roessler Lorenz components ; Coupling Strength, C = 0 (a) 0.15 0.2 0.25 0.3 0.35 0.4 0.45 Ros_1 Ros_2 Lor_1 Lor_2 Lor_3 Ros_3 Fusion Level Clustering on Roessler Lorenz components ; Coupling Strength, C = 5 (b) Figure 3-3. Dendrogram illustra tion of similarity-based time series clustering. a) Coupling strength, C = 0. The componen ts of the two systems are clearly clustered separately. b) For C = 5, the 1st and 2nd components of the Roessler system form a stronger grouping with th e components of the Lorenz system than the 3rd component of the Roessler system itself.

PAGE 74

61 matrix WR R R R. .2 1, where W is the total number of windows. For each window, the spatial configuration/arrangement of the band pass filters is changed as follows Configuration 1 (c1): Channels 1, 2, 3 & 4 are passed through BPF1 and Channels 5, 6, 7 & 8 are passed through BPF2. Configuration 2 (c2): Channels 1 & 3 are passed through BPF1 and Channels 2, 4, 5, 6, 7 & 8 are passed through BPF2. Configuration 3 (c3): Channels 1, 2, 3, 4, 6 & 7 are passed through BPF1 and Channels 5 & 8 are passed through BPF2. Configuration 4 (c4): Channels 1, 3, 6 & 7 are passed through BPF1 and Channels 2, 4, 5 & 8 are passed through BPF2. Total number of Windows W is fixed to 100 and the probabi lity of random occurrence of channel configurations is as follows: {p(c1) = 0.5, p(c2) = 0.1, p(c3) = 0.1, p(c4) = 0.3}. Applying the spectral clustering algorithm on the spatial-correlation matrix of each window Ri resulted in ) 100 8 (. 2 . . 1 2 2 1 2 1 2 1 . . 1 2 1 2 2 2 2 1 . . 1 2 1 2 2 2 2 2 . . 1 2 2 1 2 1 2 1 . . 2 2 2 1 2 2 1 1 . . 2 1 1 2 1 2 1 1 . . 2 2 2 1 2 2 1 1 . . 2 1 1 2 1 2 1x spect c1 c3 c2 c4 c4 c2 c1 . . c3

PAGE 75

62 This illustration shows how the cluster labels are assigned to the channels depending on their configuration, denoted by ci’s. Heuristically, it is evident that the channels can either be clustered in groups of four or in groups of two. For the case when the number of clusters equals four (4), the hammingdistance based clustering resulted in cluster assignment of the channels C1: {1, 3}, C2: {2, 4}, C3: {6, 7}, C4: {5, 8}. Specifying only 2 clusters resulted in C1: {1, 2, 3, 4}, C2: {5, 6, 7, 8}. Figure 3-4. Model showing the generati on of linearly dependent signals using configuration 1. We can easily see that the results with 4 clusters and 2 clusters are in perfect agreement with the construction of the spat io-temporal configurations. Channel 1and 3 are always together and belong to the same cluster, regard less of the configuration. Likewise, channels 5 and 8 always belong to the same cluster; alt hough different from the cluster to which 1 and 3 belong to. This is because the source data Z(t) is always passed

PAGE 76

63 through BPF2 to obtain channels 5 and 8, unlike ch annels 1 and 3 that are obtained by filtering Z(t) by BPF1. The channels 2 and 4 and similarly 6 and 7 change their memberships often enough and therefore can be partitioned into separate clusters. 0 200 400 600 800 1000 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Frequency (Hz) BPF1 BPF2 Figure 3-5. Plot showing the two band pa ss filters, separated by a narrow pass band -1 -0.5 0 0.5 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Real PartImaginary PartPole Zero configuration for BPF1 Z1 = 0.9 P2 = 0.5 P1 = -0.75 -1 -0.5 0 0.5 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Real PartImaginary PartPole Zero configuration for BPF2 Z1 = 0.8 P2 = 0.65 P1 = 0.7 Figure 3-6. Z-planes showing the polezero configurations for th e two band-pass filters in Fig 3-5.

PAGE 77

64 3.5 Summary Adopting a similarity-based time-series cl ustering approach, we proposed a spatiotemporal model to evaluate long term cluste rings in multi-variate time series data. The earlier proposed SOM-SI measure is used as the distance measure to evaluate the affinity between various nodes in time-series. One of the major issues in any clustering problem is determining the appropriate number of clusters. In our model we use the spectral clustering algorithm in which it relatively easie r to specify the cluster size. This is because the significant eigen values can dir ectly provide a reasona ble indication of the number of centroids in the data. Similarly, through simulations we showed that the dendrograms can also be used to get useful hints for choosing the number of clusters. Simulations on Roessler-Lorenz dyamical syst em and on a linear set up demonstrate the efficacy of our proposed approach to evalua te long term associability patterns among dynamically connected nodes.

PAGE 78

65 CHAPTER 4 APPLICATION OF SPATIO-TEMPORAL MO DEL TO EPILEPTI C INTRACRANIAL EEG 4.1 Introduction In the last chapter, we proposed a spatio -temporal model to extract groupings from long term multivariate recordings. This chap ter will focus on the application of that model on the epileptic intracranial EEG time series. The first part of the chapter will describe the details on the application of th e model and the second part will discuss the results of analyses on 11 complex-partia l seizures, from 2 epileptic patients. 4.2 Application of SOM-SI on Epileptic Intracranial EEG Data The temporal changes in the spatial struct ure of an epileptic brain was analyzed on twenty four (24) representative channels re corded bilaterally from the orbito-frontal, temporal and sub-temporal regions on the br ain. One of the fundamental requirements for analyzing the dynamics of a non linear system is to construct the state-space attractor from just a single recording of the time-ser ies. From previous studies that estimated intracranial EEG attractor size using correlation-dimension techniques [ 38-39 ], we know that the EEG state space is bounded between 3 a nd 10. This is of course assuming that the data is entirely noise-free and infinite amount of measurement data is available. The large range in dimension is due to the fact that the properties of the EEG dynamical system undergo dimension changes with changes in br ain’s physiological stat es and pathological conditions. On our intracranial EEG data, the embedding dimension (m) and the delay

PAGE 79

66 ( ), however, were chosen to be m = 10 and = 4. The parameters were compatible with other studies [ 38-39 ], performed on the same data. The following steps describe the pro cedure to track the spatio-temporal connectivity patterns in intracranial EEG data. 1. Choosing suitable embedding parameters, in tracranial EEG attractors in higher dimensional state space is constructe d using Taken’s time-delay embedding theorem [ 98 ]. On 10 second epochs, pair-wise analyses of interdependence among 24 channels are computed using the SOM-SI measure. 2. The similarity indices, from every window, are translated into an asymmetric similarity/affinity/proximity matrix. Si nce spectral cluste ring involves eigendecomposition on a symmetric affinity matrix, the asymmetric matrix is transformed to a symmetric matrix by addi ng itself to its tran spose and dividing by 2. 3. With the number of clusters (say n1) speci fied apriori, spectr al clustering on the affinity matrix results in channels be ing labeled as one of the n1 clusters. 4. Steps 1 to 3 are repeated for all the successive windows, representing 10-second stationary segments. However, the overall ab ility of the channels to associate with each other over a longer time-duration needs to be quantified. On 30 minute time segments (equal to 90, 10-second windows), pair-wise Hammingdistance based cluster-similarity matrix P is computed among all the channels. The matrix elements essentially index the proba bility of channels to group into the same cluster over a 30 minute time interval. Spectral clustering or any other clustering al gorithm on the cluster-similarity matrix P will result in final cluster memberships. The number of clusters is fixed to n2. For computing similarity indices in step 1, the epoch length of 10 seconds is chosen as a tradeoff between stationarity and sample-size requirements. Also note that the successive windows are 10 seconds apart (a lternate 10 second windows) for reasons specific to computational feasibility.

PAGE 80

67 For illustration, we present the results of SOM-SI mutual interactions on a seizure data in Fig 4-1. Pair-wise SOM-SI analysis of the interdependence between 24 channels resulted in some interesting observations Firstly, the SOM-similarity indices are temporally averaged (over 3 windows) to smoothen out the fluctuations between windows. On the smoothened interdependen ce indices, for each window and for each channel, we find the maximum driving influe nce that the channel exerts on any other channel. Over windows (in time) these ma ximum driving indices give the maximum driving ability of the particular channel of interest. The maximum driving abilities are evaluated for every channel unde r consideration and for simp licity; we show the driving abilities for just 6 channels in fig 4-1. In the inter-ictal stage, low driving ability of all the channels indicate that the channels are de-synchronized, even though they exhibit an upward trend. Synchronization go es up momentarily a few minutes pre-seizure and at the onset of the seizure, there is a sudden drop followed by a gradual in crease post-seizure. Higher degree of post-seizure global synchr onization is followed by a gradual drop, leading to the inter-ictal state. It is clear from fig.4-1 that the temporal -evolution of the interdependency values exhibit distinguishable patterns at different stages of seizure [ 30 ]. However, the information on the spatial changes among the ch annels is still missing. As we stated earlier, one of the primary objectives of this study is to determine the spatialconfigurations and their temporal changes, at various stages in a seizure data. In this context therefore, clustering can be a useful tool to en able us to quantify groupings among channels and progressively track their changes.

PAGE 81

68 We now describe step 2, with more deta ils. The channel interdependencies obtained from SOM-SI, as indicated by 1 in (3.8) represent the spatio-temporal correlation indices obtained by computing pa ir wise similarity-index among 24 channels. In spectral clustering jargon, the 1k matrix can also be interpreted as an affinity matrix representing the pair-wise distances between 24 nodes. -100 -80 -60 -40 -20 2 20 40 60 80 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time(mins)-----> Maximum averaged driving ability 0 ROF1 LOF1 RST1 LST1 RTD4 LTD3 Figure 4-1. Maximum average driving ability of each of the six (6) channels, nearly 100 minutes before and 70 minutes after Se izure-1 in patient P092. (The thin vertical bar corresponds to the time when seizure occurred (0 to 2 on the xaxis). For clarity, the box inside the figure shows a small portion of the maximum average driving ability of each of the 6 channels, baseline-offset by different scales.) Drop in synchronization followed by an abrupt increase in phase synchronization at the onset of seizure is evident. Synchronization across channels during seizure is also clearly seen However, notice that the affinity between nodes will not be symmetric, as is usually the case in spectral clustering. But the transpose transformation trick, mentioned in step 2, is applied to bring it to symmetric form Since the weighting is normalized between 0 and 1, the diagonal elements representing the a ffinity of a channel w ith itself, are coded as 1. Consequently, after spectral-clusteri ng, we have a set of labeled clusters

PAGE 82

69 representing the membership of the channels [ 32 ]. Repeating this procedure on every 10sec windows will yield a discrete-valued matrix spect similar to (3.10). Typically, the choice for the number of clusters n1 in step 3 is conditioned on the significant eigen-values (say 90 %). In our analysis, the su m of first 3 eigen-values typically ranged from 60% to 80% of the tota l variance, due to changes in seizure states. Considering this variability between epochs a nd the fact that the number of clusters need to be the same for all epochs in order to be able to determine the overall grouping in channels (using cluster-similarity matrix P ), we fixed the membership size to n1 = 3. The notion of pre-seizure state has wide ly been debated upon. A lot of studies argue against the notion of a ‘p re-seizure’ state transition. However, experimental studies using non linear dynamics have shown that th e quantitative descriptors of EEG exhibit seizure pre-cursors in the form of inter-ictal to pre-ictal state transitions. The pre-ictal transition time is not exactly known, however l iterature suggests that it has a broad range of 5 minutes to 120 minutes before seizure. Therefore in step 5, as a tradeoff between state transition periods and time-resolution, we choose 30 minutes time window to characterize both the pre-ictal and the post-ictal periods. 4.3 Results Following are the clustering results from applying the spatial-clustering procedure on different seizures. Patient P093 This patient had a history of complex partial se izures, localized in the mesial structures of the temporal lobe. Surgery revealed a lesi on in the right hippocampus (RTD electrodes) region.

PAGE 83

70 The set of 24 channels are listed below, Channels 1 to 4: {LTD3, LTD5, LTD7, LTD9} Channels 5 to 8: {RTD4, RTD6, RTD8, RTD10} Channels 9 to 12: {LST1, LST2, LST3, LST4} Channels 13 to 16: {RST1, RST2, RST3, RST4} Channels 17 to 20: {LOF1, LOF2, LOF3, LOF4} Channels 21 to 24: {ROF1, ROF2, ROF3, ROF4} Before data analysis, a validation test wa s performed to check whether application of different clustering algorithms on P would consistently result in same cluster memberships or not. For a given number of clusters n2, it turned out that all the clustering algorithms including spectral clustering pr oduced the same outputs. Therefore, we decided to choose the simple hierarchical clustering algor ithm used in the Matlab 6.5 package owing to its graphical support. Cluster-similarity matrices P indicating the probability that two channels share the same grouping in a 30-minute time segment are sh own gray-scale coded in Fig. 4-2. Preseizure analysis on 30-minute windows is show n for up to 3 hours. Similarly, the postseizure analysis is shown for the first 30 minutes The ability of the left side channels to have a higher tendency to group together compar ed to the right hemisphere channels is quite noticeable from fig 4-2. In addition, th e orbito frontal lobes seem like the only area on the brain to have a high probability of making a cross-hemisphere grouping. On the left hemisphere, the LST and the LTD channels are consistently seen to share the same clusters.

PAGE 84

71 To confirm the observations from fig 4-2, the hierarchical clustering algorithm was applied on each of those P matrices. Fig. 4-3 graphically illustrates two instances of the clustering outputs through dendr ograms. A dendrogram is strict ly defined as a binary tree with a distinguished root that has all the data items at it s leaves. Conventionally, all the leaves are shown at the same level of the draw ing. The ordering of the leaves is arbitrary. The heights of the internal nodes are related to the metric information (P here) used to form the clustering. Using a threshold of 0.4 and the average-linkage technique to determine fusion levels, clustering was perf ormed on a pre-defined number of clusters ( n2). For determining apriori the number of clusters n2, several dendrograms were visually analyzed. There seemed to be at least 3 to 4 strong groupings among channels in most of the dendrograms. For consistency, th erefore, we chose to fix the number of clusters n2 to 3 for all the analyses. Both dendrograms in fig 4-3. clearly tran slate the spatial patterns observed in the corresponding P matrices of fig 4-2. The top dendr ogram in fig 4-3. corresponds to the 2.5 hours to 3 hours time window (indicated by -5) in fig 4-2. It is easy to see that the dendrogram considers the RTD and the RST as isolated clusters due to their weak between-cluster fusion level. Si nce the number of clusters n2 is restricted to 3, all the remaining channels form a single large cluste r. Similarly, the bottom dendrogram in fig 4-3. corresponds to the P matrix indicated by -1 in fig 4-2. In this case, the RST and the RTD channels group into one cluster; also well supported by a dark patch in fig 4-2. This enables the LST/LTD channels and the LOF/RO F channels to group t ogether as separate clusters.

PAGE 85

72 Figure 4-2. Seizure 11 of patient P093: Cl uster-similarity matrices P indicating the probability that two channels share the same cluster label in a 30 minute timeinterval.

PAGE 86

73 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 ROF2 ROF4 ROF3 LOF1 LOF2 ROF1 LOF3 LOF4 LTD3 LST4 LST1 LST2 LST3 LTD5 LTD7 LTD9 RTD4 RTD6 RTD8 RTD10 RST1 RST2 RST3 RST4 Fusion LevelP093, 2.5 hours before Seizure 11 0 0.1 0.2 0.3 0.4 0.5 0.6 RTD4 RTD6 RTD8 RTD10 RST1 RST2 RST3 RST4 LTD3 LST3 LST4 LST2 LTD5 LTD7 LTD9 LST1 LOF1 LOF4 LOF2 ROF1 LOF3 ROF2 ROF4 ROF3 Fusion level P093, 0 30 minutes before Seizure 11 Figure 4-3. Dendrogram repres entation of the cluster result s in Seizure 11, P093. TOP: Dendrogram corresponding to 2.5 hours before seizure. BOTTOM: Dendrogram corresponding to the 30 minute pre-seizure period.

PAGE 87

74 The overall cluster configurati on is listed in table 4-1. Table 4-1. Spatio-temporal groupings as obtained for seizure 11 of patient P093 P093, Seizure 11 C1 C2 C3 Pre seizure, (2.5 – 3 hrs) RTD RST LTD, LST, LOF, ROF Pre seizure, (2 – 2.5 hrs) RTD, RST LOF, ROF LTD, LST Pre seizure, (1.5 – 2 hrs) RTD, RST LOF, ROF LTD, LST Pre seizure, (1 – 1.5 hrs) RTD, RST LOF, ROF LTD, LST Pre seizure, (30mins – 1 hr) RTD, RST LOF, ROF LTD, LST Pre seizure, (0 30mins) RTD, RST LOF, ROF LTD, LST Post seizure, (30mins – 1hr) RTD, RST LOF, ROF LTD, LST We summarize the spatial patterns at different time-intervals of seizure 11 as follows: 1. The LST and the LTD channels in particular exhibit a strong tendency to belong to the same group. 2. The LOF and the ROF channels form a str ong bi-lateral homol ogous connection, as seen from all the matrices in fig 4-2. 3. Relatively strong similarity can be se en between RTD and the RST channels. 4. Common observation in all the matrices is the strong similarity between the left hemisphere channels as opposed to the righ t hemisphere channels. This is reflected in the ability of LOF channels to have a higher probability of sh aring clusters with other left-hemisphere channels as seen in Fig 4-2. 5. Interestingly, no temporal changes are seen in the spatial-patterns yet. In Fig 4-3. even though the two dendrogram s share certain common attributes, for instance, channel labels within a group, observe slight dispar ities in topologies and in the order of connections. Quantifying the differe nce between topologies may possibly serve to differentiate between dendrogr ams, obtained at different states of a seizure. Statistical

PAGE 88

75 techniques, such as Double Permutation Statis tics (DPS) and mantel statistics are widely [ 50, 51 ] known to be best suited to compare two dendrograms. However, the comparisons usually test the hypothesis that the two dendrograms under comparison are not more similar than the dendrograms randomly genera ted in terms of topology, leaf positions and fusion level positions. Such statistical comparisons are particularly useful when the goal is to detect whether the spatial-correlations have significantly cha nged over time or not. We address this particular issu e of tracking the temporal chan ges in spatial connections in chapter 5. In the rest of this chapter, how ever, we analyze the spatial groupings for a larger set of seizures, from different patients. 4.4 Statistical Validation The cluster configurations observed from analyzing 30 minute segments necessitates validation. Recall that, we partia lly validated our model (until the spectral clustering stage), using synthetically coupled multivariate time sequences (non linear and linear, both). Simulations involving creati on of dynamic graphs involve multidimensional time series that continuously change clus ter memberships over time. Determining the average spatio-temporal groupings from a co llection of multi-variate time-series is relatively easier to demons trate in linear coupling cases However, non linear dynamic model constructions are extremely hard and mo stly, non-trivial. We therefore decided to pursue a posterior verification of the time-a veraged cluster groupi ngs on the intracranial EEG data, using the quas i-surrogate analysis [ 40, 66-68, 73, 94-95 ], technique. Recall that the cluster groupings obtained over 30 minute time-segments involve two steps. First step consists of applying spect ral clustering technique on the SOM-similarity indices (computed on 10 second intracranial EEG data segments). Then similar grouping patterns among channels are extracted by usi ng hierarchical cluste ring approach on the

PAGE 89

76 cluster-similarity matrices P. Specifically, 91 10-second windows (spanning 30 minutes, 3 windows per minute) are analyzed in each pass to obtain long-term groupings among channels. In order to validate this 2-step approa ch, we define our hypothesis as follows H0: The average within-cluster channel inte raction at each window (out of 91 10-second windows) is not significantly different from the corresponding between-cluster channel interactions. We propose to test this hypothesis on all the 3 ( n2) clusters separately, for every 10second window within the 30 minute period. Within-cluster interaction is computed by averaging the pair-wise similarity indices for all the channels within a cluster. For between-cluster interaction, the pair-wise interactions between 3 channels, picked randomly from each of the 3 clusters are computed. A between-cluster interaction stat istic is formed by computing the average interactions from random sel ection of 3 channels (one fr om each cluster), over a number of trials. We found that this statistic follows a quasi-normal distribu tion, implying that the within-cluster interaction value can now be compared with the mean and the variance sample estimates of the between-cluster sta tistic. Mathematically, we construct the zscore as follows ) (t b b i w i tC C C Zt t t = 1, 2,.. ,90 & i = 1, 2, 3. (4.2) where t i wC is the within-cluster interaction at time ‘ t ’, for cluster ‘ i ’. tbC is the mean and ) (t bC is the standard-deviati on of the between-cluster interaction at time ‘t’. i tZ reflects the z-score and is considered signi ficant at the 95 percen tile significance if i tZ>

PAGE 90

77 1.96 (reject H0). In table 4-2 the bolded value in each cell represents the number of windows (out of 91) having signi ficant z-scores in the 30 minute period corresponding to fig 4-2 (P093, Seizure 11). It is easy to obser ve that the null-hypothe sis H0 is rejected beyond doubt, validating the clustering results. Table 4-2. P093, Seizure 11: Over each 30 mi nute (91 samples total) window, number of times the within-cluster interaction is gr eater than between-clu ster interaction, at 95% significance level. P093, Sz 11 -5 -4 -3 -2 -1 0 (Sz) 1 C1 1 1 0.91 0.95 0.99 1 0.93 C2 0.82 0.89 0.96 0.91 0.89 0.85 0.98 C3 0.95 0.55 0.80 0.70 0.46 0.46 0.97 Seizures 4, 5, 6 & 7: Spatio-temporal clustering analyses, similar to the one described on seizure 11 were performed on several other seizures, of the same patient P093. The cluster-similarity matrices P obtained from time-intervals surrounding seizures 4 & 5 and 6 & 7 of patient P093 are shown in fig 4-5. & fig 4-6. resp ectively. Channel groupings for the same are listed in tables 4-3 & 4-4. All the 4 se izures present very consistent groupings.

PAGE 91

78 Figure 4-4. Statistical validation of the cluste ring results. In each panel, thick lines are used to represent the profiles of the three clusters in a 30 minute time interval. The thin lines are the surrogate pr ofiles indicating between-cluster interactions. Cluster veracity can be visually verified by observing that amplitudes representing within-cluster interaction for cluster profiles are mostly higher that the amplitudes repr esenting between-cluster interaction for surrogate profiles, at each time-instance. 1. Consistent to the observation in seizure 11, we observe the temporal depth and the sub-cortical regions of the left hemisphere are always bunched together. 2. Once again, the association of ROF-LOF ar eas into same cluster suggests a strong homologous connection between the orbito -frontal areas of the brain. This observation is also in agreemen t with those in seizure 11. 3. The dendrograms once again presented 4 unamb iguous clusters in the form of RST, RTD, LST/LTD and LOF/ROF. The fusi on levels, indicatin g the strength of connection between clusters, often turn out in favor of RTD and RST to be grouped separately. Owing to the fact that we have pre-defined the number of clusters to 3, the LST, LTD, LOF & ROF channels w ill consequently get grouped into one cluster. 4. Once again, temporal changes are not very evident in the spatial patterns. However, observing fig 4-5 & fig 4-6 and their co rresponding dendrograms (not shown), the fusion levels and the topology of the conn ections change with time. These changes can only be quantified using statistical tests such as Mantel test statistics or the Double Permutation Statistics (DPS).

PAGE 92

79 Patient P092 In this section, we present the summary results of the clustering analyses performed on patient P092 suffering from a lesion in the medi al temporal lobe stru ctures of the right hemisphere. Channel configuration for the patient P092 is as follows Channels 1 to 4: {LTD1, LTD3, LTD5, LTD7} Channels 5 to 9: {RTD2, RTD4, RTD6, RTD8, RTD12} Channels 10 to 13: {LST1, LST2, LST3, LST4} Channels 14 to 17: {RST1, RST2, RST3, RST4} Channels 18 to 21: {LOF1, LOF2, LOF3, LOF4} Channels 22 to 24: {ROF1, ROF2, ROF3} Note that a separate 2 dimensional 25x25 sized EEG-SOM grid was created to model the data dynamics of P092. Post spectral clusteri ng analysis on 30 minutes data segments led to some interesting observations. Fig 4-7. shows the dendrograms created for seizure segments 2 hours prior to seizure 1 and 30 minutes pre-seizure, respectivel y. As before, the number of clusters (n1) specified in the spectral-clust ering step after SOM-SI bloc k was fixed to 3. The fusion levels between most of the channel clusters is greater than 0.4, indi cating a lack of strong connectivity between regions.

PAGE 93

80 Table 4-3. Spatio-temporal groupings as obt ained for seizures 4 and 5 of patient P093. P093, Seizure 4 & 5 C1 C2 C3 Pre seizure 4, (30 – 60 mins) RTD, RST LOF, ROF LTD, LST Pre seizure 4, (0 – 30 mins) RTD, RST LOF, ROF LTD, LST Post seizure 4, (0 – 30 mins) RTD LTD, LST, LOF, ROF RST Post seizure 4, (30 mins – 1 hr) RTD LOF, ROF LTD, LST, RST Pre seizure 5, (30mins – 1 hr) RTD LTD, LST, LOF, ROF RST Pre seizure 5, (0 30mins) RTD LTD, LST, LOF, ROF RST Post-Seizure 5, (30 – 1 hr) RTD LTD, LST, LOF, ROF RST Table 4-4. Spatio-temporal groupings as obt ained for seizure 6 and 7 of patient P093. P093, Seizure 6 & 7 C1 C2 C3 Post seizure 6, (0 30 mins) RTD, RSTLTD, LST LOF, ROF Pre seizure 7, (30 mins – 1 hr) RTD, RSTLTD, LST LOF, ROF Pre seizure 7, (0 – 30 mins) RTD LTD, LST, LOF, ROF RST Post seizure 7, (0 – 30 mins) RTD LTD, LST, RST LOF, ROF Post seizure 7, (30mins – 1 hr) RTD LTD, LST, LOF, ROF RST Post seizure 7, (1 hr – 1.5 hrs) RTD LTD, LST, LOF, ROF RST

PAGE 94

81 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 (1) (0) (-2) (-4) (-1) (-3) (-5) Figure 4-5. Seizures 4 and 5 of patient P093: Cluster-simila rity matrices indicating the probability that two channels share the same cluster label in a 30 minute timeinterval.

PAGE 95

82 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 LTD3 LST1 LOF1 RTD4 RST1 ROF1 (-2) (-1) (0) (1) (2) (3) Figure 4-6. Seizures 6 and 7 of patient P093: Cluster-simila rity matrices indicating the probability that two channels share the same cluster label in a 30 minute timeinterval.

PAGE 96

83 0 0.1 0.2 0.3 0.4 0.5 LTD3 LTD5 LTD1 LTD7 LST1 LST3 LST4 LST2 LOF3 LOF4 LOF1 LOF2 ROF1 ROF2 ROF3 RTD2 RST1 RST2 RST3 RTD4 RTD6 RTD8 RTD12 RST4 Fusion Level P092, 2 hours before Seizure 1 0.1 0.2 0.3 0.4 0.5 0.6 LTD3 LTD5 LTD1 LTD7 LST1 LST2 LST3 LST4 RST4 RTD2 RST1 RST2 RST3 RTD4 RTD6 RTD8 RTD12 LOF1 LOF2 LOF3 ROF2 ROF3 ROF1 LOF4 Fusion Level P092, 30 minutes before Seizure 1 Figure 4-7. Dendrograms corresponding to P092, Seizure 1. Top: 2 hours before Seizure Bottom: 30 minutes pre-seizure. For the second level of clustering, as before, let the number of clusters n2 be fixed at 3. Cluster analysis on the 30 minutes segment 2 hours prior to seizure 1 (top dendrogram in Fig. 4-7) results in the following groups of channels: Cluster #1: LTD & LST Cluster #2: RTD & RST Cluster #3: LOF & ROF

PAGE 97

84 Observe the cluster formed from LTD & LST channels, in the dendrogram. It is made up of two sub-clusters, a large and a small clus ter. The small cluster consists of only two channels, LTD (3 & 5) and fuses with the ot her sub-cluster at a ve ry high fusion level (implying weak link). If n2 was to be increased to 4, the clustering algorithm would classify this sub-cluster as an independent cl uster. A detailed analysis on all seizures in P092 revealed a strong intra-channel correlati on (or low fusion level) between channels LTD (3 & 5) and a weak inter-channel correla tion with the rest of the channels. Surrogate analysis also confirmed the imbalance by ha ving very few rejections for the cluster consisting of LTD (3 & 5) channels. It is obvious that the average interaction (withincluster interaction) of the largest cluster would be pulled down if there are sub-clusters that have a strong within-s ub-cluster interaction, but a weak between-sub-cluster interaction. Consequently, the within-cluster interaction of the largest cluster can be expected to be as weak as or marginally better than the between-cluster interactions, leading to fewer rejectio ns of the null hypothesis H0. This problem can possibly be overcome by increasing the number of clusters to 4 or more. However, for consistency, we let the number of clusters n2 be fixed at 3 in the rest of the analyses.

PAGE 98

85 LTD1 LST1 LOF1 RTD2 RST1 ROF1 LTD1 LST1 LOF1 RTD2 RST1 ROF1 LTD1 LST1 LOF1 RTD2 RST1 ROF1 LTD1 LST1 LOF1 RTD2 RST1 ROF1 LTD1 LST1 LOF1 RTD2 RST1 ROF1 LTD1 LST1 LOF1 RTD2 RST1 ROF1 LTD1 LST1 LOF1 RTD2 RST1 ROF1 LTD1 LST1 LOF1 RTD2 RST1 ROF1 LTD1 LST1 LOF1 RTD2 RST1 ROF1 LTD1 LST1 LOF1 RTD2 RST1 ROF1 LTD1 LST1 LOF1 RTD2 RST1 ROF1 LTD1 LST1 LOF1 RTD2 RST1 ROF1 ( 3 ) ( 1 ) ( 1 ) ( 2 ) ( 0 ) ( 2 ) Figure 4-8. Cluster-similarity matrices indi cating the probability that two channels share the same cluster label in a 30 minute time-interval. a) Seizure 1 of patient P092 b) Seizure 3 of patient P092.

PAGE 99

86 Seizures 1, 3, 4, 5, 6 & 7: The spatio-temporal clustering results for se izures 1, 3, 4, 5, 6 & 7 are summarized in tables 4-5 to 4-10. Table 4-5. Spatio-temporal groupings as obtained for seizure 1 of Patient P092. P092, Seizure 1 C1 C2 C3 Pre seizure, (1.5 – 2 hrs) RTD, RSTLTD, LST (1,3,4) LOF, ROF, LST (2) Pre seizure, (1 – 1.5 hrs) RTD LST, RST, LOF, ROF, LTD (1,7) LTD (3,5) Pre seizure, (30 mins – 1 hr) RTD, RSTLTD, LST LOF, ROF Pre seizure, (0 30 mins) RTD, RSTLTD, LST LOF, ROF Post-Seizure, (0 – 30 mins) RTD, RSTLTD, LST LOF, ROF Post-Seizure, (30 – 1 hr) RTD LTD, LST, LOF, RST ROF Table 4-6. Spatio-temporal groupings as obtained for seizure 3 of Patient P092. P092, Seizure 3 C1 C2 C3 Pre seizure, (1.5 – 2 hrs) RTD LST, LTD, RST LOF, ROF Pre seizure, (1 – 1.5 hrs) RTD LST, LTD, RST LOF, ROF Pre seizure, (30mins – 1 hr) RTD, RST LST, LTD LOF, ROF Pre seizure, (0 30mins) RTD, RST LST, LTD LOF, ROF Post-Seizure, (0 – 30 mins) RTD, RST LST, LTD LOF, ROF Post-Seizure, (30 – 1 hr) RTD, RST LST, LTD LOF, ROF

PAGE 100

87 Table 4-7. Spatio-temporal groupings as obtained for seizure 4 of Patient P092. P092, Seizure 4 C1 C2 C3 Pre seizure, (1.5 – 2 hrs) RTD, RST LST, LTD LOF, ROF Pre seizure, (1 – 1.5 hrs) RTD, RST LST, LTD LOF, ROF Pre seizure, (30mins – 1 hr) RTD LST, LTD, RST LOF, ROF Pre seizure, (0 30mins) RTD, RST LST, LTD LOF, ROF Post-Seizure, (0 – 30 mins) RTD, RST LST, LTD LOF, ROF Post-Seizure, (30 – 1 hr) RTD, RST LST, LTD LOF, ROF Table 4-8. Spatio-temporal groupings as obtained for seizure 5 of Patient P092. P092, Seizure 4 C1 C2 C3 Pre seizure, (1.5 – 2 hrs) RTD LST, LTD, LOF, RST ROF Pre seizure, (1 – 1.5 hrs) RTD LST, LTD, LOF, RST ROF Pre seizure, (30mins – 1 hr) LTD (3,5) LST, LTD (1,7), RST, RTD LOF, ROF Pre seizure, (0 30mins) LTD (3,5) LST, LTD (1,7), RST, RTD, LOF ROF Post-Seizure, (0 – 30 mins) RTD, RST LST, LTD LOF, ROF Post-Seizure, (30 – 1 hr) RTD, RST LST, LTD LOF, ROF

PAGE 101

88 Table 4-9. Spatio-temporal groupings as obtained for seizure 6 of Patient P092. P092, Seizure 4 C1 C2 C3 Pre seizure, (1.5 – 2 hrs) RTD, RST LST, LTD LOF, ROF Pre seizure, (1 – 1.5 hrs) RTD LST, LTD, RST LOF, ROF Pre seizure, (30mins – 1 hr) RTD LST, LTD, RST LOF, ROF Pre seizure, (0 30mins) LTD (3,5) LST, LTD, RTD, RST LOF, ROF Post-Seizure, (0 – 30 mins) RTD, RST LST, LTD LOF, ROF Post-Seizure, (30 – 1 hr) RTD, RST LST, LTD, LOF ROF Table 4-10. Spatio-temporal groupings as obtained for seizure 7 of Patient P092. P092, Seizure 4 C1 C2 C3 Pre seizure, (1.5 – 2 hrs) RTD, RST LST, LTD LOF, ROF Pre seizure, (1 – 1.5 hrs) RTD, RST LST, LTD LOF, ROF Pre seizure, (30mins – 1 hr) RTD LST, LTD, RST LOF, ROF Pre seizure, (0 30mins) RTD, RST LST, LTD LOF, ROF From the P092 cluster results, we note the following: 1. The non-focal zone LTD has a strong coupling with the LST region. Correspondingly, strong affinity is observ ed between RTD and RST as well. These observations are consistent with the obse rvations for P093. However, unlike in P093, we also see here that LTD conn ects and disconnects wi th several other channels, depending on the seizure state. 2. As in P093, we observe an exclusiv ely strong connection between ROF-LOF regions at all stages surrounding a seizure. There are few instances where the ROF breaks into a separate group. We do not have any explanation for this drift in ROF, at this point in time.

PAGE 102

89 3. Statistics from the surrogate analyses c onfirmed the veracity of the technique in most of the cases. As pointed out earlier, discrepancies occurred in a few instances for the clusters containing LTD (3, 5) channels. Finally, we summarize the analysis on 2 pa tients and 11 complex partial seizures: 1. Contrary to the stand point that the seiz ure activity initiates in the focal zone followed by a gradual propagation to other regions, we observed that the spatial organization reflected by intracranial EEG activity exhibits either minimum or no progressive changes from the focal zone (RTD) to other zones (based on how it groups with other regions in the brain). 2. Evidence show stronger ipsi lateral connection between the LTD and LST zones compared to the connection strength betw een RTD-RST. Statistical analysis to check if a significant difference in intrahemisphere coupling strengths exists is needed. 3. We also found evidence to show a str ong cross-hemispheric activity by observing consistent groupings of the right and left or bito-frontal lobes at all seizure states. 4. On how the overall spatial networks ch ange every 30 minutes, patient P093 was seen to have qualitatively, lesser spatio-temporal changes in its P matrices than P092. It remains to be checked whethe r a significant cha nge in the spatial organization before seizure is a pre-requi site to its subsequent occurrence. (See chapter 5 for further investigation). 4.4 Summary From a clinical standpoint, one of the objectives of the a bove analysis, although qualitative, was two fold: 1) to have a better understanding of the spatial patterns in an epileptic brain and 2) to fo rge a link between the spatio-t emporal organization and the evolution of seizure states. In cognitive term s, it would be particularly interesting to know if an epileptic brain exhi bits certain predicta ble pattern differentia ls in its spatial configuration en-route to seizur es. Such investigations would help us determine if it is possible to extract markers pointi ng to occurrence of seizures. Although our first objective could be accomplished to a certain extent, analysis on 11 seizures from 2 patients failed to presen t any predictable and consistent changes, suggestive of a seizure marker. Qualitativ e comparisons between 30 minute pre-seizure

PAGE 103

90 block with 30 minute post-seizure blocks ar e very inconclusive; therefore, we can only remark to the extent that the changes in spatial groupings surrounding seizures have no consistent and predictable patt erns. Referring again to the co nsistent groupings among the brain sites regardless of the seizure states, we wonder if the grouping patterns are finger prints of an epileptic brain. In other word s, are the results compelling enough to point to an abnormality in the brain? Similar cl ustering analyses on a normal brain and a subsequent comparison with analyses from an epileptic brain might be a worthwhile effort to make these hypothetical arguments complete

PAGE 104

91 CHAPTER 5 STATISTICAL ASSESSMENT OF SPATIO-TEMPORAL DEPENDENCY CHANGES IN EPILEPTIC INTR ACRANIAL EEG 5.1 Introduction In the previous chapter, the proposed spa tio-temporal clustering model was used to analyze the groupings among critical regions in an epileptic brain. Remarkably, certain regions in the brain displayed tendencies to strongly bind with each other, regardless of the clinical states. The model is useful as an indicator of the gr oupings in the brain; however, it does not quantify the exact changes in the overall associability patterns that are speculated to cause seizures. In the next phase of this study therefore, we explore the possibility whether seizure events are coupled with changes in the f unctional relationships among certain groups of channels. The chapter is split into two inde pendent parts. The first half discusses a statistical approach for analyzing the temporal changes in the overall spatial-patterns of an epileptic intracrania l EEG. In particular, mantel statis tics are used for comparing the similarities between affinity matrices, eval uated at different seiz ure states. Significant changes of mutual interactions in 2 hours wi ndows prior to seizures are reported. The second half of the chapter determines the st atistical differences in the synchronization levels between the focal and the non-focal hemi spheres at pre-ictal and post-ictal states. 5.2 Spatio-temporal Changes It is widely believed that the dynamics of the brain, understood as a multivariate biological system, can be characterized by an alyzing its temporal state changes. Such

PAGE 105

92 changes are known to eventually lead to certa in clinical manifestations. Much of the previous studies [ 28, 35-40, 77, 82 ] in epilepsy have focused on analyzing the temporal changes associated with brain’s non linea r dynamics. Feature descriptors such as system’s complexity or the short-term Lyapunov exponents [ 35-40 ] are derived individually on each of the systems dimensi ons. Temporal dynamical changes associated with such features, however, fail to explain a system’s state associated changes in its overall spatial configurati on. Rather than studying the temporal dynamics of each dimension individually, the emphasis should be on considering all the dimensions in unison in a multi-variate scheme of things. In other words, in spatially extended systems, dynamics change both in time and in space and therefore, approaches that track the temporal changes of the spatial networks can be much more effective and helpful in our efforts to characterize cer tain clinical events. Earlier, we derived the SOM-SI measur e to define mutual interactions among various nodes in a spatially coupled multi-dimensional system. The affinity matrix representation of the SOM-SI at every time window (of length 10 seconds) provides information on the interactions among all the possible pairs of nodes in a graph. For the epileptic brain, in particular, changes associ ated with the epileptic activity will therefore be reflected by changes in the overall spat ial connectivities and we propose to quantify this by tracking the activity changes associat ed with the SOM-SI affinity matrices. 5.3 Mantel Test of Matrix Comparison Mantel tests were first developed in 1967 to correlate temporal and spatial distributions of cancer incidences [ 58 ] and since then it has been widely used as a correlation tool in various biological [ 96 ] and ecological disciplines [ 18, 50-51 ]. It is a linear correlation estimate of the relationship between two square distance matrices based

PAGE 106

93 on the degree of relationship of two sets of variables taken at the same sampling locations. In short, mantel test is essentially a statistical framework to test the consensus of two distance/proximity/affinity matrices. In the mantel test, the hypothesis is that the distances (or similarities) in matrix A are independent of the distances, for the same se t of objects, in another matrix B. In other words, we test the hypothesis that the two matrices under st udy are no more similar than they would be by chance assignment of th e labels to the rows and corresponding columns. The normal procedure to test the hypo thesis would be to compute a measure of resemblance between the values in the two uppe r (or lower) triangular parts of the square symmetric matrices under comparison and test against a random di stribution. The random distribution is constructe d by repeatedly permuting at random, the rows and corresponding columns of one of the matrices, and re-computing the statistic. Finally, the original value of the statis tic is compared with the di stribution obtained by randomly reallocating the order of the elements in one of the matrices. The statistic used for the measure of correlation between the matrices is the classical Pearson correlation coefficient N j B ij A ij N is B B s A A N r1 11 1 (5.1) where N is the number of elements in the lo wer or upper triangular part of the matrix, A is the mean for A elements and sA is the standard deviation of A elements. If the two matrices are normalized, i.e. ; B ij ij A ij ijs B B b s A A a

PAGE 107

94 we have 1 0 1 0 B As b s a and therefore (1) can be re-written as N j ij ij N ib a N r1 11 1 (5.2) Note that the coefficient measures the linear correlation coefficient; therefore if any non linear relationships exist, they will be lost. Finally, the testing procedure using the mantel test statistic is as follows. Assume two square symmetric matrices A and B of size n x n. The rows and columns in both the matrices correspond to the same objects. The first step is to compute the Pearson correlation coefficient between the corres ponding elements of the lower (or upper) triangular matrices. 1. Compute the original (or th e non-permuted) statistic rAB using (5-1) 2. Permute randomly the rows and the corres ponding columns of one of the matrices (say B) to create a new matrix B’. 3. Recompute (5-1) using A and B’ to obtain permuted statistic 'ABr 4. Repeat the steps (2) and (3 ) several times (> 500) to form a distribution of the permuted statistics. This distribution w ill be the reference distribution under the null hypothesis. 5. Assuming normal approximation on the refere nce distribution, compute the z-score by comparing the non-permuted statistic in st ep (1) with the mean and variance of the reference distribution. i.e, ''rr r z 5.4 Application to Epileptic Intracranial EEG data Recall from the 3rd chapter that the similarity indices on 10-second windows were evaluated over several hours of epileptic intracran ial EEG recording. In this study, we propose to use the mantel-test on the SOM-SI affinity matrices, evaluated at different time periods. Particularly, the emphasis will be on tracking the temporal changes of the

PAGE 108

95 spatial connections, in the intervals prior to seizure. Experimental design procedure is explained in the following steps. 1. Select 2 hours of intracranial EEG segments prior to a seizure. 2. Compute the SOM-SI affinity matric es on 10 second windows, for 2 hours preseizure (only the alternate 10 seconds are ev aluated; therefore we get a total of 180 SOM-SI affinity matrices from 2 hours data). 3. Fluctuations between successive matrices are smoothened by temporal averaging. A window size of 3 (equal to 1 minute) was used. 4. Affinity matrix corresponding to the window 2 hours prior to seiz ure is individually compared with the affinity matrices at all other times leading to seizure. For comparison, every 3rd minute was used to ensure that the eff ects of correlation between matrices due to temporal averag ing were eliminated. Note that the nullhypothesis distribution was formed from 700 randomly permuted statistics. 5. Hypothesis testing is done by checking the zscores at the 95% significance level. The mantel test statistic was tested on si x seizures of the patient P092. Fig. 5-1 illustrates the temporal changes in the spatialactivity of the channels, in the interval 2 hours prior to seizure. z scores greater than zcrit,0.95 = 1.96 indicates that the null hypothesis H0 is rejected at the 95% significance level. Rejecting H0 implies that the similarity statistic between the test-affinity matrix at time‘t’ and the reference affinity matrix (corresponding to 2 hours before seizure) is significantly di fferent than the one obtained by randomly permuting the rows and columns of the test-affinity matrix. In Fig.1 the z scores in all seizures are al most always greater than zcrit. However, in a few instances, the z ‘s exhibit a tendency to decrease gradua lly as the seizures approach. In a few other seizures, z’s seem to have a negative bump that lasts for several minutes. Reduced z’s do not necessarily imply that the correlation estimates are small. However, verification of the original similarity statis tic (r) indicated a redu ced correlation estimate

PAGE 109

96 at those points corresponding to reduced z’s. This observation dire ctly suggests that the spatial relationships in the intracranial EEG data at those points are indeed very different (less correlated) from the spatial relationships observed around 2 hours prior to seizures. 5.5 Statistical Approach As stated earlier, even though the z’s show a remarkable decrease they are still greater than the zcrit. Also, notice that the absolute values of z’s vary across seizures. These non-uniformities will render any comments regarding spatio-temporal changes prior to seizures, purely qualita tive. It is therefore absolutely necessary to quantify the temporal decrease observed in correlation es timates. We propose a simple approach to statistically verify the decreas e in the mantel statistics, z’s. The approach consists of checking whethe r the decreases observed in the mantel test procedure are statistically significant or not. Let zref be the reference z score at a time instance close to 2 hours before seizure. If , _ref t diff tz z z the null hypothesis H0 can be stated as follows Null Hypothesis H0: There is no difference in the z-scores, , and ref tz z. The alternate hypothesis H1: The z-scores , and ref tz z are not equal. Testing the hypothesis cons ists of following steps 1. Construct a distribution of z’s from samples taken during inter-ictal states. Specifically, take any 2 hour segment from the background activity of the intracranial EEG data. Construct affinity matrices as before and evaluate mantel statistics to get z-scores. Repeat the above procedure on a number of such 2 hour segments (say 40) from background record ing of the same patient and evaluate mantel test procedure for each one of them, to form an ensemble of z-scores, i s tz,where t represents the time-index (hours), s represents the segment index (s = 1, 2, 3….., 40) and i is used to denote that these zscores are computed on inter-ictal segments. For illustration, fig. 5-2 shows the smoothened iz profiles

PAGE 110

97 corresponding to the inter-ictal segments. The smoothened z profile corresponding to seizure 1 is also shown superimposed on the iz profiles. 1. If i s t i s ref i s diff tz z z, , is the difference in the z-scores at time t relative to a reference time for the sth segment, then i diff tz_ forms a distribution that is constructed from all the 40 segments and is observed to be approximately Gaussian. 2. The idea is to check if the difference obs erved in the z-scores in the segment 2 hours before seizure (zt_diff) is significantly different fr om the differences observed in the z-scores during inter-ictal states (i diff tz_). We quantify this idea as follows ) (_ _i diff t i diff t diff t tz z z z (5.3) where i diff tz_ and ) (_i diff tzare the mean and variances of, i diff tz_respectively. Zt > 1.96 indicates the difference is significant a nd so the null hypothesi s can be rejected at a 95% significance level. Since the test is one-tai led, the null hypothesis cannot be rejected if Zt< -1.96. Even though it means that the differences are significant in those cases, a negative , ref t diff tz z z implies an upward curve in fig 5-1. as the seizure is approached. Seizure 6 in fig. 5-1 can be one such instance. Fig. 5-3 presents the Zt scores as a function of time, fo r seizures 1 to 7 of the patient P092. Zt scores were computed for three different reference times. In other words, the zscores at 90 minutes, 100 minutes and 110 mi nutes prior to a seizure were used as, refz Fig 5-3. through to Figs 5-5 show profiles of Zt scores for all the three reference times. Since we are particularly interested in Zt scores > 1.96 (critical threshold), we make the following observations: 1. 6 out of the 7 seizures (the exemption be ing seizure 3) have clear time-instances where the Zt scores are greater than the critical threshold. 2. The time-instances where critical threshol ds are crossed vary from seizure to seizure. 3. The Zt profiles corresponding to the three reference times (90, 100 & 110 minutes) are mostly consistent with respect to crossing the thresholds. Occasional discrepancies are seen, perhaps due to rapi d fluctuations within the z-scores at reference times.

PAGE 111

98 From the 1st observation above, we have clear eviden ce of spatial activit y changes in the intracranial EEG data as the patient approach es seizure state (or a pre-seizure state!!). The 2nd observation tells us that the times at which the statistical changes are observed vary across seizures. This can lead us to conc lude that the seizure markers in the form of spatio-temporal pattern changes perhaps are a function of the t ype of seizures and also the pathological state of the brai n and many other variables, as well. There is also an ongoing debate on the notion of a pre-se izure state and the transition from inter-ictal state to preictal state. The above analyses may be a step ahead in providing ev idence to the existence of pre-seizure states and tran sition between states. Significan t changes in the intracranial EEG’s spatio-temporal patterns can possibly serve as one of the precursors for an impending seizure. 5.6 Inter-Hemisphere Synchronization Differences In chapter 3, we addressed the problem of quantifying mutual interactions among various regions in the brain by proposing the nearest-neighbor based SOM-SI technique. The SOM-SI is fundamentally a state space approach that uses the recurrence in dynamical states of one signal to quantify th e recurrence of dynamical states in another signal. The similarity indices, in fact, refl ect the activity shared between the neuronal networks at the sub-cortical levels.

PAGE 112

99 1.5 2 2.5 3 0 2 4 6 8 10 12 14 time (hours)ZSeizure 1 3.5 4 4.5 5 5.5 0 5 10 15 time (hours)ZSeizure 2 4.5 5 5.5 6 0 5 10 15 time (hours)ZSeizure 3 10.5 11 11.5 12 0 5 10 15 time (hours)ZSeizure 4 1 1.5 2 2.5 3 3.5 0 5 10 15 time (hours)ZSeizure 6 10 10.5 11 11.5 12 0 5 10 15 time (hours)ZSeizure 5 Figure 5-1. Quantifying the spatial change s along time using Mantel statistics. The vertical lines in the figure indicate seizure onset and termination periods. Notice that the statistical values show a slight decrease, approaching the seizures.

PAGE 113

100 1 1.5 2 2.5 3 0 5 10 15 time (hours)z Figure 5-2. Time smoothened profiles of the Mantel statistics duri ng inter-ictal states. The dark line shows the profile for seizure 1, P092. The mesial lobe temporal lobe epilepsy is mainly characterized by partial seizures that are both complex and generalized in na ture. While the generalized seizures are known to be widespread across both the hemi spheres, the complex partial seizures are more localized. The complex pa rtial seizures are known to be particularly active in the hippocampal networks that contain damage d lesions. Based on that, it is argued [ 46-47 ] that the neuronal networks co rresponding to the focal zone have a higher synchronizing ability, compared to the networks in the non-fo cal zone. It has also raised speculations on whether the interaction levels over the enti re focal-hemisphere are higher than the corresponding interaction levels in the non-focal hemisphere. A recent study by Kraskov et al. [ 46-47 ], showed that, in the inter-ictal peri od, the average synchronization in the focal hemisphere of the brain is indeed significantly higher than in the non-focal hemisphere.

PAGE 114

101 -1.5 -1 -0.5 0 -2 -1 0 1 2 3 time before seizure (hours)ZSeizure 1 Reference: 90 mins before seizure Reference:100 mins before seizure Reference:110 mins before seizure -1.5 -1 -0.5 0 -3 -2 -1 0 1 2 3 time before seizure (hours)ZSeizure 2 -1.5 -1 -0.5 0 -3 -2 -1 0 1 2 3 time before seizure (hours)ZSeizure 3 Figure 5-3. Illustrating the st atistical difference between z-sc ore at time‘t’ and z-score at a reference time, 90, 100 and 110 minutes before a seizure. Top to Bottom: profiles of z-score statistics for seiz ures 1 through 3 of patient P092.

PAGE 115

102 -1.5 -1 -0.5 0 -3 -2 -1 0 1 2 3 4 time before seizure (hours)ZSeizure 4 -1.5 -1 -0.5 0 -3 -2 -1 0 1 2 3 4 time before seizure (hours)ZSeizure 5 -1.5 -1 -0.5 0 -4 -3 -2 -1 0 1 2 3 time before seizure (hours)ZSeizure 6 Reference: 90 mins before seizure Reference:100 mins before seizure Reference:110 mins before seizure Figure 5-4. Illustrating the st atistical difference between z-sc ore at time‘t’ and z-score at a reference time, 90, 100 and 110 minutes before a seizure. Top to Bottom: profiles of z-score statistics for seiz ures 4 through 7 of patient P092.

PAGE 116

103 -1.5 -1 -0.5 0 -2 -1 0 1 2 3 time before seizure (hours)ZSeizure 7 Figure 5-5. Illustrating the st atistical difference between z-sc ore at time‘t’ and z-score at a reference time, 90, 100 and 110 minutes before a seizure. Top to Bottom: profiles of z-score statistics for seizure 7 of patient P092. In this study, we investig ate on the differences in the synchronization levels between the focal and the non-focal hemisphere s at pre-ictal and post-ictal states. Recall from the previous chapter that we observe d that a sudden increase in synchronization levels was accompanied by a gradual decrease at post-seizure, in al l the regions of an epileptic brain. Granted that the brain consistently observes an overall increase in synchronization post-seizure, it would be interesting to know whether the difference in synchronization levels between the two hemis pheres continue or cease to exist. These questions bring us to the fo llowing hypothesis formulation: Define the average interaction between th e Right Temporal Depth (RTD) and all the other regions in the right hemisphere as Ravg. Similarly, define the average interaction be tween the Left Temporal Depth (LTD) and all the other regions in the left hemisphere as Lavg. H0: Average synchronization difference between Ravg and Lavg is the same from pre-ictal to post-ictal state.

PAGE 117

104 H1: Average synchronization difference between Ravg and Lavg undergoes a change going from pre-ictal to post-ictal state. 5.7 Statistical Assessment To test the hypothesis, the average synchronization index Ravg from the focal hemisphere is compared against the Lavg, from the non-focal he misphere. As defined earlier, Ravg is obtained by averaging the mutual interaction between the right-temporal depth channels and the other channels (nam ely Right Orbito Frontal and the Right SubTemporal regions) in the right hemis phere. The same definition holds for Lavg obtained from the left hemisphere channels. Comparisons consist of checking whether pairs of Ravg and Lavg synchronization indices, obtained from a group of successive windows, come from same distributions or not. Statistical significance is checked by using the Wilcoxon signed-rank test. It is a nonparametric procedure to check the null hypothe sis that the two sets of variables come from same distribution or not. The test does not make any apriori assumptions about the distribution of the variables. However, it takes into account the magnitude of differences within pairs and gives more wei ght to pairs that show large di fferences than to pairs with small differences. The test statistic is based on the ranks of the absolute values of the differences between the two variables. An observed significance leve l is often called the p-value. This value is the basi s for deciding whether or not to reject the null hypothesis. It is probability that a statistical result as extreme as the one observed would occur if the null hypothesis were true. If the observed signi ficance level is small enough, usually less than 0.05 or 0.01, the null hypothesis is rejected.

PAGE 118

105 5.8 Experimental Results Intracranial EEG recordings consisting of 4 seizures from patient P093 and 6 seizures from patient P092 are analyzed, retros pectively. Both the patients had a history of temporal lobe epilepsy in the mesial stru ctures. Following surgic al resection, a lesion had been identified in the anterior RTD elect rodes region. In the rest of the chapter, therefore, the right hemisphe re would be called as the focal hemisphere and the RTD region as the epileptic zone. Co rrespondingly, the left hemis phere would be called as the non-focal hemisphere. The definition of existence of a pre-seizure state is very subjective. Even more, the discussion on how long before seizures is a pr e-seizure state differs very widely. Without getting into the intricacies of detecting inter-i ctal to pre-ictal transitions, we just fix the window of analysis before seizures to 60 mi nutes. Due to the fact that the arguments on post-seizure durations also have no general consensus, we trea t the entire 60 minutes of post-seizure as post-seizure state. Synchronization analysis consisted of an alyzing 60 minutes pre-seizure and 60 minutes post-seizure onset of intracranial EEG data, for each of the 10 seizures from 2 patients. For each non-overlapped 10-sec window (2000 samples, 200 samples/sec digitization rate), SOM-similarity index meas ure was applied to estimate the pair-wise synchronization index between the set of 24 ch annels. Note that, for computation reasons, only alternate 10-second windows were analyzed. 5.8.1 Experiment #1 The average synchronization indices Ravg and the Lavg are computed for every 10second window, over the entire time range of 120 minutes (60 minutes pre and 60 minutes post-seizure interval). A 10 minute m oving window is then used for statistical

PAGE 119

106 validation of synchronization indices, using a paired signed rank test. Specifically, in each sliding window of 10 minutes, the test ve rifies whether the distribution formed from computing the paired difference between Ravg and Lavg indices have zero median or not. Null hypothesis is rejected if the test produces a score outside the critical score for 95% significance. For example, z > zcrit, 0.05 (1.96 according to the table), will reject the null hypothesis, implying that the Ravg and Lavg indices have different medians. Successive windows overlap; the amount of sliding being 20 seconds. Results of a signed-rank test on one of the seizures are presented in Fig 56. Consistent with our earlier observations in previous chapter, the average synchronization indices exhibit a seeming increase in synchronization post-seizure, lasting for seve ral minutes. Results from the paired signedrank test show significant pre-seizure sync hronization differences. This seizure in particular, was a complex partial seizure implyi ng that the seizure is more hemispheric or rather, localized. The pre-seizure statistical results provide convincing evidence on the nature of localization by pointi ng out to the differences in the interactions between the two hemispheres. The situation post-seizure, however, is different. The statistical results fail to reject the null-hypothesis, for a peri od lasting for about 30 minutes post-seizure. Between the 30 minutes and 60 minutes interv al post-seizure, we again observe nullhypothesis rejections. On a larger group of seizures (10 seizures from 2 patie nts), synchronization results were analyzed using the procedure described above. Define H to be the index of nullhypothesis evaluation. H = ‘1’ indicates a rejection of null-hy pothesis and H = ‘0’ indicates insignificant evidence to do so. The average H from 10 seizures, corresponding to a 10 minute time-window are computed and plotted in fig 5-7. Clearly, the 60 minutes

PAGE 120

107 interval prior to seizure onset seems to have high rejection percen tage. It highlights the degree of difference in the way the focal z one (RTD) interacts w ith all other righthemisphere channels and the way the non-focal zone (LTD) interacts with other lefthemisphere channels. Post-seizure results exhi bit a sudden drop in H, from seizure onset to 6-8 minutes after seizure. The low rejecti on rate implies a higher degree of similarity between Ravg and Lavg. 20 minutes after seizure onse t, we again observe a gradual transition from strong similarity between Ravg and Lavg to a strong difference between them. In summary, the statistics suggest th at, in the first 20 minutes post-seizure, synchronization differences between Ravg and the Lavg cease to exist implying a global synchronization without any bi as towards hemispheres. 5.8.2 Experiment #2 In the first experiment, analysis was primarily focused on testing the synchronization differences between Ravg and Lavg during the pre-ictal and post ictal states. While the test results clearly indicate that the RTD interactions with other right hemispheric channels are similar to the corresponding LTD interactions in the post-ictal state, it will be interesting to check whether the RTD interactions with channels in the left hemisphere also yield similar results. With this objective, we set up the hypothesis as follows: Define the average interaction between th e Right Temporal Depth (RTD) and all the other regions in the right hemisphere as r avgR. Similarly, define the average interaction between the Right Temporal Depth (LTD) and all the other regions in the left hemisphere as l avgR

PAGE 121

108 H0: Average synchronization difference between r avgR and l avgR is the same from preictal to post-ictal state. H1: Average synchronization difference between r avgR and l avgR undergoes a change going from pre-ictal to post-ictal state. In tune with the design procedure of experiment #1, r avgR and l avgR synchronization profiles are computed from the SOM-SI in dices, following which the wilcoxon signed rank test is applied to test the hypothesis. Defining H as before, the average H from 10 seizures is computed and plotte d as a function of time in fig 5-8. As in experiment #1, the 60 minutes interval before seizure has a high null-hypothesis rejection percentage, suggesting the difference between the RTD inte raction with right hemisphere channels and the corresponding left hemisphere channels Interestingly, the interval post-seizure also exhibits a high rejection rate. It imp lies that significant sync hronization differences continue to exist between r avgR andl avgR, as the brain transits between seizure states. Effectively, this may mean that the na ture of the homologous interactions (or connectivities) of the focal zone are very different compared to its ipsilateral interactions (or connections), at all ictal states. 5-9 Summary We addressed the question if the streng th of connections between non-focal channels in the focal hemisphere (RST and ROF) and the channels in focal zone (RTD) differed significantly from the connections between non-focal channels in the non-focal hemisphere (LTD and LOF) and channels in non-focal zone (LTD). Hypothesis testing reported significant degree of differences in th e pre-seizure states (1 hour before seizure). However fewer rejections of null-hypothesis in the time-intervals 20 minutes post-seizure

PAGE 122

109 led us to concede that the non-focal intera ctions may be highly entrained (both in amplitude and phase) with that of the focal interactions following seizures. In a similar setup aimed at determining the functional diffe rences between (RST and ROF) with focal zone RTD and the (LST and LOF) with RTD, we observed that with reference to the focal zone, the nature of cross-hemispheric connections always differed from the intrahemispheric connections at all seizure states. -60 -40 -20 0 20 40 60 -5 -4 -3 -2 -1 0 z scores -60 -40 -20 0 20 40 60 -0.2 0 0.2 0.4 0.6 0.8 Time (minutes)Average Synchronziation index RTD / (RST+ROF) LTD / (LST+LOF) zcrit,0.05Time (minutes) Seizure onset Figure 5-6. Results of synchronization anal ysis on seizure 1 of P092, to check the differences between the focal and non-fo cal hemispheric interactions, at pre and post-ictal states. Top: Average SO M-SI synchronization indices of the right (focal) and left (non-focal) hemisphe res, 1 hour before and after seizure. Bottom: Results of Wilcoxon signed rank test show that the differences in hemispheric interactions cease immediatel y after seizure and lasts for several minutes post-seizure.

PAGE 123

110 -60 -40 -20 0 20 40 60 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 time (minutes)Average rejections of Null Hypothesis Figure 5-7. Null-hypothesis reje ctions rate averaged over 10 seizures for experiment 1 The vertical line at 0th minute indicates seizure onset. 60 minutes interval prior to seizure onset seems to have cons iderable difference between the intrahemispheric interactions. Post-seizu re results suggest a very minimal difference between the intra-hemisphere interactions, approximately in the first 20 minutes. Later than 20 minutes seems to be a transition phase from strong similarity between hemispheres to a strong difference between them.

PAGE 124

111 -60 -40 -20 0 20 40 60 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 time (minutes)Average Rejections of Null Hypothesis Figure 5-8. Null-hypothesis reje ctions rate averaged over 10 seizures for experiment 2. The vertical line at 0th minute indicates seizure onset. 60 minutes interval prior to seizure onset seems to have cons iderable difference between the intrahemispheric interactions. Post-seizu re results suggest a very minimal difference between the intra-hemisphere interactions, approximately in the first 20 minutes. Later than 20 minutes seems to be a transition phase from strong similarity between hemispheres to a strong difference between them.

PAGE 125

112 CHAPTER 6 CONCLUSIONS AND FUTURE DIRECTIONS 6-1 Conclusions With the advent of sophisticated linea r and non linear signal processing tools, application of engineering tec hniques to analyze epileptic data has taken a major step forward. A lot of advances are being seen in development of seizure detection and prevention systems. These advances have e nhanced our basic unders tanding of neural dynamics associated with seizure generation. There is sufficient evidence to believe that the brain as a functional model consists of highly complex non linear dynamics. App lication of non linear dynamical measures such as short term Lyapunov exponents (S TLmax) and correlation dimension on an epileptic brain has revealed that the comple xity of the brain reduces significantly as seizure is approached. In ot her words, the temporal dynamics of the brain progresses from a ‘severe chaotic’ state to a much le sser or rather ‘mildly chaotic’ state. Much of the analysis on temporal dynamics focuses on analyzing and characterizing the irregular beha vior of the time-signal of ei ther intracranial or scalp intracranial EEG. However, it is important to realize that the brain is a multi-dimensional system with a large set of neuronal oscill ators that are physica lly coupled together. Obviously the neurons at different spatial locations communicate with each other through synaptic potentials resulting in spike discharg es. In such an environment therefore, it would be natural to expect a spatio-temporal interaction at various patho-physiological states. Specifically, in the context of epilepsy, one can perceive that the seizure related

PAGE 126

113 abnormal spike discharges and subsequently the seizure events are accompanied by different regions in the brain communicating in a structured fashion. Therefore, it is very essential to unravel the functional connectivity of the neural networks and analyze how the structures change during seizure events. In this research, we proposed the SOM-based similarity index measure to analyze the mutual interactions among critical areas of an epileptic brain. We found that the interaction levels were very low at the inter-ictal and pre-ictal intervals. Post-seizure, where the intracranial EEG profiles are ge nerally high frequency and low amplitude waveforms, the interactions exhibited a significant positive jump. Based on the functional relationships, we analyzed long term structur al connectivity’s related to various seizure states by proposing a spatio-temporal cluste ring model. On analyzing 12 complex partial seizures from 2 patients suffering from right temporal lobe epilepsy, we found that the orbito-frontal regions always exhibit a st rong homologous connectivity while maintaining a low relationship with other regions. The left sub-temporal and the left-temporal depth regions (non-focal hemisphere) were identified to have a st rong ipsilateral connection, regardless of seizure states. Fi nally, we found that the epil eptic zone, namely the right hippocampus depth region maintained a relati vely strong connection with the right subtemporal region. Interestingly, the confi guration of the groupings between different regions always remained the same, regardless of whether the patient was in an inter-ictal, pre-ictal or post-ictal state although the inte r-region connectivity strengths seemed to vary slightly across states. In a related study, we analyzed the differe nces between connectivity strengths from left-hemisphere regions and the right-hemis phere regions to the focal zone (right

PAGE 127

114 temporal depth). Hypothesis test ing reported significant degree of differences in the preseizure states (1 hour before seizure). Howeve r fewer rejections of null-hypothesis in the time-intervals 20 minutes post-seizure led us to concede that the non-focal interactions may be highly entrained (both in amplitude and phase) with that of the focal interactions following seizures. Another aspect of this research was to analyze the temporal changes in the overall spatial-patterns of an epileptic intracranial EEG. Using Mantel statis tics, we showed that as seizure is approached, the overall spatial c onnectivity can be significantly different in comparison to the connectivity 2 hours before seizure. The exact time where the spatial bindings change can vary fr om seizure to seizure. 6-2 Future Research Directions In this dissertation, we explored the spa tio-temporal structure in an epileptic brain from a signal processing perspective. Unlike in many other clustering approaches where dynamical features extracted from the data are used as basis to determine groupings, our proposed clustering approach uses the depende ncies among the original data recordings to do the same. The SOM-based similarity index measure was proposed as a dependency metric to quantify couplings among different recordings and subsequently chosen to extract channel clusters. Comparisons of SOM-SI with other li near and non linear measures have already been articulated elsewhere [ 46 ]. However, in the context of clustering, it would be intere sting to check if any other de pendency measure would also produce same clusters as SOM-SI. Perhaps, the sensitivity of SOM-SI to non linear interactions can be confirmed by assessment of cluster groupings from linear measures. So far, because of the data size, we were constrained to analyze only on 12 seizures from 2 patients. Future effort in this direc tion would be to apply the proposed approach

PAGE 128

115 on a larger set of seizures and more patient s. In addition, since we analyzed only on complex partial seizures, it would be worthwh ile to check the cluster grouping in other types of seizures such as partial second ary generalized and sub-clinical seizures. Recall from the results that certain ch annels were always grouped together regardless of the seizure states. This raises a question if this pattern is unique to an epileptic patient and therefore be considered as a blueprint of seizures. One plausible way to answer this speculation would be to a pply the proposed clustering approach on normal subjects and then compare the differences in groupings with that of seizure patients As we have mentioned a few times earlier, one of our main objectives in this research was to develop engi neering tools to determine sp atio-temporal groupings in a multivariate epileptic brain. We proposed a si milarity-based clustering approach and used it to extract hidden structures from an epilep tic brain. One of the obvious limitations with any clustering approach is determining the optimal number of clusters. Techniques to combat the cluster size problems have been researched upon, without much success. In eigen vector based methods such as spectral clustering, cluster size can possibly be approximated to be equal to the number of eigen vectors corre sponding to significant eigen values. In multiple data sets however, the optimal cluster size need not have to be the same across different data sets rendering cluster comparisons weak. In our approach, we analyzed a large number of data sets and empirically, fixed the cluster size to 3. This may not be an efficient or a systematic appr oach to tackle the problem. Theoretic efforts are needed to develop mathematical criterion that allows us to determine a fixed cluster size, suitable to all groups of data. Beside s, exploring tools better than clustering to unravel hidden patterns in multidimensional time-sequences would be very beneficial.

PAGE 129

116 The Mantel statistics to track temporal ch anges in spatial patterns is a reasonably useful approach to distinguish pre-seizure pa tterns from those at inter-ictal stage of a seizure. So far, we have analyzed on a sma ll set of seizures and the results have looked encouraging. Assuming distinct pre-seizure patt erns exist; Mantel statistics analysis on a much larger set of seizures associated with different onset circumstances may help us in finding the general sensitivity of this technique.

PAGE 130

119 LIST OF REFERENCES 1. J. Arnhold, P. Grassberger, K. Lehnertz, and C. E. Elger, “A Robust Method for Detecting Interdependencies: Applica tion to Intracranially Recorded EEG,” Physica D vol. 134, pp. 419-430, 1999 2. L. A. Baccala, and K.Sameshima, “Overcoming the Limitations of Correlation Analysis for Many Simultaneously Processed Neural Structures,” Brain Research vol. 130, pp.33-47, 2001. 3. L. A. Baccala, and K. Sameshima, “Partial Directed Coherence: A New Concept in Neural Structure Determination,” Biological Cybernetics vol. 84, pp.463-474, 2001. 4. I. Baruchi, and E. Ben-Jacob, “Functi onal Holography of Recorded Neuronal Networks Activity,” Neuroinformatics vol. 2, pp. 333-351, 2004. 5. I. Baruchi, V. Towle, and E. Ben-J acob, “Functional Holography of Complex Networks Activity – From Cultures to the Human Brain,” Complexity vol. 10 (3), pp. 38 -51, Wiley Periodicals, 2005. 6. J. S. Bendat, and A. G. Piersol, Random Data: Analysis and Measurement Procedures, 2nd Ed. New York, John Wiley, 1986. 7. J. Bhattacharya, H. Petsche, and E. Pere da, “Interdependencies in the Spontaneous EEG while Listening to Music,” International Journal of Psychophysiology vol. 42, pp. 287-301, 2001. 8. J. Bhattacharya, E. Pereda, and H. Pets che, “Effective Detection of Coupling in Short and Noisy Bivariate Data,” IEEE Transactions on Systems, Man and Cybernetics Part B vol. 33 (1), pp. 85-95, 2003. 9. A. M. Bianchi, F. Panzica, F.Tinella, S. Franceschetti, S.Cerutti, and G. Baselli, “Analysis of Multichannel EEG Synchroni zation before and during Generalized Epileptic Seizures,” Proceedings of the 1st International IEEE EMBS Conference on Neural Engineering pp. 39-42, Capriland, 2003. 10. K. Blinowska, R. Kus, and M. Kaminski, “Granger Causality and Information flow in Multivariate processes,” Physical Review E vol. 70: 050902, pp: 01-04, 2004. 11. G. C. Carter, “Coherence and Time Delay Estimation,” Proceedings IEEE vol. 75(2), pp. 236-255, 1987.

PAGE 131

121 12. T. M. Cover, and J. A.Thomas, Elements of Information Theory New York, Wiley, 1991. 13. O. David, D. Cosmelli, and K. J.Friston, “Evaluation of Different Measures of Functional Connectivity Using a Neural Mass Model,” NeuroImage vol. 21, pp. 659-673, 2004. 14. O. David, D. Cosmelli, J. -P.Lachaux, S. Baillet, L. Garnero and J. Martinerie, “A Theoretical and Experimental Introduction to the Non-Invasive Study of the Large Scale Neural Phase Synchroni zation in Human Beings,” Invited Paper, International Journal of Computational Cognition vol. 1(4), pp. 53-77, December, 2003. 15. P. Dayan, and L. F. Abbot, Theoretical Neuroscience MIT Press, Boston, MA, 2001. 16. R. B. Duckrow, and S. S. Spencer, “Regi onal Coherence and the Transfer of Ictal Activity During Seizure Onset in the Medial Temporal Lobe,” Electroenceph Clin Neurophysiol vol. 82, pp. 415-422, 1992. 17. J. P. Eckmann, and K. D. Ruelle, “Recurrence Plots of Dynamical Systems,” Europhysics Letters, vol. 5, pp. 973-977, 1987. 18. M. -J. Fortin, M. R. T.Dale, and J, V. Hoef, “Spatial Analysis in Ecology,” Encyclopedia of Environmetrics Ed. By El-Shaarawi, A.H. and Piegorsch, W.W, vol. 4, pp. 2052-2058, 2002. 19. P. J. Franaszczuk, and G. K. Bergey, “Application of the Directed Transfer Function Method to Mesial and Latera l Onset Temporal Lobe Seizures,” Brain Topography vol. 11, pp. 13-21, 1998. 20. P. J. Franaszczuk, G. K. Bergey, and M. Kaminski, “Analysis of Mesial Temporal Seizure Onset and Propagation Using the Directed Transfer Function Method,” Electroenceph. clin. Neurophysiology vol. 91(6), pp. 413-427, 1994. 21. J. W. Freeman, “Characteristics of the Synchronization of Brain Activity Imposed by Finite Conduction Velocities of Axons,” International Journal of Bifurcation and Chaos vol. 10, pp. 2307-2322, 2000. 22. J. W. Freeman, and J. M. Barrie, “Ana lysis of Spatial Patterns of Phase in Neocortical Gamma EEGs in Rabbit,” Journal of Neurophysiology, vol. 84, pp. 1266-1278, 2000. 23. J. K. Friston, “Brain Function, Non-Li near Coupling and Neural Transients,” The NeuroScientist vol. 7 (5), pp. 406-418, 2001. 24. J. Gao, and H. Cai, “On the Structures and Quantification of Recurrence Plots,” Physics Letters A vol. 270, pp. 75-87, 2000.

PAGE 132

122 25. J. A. Garcia, J. Fdez-Valdivia, F. J. Cort ijo, and R. Molina, “A Dynamic Approach for Clustering Data,” Signal Processing, vol. 44(2), pp. 181-196, 1995. 26. J. Gotman, J. Ives, P. Gloor, A. Olivier, and L. F. Quesney, “Changes in Inter-Ictal EEG Spiking and Seizure Occurrence in Humans,” Epilepsia vol. 23, pp. 432-433, 1982. 27. C. W. J. Granger, “Investigating Caus al Relations by Econometric Models and Cross-spectral Methods,” Econometrica vol. 37(3), pp. 424-438, 1969. 28. M. A. F. Harrison, I. Osorio, M. G. Frei, S. Asuri, and Y-C. Lai, Y-C. “Correlation Dimension and Integral do not Predict Epileptic Seizures,” Chaos, vol. 15:033106, pp. 01-15, 2005. 29. A. Hegde, D. Erdogmus, Y. Rao, J. C.Principe, and J. B.Gao, “SOM-Based Similarity Index Measure: Quantifyi ng Interactions between Multivariate Structures,” Proceedings of NNSP’03, pp. 819-828, Toulouse, France Sep 2003. 30. A. Hegde, D. Erdogmus and J. C. Principe “Synchronization Analysis of Epileptic ECoG Data Using SOM-Based SI Measure,” Proceedings of EMBS 2004, pp. 952955, San Francisco, September 2004. 31. A. Hegde, D. Erdogmus, and J. C. Pr incipe, “Quantifying Spatio-Temporal Dependencies in Epileptic ECOG,” IEEE Signal Processing Sp ecial Issue Neural Co-ordination in the Brain ; A Signal Processing Perspective, vol. 85, pp. 20822100, 2005. 32. A. Hegde, D. Erdogmus, and J. C. Prin cipe, “Spatio-Temporal Clustering in Epileptic ECOG,” Proceedings of EMBS, Shanghai, September, 2005. 33. D. Hoyer, O. Hoyer, and U. Zweiner, “A New Method to Uncover Dynamic Phase Coordination and Synchronization,” IEEE Transactions on Biomedical Engineering, vol. 47 (1), pp. 68-74, 2000. 34. N. E. Huang, Z. Shen, S. R. Long, M. L. Wu, H. H. Shih, Q. Zheng, N. C. Yen, C. C. Tung, and H. H. Liu, “The Empirical Mode Decomposition and Hilbert Transform for Non-Linear and Non-St ationary Time Series Analysis,” Proc. Roy. Soc. London A, vol. 454, pp. 908-995, 1998. 35. L. D. Iasemidis, L. D. Olson, J. C. Sackellares, and R. Savit, “Time Dependencies in the Occurrence of Epileptic Seizures: a NonLinear Approach,” Epilepsy Research vol. 17, pp. 81-94, 1994. 36. L. D. Iasemidis, J. C. Principe, J. M. C zaplewski, R. L. Gilmore, S. N. Roper, and J. C. Sackellares, Spatiotemporal Transition to Epileptic Seizures: a Nonlinear Dynamical Analysis of Scalp and Intracranial EEG Recordings In: Lopes da Silva F. H., Principe J. C., Almeida L. B., eds. S patiotemporal models in biological and artificial systems. Amsterdam: IOS Press, pp: 8189, 1997

PAGE 133

123 37. L. D. Iasemidis, J. C. Sackellares, H. P. Zaveri, and W. J. Williams,“ Phase Space Topography and the Lyapunov Exponent of Electrocorticograms in Partial Seizures,” Brain Topography vol. 2, pp. 187-201, 1990. 38. L. D. Iasemidis, K. E. Pappas, J. C. Principe, and J. C. Sackellares, “SpatioTemporal Dynamics of Human Epileptic Seizures,” Proceedings of the 3rd Experimental Chaos Conference eds. R.G. Harrison et al. World Scientific, Singapore, pp. 26-30, 1996. 39. L. D. Iasemidis, P. Pardalos, J. C. Sackel lares, and D. S. Shiau, “Quadratic Binary Programming and Dynamical System Approach to Determine the Predictability of Epileptic Seizures,” Journal of Combinatorial Optimization vol. 5, pp. 09-26, 2001. 40. L. D. Iasemidis, D. S. Shiau, W. Chaovalitwongse, J. C. Sackellares, P. M. Pardalos, J. C. Principe, P. R. Carney, A. Prasad, B. Veeramani, and K. Tsakalis, “Adaptive epileptic seizure prediction system,” IEEE Transactions on Biomedical Engineering, vol. 50 (5), pp. 616-627, 2003. 41. M. Kaminski, and K. J. Blinowska, “A New Method of the Description of the Information Flow in the Brain Structures,” Biological Cybernetics vol. 65, pp. 203-210, 1991. 42. E. R. Kandel, J. H. Schaertz, and T. M. Jessell, Principles of Neural Science New York, 4th Ed., McGraw-Hill, 2000. 43. H. Kantz, and T. Schreiber, Nonlinear Time Series Analysis Cambridge University Press, Cambridge, England, 1997. 44. E. Keogh, J. Lin, and W. Truppel, “Clust ering of Time Series Subsequences is Meaningless: Implications for Past and Future Research” Proceedings of the 3rd IEEE International Conference on Data Mining (ICDM 2003), pp.115-122, Melbourne, Florida, 2003. 45. M. A. Kramer, E. Edwards, M. Soltani, M. S. Berger, R. T. Knight, and A. J. Szeri, “Synchronization Measures of Burs ting Data: Application to the Electrocorticogram of an Auditory Event Related Experiment,” Physical Review E vol. 70:011914, pp. 01-10, 2004. 46. A. Kraskov, T. Kreuz, Q. Quiroga,P. Gras sberger, F. Mormann, K. Lehnertz, and C. E. Elger, “Phase Synchronization Us ing Continuous Wavelet Transform of the EEG for Inter-Ictal Focus Localization in Mesial Temporal Lobe Epilepsy,” Epilepsia, 42(7):43, pp. 01-04, 2001. 47. A. Kraskov, T. Kreuz, Q. Quiroga, P. Gras sberger, F. Mormann, K. Lehnertz, and C. E. Elger, “Comparison of Two Phase Synchronization Analysis Techniques for Inter-Ictal Focus Lateralization in Mesial Temporal Lobe Epilepsy,” Epilepsia, vol. 43(7):48, pp. 01-04, 2002.

PAGE 134

124 48. J. P. Lachaux, E. Rodriguez, J. Mart inerie, and F. Varela, “Measuring Phase Synchrony in Brain Signals,” Human Brain Mapping vol. 8, pp. 194-208, 1999. 49. A. Laird, J. Carew, and M. E. Meyera nd, “Analysis of the Instantaneous Phase Signal of a FMRI Time-Series Vi a The Hilbert Phase Transform,” Conference Record of the Asilomar Conference on Signals, Systems and Computers v 2, pp. 1677-1681, 2001. 50. F. –J. Lapointe, and P. Legendre, “A Sta tistical Framework to Test the Consensus Among Additive Trees (Cladograms),” Systematic Biology, vol. 41 (2), pp. 158171, 1992. 51. P. Legendre, F. -J. Lapointe, and P. Casgrain,“Modeling Brain Evolution from Behavior: A Permutational Regression Approach,” Evolution vol. 48 (5), pp. 14871499, 1994. 52. K. Lehnertz, R. G. Andrzejak, J. Arnhold, T. Kreuz F. Mormann, C. Rieke, G. Widman, and C. E. Elger, “Nonlinear EEG Analysis in Epilepsy: Its Possible Use for Interictal Focus Localization, Se izure Anticipation, and Prevention, ” Journal of Clinical Neurophysiology vol. 18 (3), pp. 209-222, 2001. 53. K. Lehnertz, F. Mormann, T. Kreuz, R. G. Andrzejak, J. Arnhold, C. Rieke, P. David, and C. E. Elger, “Seizure Prediction By Nonlinear EEG Analysis, ” IEEE Engineering in Medicine and Biology Magazine pp. 57-63, 2003. 54. K. Lehnertz, “Non-Linear Time Series Anal ysis of Intracranial EEG Recordings in Patients with Epilepsy an Overview,” International Journal of Psychophysiology vol. 34, pp. 45-52, 1999. 55. M. Li, J. H. Badger, X. Chen, S. Kwong, P. Kearney, and H. Zhang, “An Information-Based Sequence Distance and its Application to Whole Mitochondrial Genome Phylogeny,” Bioinformatics vol. 17 (2), pp.149-154, 2001. 56. F. H. Lopes da Silva, W. Blanes, S. N. Kalitzin, J. Parra, P. Suffczynski, and D. N. Velis, “Dynamical Diseases of Brain Sy stems: Different Routes to Epileptic Seizures (Invited Paper),” IEEE Transactions on BioMedical Engineering vol. 50, pp. 540-549, 2003. 57. J. Malik, S. Belongie, T. Leung, and J. Shi, “Contour and Texture Analysis for Image Segmentation,” Interna tional Journal of Computer Vision, 43(1), pp. 7-27, 2001. 58. N. Mantel, “The Detection of Disease Clustering and a Generalized Regression Approach,” Cancer Research, vol. 27, pp. 209-220, 1967. 59. N. Marwan, and J. Kurths, “Nonlinear An alysis of Bivariate Data with Cross Recurrence Plots,” Physics Letters A vol. 302, pp. 299-307, 2002.

PAGE 135

125 60. F. Mormann, R. G. Andrzejak, T. Kreuz, C. Rieke, P. David, C. E. Elger, and K. Lehnertz, “Automated Detection of a Pr eseizure State Based on a Decrease in Synchronization in Intracrani al Electroencephalogram Recordings from Epilepsy Patients,” Physical Review E vol. 67:02192, pp. 01-10, 2003. 61. F. Mormann, T. Kreuz, R. G. Andrzejak, P. David, K. Lehnertz and C. E. Elger, “Epileptic Seizures are Preceded by a Decrease in Synchronization,” Epilepsy Research vol. 53, pp. 173-185, 2003. 62. F. Mormann, T. Kreuz, C. Rieke, R. G. Andrzejak, A. Kraskov, P. David, C. E. Elger, and K. Lehnertz, “On the Pred ictability of Epileptic Seizures,” Clinical Neurophysiology vol. 116, pp. 569-587, 2005. 63. F. Mormann, F., K. Lehnertz, P. David, and C. E. Elger, “Mean Phase Coherence as a Measure for Phase Synchronization and its Application to the EEG of Epileptic Patients,” Physica D vol. 144, pp. 358-369, 2000. 64. A. Y. Ng, M. I. Jordan, and Y. Weiss, “On Spectral Clustering: Analysis and an Algorithm,” Advances in Neural Information Processing Systems 14, pp. 849-856, 2002. 65. A. Ossadtchi, S. Baillet, J. C. Mosher, and R. M. Leahy, “Using Mutual Information to Select Event Related Components in ICA,” Proceedings of International Conference on Biomagnetism Finland, 2000. 66. M. Palus, “Detecting Phase Synchronization in Noisy Systems,” Physics Letters A vol. 225, pp. 341-351, 1997. 67. M. Palus, “Kolmogorov Entropy from Time Series Using Information-Theoretic Functionals,” Neural Networks World vol. 3, pp. 269 292, 1997. 68. M. Palus, V. Komarek, T. Prochazka, Z. Hrncir, and K. Sterbova, “Synchronization and Information Flow in EEGs of Epileptic Patients,” IEEE Engineering in Medicine and Biology pp. 65-71, 2001. 69. L. M. Pecora, and T. L. Caroll, “Synchronization in Chaotic Systems,” Physics Review Letter vol. 64, pp. 821 830, 1990. 70. B. Pompe, “A Tool to Measure Dependencies in Data Sequences,” Advanced Mathematical Tools in Metr ology II, World Scientific Singapore, pp. 81-99, 1996. 71. B. Pompe, “Measuring Statistical Dependencies in a Time Series,” Journal of Statistical Physics vol. 73, pp. 587-610, 1993. 72. B. Pompe, P. Blidh, D. Hoyer, and M. Eiselt, “Using Mutual Information to Measure Coupling in the Ca rdiorespiratory System,” Proceedings of the Engineering in Medicine and Biology Magazine, IEEE vol. 17 (6), pp. 32-39, 1998.

PAGE 136

126 73. D. Prichard, and J. Theiler, “Generati ng Surrogate Data for Time Series with Several Simultaneously Measured Variables,” Physical Review Letters vol. 73 (7), pp. 951-954, 1994. 74. J. C. Principe, N. R. Euliano, and W. C. Lefebvre, Neural and Adaptive Systems: Fundamentals through Simulations John Wiley & Sons, 2000. 75. J. G. Proakis, Digital Communications 4th edition, McGraw Hill, 2000. 76. Q. R Quiroga, J. Arhnold, and P. Grassberger, “Learning Driver-Response Relationships from Synchronization Patterns,” Physical Review E vol. 61, pp. 5142-5148, 2000. 77. Q. R. Quiroga, J. Arnold, K. Lehnertz, and P. Grassberger, “Kulback-Leibler and Renormalized Entropies: Applications to Electroencephalograms of Epilepsy Patients,” Physical Review E vol. 62 (6), pp. 8380 – 8386, 2000. 78. Q. R Quiroga, A. Kraskov, and T. Kreuz, “Performance of different synchronization measures in real data: a case study on EEG signals” Physica Review E vol. 65:041903, pp. 01 – 14, 2002. 79. Q. R Quiroga, T. Kreuz, and P. Grassber ger, “Event Synchronization: A Simple and Fast Method to Measure Synchr onicity and Time Delay Patterns,” Physica Review E vol. 66:041904, pp: 01-06, 2002. 80. Q. R. Quiroga, O. W. Sackowitz, E. Basar, and M. Schurmann, “Wavelet Transform in the Analysis of the freq uency Composition of Evoked Potentials,” Brain Research Protocols vol. 8, pp. 16-24, 2001. 81. M. L. V. Quyen, C. Adam, M. Baulac, J. Martinerie, and F.J. Varela, “Nonlinear Interdependencies of EEG Signals in Hu man Intracranially Recorded Temporal Lobe Seizures,” Brain Research vol. 792, pp. 2440, 1998. 82. M. L. V. Quyen, J. Martinerie, V. Navarro, P. Boon, M. D. Have, C. Adam, B. Renault, F. Varela, and M. Baulac, “Anticipation of Epileptic Seizures from Standard EEG Recordings,” The Lancet vol. 357, pp. 183-188, 2001. 83. M. L. V. Quyen, J. Foucher, J. Lachaux, E. Rodriguez, A. Lutz, M. Francisco,and J. Varela, “Comparison of Hilbert Transf orm and Wavelet Methods for Analysis of Neuronal Synchrony,” Journal of Neuroscience Methods vol. 111, pp. 83-98, 2001. 84. E. Rodriguez, N. George, J. –P. Lachaux, J. Martinerie, B. Renault, and F.J Varela, “Perception’s Shadow: Long-Distance Sync hronization of Human Brain Activity,” Nature vol. 397, pp. 430-437, 1999.

PAGE 137

127 85. M. Rosenblum, A. Pikovsky, J. Kurths, C. Schaefer, and P. A. Tass, Phase Synchronization: from Theory to Data Analysis. In F. Moss and S. Gielen, editors, Handbook of Biological Physics, page 297, Elsevier Science Foundation, Amsterdam, 2001. 86. M. Rosenblum, A. Pikovsky, and J. Kurt hs, “Phase Synchronization of Chaotic Oscillators,” Physics Review Letter vol. 76, pp. 1804-807, 1996. 87. M. G. Rosenblum, A. S. Pikovsky, a nd J. Kurths “From Phase to Lag Synchronization in Couple d Chaotic Oscillators,” Physics Review Letter vol. 76, pp. 1804-1807, 1996. 88. N. F. Rulkov, M. M. Sushchik, L. S. Tsimring, and H. D. I. Abarbanel, “ Generalized Synchronization of Chaos in Directionally Coupled Chaotic Systems”, Physics Review E vol. 51(2), pp. 980-994, 1995. 89. S. Haykin, Neural Networks: A Comprehensive Foundation 2nd edition, Prentice Hall, 1999. 90. V. J. Samar, A. Bopardikar, R. Rao, and K. Swartz, “Wavelet Analysis of Neuroelectric Waveforms: A Conceptual Tutorial,” Brain and Language vol. 66, pp. 07-60, 1999. 91. K. Sameshima, and L. A. Baccala, “Using Partial Directed Coherence to Describe Neuronal Ensemble Interactions,” Journal of Neuroscience Methods vol. 94, pp. 93-103, 1999. 92. S. J. Schiff, P. So, T. Chang, R. E. Burke, and T. Sauer, “Detecting Dynamical Interdependence and Generalized Synchrony through Mutual Prediction in a Neural Ensemble,” Physics Review E vol. 54, pp. 6708, 1996. 93. A. Schmitz, “Measuring Statistical Depe ndence and Coupling of Subsystems,” Physical Review E vol. 62, pp: 7508-7514, 2000. 94. T. Schreiber, “Measuring Information Transfer,” Physics Review Letters vol. 85, pp. 461-470, 2000. 95. T. Schreiber, and A. Schmitz, “Improved Surrogate Data for Nonlinearity Tests,” Physical Review Letters vol.77 (4), pp. 635-638, 1996. 96. W. D. Shannon, M. A. Watson, A. Perry, and K. Rich, “Mantel Statistics to Correlate Gene Expression Levels from Mi croarrays with Clinical Covariates,” Genetic Epidemiology vol. 23, pp. 87-96, 2002. 97. C. J. Stam, M. Breakspear, A. M. C. Walsum, and B. W. Dijk, “Nonlinear Synchronization in EEG and Whole-head MEG recordings of Healthy Subjects,” Human Brain Mapping vol. 19 (2), pp. 63-78, 2003.

PAGE 138

128 98. F. Takens, Detecting Strange Attractors in Turbulence In: Rand DA, Young LS, eds. Dynamical systems and turbulence. Lecture notes on mathematics. vol. 898, Berlin: Springer, pp. 366-381, 1981. 99. P. Tass, M. G. Rosenblum, J. Weule, J. Kurths, A. Pikovsky, and J. Volkmann, “Detection of n:m Phase Locking fr om Noisy Data: Application to Magnetoencephalography,” Physics Review Letter vol. 81, pp. 3291-3294, 1998. 100. Y. Weiss, “Segmentation Using EigenVectors: A Unifying View,” International Conference on Computer Vision pp: 975-982, 1999. 101. C. S. Burrus, R. A. Gopinath, and H. Guo, Introduction to Wavelets and Wavelets Transforms A Primer, Prentice-Hall Interna tional Editions, Engions, Clifs, NJ, 1998.

PAGE 139

129 BIOGRAPHICAL SKETCH Anant graduated from Mysore University, India, with a bachelors of Engineering degree in Electronics and Co mmunication engineering (1998). For a brief period, he worked in the personal communications sector at the Bangalore design center in Motorola India Ltd., Later, during master’s research, he worked on evoked potentials at the BioSignal Analysis Laboratory at the University of Houston. Currently he is pursuing his doctorate degree in electrical e ngineering at the University of Florida. His research at the Computational Neuroengineering Laboratory involves investigating the functional connectivities between brain areas at various states of epile ptic seizures. His research interests are in biomedical signal pro cessing, broad band communication technologies and machine learning networks adapted from in formation theoretic le arning principles.


Permanent Link: http://ufdc.ufl.edu/UFE0013522/00001

Material Information

Title: Spatio-Temporal Dependency Analysis of Epileptic Intracranial Electroencephalograph
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0013522:00001

Permanent Link: http://ufdc.ufl.edu/UFE0013522/00001

Material Information

Title: Spatio-Temporal Dependency Analysis of Epileptic Intracranial Electroencephalograph
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0013522:00001


This item has the following downloads:


Full Text












SPATIO-TEMPORAL DEPENDENCY ANALYSIS OF EPILEPTIC INTRACRANIAL
ELECTROENCEPHALOGRAPH















By

ANANT HEGDE


A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA


2006

































Copyright 2006

by

ANANT HEGDE
































This dissertation is dedicated to my family, teachers and friends for their enduring love,
support and friendship.















ACKNOWLEDGMENTS

It has been a privilege being Dr. Jose Principe's student. His constant guidance,

encouragement and technical directions helped me tremendously in approaching research

challenges. The fact that he has infinite time for his students makes him very special in

my eyes. No wonder his enthusiasm, his quest for excellence and his abilities to inspire

and motivate have taken him a long way. My experience as a researcher under his

supervision has been very rich and fulfilling and I would like to thank him for teaching

me how to research, for giving me enough freedom to explore problems independently

and also for his valuable advice whenever I needed it most.

I would like to thank Dr. Chris Sackellares for being on my committee and offering

me plentiful guidance in research. It is always a privilege to know someone like Dr. John

Harris. I thank him for his support, both moral and technical, and also for serving in my

committee. I also thank Dr. Mingzhou Ding for being on my committee and his valuable

feedback in my research.

It is seldom that someone working in CNEL would not know or have heard of Dr.

Deniz Erdogmus. I consider myself lucky having been around when Deniz was here.

Besides being a great friend, he offered relentless guidance all along my research and his

contributions in this accomplishment are as precious as anyone else's. I would also like

thank Dr. Deng Shiau and Dr. Yadunandana Rao for their technical feedback.

I also wish to extend my acknowledgements to all the members of CNEL who have

been primarily responsible for my fruitful stay in the lab. In particular, I would like to









cherish those productive discussions with Mustafa Can, Rui Yan, Jianwu Xu, Puskal,

Sudhir, Han, Kyu-hwa and Jack. Overall, my stay in CNEL has been very rewarding and

if I look back a few years from now, I am sure I'll say 'No doubt, those were my

formative years.' I would like to extend my gratitude to the always cheerful Ellie

Goodwin for her golden words of wisdom. Her ability to get things done was truly

remarkable. I would also like to acknowledge Linda Kahila and Shannon for their

extensive support and assistance during my stay at UF. Kartikeya Tripati has been very

inspiring all along and I am very grateful to his constant support.

I would like to thank my family and friends for their constant love and

encouragement. They have allowed me to pursue whatever I wanted in life. Without their

guidance and affection, it would have been impossible for me to be where I am today.

Lastly, I would like to thank my life partner Madhoo for making my life so colorful

and for being on my side whenever I need her. Her everlasting love has made me a better

individual.















TABLE OF CONTENTS

page


A C K N O W L E D G M E N T S ................................................................................................. iv

LIST OF TABLES ................................................................ ........ viii

LIST OF FIGURES ......... ......................... ...... ........ ............ ix

ABSTRACT .............. ..................... .......... .............. xii

CHAPTER

1 SPATIO-TEMPORAL INTERACTIONS ................................ ............... 1

1.1 Introdu action ..............................1..........................
1.2 Different Kinds of Synchronization.................. ............... ..... ...............3
1.3 Second Order Synchronization M easures.................................... .....................4
1.3.1 C ross-C correlation ................................ .... .......... .............. .. ......5
1.3.2 Partial D directed Coherence (PD C) ........................................ ....................6
1.3.3 Synthetic sim ulation ......................................... ...................... ...... .. .. .. 8
1.4 M mutual Inform action ................................................ ........ .. ............ 11
1.5 Phase Synchronization.................... .......... ............................. 12
1.5.1 Hilbert Transformation...................... ........... ................... 13
1.5.2 Empirical Mode Decomposition (EMD)..............................................15
1.5.3 W avelet decomposition ......... .. ...................................... .......... 18
1.5.4 Indexing Phase Synchronization ..................................... ............... 19
1.6 Non-Linear synchronization measures ............. ............................. ....... ........... 22
1.6.1 Signal Reconstruction.................. ......... .......................... 22
1.6.2 Non-Linear Synchronization Measures ........... ...................................23
1.7 Objectives and Author's Contribution ............... .......................................... 25

2 SELF-ORGANIZING-MAP BASED SIMILARITY INDEX MEASURE ..............27

2 .1 Introduction ................... ............................................... ................. 27
2.2 The Similarity Index Technique .............................................. ..... .......... 27
2.3 SOM-based Similarity Index (SOM-SI)............................................................28
2.4 Simulation Results .................. ........................................ .... ........ 31
2.4.1 R osseler-Lorenz sim ulation.................................... ........................ 31
2.4.2 Intracranial EE G Sim ulation ........................................... .....................35
2.5 Epileptic Intracranial EEG Data Description ....................................... .......... 38









2.6 Testing the application of SOM-SI on Epileptic Intracranial EEG...................38
2.7 Statistical comparison between SI and SOM-SI............................... .............40
2.8 Testing the Robustness of SOM-SI on Multiple SOMs ....................................42
2 .9 S u m m ary ............. ............ .. ................................................... 4 6

3 SPATIO-TEMPORAL CLUSTERING MODEL ............................................... 49

3 .1 Intro du action .................. ...................................... ......... ................ 4 9
3.2 M odel for Spatio-temporal clustering........................................ ............... 50
3.3 Temporal Quantification Using Hamming Distance ........................................53
3.4 Sim ulations w ith Synthetic D ata ........................................ ....... ............... 55
3.4.1 R oessler- Lorenz System ........................................ ....... ............... 55
3.4.2 L inear M odel ...................................... .............................58
3 .5 S u m m ary ................................4...........................

4 APPLICATION OF SPATIO-TEMPORAL MODEL TO EPILEPTIC
IN TR A C R A N IA L E E G ................................................................... .....................65

4.1 Introduction ................ ....................... .................. .................. .....65
4.2 Application of SOM-SI on Epileptic Intracranial EEG Data ............................65
4 .3 R e su lts................................ ........................................................ ............... 6 9
4 .4 Statistical V alidation ......... ........................................................ .. .... ...... 75
4 .4 S u m m a ry .....................................................................................................8 9

5 STATISTICAL ASSESSMENT OF SPATIO-TEMPORAL DEPENDENCY
CHANGES IN EPILEPTIC INTRACRANIAL EEG...........................................91

5.1 Introduction ..................................................................... 91
5.2 Spatio-tem poral Changes...................... .... ................................. ............... 91
5.3 Mantel Test of Matrix Comparison .................. ..........................92
5.4 Application to Epileptic Intracranial EEG data......... ............ ............... 94
5.5 Statistical A pproach................................................... .. ...... .......... .. ..96
5.6 Inter-Hemisphere Synchronization Differences ....................................... 98
5 .7 Statistical A ssessm ent............................................ ....................................... 104
5.8 E xperim ental R results ........................................................................... ...... 105
5.8.1 Experim ent #1 ...................... ................ ................... .. .... .. 105
5.8.2 Experiment #2 ............. ...... .............................. .......... 107
5 -9 S u m m ary ...................................................... ................ 10 8

6 CONCLUSIONS AND FUTURE DIRECTIONS .............................112

6 -1 C o n clu sio n s ...................................... ......... ............... ................ 1 12
6-2 Future R research D directions ......................................................... .................114

L IST O F R EFER EN CE S ......... ............................................................ ........... ........ .. 119

BIOGRAPH ICAL SKETCH .............. ......................... ................... ............... 129















LIST OF TABLES


Table p

2-1 Coupling strength between pairs of signals, the normalized similarity index and
the original similarity index between them........... ................. ..... ............37

2-2 Quantitative comparisons between the SOM-SI profiles obtained from SOM-1
and SOM-2. LOF3 & LOF4 data was projected on each of the SOMs and then
the SOM-SI measure was applied to analyze the dependency of LOF3 on LOF4. .46

2-3 Summary of the comparisons between the SOM-SI profiles from SOM-1 and
SOM-2. Each row represents the statistics (mean and variance) of pair-wise
SOM-SI analyses of the epileptic ECOG data from 6 channels (15
com binations). .........................................................................48

3-1 Cluster groupings for the simulation example ............................................. ........... 58

4-1 Spatio-temporal groupings as obtained for seizure 11 of patient P093 ..................74

4-2 P093, Seizure 11: Over each 30 minute (91 samples total) window, number of
times the within-cluster interaction is greater than between-cluster interaction, at
95% significance level. .................................................................... ..................77

4-3 Spatio-temporal groupings as obtained for seizures 4 and 5 of patient P093. ........80

4-4 Spatio-temporal groupings as obtained for seizure 6 and 7 of patient P093............80

4-5 Spatio-temporal groupings as obtained for seizure 1 of Patient P092. ..................86

4-6 Spatio-temporal groupings as obtained for seizure 3 of Patient P092. ....................86

4-7 Spatio-temporal groupings as obtained for seizure 4 of Patient P092. ....................87

4-8 Spatio-temporal groupings as obtained for seizure 5 of Patient P092. ....................87

4-9 Spatio-temporal groupings as obtained for seizure 6 of Patient P092. ....................88

4-10 Spatio-temporal groupings as obtained for seizure 7 of Patient P092. ....................88















LIST OF FIGURES


Figure p

1-1 Coupling diagram showing the direction of interactions between nodes, xl, x2
and x 3 ............................................................................................................ 8

1-2 Matrix layout plots for PDC, describing the interactions in example simulation......9

1-3 Synthetic Sinusoidal signal X(t) and its components xl, x2 and x3. X(t) also
consists of noise w -N (0, 0.1) ................................................................................... 17

1-4 Decomposition of sinusoids using EM D .................................... ......... ............... 17

1-5 Hilbert spectrum constructed from the IMF functions .........................................18

1-6 Decomposition of sinusoids using wavelets ......................................................19

1-7 Hilbert spectrum constructed from the wavelet reconstructed narrow-band
sig n a ls ............. ......... .. .. ......... .. .. .................................................. 2 0

2-1 Phase-space trajectories of the Rosseler-Lorenz system for various coupling
strengths. ............................................................................33

2-2 Illustrating the dependency relationships between the Rosseler (driver) and the
Lorenz system (response) as a function of coupling strength................................34

2-3 Original EEG signals X and Y. for the EEG simulation example .........................34

2-4 The nonlinearly mixed (synthetic) EEG signals V, Z, and W...............................34

2-5 Schematic diagram of the coupling function between the signals........................35

2-6 The phase-space trajectory of the training EEG signal (projected in two
dimensions) and the weights of the trained EEG-SOM (circles). ............................36

2-7 Diagram of the depth and subdural electrode montage in an epileptic brain...........39

2-8 The phase-space trajectory of the training intracranial EEG signal (solid lines)
superimposed on the weights (dots) of the trained 25x25 EEG-SOM grid. ...........40

2-9 A qualitative illustration of the accuracy of the 25x25 EEG-SOM.. ......................41















2-10 SOM-SI results showing the dependencies between two signals. X and Y
correspond to channels RTD4 and RST1, respectively.......................................42

2-11 Experiment setup to compare SOM-Similarity Indices obtained from 2 separate
m ap s. ............................................................................... 44

2-12 Time intervals from all the seizures used as test data (not drawn to scale) ............44

2-13 Comparing interdependencies between channels LOF3 and LOF4.......................47

3-1 Block diagram to extract spatio-temporal groupings information in Multivariate
E E G stru ctu res............................. .................................................... ............... 50

3-2 Time profiles of the Roessler and the Lorenz components. a) Coupling strength,
C = 0 and b) C = 5 .....................................................................57

3-3 Dendrogram illustration of similarity-based time series clustering, a) Coupling
strength, C = 0. .........................................................................60

3-4 Model showing the generation of linearly dependent signals using configuration
1..... ......... .................................. ........ 62

3-5 Plot showing the two band pass filters, separated by a narrow pass band ..............63

3-6 Z-planes showing the pole- zero configurations for the two band-pass filters in
F ig 3-5 ...............................................................................63

4-1 Maximum average driving ability of each of the six (6) channels, nearly 100
minutes before and 70 minutes after Seizure-1 in patient P092 .............. ..............68

4-2 Seizure 11 of patient P093: Cluster-similarity matrices P indicating the
probability that two channels share the same cluster label in a 30 minute time-
interv al ............................ ...................... ..... ............. ................. 72

4-3 Dendrogram representation of the cluster results in Seizure 11, P093 ...................73

4-4 Statistical validation of the clustering results. In each panel, thick lines are used
to represent the profiles of the three clusters in a 30 minute time interval ............78

4-5 Seizures 4 and 5 of patient P093: Cluster-similarity matrices indicating the
probability that two channels share the same cluster label in a 30 minute time-
interv al ................. .................................... ........................... 8 1















4-6 Seizures 6 and 7 of patient P093: Cluster-similarity matrices indicating the
probability that two channels share the same cluster label in a 30 minute time-
interval .......... .............. ............................................ ................ 82

4-7 Dendrograms corresponding to P092, Seizure 1.............................................83

4-8 Cluster-similarity matrices indicating the probability that two channels share the
same cluster label in a 30 minute time-interval................................. ............... 85

5-1 Quantifying the spatial changes along time using Mantel statistics.........................99

5-2 Time smoothened profiles of the Mantel statistics during inter-ictal states...........100

5-3 Illustrating the statistical difference between z-score at time't' and z-score at a
reference time, 90, 100 and 110 minutes before a seizure. Top to Bottom:
profiles of z-score statistics for seizures 1 through 3 of patient P092.................. 101

5-4 Illustrating the statistical difference between z-score at time't' and z-score at a
reference time, 90, 100 and 110 minutes before a seizure. Top to Bottom:
profiles of z-score statistics for seizures 4 through 6 of patient P092.................. 102

5-5 Illustrating the statistical difference between z-score at time't' and z-score at a
reference time, 90, 100 and 110 minutes before a seizure. Top to Bottom:
profiles of z-score statistics for seizure 7 of patient P092.................................... 103

5-6 Results of synchronization analysis on seizure 1 of P092, to check the
differences between the focal and non-focal hemispheric interactions, at pre and
post-ictal states.. ..................................................................... 109

5-7 Null-hypothesis rejections rate averaged over 10 seizures for experiment 1.........110

5-8 Null-hypothesis rejections rate averaged over 10 seizures for experiment 2.........111















Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy

SPATIO-TEMPORAL DEPENDENCY ANALYSIS OF EPILEPTIC INTRACRANIAL
ELECTROENCEPHALOGRAPH

By

Anant Hegde

May 2006

Chair: Jose C. Principe
Major Department: Electrical and Computer Engineering

There is evidence that the mechanisms leading to epileptic seizures can be

understood by continuously tracking the ongoing spatio-temporal mappings in the brain.

However, the exact spatio-temporal changes leading to epileptic seizures, although

widely studied, are not yet well understood. Previous studies have mostly focused on

individually tracking the dynamical changes along time. However, approaches to

simultaneously explain a system's temporal changes in its overall spatial configuration

can be much more effective in efforts to characterize epileptic events.

In this dissertation, we propose a simple statistical approach to quantify the

temporal changes in spatial patterns of an intracranial EEG. Previously, we developed a

non-linear synchronization measure, called the SOM-Similarity Index, to quantify mutual

associations between various brain regions. We propose to apply the Mantel test statistics

on the SOM-similarity indices to track the temporal changes of the spatial patterns.









Statistical comparisons between inter-ictal and pre-ictal states suggest significant changes

in the spatial connectivity prior to a seizure.

Another aspect of this study investigates the regional groupings in the intracranial

EEG spatial networks at different ictal-states. We develop a model that will allow us to

study the various cluster patterns in an epileptic brain. Studies on 10 seizures from 2

patients reveal strong connections between the left-sub-temporal and the left-temporal

depth areas. In addition, strong homologous connectivity is found in the orbito-frontal

regions.

In the third aspect of this study, we investigate the differences in the

synchronization levels between the focal and the non-focal hemispheres at pre-ictal and

post-ictal states. Statistical tests confirm the existence of significant differences at pre-

ictal states followed by a strong entrainment, 20 minute post-ictal period.














CHAPTER 1
SPATIO-TEMPORAL INTERACTIONS

1.1 Introduction

The word "interaction" in a system can possibly be referred to as an exchange of

information between its components. A system can be understood to be made up of a

large number of independent components, each acting as an individual sub-system. These

sub-systems generally possess an independent environment of their own. However, in a

synergetic environment, they also communicate or share information with their

counterparts in such a way that the cumulative effect is responsible for the overall

functioning of the system.

Even though the notion of a system is the same in any discipline, the specific

interpretation or representation of the same differs and is usually tailored to suit our

understanding. In physics or engineering, a system is usually considered to be made up of

a set of well connected oscillators that are assumed to have independent rhythms

specified over periods of time and expressed as time-sequences. Harmonics induced by

the oscillator rhythms encode the state of the oscillator and perhaps the system as well.

As a result of coupling with other oscillators in an interactive environment, a signal

recorded from a transducer can be expressed in a number of ways: a) as a superposition

of oscillations generated individually with oscillations shared from other oscillators or b)

as a complex interactive process induced by phase, frequency or amplitude variations

between oscillators in the spatial neighborhood. Thus the process where oscillators share









complex dynamics with their counterparts can be characterized as a state of mutual-

harmony or synchronization. Referring to the degree of interaction, it is possible that it be

weak or strong, but is nevertheless capable of triggering an event. While a lot of

enthusiasm is placed on understanding the exact nature of synchronization, equal

importance has been devoted to relate the event manifestations caused due to

synchronization changes. Finally, synchronization analysis of the time-recordings of a

physical system allows us to probe the underlying hidden dynamics involved in creating

the system.

Synchronization changes in a system can be across both space and time.

Particularly, in a multi-variate system, understanding the interactions among its various

nodes, whose behavior can be represented along time as time-sequences, presents

numerous challenges. One of the key aspects of tightly coupled systems with spatial

extent is their ability to interact both across space and time, which complicates the

analysis greatly. In biological systems such as the central nervous system, this difficulty

is compounded by the fact that the components of interest have nonlinear complicated

dynamics that can dictate overall changes in the system behavior [22]. In summary, it is

important to realize that the spatial structures also undergo functional changes that can

dictate the overall behavior of a system.

With the advent of sophisticated tools, multivariate time-signal analysis has

become the forefront of research in discrete time signal processing [66-72, 85-87]. Time-

sequences recorded from real-world systems encode huge amounts of information that

are both temporally and spatially extended. Especially in neural systems, synchronization

between cortical structures or local-field potentials recorded at the tissue level is believed









to be one of the prominent indicators of physiological and pathological events [21-22]. In

the next section therefore, we provide a brief literature overview of the efforts in this area

and then in the last part of this chapter, we try to address the fundamental question related

to epileptic seizures that eventually served as a motivation for us to pursue this study.

1.2 Different Kinds of Synchronization

In the previous section, we introduced the notion of coupling and abstractly

described it as a process of information exchange between components of a system.

Exactly how information is shared is a question that remains ambiguous. Previous studies

on multi-variate time series analyses have resulted in development of a wide array of

signal-processing tools for quantification of coupling in systems. However, the general

consensus on the best approach to quantify this phenomenon is largely uncertain.

Synchronization can be conveniently classified as linear or non linear [13]. Linear

synchronization assumes that the processes are Gaussian distributed, stationary, ergodic

and therefore the couplings can be expressed using just the first and second order

moments. While some of the tools quantify the correlations by exploiting the time

parameters such as amplitude or phase of the data, the others operate in the frequency

domain to achieve the same [11, 60-63]. Recently, multi-resolution techniques have

facilitated significant progress in exploiting the time and frequency parameters jointly

[90].

One of the major advantages of assuming linearity is the fact that the system

interactions can be quantitatively analyzed in great depth. However, assuming a linear

model can be arguably very weak and in practice, it is very restrictive to assume that the

functional mapping between two coupled systems is sensitive to the differences in first

and second order moments. Nonlinear coupling can produce changes in the higher order









moments of stochastic signals without affecting appreciably the lower order moments,

which raises concerns on the general applicability of these methods.

Most real world systems are evidently non-linear processes and their functional

interactions are known to be dynamical in nature. Unlike in the linear case, it is hard to

determine the exact functional mapping and therefore without any apriori assumptions,

only the degree to which the two systems interact can be quantified [43]. This is exactly

how most of non linear tools operate. These tools, in general, assume that if two systems

are related, states of one system should be described by the states of another. Even

though these techniques have a wider sense of applicability, they still suffer from the lack

of analytical transparency, i.e. their inability to analyze the systems in depth. Another

major setback is the huge demand in the computational complexity involved.

In the literature, synchronization between systems has also been distinguished as

identical synchronization, phase synchronization and generalized synchronization. As far

as the tools are concerned, they can be classified into two categories, 1) bi-variate and 2)

multi-variate. As stated earlier, these classifications are made to suit the hypothetical

definition of what synchronization is. In the next section, we provide a brief literature

overview of the various signal processing tools designed to quantify the degree of

synchronization between systems.

1.3 Second Order Synchronization Measures

The literature on quantification of the spatio-temporal couplings in multi-

dimensional dynamical systems is rich and abundant, and for simplicity, it can be grossly

divided into linear and non linear. This section reviews some of the widely used second

order linear tools in this effort. We will discuss the common goals among these measures,

their approaches and the advantages and trade-offs along the way.









1.3.1 Cross-Correlation

Cross-correlation analysis is one of the earliest and most relied-on linear

techniques [6]. It measures the linear coupling between two signals and by design, is

symmetric.

Consider an identically distributed stationary stochastic process X. Assuming

ergodicity, represent a single time recording of the process as the discrete time-signal

x(n), n = 1,2,... M, where Mis the number of total number of samples in the signal.

Cross-correlation between two discrete time signals x(n) and y(n) can be estimated as a

function of lags, r, as follows:

1 Nr
R()- r= -N- 1,...,-1,,1.....N -1 (1.1)
(N r)

Due to the symmetric property, Rx (r) = (-) If u, uy are the mean-estimates of

x(n) and y(n) respectively, the cross-covariance between x(n) and y(n) can be estimated

as:

C, (r) = R (r) i (1.2)

Again, it is obvious that Cy (r) is symmetric. Cross-correlation coefficient between x(n)

and y(n) can now be defined as the cross-covariance normalized by the product of the

square root of the variances of the two signals. Mathematically,

C (r)
)3 (= (1.3)
~x y


where -x, i-y are the standard-deviations of x(n) and y(n) respectively. p,3 is

normalized between -1 and +1. The bounds represent strongest correlation (positive and

negative) while 0 represents independence.









Similarly, define X(o)) = F(x(n)) and Y(o)) = F(y(n)) as the Fourier transform

equivalents of x(n) and y(n) respectively. If the cross-spectrum C, (o) and the auto-

spectrums C, (o) and C, (o) are defined accordingly, normalized cross-coherence can

be mathematically represented as


^ (0) = (1.4)
c, (w).C, (w)


The cross-coherence quantifies the degree of coupling between X and Y at each

frequency. Again, the bounds are the same as in cross-correlation coefficient.

One of the problems with cross-correlation and the cross-coherence measures is

that they are symmetric. Symmetry fails to quantify the directional information and

therefore allows us to know only the average amount of information exchanged between

X and Y. Also, as with most other techniques, pair-wise computation is a limitation with

correlation/coherence analysis as well.

1.3.2 Partial Directed Coherence (PDC)

Efforts to improve on the correlation technique resulted in development of

multivariate auto-regressive modeling techniques such as directed coherence, directed

transfer functions (DTF) [19-20], and partial directed coherence (PDC) [2-3, 91]. It is a

frequency domain technique based on the concept of Granger-causality [27] which says

that an observed time series x,(n) causes another series x,(n), if knowledge of xj(n) 's past

significantly improves prediction ofx,(n). Along the same lines, the reverse case may or

may not be true. One way to make a quantitative assessment of the amount of linear

interaction and the direction of interaction among multiple time-series would be by

translating the concept of Granger-causality into a mathematical model such as a multi-









variable autoregressive model (MVAR). This will enable us to check the number of

predictable variations arising from using the knowledge of the past samples of another

time-series to predict the current sample of the time-series of interest.

The partial directed coherence (PDC) from to i at a frequencyfis given by

Ai (f)
~, (f)= a (f (1.5)
af (f)a (f)


1- a (r)e-jr,=
where A ) (f) = (1.6)
f a. (r)e-j2r, otherwise
r=\

a# (r) are the multivariate auto-regressive (MAR) coefficients at lag r, obtained by

finding the least-squares solution of the MAR model

p-1
x= A'x(p-r)+n (1.7)
r=0

Sr r r x (p r)
1a1 a"2 alN (p
x2 ( -r)
Here, A'r = and x(p-r)=

aN1 aN2 aNN_ XN(p r)

p represents the depth of the AR model, r represents the delay and n represents the

prediction error or the noise. Note that rij (f) represents the relative coupling strength of

the interaction of a given signal source with regard to signal i as compared to all ofj 's

interactions to other signals. It turns out that the PDC is normalized between 0 & 1 at all

frequencies. When i=j, the PDC represents the influence of a signal's past on its current

state. The PDC also differentiates between direct and indirect interactions among








multiple time series. Evidently, the MVAR models are better than most bivariate

measures because they use information simultaneously from all the channels, and thus are

able to unambiguously distinguish between direct and indirect causal connectivity

between nodes [41].

1.3.3 Synthetic simulation
Suppose that three simultaneously observed time-series are generated as follows

x, (n) = 0.1 sin(2fln) + sin(2zfn) + sin(2f3n) + w (n)
x2 (n) = 0.75x, (n 2) 0.1x2, (n -1)+ w (n) (1.8)
x3 (n) = 0.36x2 (n- 4) + Cx, (n -1) + w3(n)


where f =0.f,f 2=0.2f, f3 0.44f, sampling frequency f 100Hz; Coefficient C denotes

the coupling strength between xl and x3. w1, w2 & w3 are zero-mean uncorrelated white

noise processes with identical variances equal to 0.5.

One can observe from (1-5) that x2 is influenced by xl and x3 is influenced by x2 and can

possibly be influenced by x, as well, depending on the coupling strength C. Therefore, x3

can be influenced by x] directly or indirectly or both directly and indirectly. The coupling

state diagram is plotted in fig. 1-1.


0.-






Figure 1-1. Coupling diagram showing the direction of interactions between nodes, xl, x2
and x3.











PDO









1 1 2 0 0 1 2 3


--- 01 O1 2 O.
1=2 05 0.5 0.5




0 1 1


r= .=12 l--
0 1 20-S O1 2 0

(a)
PDC



o 11111 MBI ^nilnllfiffil llmmm0:]im mltafflff
0 1 2 3 0 1 2 0 1 2 3




0 1 2 3 1 2 3 0 1 2 3
1 1 1 --
0.5 0. 0.5



0 1 2 0 1 2 0 1 2 3
.= .=-2 =-3

(b)
Figure 1-2. Matrix layout plots for PDC, describing the interactions in example
simulation (1-8). 7tij indicates the influence of xj on xi at frequency f. The x-
axis limits range from 0 to pi, and the y-axis limits range from 0 to 1. a) Case
when C = 0.3 and b) C = 0.

Corresponding PDC matrix plots for the case when C = 0.3 are shown in fig 1-2a.


The plots clearly suggest a dependency of x2 and x3 on xl. Dependency of x3 on x2 is also


prominent from the plot. All the other plots (excluding the diagonal plots) indicate some


minor influences in the reverse direction, i.e, x2 to xl etc. However, a suitable statistical


threshold may clearly indicate an unambiguous influence among the nodes. The PDC









results seem to be in agreement, at least qualitatively, with the coupling diagram in Fig 1-

1.

For the case when there is no direct coupling between xl and x3, i.e, C=0, the plots

are described in fig. 1-2b. In this case, PDC indicates a much lesser influence of xj on x3,

as expected from its design.

Qualitatively, the PDC plots are reasonably good indicators of directional

dependence between multi-variate structures. However, it is difficult to exactly mark the

directions without setting an appropriate threshold. Quantitatively, the information on the

exact amount of coupling between two signals at a particular frequency as determined

by PDC, is ambiguous. In the simulation example, we created x2 and x3 using direct and

indirect influences from x1. Since x, is the driver and also is composed of multiple

sinusoids (f,f2 &f3), it is natural to expect the frequency-spectrum of x2 and x3 to mainly

consist of these sinusoidal frequencies. This also implies that the coherence between

these time-structures will be strong at these pole locations. However, the PDC plots fail

to provide that information, as we can see from fig. 1-2a and 1-2b. The coherence plots

are either smeared across all the frequencies or exhibit high values at frequencies other

than the expected frequency locations.

The MVAR approaches have been used to determine the propagation of epileptic

intracranial EEG activity in temporal lobe and mesial seizures [2-3, 10, 19-20]. However,

these models strictly require that the measurements be made from all the nodes, or the

directional relationships could be ambiguous. In addition, there is no clear evidence of

causality relationships among the cortical structures and as suggested in Rodriguez et al.

[84], the nature of synchronization is mostly instantaneous or without any detectable









delay. The applicability of MVAR measures is constrained to the linear dependencies

even though there are strong evidences pointing to non linear patterns of interactions

[23].

1.4 Mutual Information

Linear measures are typically restricted to measure statistical dependencies up to

the 2nd order. If signals are Gaussian distributed, the 2nd order statistics are sufficient to

capture all the information in the data. However, in practice, most signals are either non-

Gaussian or quassi-Gaussian rendering the linear statistical measures inadequate. As an

alternative, information theory measures [11-12, 66-68, 94-95] have been widely claimed

as acceptable choices to exploit the higher order statistical dependencies between signals.

A useful quantity is the mutual information approach that measures the Kullback-Liebler

(K-L) divergence between the joint probability and the product of the marginal

probabilities [11].

Consider that the vector samples x = {X1,....,XN} from an information source in a d-

dimensional signal space. Suppose that these samples are drawn from the distribution

f(x). Similarly, let the vector samples y = {y1,....,yN} form another information source in

a d-dimensional signal space belongs to distribution g(x). By definition, mutual-

information between x and y can be computed as follows


I(x;y) = JJ (x,y)log f(,y)dxdy
x,y f-(x)fy(y)
= Jfy (X, ) log fy (x,y)dxdy- ff(x)dx- f, (y)dy (1.9)
x,y x y
= H(x)+ H(y)- H(x,y)

where fx (x, y) is the joint probability event between x and y, and fx (x), f (y) are the

marginal probabilities, respectively. H(x) can be regarded as the average measure of









uncertainty and in the information theory literature is also called as Shannon's marginal

entropy, in x. Similarly H(y) is Shannon's marginal entropy in y, while their joint entropy

is denoted by H(x,y).

MI has been successfully used to quantify statistical couplings in biological

applications such as newborn cardio-respiratory systems [70-72], event related potentials

(ERPs) [65] and epileptic seizures [66-68]. One of the problems with this approach,

however, is that it requires large data sets for probability density estimation making it

computationally unviable. Efficient algorithms such as generalized mutual information

function (GMIF) have been recently proposed [70-72] to enable speedy computation

using very small data sequences.

1.5 Phase Synchronization

Conceptually, if the rhythms of one signal are in harmony with that of the other, the

two signals are known to be phase locked [99]. Phase synchronization can therefore be

defined as the degree to which two signals are phase locked.

Most of real world signals have broad spectra. For example, EEG signal recordings

are usually in the range of 0.1 to 1000 Hz even though they are usually band pass filtered

between 0.1 and 80 Hz since a major portion of the energy is contained in that spectrum.

The EEG can be classified roughly into five (5) different frequency bands, namely the

delta (0-4 Hz), theta (4-8 Hz), alpha (8-12 Hz), Beta (12-16 Hz) and the Gamma (16-80

Hz) frequency bands. Freeman [21] demonstrated evidence of phase locking between

EEG frequency bands across different regions of the brain, leading to certain clinical

events such as evoked potentials. Similarly, it is also believed that phase synchronization

across narrow frequency EEG bands, pre-seizure and at the onset of seizure [81-83] may

provide useful hints of the spatio-temporal interactions in epileptic brain.









The procedure normally adopted to measure phase synchronization consists of

initially tracking the instantaneous phases of the two signals of interest. This is followed

by computing their differences called the relative phase. Quantifying the distribution of

the relative phase will yield an index describing the degree to which the two signals are

phase locked. Hilbert transform [33, 49, 75, 81-83] is a widely used signal-transformation

trick to compute the instantaneous parameters of a time-signal. The following subsection

briefly describes the steps involved.

1.5.1 Hilbert Transformation

Consider a real-valued narrow-band signal x(t) concentrated around frequencyfc.

Define x(t)as

1
x(t)= x(t)*- (1.10)

x(t) can be viewed as the output of the filter with impulse response

1
h(t)=- -o < t < o (1.11)
2-t
excited by an input signal x(t). We call this filter a Hilbert transformer.

Next, we construct an analytical signal

z(t) = x(t)+ jx(t) (1.12)
which contains only positive frequencies. Our objective is to extract the instantaneous

parameters of the signal such as instantaneous amplitude, instantaneous frequency and

the instantaneous phase. z(t) can also be expressed in Cartesian co-ordinates as

z(t) = A(t)e(t) (1.13)
where A(t) is the instantaneous amplitude and 0(t) is the instantaneous angle of the

signal x(t). 0(t) can be expressed as


0(t) = cos((0 (t) + O(t))


(1.14)









(p(t) is the instantaneous phase of the signal and o)0 (t) is the instantaneous frequency of

the signal x(t).

Hilbert transforms are accurate only when the signals have a narrow-band

spectrum, which is unrealistic for most real-world signals. Pre-processing of the signal

such as decomposing it into narrow frequency bands is needed before we apply Hilbert

transformation to compute the instantaneous parameters. There are several ways of

breaking a signal into its sub-band components. Previous work [48, 76-80] to achieve

signal decomposition either points to a bank of digital FIR filters or the sub-band coding

using wavelet transforms. In the former, a multitude of conventional FIR filters are

designed to cater to the desired pass-bands using the Hamming or the Hann window. It is

well known from the signal processing literature that the transient events of the signal,

such as the spike activity and the sudden transitions that quite often exist in real signals,

are not adequately captured using the conventional FIR filters. In other words, it is

essential to capture the changing events in the signal locally at the time of their

occurrence. If the signal is non-stationary, the requirement for locality is certainly not

met with conventional FIR filters. The Wavelet filters [101] have shown better capability

of representing a signal in both the time and frequency domain and also capture the

transient details of the signal more effectively than the conventional digital FIR filters.

However, in both the techniques the decomposition does not always lead to narrow-band

signals. In other words, the decomposed signals do not conform to the definitions of a

narrow-band signal, even though they exist in a very narrow band frequency range.

Certain conditions need to be met to define a meaningful instantaneous frequency on a

narrow-band signal. They are as follows:









* the signal needs to be symmetric with respect to the local zero mean, and

* the signal should have the same number of zero crossings and extrema.

1.5.2 Empirical Mode Decomposition (EMD)

EMD is another preprocessing technique [34] in which the data are decomposed

into narrow band components called the intrinsic mode functions (IMF). Consider a

broad-band real-valued signal s(t). Using an iterative sifting process [34], the IMF

components are extracted from the signal s(t), iteratively. The first step in the sifting

process involves identifying the local minima and the maxima of s(t), and then using an

interpolation filter to combine the extremes. The local mean of the signal mi is computed

using the local extrema information and subtracted from s(t) i.e,

s(t)- m, = h,, (1.15)

The new signal is checked for the above 2 conditions. If it did not satisfy them, the

whole sifting process is repeated on hn1 to obtain say h12. This sifting process is

iteratively repeated until an IMF component, say hi, satisfying the two criteria is

obtained. The whole sifting process is now repeated on the residue signal si(t) obtained

from subtracting the IMF component hi from the original signal s(t) i.e,

s,(t) = s(t) -h, (1.16)

Eventually, all the successive IMFs satisfy the narrow band constraints. The second

condition is particularly essential in order to ensure that the instantaneous frequency does

not have any undesired fluctuations, introduced by asymmetric wave forms. Contrary to

previously discussed methods, the EMD derives IMFs or the basis functions directly and

adaptively from the data. This characteristic makes it particularly useful to detect local









changes in the signal. The basis functions can be empirically shown to be orthonormal

and complete.

We demonstrate the ability of the EMD to decompose a signal into its narrow-band

components by presenting a simple synthetic simulation. Consider a non-stationary

stochastic process X(t) consisting of a combination of sinusoidal components and noise

(constructed as shown in fig. 1-3). Components xl and x3 consist of 100 Hz (f1) and 450

Hz (f3) sinusoid oscillations respectively and span over the entire duration of the signal

X(t). x2 is a sinusoidal component of frequency 260 Hz (f) spanning only for some

period of the entire signal's duration. The signals were sampled at 1024 Hz (f)

frequency. Decomposition of X(t) into its narrow band components using EMD is shown

in fig 1-4.

The IMF components (fig 1-4a) and their corresponding frequency spectrum (fig 1-

4b) reveal a few things and can be listed as follows

1. All three frequencies are well separated by the EMD and also the middle panel in
fig 1-4a shows how the frequency component f2 exists only between the time-
duration that it was created in.

2. There exists leakage of f2 and f3 into IMF 1, but by looking at the power spectrum,
it is easy to observe that the energy at these frequencies is very low compared to fl.

3. Peak to peak amplitude of the extracted IMFs is almost half of the peak to peak
amplitudes in the original signal; however, this is just a scaling problem and can be
easily resolved by scaling up each of the IMFs.

4. The Hilbert spectrum, computed using the information from the instantaneous
parameters of the IMFs, reveals the joint time-frequency structure of the signal. The
extraneous frequencies are either created during the first transient period or due to
the presence of noise. Energy at these frequencies is below OdB.









17





X2















0


-2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
X 2



0


-2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1


X=x1 + X

+ X3 +


0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9


Figure 1-3. Synthetic Sinusoidal signal X(t) and its components xl, x2 and x3. X(t) also

consists of noise w-N(0, 0.1)


IMF components
2


0


-2
0 0.5 1
2





-2 0.5 1
O 0.5 1


0.5
time (sec)


Frequency spectrum


200



0


500


1000


400


200



0 500 1000
400 5


200



0 500 1000
frequency (Hz)


Figure 1-4. Decomposition of sinusoids using EMD. a) IMF components obtained from

successive decomposition of x(t) using the EMD b) Corresponding frequency

spectrum of the IMF components. IMF1 predominantly consists of 450 Hz,

IMF2 consists of 260 Hz and IMF3 consists of 100 Hz, as desired.


IMF1







IMF2








IMF3


-u,


ir










1.5.3 Wavelet decomposition

As described in previous sections, one of the many ways to achieve narrow-band

decomposition is by using ortho-normal wavelet filters [101]. For the same synthetic

signal X(t), we apply the '7 length Daubechies' filter functions, at various scales, to

achieve signal decomposition (fig 1-6).


IUU LZUU .UU 4UU UUU .UU
Time Sample


CUU U~U -UU IUUU


Figure 1-5. Hilbert spectrum constructed from the IMF functions. The spectrum reveals
the presence of 100 Hz and 450 Hz along the entire signal duration. Presence
of 260 Hz (corresponding to 0.26 on the Y-axis) component between the
200th and 800th sample (corresponding to approximately 0.2 and 0.8 second)
is also clearly visible from the spectrum plot. The leakage frequencies are
mainly due to the transient oscillations created in the first few samples of the
IMFs and also due to noise.

The following observations indicate the superior ability of EMD compared to wavelets in

a non-stationary scenario.

Two spurious frequency components of reasonably high energy, namely 150 Hz
and 240 Hz are created by wavelet analysis and synthesis (reconstruction) filters.

5. Qualitative comparison of Hilbert spectrum plots (fig 1-7 and fig 1-5) suggest a
larger smearing of the frequencies in wavelets compared to EMD.











6. The amount of energy leakage from the 260 Hz component into the first panel (fig.
1-4b) is considerably higher than the corresponding amount of energy leakage in
IMF (first panel in fig. 1-4b).


Wavelet reconstruction 10Frequency spectrum
5 1000 spectrum

0 500

-5 0
0 0.5 10 500 1000
5 1000

0 0.5 0 500 1000


0 0.5 1 0 500 1000
5 1000

3 0 500


0 0.5 1 0 500 1000



Figure 1-6. Decomposition of sinusoids using wavelets. Left) Wavelet reconstructed
components obtained from using Daubechies 7-scaled filters. Right)
Corresponding frequency spectrum of the wavelet reconstructed signals. W1
predominantly consists of 450 Hz, but also has a major portion of 240 Hz. W2
has more power concentrated in 260 Hz and 150 Hz frequencies and IMF3
consists of 100 Hz, as desired. Presence of spurious frequencies such as 240
Hz and 150 Hz are revealed using wavelet-decomposition.

1.5.4 Indexing Phase Synchronization

As we saw earlier, the Hilbert transformation is a useful time domain to a time-

domain transformation technique to compute the instantaneous phase of a signal. Due to

the narrow band constraints, any wide-band signal needs to be broken down into its

narrow band components and we saw that EMD and wavelet filters can be used as

potential pre-processors to achieve the same. Once a meaningful estimate of the

instantaneous phases is obtained, we need to determine the degree of phase locking by

quantifying them with a suitable metric.











Hilbert spectrum of wavelet reconstructed signals


0 45
0,4

- 035
C

o 025
.-
E 02
z0
0.15


01
0 US


100 200 300 400 500 600 700 800 900 1000
Time Sample

Figure 1-7. Hilbert spectrum constructed from the wavelet reconstructed narrow-band
signals. The spectrum reveals the presence of 100 Hz and 450 Hz along the
entire signal duration. Frequencies in the range of 200 Hz and 350 Hz
(corresponding to 0.26 and 0.35 on the Y-axis) are smeared between the 200th
and 800th sample (corresponding to approximately 0.2 and 0.8 second).
Presence of 260 Hz is not very clear from the spectrum plot.

Let (x (t) denote the instantaneous phase of one of the narrowband components of


x(t). Similarly, let op (t) be the instantaneous phase of the corresponding narrowband


component of y(t). Relative phase 0(t)can be computed as the difference between (x (t)


and py (t) i.e, 0(t) = ( (t) (p,(t), t = 1, 2,......, N. The relative phase or the phase-


difference 0(t) forms a distribution and a number of metrics can be used to compute an


index of phase synchronization by using the distribution of 0(t). Two popular metrics


exist in the literature [60-63, 81-83, 85, 99]. The first one expresses the phase


synchronization index as a deviation in the distribution of 0(t) from the uniform


distribution Odff [85, 99]. If we use the Shannon's definition for entropy, comparison of


the entropy H of 0(t) with entropy of a uniform phase distribution H,,, will result in,


H H
H X=- (1.17)
unif


V


0.5









Note that = 0 for complete lack of synchronization and is 1 for perfect phase-locking.

Another metric commonly used synchronization index is the phase locking value (PLV),

also called as mean phase coherence [81-83]. It is given by,

z= (ej() (1.18)

The bounds are the same as in previous index.

EMD is a useful pre-processing technique for decomposing a time-series into its

narrow band components. However, it does not provide any design parameter using

which a signal can be filtered into specific bands of interest. In phase synchrony analysis

where instantaneous phase is used as the synchrony metric to evaluate phase-locking

between two time-series, it is very important that the narrow band components of the

time-series are in sufficiently narrow (such as 2 to 4 Hz pass band range) spectra. While

it is possible to apriori design a conventional FIR filter catering to a desired pass band

and stop band, the EMD does not provide the luxury of choosing the desired frequency

bands. Due to the fact that EMD adaptively derives its basis functions based on the data,

the IMFs of two different signals can easily lie in different bands. Hence, another

limitation associated with using EMD is that, it is not possible to compare between the

IMFs of two signals.

Unlike the MVAR model based approaches, the applicability of phase-synchrony

measures is restricted to bi-variate time-series [10]. We also saw that both the

synchronization indices are symmetric and are therefore incapable of providing the

directionality of dependency. Also they are restricted to identifying only the phase

locking between time-signals. Synchronization between two systems is a general

phenomenon and can broadly be defined as the functional relationship governing the









oscillation of the two systems [88-92]. The function could be linear invertible, linear non-

invertible, non linear invertible and non linear non-invertible. Phase synchronization,

being only a subset of the general synchronization definition, fails to address the overall

functional relationship between systems.

1.6 Non-Linear synchronization measures

A lot of studies in the last two decades, to quantify nonlinear dependencies [88-93]

between two or more signals have met with reasonable success. As pointed out earlier,

most real world systems can be characterized as non linear processes [52-54] with their

samples dynamically evolving from interacting with past samples of the same process as

well as from the samples in the spatial neighborhood. Functional non linear mappings are

generally very hard to characterize and the fact that the recordings are embedded in noise

makes it harder.

A common procedure in any non linear measure is to construct a system map (or

the Poincare map) by embedding a scalar time series into a high dimensional Euclidean

space and building a trajectory by connecting the time points in that space [98]. This is

followed by computations that are generally directed towards determining the mean

distances between the mutual neighbors in the phase space of the signals. The

synchronization indices indicate the degree to which the state of one system is a function

of another and vice-versa.

1.6.1 Signal Reconstruction

Taken's embedding theorem states that under certain conditions, it is possible to

construct a finite dimensional dynamical system from a scalar time series {x, }I, where

xt= g(ut). utis the actual dynamical system, residing in dimension say 'd'. Taken [98]









postulated that, using time-delay embedding the exact dynamics of a system can be

reconstructed, while also preserving its topological characteristics. The reconstruction of

the attractor is done from a finite time series of the observation of a single variable x(t).

Thus a multi-dimensional embedding space is constructed from the time series data, and a

point in it represents the state of the system at a given time.

x(t,)
x(t, r)
x(t 2r)
x=< (1.19)


x(t, -(d- 2)r)
x(t (d )r)


Note that time-delay reconstruction relies heavily on the embedding parameters such as

the embedding dimension (d) and the delay time (z). Nevertheless, these parameters are

important since they convey a lot of information about the dynamical properties of the

underlying processes. The embedding dimension d can be selected by the false-nearest

neighbor criteria [98] and the delay time T can be estimated by the method of mutual-

information [98], or by measuring zeros of autocorrelation.

1.6.2 Non-Linear Synchronization Measures

Eckmann et al. [17] proposed the method of recurrence plots (RPs) that represents

the recurrence of states in the phase-space trajectory of a chaotic signal. Since chaotic

systems are non linear in nature, this method has been fairly successful in detecting

bifurcations and non-stationarities in time sequences. Cross-recurrence plot (CRP) is an

extension of the RP idea to multi-dimensional time signals [24]. The CRP has found

applicability in describing the time-dependency between multiple time-series recorded









from multiple locations. However, the lack of quantitative information and the

computational complexity makes it tedious for analyzing large sets of data. One of the

other drawbacks of this measure is that it fails to indicate the direction of information

flow (or influence). Variants of recurrence plots [59] are used to measure recurrence of

states in the phase-space between two chaotic signals. However, these plots are difficult

to analyze due to the lack of quantitative measures, in addition to their computational

complexity.

Several other non linear measures were proposed following the definition of

generalized synchronization by Rulkov and co-workers [88-92] Mutual predictabilities

were defined and studied extensively by Schiff et al. [92] and Le Van Quyen et al. [81].

Principally, if two systems are functionally related, observing the states of one system

should help predict the states of the other. A statistic constructed on this idea will indicate

the relative strength of coupling as well as the direction of coupling between any two

systems.

Recently, Arnhold et al. [1, 76] introduced the similarity-index technique (SI) to

measure such asymmetric dependencies between time-sequences. Conceptually, this

method relies on the assumption that if there is a dependency between two signals, the

neighboring points in time will also be neighboring points in state space. In other words,

if there is a functional dependency between two signals, then the recurrence of dynamics

of one signal will mean the recurrence of dynamics in the other signal.

Arnhold et al. [1, 76] propose to search for recurrence in the signals' state-space,

which poses an enormous computational burden, especially for large data sets. In this

dissertation, we propose a self-organizing map (SOM) based improvement to this method









to reduce computational complexity, while maintaining accuracy. This is achieved by

mapping the embedded data from signals onto a quantized output space through a SOM

specialized on these signals, and utilizing the activation of SOM neurons to infer about

the influence directions between the signals, in a manner similar to the original SI

technique. The approach not only facilitates real-time computation but it also enables us

to derive long-range synchronization patterns.

1.7 Objectives and Author's Contribution

Mechanisms involved in the synchronization among neural populations have

received a lot of attention in recent years [4-5, 8-9, 14-16, 42, 60-63, 97] It is agreed that

even in the normal physiological states, the neuron ensembles at the cortical structures

always interact and undergo a certain degree of synchronization. The Response to

sensory stimulations in the form of cohesive firing of neurons is known to lead to certain

transient changes. A larger question, however, that is still debated upon is whether

elevated synchrony in firing in a spatial neighborhood and its subsequent propagation to

other areas in the brain can lead to pathological situations such as seizures, stroke or

Alzheimer's.

Epileptic studies on adults and neo-nates have been widely researched [26, 35-40,

60-62, 82]. Researchers from different disciplines have contributed immensely in

understanding the various neurological aspects leading to seizures. Synchronization and

de-synchronization between the cortical networks are believed to be one of the plausible

reasons for this neurological disorder. It is soon becoming clear that epilepsy is a

dynamical disease [56], i.e. the macroscopic spatio-temporal dynamical across different

regions of the brain are consistent with rapid, sometimes gradual and often very subtle

nonlinear dynamical interactions. It is believed that synchronization occurs due to both









local and global discharges of the neurons. From the epilepsy perspective, quantifying the

changes in spatio-temporal interactions could potentially lead to the development of

seizure-warning systems. This quantification would also help us identify the regions that

actively participate during epileptic seizures.

In this dissertation, we apply the SOM-SI measure on epileptic seizure data to

investigate i) the temporal evolution of dependencies at different stages of a seizure, ii)

temporal changes in spatial connectivity of signals as a seizure is approached and iii)

spatial-temporal clustering of channels at different stages of seizure.

In Chapter 2, we present the SOM-based SI measure to quantify spatio-temporal

patterns. The SOM based SI approach to detect spatial interdependencies can be

generalized to any multi-variate dynamical system. Initially we present this idea from a

general perspective and then proceed to demonstrate the robustness brought about by the

SOM-improvisation through synthetic simulations. In chapter 3, we propose a spatio-

temporal clustering model that effectively utilizes similarity indices to achieve time-

series clustering. From an application standpoint however, the goal of this thesis is to

identify spatio-temporal connectivity's in epileptic seizures. In chapter 4 we present the

application of our model to epileptic intracranial EEG data and analyze the clustering

results on a large set of seizures. In chapter 5, we investigate the hypothesis that the

spatial patterns in intracranial EEG undergo significant changes as the brain progresses

towards a seizure state. This chapter also investigates if the dependencies in the focal

hemisphere are significantly different from the non-focal hemisphere, at pre-seizure and

post-seizure states. Chapter 6 presents a detailed discussion on conclusions and possible

directions for future research.














CHAPTER 2
SELF-ORGANIZING-MAP BASED SIMILARITY INDEX MEASURE

2.1 Introduction

In this chapter, we will mainly discuss the Similarity Index (SI) measure proposed

by Amhold et al. [1] and the proposed SOM based improvement to the SI measure.

Simulations with synthetic data are used to demonstrate the applicability and robustness

of the proposed SOM based SI measure.

2.2 The Similarity Index Technique

Synchronization between two signals Xand Yis understood, in principle, as a

functional mapping between Xand Y. In most neuro-biological systems especially, it is

realistic to expect the functional mappings to be nonlinear. The attractors of functionally

synchronous systems are related such that, the trajectory of one system partly describes

the trajectory of another and vice-versa. In unilateral coupling, if the trajectory of Y is

influenced by the trajectory of X, then Yis said to be functionally dependent on X

However, in instances where bi-directional relationships exist, the trajectory of Y can also

partly influence the trajectory of X In such a case, it is important to quantify two aspects

of synchronization: direction and strength.

Using the principles of generalized synchronization defined by Rulkov [88] and,

Arnold et al. [1] recently formulated a synchronization measure to characterize functional

relationships between non linearly coupled systems. Assume that Xand Y are two time

series generated by a system, which are embedded into two vector signals in time using

delays. N(X Y) is defined as the average dependency of Xon Y and it can be written as,










N(X ) N R(X)RXY) (2.1)
-Y
N, = R R"(X)
where R'(X) is the average Euclidean distance between the state-vector of X, and the

remaining state-vectors in X Y-conditioned Euclidean distance R'(X Y) measures the

average Euclidean distance between X, and the vectors in Xwhose corresponding time-

partners are the k-nearest neighbors of Y,. This measure takes values in [0, 1], where 0

implies no coupling and 1 implies perfect synchronization. In other words, 1 suggests that

recurrence of a state in Y implies a recurrence in X [1, 7, 8, 45, 78]. On the same

principles, N(X | Y) = 0 implies complete independence between Xand Y. By design, SI

can quantify nonlinear dependencies.

Similarly, it is possible to quantify the average dependence of Y on X by

1 R (Y)-Rn (Y I X)
N(Y I X) = R(Y) R"(Y ) (2.2)
N,_0 R"(Y)
Comparing N(X Y) and N(Y\X), we can determine which signal is more dependent on the

other. The difficulty with this approach is that at every time instant, we must search for

the k nearest neighbors of the current embedded signal vectors among all N sample

vectors; this process requires O(N2) operations. This high complexity hinders real-time

implementation and analysis. In addition, the measure depends heavily on the free

parameters, namely, the number of nearest neighbors and the neighborhood size s. The

neighborhood size s needs to be adjusted every time the dynamic range of the windowed

data changes.

2.3 SOM-based Similarity Index (SOM-SI)

The SOM-based SI algorithm is aimed at reducing the computational complexity of

the SI technique. The central idea is to create a statistically quantized representation of

the dynamical system using a SOM [74, 89]. A SOM is a neural-network in which spatial









patterns from the input-space are mapped onto an ordered output space consisting of a set

of neurons, called processing elements (PE). Thus each neuron in the SOM, based on its

location on the map, compactly models different features/dynamics of the input.

For best generalization, the map needs to be trained to represent all possible states

of the system (or at least with as much variation as possible). As an example, if we were

to measure the dependencies between EEG signals recorded from different regions of the

brain, it is necessary to create a SOM that represents the dynamics of signals collected

from all channels. The SOM can then be used as a prototype to represent any signal

recorded from any spatial location on the brain, assuming that the neurons of the SOM

have specialized in the dynamics from different regions.

One of the salient features of the SOM is topology preservation [74, 89]; i.e., the

neighboring neurons in the feature space correspond to neighboring states in the input

data. In the application of SOM modeling to the similarity index concept, the topology

preserving quality of the SOM feature of the SOM will be of added advantage, because of

the fact that the neighboring neurons in the feature space will now correspond to

neighboring states in the input data.

Assume that X and Y are two time series generated by a system, which are

embedded into two vector signals in time using delays. Define the activation region of a

neuron in the SOM as the set of all input vectors (the embedded signal vectors) for which

the neuron is the winner based on some distance metric (Euclidean in most cases). Let X,

be the set of time indices of input vectors x, that are in the activation region of the winner

neuron corresponding to the input vector x, at time n. Similarly define the set Y,.

Then the procedure to estimate the directed SOM-SI between X and Y is as follows









7. Train a SOM using embedded vectors from both X and Y as the input.

8. At time n, find Wf the winner neuron for vector xn, and find WY the winner
neuron for vector yn.

9. To find R"(X), compute the average Euclidean distance between Wx and all the
other winner neurons in the SOM. Similarly, compute R"(Y).

10. Determine the sets Xn and Y, for Wx and WY respectively.

11. Determine the nearest neurons WY corresponding to vectors y,, where j e Xn
Determine the nearest neurons Wn, corresponding to vectors y where j Yn.

12. Calculate R (X Y) = (1 / q) = W Wnxj II, where q is the number of elements in Xn.
CalculateRn (Y IX) = (1 / q)X =1 I W -wy I where q is the number of elements of
Yn.

13. Compute the ratios,

Sn"(X Y)= n (X)- Rn (X I Y))/ Rn (X) (2.3)

Nn (YX) =n(Y)-R n(Y X) /Rn(Y) (2.4)

14. Find interdependencies N(X I Y) and N(Y IX) as the average of N"(X I Y) and
Nn(Y X) overall n.

15. Compute the SOM-SI as the difference,

x = N(Y X)- N(X IY).

Positive values of indicate that influence of X on Y is more than the influence of Y on X,

while negative values indicate the opposite. Higher magnitude of indicates a stronger

coupling of the signals [29].

The computational savings of the SOM approach is an immediate consequence of

the quantization of the input (signal) vector space. The nearest neighbor search involves

O(NM) operations as opposed to O(N2) in the original SI, where Mis the number of PEs.









Traditionally M
compared to original SI.

2.4 Simulation Results

In this section, we demonstrate the viability of the SOM-based similarity index

approach in determining couplings and influence directions between synthetic and real

signals. One case study considers a coupled Rosseler-Lorenz system (as described in [7-

8, 76, 78]), and the other considers real EEG signals.

2.4.1 Rosseler-Lorenz simulation

The same Rosseler-Lorenz example used by Quiroga et al. [76] is used here. A

synthetic nonlinear dependency between a Rosseler (X) and a Lorenz (Y) system is

created by having the second state of the Rosseler system drive the Lorenz system in the

following manner

Roessler Lorenz
ii = -6{x, + x,} y = 10(-y, + y,)
X2 =6{x1 + 0.2x2} Jy =28y1 -y -y y3 y +Cx2 (2.5)
3 =6{0.2 + x3(x, 5.7))} 3 =y1y, y3
3
where C is the coupling strength. To be realistic, measurement noise, up to 30dB SNR are

added to both the systems. Two SOMs, one corresponding to x2 component of the

Rosseler system and the other, corresponding to the y2 component of the Lorenz system

(a new SOM for each value of C), were trained separately on embedded data, using an

embedding delay of 0.3 time-units and embedding dimension of 4. Note that separate

SOMs were trained because of the fact that the dynamics of the two systems were non-

identical. However, in general, while working with signals generated from the same

system, it suffices to train a SOM that represents all possible dynamics of that system.

SOM neurons overlapped with the phase-space dynamics of the Rosseler system and the









Lorenz system (for different C values) are shown in Fig. 2-1. Each SOM is an 8x8

rectangular grid, and is trained on a set of 4000 samples using a Gaussian neighborhood

function for about 400 iterations. The neighborhood radius (standard deviation of the

Gaussian neighborhood function) is exponentially annealed starting from an initial value

of 4 with a time constant of 100. The step size is also annealed exponentially from 0.08

using the same time constant.

Using the SOM-SI approach, the normalized indices Z are calculated for coupling

strengths of C= 0, 0.5, 1, 2, 3 and 5. The robustness to varying SNR is demonstrated in

the test data [31], as shown in Fig 2-2. Firstly, we observe that the asymmetric difference

in the interdependencies between the Rosseler and the Lorenz system increases as the

coupling (C) is increased (irrespective of the noise level). Despite the fact that the

training was done with 30dB SNR level, the absolute values of the dependencies do not

change much for SNR up to around 20dB, indicating noise-robustness. However, a

drastic decrease in the absolute dependency values for 10dB and OdB SNR situations

suggest that the structural relationship between the two systems is being largely

overwhelmed by the stochastic noise. In summary, we see that the idea of SOM modeling

of system's dynamics not only reduces computational complexity but also preserves

neighborhood mappings during noisy perturbations, thus providing enhanced noise

robustness.

















C=-0 R~euflrSyam



*



a ,

o. '0
0 0* L'1

o O
5 1 05 1
0 -i
o-`p"-0
* ; .


C *0.5. Lorn Sytm


o 0


C 3. Lornz Sy em


0O




0 1
o e

- 0 o


C 0. Lonnz Sy*m






0 0
0 *



0 0





o o


C -2, LoMI Sysm


0








,, o
oii *


12 a *2













(d)



C .5, Lorenz Syntm



o

0
2 o
I.s o ; <

S* o


*o o 0
100






a 0




S 45 0 1 14 2 L
a3 r 1 ->* n n i i* 2 l


Figure 2-1. Phase-space trajectories of the Rosseler-Lorenz system for various coupling

strengths a) Rosseler b) Lorenz (C=0) c) Lorenz (C=0.5) d) Lorenz (C=2) e)

Lorenz (C=3) f) Lorenz (C=5). The SOM weights (circles) for each signal are

superimposed on the trajectory.


1 0 1 2


e


*4


a
a


S1 0 1 2


3 21
















0.6

0.4

0.2,
SNR = 30dB

0 2 4 6



0.6 SNR = 10dB

0.4

0.2


0.6

0.4

0.2
SNR= 20dB
0 2 4 E
Coupling Strength (C)


0.6 SNR = OdB

-*- N(Roess Lor)
0.4 -e- N(LorlRoess)

0.2


0 2 4 6 0 2 4 6


Figure 2-2. Illustrating the dependency relationships between the Rosseler (driver) and

the Lorenz system (response) as a function of coupling strength.


Figure 2-3. Original EEG signals X and Y. for the EEG simulation example


v z w











W Iw II M I I I





Figu0 12-4. The nonlinearly mixed (s ) EEG snals V, Z, and W. 0 1 25




Figure 2-4. The nonlinearly mixed (synthetic) EEG signals V, Z, and W.




















Figure 2-5. Schematic diagram of the coupling function between the signals.

2.4.2 Intracranial EEG Simulation

In this example [30], intracranial EEG signals recorded from epileptic patients are

considered. Clinically, it is useful to know or find out the direction of information flow

through intracranial EEG signals. This analysis may help locate the epileptic foci of the

seizures as well as providing means of predicting them.

In this signal, since the intracranial EEG measurements are generated by a closed

system (the brain), we assume that the dynamical statistics of the signals observed at

different channels can be modeled using a single SOM, unlike the previous synthetic

example, where a separate SOM was used for each signal. The SOM, however, must be

trained using data that represents all possible psycho-physiological states that the

intracranial EEG signals might exhibit. In the case of an epilepsy patient, these include

pre-ictal, itcal and post-ictal states, in addition to the inter-ictal state.

In this example, intracranial EEG signals collected from two different patients at

different locations (labeled X and Y) are used. A synthetic nonlinear functional

relationship with influence direction from Xto V, W, and Z is created according to (2.6).

Care is taken in choosing these functions to make sure that the synthetic intracranial EEG









mixtures in (2.6) exhibit the characteristics of real EEG signals. The signals are shown in

fig. 2-3. This is achieved by verifying that the time structure

EEG SOM


4 .: '- o.-.'a.- "'
2 .*" o. ..




-2 -4 -2 0 .- -., 2 4 6




x(n)
Figure 2-6. The phase-space trajectory of the training EEG signal (projected in two
dimensions) and the weights of the trained EEG-SOM (circles).

and the power spectra of these signals are consistent with that of an EEG signal.
v(n) y(n) cx(n 2)x(n 3)

z(n) cv(n 2) cxzx(n 2) (2.6)

w(n) 0.05r(n)cx(n 2) v(n)






Here, x(n) and y(n) denote the original time sequences and v(n), w(n), and z(n) denote the

synthetic signals driven by the two original signals. In addition, r(n) is a zero-mean unit-

variance Gaussian noise term. The synthetic intracranial EEG signals are shown in Fig. 2-

4 and the flow diagram representing the relationships in (2.6) is depicted in Fig. 2-5.

A lOx10 rectangular SOM (referred to as EEG-SOM) is trained using 3000

samples of embedded intracranial EEG data (with an embedding dimension of 10 and
Figurembedding delay of 30ms) 2-6. The phase-space trajectory of the training EEG signal (projected in two
dimensions) and theweights of the trained M are shown in Fig. 2-6. The normal EG state is represented

andby the smapowller spectra of these signals are consistent withe training data), whereas the
v(n) = y(n) + cx(n 2)x(n 3)
z(-) = c~v( 2)+ c,,x( 2) (2.6)
w(n) = 0.05r(n) + cwz(n 2) + Cywv(n)

Here, x(n) and y(n) denote the original time sequences and v(n), w(n), and z(n) denote the

synthetic signals driven by the two original signals. In addition, r(n) is a zero-mean unit-

variance Gaussian noise term. The synthetic intracranial EEG signals are shown in Fig. 2-

4 and the flow diagram representing the relationships in (2.6) is depicted in Fig. 2-5.

A 10x10 rectangular SOM (referred to as EEG-SOM) is trained using 3000

samples of embedded intracranial EEG data (with an embedding dimension of 10 and

embedding delay of 30ms). The phase-space trajectory of the training data and the

weights of the trained SOM are shown in Fig. 2-6. The normal EEG state is represented

by the smaller amplitude activity (the dominant portion of the training data), whereas the









larger amplitudes correspond to the spiky, sharp, and slow wave activity formed during

the ictal state of the brain, or to artifacts formed due to muscle movements, etc. After

training the SOM, the normalized similarity index between the original signals Xand Y is

evaluated to verify whether these intracranial EEG signals are indeed independent or not.

The dependencies between X, V, W, and Z are also evaluated using the SOM to

calculate the normalized similarity index. The results are summarized in Table 2-1, where

both the coupling strength and the estimated similarity index between pairs of signals are

presented.

Table 2-1. Coupling strength between pairs of signals, the normalized similarity index
and the original similarity index between them.

Coupling Strength S(X, Y) = S(X I Y) S(Y I X)

c,= 1 -0.1668 X->V -0.1112 X->V

c, 2 -0.0756 x- >z -0.0612 x->z

cz -0.8 0.0901 z->v 0.0772 Z->V

cw = 1 -0.213 V- > W -0.336 V- >

c, = 0.3 -0.1225 z->W -0.1044 z->W



The results obtained from the SOM-based SI measure and the original SI measure

(in Table 2-1) is in perfect agreement. We conclude that X influences V and Z, V

influences Z and W, and W influences Z. Comparing these with the flow diagram in Fig.

2-5, it is seen that all directional couplings are consistent with the true construction

except for the relationship between V and Z. Possibly, this discrepancy is due to some

cancellations between the couplings from X and from V. Also, we can see that V is

exclusively constructed from the X and the Y components and does not have any









independent oscillations of its own, unlike W. These results indicate that the similarity

index approach might not produce results that are consistent with what one would expect

from the equations (if these are known) when the coupling diagram has closed loops.

2.5 Epileptic Intracranial EEG Data Description

Intracranial EEG signals were recorded from hippocampus, sub-temporal and

frontal-cortex structures of epileptic patients having a history of complex-partial, partial-

secondary and sub-clinical seizures of temporal-lobe focus, using bilaterally and

surgically implanted electrodes (Fig 2-7). Using amplifiers with an input range

of+ 0.6mv the recorded signals were converted to narrow-band using an anti-aliasing

filter with a cut-off range between 0.1Hz and 70Hz. Using an analog-to-digital converter

with 10-bit quantization precision, the narrow-band signals were sampled/digitized at 200

samples/sec. Measurements involved recording intracranial EEGs from multiple sensors

(28 to 32, with common reference channels) and the recordings spanned over 6

continuous days. A total of 55 seizures, of temporal lobe onset, were recorded from 5

patients, in the range of 6 to 18 seizures for each patient.

2.6 Testing the application of SOM-SI on Epileptic Intracranial EEG

In this section, we demonstrate the utility of SOM-SI in epileptic intracranial EEG

analysis and statistically verify accuracy of the results with the original SI. A 25x25

EEG-SOM is trained using 3000 input vectors constructed by embedding (dimension, m

= 10, delay, r = 30ms) intracranial EEG signals collected from various regions such as

temporal, sub-temporal, and orbit frontal, of an epilepsy patient. The EEG-SOM needs to

represent all possible EEG-dynamics, so the training data must include samples from the

inter-ictal, ictal, and the pre-ictal states of the patient. Fig. 2-8 shows the phase-space

trajectory of the data and the PEs of the EEG-SOM in two-dimensions. The normal











EEG state is represented by the smaller amplitude activity (the dominant portion of the

training data), whereas the larger amplitudes correspond to the spiky, sharp and slow

wave activity,


OIRBTOFRONTAL
{RODF



RIGIH'T
SIBTEMPORAL




TEMPORAL
(BrrD)


LEVI
ORBITOFRONTAL
(ILOF)



LErr




LEFT
FLMPTOD L


Figure 2-7. Diagram of the depth and subdural electrode montage in an epileptic brain.
Electrode strips are placed over the left orbitofrontal (LOF), right orbitofrontal
(ROF), left subtemporal (LST), right subtemporal cortex (RST). Depth
electrodes are placed on the left temporal depth (LTD) and right temporal
depth (RTD), to record hippocampus EEG activity.

mostly formed during the ictal state. We note that the distribution of the neurons is sparse

in the higher amplitude region because of the density matching property of the SOM.

To ensure generalization of the SOM, a test set of intracranial EEG signals were

quantized by the trained SOM. As seen from Fig. 2-9, the quantization results

successfully approximate the dynamics of the test data set (projected in one-dimension).

The correlation coefficient between the two signals was found to be 90.1%. For the most










part, the correlation coefficient was between 80% and 95%. Note that the amplitude

errors are higher in the larger amplitude regions corresponding to spike and slow waves.

This is expected because of the sparse distribution of the neurons in the higher amplitude

regions. These errors can be compensated by using a larger SOM grid (> 25x25), but

since the dynamics of the data are more important for the neighborhood information in

the SI measure and computational complexity will be an issue, we chose not to increase

the SOM grid size.

2.7 Statistical comparison between SI and SOM-SI

Next, we quantify the accuracy of the SOM-SI measure relative to the original SI

measure by comparing their results. SI values were calculated on nearly 39 minutes of

data corresponding to a pair of signals obtained from the right temporal (RTD4) and the

right sub temporal depth (RST1) electrodes, corresponding to patient P093. The entire

6
.. .... .._..

S "- -. -






-66 -4 -2 0 2 4 6
4 -.. .- .
-6 ------- ------- ------ '- --------------



X1
-6 -4 -2 0 2 4 6

Figure 2-8. The phase-space trajectory of the training intracranial EEG signal (solid
lines) superimposed on the weights (dots) of the trained 25x25 EEG-SOM
grid.

interval of 39 minutes data was segmented into 230 non-overlapping windows of 10

seconds each. Fig.2-10 shows the interdependency values of both the measures. It is easy

to see that the results from both the measures are in agreement to a large extent. There are

also subtle differences, which need to be quantified using statistical tests.



















200 300 400 500 600 700 800 900 1000
time (samples)
Figure 2-9. A qualitative illustration of the accuracy of the 25x25 EEG-SOM. Sample
test signal (solid line) overlapped on the SOM-reconstructed output (dash and
dot line). In this case, the correlation coefficient = 90.1%.

The comparison will be two-fold: (i) identify if the number of windows in which the

predicted directions of influence differ is significant or not, (ii) given time instances

where both measures agree on the direction, check if significant differences exist in

predicted strengths of influence. Assuming that SOM-SI and SI measure values come

from normal populations, we use the two-sided paired t-test to investigate the extent of

disagreement between the two methods. The test was performed at a significance level of

ca=0.05, over a size of 138 randomly selected samples out of the 230 available samples.

Null hypothesis: Ho: ,td = '(Zomsi-XZs) = 0

Alternate hypothesis: Hi: ,d = /,(ZsoX s-Zst) f 0

Paired t-test is chosen, because the observation in window 1 of original SI is

obtained under similar conditions as the window 1 of SOM-based SI, and hence, the data

may be said to occur in pairs. In this case, texp was found as -0.9441, whereas tcrit=t(o.o5),

137=1.960. Since texp
HO. This was also the case in 20 other comparisons made using different electrode pairs

from different patients. Therefore, we conclude that statistically the SOM-SI measure,

computed with a 25x25 grid SOM, performs as well as the original SI measure.













I0 .'-_-_-_-_- -I - I -__ _
0.4




-030 -25 -20 -15 -10 -5 0 5 10
Time ( minutes )-------------->
SOriginal Similarity Index
0)I





N(XIY)
-0.2- N(YIX)
-30 -25 -20 -15 -10 -5 0 5 10
Time ( minutes ) ----------------->

Figure 2-10. SOM-SI results showing the dependencies between two signals. X and Y
correspond to channels RTD4 and RST1, respectively. The time instant '0'
corresponds to the seizure onset (top). Results produced by the original SI
(bottom).

2.8 Testing the Robustness of SOM-SI on Multiple SOMs

The previous simulations successfully demonstrated the accuracy of the SOM-

based measure by statistically comparing its results with results from the original SI

measure. For application on seizures especially, a 25 x 25 SOM grid was trained to

embed all the dynamical states of the EEG attractor. SOM being one of the most

important elements of this improvised measure, one of the pre-requisites of this approach

is to ensure the following: a) that for data modeling purposes, the training set captures the

essence of variance found in the dynamics of the ictal states from all the channels, for a

given patient and b) the similarity indices computed using the SOM's processing

elements (PEs) are independent of the SOM and the corresponding training data set. In

other words, pair-wise similarity indices computed on two separate SOMs should be

significantly close to each other, if not equal.









We performed a simple test to check the latter point before proceeding with

extensive data analysis. From the multivariate intracranial EEG data samples of the

patient P093, two separate training sets were constructed. One of the training sets

consisted of portions of data sampled from the inter-ictal, ictal, pre and post-ictal states of

seizures 1 & 2. The other training set consisted of data portions picked around seizure 4

& 5. Using the same normalization procedures on both the sets and with the same set of

training parameters as before, two separate SOMs (called as SOM-1 and SOM-2 for

convenience) were trained. Following that, the SOM-similarity indices were obtained

from pair-wise analysis of interdependence among channels chosen from the ROF and

LOF regions of the brain, as illustrated in Fig 2-11.

The test data from the ROF and the LOF regions were picked from intervals

surrounding seizures 6 & 7, seizures 4 & 5 and seizure 11, respectively. Fig 2-12 shows

the exact time-intervals of the test data. The similarity index profiles {N1 (X|Y)t and {N2

(X|Y)t obtained from computing the SOM-SI on large intervals (say time t = 1,...,T) of

seizure data are quantitatively compared using the classical correlation coefficient and

error-percentage as the comparison metrics. The error-percentage is computed as follows


e}= 100" *N' (X I Y) Nz (X I Y),_
{ le} 100*{ x }(2.7)
N' (X I Y),

where N(XIY) is the normalized interdependency of X on Y. Note that the notations X

are Y are used to denote the two channels of interest. Normalized error e quantifies the

percentage difference between the interdependency values from SOM-2 and SOM-1,

keeping interdependency value from SOM-1 as the reference. From the error population,

the fraction of the absolute error values less than 20% and the fraction less than 10% are











computed to determine the degree of dependence of the SOM-SI measure on the data

used to train a SOM.





SSOMM-1 --






.... SOM-2




Figure 2-11. Experiment setup to compare SOM-Similarity Indices obtained from 2
separate maps.





1 hour 3 .urntes 1 hour 56 mirutes
1 hainr 48
minutes
Sz4 Sz5




17 98 minutes Appix 117 minutes

Sz6
S.7



3 hours 53 nuizues

S z 11

Figure 2-12. Time intervals from all the seizures used as test data (not drawn to scale)

For illustration, the results from analyzing the interdependency of LOF3 on LOF4

on various seizures are shown in fig 2-12. The histograms correspond to the error

ensembles obtained from analyzing over long seizure intervals. Qualitatively, the

superimposed traces in fig 2-13 indicate the extent of agreement or disagreement between

the SOM-SI profiles. Table 2-2 compiles a summary of the agreement between the SOM-









SI profiles for about 13 hours of intracranial EEG data from 5 seizures. A large fraction

of errors less than 20%, supported by a high correlation coefficient between the two

SOM-SI profiles suggests that there was very little disparity between the SOM-SI profiles

from SOM-1 and SOM-2. Besides, the high percentages also suggest the intracranial

EEG data dynamics do not vary drastically from one seizure to another, and therefore the

two SOM models produced the almost identical SI results. This finding will consequently

support (or reinforce) our original belief that a well-trained SOM and a well picked

training data set is sufficient to carry out independency analysis on all the seizures of a

patient.

Overall, pair-wise analyses of the interdependency among 6 channels (C6 = 15

combinations) on 5 seizures of P093 was performed on SOM-1 and SOM-2. The average

correlation coefficient and the error results between the SOM-SI profiles are shown in

Table 2-3.

Results from table 2-3 indicate that around 80% of the times, the differences

between the SOM-SI results are less than 20%. This is not surprising considering that the

differences are measured in percentages (2.7) and therefore even small discrepancies in

the case of small dependency values can appear magnified. In addition, we also speculate

that the discrepancies could be the outcome of the two SOMs being trained in an identical

fashion instead of being fine-tuned to obtain the lowest reconstruction error in each. In









Table 2-2. Quantitative comparisons between the SOM-SI profiles obtained from SOM-1
and SOM-2. LOF3 & LOF4 data was projected on each of the SOMs and
then the SOM-SI measure was applied to analyze the dependency of LOF3 on
LOF4.

P093, Correlation Fraction of error Fraction of error

Interdependency Coefficient (%) less than 20% less than 10%

N(LOF31 LOF4)

Seizure 6 & 7 95.74 0.8504 0.5597

Seizure 4 & 5 98.45 0.9234 0.7543

Seizure 11 91.59 0.6452 0.3614



general, if the SOMs can be designed to obtain the lowest reconstruction error, by

iteratively choosing the best sets of parameters, a slight improvement in the performances

can be easily achieved. But as it stands, a slight discrepancy can nevertheless be always

expected although it may have very little impact in the overall scheme of analysis.

2.9 Summary

The similarity index measure determines directional dependencies between two

signals using the basic assumption that two related signals will have similar recurrences

of the embedded state vector. This method has high computational complexity in terms of

the number of samples, since a search for nearest neighbors must be performed in the

phase-space of the signal. In this dissertation, we proposed a SOM-based approach as an

improvement over the original SI measure, for detecting functional dependencies among

multivariate structures. This approach reduces the computational complexity drastically

by exploiting the accurate quantization properties of the SOM in representing the

dynamics of the signal in the phase space. Another advantage of the SOM-based








47



approach is that the difficulties that the original similarity index approach encounters in


handling nonstationary data (such as the necessity to tweak parameters) are eliminated by


training the SOM using samples from various regimes of the nonstationary system. On

P093, Seizure 67 Patient P093, Seizure 6 & 7


i 0" l3 I
0.2


E0.1
O


Patient P093, Seizure 4 & 5


Patient P093, Seizure 11
0.4 -


.0.3


0.2

0.1

0

0

-0.1 time (minutes)
0 50 100 150 200 250


100


80


60


40


20

- o0


50 100


Error (%)


Figure 2-13. Comparing interdependencies between channels LOF3 and LOF4. Left:
SOM-similarity profiles from the output of SOM-1 and SOM-2 are
superimposed. Right: Histogram of the errors in %. Top: Seizure 4 & 5
Middle: Seizure 6 & 7. Bottom: Seizure 11.


60

50

40

30


error (%)


50 100


50


5(


Error (%)









Table 2-3. Summary of the comparisons between the SOM-SI profiles from SOM-1 and
SOM-2. Each row represents the statistics (mean and variance) of pair-wise
SOM-SI analyses of the epileptic intracranial EEG data from 6 channels (15
combinations).

Correlation Fraction of error Fraction of error
P093
Coefficient (%) less than 20% less than 10%


Seizure 6 & 7 94.32 + 2.85 0.79+ 0.1 0.54 + 0.12


Seizure 4 & 5 97.46 1.08 0.91 0.06 0.73 + 0.12


Seizure 11 93.24 + 2.06 0.71 0.08 0.41 0.07




the other hand, the SOM-based approach might suffer from inaccuracy if the quantization

is severe. Therefore, the size of the SOM could be decided by a trade-off between

representation accuracy and computational complexity. Through simulations, we also

demonstrated the noise-robustness features of the measure.














CHAPTER 3
SPATIO-TEMPORAL CLUSTERING MODEL

3.1 Introduction

In Chapter 2, we proposed the SOM-SI measure [29-31] and demonstrated SOM's

resourcefulness as a model infrastructure for computational purposes. Often, the time-

sequences in a multi-variate system share similar information that is reflected in their

interactive or synchronization abilities. By definition the word similar could mean that

the information shared among a set of channels is stronger than the information they

share with other channels. Such spatial similarities could possibly be momentary up to a

few seconds or could even stretch to several minutes or hours. As we postulated earlier,

spatio-temporal similarity changes could be one of the driving factors to trigger certain

events in biological systems. From a seizure point of view, we suspect that analyzing the

temporal changes in channel similarities could reveal some interesting aspects about the

epileptic brain.

Similarity-based time-series clustering [25, 44] is a well researched topic in the

area of dynamical graph theory. It is an extremely useful approach to characterize the

spatial groupings in time-sequences. Similar time sequences are typically grouped based

on their mutual interactions. In this study, using the SOM-SI as a computational tool to

derive the distance/similarity/proximity matrix, we propose a clustering model to

dynamically analyze the spatio-temporal groupings in muti-variate time sequences.

Simple toy simulations are presented to validate the robustness of the model. Since our









goal is specific to intracranial EEG, we use this model in chapter 4 to track the seizure-

related temporal changes in groupings.

3.2 Model for Spatio-temporal clustering

In this section, we propose a clustering approach to extract information on spatio-

temporal distribution of multivariate time-measurements. A 3-fold approach consisting of

spatial-discretization of the data using spectral-clustering technique [57, 64, 100],

temporal quantification using hamming distance, followed by application of another

clustering technique, is presented in Fig 3-1. The rational will become apparent during

the explanation.

Multichanmel
jata Time-delay Ol-I Spal-
-------Eniiddhg hg





Temponal Clusteriz
r Quantification
usiK Harmiir
Distarce

Figure 3-1. Block diagram to extract spatio-temporal groupings information in
Multivariate EEG structures

Spectral clustering is one of the many clustering methods that use subspace

decomposition on data-derived affinity matrix to achieve data-clustering. Using kernel

methods, the data samples are projected onto a higher dimensional space where the

discriminant analysis is much easier. Projecting the data onto a feature space results in

tightly formed clusters such that the between cluster entropy is maximized and within

cluster entropy is minimized. Spectral Clustering is inspired by the normalized cut theory

in computer science where the distance between the nodes/data samples is interpreted as









an affinity matrix on which subspace decomposition yield membership labels of the

nodes/data samples.

In our study, we use the standard spectral clustering algorithm by Ng et al. [64] to

spatially cluster the similarity-indices obtained by the SOM-SI technique. The algorithm

is outlined below,

1. Assuming N to be the number of nodes S ={si, s2, .... SN} and k as a priori
specified cluster size, form a similarity matrix A C R N using a distance metric.

2. Define D to be the diagonal matrix whose (i,i)th element is the sum ofA's ith row
and construct the matrix L D-/2A D-'/2

3. Find xl, x2, ...Xk, the k largest eigen-vectors of L and form the matrix E by
stacking [xl, x2, ....Xk] C RNxK the eigen-vectors in columns.

4. Form the matrix Y from X by renormalizing each of X's rows to have unit length.

5. Treating each row in Y as a point in Rk, cluster them into k clusters using k-means
or any other distortion minimization algorithm.

6. Finally, assign original node (si) to cluster j if and only if row i of the matrix Y was
assigned to cluster j.



Pair-wise evaluation of SOM-SI measure on all the possible combinations ((CN2,

where N is assumed to be the number of channels) of a portion of a multi-variate time

series leads to k = 2*(CN2) similarity indices in the bounds of [0, 1]. k is multiplied by 2

because of the asymmetric nature of the SOM-SI measure. If we imagine the time-series

as various inter-connected nodes in a multi-dimensional graph, the SOM-SI similarity

indices represent the affinity or rather, the weights of the connection, between those

nodes. Therefore, we can translate them into a square matrix of size Nx N, where N is the

# of channels. Since the weighting is normalized between 0 and 1, the diagonal elements

representing the affinity of a channel with itself, are coded as 1.









However, to be able to perform spectral-decomposition on an affinity matrix, Ng's

algorithm [64] requires that the affinity matrix be square and symmetric in nature. This is

because the eigen decomposition yields orthogonal column vectors (also called

eigenvectors) only if the projection matrix is square-symmetric. The asymmetric matrix

can be transformed to a symmetric matrix by adding it to its transpose and dividing each

entry by 2. Mathematically, if K, represents the asymmetric affinity matrix, the

transformation can be represented as,

(r, + kk)
3- = (3.1)
2

The transformed affinity matrix now represents the average information exchanged

between all pairs of channels. This implies that we do not have the luxury of using the

asymmetric nature of the dependencies to create a membership grouping among channels.

In epileptic study however, we will see that the averaging of the dependency information

will not affect channel dependencies because of the earlier observation that there is no

major difference in the driving and receiving information of the channels.

On the transformed affinity matrix 3, eigen decomposition can be done, followed

by K-means clustering. Consequently, we have a set of labeled clusters (number of

clusters is set, based on significant eigen values) representing the membership of the

channels.

If the above procedure is repeated over consecutive time windows (overlapped or

non-overlapped), channel groupings obtained on each time-window can be arranged as in

a matrix form as in (3.2).









3 2 2 3 1
1 2 2 3 2
Aspect = ... (3.2)


3 1 2 1 2

To characterize the average clustering of the channels over a longer period of time,

we propose another, albeit simple, hierarchical clustering approach that uses hamming

distance to derive the proximity matrix.

3.3 Temporal Quantification Using Hamming Distance

We showed in the previous section that the multivariate time-series can be grouped

by using similarity-based clustering techniques such as spectral clustering. The spectrally

clustered labels specify the groups of channels exhibiting high degree of within-cluster

similarities and low between-cluster similarities. Often in applications such as epileptic

EEG analyses, it is important to identify channel groupings over a longer period of time.

The epileptic seizures, in particular, are characterized by dynamic states (inter-ictal, ictal,

pre-ictal & post-ictal) that are known to possess both local and global spatio-temporal

clusters. Channels associate and de-associate in time; however, depending on the psycho

physiological state of the brain, certain groups of channels might have a higher likelihood

of sharing same channel labels thus forging a long-term association. In epileptic

intracranial EEG, identifying such state-dependent clusters may provide us with useful

insights on the evolution of brain patterns during seizure states.

State dependent connections can be quantified by clustering rows of the aspect

matrix that are similar with each other over a longer time interval, say T. Mutual

Information is one of the widely used metrics to quantitatively evaluate the similarity

between two discrete processes. However, in our study, using mutual information to









assess the degree of similarity between two channels/rows can have serious drawbacks.

For any two channels, two pairs at different time-instances, having same elements (say 3

& 1) need not bear any semblance. This is because the cluster labels are arbitrarily

labeled at each time-instance and do not have a fixed representation across time.

Technically, it implies that the probability term P(,-,,J=, i v j, where i,j are the

channels/rows in the Kicec matrix and u, v are the channel labels/memberships, can be

skewed unless u = v.

In this context, a very simple statistic to determine the degree of similarity would

be to compute the relative frequency of any two channels sharing the same

labels/groupings. In other words, in a time window of length T, we check the average

number of times when the two channels of interest, share the same cluster label.

In an algebraic context, the above operation is equivalent to computing pair-wise

hamming distance in a time window T. Of course, similarity would be obtained by

subtracting the hamming distance from 1. That is,

If dham is the hamming-distance between channels 'i' & 'j', similarity in probabilistic

terms is obtained as,

pim= 1 dham (3.3)

Thus, computing the pair-wise similarity for all i &j combinations will result in a P

matrix of size Nx N (Nis the number of channels). For convenience, we will call the

matrix P the cluster-similarity matrix in all our future references.

Hierarchical clustering on the cluster-similarity matrix P will yield information on

the cluster groupings over a time T. In the context of intracranial EEG data, clustering









will thus enable us to know the groups of channels that have similar behavioral structure

in the brain, over a wider time interval.

3.4 Simulations with Synthetic Data

This section presents simple simulations to demonstrate the validity of hierarchical

clustering approach to cluster a multivariate time-series data based on their mutual

interactions.

3.4.1 Roessler Lorenz System

Recall the Rosseler and the Lorenz system example from section 2.3.1. The idea is

to first generate non linearly coupled dynamical signals in each system and identify

mutual interactions using the SOM-SI measure. A hierarchical clustering algorithm such

as spectral clustering is then applied to group the time-series into pre-determined number

of clusters.

For the zero-coupling case (C = 0), 5000 data-samples from the non-transient

portion of each of the 6 components are extracted as shown in Fig 3-2a. Observe that the

time profiles of the first and the second components of each system have a strong phase-

similarity between them. The dynamics of each system evolve as a result of non linear

interactions within their components.

Ideally, only two SOMs suffice, each trained to represent the Rosseler and the

Lorenz system individually. However, if non linear coupling is introduced between the

two systems, the attractor dynamics change depending on the coupling strength. This

situation will force us to have individually SOMs to represent each dimension of a

system. Therefore, without any loss in generality, we train individual SOMs for each

component of the Roessler and the Lorenz system. Specific to this simulation, each SOM

was a 2-d grid of size 8 x 8 processing elements.









Affinity between the components is obtained by pair-wise computing the SOM-

similarity indices. Further, groupings among the 6 time-profiles are obtained by applying

spectral clustering algorithm. Spectral clustering, for a cluster-size n = 2, clear

distinguishes between Rosseler and the Lorenz components. For visual purposes, we

present the dendrogram (as in Fig 3-3a) generated in MATLAB 7.0. As expected, the

Rosseler and the Lorenz components form separate individual groupings. The fusion

levels also support the fact that the coupling strength C is very low (in this case = 0), by

clearly suggesting weak between-component interaction.

The next case consists of establishing a coupling (C = 5) from the X2, the 2nd

component of the Roessler to Y2, the 2nd component of the Lorenz system. Fig 3-2b

shows the time-traces of the Roessler and the Lorenz components. It is clear that Y2, (and

therefore the Y1, thelst component as well) has an altered time-profile. In addition, since

Y2 and Y3 of the Lorenz system also interact in a non linear fashion, we observe a clear

change in the time-dynamics of Y3. The procedure of applying SOM-SI to compute the

affinity matrix followed by spectral clustering is repeated. Dendrogram is shown in Fig.

3-3b, to graphically illustrate the cluster organization. Firstly, we can observe that the

topology of this dendrogram (for C = 5) is drastically different from the topology of

dendrogram for C = 0. We no longer see two clearly separated Rosseler and Lorenz

groupings. The fusion levels among the Lorenz components is still pretty much at the

same level as they were for C = 0. However, the fusion levels between the combined

Lorenz components and X1, X2 combined is reduced considerably, reflecting strong

coupling introduced earlier. The topology-difference between the two dendrograms in fig.








57



3-3b suggests that due to the non linear coupling introduced in (2.5) through the coupling


constant C, the dynamics of X2 & Xi are more closer to the dynamics to Y than to X3.


Roessler






500 1000 1500 2000 2500 3000 3500 4000 4500 5000




2- I


500 1000 1500 2000 2500 3000 3500 4000 4500 5000




0 .. .. .. .. .. .... .. . .. .. .. . .. .


500 1000 1500 2000 2500 3000 3500 4000 4500 5000


Lorenz


20 0 00 00 20 20 00 00 00 00
500 1000 1500 2000 2500 3000 3500 4000 4500 5000


time (samples)


Roessler






500 1000 1500 2000 2500 3000 3500 4000 4600 5000







500 1000 1500 2000 2500 3000 3500 4000 4500 5000







500 1000 1500 2000 2500 3000 3500 4000 4500 5000


time (samples)


Figure 3-2. Time profiles of the Roessler and the Lorenz components. a) Coupling
strength, C = 0 and b) C = 5









Table 3-1 shows the cluster groupings, obtained as a function ofn for various

coupling constants. While it is hard to obtain any distinction among the groupings, with n

= 4, the difference is clear for n = 2. Spectral Clustering, consistent with the model

construction and the observation in dendrogram, groups the X1, X2, Y1, Y2 & Y3,

components separately from X3.

Table 3-1. Cluster groupings for the simulation example.


3.4.2 Linear Model

Consider a set of colored noise data obtained as the band-pass filtered output of an

additive white-Gaussian stochastic process of zero mean and unit-variance (as shown in

Fig. 3-2). Each of the outputs can be denoted by s,(t), i = 1,2,.....,8 & t = 1,2,...,T. We

construe x,(t) as the observation measurements, sampled from eight (8) different channels

of a spatio-temporal system, over an interval T. Observe that x,(t) is constructed by

adding Gaussian white noise (of variance between 0.1 and 3) to s,(t).

xi(t) = si(t) + ni(t) (3.4)


C=0 C=5

# of clusters, n = 2 (Xi, X2, X3), (Xi, X2, Y1, Y2 ,Y3),

(Y1, Y2 ,Y3) (X3)

# of clusters, n = 3 (Xi, X2, X3), (Xi, X2), (X3),

(Yi, Y2), (Y3) (Y1, Y2, Y3)

# of clusters, n = 4 (Xi, X2), (X3), (Xi, X2), (X3),

(Yi, Y2), (Y3) (Y1, Y2), (Y3)









Design of noise n,(t) is such that, a variance higher than 3 will destroy the structure

ofx,(t) to the extent that the signals become completely un-correlated with each other.

Even though the two band-pass filters denoted by BPF1 and BPF2 respectively (as shown

in fig. 3-3 and fig. 3-4) share different pass-bands, there is a significant amount of

overlap as seen by the narrow separation of their pass-bands. It should also be noted from

fig. 3-4 that some of the output measurements x,(t) can share the same band-pass filters.

However the additive noise component n,(t) can introduce some stochastic differences in

their linear structures. Measurements x,(t) that share the same band-pass filter can loosely

be considered analogous to time-structures having reasonably strong linear-interactions.

This is because information exchanged between two time-series can sometimes cause

their phase dynamics to be similar.

We compute pair-wise cross-correlation indices (1) among all the channels to

quantify the interactions or rather, the exact amount of information exchanged among the

spatio-temporal measurements x,(t). Due to the fact that the signals are stochastic, linear

and do not possess any complex or chaotic phase dynamics, we felt that the cross-

correlation computation would suffice to indicate the amount of their pair-wise

interactions. The cross-correlation indices as denoted by

1 1 1
11 1,2 1,8
R2,1 R2,2 2,8
R=


R1 R2 R
8,1 8,2 8,8

denote the affinity matrix for the time series corresponding to say window 1. Repeating

this process over a number of windows Wwill result in the spatio-temporal affinity











Clustering on Roessler Lorenz components ; Coupling Strength, C = 0


0.2 0.3 0.4 0.5 0.6 0.7 0.8
Fusion Level


Clustering on Roessler Lorenz components ; Coupling Strength, C = 5


Fusion Level


(b)

Figure 3-3. Dendrogram illustration of similarity-based time series clustering, a)
Coupling strength, C = 0. The components of the two systems are clearly
clustered separately. b) For C = 5, the 1st and 2nd components of the Roessler
system form a stronger grouping with the components of the Lorenz system
than the 3rd component of the Roessler system itself.









matrix R =R1 R2 R where Wis the total number of windows. For each

window, the spatial configuration/arrangement of the band pass filters is changed as

follows

Configuration 1 (cl): Channels 1, 2, 3 & 4 are passed through BPF1 and

Channels 5, 6, 7 & 8 are passed through BPF2.

Configuration 2 (c2): Channels 1 & 3 are passed through BPF1 and

Channels 2, 4, 5, 6, 7 & 8 are passed through BPF2.

Configuration 3 (c3): Channels 1, 2, 3, 4, 6 & 7 are passed through BPF1 and

Channels 5 & 8 are passed through BPF2.

Configuration 4 (c4): Channels 1, 3, 6 & 7 are passed through BPF1 and

Channels 2, 4, 5 & 8 are passed through BPF2.

Total number of Windows Wis fixed to 100 and the probability of random occurrence of

channel configurations is as follows: {p(ci) = 0.5, p(c2) = 0.1,p(c3) = O., p(c4) = 0.3}.

Applying the spectral clustering algorithm on the spatial-correlation matrix of each

window Rk resulted in

1 2 1 2 1 1 2 1
1 2 2 1 2 2 2 1
1 2 1 2 1 1 2 1
1 2 2 1 2 2 2 1
aspect 2 1 2 1 2 2 1 2

2 2 2 2 1 2 1 1
2 2 2 2 1 2 1 1
2 1 2 1 2 2 1 2
T (8x100)


Cl C3 C2 C4 C4 C2 C1 C3










This illustration shows how the cluster labels are assigned to the channels depending on

their configuration, denoted by ci's. Heuristically, it is evident that the channels can

either be clustered in groups of four or in groups of two. For the case when the number

of clusters equals four (4), the hamming-distance based clustering resulted in cluster

assignment of the channels

C1: {1, 3), C2: {2, 4), C3: {6, 7), C4: {5, 8).

Specifying only 2 clusters resulted in

C1: {l, 2, 3, 4), C2: {5, 6, 7, 8).








BPF,

SBPF-----

BPF, ---3r ---
BPF -


BPF,





Figure 3-4. Model showing the generation of linearly dependent signals using
configuration 1.

We can easily see that the results with 4 clusters and 2 clusters are in perfect

agreement with the construction of the spatio-temporal configurations. Channel land 3

are always together and belong to the same cluster, regardless of the configuration.

Likewise, channels 5 and 8 always belong to the same cluster; although different from the

cluster to which 1 and 3 belong to. This is because the source data Z(t) is always passed







63


through BPF2 to obtain channels 5 and 8, unlike channels 1 and 3 that are obtained by


filtering Z(t) by BPF1. The channels 2 and 4 and similarly 6 and 7 change their


memberships often enough and therefore can be partitioned into separate clusters.


Figure 3-5. Plot showing the two band pass filters, separated by a narrow pass band


Pole Zero configuration for BPF1


Pole Zero configuration for BPF2


8 /
6 /
4

P1 -0.7 P,= 0.65
Z, 01.



6
8

-I
-1 -0.5 0 0.5 1
Real Part


Figure 3-6. Z-planes showing the pole- zero configurations for the two band-pass filters
in Fig 3-5.









3.5 Summary

Adopting a similarity-based time-series clustering approach, we proposed a spatio-

temporal model to evaluate long term clusterings in multi-variate time series data. The

earlier proposed SOM-SI measure is used as the distance measure to evaluate the affinity

between various nodes in time-series. One of the major issues in any clustering problem

is determining the appropriate number of clusters. In our model we use the spectral

clustering algorithm in which it relatively easier to specify the cluster size. This is

because the significant eigen values can directly provide a reasonable indication of the

number of centroids in the data. Similarly, through simulations we showed that the

dendrograms can also be used to get useful hints for choosing the number of clusters.

Simulations on Roessler-Lorenz dyamical system and on a linear set up demonstrate the

efficacy of our proposed approach to evaluate long term associability patterns among

dynamically connected nodes.














CHAPTER 4
APPLICATION OF SPATIO-TEMPORAL MODEL TO EPILEPTIC INTRACRANIAL
EEG

4.1 Introduction

In the last chapter, we proposed a spatio-temporal model to extract groupings from

long term multivariate recordings. This chapter will focus on the application of that

model on the epileptic intracranial EEG time series. The first part of the chapter will

describe the details on the application of the model and the second part will discuss the

results of analyses on 11 complex-partial seizures, from 2 epileptic patients.

4.2 Application of SOM-SI on Epileptic Intracranial EEG Data

The temporal changes in the spatial structure of an epileptic brain was analyzed on

twenty four (24) representative channels recorded bilaterally from the orbito-frontal,

temporal and sub-temporal regions on the brain. One of the fundamental requirements for

analyzing the dynamics of a non linear system is to construct the state-space attractor

from just a single recording of the time-series. From previous studies that estimated

intracranial EEG attractor size using correlation-dimension techniques [38-39], we know

that the EEG state space is bounded between 3 and 10. This is of course assuming that the

data is entirely noise-free and infinite amount of measurement data is available. The large

range in dimension is due to the fact that the properties of the EEG dynamical system

undergo dimension changes with changes in brain's physiological states and pathological

conditions. On our intracranial EEG data, the embedding dimension (m) and the delay









(r), however, were chosen to be m = 10 and r = 4. The parameters were compatible

with other studies [38-39], performed on the same data.

The following steps describe the procedure to track the spatio-temporal

connectivity patterns in intracranial EEG data.

1. Choosing suitable embedding parameters, intracranial EEG attractors in higher
dimensional state space is constructed using Taken's time-delay embedding
theorem [98]. On 10 second epochs, pair-wise analyses of interdependence among
24 channels are computed using the SOM-SI measure.

2. The similarity indices, from every window, are translated into an asymmetric
similarity/affinity/proximity matrix. Since spectral clustering involves eigen-
decomposition on a symmetric affinity matrix, the asymmetric matrix is
transformed to a symmetric matrix by adding itself to its transpose and dividing by
2.

3. With the number of clusters (say nl) specified apriori, spectral clustering on the
affinity matrix results in channels being labeled as one of the nl clusters.

4. Steps 1 to 3 are repeated for all the successive windows, representing 10-second
stationary segments. However, the overall ability of the channels to associate with
each other over a longer time-duration needs to be quantified.

On 30 minute time segments (equal to 90, 10-second windows), pair-wise Hamming-

distance based cluster-similarity matrix P is computed among all the channels. The

matrix elements essentially index the probability of channels to group into the same

cluster over a 30 minute time interval.

Spectral clustering or any other clustering algorithm on the cluster-similarity matrix

P will result in final cluster memberships. The number of clusters is fixed to n2.

For computing similarity indices in step 1, the epoch length of 10 seconds is chosen as a

tradeoff between stationarity and sample-size requirements. Also note that the successive

windows are 10 seconds apart (alternate 10 second windows), for reasons specific to

computational feasibility.









For illustration, we present the results of SOM-SI mutual interactions on a seizure

data in Fig 4-1. Pair-wise SOM-SI analysis of the interdependence between 24 channels

resulted in some interesting observations. Firstly, the SOM-similarity indices are

temporally averaged (over 3 windows) to smoothen out the fluctuations between

windows. On the smoothened interdependence indices, for each window and for each

channel, we find the maximum driving influence that the channel exerts on any other

channel. Over windows (in time) these maximum driving indices give the maximum

driving ability of the particular channel of interest. The maximum driving abilities are

evaluated for every channel under consideration and for simplicity; we show the driving

abilities for just 6 channels in fig 4-1. In the inter-ictal stage, low driving ability of all the

channels indicate that the channels are de-synchronized, even though they exhibit an

upward trend. Synchronization goes up momentarily a few minutes pre-seizure and at the

onset of the seizure, there is a sudden drop followed by a gradual increase post-seizure.

Higher degree of post-seizure global synchronization is followed by a gradual drop,

leading to the inter-ictal state.

It is clear from fig.4-1 that the temporal-evolution of the interdependency values

exhibit distinguishable patterns at different stages of seizure [30]. However, the

information on the spatial changes among the channels is still missing. As we stated

earlier, one of the primary objectives of this study is to determine the spatial-

configurations and their temporal changes, at various stages in a seizure data. In this

context therefore, clustering can be a useful tool to enable us to quantify groupings

among channels and progressively track their changes.










We now describe step 2, with more details. The channel interdependencies obtained

from SOM-SI, as indicated by Ki in (3.8) represent the spatio-temporal correlation

indices obtained by computing pair wise similarity-index among 24 channels. In spectral

clustering jargon, the k, matrix can also be interpreted as an affinity matrix representing

the pair-wise distances between 24 nodes.


ROF1
Z 0.9
LOF1
S0.8
RST1
-5 0.7 r
0.6 LST1
0.5 ^ RTD4
S0.4
E 0.3 LTD3

0.20
0.1

00 -80 -60 -40 -20 2 20 40 60 80
Time(mins)----->

Figure 4-1. Maximum average driving ability of each of the six (6) channels, nearly 100
minutes before and 70 minutes after Seizure-1 in patient P092. (The thin
vertical bar corresponds to the time when seizure occurred (0 to 2 on the x-
axis). For clarity, the box inside the figure shows a small portion of the
maximum average driving ability of each of the 6 channels, baseline-offset by
different scales.) Drop in synchronization followed by an abrupt increase in
phase synchronization at the onset of seizure is evident. Synchronization
across channels during seizure is also clearly seen

However, notice that the affinity between nodes will not be symmetric, as is usually

the case in spectral clustering. But the transpose transformation trick, mentioned in step

2, is applied to bring it to symmetric form. Since the weighting is normalized between 0

and 1, the diagonal elements representing the affinity of a channel with itself, are coded

as 1. Consequently, after spectral-clustering, we have a set of labeled clusters









representing the membership of the channels [32]. Repeating this procedure on every 10-

sec windows will yield a discrete-valued matrix specK t similar to (3.10).

Typically, the choice for the number of clusters nl in step 3 is conditioned on the

significant eigen-values (say 90 %). In our analysis, the sum of first 3 eigen-values

typically ranged from 60% to 80% of the total variance, due to changes in seizure states.

Considering this variability between epochs and the fact that the number of clusters need

to be the same for all epochs in order to be able to determine the overall grouping in

channels (using cluster-similarity matrix P), we fixed the membership size to nl = 3.

The notion of pre-seizure state has widely been debated upon. A lot of studies

argue against the notion of a 'pre-seizure' state transition. However, experimental studies

using non linear dynamics have shown that the quantitative descriptors of EEG exhibit

seizure pre-cursors in the form of inter-ictal to pre-ictal state transitions. The pre-ictal

transition time is not exactly known, however literature suggests that it has a broad range

of 5 minutes to 120 minutes before seizure. Therefore in step 5, as a tradeoff between

state transition periods and time-resolution, we choose 30 minutes time window to

characterize both the pre-ictal and the post-ictal periods.

4.3 Results

Following are the clustering results from applying the spatial-clustering procedure

on different seizures.

Patient P093

This patient had a history of complex partial seizures, localized in the mesial structures of

the temporal lobe. Surgery revealed a lesion in the right hippocampus (RTD electrodes)

region.









The set of 24 channels are listed below,

Channels 1 to 4: {LTD3, LTD5, LTD7, LTD9}

Channels 5 to 8: {RTD4, RTD6, RTD8, RTD10}

Channels 9 to 12: {LST1, LST2, LST3, LST4}

Channels 13 to 16: {RST1, RST2, RST3, RST4}

Channels 17 to 20: {LOF1, LOF2, LOF3, LOF4}

Channels 21 to 24: {ROF1, ROF2, ROF3, ROF4}

Before data analysis, a validation test was performed to check whether application

of different clustering algorithms on P would consistently result in same cluster

memberships or not. For a given number of clusters n2, it turned out that all the clustering

algorithms including spectral clustering produced the same outputs. Therefore, we

decided to choose the simple hierarchical clustering algorithm used in the Matlab 6.5

package owing to its graphical support.

Cluster-similarity matrices P indicating the probability that two channels share the

same grouping in a 30-minute time segment are shown gray-scale coded in Fig. 4-2. Pre-

seizure analysis on 30-minute windows is shown for up to 3 hours. Similarly, the post-

seizure analysis is shown for the first 30 minutes. The ability of the left side channels to

have a higher tendency to group together compared to the right hemisphere channels is

quite noticeable from fig 4-2. In addition, the orbito frontal lobes seem like the only area

on the brain to have a high probability of making a cross-hemisphere grouping. On the

left hemisphere, the LST and the LTD channels are consistently seen to share the same

clusters.









To confirm the observations from fig 4-2, the hierarchical clustering algorithm was

applied on each of those P matrices. Fig. 4-3 graphically illustrates two instances of the

clustering outputs through dendrograms. A dendrogram is strictly defined as a binary tree

with a distinguished root that has all the data items at its leaves. Conventionally, all the

leaves are shown at the same level of the drawing. The ordering of the leaves is arbitrary.

The heights of the internal nodes are related to the metric information (Phere) used to

form the clustering. Using a threshold of 0.4 and the average-linkage technique to

determine fusion levels, clustering was performed on a pre-defined number of clusters

(n2). For determining apriori the number of clusters n2, several dendrograms were

visually analyzed. There seemed to be at least 3 to 4 strong groupings among channels in

most of the dendrograms. For consistency, therefore, we chose to fix the number of

clusters n2 to 3 for all the analyses.

Both dendrograms in fig 4-3. clearly translate the spatial patterns observed in the

corresponding P matrices of fig 4-2. The top dendrogram in fig 4-3. corresponds to the

2.5 hours to 3 hours time window (indicated by -5) in fig 4-2. It is easy to see that the

dendrogram considers the RTD and the RST as isolated clusters due to their weak

between-cluster fusion level. Since the number of clusters n2 is restricted to 3, all the

remaining channels form a single large cluster. Similarly, the bottom dendrogram in fig

4-3. corresponds to the P matrix indicated by -1 in fig 4-2. In this case, the RST and the

RTD channels group into one cluster; also well supported by a dark patch in fig 4-2. This

enables the LST/LTD channels and the LOF/ROF channels to group together as separate

clusters.














LTD
LST1 -
R T DI ,--

R ST1 .--- ---- -L---- ---
ROF1 -------- r--------

LTND LST1 LOF1 RTD4 RST1 ROF1
LTN DS
LST1 ,. .- ..
I I
LOF1 ---- .----I----
RTD4 .., ,- "am- ,--
RSTi j----I -- L --- --- ----
ROF1RS1 ,--- .__L, ......
ROM.

LTDS LST1 LO1 RTD4 RST1 ROF1
LTDSN
LST1 ... .... --- H
LOF1 --.. ,-....
RTD~ ___ (-1)
RST1C----L-
RST1 J---- ---- --- ---
ROF1

LTDS LIT1 LOF1 RTD4 RST1 ROF1
LT- -
LST1 ---- L ---'- ---
I I
LO1 ----- --------
RTD4 ------.. --- -I--
RST1 ---- L --- ---- --- ----
ROMF1 ---- L- --. r -- -- ....

LT M LST1 LOF1 RT [4 RST1 ROF1


LTOB
LS1L...... I.... J....
LST1 -- -
LOF1 ------ ----I----
RTD-
RST1 ---- ----. --- -
I I II
ROF1 --------

LTD3 LST1 LOF1 RTO4 RST1 ROF1
LTIi-
LST1 ........ ....
LOF1 .-- ----- ----I----
RTI z.
RST1 ---- J ---- L ---
ROF1 ---- ---- --- -- -


LTD[ LST1 LOF1 RTD4 RST1 ROF1
LTC ST1
LST1 ------
LOF1 ---- --
RTD4 ----.-----r -- F---
RST1 ---- ----L--- ---


LTDS LST1 LOF1 RTO4 RST1 ROF1


Figure 4-2. Seizure 11 of patient P093: Cluster-similarity matrices P indicating the
probability that two channels share the same cluster label in a 30 minute time-
interval.


(-2)













































ROF3
ROF4
ROF2
LOF3
ROF1
LOF2
LOF4
LOF1
LST1
LTD9
LTD7
LTD5
LST2
LST4
LST3
LTD3
RST4
RST3
RST2
RST1
RTD10
RTD8
RTD6
RTD4


P093, 2.5 hours before Seizure 11


0.3 0.4
Fusion Level


P093, 0 30 minutes before Seizure 11



















-----------
----------

I I-I------


0 0.1 0.2 0.3 0.4
Fusion level


0.5 0.6


Figure 4-3. Dendrogram representation of the cluster results in Seizure 11, P093. TOP:
Dendrogram corresponding to 2.5 hours before seizure. BOTTOM:
Dendrogram corresponding to the 30 minute pre-seizure period.









The overall cluster configuration is listed in table 4-1.

Table 4-1. Spatio-temporal groupings as obtained for seizure


11 of patient P093


P093, Seizure 11 C1 C2 C3

Pre seizure, (2.5 3 hrs) RTD RST LTD, LST, LOF, ROF

Pre seizure, (2 2.5 hrs) RTD, RST LOF, ROF LTD, LST

Pre seizure, (1.5 2 hrs) RTD, RST LOF, ROF LTD, LST

Pre seizure, (1 1.5 hrs) RTD, RST LOF, ROF LTD, LST

Pre seizure, (30mins 1 hr) RTD, RST LOF, ROF LTD, LST

Pre seizure, (0 30mins) RTD, RST LOF, ROF LTD, LST

Post seizure, (30mins lhr) RTD, RST LOF, ROF LTD, LST



We summarize the spatial patterns at different time-intervals of seizure 11 as follows:

1. The LST and the LTD channels in particular, exhibit a strong tendency to belong to
the same group.

2. The LOF and the ROF channels form a strong bi-lateral homologous connection, as
seen from all the matrices in fig 4-2.

3. Relatively strong similarity can be seen between RTD and the RST channels.

4. Common observation in all the matrices is the strong similarity between the left
hemisphere channels as opposed to the right hemisphere channels. This is reflected
in the ability of LOF channels to have a higher probability of sharing clusters with
other left-hemisphere channels, as seen in Fig 4-2.

5. Interestingly, no temporal changes are seen in the spatial-patterns yet.



In Fig 4-3. even though the two dendrograms share certain common attributes, for

instance, channel labels within a group, observe slight disparities in topologies and in the

order of connections. Quantifying the difference between topologies may possibly serve

to differentiate between dendrograms, obtained at different states of a seizure. Statistical









techniques, such as Double Permutation Statistics (DPS) and mantel statistics are widely

[50, 51] known to be best suited to compare two dendrograms. However, the comparisons

usually test the hypothesis that the two dendrograms under comparison are not more

similar than the dendrograms randomly generated in terms of topology, leaf positions and

fusion level positions. Such statistical comparisons are particularly useful when the goal

is to detect whether the spatial-correlations have significantly changed over time or not.

We address this particular issue of tracking the temporal changes in spatial connections in

chapter 5. In the rest of this chapter, however, we analyze the spatial groupings for a

larger set of seizures, from different patients.

4.4 Statistical Validation

The cluster configurations observed from analyzing 30 minute segments

necessitates validation. Recall that, we partially validated our model (until the spectral

clustering stage), using synthetically coupled multivariate time sequences (non linear and

linear, both). Simulations involving creation of dynamic graphs involve multidimensional

time series that continuously change cluster memberships over time. Determining the

average spatio-temporal groupings from a collection of multi-variate time-series is

relatively easier to demonstrate in linear coupling cases. However, non linear dynamic

model constructions are extremely hard and mostly, non-trivial. We therefore decided to

pursue a posterior verification of the time-averaged cluster groupings on the intracranial

EEG data, using the quasi-surrogate analysis [40, 66-68, 73, 94-95], technique.

Recall that the cluster groupings obtained over 30 minute time-segments involve two

steps. First step consists of applying spectral clustering technique on the SOM-similarity

indices (computed on 10 second intracranial EEG data segments). Then similar grouping

patterns among channels are extracted by using hierarchical clustering approach on the









cluster-similarity matrices P. Specifically, 91 10-second windows (spanning 30 minutes,

3 windows per minute) are analyzed in each pass to obtain long-term groupings among

channels.

In order to validate this 2-step approach, we define our hypothesis as follows

Ho: The average within-cluster channel interaction at each window (out of 91 10-second

windows) is not significantly different from the corresponding between-cluster channel

interactions.

We propose to test this hypothesis on all the 3 (n2) clusters separately, for every 10-

second window within the 30 minute period.

Within-cluster interaction is computed by averaging the pair-wise similarity indices

for all the channels within a cluster. For between-cluster interaction, the pair-wise

interactions between 3 channels, picked randomly from each of the 3 clusters are

computed. A between-cluster interaction statistic is formed by computing the average

interactions from random selection of 3 channels (one from each cluster), over a number

of trials. We found that this statistic follows a quasi-normal distribution, implying that the

within-cluster interaction value can now be compared with the mean and the variance

sample estimates of the between-cluster statistic. Mathematically, we construct the z-

score as follows

c; -(cb
Z t = 1, 2,.. ,90 & i = 1, 2, 3. (4.2)
o-(Cbt(

where C,, is the within-cluster interaction at time '', for cluster 'i'. Cb ) is the mean


and c(Cb,) is the standard-deviation of the between-cluster interaction at time t'. Z,

reflects the z-score and is considered significant at the 95 percentile significance if Z' >









1.96 (reject Ho). In table 4-2 the bolded value in each cell represents the number of

windows (out of 91) having significant z-scores in the 30 minute period corresponding to

fig 4-2 (P093, Seizure 11). It is easy to observe that the null-hypothesis HO is rejected

beyond doubt, validating the clustering results.

Table 4-2. P093, Seizure 11: Over each 30 minute (91 samples total) window, number of
times the within-cluster interaction is greater than between-cluster interaction,
at 95% significance level.

P093, Sz 11 -5 -4 -3 -2 -1 0 (Sz) 1

C1 1 1 0.91 0.95 0.99 1 0.93

C2 0.82 0.89 0.96 0.91 0.89 0.85 0.98

C3 0.95 0.55 0.80 0.70 0.46 0.46 0.97



Seizures 4, 5, 6 & 7:

Spatio-temporal clustering analyses, similar to the one described on seizure 11

were performed on several other seizures, of the same patient P093. The cluster-similarity

matrices P obtained from time-intervals surrounding seizures 4 & 5 and 6 & 7 of patient

P093 are shown in fig 4-5. & fig 4-6. respectively. Channel groupings for the same are

listed in tables 4-3 & 4-4. All the 4 seizures present very consistent groupings.














0.1 0.1
0 10 20 so O 40 0

0.2 _0.2 -I_-.- &] -

0.2 ____________ __ 0.2
G0 0
so ro m nuB 100 110




0.2 C2-


01EO 150 200



Figure 4-4. Statistical validation of the clustering results. In each panel, thick lines are
used to represent the profiles of the three clusters in a 30 minute time interval.
The thin lines are the surrogate profiles indicating between-cluster
interactions. Cluster veracity can be visually verified by observing that
amplitudes representing within-cluster interaction for cluster profiles are
mostly higher that the amplitudes representing between-cluster interaction for
surrogate profiles, at each time-instance.

1. Consistent to the observation in seizure 11, we observe the temporal depth and the
sub-cortical regions of the left hemisphere are always bunched together.

2. Once again, the association of ROF-LOF areas into same cluster suggests a strong
homologous connection between the orbito-frontal areas of the brain. This
observation is also in agreement with those in seizure 11.

3. The dendrograms once again presented 4 unambiguous clusters in the form of RST,
RTD, LST/LTD and LOF/ROF. The fusion levels, indicating the strength of
connection between clusters, often turn out in favor of RTD and RST to be grouped
separately. Owing to the fact that we have pre-defined the number of clusters to 3,
the LST, LTD, LOF & ROF channels will consequently get grouped into one
cluster.

4. Once again, temporal changes are not very evident in the spatial patterns. However,
observing fig 4-5 & fig 4-6 and their corresponding dendrograms (not shown), the
fusion levels and the topology of the connections change with time. These changes
can only be quantified using statistical tests such as Mantel test statistics or the
Double Permutation Statistics (DPS).









Patient P092

In this section, we present the summary results of the clustering analyses performed

on patient P092 suffering from a lesion in the medial temporal lobe structures of the right

hemisphere.

Channel configuration for the patient P092 is as follows

Channels 1 to 4: {LTD1, LTD3, LTD5, LTD7}

Channels 5 to 9: {RTD2, RTD4, RTD6, RTD8, RTD12}

Channels 10 to 13: {LST1, LST2, LST3, LST4}

Channels 14 to 17: {RST1, RST2, RST3, RST4}

Channels 18 to 21: {LOF1, LOF2, LOF3, LOF4}

Channels 22 to 24: {ROF1, ROF2, ROF3}

Note that a separate 2 dimensional 25x25 sized EEG-SOM grid was created to model the

data dynamics of P092. Post spectral clustering analysis on 30 minutes data segments led

to some interesting observations.

Fig 4-7. shows the dendrograms created for seizure segments 2 hours prior to

seizure 1 and 30 minutes pre-seizure, respectively. As before, the number of clusters (nl)

specified in the spectral-clustering step after SOM-SI block was fixed to 3. The fusion

levels between most of the channel clusters is greater than 0.4, indicating a lack of strong

connectivity between regions.










Table 4-3. Spatio-temporal groupings as obtained for seizures 4 and 5 of patient P093.

P093, Seizure 4 & 5 C1 C2 C3

Pre seizure 4, (30 60 mins) RTD, RST LOF, ROF LTD, LST

Pre seizure 4, (0 30 mins) RTD, RST LOF, ROF LTD, LST

Post seizure 4, (0 30 mins) RTD LTD, LST, RST

LOF, ROF

Post seizure 4, (30 mins 1 hr) RTD LOF, ROF LTD, LST,

RST

Pre seizure 5, (30mins 1 hr) RTD LTD, LST, RST

LOF, ROF

Pre seizure 5, (0 30mins) RTD LTD, LST, RST

LOF, ROF

Post-Seizure 5, (30 1 hr) RTD LTD, LST, RST

LOF, ROF



Table 4-4. Spatio-temporal groupings as obtained for seizure 6 and 7 of patient P093.

P093, Seizure 6 & 7 C1 C2 C3

Post seizure 6, (0 30 mins) RTD, RST LTD, LST LOF, ROF

Pre seizure 7, (30 mins 1 hr) RTD, RST LTD, LST LOF, ROF

Pre seizure 7, (0 30 mins) RTD LTD, LST, LOF, ROF RST

Post seizure 7, (0 30 mins) RTD LTD, LST, RST LOF, ROF

Post seizure 7, (30mins 1 hr) RTD LTD, LST, LOF, ROF RST

Post seizure 7, (1 hr 1.5 hrs) RTD LTD, LST, LOF, ROF RST













LTD3
LTD3 -. - -
LST1 -
LOF1
RTD4 -- -
RST1
ROF1
LTD3 LST1 LOF1 RTD4 RST1 ROF1
LTD3
LST1
LOF1 N T
RTD4
RST1
ROF1
LTD3 LST1 LOF1 RTD4 RST1 ROF1
LTD3

LOF1 N-
RTD4 -
RST1 -
ROF1 -
LTD3 LST1 LOF1 RTD4 RST1 ROF1
LTD3
LSTI
LOF1
RTD4
RST1 .. .
ROF1
LTD3 LST1 LOF1 RTD4 RST1 ROF1


LTD3
LST1
LOF1
RTD4
RST1
ROF1
LTD3


LTD3
LSTI
LOF1
(-3) RTD4
RST1
ROF1
LTD3
LTD3
LSTl
LOF1
(-1) RTD4
RST1
ROF1
LTD3


;Zr4:- (4)

LST1 LOF1 RTD4 RST1 ROF1



(-2)


LST1 LOF1 RTD4 RST1 ROF1



L(0)



LST1 LOF1 RTD4 RST1 ROF1


Sz4 Sz5


Figure 4-5. Seizures 4 and 5 of patient P093: Cluster-similarity matrices indicating the
probability that two channels share the same cluster label in a 30 minute time-
interval.


I ,,





















LTD3 LST1 LOF1 RTD4 RST1 ROF1


LTD3 LST1 LOF1 RTD4 RST1 ROF1


LTD3 LST1 LOF1 RTD4 RST1 ROF1


LTD3
LST1
LOF1
(-2) RTD4
RST1
ROF1

LTC

LTD3
LST1
LOF1
(0) RTD4
RST1
ROF1

LTC

LTD3
LST1

(2) RTD4
RST1
ROF1


D3 LST1 LOF1 RTD4 RST1 ROF1


D3 LST1 LOF1 RTD4 RST1 ROF1


a-a-


LTD3 LST1 LOF1 RTD4 RST1 ROF1


30 to-60 '
m6 uZ
Sz6 Sz7


Figure 4-6. Seizures 6 and 7 of patient P093: Cluster-similarity matrices indicating the
probability that two channels share the same cluster label in a 30 minute time-
interval.


(3)


i










P092, 2 hours before Seizure 1


0 0.1 0.2 0.3 0.4 0.5
Fusion Level

P092, 30 minutes before Seizure 1


Fusion Level


Figure 4-7. Dendrograms corresponding to P092, Seizure 1. Top: 2 hours before Seizure
Bottom: 30 minutes pre-seizure.

For the second level of clustering, as before, let the number of clusters n2 be fixed at 3.

Cluster analysis on the 30 minutes segment 2 hours prior to seizure 1 (top dendrogram in

Fig. 4-7) results in the following groups of channels:

Cluster #1: LTD & LST

Cluster #2: RTD & RST

Cluster #3: LOF & ROF









Observe the cluster formed from LTD & LST channels, in the dendrogram. It is made up

of two sub-clusters, a large and a small cluster. The small cluster consists of only two

channels, LTD (3 & 5) and fuses with the other sub-cluster at a very high fusion level

(implying weak link). If n2 was to be increased to 4, the clustering algorithm would

classify this sub-cluster as an independent cluster. A detailed analysis on all seizures in

P092 revealed a strong intra-channel correlation (or low fusion level) between channels

LTD (3 & 5) and a weak inter-channel correlation with the rest of the channels. Surrogate

analysis also confirmed the imbalance by having very few rejections for the cluster

consisting of LTD (3 & 5) channels. It is obvious that the average interaction (within-

cluster interaction) of the largest cluster would be pulled down if there are sub-clusters

that have a strong within-sub-cluster interaction, but a weak between-sub-cluster

interaction. Consequently, the within-cluster interaction of the largest cluster can be

expected to be as weak as or marginally better than the between-cluster interactions,

leading to fewer rejections of the null hypothesis Ho. This problem can possibly be

overcome by increasing the number of clusters to 4 or more. However, for consistency,

we let the number of clusters n2 be fixed at 3 in the rest of the analyses.










85






LTD1
LST1 *
LOFI --- --
RTD2 (-3)
RST1
ROFM
LTD1 LSTI LOF RTD2 RST1 ROM


LTD1
LST1
LOF1 r )
RTD2 (-1)
RST1


LTD1 LST1 LOF1 RTD2 RST1 ROF


LTD1
LST1
LOF1 --- -- (0)
RTD2 -
RST1
ROF
LTD1 LST1 LOF1 RTD2 RST1 ROF

LTD1
LST1 -
LOF1
RT2 -- -- (2)
RST1

RO F1- - -
LTD1 LST1 LOF1 RTD2 RSTI ROM


-30 to -60

itiin5


Oto -33

iTiiwS


LIU' L I U
L G I- --LI-' ----

141U2 L-- r HII U2 --

--- -T- L LMG .1 I1 7I7-1
L I U' LISl' LGI- H I 2 H1yl II1"G LIU' LYI' LGI- HIU2 UY", 14GI-1


LI '
Lll.
LCIr -------------~---

141.u --r---.- ( )
1 U41' L----r--- ---

LIU' L1I" LGI- I- IU2 II' I GI-'


LI U' L I U'
L I I I I.-- LL r
LCH^ ----L- '-----3""-- LCH' ---- -- -----1.._
141U2 --- r---l--- 1-- 1,I .12 -- -----
14I: f1 UI


LIU' L 1' LGI' 14UJ 1 I' 1G LIU' L I' LGH- 1IU2 1, 1' 1 I-





Figure 4-8. Cluster-similarity matrices indicating the probability that two channels share

the same cluster label in a 30 minute time-interval, a) Seizure 1 of patient

P092 b) Seizure 3 of patient P092.


(-2)


LIU'
L I 4-- 1 - -.....
141U2 ---- -------

I. Y ----. .--.-.--.--- -. .

L IU' LI' LGI-' PII.2 14 11' 1GI-'









Seizures 1, 3, 4, 5, 6 & 7:


The spatio-temporal clustering results for seizures 1, 3, 4, 5, 6 & 7 are summarized

in tables 4-5 to 4-10.


Table 4-5. Spatio-temporal groupings as obtained for seizure 1


of Patient P092.


P092, Seizure 1 C1 C2 C3

Pre seizure, (1.5 2 hrs) RTD, RST LTD, LST (1,3,4) LOF, ROF, LST (2)

Pre seizure, (1 1.5 hrs) RTD LST, RST, LOF, LTD (3,5)

ROF, LTD (1,7)

Pre seizure, (30 mins 1 hr) RTD, RST LTD, LST LOF, ROF

Pre seizure, (0 30 mins) RTD, RST LTD, LST LOF, ROF

Post-Seizure, (0 30 mins) RTD, RST LTD, LST LOF, ROF

Post-Seizure, (30 1 hr) RTD LTD, LST, ROF

LOF, RST



Table 4-6. Spatio-temporal groupings as obtained for seizure 3 of Patient P092.

P092, Seizure 3 C1 C2 C3

Pre seizure, (1.5 2 hrs) RTD LST, LTD, RST LOF, ROF

Pre seizure, (1 1.5 hrs) RTD LST, LTD, RST LOF, ROF

Pre seizure, (30mins 1 hr) RTD, RST LST, LTD LOF, ROF

Pre seizure, (0 30mins) RTD, RST LST, LTD LOF, ROF

Post-Seizure, (0 30 mins) RTD, RST LST, LTD LOF, ROF

Post-Seizure, (30 1 hr) RTD, RST LST, LTD LOF, ROF










Table 4-7. Spatio-temporal groupings as obtained for seizure 4 of Patient P092.

P092, Seizure 4 C1 C2 C3

Pre seizure, (1.5 2 hrs) RTD, RST LST, LTD LOF, ROF

Pre seizure, (1 1.5 hrs) RTD, RST LST, LTD LOF, ROF

Pre seizure, (30mins 1 hr) RTD LST, LTD, RST LOF, ROF

Pre seizure, (0 30mins) RTD, RST LST, LTD LOF, ROF

Post-Seizure, (0 30 mins) RTD, RST LST, LTD LOF, ROF

Post-Seizure, (30 1 hr) RTD, RST LST, LTD LOF, ROF


Table 4-8. Spatio-temporal groupings as


obtained for seizure 5 of Patient P092.


P092, Seizure 4 C1 C2 C3

Pre seizure, (1.5 2 hrs) RTD LST, LTD, ROF

LOF, RST

Pre seizure, (1 1.5 hrs) RTD LST, LTD, ROF

LOF, RST

Pre seizure, (30mins 1 hr) LTD (3,5) LST, LTD (1,7), LOF, ROF

RST, RTD

Pre seizure, (0 30mins) LTD (3,5) LST, LTD (1,7), ROF

RST, RTD, LOF

Post-Seizure, (0 30 mins) RTD, RST LST, LTD LOF, ROF

Post-Seizure, (30 1 hr) RTD, RST LST, LTD LOF, ROF