UFDC Home  myUFDC Home  Help 



Full Text  
SPATIOTEMPORAL 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, Kyuhwa 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 SPATIOTEMPORAL 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 rossC 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 NonLinear synchronization measures ............. ............................. ....... ........... 22 1.6.1 Signal Reconstruction.................. ......... .......................... 22 1.6.2 NonLinear Synchronization Measures ........... ...................................23 1.7 Objectives and Author's Contribution ............... .......................................... 25 2 SELFORGANIZINGMAP BASED SIMILARITY INDEX MEASURE ..............27 2 .1 Introduction ................... ............................................... ................. 27 2.2 The Similarity Index Technique .............................................. ..... .......... 27 2.3 SOMbased Similarity Index (SOMSI)............................................................28 2.4 Simulation Results .................. ........................................ .... ........ 31 2.4.1 R osselerLorenz 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 SOMSI on Epileptic Intracranial EEG...................38 2.7 Statistical comparison between SI and SOMSI............................... .............40 2.8 Testing the Robustness of SOMSI on Multiple SOMs ....................................42 2 .9 S u m m ary ............. ............ .. ................................................... 4 6 3 SPATIOTEMPORAL CLUSTERING MODEL ............................................... 49 3 .1 Intro du action .................. ...................................... ......... ................ 4 9 3.2 M odel for Spatiotemporal 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 SPATIOTEMPORAL MODEL TO EPILEPTIC IN TR A C R A N IA L E E G ................................................................... .....................65 4.1 Introduction ................ ....................... .................. .................. .....65 4.2 Application of SOMSI 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 SPATIOTEMPORAL DEPENDENCY CHANGES IN EPILEPTIC INTRACRANIAL EEG...........................................91 5.1 Introduction ..................................................................... 91 5.2 Spatiotem 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 InterHemisphere 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 62 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 21 Coupling strength between pairs of signals, the normalized similarity index and the original similarity index between them........... ................. ..... ............37 22 Quantitative comparisons between the SOMSI profiles obtained from SOM1 and SOM2. LOF3 & LOF4 data was projected on each of the SOMs and then the SOMSI measure was applied to analyze the dependency of LOF3 on LOF4. .46 23 Summary of the comparisons between the SOMSI profiles from SOM1 and SOM2. Each row represents the statistics (mean and variance) of pairwise SOMSI analyses of the epileptic ECOG data from 6 channels (15 com binations). .........................................................................48 31 Cluster groupings for the simulation example ............................................. ........... 58 41 Spatiotemporal groupings as obtained for seizure 11 of patient P093 ..................74 42 P093, Seizure 11: Over each 30 minute (91 samples total) window, number of times the withincluster interaction is greater than betweencluster interaction, at 95% significance level. .................................................................... ..................77 43 Spatiotemporal groupings as obtained for seizures 4 and 5 of patient P093. ........80 44 Spatiotemporal groupings as obtained for seizure 6 and 7 of patient P093............80 45 Spatiotemporal groupings as obtained for seizure 1 of Patient P092. ..................86 46 Spatiotemporal groupings as obtained for seizure 3 of Patient P092. ....................86 47 Spatiotemporal groupings as obtained for seizure 4 of Patient P092. ....................87 48 Spatiotemporal groupings as obtained for seizure 5 of Patient P092. ....................87 49 Spatiotemporal groupings as obtained for seizure 6 of Patient P092. ....................88 410 Spatiotemporal groupings as obtained for seizure 7 of Patient P092. ....................88 LIST OF FIGURES Figure p 11 Coupling diagram showing the direction of interactions between nodes, xl, x2 and x 3 ............................................................................................................ 8 12 Matrix layout plots for PDC, describing the interactions in example simulation......9 13 Synthetic Sinusoidal signal X(t) and its components xl, x2 and x3. X(t) also consists of noise w N (0, 0.1) ................................................................................... 17 14 Decomposition of sinusoids using EM D .................................... ......... ............... 17 15 Hilbert spectrum constructed from the IMF functions .........................................18 16 Decomposition of sinusoids using wavelets ......................................................19 17 Hilbert spectrum constructed from the wavelet reconstructed narrowband sig n a ls ............. ......... .. .. ......... .. .. .................................................. 2 0 21 Phasespace trajectories of the RosselerLorenz system for various coupling strengths. ............................................................................33 22 Illustrating the dependency relationships between the Rosseler (driver) and the Lorenz system (response) as a function of coupling strength................................34 23 Original EEG signals X and Y. for the EEG simulation example .........................34 24 The nonlinearly mixed (synthetic) EEG signals V, Z, and W...............................34 25 Schematic diagram of the coupling function between the signals........................35 26 The phasespace trajectory of the training EEG signal (projected in two dimensions) and the weights of the trained EEGSOM (circles). ............................36 27 Diagram of the depth and subdural electrode montage in an epileptic brain...........39 28 The phasespace trajectory of the training intracranial EEG signal (solid lines) superimposed on the weights (dots) of the trained 25x25 EEGSOM grid. ...........40 29 A qualitative illustration of the accuracy of the 25x25 EEGSOM.. ......................41 210 SOMSI results showing the dependencies between two signals. X and Y correspond to channels RTD4 and RST1, respectively.......................................42 211 Experiment setup to compare SOMSimilarity Indices obtained from 2 separate m ap s. ............................................................................... 44 212 Time intervals from all the seizures used as test data (not drawn to scale) ............44 213 Comparing interdependencies between channels LOF3 and LOF4.......................47 31 Block diagram to extract spatiotemporal groupings information in Multivariate E E G stru ctu res............................. .................................................... ............... 50 32 Time profiles of the Roessler and the Lorenz components. a) Coupling strength, C = 0 and b) C = 5 .....................................................................57 33 Dendrogram illustration of similaritybased time series clustering, a) Coupling strength, C = 0. .........................................................................60 34 Model showing the generation of linearly dependent signals using configuration 1..... ......... .................................. ........ 62 35 Plot showing the two band pass filters, separated by a narrow pass band ..............63 36 Zplanes showing the pole zero configurations for the two bandpass filters in F ig 35 ...............................................................................63 41 Maximum average driving ability of each of the six (6) channels, nearly 100 minutes before and 70 minutes after Seizure1 in patient P092 .............. ..............68 42 Seizure 11 of patient P093: Clustersimilarity matrices P indicating the probability that two channels share the same cluster label in a 30 minute time interv al ............................ ...................... ..... ............. ................. 72 43 Dendrogram representation of the cluster results in Seizure 11, P093 ...................73 44 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 45 Seizures 4 and 5 of patient P093: Clustersimilarity matrices indicating the probability that two channels share the same cluster label in a 30 minute time interv al ................. .................................... ........................... 8 1 46 Seizures 6 and 7 of patient P093: Clustersimilarity matrices indicating the probability that two channels share the same cluster label in a 30 minute time interval .......... .............. ............................................ ................ 82 47 Dendrograms corresponding to P092, Seizure 1.............................................83 48 Clustersimilarity matrices indicating the probability that two channels share the same cluster label in a 30 minute timeinterval................................. ............... 85 51 Quantifying the spatial changes along time using Mantel statistics.........................99 52 Time smoothened profiles of the Mantel statistics during interictal states...........100 53 Illustrating the statistical difference between zscore at time't' and zscore at a reference time, 90, 100 and 110 minutes before a seizure. Top to Bottom: profiles of zscore statistics for seizures 1 through 3 of patient P092.................. 101 54 Illustrating the statistical difference between zscore at time't' and zscore at a reference time, 90, 100 and 110 minutes before a seizure. Top to Bottom: profiles of zscore statistics for seizures 4 through 6 of patient P092.................. 102 55 Illustrating the statistical difference between zscore at time't' and zscore at a reference time, 90, 100 and 110 minutes before a seizure. Top to Bottom: profiles of zscore statistics for seizure 7 of patient P092.................................... 103 56 Results of synchronization analysis on seizure 1 of P092, to check the differences between the focal and nonfocal hemispheric interactions, at pre and postictal states.. ..................................................................... 109 57 Nullhypothesis rejections rate averaged over 10 seizures for experiment 1.........110 58 Nullhypothesis 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 SPATIOTEMPORAL 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 spatiotemporal mappings in the brain. However, the exact spatiotemporal 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 nonlinear synchronization measure, called the SOMSimilarity Index, to quantify mutual associations between various brain regions. We propose to apply the Mantel test statistics on the SOMsimilarity indices to track the temporal changes of the spatial patterns. Statistical comparisons between interictal and preictal 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 ictalstates. 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 leftsubtemporal and the lefttemporal depth areas. In addition, strong homologous connectivity is found in the orbitofrontal regions. In the third aspect of this study, we investigate the differences in the synchronization levels between the focal and the nonfocal hemispheres at preictal and postictal states. Statistical tests confirm the existence of significant differences at pre ictal states followed by a strong entrainment, 20 minute postictal period. CHAPTER 1 SPATIOTEMPORAL 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 subsystem. These subsystems 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 timesequences. 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 timerecordings 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 multivariate system, understanding the interactions among its various nodes, whose behavior can be represented along time as timesequences, 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 timesignal analysis has become the forefront of research in discrete time signal processing [6672, 8587]. Time sequences recorded from realworld systems encode huge amounts of information that are both temporally and spatially extended. Especially in neural systems, synchronization between cortical structures or localfield potentials recorded at the tissue level is believed to be one of the prominent indicators of physiological and pathological events [2122]. 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 multivariate time series analyses have resulted in development of a wide array of signalprocessing 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, 6063]. Recently, multiresolution 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 nonlinear 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) bivariate and 2) multivariate. 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 spatiotemporal 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 tradeoffs along the way. 1.3.1 CrossCorrelation Crosscorrelation analysis is one of the earliest and most reliedon 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 timesignal x(n), n = 1,2,... M, where Mis the number of total number of samples in the signal. Crosscorrelation 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 meanestimates of x(n) and y(n) respectively, the crosscovariance 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. Crosscorrelation coefficient between x(n) and y(n) can now be defined as the crosscovariance 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, iy are the standarddeviations 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 crossspectrum C, (o) and the auto spectrums C, (o) and C, (o) are defined accordingly, normalized crosscoherence can be mathematically represented as ^ (0) = (1.4) c, (w).C, (w) The crosscoherence quantifies the degree of coupling between X and Y at each frequency. Again, the bounds are the same as in crosscorrelation coefficient. One of the problems with crosscorrelation and the crosscoherence 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, pairwise 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 autoregressive modeling techniques such as directed coherence, directed transfer functions (DTF) [1920], and partial directed coherence (PDC) [23, 91]. It is a frequency domain technique based on the concept of Grangercausality [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 timeseries would be by translating the concept of Grangercausality 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 timeseries to predict the current sample of the timeseries 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)ejr,= where A ) (f) = (1.6) f a. (r)ej2r, otherwise r=\ a# (r) are the multivariate autoregressive (MAR) coefficients at lag r, obtained by finding the leastsquares solution of the MAR model p1 x= A'x(pr)+n (1.7) r=0 Sr r r x (p r) 1a1 a"2 alN (p x2 ( r) Here, A'r = and x(pr)= 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 timeseries 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 zeromean uncorrelated white noise processes with identical variances equal to 0.5. One can observe from (15) 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. 11. 0. Figure 11. 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 20S 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 12. Matrix layout plots for PDC, describing the interactions in example simulation (18). 7tij indicates the influence of xj on xi at frequency f. The x axis limits range from 0 to pi, and the yaxis 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 12a. 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. 12b. 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 multivariate 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 frequencyspectrum of x2 and x3 to mainly consist of these sinusoidal frequencies. This also implies that the coherence between these timestructures will be strong at these pole locations. However, the PDC plots fail to provide that information, as we can see from fig. 12a and 12b. 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 [23, 10, 1920]. 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 quassiGaussian rendering the linear statistical measures inadequate. As an alternative, information theory measures [1112, 6668, 9495] 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 KullbackLiebler (KL) 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 ddimensional 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 cardiorespiratory systems [7072], event related potentials (ERPs) [65] and epileptic seizures [6668]. 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 [7072] 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 (04 Hz), theta (48 Hz), alpha (812 Hz), Beta (1216 Hz) and the Gamma (1680 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, preseizure and at the onset of seizure [8183] may provide useful hints of the spatiotemporal 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, 8183] is a widely used signaltransformation trick to compute the instantaneous parameters of a timesignal. The following subsection briefly describes the steps involved. 1.5.1 Hilbert Transformation Consider a realvalued narrowband 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) 2t 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 coordinates 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 narrowband spectrum, which is unrealistic for most realworld signals. Preprocessing 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 subband components. Previous work [48, 7680] to achieve signal decomposition either points to a bank of digital FIR filters or the subband coding using wavelet transforms. In the former, a multitude of conventional FIR filters are designed to cater to the desired passbands 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 nonstationary, 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 narrowband signals. In other words, the decomposed signals do not conform to the definitions of a narrowband 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 narrowband 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 broadband realvalued 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 narrowband components by presenting a simple synthetic simulation. Consider a nonstationary stochastic process X(t) consisting of a combination of sinusoidal components and noise (constructed as shown in fig. 13). 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 14. The IMF components (fig 14a) 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 14a 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 timefrequency 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 13. Synthetic Sinusoidal signal X(t) and its components xl, x2 and x3. X(t) also consists of noise wN(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 14. 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 narrowband decomposition is by using orthonormal 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 16). IUU LZUU .UU 4UU UUU .UU Time Sample CUU U~U UU IUUU Figure 15. 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 Yaxis) 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 nonstationary 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 17 and fig 15) 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. 14b) is considerably higher than the corresponding amount of energy leakage in IMF (first panel in fig. 14b). 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 16. Decomposition of sinusoids using wavelets. Left) Wavelet reconstructed components obtained from using Daubechies 7scaled 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 waveletdecomposition. 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 wideband signal needs to be broken down into its narrow band components and we saw that EMD and wavelet filters can be used as potential preprocessors 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 17. Hilbert spectrum constructed from the wavelet reconstructed narrowband 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 Yaxis) 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 [6063, 8183, 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 phaselocking. Another metric commonly used synchronization index is the phase locking value (PLV), also called as mean phase coherence [8183]. It is given by, z= (ej() (1.18) The bounds are the same as in previous index. EMD is a useful preprocessing technique for decomposing a timeseries 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 phaselocking between two timeseries, it is very important that the narrow band components of the timeseries 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 phasesynchrony measures is restricted to bivariate timeseries [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 timesignals. Synchronization between two systems is a general phenomenon and can broadly be defined as the functional relationship governing the oscillation of the two systems [8892]. The function could be linear invertible, linear non invertible, non linear invertible and non linear noninvertible. Phase synchronization, being only a subset of the general synchronization definition, fails to address the overall functional relationship between systems. 1.6 NonLinear synchronization measures A lot of studies in the last two decades, to quantify nonlinear dependencies [8893] 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 [5254] 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 viceversa. 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 timedelay 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 multidimensional 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 timedelay 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 falsenearest 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 NonLinear Synchronization Measures Eckmann et al. [17] proposed the method of recurrence plots (RPs) that represents the recurrence of states in the phasespace trajectory of a chaotic signal. Since chaotic systems are non linear in nature, this method has been fairly successful in detecting bifurcations and nonstationarities in time sequences. Crossrecurrence plot (CRP) is an extension of the RP idea to multidimensional time signals [24]. The CRP has found applicability in describing the timedependency between multiple timeseries 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 phasespace 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 coworkers [8892] 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 similarityindex technique (SI) to measure such asymmetric dependencies between timesequences. 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' statespace, which poses an enormous computational burden, especially for large data sets. In this dissertation, we propose a selforganizing 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 realtime computation but it also enables us to derive longrange 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 [45, 89, 1416, 42, 6063, 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 neonates have been widely researched [26, 3540, 6062, 82]. Researchers from different disciplines have contributed immensely in understanding the various neurological aspects leading to seizures. Synchronization and desynchronization 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 spatiotemporal 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 spatiotemporal interactions could potentially lead to the development of seizurewarning systems. This quantification would 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 evolution of dependencies at different stages of a seizure, ii) temporal changes in spatial connectivity of signals as a seizure is approached and iii) spatialtemporal clustering of channels at different stages of seizure. In Chapter 2, we present the SOMbased SI measure to quantify spatiotemporal patterns. The SOM based SI approach to detect spatial interdependencies can be generalized to any multivariate dynamical system. Initially we present this idea from a general perspective and then proceed to demonstrate the robustness brought about by the SOMimprovisation 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 spatiotemporal 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 nonfocal hemisphere, at preseizure and postseizure states. Chapter 6 presents a detailed discussion on conclusions and possible directions for future research. CHAPTER 2 SELFORGANIZINGMAP 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 neurobiological 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 viceversa. 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 bidirectional 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 statevector of X, and the remaining statevectors in X Yconditioned Euclidean distance R'(X Y) measures the average Euclidean distance between X, and the vectors in Xwhose corresponding time partners are the knearest 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 realtime 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 SOMbased Similarity Index (SOMSI) The SOMbased 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 neuralnetwork in which spatial patterns from the inputspace 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 SOMSI 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 SOMSI 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 SOMbased similarity index approach in determining couplings and influence directions between synthetic and real signals. One case study considers a coupled RosselerLorenz system (as described in [7 8, 76, 78]), and the other considers real EEG signals. 2.4.1 RosselerLorenz simulation The same RosselerLorenz 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 timeunits 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 phasespace dynamics of the Rosseler system and the Lorenz system (for different C values) are shown in Fig. 21. 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 SOMSI 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 22. 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 noiserobustness. 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 21. Phasespace trajectories of the RosselerLorenz 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 22. Illustrating the dependency relationships between the Rosseler (driver) and the Lorenz system (response) as a function of coupling strength. Figure 23. Original EEG signals X and Y. for the EEG simulation example v z w W Iw II M I I I Figu0 124. The nonlinearly mixed (s ) EEG snals V, Z, and W. 0 1 25 Figure 24. The nonlinearly mixed (synthetic) EEG signals V, Z, and W. Figure 25. 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 psychophysiological states that the intracranial EEG signals might exhibit. In the case of an epilepsy patient, these include preictal, itcal and postictal states, in addition to the interictal 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. 23. 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 26. The phasespace trajectory of the training EEG signal (projected in two dimensions) and the weights of the trained EEGSOM (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 zeromean 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. 25. A lOx10 rectangular SOM (referred to as EEGSOM) is trained using 3000 samples of embedded intracranial EEG data (with an embedding dimension of 10 and Figurembedding delay of 30ms) 26. The phasespace trajectory of the training EEG signal (projected in two dimensions) and theweights of the trained M are shown in Fig. 26. 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 zeromean 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. 25. A 10x10 rectangular SOM (referred to as EEGSOM) is trained using 3000 samples of embedded intracranial EEG data (with an embedding dimension of 10 and embedding delay of 30ms). The phasespace trajectory of the training data and the weights of the trained SOM are shown in Fig. 26. 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 21, where both the coupling strength and the estimated similarity index between pairs of signals are presented. Table 21. 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 SOMbased SI measure and the original SI measure (in Table 21) 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. 25, 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, subtemporal and frontalcortex structures of epileptic patients having a history of complexpartial, partial secondary and subclinical seizures of temporallobe focus, using bilaterally and surgically implanted electrodes (Fig 27). Using amplifiers with an input range of+ 0.6mv the recorded signals were converted to narrowband using an antialiasing filter with a cutoff range between 0.1Hz and 70Hz. Using an analogtodigital converter with 10bit 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 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 SOMSI on Epileptic Intracranial EEG In this section, we demonstrate the utility of SOMSI in epileptic intracranial EEG analysis and statistically verify accuracy of the results with the original SI. A 25x25 EEGSOM 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, subtemporal, and orbit frontal, of an epilepsy patient. The EEGSOM needs to represent all possible EEGdynamics, so the training data must include samples from the interictal, ictal, and the preictal states of the patient. Fig. 28 shows the phasespace trajectory of the data and the PEs of the EEGSOM in twodimensions. 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 27. 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. 29, the quantization results successfully approximate the dynamics of the test data set (projected in onedimension). 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 SOMSI Next, we quantify the accuracy of the SOMSI 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 28. The phasespace trajectory of the training intracranial EEG signal (solid lines) superimposed on the weights (dots) of the trained 25x25 EEGSOM grid. interval of 39 minutes data was segmented into 230 nonoverlapping windows of 10 seconds each. Fig.210 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 29. A qualitative illustration of the accuracy of the 25x25 EEGSOM. Sample test signal (solid line) overlapped on the SOMreconstructed output (dash and dot line). In this case, the correlation coefficient = 90.1%. The comparison will be twofold: (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 SOMSI and SI measure values come from normal populations, we use the twosided paired ttest 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 = '(ZomsiXZs) = 0 Alternate hypothesis: Hi: ,d = /,(ZsoX sZst) f 0 Paired ttest is chosen, because the observation in window 1 of original SI is obtained under similar conditions as the window 1 of SOMbased 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 SOMSI 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 210. SOMSI 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 SOMSI 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 prerequisites 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, pairwise 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 interictal, ictal, pre and postictal 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 SOM1 and SOM2 for convenience) were trained. Following that, the SOMsimilarity indices were obtained from pairwise analysis of interdependence among channels chosen from the ROF and LOF regions of the brain, as illustrated in Fig 211. 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 212 shows the exact timeintervals of the test data. The similarity index profiles {N1 (XY)t and {N2 (XY)t obtained from computing the SOMSI on large intervals (say time t = 1,...,T) of seizure data are quantitatively compared using the classical correlation coefficient and errorpercentage as the comparison metrics. The errorpercentage 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 SOM2 and SOM1, keeping interdependency value from SOM1 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 SOMSI measure on the data used to train a SOM. SSOMM1  .... SOM2 Figure 211. Experiment setup to compare SOMSimilarity 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 212. 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 212. The histograms correspond to the error ensembles obtained from analyzing over long seizure intervals. Qualitatively, the superimposed traces in fig 213 indicate the extent of agreement or disagreement between the SOMSI profiles. Table 22 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 SOMSI profiles suggests that there was very little disparity between the SOMSI profiles from SOM1 and SOM2. 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 welltrained SOM and a well picked training data set is sufficient to carry out independency analysis on all the seizures of a patient. Overall, pairwise analyses of the interdependency among 6 channels (C6 = 15 combinations) on 5 seizures of P093 was performed on SOM1 and SOM2. The average correlation coefficient and the error results between the SOMSI profiles are shown in Table 23. Results from table 23 indicate that around 80% of the times, the differences between the SOMSI 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 finetuned to obtain the lowest reconstruction error in each. In Table 22. Quantitative comparisons between the SOMSI profiles obtained from SOM1 and SOM2. LOF3 & LOF4 data was projected on each of the SOMs and then the SOMSI 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 phasespace of the signal. In this dissertation, we proposed a SOMbased 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 SOMbased 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 213. Comparing interdependencies between channels LOF3 and LOF4. Left: SOMsimilarity profiles from the output of SOM1 and SOM2 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 23. Summary of the comparisons between the SOMSI profiles from SOM1 and SOM2. Each row represents the statistics (mean and variance) of pairwise SOMSI 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 SOMbased approach might suffer from inaccuracy if the quantization is severe. Therefore, the size of the SOM could be decided by a tradeoff between representation accuracy and computational complexity. Through simulations, we also demonstrated the noiserobustness features of the measure. CHAPTER 3 SPATIOTEMPORAL CLUSTERING MODEL 3.1 Introduction In Chapter 2, we proposed the SOMSI measure [2931] and demonstrated SOM's resourcefulness as a model infrastructure for computational purposes. Often, the time sequences in a multivariate 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, spatiotemporal 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. Similaritybased timeseries 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 timesequences. Similar time sequences are typically grouped based on their mutual interactions. In this study, using the SOMSI as a computational tool to derive the distance/similarity/proximity matrix, we propose a clustering model to dynamically analyze the spatiotemporal groupings in mutivariate 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 Spatiotemporal clustering In this section, we propose a clustering approach to extract information on spatio temporal distribution of multivariate timemeasurements. A 3fold approach consisting of spatialdiscretization of the data using spectralclustering technique [57, 64, 100], temporal quantification using hamming distance, followed by application of another clustering technique, is presented in Fig 31. The rational will become apparent during the explanation. Multichanmel jata Timedelay OlI Spal Eniiddhg hg Temponal Clusteriz r Quantification usiK Harmiir Distarce Figure 31. Block diagram to extract spatiotemporal groupings information in Multivariate EEG structures Spectral clustering is one of the many clustering methods that use subspace decomposition on dataderived affinity matrix to achieve dataclustering. 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 similarityindices obtained by the SOMSI 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 eigenvectors of L and form the matrix E by stacking [xl, x2, ....Xk] C RNxK the eigenvectors 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 kmeans 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. Pairwise evaluation of SOMSI measure on all the possible combinations ((CN2, where N is assumed to be the number of channels) of a portion of a multivariate 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 SOMSI measure. If we imagine the timeseries as various interconnected nodes in a multidimensional graph, the SOMSI 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 spectraldecomposition 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 squaresymmetric. 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 Kmeans 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 nonoverlapped), channel groupings obtained on each timewindow 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 timeseries can be grouped by using similaritybased clustering techniques such as spectral clustering. The spectrally clustered labels specify the groups of channels exhibiting high degree of withincluster similarities and low betweencluster 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 (interictal, ictal, preictal & postictal) that are known to possess both local and global spatiotemporal clusters. Channels associate and deassociate 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 longterm association. In epileptic intracranial EEG, identifying such statedependent 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 timeinstances, having same elements (say 3 & 1) need not bear any semblance. This is because the cluster labels are arbitrarily labeled at each timeinstance 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 pairwise 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 hammingdistance between channels 'i' & 'j', similarity in probabilistic terms is obtained as, pim= 1 dham (3.3) Thus, computing the pairwise 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 clustersimilarity matrix in all our future references. Hierarchical clustering on the clustersimilarity 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 timeseries 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 SOMSI measure. A hierarchical clustering algorithm such as spectral clustering is then applied to group the timeseries into predetermined number of clusters. For the zerocoupling case (C = 0), 5000 datasamples from the nontransient portion of each of the 6 components are extracted as shown in Fig 32a. 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 2d grid of size 8 x 8 processing elements. Affinity between the components is obtained by pairwise computing the SOM similarity indices. Further, groupings among the 6 timeprofiles are obtained by applying spectral clustering algorithm. Spectral clustering, for a clustersize n = 2, clear distinguishes between Rosseler and the Lorenz components. For visual purposes, we present the dendrogram (as in Fig 33a) 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 betweencomponent 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 32b shows the timetraces of the Roessler and the Lorenz components. It is clear that Y2, (and therefore the Y1, thelst component as well) has an altered timeprofile. In addition, since Y2 and Y3 of the Lorenz system also interact in a non linear fashion, we observe a clear change in the timedynamics of Y3. The procedure of applying SOMSI to compute the affinity matrix followed by spectral clustering is repeated. Dendrogram is shown in Fig. 33b, 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 topologydifference between the two dendrograms in fig. 57 33b 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 32. Time profiles of the Roessler and the Lorenz components. a) Coupling strength, C = 0 and b) C = 5 Table 31 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 31. Cluster groupings for the simulation example. 3.4.2 Linear Model Consider a set of colored noise data obtained as the bandpass filtered output of an additive whiteGaussian stochastic process of zero mean and unitvariance (as shown in Fig. 32). 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 spatiotemporal 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 uncorrelated with each other. Even though the two bandpass filters denoted by BPF1 and BPF2 respectively (as shown in fig. 33 and fig. 34) share different passbands, there is a significant amount of overlap as seen by the narrow separation of their passbands. It should also be noted from fig. 34 that some of the output measurements x,(t) can share the same bandpass filters. However the additive noise component n,(t) can introduce some stochastic differences in their linear structures. Measurements x,(t) that share the same bandpass filter can loosely be considered analogous to timestructures having reasonably strong linearinteractions. This is because information exchanged between two timeseries can sometimes cause their phase dynamics to be similar. We compute pairwise crosscorrelation indices (1) among all the channels to quantify the interactions or rather, the exact amount of information exchanged among the spatiotemporal 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 pairwise interactions. The crosscorrelation 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 spatiotemporal 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 33. Dendrogram illustration of similaritybased 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 spatialcorrelation 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 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: {l, 2, 3, 4), C2: {5, 6, 7, 8). BPF, SBPF BPF, 3r  BPF  BPF, Figure 34. 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 spatiotemporal 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 35. 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 36. Zplanes showing the pole zero configurations for the two bandpass filters in Fig 35. 3.5 Summary Adopting a similaritybased timeseries clustering approach, we proposed a spatio temporal model to evaluate long term clusterings in multivariate time series data. The earlier proposed SOMSI measure is used as the distance measure to evaluate the affinity between various nodes in timeseries. 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 RoesslerLorenz 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 SPATIOTEMPORAL MODEL TO EPILEPTIC INTRACRANIAL EEG 4.1 Introduction In the last chapter, we proposed a spatiotemporal 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 complexpartial seizures, from 2 epileptic patients. 4.2 Application of SOMSI 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 orbitofrontal, temporal and subtemporal regions on the brain. One of the fundamental requirements for analyzing the dynamics of a non linear system is to construct the statespace attractor from just a single recording of the timeseries. From previous studies that estimated intracranial EEG attractor size using correlationdimension techniques [3839], we know that the EEG state space is bounded between 3 and 10. This is of course assuming that the data is entirely noisefree 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 [3839], performed on the same data. The following steps describe the procedure to track the spatiotemporal connectivity patterns in intracranial EEG data. 1. Choosing suitable embedding parameters, intracranial EEG attractors in higher dimensional state space is constructed using Taken's timedelay embedding theorem [98]. On 10 second epochs, pairwise analyses of interdependence among 24 channels are computed using the SOMSI 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 10second stationary segments. However, the overall ability of the channels to associate with each other over a longer timeduration needs to be quantified. On 30 minute time segments (equal to 90, 10second windows), pairwise Hamming distance based clustersimilarity 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 clustersimilarity 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 samplesize 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 SOMSI mutual interactions on a seizure data in Fig 41. Pairwise SOMSI analysis of the interdependence between 24 channels resulted in some interesting observations. Firstly, the SOMsimilarity 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 41. In the interictal stage, low driving ability of all the channels indicate that the channels are desynchronized, even though they exhibit an upward trend. Synchronization goes up momentarily a few minutes preseizure and at the onset of the seizure, there is a sudden drop followed by a gradual increase postseizure. Higher degree of postseizure global synchronization is followed by a gradual drop, leading to the interictal state. It is clear from fig.41 that the temporalevolution 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 SOMSI, as indicated by Ki in (3.8) represent the spatiotemporal correlation indices obtained by computing pair wise similarityindex among 24 channels. In spectral clustering jargon, the k, matrix can also be interpreted as an affinity matrix representing the pairwise 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 41. Maximum average driving ability of each of the six (6) channels, nearly 100 minutes before and 70 minutes after Seizure1 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, baselineoffset 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 spectralclustering, 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 discretevalued matrix specK t similar to (3.10). Typically, the choice for the number of clusters nl in step 3 is conditioned on the significant eigenvalues (say 90 %). In our analysis, the sum of first 3 eigenvalues 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 clustersimilarity matrix P), we fixed the membership size to nl = 3. The notion of preseizure state has widely been debated upon. A lot of studies argue against the notion of a 'preseizure' state transition. However, experimental studies using non linear dynamics have shown that the quantitative descriptors of EEG exhibit seizure precursors in the form of interictal to preictal state transitions. The preictal 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 timeresolution, we choose 30 minutes time window to characterize both the preictal and the postictal periods. 4.3 Results Following are the clustering results from applying the spatialclustering 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. Clustersimilarity matrices P indicating the probability that two channels share the same grouping in a 30minute time segment are shown grayscale coded in Fig. 42. Pre seizure analysis on 30minute 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 42. In addition, the orbito frontal lobes seem like the only area on the brain to have a high probability of making a crosshemisphere 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 42, the hierarchical clustering algorithm was applied on each of those P matrices. Fig. 43 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 averagelinkage technique to determine fusion levels, clustering was performed on a predefined 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 43. clearly translate the spatial patterns observed in the corresponding P matrices of fig 42. The top dendrogram in fig 43. corresponds to the 2.5 hours to 3 hours time window (indicated by 5) in fig 42. It is easy to see that the dendrogram considers the RTD and the RST as isolated clusters due to their weak betweencluster 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 43. corresponds to the P matrix indicated by 1 in fig 42. In this case, the RST and the RTD channels group into one cluster; also well supported by a dark patch in fig 42. 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 jI  L    ROF1RS1 , .__L, ...... ROM. LTDS LST1 LO1 RTD4 RST1 ROF1 LTDSN LST1 ... ....  H LOF1 .. ,.... RTD~ ___ (1) RST1CL 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 42. Seizure 11 of patient P093: Clustersimilarity 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 II 0 0.1 0.2 0.3 0.4 Fusion level 0.5 0.6 Figure 43. 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 preseizure period. The overall cluster configuration is listed in table 41. Table 41. Spatiotemporal 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 timeintervals 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 bilateral homologous connection, as seen from all the matrices in fig 42. 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 lefthemisphere channels, as seen in Fig 42. 5. Interestingly, no temporal changes are seen in the spatialpatterns yet. In Fig 43. 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 spatialcorrelations 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 spatiotemporal groupings from a collection of multivariate timeseries is relatively easier to demonstrate in linear coupling cases. However, non linear dynamic model constructions are extremely hard and mostly, nontrivial. We therefore decided to pursue a posterior verification of the timeaveraged cluster groupings on the intracranial EEG data, using the quasisurrogate analysis [40, 6668, 73, 9495], technique. Recall that the cluster groupings obtained over 30 minute timesegments involve two steps. First step consists of applying spectral clustering technique on the SOMsimilarity indices (computed on 10 second intracranial EEG data segments). Then similar grouping patterns among channels are extracted by using hierarchical clustering approach on the clustersimilarity matrices P. Specifically, 91 10second windows (spanning 30 minutes, 3 windows per minute) are analyzed in each pass to obtain longterm groupings among channels. In order to validate this 2step approach, we define our hypothesis as follows Ho: The average withincluster channel interaction at each window (out of 91 10second windows) is not significantly different from the corresponding betweencluster 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. Withincluster interaction is computed by averaging the pairwise similarity indices for all the channels within a cluster. For betweencluster interaction, the pairwise interactions between 3 channels, picked randomly from each of the 3 clusters are computed. A betweencluster 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 quasinormal distribution, implying that the withincluster interaction value can now be compared with the mean and the variance sample estimates of the betweencluster 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 withincluster interaction at time '', for cluster 'i'. Cb ) is the mean and c(Cb,) is the standarddeviation of the betweencluster interaction at time t'. Z, reflects the zscore and is considered significant at the 95 percentile significance if Z' > 1.96 (reject Ho). In table 42 the bolded value in each cell represents the number of windows (out of 91) having significant zscores in the 30 minute period corresponding to fig 42 (P093, Seizure 11). It is easy to observe that the nullhypothesis HO is rejected beyond doubt, validating the clustering results. Table 42. P093, Seizure 11: Over each 30 minute (91 samples total) window, number of times the withincluster interaction is greater than betweencluster 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: Spatiotemporal clustering analyses, similar to the one described on seizure 11 were performed on several other seizures, of the same patient P093. The clustersimilarity matrices P obtained from timeintervals surrounding seizures 4 & 5 and 6 & 7 of patient P093 are shown in fig 45. & fig 46. respectively. Channel groupings for the same are listed in tables 43 & 44. 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 44. 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 betweencluster interactions. Cluster veracity can be visually verified by observing that amplitudes representing withincluster interaction for cluster profiles are mostly higher that the amplitudes representing betweencluster interaction for surrogate profiles, at each timeinstance. 1. Consistent to the observation in seizure 11, we observe the temporal depth and the subcortical regions of the left hemisphere are always bunched together. 2. Once again, the association of ROFLOF areas into same cluster suggests a strong homologous connection between the orbitofrontal 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 predefined 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 45 & fig 46 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 EEGSOM 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 47. shows the dendrograms created for seizure segments 2 hours prior to seizure 1 and 30 minutes preseizure, respectively. As before, the number of clusters (nl) specified in the spectralclustering step after SOMSI 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 43. Spatiotemporal 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 PostSeizure 5, (30 1 hr) RTD LTD, LST, RST LOF, ROF Table 44. Spatiotemporal 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 45. Seizures 4 and 5 of patient P093: Clustersimilarity 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 aa LTD3 LST1 LOF1 RTD4 RST1 ROF1 30 to60 ' m6 uZ Sz6 Sz7 Figure 46. Seizures 6 and 7 of patient P093: Clustersimilarity 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 47. Dendrograms corresponding to P092, Seizure 1. Top: 2 hours before Seizure Bottom: 30 minutes preseizure. 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. 47) 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 subclusters, a large and a small cluster. The small cluster consists of only two channels, LTD (3 & 5) and fuses with the other subcluster at a very high fusion level (implying weak link). If n2 was to be increased to 4, the clustering algorithm would classify this subcluster as an independent cluster. A detailed analysis on all seizures in P092 revealed a strong intrachannel correlation (or low fusion level) between channels LTD (3 & 5) and a weak interchannel 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 subclusters that have a strong withinsubcluster interaction, but a weak betweensubcluster interaction. Consequently, the withincluster interaction of the largest cluster can be expected to be as weak as or marginally better than the betweencluster 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 7I71 L I U' LISl' LGI H I 2 H1yl II1"G LIU' LYI' LGI HIU2 UY", 14GI1 LI ' Lll. LCIr ~ 141.u r. ( ) 1 U41' Lr  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  rl 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 48. Clustersimilarity matrices indicating the probability that two channels share the same cluster label in a 30 minute timeinterval, 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 spatiotemporal clustering results for seizures 1, 3, 4, 5, 6 & 7 are summarized in tables 45 to 410. Table 45. Spatiotemporal 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 PostSeizure, (0 30 mins) RTD, RST LTD, LST LOF, ROF PostSeizure, (30 1 hr) RTD LTD, LST, ROF LOF, RST Table 46. Spatiotemporal 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 PostSeizure, (0 30 mins) RTD, RST LST, LTD LOF, ROF PostSeizure, (30 1 hr) RTD, RST LST, LTD LOF, ROF Table 47. Spatiotemporal 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 PostSeizure, (0 30 mins) RTD, RST LST, LTD LOF, ROF PostSeizure, (30 1 hr) RTD, RST LST, LTD LOF, ROF Table 48. Spatiotemporal 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 PostSeizure, (0 30 mins) RTD, RST LST, LTD LOF, ROF PostSeizure, (30 1 hr) RTD, RST LST, LTD LOF, ROF 