# **IDENTIFYING THE EPILEPTIC NETWORK**

**Topic Editors Mark Holmes and Don Tucker**

#### *FRONTIERS COPYRIGHT STATEMENT*

© Copyright 2007-2014 Frontiers Media SA. All rights reserved.

All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.

The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.

Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.

Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.

As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.

All copyright, and all rights therein, are protected by national and international copyright laws.

The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use.

**ISSN** 1664-8714 **ISBN** 978-2-88919-296-0 **DOI** 10.3389/978-2-88919-296-0

### *ABOUT FRONTIERS*

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

### *FRONTIERS JOURNAL SERIES*

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing.

All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

### *DEDICATION TO QUALITY*

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view.

By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

### *WHAT ARE FRONTIERS RESEARCH TOPICS?*

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area!

Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

## **IDENTIFYING THE EPILEPTIC NETWORK**

Topic Editors: **Mark Holmes,** University of Washington, USA **Don Tucker,** Electrical Geodesics, Inc. and the University of Oregon, USA

Image taken from: Song J, Tucker DM, Gilbert T, Hou J, Mattson C, Luu P and Holmes MD (2013) Methods for examining electrophysiological coherence in epileptic networks. Front. Neurol. 4:55. doi: 10.3389/fneur.2013.000550055

# Table of Contents

### *04 Identifying the Epileptic Network* Mark D. Holmes and Don M. Tucker


Robert Todd Constable, Dustin Scheinost, Emily S. Finn, Xilin Shen, Michelle Hampson, F. Scott Winstanley, Dennis D. Spencer and Xenophon Papademetris

*45 Local Functional Connectivity as a Pre-Surgical Tool for Seizure Focus Identification in Non-Lesion, Focal Epilepsy*

K. E. Weaver, W. A. Chaovalitwongse, E. J. Novotny, A. Poliakov, T. G. Grabowski and J. G. Ojemann

*59 Computer-Aided Diagnosis and Localization of Lateralized Temporal Lobe Epilepsy Using Interictal FDG-PET*

Wesley T. Kerr, Stefan T. Nguyen, Andrew Y. Cho, Edward P. Lau, Daniel H. Silverman, Pamela K. Douglas, Navya M. Reddy, Ariana Anderson, Jennifer Bramen, Noriko Salamon, John M. Stern and Mark S. Cohen

*73 Focal Peak Activities in Spread of Interictal-Ictal Discharges in Epilepsy With Beamformer MEG: Evidence for an Epileptic Network?*

Douglas F. Rose, Hisako Fujiwara, Katherine Holland-Bouley, Hansel M. Greiner, Todd Arthur and Francesco T. Mangano


Ceon Ramon and Mark D. Holmes

## Identifying the epileptic network

#### *Mark D. Holmes1 \* and Don M. Tucker <sup>2</sup>*

*<sup>1</sup> Neurology/Regional Epilepsy Center, Harborview Medical Center, University of Washington, Seattle, WA, USA*

*<sup>2</sup> Electrical Geodesics Inc., Eugene, OR, USA*

*\*Correspondence: mdholmes@u.washington.edu*

#### *Edited by:*

*Jorge Asconape, Loyola University, USA*

Progress in characterizing the functional networks of the normal human brain is now rapid, with evidence from both regional correlational patterns from functional MRI and fiber tractography from diffusion MRI. Increasingly, the tools of cerebral network analysis are being applied to understand the derangement of specific cortical and subcortical networks in epileptic disorders. In this approach, the clinical manifestations of epilepsy are viewed as the consequence of the pathologies of network dynamics and functional connectivity that may involve abnormal network pathways. Importantly, concepts of epileptic networks are supplanting the older, and more simplistic, notion that epileptic seizures must be either "focal" (or partial) or "generalized" in nature. Rather, seizures can be understood to result from the paroxysmal and pathological activation of specific neuronal connections. The characteristics of these may not fit with conventional assumptions, and could include widespread and bilateral involvement during seizures which classically are considered as focal, or could involve restricted cortical/subcortical regions during some seizures that are typically considered as generalized in nature. We believe that identifying patient-specific epileptic networks will provide critical insights into epilepsy syndromes, and more importantly, these insights will lead the way to novel forms of treatment for affected individuals.

Technological improvements in several fields have contributed to the tools applied to understanding epileptic networks, particularly in neuroimaging (MRI, FDG-PET, fMRI), and in electromagnetic recordings (dense array EEG, MEG). Investigators are also finding that combining these methodologies may have a synergistic effect in regard to enhancing our understanding of the involved cortical networks. In this volume we have assembled contributions from an international group of investigators, each of whom has approached the problem of identifying the epileptic network from somewhat different perspectives. The unifying theme in all cases is the question of how the application of a specific technology, or a simultaneous combination of technologies, may enhance our insight into the recognition of the epileptogenic zone in the resting state.

This book opens with a chapter by Stefan and Lopes da Silva (1), who review the evidence for the concept of epileptic networks. These authors discuss the structure and dynamics of cortical networks, describe how these connections can be analyzed through linear and non-linear methodologies, and outline the dynamics of neuronal networks in the context of combined EEG/MEG and EEG/fMRI signals analysis. They conclude that the resulting network analysis has clear relevance to understanding the nature of seizures occurring with focal cortical dysplasia and with temporal lobe epilepsy. Remarkably, they suggest that an absence seizure, often considered the prototypical generalized seizure, is actually a fast-spreading localized event.

In a similar vein, Leite et al. (2) propose a novel method for linkage of EEG and fMRI signals in network analysis by describing in their report a "transfer function" between these divergent measures. They perform independent component analysis of EEG and extract metrics that express models of EEG-fMRI function from resulting time courses. These metrics are then used to predict fMRI activity and thus the brain regions associated with epileptic activity. The authors illustrate the methodology in a proof of concept report on the application of this function to fMRI-EEG data obtained during both ictal and interictal states in one subject with a hypothalamic hamartoma.

In the next two chapters, by Constable et al. (3) and Weaver et al. (4) the focus is on using resting state fMRI to assess functional connectivity in the human brain, and how this approach can be applied to epilepsy. These two groups describe the functional reorganization that occurs in epilepsy, and the potential that connectivity measures have in identifying a network of seizure-generating tissues. Both groups stress the importance of focal connectivity measures as adjunctive tools in the identification of the epileptogenic zone in patients with refractory epilepsy who are being considered for resective surgery.

On the other hand, Kerr et al. (5) find that the interictal FDG-PET, by visualization of the metabolic changes that take place across the whole brain in epilepsy patients, offer another method to observe abnormal brain networks in the resting state. These authors report that in temporal lobe epilepsy, examination of patterns of metabolic dysfunction may assist in lateralizing the onset of seizures. They report on the development of a computerized assisted diagnostic tool for implementing the metabolic analysis in clinical practice.

Rose et al. (6) studied simultaneous MEG-EEG activity in a series of children with refractory epilepsy. They studied the MEG signals throughout the brain using a beamformer algorithm, and they determined virtual MEG spike locations with a spike detection program. Comparisons of the MEG results with intracranial EEG recordings were conducted both for EEG spikes and for the onset and spread of seizures. By demonstrating similarities with the invasive electrographic findings, the authors conclude that the pattern of interictal MEG findings has the potential to define the distribution of the epileptic network, thereby providing a noninvasive method to analyze abnormal neuronal connections.

Yamazaki et al. (7) have pioneered the ability to simultaneously record 256 channel dense EEG (dEEG) and invasive subdural EEG recordings in temporal lobe epilepsy, thus helping to establish the validity of dEEG recordings. In their chapter in this volume, Yamazaki et al. (7) extend this work to cases of neocortical epilepsy by demonstrating that dEEG, by covering the whole head with sufficient sensor density, can reliably localize epileptiform discharges when compared to invasive studies.

The final two chapters concern the application of analytic techniques to examine abnormal synchronization of the interictal dEEG data to establish the presumptive epileptogenic zone. Song et al. (8) discuss the use of coherence measures in the examination of interictal spikes to determine the extent and distribution of epileptic networks. In their contribution, Ramon and Holmes (9) provide evidence that brief segments of interictal dEEG, free of classical epileptiform patterns, nevertheless may contain stable markers that reveal the likely epileptic network. These markers are identified through analysis of localized patterns of phase synchronization and cross-frequency coupling that appear specific to the epileptogenic region as proven by later intracranial recordings.

The topics covered in this volume present an introduction to the study of identifying epileptic networks. They are only a sample of the many current approaches to cerebral network analysis that could be applied to epilepsy. Nevertheless, we are hopeful that the material presented here will provide encouragement for additional work to clarify – and treat – the pathological dynamics of human cerebral networks in epilepsy.

#### **References**


#### *Received: 06 June 2013; accepted: 17 June 2013; published online: 01 July 2013. Citation: Holmes MD and Tucker DM (2013) Identifying the epileptic network. Front. Neurol. 4:84. doi: 10.3389/fneur.2013.00084*

*This article was submitted to Frontiers in Epilepsy, a specialty of Frontiers in Neurology. Copyright © 2013 Holmes and Tucker. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*

## Epileptic neuronal networks: methods of identification and clinical relevance

#### **Hermann Stefan<sup>1</sup>\* and Fernando H. Lopes da Silva2,3**

<sup>1</sup> Department of Neurology, University Hospital Erlangen, Erlangen, Bavaria, Germany

<sup>2</sup> Centre of Neuroscience, Swammerdam Institute for Life Sciences, University of Amsterdam, Amsterdam, Netherlands

<sup>3</sup> Department of Bioengineering, Instituto Superior Técnico, Lisbon Technical University, Lisbon, Portugal

#### **Edited by:**

Mark Holmes, University of Washington, USA

#### **Reviewed by:**

Andreas Schulze-Bonhage, University Hospital Freiburg, Germany Don Tucker, University of Oregon, USA

Marino M. Bianchin, Universidade Federal do Rio Grande do Sul, Brazil

#### **\*Correspondence:**

Hermann Stefan, Neurological Clinic, Friedrich-Alexander University Hospital Erlangen, 10 Schwabachanlage, 91054 Erlangen, Bavaria, Germany. e-mail: hermann.stefan@ uk-erlangen.de

The main objective of this paper is to examine evidence for the concept that epileptic activity should be envisaged in terms of functional connectivity and dynamics of neuronal networks. Basic concepts regarding structure and dynamics of neuronal networks are briefly described. Particular attention is given to approaches that are derived, or related, to the concept of causality, as formulated by Granger. Linear and non-linear methodologies aiming at characterizing the dynamics of neuronal networks applied to EEG/MEG and combined EEG/fMRI signals in epilepsy are critically reviewed. The relevance of functional dynamical analysis of neuronal networks with respect to clinical queries in focal cortical dysplasias, temporal lobe epilepsies, and "generalized" epilepsies is emphasized. In the light of the concepts of epileptic neuronal networks, and recent experimental findings, the dichotomic classification in focal and generalized epilepsy is re-evaluated. It is proposed that so-called "generalized epilepsies," such as absence seizures, are actually fast spreading epilepsies, the onset of which can be tracked down to particular neuronal networks using appropriate network analysis. Finally new approaches to delineate epileptogenic networks are discussed.

**Keywords: epileptic networks, neurophysiological classification, MEG/EEG, basic concepts, clinical approaches**

#### **EPILEPTIC NEURONAL NETWORKS: BASIC CONCEPTS, STRUCTURE, AND DYNAMICS**

Seminal descriptions of neuronal networks in which neurons are the elementary units that transmit signals through synaptic contacts were performed by Ramón y Cajal (1894). The concept of neuronal networks has occupied a prominent role in the Neurosciences since. Research into how neuronal networks are interconnected forming the wiring structure of the brain has been a constant thread along the years. An important question has been to find rules that link structural connectivity of neuronal networks with information flow and processing in such networks. A model of how information may flow in cortical networks was proposed by Abeles (1991), Abeles and Gerstein (1998) who introduced the concept of "synfire chains" meaning synchronous working chains of neurons in a network, i.e., sets of interconnected neurons that participate in common tasks. He elaborated this concept further in what he called Corticonics where insights from anatomical and physiological studies are combined with mathematical and computer modeling to obtain quantitative descriptions of cortical functions. These notions have been explored in modern neural network modeling. At the level of the organization of the whole brain the network concept has been extended, among others by Mesulam (1990), describing "local networks" (engaged in modality-specific processing such as analysis of shape, spatial location, and object identification in the visual modality) and "large-scale networks" that incorporate numerous parallel lines of communication with multiple cross-links, enabling integrative processing. The complexity of the organization of these networks

of the brain has been compared with that of other large-scale networks, such as the World-Wide Web, the Internet, social networks, or metabolic networks (Jeong et al., 2000), and has been the object of similar mathematical analyses based on topological properties; graph analysis is an example of this approach. In this way the notions used in these mathematical analyses have been adopted in the description of neuronal networks of the brain, and terms as "nodes" and "hubs" have entered the field of the neurosciences. Thus concepts from graph theory are being used to represent neuronal networks: a neuron is denominated a node and a neuronal network consists of nodes connected by links or edges; highly connected nodes are called hubs and an uninterrupted sequence of links forms a path; questions such as how the flow of information takes place from a node to another, can be analyzed by finding the possible paths in a graph. The structure of these networks deviates from random; this structure has some properties of "small-world" networks where any two nodes may be connected by short paths and where a few "hubs" may dominate the whole connectivity of the network (Watts and Strogatz, 1998). In such a structure the number of links of a node mayfollow a power-law which characterizes the so-called "scale-free" networks which are heterogeneous. Small world networks are hypothesized to optimize rapid synchronization transfer creating a balance between local processing and global integration (Meador, 2011). We should note that in current applications of these concepts in the field of neurophysiology, the "nodes" correspond simply to the sites where signals are recorded, be it using EEG, MEG, or functional magnetic resonance imaging (fMRI), and not to the neuronal elements as such. This implies

that there is an enormous distance between those "nodes" and the neuronal reality.

A particular feature of some neuronal networks is that these are interconnected by means of re-entrant connections, i.e., that some nodes tend to receive connections from other nodes to which they project by relatively short paths, some of which have been well characterized both anatomically and physiologically such as the cortico-thalamic-cortical system (Steriade, 2001), and the entorhinal – hippocampal-entorhinal system (Kloosterman et al., 2004). A few hubs may dominate the whole connectivity of the network, and if these hubs would represent neuronal features with a high degree of excitability, the latter may be rapidly distributed throughout the whole network. Furthermore if the network possesses re-entrant properties an even larger network may display this high degree of excitability with complex dynamics.

Whereas the historical approach to understand higher level operating principles in the brain was to consider it subdivided into anatomical regions with local functional properties, the current approach, inspired by the theoretical analysis of complex networks, as described above, is to emphasize networks interactions and connectivity at short and long range.

These theoretical considerations provide a convenient approach to better understand the pathophysiology of epilepsies. In this context, however, it is essential to go forward from the description of network connectivity and structural properties, as presented above, to the dynamic dimension, i.e., to the study of the activity of the networks as function of time. A fundamental characteristic of these neuronal networks is that their dynamics are essentially non-linear given the non-linear transfer properties of neuronal elements.

The dynamics of the population of neurons that constitute the neuronal networks can be considered at different spatial scales: microscopic, mesoscopic, and macroscopic. The last two levels are particularly relevant with respect to the dynamics of EEG or MEG signals in general, and in the case of epileptic activity in particular. A variety of molecular processes at the microscopic (cellular) level may lead to changes in the stability of neuronal networks causing epileptiform seizures, that become manifest at the mesoscopic and macroscopic levels. A generalized concept is that seizures occur in strongly coupled neuronal networks due to a shift in the dynamical balance between excitatory and inhibitory processes with a predominance of the former. In terms of the mathematical theory of complex non-linear systems we may state that such networks display bistability, i.e., they may feature two stable operational states that may exist simultaneously for the same set of system parameters. One of these states is the normal, inter-ictal state, and the other is the epileptic or ictal state of the network. The transition between the two states is called a dynamical bifurcation (Lopes da Silva et al., 2003). In epileptic brain certain networks have abnormal parameters at the molecular and cellular levels, due to genetic or to acquired pathogenic factors, rendering some essential parameters, that control network stability, extremely vulnerable to the influence of exogenous and endogenous factors, such that this kind of bifurcations may occur easily. In this way abnormal oscillations and other events, such as epileptiform spikes, may occur in hubs of these neuronal networks with abnormal parameters.

At the local neuronal network level, some hubs constituted by neurons and associated glia constitute oscillatory systems that became increasingly coupled at the transition to a seizure, thereby recruiting more distant neuronal networks, constituting complex oscillatory circuits, which can be recognized by EEG or MEG recordings (Zhang et al., 2011). Accordingly, circuits of this kind have been described in several forms of epilepsy, such as in the thalamocortical system involved in Absence epilepsies (Meeren et al., 2002, 2005; Suffczynski et al., 2004), and also in several other forms of epilepsy as discussed by Halasz (2010) for rolandic epilepsy (inner part of sylvian fissure), Landau Kleffner syndrome (perisylvian opercular structure and/or posterior part of first temporal convolution), electrical status epilepticus in sleep (perisylvian area, bilateral widespread involvement of cortical mantle, thalamic mediodorsal nucleus), Lennox–Gastaut syndrome (diffuse bi-synchronous epileptogenic system and cortical excitation with augmented cortico-thalamic oscillations), nocturnal frontal lobe epilepsy (frontal medial and orbital surfaces). Also McIntyre and Gilby (2008) described in various models of temporal lobe epilepsy the recruitment of the parahippocampal cortices including piriform, perirhinal, and entorhinal cortex in addition to the hippocampus proper. Along the same line of thought Spencer (2002) put forward the concept of human epilepsy as a disorder of large neural networks, and Avanzini et al. (2012) proposed the term "system epilepsies" to describe some types of epilepsy that depend on the dysfunction of specific functional neural systems. Clinical and network analytical studies are required to advance detection of such dysfunctional specific systems, and characterize more precisely their abnormal structure and dynamics.

#### **EPILEPTIC NEURONAL NETWORKS: CAUSALITY, LINEAR AND NON-LINEAR METHODS, AND NEW APPROACHES**

Epileptic conditions have to be characterized on the basis of clinical evidence, but a comprehensive analysis of the brain systems responsible for epileptic manifestations resorts to neuroimaging techniques that may reveal structural abnormalities, and to the EEG, MEG, and fMRI which can reveal the underlying dynamics. Here we concentrate on general aspects of methodologies aiming at characterizing epileptogenic networks. This implies *functional connectivity mapping* to determine the dynamics of epileptiform activities displayed as patterns of interactions between anatomically connected neural nodes responsible for these abnormal activities. Current methodologies allow a direct evaluation of correlations between EEG seizure activities, their propagation dynamics on the one hand, and the evolution of clinical signs on the other, observed using combined EEG-video recording. This is nicely illustrated in the case of patients with intractable Jacksonian seizures in whom intra-cranial EEG recordings (iEEG) were made in order to assess the indication of surgery (Akiyama et al., 2011). During the Jacksonian seizures High Frequency Oscillations (HFOs > 40 Hz) started in the sensory cortex and propagated to the motor cortex when ictal motor signs occurred, but as the seizure progressed ictal HFO spread or reverberated into the rolandic region; further when the seizure became secondarily generalized the ictal HFOs were limited to the Rolandic region (**Figure 1**).

Considering that the main objective of this paper is to examine evidence for the concept that epileptic activity should be envisaged in terms of functional connectivity and dynamics of neuronal networks, we emphasize here approaches that are derived, or related to the concept of causality, as formulated by Granger (1998) in econometrics. In short, according to Ganger causality an observed time-series *x*(*t*) can be considered the cause of another series *y*(*t*) if knowledge of the past values of *x*(*t*) improves the prediction of *y*(*t*).

The"directed transfer function" (DTF) extends Granger causality to multichannel EEG/MEG data (Kaminski and Blinowska, 1991) and has been applied to estimate functional connectivity in epilepsy (Franaszczuk and Bergey, 1998) and more recently by Dai et al. (2012) as discussed more in detail below. It should be noted, however, that DTF represents a linear combination of causal relations, not only along direct pathways, but also along indirect pathways. This led to the development of another measure, the

so-called "direct DTF" (dDTF), which emphasizes direct associations over indirect ones (Korzeniewska et al., 2003). The latter, however, has the limitation that it needs relatively long signal epochs to be estimated reliably. To minimize the effects of nonstationary behavior of EEG/MEG signals, several methods have been developed, among which the method of short-time DTF (SDTF). By combining dDTF and SDTF new measures were proposed: the SdDTF which estimates direct causal influences between signals, not mediated by other signals, in short-time epochs, and the Event-related causality (ERC) which estimates event-related changes in SdDTF (Korzeniewska et al., 2008). This has been applied mainly in the analysis of task-related changes in EEG or MEG signals, particularly in the high frequency gamma range during cognitive tasks.

It should be added that the methodologies described above are based on the assumption that the relationships between EEG/MEG signals are linear; this may be an acceptable approximation in many cases, although during epileptic seizures it is doubtful whether the linear assumption always holds. Therefore non-linear methods were developed with the objective of estimating the coupling between different EEG/MEG signals in general. A first group of these methods was based on the estimation of mutual information (Mars and Lopes da Silva, 1983) and on non-linear regression (Lopes da Silva et al., 1993; Wendling et al., 2001; Kalitzin et al., 2007) applied to EEG or MEG signals. A second group of methods, was based on tools imported from the field of non-linear dynamical systems and chaos theory (Lehnertz, 1999; Iasemidis, 2003). Related to this class of methods two other may be distinguished, namely: phase synchronization (PS) methods (Bhattacharya et al., 2001; Rosenblum et al., 2004), generalized synchronization (GS) methods (Arnhold et al., 1999; Stam and van Dijk, 2002), and more recently directed phase lag index (dPLI; Stam and van Straaten, 2012). The former estimate the instantaneous phase of each signal and then compute a quantity based on co-variation of extracted phases to determine the degree of coupling between signals. GS methods also consist of two steps: the reconstruction of state space trajectories from time-series signals and the computation of a similarity index on reconstructed trajectories. The dPLI characterizes spatial and temporal patterns of phase relations in functional brain networks.

A study of Wendling et al. (2009) is particularly interesting because these authors compared directly several methods of estimating functional connectivity between EEG/MEG signals, namely (a) linear (cross-correlation, cross-spectral analysis – coherence and phase), (b) non-linear regression (mutual information, h2 association index), (c) PS, and (d) GS, applied to a well defined data set. To make the comparison these authors built computer models of interconnected neuronal networks with defined coupling parameters, that can generate oscillatory activity typical of epileptic seizures. This model-based methodology allows establishing at will the degree of coupling between the different neuronal networks that generate the EEG signals. In this way the coupling between the signals of different networks can be estimated and the computed values obtained using different methods can be directly compared among themselves, and with the values of the coupling parameters established *a priori*. This comparison revealed that there was no "ideal" method, i.e., none of the methods performed better than all the other ones in all nine studied situations. Nevertheless, regression methods (linear or non-linear) showed sensitivity to the coupling parameter in all tested models with average or good performance, what leads to the conclusion that these are robust, and it is advisable to first apply these regression methods in order to characterize functional brain connectivity, under normal or pathologic conditions. In any case it is useful in practice to compare the results of different measures to get more reliable estimates of the coupling of interest. **Figure 2** shows examples of the results of the application of nonlinear regression analysis to intra-cerebral EEG signals identifying network associations around the moment of seizure onset.

In the last decade new techniques have entered the field based on the application of MRI: namely fMRI, particularly in conjunction with EEG, and diffusion-based tractography imaging (DTI). Regarding fMRI, and in the context of determining pathways of propagation of epileptic activity in neuronal networks, Dynamic Causal Modeling (DCM) applied to the interpretation of hemodynamic signals (BOLD) is being extensively used to determine the patterns of interaction between different neuronal networks (Friston et al., 2003).

In an animal experimental model of absence epilepsy, this methodology has been integrated with associated EEG signals (David et al., 2008). In this case the performances of DCM and Granger causality estimates were compared, showing approximately similar results. In human epilepsy these methodologies were applied recently in epileptic patients with Hypothalamic Hamartomas and were able to yield plausible estimates of seizure propagation pathways (Murta et al., 2012). DTI is based on the principle of anisotropic diffusion of water molecules in white matter tracts throughout brain tissue. In a study of children with temporal lobe epilepsy displaying spikes over the Rolandic region identified in the MEG, the hypothesis that the latter occurred due to activity propagating along neural aberrant pathways connecting the temporal lobe and the Rolandic cortex appeared plausible according to the DTI analysis (Bhardwaj et al., 2010).

#### **EPILEPTIC NEURONAL NETWORKS: CLINICAL QUERIES AND PRACTICAL RELEVANCE**

#### **NETWORKS IN FOCAL CORTICAL DYSPLASIAS AND OTHER LESION-RELATED EPILEPSIES**

To put in evidence functional dynamics of neuronal networks engaged in epileptic seizure activity the study of Focal Cortical Dysplasias (FCD), Dysembryoplastic NeuroEpilethelial Tumors (DNET), and Periventricular Nodular Heterotopias (PNH), which frequently are associated with pharmaco-resistant epilepsy, is particularly enlightening. These lesions may be synaptically connected with other neuronal networks, such that the epileptic activity may propagate along the connecting pathways constituting an "epileptogenic network" (Aubert et al., 2009), or otherwise may stay confined to the region of, and around, the FCD lesion. Interestingly, anatomical alterations in tissue microstructure adjacent to some FCDs were detected using DTI-MR imaging. These overlapped with the localization of clusters of equivalent dipoles of epileptiform spikes (Widjaja et al., 2009).

Therefore it is most relevant to determine the functional organization of these epileptogenic networks, since this may give useful indications for a possible surgical intervention and the corresponding prognosis. Different methods have been applied to estimate the functional connectivity of neuronal networks in these cases. Using depth EEG registrations (stereoencephalography) functional analysis of multiple EEG signals was performed using non-linear regression (h2 association index) by Valton et al. (2008), and by computing the so-called "Epileptogenicity index" (Aubert et al., 2009). In the former study Valton et al. (2008) applied this methodology to analyze depth seizure EEG signals in a patient with bilateral PNH. Non-linear regression analysis revealed a large epileptogenic network extending beyond the PNH and involving remote cortical structures. The fact that this is a vast epileptogenic network may account for surgical failures in patients with this kind of heterotopias.

Aubert et al. (2009) used depth electrodes (stereoencephalography) and the "Epileptogenicity index" from Bartolomei et al. (2008). This index is based on the spectral content of

high frequency components [beta: (12.4–24 Hz) + gamma: (24– 90 Hz)] relative to lower frequency components. It accounts for both the propensity of a brain area to generate high frequency oscillations and the time for this area to get involved in the seizure. In this study depth EEG signals of several cortical areas were analyzed; the authors found in one group of patients (31% of the patients) one strictly localized epileptogenic zone, while a second group (61%) displayed more than one epileptogenic network. The success of surgery was related to the extent of the epileptogenic network. This implies that it is important to determine with the highest precision the extent of the epileptogenic network while planning a surgical intervention.

This study illustrates the clinical relevance of making a functional analysis of the neuronal networks associated with an epileptogenic lesion using quantitative methods of EEG signal analysis, and that the "epileptogenicity index" may be a useful tool even in cases where a relatively small number of EEG signals are available. In the case of Mesial Temporal Lobe Epilepsies (MTLE) extended networks outside the temporal lobe may be involved what led Ryvlin and Kahane (2005) to coin the denomination "Temporal-plus Epilepsy" where epileptiform activity appears in multiple brain lobes in addition to the temporal lobe.

The importance of differentiated "primary and secondary" inter-ictal spike activity for the identification of ictal and

propagated epileptic activity was pointed out by Badier and Chauvel (1995).

The determination of the extent of an epileptogenic network may also be obtained from recordings of inter-ictal spikes, either using EEG or MEG, since it is known that the propagation patterns of spikes may give valuable information about neural networks associated with epilepsy (Spencer, 2002) and to the outcome of epilepsy surgery (Alarcon et al., 1997; Hufnagel et al., 2000; Schulz et al., 2000). Indeed inter-ictal EEG (reviewed in Rodin et al., 2009) and MEG spikes (Tanaka et al., 2010) may be considered biomarkers of epileptogenic networks; more recently, in the same context, a novel kind of biomarker – High Frequency Oscillations, i.e., HFOs (ripples and fast ripples), – has raised a lot of interest and its role is the object of investigation, reviewed by Jacobs et al. (2012). In order to perform this kind of analyses using scalp EEG and/or MEG recordings it is necessary to combine information about electric and/or magnetic fields with anatomical information obtained from MRI. In this way spatio-temporal analysis of EEG/MEG can yield useful information about the propagation of inter-ictal spikes in epileptogenic networks. The validation of this approach was systematically assessed in the study of Tanaka et al. (2010) who compared the results of the spatio-temporal source analysis of EEG and MEG of inter-ictal spikes with data obtained by way of iEEG recordings made extraoperatively in 10 patients.

This study showed that the analysis based on MEG spikes yield a very similar propagation pattern as observed in iEEG, better than EEG data. It should be noted, however, that the MEG signals were obtained from 102 sites whereas for the EEG only 70 recording sites were used, what may account, at least partially, for the difference in performance between the two sets of recordings.

A clinical example where the value of inter-ictal MEG spikes is illustrated is the case study of a patient with PNH where Magnetic Source Imaging (MSI), i.e., a combination of MEG and co-registered MRI, was applied, and used to guide intraoperative electrocorticography. Thereafter intra-cerebral depth electrodes and subdural strips were implanted guided by the MSI data, which revealed two separated zones of spike activity. A cluster analysis of electrographic recordings of spikes showed two clusters based on source localization using equivalent dipole source models and a realistic volume conductor model; the main cluster arose from the temporal neocortex and not from the PNH area. This led to the assumption that the networks of the temporal cortex distant from the lesion might act as primary epileptogenic network. Accordingly a resection of the temporo-occipital neocortical tissue including the main spike focal area was performed, with excellent postoperative seizure control (**Figures 3A,B**; Stefan et al., 2007). This led to the definition of the epileptogenic network, including the heterotopia and overlying neocortex, what was essential to assure a positive surgical outcome. Another study (Dai et al., 2012) used also inter-ictal Magnetoencephalography (MEG) to identify epileptogenic networks; these authors estimated cortical sources of spikes using realistic head models, computed the directional connectivity between those sources in patients with pharmaco-resistant epilepsies (different kind of pathologies)

who were pre-surgically investigated using intra-cranial electrodes under the guidance of the MEG data, and whose epileptogenic zone was later removed. They found a good overlap between the primary sources of epileptiform spikes and the epileptogenic zones that were later surgically resected, in line with the results obtained by Tang et al. (2003), Shiraish et al. (2005), Stufflebeam et al. (2009), and Wang et al. (2012). In the context of the present discussion on the importance of network analysis in epilepsy, the study of Dai et al. (2012) goes further than previous investigations because these authors applied the DTF, introduced above, using an open-source software package (eConnectome). In this way they estimated the direction of the propagation of the epileptiform activity along the interconnected neuronal networks as illustrated in **Figure 4**.

#### **NETWORKS IN GENERALIZED EPILEPSIES**

In the light of the concepts of epileptic neuronal networks and recent experimental findings the dichotomic classification in focal and generalized epilepsy has to be re-evaluated. Indeed it is necessary to reassess the role of epileptic networks in the so-called generalized epilepsies.

#### **Are generalized epilepsies actually fast spreading focal epilepsies?**

The concept of "primarily generalized epilepsy," as for example in Childhood Absence Epilepsy (CAE), implies that all brain regions simultaneously would generate spike-and-wave discharges (SWDs) and the associated seizure, classically considered to be triggered by some central process associated with the diffuse cortico-thalamic system according to the "centrencephalic" concept of Jasper and Penfield (1954), or by the interplay between

**FIGURE 3 | (A)** (a) Streamtube visualization of the right optic radiation based on diffusion tensor imaging. (b) For navigation, a three-dimensional object representing the optic radiation (wrapping the individual fibers) and two distinct MSI foci (red) are generated. (c) Relation of optic radiation (visualized as streamlines) to MSI foci. (d–f) Sagittal/coronal/axial view of T1-weighted images with registered DTI and MSI data. Localization of focal epileptic activity is below the optic Tract (Stefan et al., 2007).

(Continued)

**FIGURE 3 | (B)** (a) Intraoperative electroencephalographic recordings with platinum-electrodes electrode position confirmed by intraoperative T1- and T2-weighted high-field-MR imaging. (b) MSI guided electrode implantation of intra-cerebral depth and subdural electrodes; spike activity in lateral cortex and periventricular heterotopia, the

polyspikes during intraoperative ECoG. (c) Intraoperative MR-imaging after cortical resection of MSI-focus No. II (with platinum-electrodes still in situ).

thalamus and cortex as in the corticoreticular hypothesis of Gloor (1968). There is growing evidence, however, that this is a too simplified view, and that in these cases there is a cortical frontal neuronal network where the onset of the SWDs is located. Thus the question may be formulated in a simple way: are SWDs manifestations of a generalized or of a focal process? In the same vein

a critical analysis of the classification of epilepsies recommend "abandoning these terms as overall classification categories into which all epilepsies must fit" (Berg and Scheffer, 2011).

One important hindrance in solving this controversy is that many studies have used inadequate methodologies. Most of the classic observations of "generalized SWDs" (GSWDs) were based on visual inspection of EEG recordings on paper at relatively high speed, as even reported recently (Koutroumanidis et al., 2012). Although the human condition of CAE may differ in some respect from genetic rodent models of Absence epilepsy [Genetic Absence Epileptic Rats of FStrasbourg and Wag/Rij rats (GAERS) =WistarAlbino Glaxofrom Rijswijk] the detailed observation of the evolution of the typical SWDs in the latter, using appropriate techniques, allowed Meeren et al. (2002) to identify the focal cortical onset of SWDs as the peri-oral region of the somatosensory cortex (**Figure 5**). A crucial point is that these local SWDs propagate very quickly throughout the cortex and to the thalamus at the millisecond scale. This fast propagation can only be accurately determined using appropriate analytical methods, such as non-linear regression analysis, at the very beginning of the burst of SWDs; within a few hundreds of milliseconds the propagation to other cortical areas and to the thalamus feeds back to the initial cortical area, what confounds any possibility of determining later where the SWDs had started.

These observations have two important methodological consequences: (1) any method of analysis needs to be reliably applicable to very short signal epochs in the order of <500 ms: (2) the sampling both in time and space has to be very high, at the millimeter and millisecond scales.

This cortical focal onset of SWDs first identified in the cortex of Wag/Rij rat was subsequently further characterized by means of intracortical and intracellular recordings by Polack et al. (2007, 2009) and Pinault (2003) in GAERS, what is in line with the previous observations that the appearance of SWDs in the electrocorticogram precede the corresponding discharges in the thalamus by Seidenbecher et al. (1998). The hypothesis that absence seizures have a focal cortical origin in WAG/Rij rats is also supported by the abolition of SWDs after pharmacological inhibition of the focal region (Sitnikova and van Luijtelaar, 2004; van Luijtelaar and Sitnikova, 2006) and by a series of studies showing abnormal excitability of the cortical focal region (see also Lüttjohann et al., 2011). Furthermore an integrated fMRI/EEG study in GAERS confirmed also the leading role of the somatosensory cortical focal area in SWDs, where the hemodynamic signals were analyzed using Granger causality and DCM along with local field potentials (David et al., 2008).

#### **Can the cortical focal origin of SWDs in rat absence models be extrapolated to human patients?**

Early classic observations of Bancaud et al. (1974), Rodin et al. (1994), Niedermeyer et al. (1969), and Niedermeyer (1972) suggested that in CAE there was a focal onset of SWDs in the frontal cortex. In the same line the temporal analysis of ictal absence EEG signals revealed a rapid motor involvement in cranio-caudal direction from the ocular/peri-oral regions to the extremities indicating dynamic propagation in a network involving the frontal lobe and motoric system (Stefan, 1982; Stefan and Carter Snead, 1997).

Notwithstanding these and other observations also pointing out in the direction that SWDs in absences are not primarily

generalized, this feature was overlooked in the Epilepsy community probably due to the overwhelming weight that the term " generalized" carries until now in the classification of epileptic seizures. Technical advances of EEG technology (Holmes, 2008), however, gave a new thrust to the search for possible sources of SWDs in the human cortex in patients with CAE. An important advance in this respect was achieved by Holmes et al. (2004, 2010) who recorded scalp EEG signals with a dense-array, 256-channel system in patients with "primary generalized epilepsy," absence seizures, and carried out source analysis estimating equivalent dipole distributions smoothed by linear inverse methods (LORETA, Pascual-Marqui et al., 2002). The onset of the slow component of SWDs was located at the frontal cortex and the spike at the frontopolar region of orbital frontal lobe. Furthermore a similar analysis (Holmes et al., 2010) in patients with Juvenile Myoclonic Epilepsy (JME) revealed sources in orbitofrontal/medial frontopolar cortex in all patients examined, and in half of the patients sources in basalmedial temporal lobe sources were found. Thus these authors conclude from these observations that JME is "not generalized in the sense of bilaterally diffuse onset."

Further clinical studies in patients with CAE, JME, and epilepsy with generalized tonic-clonic seizures (GTC) by Stefan et al. (2009) using MEG/EEG, demonstrated regional activations of the fronto polar medial cortex and rapid involvement of other brain areas (**Figures 6A,B**). The distribution of SWDs in absence patients appeared to involve a prefrontal-insular-thalamic network, whereas in patients with myoclonic components the dominant networks were (pre)motorinsular-thalamic.

Using 204-channel MEG recordings in patients with juvenile absence epilepsy (JAE) and applying dynamic statistical parameter mapping (sSPM) Sakurai et al. (2010) found that the onset of SWDs corresponded to focal cortical activation with secondary activation of posterior cingulate and precuneus, brain structures that belong to the "default mode network (DMN)," at least in some patients.

Since in the great majority of patients with apparently generalized Spike-and-Wave (GSW) and "Idiopathic Generalized Epilepsies" no iEEG recordings are made, several investigators have attempted to determine the role of the thalamus and the cortex in the generation of these discharges using imaging methods, namely fMRI associated with EEG. Using this methodology Salek-Haddadi et al. (2003) studied a patient with IGE and frequent absences and found that generalized SWDs were time-locked with bilateral BOLD activation of the thalamus and cortical deactivation most prominent in the frontal cortex. A similar study by Aghakhani et al. (2004), in a group of patients with IGE and GSWDs described bilateral activation in thalamic regions associated with SWDs, while in the cortex deactivations were observed.

Moeller et al. (2009)investigated not only GSWDs in patients with absence seizures but also triggered photoparoxysmal responses (PPRs) using EEG based coherent source dynamic imaging and BOLD signals. Partial directed coherence analysis indicated that the thalamus appeared to act as "pacemaker" of GSWDs in absence seizures; in contrast PPRs could be accounted for by an activation of the occipital cortex that propagates along cortical networks to frontal areas.

Carney et al. (2012) used also EEG-fMRI to further investigate the role of the frontal cortex in absence seizures. They identified two major patterns of frontal cortical BOLD signal change following onset of SWDs using event-related time course analyses, in the dorsolateral prefrontal cortex (DLPF): one group showed a pronounced and prolonged positive cortical BOLD signal change, whereas another group displayed a less pronounced BOLD increase followed by a predominant negative BOLD signal change. Similar patterns were found in the midline and lateral parietal cortex, caudate, and thalamus. They report also evidence of BOLD signal changes that precede the Spike-and-Wave onset, particularly in the mesial frontal cortex, parietal cortex, and precuneus, but not in the thalamus.

These authors suggest that the differences in frontal cortical BOLD associations with onset of absences may have phenotypic implications. This implies that group-averaged data have to be interpreted with caution, and that individual recordings have always to be examined.

A general comment about the findings of these studies in the light of the experimental evidence obtained in animal models of absences, is that the dynamics of the onset and early propagation of SWDs from the cortical onset zone to the thalamus takes place at a very fast pace. Only during a couple of hundreds of milliseconds there is a sustained flow of SWD signals from cortex to thalamus; after that the cortico-thalamo-cortical loop is entrained in the oscillations and it is not possible anymore to identify where the

onset was situated. Therefore BOLD signals, that have much longer time constants don't have the adequate time resolution to catch the dynamical onset event of SWDs. The changes in BOLD signals put in evidence by all the studies referred to above, however, are consistent: during a burst of SWDs there is activation of thalamus and deactivation of cortical areas, what represents the steady-state condition of the underlying neuronal networks during a SWD burst, but are not informative regarding the identification of the onset dynamics of SWDs.

activity). Numbers correspond to time points in the waveform. **(A)** Activity propagated from left frontal mesial area, to larger frontal areas, including

#### **Does network analysis reveal abnormalities of the inter-ictal state in patients with IGE?**

Children with CAE can display besides SWDs associated with absences also SWDs that may be clinically silent. Li et al. (2009) using EEG-fMRI analyzed both types of GSWDs, inter-ictal and ictal and reported that both types were associated with changes of BOLD signal in the basal ganglia-thalamo-cortical loop, but whereas the ictal type showed widespread and symmetrical deactivation in the cortex, the inter-ictal type showed predominant cortical activation. The authors advance the hypothesis that the cortical deactivation would be the substrate for the abrupt loss of consciousness of the absences. Interestingly Luo et al. (2011) investigated, using fMRI, the inter-ictal state also in patients with absences but taking care of avoiding that SWDs were present during the recording, i.e., in the resting state. Cross-correlation functional connectivity analysis revealed decreased integration within the DMN in the absence epilepsy patients, in particular a decrease of functional connectivity among the frontal, parietal, and temporal lobe. It is to be investigated whether this DMN abnormalities may be related to cognitive impairments in these patients.

to color in this figure legend, the reader is referred to the web version of

the article; Stefan et al., 2009).

#### **NEW APPROACHES TO DELINEATE EPILEPTOGENIC NETWORKS ALSO WITH REGARD TO GUIDING THERAPY**

In the previous sections we considered how various modern methodologies can be applied to improve the delineation of the epileptogenic zone having in mind that the latter should be envisaged as an *epileptogenic neuronal network*, and we discussed several typical clinical cases. In the modern conceptualization of epilepsy one has to evolve from "zones to networks" paraphrasing Laufs (2012). The identification of epileptogenic networks is, of course, of paramount importance in order to improve planning of resective surgery, and also to guide targeted therapies with the aim of controlling abnormal activity in relevant hubs and nodes within the epileptogenic network.

In the context of surgical planning in patients with medically intractable epilepsies one crucial aspect is to optimize the placement of depth electrodes (iEEG), such that the brain space of interest may be appropriately explored. Currently the possibilities offered in this respect by combining EEG and fMRI are the object of active research (Vulliemoz et al., 2010). A topical issue is whether fMRI may have a relevant added-value to classic EEG scalp recordings regarding the analysis of epileptogenic networks. Of course fMRI being non-invasive would represent an important practical tool in this context. Also in association with iEEG, fMRI might be valuable since it permits to extend the scope of the search for biomarkers of epileptic activity to the whole brain, compensating the spatial sampling limitations of iEEG. It is important to emphasize that EEG-fMRI studies should take into account not only the question of localization within brain space of "hot spots"of epileptiform activities but also the dynamic features of these activities, i.e., the flows of propagation and the corresponding time delays. In many investigations these dynamic aspects are still too little explored.

Some promising findings, however, have been reported. In a study of patients with focal epilepsies Jacobs et al. (2009) demonstrated the occurrence of BOLD changes associated with inter-ictal spikes recorded at the scalp that preceded the latter by a few seconds. This early BOLD response may be interpreted as resulting from changes in neuronal activity in epileptogenic neuronal networks situated deep in the brain that are not reflected at the level of the scalp EEG.

The study of Fahoum et al. (2012) addressed the same question by investigating the distribution of cortical and subcortical hemodynamic changes associated with inter-ictal spikes (IEDs) recorded in scalp EEG in patients with different epileptic conditions (temporal lobe-TLE, frontal lobe-FLE, and posterior quadrant-PQE epilepsies) using a similar EEG-fMRI approach. These authors modeled the BOLD response to the IEDs using the timing of the epileptiform events as regressor convolved with a series of Hemodynamic Response Functions (HRFs) consisting of gamma functions peaking at successive delays. Without going here into details the main findings of the analyses showed widespread clusters of activation and deactivation in TLE and FLE patients, while in PQE only deactivations clusters were found, that reached brain areas outside the presumed epileptogenic zone. The largest activations both in TLE and FLE patients were found bilaterally in mid-cingulate gyri. All patient groups showed deactivations of DMN regions, particularly in TLE patients (inferior parietal lobules, posterior cingulate cortex, and precuneus bilaterally). The involvement of the mid-cingulate gyri likely reflects the rapid propagation of epileptic activity from sources in temporal and frontal areas. The pathophysiological significance of this finding is, however, not yet clear. Nonetheless this may be a pointer to further investigate whether these networks involving the cingulated cortex may be interesting targets for therapeutic interventions.

Simultaneous intra-cranial recordings with an appropriate spatial sampling, however, are necessary to clarify these functional relationships revealed by the scalp IEDs – fMRI analyses. In a study of patients with FCDs Thornton et al. (2011) made a comparison between EEG-fMRI signals associated with inter-ictal epileptiform discharges (IEDs) in iEEG recordings in order to delineate in a classic way the seizure onset zone (SOZ). These authors studied also the surgical outcome of these patients. About 5 of 11 patients showed IED-related BOLD signals that were concordant with the electrophysiologically determined SOZs, that in these patients had a limited extent. Six of 11, however, did not display this concordance and the IED-related BOLD signals revealed more widespread epileptogenic regions in comparison with the extent of the SOZ delineated based only on iEEG data. Most interesting the five former patients had a good surgical outcome, but this was not the case for the latter 6. These findings suggest that EEG-fMRI may be useful to identify patients with extensive epileptogenic networks that extend beyond those delineated using only classic iEEG. This information may contribute to making decisions concerning surgical resections more appropriately.

Thus EEG/fMRI studies may be helpful in order to plan more efficiently iEEG recordings and to estimate more accurately the extent of epileptogenic neuronal networks. One critical note should be added: many EEG-fMRI studies tend to focus mainly on localizations rather than on dynamics. In general these studies don't take into account time-delays between different hubs and nodes within extensive epileptogenic networks. In order to obtain this information one has to resort to electrophysiological measurements, given the low time resolution of fMRI signals. To obtain a comprehensive picture of epileptogenic networks it is essential to uncover the dynamics of the propagation of epileptic activity in these networks.

It is interesting to note that these techniques are being applied also to find whether the pattern of functional connectivity of neuronal networks may yield relevant information about epileptogenic networks during inter-ictal resting states. Using non-linear correlation analysis Bettus et al. (2011) showed in TLE patients a general increase of iEEG signal interdependencies (specific for the beta frequency band) in regions affected by electrical epileptiform abnormalities relative to non-affected areas, whereas the opposite pattern was found for functional connectivity measured using fMRI signals. This latter finding may be due to anomalies in metabolism and in neurovascular coupling (blood-brain-barrier permeability) in epileptogenic networks that may be affected in TLE. Such anomalies have been also shown in animal models of TLE and in human patients by means of MR diffusion imaging – tractography (Yogarajah et al., 2008). Bettus et al. (2011) suggest that this increase in beta frequency interdependencies in epileptogenic networks"could be a reliable pathological marker of epileptic processes." This needs further confirmation.

With respect to locally targeted therapies experimental work in animal models reveals some interesting novel perspectives for future therapeutic interventions, for example with the purpose of delivering anticonvulsants locally into specific hubs of epileptogenic networks as proposed by Löscher and collaborators (Bröer et al., 2012). Strategies for neuromodulation aiming at the therapeutic control of epileptogenic networks are being considered besides local drug delivery, such as local electrical stimulation, transcranial magnetic stimulation, stem cells transplantation, and gene therapy (see review Al-Otaibi et al., 2011). Regarding deep electrical brain stimulation (DBS) Kahane and Depaulis (2010) stress the importance of gaining a better understanding of the functional properties of epileptogenic neuronal networks in which seizures originate and propagate, as much as of the mechanisms by which neurostimulation works, in order to define the types of DBS that may be effective. The need of acquiring a comprehensive insight in these functional epileptogenic networks applies to all strategies to develop novel targeted therapies.

#### **CONCLUSION**

The interest on neuronal network analysis in epilepsy has gained strength with the use of high resolution recording techniques and signal analytical methodologies that opened up the possibility of studying dynamic brain states with high resolution both in space and time. In addition to invasive EEG recordings, non-invasive recording techniques like MEG/EEG and EEG-fMRI are being used more and more to identify network involvements in various epileptic conditions and age dependent syndromes, both structurally and dynamically. In some cases diverse techniques are used in a complementary way. A typical example is the study ofVaudano et al. (2012) who made ictal MEG and EEG-fMRI recordings in a patient with reading epilepsy. Using this information combined with DCM, they found evidence for a causal link between activity in the left piriform cortex and the seizures elicited by reading. More research is needed to integrate this kind of findings with the

#### **REFERENCES**


C., Juler, J., Polkey, C. E., et al. (1997). Origin and propagation of interictal discharges in the acute electrocorticogram. Implications for pathophysiology and surgical treatment of temporal lobe epilepsy. *Brain* 120(Pt 12), 2259–2282.


detection of relevant biomarkers of epileptogenic zones, that have been recently described such as (fast)ripples (HFOs, Jacobs et al., 2012; Jefferys et al., 2012) in order to refine systematic network analyses in all sorts of epileptic patients.

In the last decade the development of new methodologies to analyze the dynamics of neuronal networks has gained momentum, and has yield a wide range of computer tools that are being tested in clinical and experimental environments. Using such tools, many of which were discussed above, new insights are emerging.

One of these is that the concept that some types of epilepsy are "generalized" is outdated. Even in cases of IGE and CAE there is now compelling evidence that the seizures start in a well defined brain area and spread at great speed to connected brain areas recruiting specific neuronal networks into typical oscillatory behavior. This conclusion is now supported by high resolution human EEG and MEG data, by animal experimental electrophysiological and fMRI data. Thus the dividing line between "generalized" and "focal" epilepsies recedes on the face of new relevant neuronal network analyses. In the classification of seizures, it is important, in addition to the localization of seizure onset, to identify cortical and subcortical network involvement and the speed of propagation. This implies that a neurophysiological seizure classification approach should include the identification of the predominant involved networks,within one hemisphere or involving both. This conclusion is in line with Laufs' (2012) conclusion that the concept that epilepsies should be considered as resulting from disturbed network interactions that implies "multitargeted treatments." A major challenge will be to perform connectivity analysis in order to identify the primary epileptogenic source, or "hub," in complex multifocal epilepsies in order to optimize tailoring surgical resections and/or targeted therapies in epileptic patients.

#### **ACKNOWLEDGMENTS**

Supported by DFG STE 15-1.


epilepsies: entering the 21st century. *Epilepsia* 52, 1058–1062.


*J. Neurosci. Methods* 30, 125, 195–207.


N. M., et al. (2003). Consistency of interictal and ictal onset localization using magnetoencephalography in patients with partial epilepsy. *J. Neurosurg.* 98, 837–845.


electroencephalographic signals in the time and frequency domain: identification of epileptogenic networks in partial epilepsy. *Philos. Trans. A Math. Phys. Eng. Sci.* 367, 297–316.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 17 October 2012; paper pending published: 10 December 2012; accepted: 24 January 2013; published online: 01 March 2013.*

*Citation: Stefan H and Lopes da Silva FH (2013) Epileptic neuronal networks: methods of identification and clinical relevance. Front. Neurol. 4:8. doi:10.3389/fneur.2013.00008*

*This article was submitted to Frontiers in Epilepsy, a specialty of Frontiers in Neurology.*

*Copyright © 2013 Stefan and Lopes da Silva. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, providedthe original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*

#### **Marco Leite1,2, Alberto Leal 3,4 and Patrícia Figueiredo1,2\***

<sup>1</sup> Department of Bioengineering, Instituto Superior Técnico, Technical University of Lisbon, Lisbon, Portugal

<sup>3</sup> Centro de Investigação e Intervenção Social, Lisbon, Portugal

<sup>4</sup> Department of Neurophysiology, Centro Hospitalar Psiquiátrico de Lisboa, Lisbon, Portugal

#### **Edited by:**

Mark Holmes, University of Washington, USA

#### **Reviewed by:**

Jose F. Tellez-Zenteno, University of Saskatchewan, Canada Don Tucker, Electrical Geodesics, Inc. and the University of Oregon, USA

#### **\*Correspondence:**

Patrícia Figueiredo, Department of Bioengineering, Institute for Systems and Robotics, Instituto Superior Técnico, Av. Rovisco Pais, 1, 1049-001 Lisboa, Portugal. e-mail: patricia.figueiredo@ist.utl.pt

Simultaneous electroencephalogram (EEG)-functional Magnetic Resonance Imaging (fMRI) recordings have seen growing application in the evaluation of epilepsy, namely in the characterization of brain networks related to epileptic activity. In EEG-correlated fMRI studies, epileptic events are usually described as boxcar signals based on the timing information retrieved from the EEG, and subsequently convolved with a hemodynamic response function to model the associated Blood Oxygen Level Dependent (BOLD) changes. Although more flexible approaches may allow a higher degree of complexity for the hemodynamics, the issue of how to model these dynamics based on the EEG remains an open question. In this work, a new methodology for the integration of simultaneous EEG-fMRI data in epilepsy is proposed, which incorporates a transfer function from the EEG to the BOLD signal. Independent component analysis of the EEG is performed, and a number of metrics expressing different models of the EEG-BOLD transfer function are extracted from the resulting time courses.These metrics are then used to predict the fMRI data and to identify brain areas associated with the EEG epileptic activity.The methodology was tested on both ictal and interictal EEG-fMRI recordings from one patient with a hypothalamic hamartoma. When compared to the conventional analysis approach, plausible, consistent, and more significant activations were obtained. Importantly, frequency-weighted EEG metrics yielded superior results than those weighted solely on the EEG power, which comes in agreement with previous literature. Reproducibility, specificity, and sensitivity should be addressed in an extended group of patients in order to further validate the proposed methodology and generalize the presented proof of concept.

**Keywords: BOLD, EEG-fMRI, epilepsy, ICA, heuristic**

#### **INTRODUCTION**

Over the years, the electroencephalogram (EEG) has been the tool of choice for the diagnosis and characterization of epilepsy. With the possibility to acquire the EEG simultaneously with functional Magnetic Resonance Imaging (fMRI), studies of Blood Oxygen Level Dependent (BOLD) signals correlated with epileptic activity proliferated (Ives et al., 1993; Hoffmann et al., 2000; Lemieux et al., 2001; Salek-Haddadi et al., 2006; LeVan and Gotman, 2009). Despite its potential for the localization of epileptogenic brain networks in patients with drug-resistant focal epilepsy undergoing evaluation for surgical treatment, simultaneous EEG-fMRI has yet to reach its full potential in clinical practice. One factor contributing to this state of affairs is the lack of sensitivity in the identification of hemodynamic changes associated with the EEG epileptiform discharges in a significant number of studies (Aghakhani et al., 2006; Salek-Haddadi et al., 2006; Gotman, 2008; Grouiller et al., 2011). Although technical difficulties related with data acquisition and artifact correction may in part explain such results, probably most important are the limitations related with the remaining conceptual and methodological challenges associated with the integration of the two types of signals.

Although a great amount of both experimental and theoretical work has been dedicated to the clarification of the relationship between neural activity and associated hemodynamic changes, neurovascular coupling mechanisms remain an active area of research (Rosa et al., 2010a). The most consensual evidence comes from the recording of electrical activity using micro-electrodes implanted in the cortex of (non-human) animals, simultaneously with fMRI, indicating that the BOLD signal reflects mostly slow, post-synaptic input activity measured by local field potentials (LFPs), rather than fast, spiking output activity measured by single/multi-unit activity (S/MUA; Logothetis et al., 2001). In humans, a growing number of simultaneous EEG-fMRI studies on healthy subjects as well as epilepsy patients have now been reported (Goldman et al., 2002; Laufs et al., 2003, 2006; Moosmann et al., 2003; de Munck et al., 2009), and biophysical models of the neurovascular coupling have been proposed (Riera et al., 2006, 2007). Overall, reports in the literature do not provide a clear picture of the link between EEG and BOLD signals. In particular, contradictory results have been presented regarding the dependency of BOLD changes on the EEG power and spectral profiles. These include, for example, positive and negative BOLD

<sup>2</sup> Institute for Systems and Robotics, Lisbon, Portugal

correlations with specific frequency band power changes in the human EEG (de Munck et al., 2009), BOLD decoupling from LFP power in mice (Ekstrom, 2010), and negative BOLD associated with large increases in LFP and MUA during seizures also in mice (Schridde et al., 2008). Rosa et al. (2010b) addressed this topic by comparing different models of the transfer function between EEG and BOLD signals, in the prediction of fMRI data, in a visual stimulation experiment with human healthy subjects. The models explored included the EEG total power (TP; Wan et al., 2006); linear combinations of the power from different frequency bands (Goense and Logothetis, 2008); and several variations of a heuristic model proposed by Kilner et al. (2005)in which BOLD changes are assumed to be proportional to the root mean square frequency (RMSF) of the EEG spectrum. The results obtained showed a clear superiority of the RMSF metrics in predicting the BOLD signal, when compared to power-weighted metrics.

In most epilepsy EEG-fMRI studies, the goal is to find brain networks exhibiting hemodynamic changes associated with interictal and/or ictal activity, which are expected to be correlated with the epileptogenic areas (Hoffmann et al., 2000; Salek-Haddadi et al., 2006; Marques et al., 2009). Both ictal and interictal events are identified on the EEG trace and are then used to define regressors of interest in a general linear model (GLM) analysis of the fMRI data. Interictal spikes are usually described as stick functions and ictal activity as boxcar signals between seizure onset and offset, eventually sub-divided into up to three phases: early ictal EEG, clinical onset, and late ictal EEG. Both types of events are then convolved with a model of the hemodynamic response function (HRF; Tyvaert et al., 2008; Moeller et al., 2010; Thornton et al., 2010). Independent component analysis (ICA) of the fMRI data has also been performed in order to identify interictal/ictal BOLD patterns without resorting to the EEG (LeVan et al., 2010; Thornton et al., 2010). When the EEG accurately reflected the seizure onset, the GLM approach yielded activations concordant with the ictal onset zone; otherwise ICA still gave valuable insight on the ictal hemodynamics. Importantly, the question of whether or not the neurovascular coupling is preserved from healthy to disease conditions has also been investigated, by allowing variations in the HRF (Grouiller et al., 2010). It is in principle possible to achieve a higher degree of complexity for the epileptic activity hemodynamics by using more flexible HRFs, or completely model-free approaches such as ICA. However, such approaches do not address the issue of how to model these dynamics based on the information available in the EEG data (Lemieux, 2008).

In this paper, we address the issue of modeling the BOLD dynamics associated with epileptic EEG activity by extracting metrics from the EEG spectrum expressing different models of the transfer function between neuronal and hemodynamic signals. A methodology is proposed for the analysis of EEG-fMRI data in epilepsy, consisting of ICA decomposition of the EEG followed by component selection based the reproducibility across different acquisition runs, Morlet wavelet spectral analysis, and EEG metric extraction. The resulting time courses are convolved with a canonical HRF and used as regressors of interest in a GLM analysis of the fMRI data. The proposed methodology is applied to the study of a patient with epilepsy associated with a hypothalamic hamartoma. The different EEG-fMRI transfer functions are compared with each other, as well as with a conventional GLM methodology based on the identification of ictal and interictal events on the EEG by the neurophysiologist, and also with an fMRI-ICA approach.

#### **MATERIALS AND METHODS PATIENT CHARACTERIZATION**

We focused on the simultaneous EEG-fMRI data recorded from a 2-year-old patient with a giant hypothalamic hamartoma suffering from gelastic epilepsy, as part of the pre-surgical evaluation under the Program for Epilepsy Surgery of the Hospital Center of West Lisbon. This case study has been previously described (Leal et al., 2009). This patient was studied in a run of 30 patients, from which only five had ictal events during the scanning sessions. From these five, only this one patient had an EEG trace clear from movement related MR artifacts triggered by the beginning of the seizures and therefore presented sufficient data quality for subsequent EEG quantification and was selected for this study.

Seizures occurred more than 50 times per day and typically lasted for 20–30 s, involving almost exclusively the left hemisphere. The clinical manifestations of the seizures consisted of slowing of motor activity, variable interruption of consciousness, eyelid rhythmic movements with bilateral nystagmus to the right, and occasionally gelastic laughter. After the acquisition of the EEG-fMRI data, the patient underwent a two-stage hamartoma disconnection surgery, 1 year after which the seizures were reduced to 1–3 per day.

The EEG interictal activity demonstrated a persistent slowwave abnormality over the left temporal-occipital area, associated with abundant spike activity with occasional contralateral propagation. Abundant spike activity also occurred over the left hemisphere frontal lobe. Topological analysis of the interictal spikes presented a spatial stationarity for the posterior spikes,whereas the frontal ones changed significantly in configuration from an occipital dipolar potential at spike onset to a dipolar frontal potential at spike peak. The ictal EEG pattern was very monotonous and consisted of early diffuse desynchronization, followed by the build-up of spike activity over the left occipital and temporal areas and, in the later stages of the seizure, over the frontal area. Occasionally secondary propagation of spike activity to the right temporal areas occurred.

#### **SIMULTANEOUS EEG-fMRI ACQUISITION**

The EEG was recorded using an MR-compatible 37-channel system (Maglink, Neuroscan, Charlotte, NC, USA) with two ECG channels, using a sampling rate of 1000 Hz with a bandwidth DC-250 Hz and reference at electrode FCz. The imaging was performed on a 1.5 T MRI scanner (GE Cvi/NVi). Six fMRI runs were collected using a gradient-echo echoplanar imaging (EPI) sequence, with TR = 2.275 s, 3.75 mm × 3.75 mm × 5.00 mm voxel size and a total of 154 volumes (the first four were subsequently rejected in order to discard T<sup>1</sup> relaxation unstationarities). A T1-weighted structural image was also acquired using a 3D spoiled gradient recovery (SPGR) sequence, with 0.94 mm × 0.94 mm in-plane resolution and 0.6 mm slice thickness. During the scanning session, the patient was administered light anesthesia with Sevoflurane at 1% (Abbott Laboratories, Abbot Park, IL, USA), through mask, as established by the MRI protocol for small children and

uncooperative patients. A 5 min EEG recording was also performed outside the scanner and before anesthesia,in order to allow for cross-validation of the simultaneous EEG-fMRI recordings in terms of the presence of MR related artifacts and the effects of anesthesia. Periods of interictal and ictal activity were identified on the artifact-corrected EEG by the neurophysiologist. For runs 2 and 3, two and five ictal events occurred, respectively. For the remainder runs, only interictal spikes were detected on the EEG. The visual inspection of the EEG inside and outside the scanner revealed a slight increase in both beta and theta background activity under anesthesia, but there was no significant change in the morphology of interictal spikes.

### **EEG ANALYSIS**

#### **Pre-processing**

The EEG pre-processing was executed using the EEGLAB toolbox (Delorme and Makeig, 2004). Firstly, the EEG traces were visually inspected for the presence and consequent rejection of bad channels, associated with poor contacts. A 2 Hz high-pass filter was then applied so as to remove baseline drifts from the signal. The times of occurrence of the gradient artifacts associated with the acquisition of each fMRI slice were automatically identified on the EEG signals, by software developed in-house. The fMRI gradient artifact correction algorithm, FMRIB's FASTR (Niazy et al., 2005), was then applied using the default parameters. For the pulse artifact removal, FMRIB's QRS complex identification algorithm (Niazy et al., 2005) was first applied to the ECG channels. An optimal basis set of three principal components was then employed for pulse artifact removal of the data, after low-pass filtering at 45 Hz and down-sampling to 100 Hz to improve manageability.

#### **ICA decomposition**

In an attempt to separate out the activity of interest in the EEG, the pre-processed data were decomposed by ICA using the *infomax* algorithm as implemented in EEGLAB (Delorme and Makeig, 2004). The reference channel was arbitrarily kept as the one chosen by the electrophysiologist during the acquisition. Although the referencing method for the EEG channels does not affect the final IC time courses, because the reference channel is linearly separable from the data, it will affect the IC's scalp topographies. Nevertheless, the reference channel was kept the same (FCz) throughout the acquisition of all six runs, so comparisons between the scalp topographies of components obtained in different runs are still possible.

A reproducibility analysis of the ICs was performed in order to identify the associated topographies that were consistent across the six acquisition runs. *IC groups* were built by fixing each IC of each run and selecting, for each of the other runs, the IC that was the most spatially correlated with it. A total of 6 runs × 25 ICs = 150 IC groups were hence generated, each composed of a string of six ICs. An IC group was considered to be *consistent* if the same associated string was generated by one IC of each run, and therefore was repeated six times in the total set of 150 IC groups. The topographies and time courses of the consistent IC groups were then visually inspected for the identification and consequent rejection of artifact related ICs, dominated by residual gradient artifacts, bad channels, or eye blink/movement artifacts.

#### **EEG metrics**

The spectral profiles of the selected ICs were obtained by timefrequency analysis through convolution of the respective time courses with Morlet wavelets. The subset of metrics found to be more relevant in Rosa et al. (2011) were applied here: mean frequency (MF), RMSF, un-normalized mean frequency (uMF) un-normalized mean square frequency (uRMSF), and TP.

The difference between frequency averaging measures (RMSF or MF) was found to have a negligible effect on BOLD signal prediction; hence the results emerging from the MF and uMF metrics will be omitted, as they were not significantly different from those of the RMSF and uRMSF metrics, respectively.

#### **fMRI ANALYSIS**

The fMRI data were analyzed using FSL<sup>1</sup> , including: (1) preprocessing; (2) GLM; and (3) ICA.

Pre-processing consisted on: motion correction using MCFLIRT (Jenkinson et al., 2002); slice timing correction using (Hanning-windowed) sinc interpolation to shift each time-series by an appropriate fraction of a TR relative to the middle of the TR period; brain extraction using BET (Smith, 2002); temporal highpass filtering rejecting periods above 100 s; and spatial Gaussian filtering with FWHM = 8 mm.

Two GLM analyses were specified in order to identify BOLD signal changes associated with: (1) the *EEG metrics* for each consistent IC group (as defined in the previous section); and (2) the epileptic activity identified by the neurophysiologist on the pre-processed EEG traces, where boxcar signals were used to describe the periods of ictal activity and stick functions were used to describe interictal spikes (from now on referred to as *EA*). Each of these variables of interest (*EEG metrics* and *EA*) was convolved with a canonical HRF (Friston et al., 1998), and the time and dispersion derivatives were also included to account for some degree of variability in the HRF shape across the patient's brain. This resulted in a set of three regressors for each variable. The final GLM regressors were obtained by re-sampling the resulting time courses to match the middle of the acquisition time period of each fMRI volume. Six motion parameters were also included as regressors of no interest, in order to account for residual motion-related signal jitter not removed by the motion correction procedure. The GLM's were fitted to the data using the FILM algorithm (Woolrich et al., 2001) *F* tests were then applied to each estimated parameter contrast, resulting in *Z* (Gaussianized *F*) statistic maps. These were thresholded using a clustering procedure, whereby each cluster is determined by a voxel *Z* > 2.3 and a (corrected) cluster significance threshold *p* = 0.05.

An ICA of each fMRI run was also performed using Probabilistic ICA as implemented in MELODIC (Multivariate Exploratory Linear Decomposition into Independent Components) Version 3.09, part of FSL (FMRIB's Software Library (see text footnote 1); Beckmann and Smith, 2004). Pre-processing included voxelwise de-meaning of the data and normalization of the voxel-wise variance. Estimated component maps were divided by the standard deviation of the residual noise and thresholded by fitting a mixture model to the histogram of intensity values.

<sup>1</sup>www.fmrib.ox.ac.uk/fsl

#### **MODEL COMPARISONS**

The comparison between EEG metrics was performed through GLM analysis followed by inferences based on *F* tests, for each consistent IC group of interest and for each acquisition run. We first considered a single GLM comprising the three EEG metrics (TP, RMSF, and uRMSF) with three regressors of interest each (canonical HRF, time, and dispersion derivatives), totaling nine regressors of interest, in one three-way comparison. Three GLM's contrasting pairs of EEG metrics (TP vs. RMSF, TP vs. uRMSF, and RMSF vs. uRMSF) were also considered, totaling six regressors of interest each, in three two-way comparisons. *F* tests were computed for each set of three regressors regarding each metric, as well as for the whole set of metrics. The resulting *Z* statistical maps were inspected for their volume, i.e., the total number of voxels exhibiting significant EEG metricrelated BOLD changes. This was used as a quantitative measure of the performance of each EEG metric in terms of BOLD prediction, and the analysis was repeated for each consistent IC group.

The comparison between consistent IC groups of interest was performed in a similar way, for each acquisition run, applying the selected EEG metric. Here, we considered a single GLM comprising the three consistent IC groups of interest (I, II, and V, as will be shown in the Results) with three regressors of interest each (canonical HRF, time, and dispersion derivatives), totaling nine regressors of interest, in one three-way comparison.

In order to assess the plausibility and consistency of the EEG metric-derived BOLD maps, the EEG metric/IC group combination yielding the largest maps (and hence best at predicting BOLD signal changes) was selected for subsequent comparison with the GLM analysis using the EEG EA regressors and also with the ICA analyses. The consistency between each two regressors was measured as their temporal correlation. The consistency between two maps was measured as their spatial overlap, i.e., the ratio of the volume of the map intersection with the volume of the map union.

### **RESULTS**

In this section, the results obtained through the proposed methodology will be presented.

#### **EEG ANALYSIS**

An exerpt of the EEG trace obtained after pre-processing is presented in **Figure 1**. The data appears to be clear from MR related artifacts and the ictal activity can be clearly identified.

The consistent IC groups, obtained as a function of the reproducibility of the associated IC topographies across runs, are presented in **Figure 2**. The IC groups I, II, and V were considered artifact free and were kept for further analysis, while the remainder were rejected. The IC groups I and II exhibit clear dipolar configurations in the left hemisphere, compatible with the frontal and occipital/parietal patterns of the patient's interictal and ictal EEG. Given its predominantly frontal topography, particular attention

was given to IC group II in order to confirm that ocular movement artifacts did not dominate the IC time course. Topography group V exhibits a more diffuse configuration difficult to interpret.

The spectral profiles obtained with Morlet wavelet decomposition for one IC (IC10) in run 3, as well as the EEG channel with highest absolute weight for this component (P3), are shown in **Figure 3**. Spectral changes associated with the ictal events identified by the neurophysiologist are visible in both spectrograms; however, these changes are clearer in the spectrogram obtained for the IC in comparison to the spectrogram of channel P3, or any other single EEG channel (data not shown). This observation suggests that ICA decomposition is capable of separating out the ictal activity spread over several EEG channels into a limited number of components.

#### **fMRI ANALYSIS**

The BOLD statistical maps obtained using each EEG metric and consistent IC group were generally consistent with each other, across runs and also with the patient's seizure semiology, but differed considerably in terms of the number of voxels showing statistically significant EEG metric-related BOLD changes. Firstly, the results of the comparison between the fMRI analysis using the different EEG metrics and consistent IC groups will be presented. A single metric and IC group will then be selected for comparison with the EA GLM and fMRI-ICA analyses.

#### **EEG metric comparisons**

The three-way analysis did not show meaningful activations to allow for the comparison of the individual EEG metrics; however, the two-way comparisons yielded a clear superiority of the RMSF metric when compared to the TP and uRMSF metrics, for every IC group and run, in terms of the number of voxels showing statistically significant EEG metric-related BOLD changes, as shown in **Figure 4**.

The three-way comparison regarding the IC group showed a significant superiority of IC group II, in terms of the number of voxels exhibiting statistically significant EEG metric-related BOLD changes, as shown in **Figure 5**. For five out of the six runs, the statistical maps yielded by IC group II had larger volumes of significant voxels than those yielded by other ICs. The EEG RMSF metric of IC group II was found to be the best at predicting BOLD changes and it will now be compared with the results of the EA GLM analysis as well as with those of the fMRI-ICA approach.

The BOLD statistical maps obtained using the RMSF metric for IC groups I and V are shown in **Figures 6** and **7**.

#### **EEG metric vs. EA GLM analysis**

The comparison between the fMRI results obtained using the EEG RMSF metric and the EA GLM analysis, for IC group II, are summarized in **Figure 8**. For each run, the regressors associated with the canonical HRF are plotted alongside with the *Z* (Gaussianized *F*) statistic maps obtained with both methods. The temporal correlation between the regressors and the spatial overlap between the respective maps are also presented. The consistency between the two methodologies is generally low for all runs (correlation < 0.3 and overlap < 3), with only run 3 exhibiting an overlap between maps above 10%. Interestingly, this is also the run during which more ictal events occurred.

In general, the maps obtained using the EEG metric approach were consistent across runs and also with the patient's seizure semiology. In fact, clusters were found on the left occipital/parietal and left frontal lobes, consistently with the observation on the EEG of spike buildups in left occipital/parietal channels followed by

the left frontal channels, after a generalized desynchronization. In run 3, the EA approach was also capable of identifying these brain areas, but the corresponding statistical maps were more significant and extensive for the EEG metric approach. Furthermore, with the EEG metric approach, clusters involving left thalamic, left hippocampal, and left frontal ventral areas, as well as the hamartoma itself, were also found, which are generally consistent with the brain network associated with seizures originating in a hypothalamic hamartoma. For the remaining runs, in contrast with the EEG metric approach, the fMRI statistical maps obtained with the EA approach were in general not consistent across runs nor with the expected epileptic network.

#### **EEG metric vs. fMRI-ICA approach**

The results regarding the EEG metric GLM analysis and the fMRI-ICA decomposition of all the runs, for IC group II, are presented in **Figure 9**. The fMRI-ICs presented are those that yielded the highest spatial map overlap with the results of the corresponding EEG metric analysis. Map spatial overlaps and time course temporal correlations between EEG metric GLM and fMRI-ICA approaches were generally higher than those found between EEG metric GLM and EA GLM approaches (**Figure 8**). Moreover, in contrast with the EA approach, with fMRI-ICA, consistency across runs was also observed for the brain areas expected in the patient's epileptic network. Interestingly, runs where only interictal activity was recorded (1, 4, 5, and 6) yielded maps consistent with runs with ictal activity (2 and 3).

#### **DISCUSSION**

We proposed a new methodology for BOLD signal prediction in EEG-correlated fMRI studies in epilepsy, by incorporating a model of the EEG-BOLD transfer function. Specifically, independent components of the EEG associated with consistent topographies were translated into BOLD signal predictions by a set of modelbased metrics. Interestingly, we found that increases in the MF of

the EEG were better than its power at predicting BOLD increases, in support of the heuristic proposed in (Kilner et al., 2005) and in agreement with the results obtained in a visual stimulation experiment with healthy subjects (Rosa et al., 2010b). Moreover, such EEG-based metrics were found to improve detection sensitivity compared with conventional approaches to EEG-fMRI data analysis in epilepsy.

The heuristic proposed by Kilner and colleagues puts forward the MF (in an RMS sense) of the EEG spectrum as a surrogate of the neuronal activity eliciting the BOLD signal, following a broad physiological inspiration: the BOLD eliciting signal is assumed to be proportional to the electric work dissipated by the ionic currents across the cell membranes; this can in turn be shown to be proportional to the RMSF of the LFP/EEG spectrum if the covariance of the membrane potentials of the cells are assumed constant (Kilner et al., 2005). Although more detailed models of EEG-fMRI integration have been presented in the literature (Riera et al., 2006, 2007), the heuristic benefits from its simplicity in application and interpretability. The dependency of the BOLD signal on the EEG spectral profile has often been experimentally reported in studies of spontaneous or evoked fluctuations of brain rhythms, and the results are most frequently found to be in concordance with Kilner's heuristic predictions (Goldman et al., 2002; Gonçalves et al., 2006; Laufs et al., 2006; de Munck et al., 2009; Michels et al., 2010; Zumer et al., 2010; Khursheed et al., 2011). Rosa et al. (2010b) have explicitly employed the heuristic to analyze EEG-fMRI data collected during a visual flicker stimulation task in healthy subjects, and showed that BOLD signal decreases were indeed associated with changes in the EEG spectral profile, namely its RMSF, which did not arise from power changes in one specific frequency band. In epilepsy, low-frequency slow-wave activity increases have been shown to be associated with BOLD decreases

(Archer et al., 2003) while high-frequency spike and wave discharges have been shown to be associated with BOLD increases (Krakow et al., 2001; Hamandi et al., 2004). Although these finding are in agreement with the heuristic, our study is the first one to test it explicitly on EEG-fMRI data of epileptic activity.

The EEG metric-related BOLD change maps were consistent with the ones obtained using the epileptic activity regressors defined by the neurophysiologist, whenever ictal activity was identified on the EEG. For the runs on which no ictal events were detected on the EEG, the proposed methodology was able to identify the same brain network that was involved in the ictal runs. This was achieved by regressors based on the same IC scalp topography as those from the ictal runs, suggesting the presence of underlying epileptic activity in the identified epileptic network, which only occasionally manifests itself with an ictal character. In contrast with the proposed EEG metric based approach, however, the conventional analysis of the interictal runs yielded inconsistent or no results for the same statistical significance threshold. These findings suggest that the interictal events detected on the EEG may not fully reflect the activity of the underlying brain network, while the selected EEG metrics may be more powerful in depicting it. The same network was also identified by fMRI-ICA on both ictal and interictal runs, which further supports this idea. However, the fully data-driven method of fMRI-ICA lacks an explanatory model for the data, in contrast with the proposed methodology, which is based on modeling the EEG-BOLD transfer function.

For the purpose of verifying the specific epileptic character of the identified brain networks, their consistency across runs and their plausibility as the epileptic brain network underlying the patient's seizure semiology were considered. The overlap of the corresponding BOLD change maps across runs was very high, both for the EEG metrics GLM and the fMRI-ICA approaches.

regressors (bottom of graphs) and the map overlap (top right of maps) are shown.

Moreover, these maps were in good agreement with the brain network known to be involved in the epileptic activity of this patient, which includes the hamartoma, as well as left hemisphere hypothalamus, hippocampus, parietal–occipital lobe, cingulate gyrus, and dorsal–lateral frontal lobe (Leal et al., 2009). Clusters in the left parietal–occipital and frontal lobes are consistent with the observation on the patient's EEG of spike buildups in left parietal– occipital channels followed by the left frontal channels, after a generalized desynchronization. Clusters involving left thalamus and hippocampus, as well as the hamartoma itself, are generally consistent with the brain network associated with seizures originating in a hypothalamic hamartoma (Leal et al., 2003). Future work should be carried out in order to further validate the networks identified by our methodology. On a first approach, the

topographies of the selected ICs could be compared with those obtained by ICA decomposition of the EEG performed outside the MRI scanner. Ultimately, intra-cranial EEG recordings (unavailable for this patient) must be used in order to achieve a conclusive validation.

An ICA decomposition of the EEG was used here with the purpose of separating out the activity of interest into a univariate timecourse, which was expected to exhibit a consistent and meaningful topography. ICA is a popular technique for the removal of muscular activity or eye movement artifacts in EEG data processing (Vigario, 1997). Because of the large amplitude of interictal epileptic activity and the fact that its sources can generally be assumed to be spatially stationary, ICA has also proved to be useful in the separation and identification of such activity (Kobayashi

et al., 1999; Urrestarazu et al., 2006; Marques et al., 2009; Formaggio et al., 2011). However, ICA decomposition of EEG ictal events that involve spatial propagation may be questionable in the sense that the spatial stationarity assumption of the EEG sources is not verified. This is in fact the case in our study, where ICA was applied to the EEG recorded during seizures exhibiting a spatial propagation pattern associated (Leal et al., 2009). The results obtained may reflect this issue to some extent, as no single IC isolated *per se* all of the seizure dynamics. Nevertheless, ICA was useful for the separation of local approximately stationary activity, giving the ICs higher signal-to-noise ratio for the signals of interest, by separating them from activity in neighboring brain regions and also from residual artifacts not fully corrected by the EEG pre-processing as observable in the spectral analyses and scalp topographies.

Other approaches for the extraction of a univariate EEG time course, representative of the epileptic activity, have been proposed in the literature, namely continuous Electrical Source Imaging (Vulliemoz et al., 2010) and weighted averaging of selected ICs (Formaggio et al., 2011). The former corresponds to the projection of the recorded EEG data into the space of an EEG source estimation, which provides a more informed way of obtaining the EEG source signal than ICA. However, this technique is more prone to include artifacts in the resultant time courses when compared to ICA, as the latter automatically rejects channels contaminated with relevant artifacts. Regarding the averaging of selected EEG IC time courses, there is no straightforward way to compute averaging weights. Furthermore, when averaging ICs with discrepant topographies, one incurs in the risk of mixing sources with significantly different dynamics and loosing the meaning of the ICA source separation.

A limitation of the proposed approach is the bias of the EEG measures toward superficial cortical activity, which possibly precludes the identification of hemodynamic changes associated with deep brain activity. This limitation is however common to all EEGfMRI "integration through prediction" approaches. Nonetheless, the particular syndrome of epilepsy associated with hypothalamic hamartomas is rather well described in the literature, with clear evidence for the hypothalamic hamartoma being the seizure focus, and this region was in fact found in our work in five out of six EEG-fMRI datasets. Although the proposed EEG-BOLD transfer functions are not specific to epileptic activity, this specificity is achieved in the presented methodology by selecting the EEG topography used to extract the metric based on an ICA procedure. Nevertheless, a possible limitation of the proposed transfer functions is that they are not specific to the start of the seizure, and hence to the epileptogenic focus. However, our aim was the description of the full epileptic network, leaving the specific identification of the seizure focus and seizure propagation dynamics for other, related lines of research (Murta et al., 2012).

The study presented here focused on data from a single patient with the aim of providing a proof of concept for the potential usefulness of the proposed methodology for the identification of epileptic networks. The choice of this patient was based on the quality of the EEG data that could be achieved due to the absence of movement associated with the beginning of the

#### **REFERENCES**


in dissociation. *Brain Res. Rev.* 62, 233–244.


seizures. Moreover, epilepsy cases associated with hypothalamic hamartomas are relatively stereotypical in terms of the electrophysiological patterns of seizure propagation, which makes the interpretation of the results relatively more robust in comparison to other types of epilepsy. In general, the clinical utility of EEGfMRI is yet to be established. Nevertheless, in this case study, the results of the EEG-fMRI investigation indicate possible alternative treatment approaches involving the surgical interruption of the seizure propagation pathways (Leal et al., 2009; Murta et al., 2012). The proposed methodology should now be applied to an extended group of patients, in order to generalize the proof of concept presented here and to further validate it. Reproducibility, specificity, and sensitivity should then also be addressed.

In conclusion, we presented a new approach for EEG-fMRI integration in the field of epilepsy, which incorporates and tests different models of the transfer function between EEG and BOLD signals, hence allowing better predictions of the hemodynamic changes associated with epileptic activity. This work therefore provides a contribution to our understanding of the link between EEG and BOLD signals as well as for improving the yield of EEG-fMRI studies in epilepsy.

#### **ACKNOWLEDGMENTS**

We acknowledge the Portuguese Science Foundation (FCT) for financial support through Project PTDC/SAU-BEB/65977/2006, Project PTDC/SAU-ENB/112294/2009, and Project PEst-OE/EEI/ LA0009/2011.

activity by simultaneous electroencephalography and functional magnetic resonance imaging. *Brain* 134, 2867–2886.


linear registration and motion correction of brain images. *Neuroimage* 17, 825–841.


memory task: modulations in low and high frequency bands. *PLoS One* 5:e10298. doi: 10.1371/journal.pone.0010298


Continuous EEG source imaging enhances analysis of EEG-fMRI in focal epilepsy. *Neuroimage* 49, 3219–3229.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 16 October 2012; accepted: 04 January 2013; published online: 25 January 2013.*

*Citation: Leite M, Leal A and Figueiredo P (2013) Transfer function between EEG and BOLD signals of epileptic activity. Front. Neur. 4:1. doi: 10.3389/fneur.2013.00001*

*This article was submitted to Frontiers in Epilepsy, a specialty of Frontiers in Neurology.*

*Copyright © 2013 Leite, Leal and Figueiredo. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*

## Potential use and challenges of functional connectivity mapping in intractable epilepsy

**Robert Todd Constable1,2,3,4\*, Dustin Scheinost 1,3, Emily S. Finn<sup>4</sup> , Xilin Shen<sup>1</sup> , Michelle Hampson<sup>1</sup> , F. ScottWinstanley <sup>2</sup> , Dennis D. Spencer <sup>2</sup> and Xenophon Papademetris 1,3**

<sup>1</sup> Department of Diagnostic Radiology, Yale School of Medicine, New Haven, CT, USA

<sup>2</sup> Department of Neurosurgery, Yale School of Medicine, New Haven, CT, USA

<sup>3</sup> Department of Biomedical Engineering, Yale University, New Haven, CT, USA

4 Interdepartmental Neuroscience Program, Yale University, New Haven, CT, USA

#### **Edited by:**

Mark Holmes, University of Washington, USA

#### **Reviewed by:**

Kaspar Schindler, University of Bern, Switzerland Andreas Schulze-Bonhage, Universitätsklinikum Freiburg, Deutschland

#### **\*Correspondence:**

Robert Todd Constable, School of Medicine, Yale University, 300 Cedar Street, The Anlyan Center, N132, New Haven, CT 06520, USA. e-mail: todd.constable@yale.edu

This review focuses on the use of resting-state functional magnetic resonance imaging data to assess functional connectivity in the human brain and its application in intractable epilepsy. This approach has the potential to predict outcomes for a given surgical procedure based on the pre-surgical functional organization of the brain. Functional connectivity can also identify cortical regions that are organized differently in epilepsy patients either as a direct function of the disease or through indirect compensatory responses. Functional connectivity mapping may help identify epileptogenic tissue, whether this is a single focal location or a network of seizure-generating tissues. This review covers the basics of connectivity analysis and discusses particular issues associated with analyzing such data. These issues include how to define nodes, as well as differences between connectivity analyses of individual nodes, groups of nodes, and whole-brain assessment at the voxel level. The need for arbitrary thresholds in some connectivity analyses is discussed and a solution to this problem is reviewed. Overall, functional connectivity analysis is becoming an important tool for assessing functional brain organization in epilepsy.

**Keywords: epilepsy, functional connectivity, connectome, fMRI, network theory, graph theory, surgical planning**

#### **INTRODUCTION**

Functional connectivity in human neuroscience refers to the synchrony of activity in anatomically distinct regions of the brain: if two areas are highly correlated in their activity over time, they are considered functionally connected. As measured by fMRI, functional connectivity relies on the blood oxygenation level dependent (BOLD) contrast mechanism (Ogawa et al., 1990, 1992), the same as that used in traditional task-based functional MRI studies. But, rather than examining changes in response to specific input stimuli in a block design or event-related paradigm, connectivity mapping can extract information from correlations in the fMRI time-course data while the subject is at rest, in the absence of any externally imposed task. Acquiring connectivity data in the resting state makes connectivity analysis easily adaptable to clinical scanning as it requires no subject participation other than to remain still in the magnet and is thus just like any clinical imaging study. It can also be used intraoperatively in anesthetized patients.

Task-based fMRI provides exquisite maps of functional regions differentially involved in the execution of a specific task and is the foundation of most functional brain imaging research. Task-based fMRI as it is used in surgical planning in epilepsy is chiefly focused on identifying eloquent cortex that must be spared in a surgical intervention,while other techniques are used to identify the epileptogenic tissue to be resected. Task-based fMRI is not well suited for identifying abnormal cortical or subcortical tissue throughout the brain, because only a very small number of regions typically show differential activation to any given task and the number of tasks that can be run in a reasonable time with sufficient statistical power is low. Task-based fMRI has not seen widespread clinical use (outside of mapping for surgical planning) because of limitations on the number of tasks, difficulties associated with presenting tasks and/or training subjects on the task in a busy clinical MR center, and the lack of whole-brain assessment from such studies.

Functional connectivity data is different from task-based fMRI data in that it does not provide information as to which areas of the brain are differentially involved in the execution of a specific task, but instead provides a more basic measure reflecting how different brain areas are functionally connected to one another. As such, functional connectivity studies are potentially more appropriate for whole-brain surveys of functional abnormalities such as the clinical challenge in epilepsy of identifying seizure-generating tissue elements or networks. The hope is that examining the network properties of the seizure-generating tissue will lead to a better understanding of the impact of epilepsy on the functional organization of the brain, assist in understanding comorbidities, and facilitate surgical planning, which in turn could lead to improved surgical outcome through identification and resection of the critical node(s) in the seizure-generating networks(s). Connectivity data could, for example, be used to identify brain areas to target with invasive intracranial recording electrodes. Resting-state fMRI is not limited to identifying tissue to be resected but also has the potential to predict cognitive change following resection if the resection plan involves specific nodes in functional networks. Research in all of these areas is ongoing and these topics are discussed in further detail in the sections that follow.

Functional connectivity mapping was first described by Biswal et al. (1995) in a seminal paper, revealing correlations between brain regions that are synchronized in time, through spontaneous fluctuations of activity, while the brain regions themselves may be quite distant. Brain regions with high resting-state temporal correlations are thought to be involved in the same network and the more well-formed the network, the stronger the connectivity. Spontaneous fluctuations lead to changes in the BOLD signal and most of the work to date has focused on very low-frequency (<0.1 Hz) fluctuations.

Spontaneous fluctuations have also been studied with both surface EEF and intracranial recordings (Towle et al., 1999;Andrzejak et al., 2001, 2006; Schevon et al., 2007; Ortega et al., 2008a,b; Lehnertz et al., 2011; Palmigiano et al., 2012) and local differences in the connectivity between nodes in the seizure focus area have been observed. For a review discussing both increases and decreases in synchronization and their role in epilepsy please see Jiruska et al. (2013). The EEG work, it should be noted, has primarily focused on higher frequency fluctuations in the 1–40 Hz range although several groups now are beginning to target much lower frequency oscillations.

After the initial paper by Biswal et al. (1995), much effort was focused over the next decade on verifying that the correlations observed between brain regions indeed reflect true functional connections and not simply noise correlations due to physiological fluctuations. Without doubt, care must be taken to minimize the effects of physiologic noise as there are clear correlations associated simply with respiratory and cardiac signal fluctuations; however, the current consensus is that meaningful functional connections are found in continuously recorded BOLD fMRI data. Since about 2005, the field of functional connectivity mapping has been expanding rapidly through novel approaches to analysis, noise removal methodology, and applications to clinical populations and basic neuroscience problems. In particular, the field is expanding to more in-depth analyses that move beyond simply examining correlations between two or a small handful of brain regions to capturing network information from a large array of nodes across the brain.

For most functional connectivity studies, BOLD fMRI data is collected with the subject in the resting state with eyes either open or closed, and the patient is not required to perform a task. It is also possible to obtain connectivity data in the presence of a task and/or after brain-state manipulations, but it should be noted that the connectivity patterns can be slightly modified by task or brain state. Functional connectivity patterns have been shown to be sensitive to brain state as well as behavioral variables (Hampson et al., 2006a,b; Johnson et al., 2006; Rogers et al., 2010; Bonelli et al., 2012; Cole et al., 2012), and to vary with development (Fair et al., 2008; Schafer et al., 2009; Myers et al., 2010) and age (Dosenbach et al., 2010; Hampson et al., 2012).

Connectivity changes have also been reported in several clinical populations (Quigley et al., 2001; Lowe et al., 2002; Irwin et al., 2004; Saini et al., 2004; Haas et al., 2006; Waites et al., 2006; Hoffman et al., 2007; Schafer et al., 2009; Wang et al., 2009; Freilich and Gaillard, 2010;Myers et al., 2010;Bai et al., 2011;Killory et al., 2011; Zhang et al., 2011a,b; Bagshaw and Cavanna, 2012; de Groot et al., 2012) and there are now more than 100 publications related to

connectivity measures in epilepsy patients. There is evidence that correlations between time-varying BOLD signals reflect intrinsic functional connections in that they are present when subjects are both awake and under anesthesia (Vincent et al., 2007; Martuzzi et al., 2010, 2011) and they are highly reproducible (Shehzad et al., 2009). Overall, resting-state functional connectivity mapping has significant potential to reveal the functional organization of the brain and how it may be altered in different diseases or disorders.

This review begins with an introduction to resting-state functional connectivity, including basic data collection and processing steps, the type of information that can be obtained, and various means of application. The particulars associated with restingstate methodology as applied in epilepsy are considered including strategies for voxel-based and region-of-interest (ROI) based analyses. Issues associated with selecting a connectivity threshold or avoiding such thresholds, ROIs, and the sensitivity of the method to choice of ROI are considered in addition to the emerging field of connectivity-based parcelation for identifying minimal functional subunits for nodal analysis. Finally we provide a number of early clinical results reflecting the potential of functional connectivity data to contribute to the clinical management of epilepsy.

#### **THE BASICS OF FUNCTIONAL CONNECTIVITY DATA COLLECTION**

Resting-state data is typically collected using a gradient-echo echo planar imaging (EPI) pulse sequence that provides whole-brain coverage with a temporal resolution of 3 s or less at a field strength of 3 T. Current state-of-the-art EPI involves the use of a multiband/multi-plexed EPI sequences (Feinberg et al., 2010) that can provide whole-brain coverage with 2 mm<sup>3</sup> voxel dimensions and a repetition time (TR) of less than 1 s. Short TR in fMRI generally provides better statistical power (Constable and Spencer, 2001) and eases removal of physiological noise in connectivity data (Lowe et al., 1998). A minimum of approximately 5 min of such data is required to obtain connectivity maps with reasonable signal-to-noise ratios (SNRs). Such an acquisition easily fits within the constraints of a clinical MRI study, although in general acquiring more data is advantageous.

Initial post-processing steps involve motion correcting the data (to ensure all volume acquisitions are aligned through time), and removing time-course signals that may be of no interest including the mean global signal through time, components associated with cerebrospinal fluid and white-matter signals, and the variables describing subject motion (the three translation and three rotation directions). Temporal drift terms are often regressed from the data and a low pass filter (<0.1 Hz) is applied.

The extent to which regressing the temporal fluctuations associated with motion is effective has been a topic of much discussion in the literature (Weissenbacher et al., 2009) and there is evidence for residual motion effects after correction (Power et al., 2012; Satterthwaite et al., 2012, 2013; Van Dijk et al., 2012). The issue of whether or not to remove the global signal mean also remains a hot topic as this can lead to the introduction of negative correlations in the data, making interpretation difficult (Fox et al., 2009; Murphy et al., 2009; Hampson et al., 2010).

Following the steps in the chosen preprocessing pipeline, functional connectivity is measured as the temporal correlation between the signal from any pair of voxels or ROIs. The following sections discuss the pros and cons of each approach in detail.

#### **CONNECTIVITY MEASURES**

Resting-state connectivity mapping has been used to examine functional connections between cortical regions since the first presentation of the method by Biswal et al. (1995). In the sections that follow, we discuss how connectivity analyses can be applied to epilepsy for examining alterations in normal brain networks as a function of disease, or as a function of surgical intervention, as well as how connectivity mapping can be used to detect segments of tissue with abnormal functional connectivity with the goal of identifying the seizure-generating foci or network.

#### **ROI-TO-WHOLE-BRAIN**

Defining a ROI and performing ROI-to-whole-brain connectivity analysis is probably the most common approach to examining connectivity in the brain; this was the original method introduced by Biswal et al. (1995). Such an approach is motivated when an investigator is interested in a particular brain region (the seed ROI) and wishes to examine what other brain regions the seed is connected to, as well as how such connections vary between healthy controls and a patient group, or preand post-surgical intervention (see **Figure 1** for an example). In epilepsy this approach has been used to examine changes in language networks before and after anterior temporal lobe resection (Bonelli et al., 2012). In another study, Pereira et al. (2010) examined connectivity between the left and right hippocampi in mesial temporal sclerosis patients and found that relative to control subjects, patients with left hippocampal sclerosis showed a larger decrease in functional connectivity than patients with right hippocampal sclerosis in subjects with left-hemispheric language dominance. Morgan et al. (2012) performed a similar study using each of the left and right hippocampi as seeds for a seed-to-whole-brain analysis. This work revealed that the connectivity between the right hippocampus and the ventral lateral nucleus of the right thalamus could distinguish between seizurefree patients with left temporal lobe epilepsy (TLE) and right TLE patients and that in general, connectivity was greater in the seizure-free patient group with left TLE compared to the healthy controls.

In another study of the medial temporal lobe, Pittau et al. (2012) selected four manually drawn ROIs and examined differences in connectivity between 23 patients with mesial TLE and compared these to 23 age- and gender-matched controls in an ROI-to-whole-brain analysis. They found that patients with right MTLE had decreased connectivity between right amygdala and right hippocampus and the brain areas associated with the default mode network, some prefrontal regions, and contralateral mesial temporal structures. The left MTLE patients showed decreased connectivity between the left amygdala and left hippocampus to the default mode network, contralateral hippocampus, and bilateral limbic prefrontal regions.

Moving outside the temporal lobe, a study by Killory et al. (2011) performed ROI-to-whole-brain analysis using eight ROIs: three defined used a functional task-based paradigm and another five spherical ROIs based on previous literature. They then compared the resulting maps from patients with childhood absence epilepsy to those of control subjects. An overall decrease in connectivity in the attention nodes was observed in the absence patients, who demonstrated impaired connectivity in the insula/frontal operculum and medial frontal nodes relative to healthy control subjects.

In a study that used simultaneous EEG/fMRI to define the initial seed regions,Negishi et al. (2011) compared connectivity to the seed within the same hemisphere and the contralateral hemisphere to derive a connectivity laterality index. Results showed decreased laterality of functional connectivity in the patients that were seizure-free after surgery compared to those that had recurrent seizures after surgery.

#### **ROI-TO-ROI**

Another approach to connectivity analysis is to define a set of ROIs and to measure functional connectivity between all possible pairs of ROIs. This has been applied in numerous studies of epilepsy and some of these are summarized below. Once multiple ROIs are considered, network properties for these ROIs can be measured as discussed in the section on network theoretic measures below.

In a study by Bettus et al. (2010), the investigators chose five ROIs in each hemisphere – primarily in medial temporal lobe

**FIGURE 1 | Seed-to-whole-brain connectivity mapping**. A reference BOLD signal time-course from a seed in Broca's area (top row) is correlated with the BOLD time-course for all gray-matter voxels in the image, revealing areas to which the seed is functionally connected (bottom row; hot colors = strong positive correlation, cool colors = weak or negative correlation). Such

seed-to-whole-brain connectivity maps may be examined for a single epilepsy patient relative to a group of healthy control subjects, or for a group of patients with a similar pathology to a group of healthy control subjects, or correlated with behavioral measures. This approach and the results it produces are highly dependent upon how the initial seed region is chosen.

regions, defined using the Pick atlas tool in SPM (Maldjian et al., 2003) – and examined connectivity between each pair of ROIs and compared this to data from 36 control subjects. They reported decreased connectivity in mesial TLE patients relative to control subjects, with the largest decreases on the ipsilateral side and some increased connectivity on the contralateral side. In a similar study, Morgan et al. (2011) defined hippocampal ROIs using both structural and functional information and revealed a relationship between functional connectivity and causal influence of the left and right hippocampi that varied with the duration of disease in a group of 19 mesial TLE patients. Looking at a broader network of nodes in a subsequent paper, Morgan et al. (2012) compared pre-operative fMRI connectivity data in two groups of patients who were postoperatively sorted into seizure-free and those with recurring seizures. They showed that the connectivity between the right hippocampus and the ventral lateral nucleus of the right thalamus could distinguish between seizure-free patients with left TLE from those with right TLE with high sensitivity and specificity. The patients with recurring seizures after surgery generally exhibited different connectivity values in this network from those that were seizure-free.

In a study of 11 children with intractable epilepsy, Widjaja et al. (2013) used independent components analysis to identify the default mode network and then calculated the functional connectivity between different pairs of ROIs in this network. They found reduced connectivity in the default mode network in children with medically refractory epilepsy.

Using ROI-to-ROI connectivity analyses with a number of predefined ROIs, Bai et al. (2011)investigated functional connectivity in a group of childhood absence epilepsy patients. In a unique twist on the ROI approach, they also examined ROI homologs in a cross-hemispheric study where the ROIs in this case were individual voxels. The primary finding was increased interhemispheric connectivity in the lateral orbitofrontal cortex in the patient group relative to the healthy control subjects.

A pattern classification approach was applied to another ROIto-ROI study in epilepsy by Zhang et al. (2011a). This work examined connectivity patterns between nodes in both hemispheres and found, much like the Negishi et al. (2011) study described in the previous section, that asymmetry in functional connectivity could correctly distinguish epilepsy patients from healthy controls with 82.5% specificity and 85% sensitivity.

#### **NETWORK THEORETIC MEASURES**

While ROI-to-whole-brain and ROI-to-ROI analyses can reveal much about functional connectivity, such approaches do not take the next step of considering nodes in the context of a functional network and the properties of that network as a whole. Recently there has been an explosion of interest in applying network theory (see Achard et al., 2006; Bullmore and Sporns, 2009; Bressler and Menon, 2010; Hagmann et al., 2010; Rubinov and Sporns, 2010) to the analysis of functional connectivity data in order to characterize brain connections at both the nodal and network level. This allows for the observation of the networks associated with specific brain functions and generally moves fMRI from identification of individual nodes to systems involved in the execution of a task. Such theory, when applied to resting-state fMRI data, can characterize the topology of normal networks in the brain and, by extension, identify abnormal patterns of connectivity.

Network theory measures have been applied in epilepsy to EEG or intracranial EEG data (Ponten et al., 2009, 2010; van Dellen et al., 2009; Varotto et al., 2012), MEG data (Chavez et al., 2010), and cortical thickness correlations (Bernhardt et al., 2011), but to our knowledge only one published study to date has specifically examined network properties using fMRI data (Zhang et al., 2011a; although in the voxel-based approaches section below we include a small study by Stufflebeam that uses the graph-theory measure of degree to attempt to identify epileptogenic tissue). In this work, a total of 36 ROIs were defined using both functional (language and motor regions) and anatomic definitions [Brodmann's (1909) areas]. Five network measures – degree, strength, clustering coefficient, closeness, and betweenness centrality (Wasserman and Faust, 1994) – were calculated for these nodes in individual subjects as well as for input into a classification strategy to determine if such measures could distinguish medial TLE patients from healthy control subjects. The results confirmed that network measures such as these can indeed aid in classifying epilepsy patients and that the epilepsy process is associated with changes in network-level functional brain organization.

Taken together, independent of whether the analysis is ROI-towhole-brain, ROI-to-ROI, or at the level of network properties, these studies suggest that functional connectivity has potential for identifying disrupted circuitry as a function of disease and, perhaps more importantly, predicting outcomes from surgical intervention. In the sections that follow we outline some important issues encountered in these studies related to defining ROIs and determining thresholds for connectivity.

#### **ROI-BASED APPROACHES AND THE PROBLEM OF ROI DEFINITION**

As described above, many connectivity analyses require the *a priori* definition of at least one ROI. While not often highlighted in previous publications, the choice of the seed ROI(s) and how the exact boundaries of that ROI are defined is critical. If an ROI contains multiple time-courses then the average time-course from the ROI may not properly represent any of the time-courses within an ROI and the results may be completely erroneous. Further, varying the spatial definition of the seed can substantial changes in results. This is easily highlighted by observing that in a typical ROI-to-whole-brain connectivity map (e.g., see **Figure 1**), there are often sharp transitions from positive to negative correlations; hence, moving the seed can result in a very different map.

As is evident in the papers already discussed, numerous approaches to defining ROIs have been used to date. Task-based fMRI has been used to define specific functional circuits within which connectivity can be analyzed (Frings et al., 2009; Bonelli et al., 2012). This approach, however, suffers from the limitations of task-based fMRI studies in general in that only a very limited number of ROIs are activated by a task and thus wholebrain assessment of connectivity is not possible using such definitions. Another approach has been to use independent component analyses (ICA) (McKeown et al., 1998) to delineate brain regions (Luo et al., 2012; Mankinen et al., 2012) but these have typically identified only a very limited number of networks – often, for example, fewer than 10.

Anatomic ROI definitions have also been used extensively (Crespo-Facorro et al., 1999; Tzourio-Mazoyer et al., 2002; Makris et al., 2005; Shattuck et al., 2008; Zhang et al., 2011a). Such definitions are ideal in structures that are well-defined anatomically (such as the hippocampus) but are difficult in areas such as the frontal and parietal cortices, and therefore the risk of mixing temporal signals into heterogeneous ROIs in these regions is high. Many investigators have used small spherical ROI placements (Shehzad et al., 2009; Bai et al., 2011; Bettus et al., 2011; Killory et al., 2011; Koyama et al., 2011); in this case the risk of mixing different functional time-courses decreases with the size of the defined sphere but is not eliminated. Many investigators have arbitrarily parcelated the cortex into anywhere from 100 to 1000 nodes, but again, with such an approach, the node definitions may not necessarily reflect true functional boundaries.

An emerging area of investigation involves performing wholebrain parcelation based on the time-courses themselves (van den Heuvel et al., 2008; Shen et al., 2010, 2013; Craddock et al., 2012). This approach appears very promising because it can provide minimal functional subunits with uniform time-courses within each unit. An example obtained using the approach of Shen et al. (2013) is shown in **Figure 2**. Both Shen et al. (2013) and Craddock et al. (2012) have shown that ROIs extracted from these parcelations had higher functional homogeneity than anatomically defined ROIs and thus were more relevant for fMRI connectivity analyses. This parcelation approach using connectivity data itself appears to solve the problem of providing whole-brain ROI definitions for meaningful connectivity analysis. The next problem is how to apply such an approach to a patient population or to a group of patients. For example, if one generates a parcelation from healthy control subjects to investigate differences in network properties between control subjects and epilepsy patients, any results could be interpreted as due to actual differences in network properties, or a mismatch in the ideal functional boundaries for the parcelation nodes that come from imposing a control-derived parcelation on patients. This latter question can be addressed by directly comparing a parcelation derived from the patient or patient group to the parcelation derived from the healthy control subjects.

Thus, while it has been difficult to define ROIs for functional connectivity analysis, these new parcelation approaches appear to have solved this problem. Another approach, however, is to avoid the ROI problem completely by moving to voxel-level connectivity analysis, in which each voxel in the gray matter is treated as an individual node and summary statistics on the network or connectivity properties of each voxel are obtained. This approach is described in the next section.

#### **VOXEL-BASED CONNECTIVITY ANALYSIS**

Voxel-level network analyses have been developed that can provide insight into the functional connectivity of individual tissue elements, and a number of important studies have been published using such approaches (Buckner et al., 2009; Martuzzi et al., 2011; Stufflebeam et al., 2011; Scheinost et al., 2012). Voxel-level analyses have the benefit of not requiring *a priori* definition of a ROI, instead treating each voxel as a node in a network analysis. In this approach a network measure such as *degree* can be calculated for each voxel. (*Degree* is the number of connections to a voxel above some arbitrary correlation threshold, i.e., *r* > *t*, where *t* = 0.25). Such a degree map, as shown in **Figure 3** below, provides then for the first time a gray-scale contrast reflecting the functional connectivity of each tissue element and is a potentially powerful approach to identifying regions that have abnormal functional connectivity on a whole-brain level. Such degree maps can also be used to define nodes for further ROI-to-ROI network analyses and/or to compare individual patients with control-group data.

Such a voxel-based approach has been used to study brain changes in Alzheimer's disease (Buckner et al., 2009), to examine the effects of anesthetics on the human brain (Martuzzi et al.,

**FIGURE 2 | A map of reproducible (across 79 subjects) functional subunits identified using resting-state connectivity data**. These functional subunits are ideal nodes for connectivity analyses as they have highly uniform time-courses for each voxel within a given node by definition. This parcelation approach (Shen et al., 2010, 2013) provides a solution to the ROI definition problem in connectivity or network analysis of the brain.

**subjects) with contrast reflecting the functional connectivity of each voxel as measured by the network measure of degree (brighter colors**

and compared to control-group data to isolate tissue elements with abnormal functional connectivity.

2011) and to identify epileptogenic tissue in epilepsy (Stufflebeam et al., 2011).

There are multiple approaches to calculating degree; three simple schematics are shown in **Figure 4** to illustrate the different approaches. In these examples the voxel-based degree measure can be calculated for connections encompassing the entire brain (whole-brain connectivity), within the same hemisphere as the voxel (ipsilateral connectivity) or spanning connections to the other hemisphere (contralateral connectivity). In each case the calculation is performed for each gray-matter voxel in the brain and the intensity of each voxel then reflects its number of connections above a predetermined correlation threshold.

An example of such an approach applied in epilepsy is shown in **Figure 5** below where the three *degree* measures are shown for a single patient relative to 20 healthy control subjects. The idea behind this approach is that seizure-generating tissue may have altered functional connectivity either because of the epilepsy processes themselves or because the tissue is not functioning normally. Generally, in all three measures, we observe widespread decreases along with some increases in functional connectivity for the patient relative to the controls. Much needs to be learned about this approach, however, as there are subtle differences in the ipsilateral and contralateral measures that may provide relevant information for specific pathologies. In the case shown in **Figure 5** the patient had suspected right medial TLE and indeed there is a large region of decreased functional connectivity in the right temporal lobe relative to healthy control subjects across all three measures. However, there are many other regions that also show differences in functional connectivity relative to control subjects and it is currently an open question whether these differences reflect part of the epilepsy network, brain reorganization of functional subunits as a result of having epilepsy, or some other neurophysiological mechanism.

Given some clinical consensus via more conventional measures (surface and/or video EEG, other clinical data, and/or invasive recordings) there is often some indication that the patient has foci in a particular region or lobe. In the next example, shown in **Figure 6** below, the right hippocampus was suspected to be part of the seizure-generating network and indeed the ipsilateral degree measure showed decreased functional connectivity in this region for the patient (see **Figure 6A**). (As an aside, we note that the ipsilateral degree measure is perhaps the most straightforward and easiest to interpret, whereas the whole-brain and contralateral measures are more difficult to understand in terms of the lateralization of the source of the problem.) As in the previous example, many other regions also show alterations in connectivity and some of these may be part of the seizure-generating network. Still, with *a priori* clinical information about the suspected site of seizure onset, it is reasonable to use the right hippocampus as a seed ROI for an ROI-to-whole-brain analysis. The nodes connected to

**FIGURE 5 | Functional connectivity difference maps for a single epilepsy patient versus a group of healthy control subjects (red** = **connectivity increased in patient relative to controls; blue** = **connectivity decreased in patient relative to controls).** Some of these regions are consistent with the

interactions with different pathological variants. increased further by examining the distribution of degree values across the entire range of thresholds; this also addresses the second

seizure onset zones but clearly a large number of regions show altered connectivity. There are subtle differences between these three measures and more work is needed to determine the sources of these differences and the

this right hippocampal region are shown in **Figure 6B** with hot colors reflecting strong connectivity to the seed region. Combining the ipsilateral connectivity map in **Figure 6A** with the map showing connections to the suspected focal region in **Figure 6B** using a logical AND operation, yields a much more circumscribed map **Figure 6C**, which highlights regions that have both altered functional connectivity and are part of the same network. (The most superior slices in the brain are not shown in the bottom two panels as there are no regions in these slices that satisfy both constraints.)

Yet, despite their potential, voxel-based analyses of connectivity have not seen widespread application primarily for two reasons: (1) a lack of sensitivity, and (2) vulnerability to threshold effects.

The first problem is partially solved by higher-field magnets of 3 T or more as well as the more recent move to ultrafast pulse sequences such as the multi-band/multi-plexed sequence developed by Feinberg et al. (2010). In addition, sensitivity can be problem, as described in the next section. The second problem, the need for an arbitrary threshold to calculate the degree measure for each voxel, is unique to the application of network theory measures to fMRI data. In most network theory applications, such as tracking the transmission of a disease or analyzing friend links on Facebook, the decision tree is binary – either the disease is transferred or it is not; either two people are friends or they are not. However, in the case of fMRI connectivity, since connectivity is measured as a correlation, values may span a continuum from −1 to 1. Thus, an arbitrary threshold is typically invoked to decide if two regions are connected or not. In general there is no principled way to choose this threshold and the results can change dramatically at different thresholds, which is a major issue for applying network or graph-theory measures (Rubinov and Sporns, 2010) to fMRI data. To get around this

**FIGURE 6 | (A)** Ipsilateral connectivity for a single patient versus a group of healthy control subjects showing multiple regions with altered connectivity (warm colors = higher connectivity for patient; cool colors = higher connectivity for controls). **(B)** In an individual analysis of the patient data, selecting the right hippocampus, which was the onset zone suspected from clinical data, as a seed region reveals a network of nodes with high vs. low connectivity to the right hippocampus. In **(C)** the intersection of the maps in **(A)** and **(B)** highlights regions that have abnormal connectivity compared to controls AND are part of the same network involving the right hippocampus.

arbitrary threshold choice, Scheinost et al. (2012) presented a new approach that examines the entire connectivity distribution curve across all correlation thresholds from 0 to 1. (The absolute value of the correlations can be taken, therefore taking negative correlations into account, or the positive and negative correlations can be treated in separate analyses). This new approach, referred to as the intrinsic connectivity distribution, is described in the section that follows.

**FIGURE 7 | Regions of high intrinsic connectivity for left temporal lobe epilepsy (LTE) patients (n** = **10) shown at different statistical levels**. At p < 0.05 corrected (top row) the ICD map (left) highlights significant patterns of increased connectivity in the left hippocampus whereas the degree map (right) does not show significantly different connectivity in this region. At a less stringent statistical threshold (bottom row), the degree map shows the same areas at p < 0.05 uncorrected (right) as the ICD maps (left) at either statistical threshold. This suggests that intrinsic connectivity is capable of detecting abnormal connectivity in regions of epileptogenic tissue and that the ICD approach yields higher sensitivity and reproducibility.

#### **INTRINSIC CONNECTIVITY DISTRIBUTION ANALYSIS**

By characterizing the entire *degree* curve for any of the three voxelbased connectivity metrics described above (whole-brain, ipsilateral, and contralateral), the need for a specific but arbitrary threshold is eliminated, providing both a more consistent measure of functional connectivity and a more sensitive measure of individual differences. Summarizing the entire degree distribution curve captures information about all connections, weak to strong. This parametric estimation approach therefore enables the interpretation of any differences (between an epilepsy patient and a control group, for example) to take the form of "more stronger connections," "more weaker connections,"or"more even spread of connections." An example in a group of 10 TLE patients comparing degree at a single threshold (*r* > 0.25) and ICD is shown in **Figure 7**.

The ICD approach improves on the earlier voxel-based degree measure through increased sensitivity as well as through more stable interpretation of the results. To date no comprehensive study has been published using these emerging methods, but our lab as well as others are currently applying this to the problem of localizing seizure foci or networks in epilepsy and no doubt more studies will appear shortly.

#### **SUMMARY**

In summary, resting-state functional connectivity as measured by BOLD fMRI reflects intrinsic connections in the brain, providing insight into how the brain is wired and how such wiring may be altered in disease or through surgical intervention. Like task-based fMRI, functional connectivity measures can be altered by task, drug or brain state, but unlike task-based fMRI, which reflects small changes in activity superimposed upon a very high baseline activity level, functional connectivity measures intrinsic functional connections that exist even in the presence of deep anesthesia.

The measures of intrinsic connectivity as described above, particularly the voxel-based degree measure, provide for the first time a contrast mechanism in MRI based on function rather than anatomy. Since this approach requires no task, the data is acquired much in the same way that anatomic MRI data is acquired and thus these measures can be easily incorporated into clinical diagnostic radiology departments.

Furthermore, while task-based fMRI provides functional information on only a very limited number of cortical regions (i.e., those few regions differentially activated by the task), the voxel-based connectivity measures described above provide functional information on the whole-brain without the need for a task. This makes it ideal for investigating the effects of pharmacological agents and the impact of diseases and specific pathologies on the brain without the need for *a priori* knowledge or selection of ROI. The use of functional connectivity methods in epilepsy will undoubtedly increase as we learn more about the functional changes that occur with epilepsy and/or as a function of surgical intervention.

#### **ACKNOWLEDGMENTS**

Funding support from the NIH NIBIB and NINDS is gratefully acknowledged.

#### **REFERENCES**


epileptic networks assessed by intracerebral EEG and BOLD signal fluctuations. *PLoS ONE* 6:e20071. doi:10.1371/journal.pone.0020071


to a right hemisphere language center in prematurely born adolescents. *Neuroimage* 51, 1445–1452.


activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. *Neuroimage* 15, 273–289.


temporal lobe epilepsy. *Ann. Neurol.* 59, 335–343.


Altered functional-structural coupling of large-scale brain networks in idiopathic generalized epilepsy. *Brain* 134, 2912–2928.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 21 February 2013; paper pending published: 20 March 2013; accepted: 11 April 2013; published online: 22 May 2013.*

*Citation: Constable RT, Scheinost D, Finn ES, Shen X, Hampson M, Winstanley FS, Spencer DD and Papademetris X (2013) Potential use and challenges of functional connectivity mapping in intractable epilepsy. Front. Neurol. 4:39. doi: 10.3389/fneur.2013.00039*

*This article was submitted to Frontiers in Epilepsy, a specialty of Frontiers in Neurology.*

*Copyright © 2013 Constable, Scheinost, Finn, Shen, Hampson, Winstanley, Spencer and Papademetris. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*

## Local functional connectivity as a pre-surgical tool for seizure focus identification in non-lesion, focal epilepsy

#### **K. E.Weaver 1,2\*,W. A. Chaovalitwongse1,2,3, E. J. Novotny 2,4, A. Poliakov <sup>5</sup> ,T. G. Grabowski 1,2,6 and J. G. Ojemann2,7,8**

<sup>1</sup> Department of Radiology, University of Washington, Seattle, WA, USA


<sup>8</sup> Neurosurgery, Seattle Children's Hospital, Seattle, WA, USA

#### **Edited by:**

Mark Holmes, University of Washington, USA

#### **Reviewed by:**

Sandrine DeRibaupierre, University of Western Ontario, Canada Mario A. Vanegas, Instituto Nacional de Neurologia y Neurocirugia, Mexico

#### **\*Correspondence:**

K. E. Weaver, Department of Radiology, University of Washington, 1959 NE Pacific Street, Box 357115, Seattle, WA 98195, USA. e-mail: weaverk@uw.edu

Successful resection of cortical tissue engendering seizure activity is efficacious for the treatment of refractory, focal epilepsy. The pre-operative localization of the seizure focus is therefore critical to yielding positive, post-operative outcomes. In a small proportion of focal epilepsy patients presenting with normal MRI, identification of the seizure focus is significantly more challenging. We examined the capacity of resting state functional MRI (rsfMRI) to identify the seizure focus in a group of four non-lesion, focal (NLF) epilepsy individuals. We predicted that computing patterns of local functional connectivity in and around the epileptogenic zone combined with a specific reference to the corresponding region within the contralateral hemisphere would reliably predict the location of the seizure focus. We first averaged voxel-wise regional homogeneity (ReHo) across regions of interest (ROIs) from a standardized, probabilistic atlas for each NLF subject as well as 16 ageand gender-matched controls. To examine contralateral effects, we computed a ratio of the mean pair-wise correlations of all voxels within a ROI with the corresponding contralateral region (IntraRegional Connectivity – IRC). For each subject, ROIs were ranked (from lowest to highest) on ReHo, IRC, and the mean of the two values. At the group level, we observed a significant decrease in the rank for ROI harboring the seizure focus for the ReHo rankings as well as for the mean rank. At the individual level, the seizure focus ReHo rank was within bottom 10% lowest ranked ROIs for all four NLF epilepsy patients and three out of the four for the IRC rankings. However, when the two ranks were combined (averaging across ReHo and IRC ranks and scalars), the seizure focus ROI was either the lowest or second lowest ranked ROI for three out of the four epilepsy subjects. This suggests that rsfMRI may serve as an adjunct pre-surgical tool, facilitating the identification of the seizure focus in focal epilepsy.

**Keywords: resting state fMRI, functional connectivity, non-lesion, focal epilepsy, ReHo, contralateral, pre-operative evaluation, epilepsy surgery**

#### **INTRODUCTION**

Current standards of care for the treatment of pharmacoresistant, focal epilepsy includes the surgical resection of epileptogenic cortex. Typically, the tissue targeted for resection encompasses an extended area around the seizure focus believed to be involved in the propagation of epileptiform discharges, generally referred to as the epileptogenic zone (Rosenow and Lüders, 2001; Laufs, 2012). The benefits of epilepsy surgery have clearly been established. Numerous prospective as well as longitudinal studies have shown that higher rates of seizure freedom, improved quality of life, and decreased long-term remission rates are associated with successful surgical intervention (Wiebe et al., 2001; Spencer and Huh, 2008; de Tisi et al., 2011).

The precise localization of the seizure focus and the extended epileptogenic zone is therefore critical to yielding positive, postoperative outcomes. Pre-surgical evaluations aimed at identifying the seizure focus are comprised of any number of interdisciplinary approaches, including electrophysiological investigations [e.g., electroencephalography and less frequently sub-dural electrophysiology such as electrocorticography (ECoG) or stereoelectroencephalography], traditional neuropsychological evaluation, modern structural [e.g., structural magnetic resonance imaging (MRI)], metabolic [e.g., [18F]fluoro-2-deoxy-glucose positron emission tomography (FDG-PET)], and functional imaging based approaches (e.g., functional MRI). The typical clinical evaluation identifies sites of pathology from structural-based MR scans and

probes surrounding tissue for epileptogenic potential using a combination of the aforementioned modalities. However, in a few patients (i.e., approximately 25% of all qualifying surgical candidates), structural imaging is normal (i.e., an absence of qualitative, gross pathology – Duncan, 2010). In these non-lesional cases, seizure localization presents an additional challenge and clinicians must rely more heavily on alternative approaches (Siegel et al., 2001; Jayakar et al., 2008).

[18F]fluoro-2-deoxy-glucose positron emission tomography, which has traditionally been a widely used pre-surgical evaluative tool, plays a particularly important role in the absence of identified structural abnormalities (Mauguière and Ryvlin, 2004). In cases of refractory, non-lesional epilepsy, identification of a focal area of hypometabolism may reflect candidate seizure focus sites. It is not uncommon however to find hypometabolic regions outside the suspected region of interest (ROI). Thus, FDG-PET hypometabolic regions are frequently used to guide ECoG recordings. During these studies it is often noted that the extent of abnormal hypometabolic regions overlaps with the ictal onset zones and in many cases these areas are substantially larger than and overlap with electrodes displaying interictal epileptic discharges (IEDs) (Duncan, 2010). Moreover, overlapping sites of hypometabolism are commonly lateralized to one hemisphere. For example, when classified by seizure-freedom rates at a 12-month follow-up, quantitative comparisons of FDG uptake rates of the hypometabolic regions relative to the contralateral side showed high accuracy (∼80%) in identifying the hemisphere harboring the epileptogenic focus (Won et al., 1999).

In recent clinical research studies, fMRI has been shown to be a reliable complementary study to FDG-PET. For example, resected cortex displaying pre-operative evoked BOLD signal activations highly concordant with simultaneously EEG recorded IEDs was associated with a greater probability of post-operative seizure freedom (Thornton et al., 2010). It was noted that the greater the degree of overlap between resected tissue and the spread of IED correlated BOLD signal across a region, the greater the probability of long-term seizure freedom. Based on this good concordance, the authors suggested that use of simultaneously acquired EEGfMRI maybe"a useful adjunct"during the pre-operative evaluation of epileptogenic cortex, particularly in the absence of identified pathology (Zijlmans et al., 2007; Thornton et al., 2010). Despite the major advantages of simultaneous EEG-fMRI during preoperative evaluation, it is not readily available in the clinical setting.

One promising application of BOLD fMRI that may aid seizure focus localization and is now commonly available in the clinical setting is resting state functional MRI (rsfMRI) functional connectivity (fc) (Fox and Raichle, 2007; Biswal et al., 2010). This method calculates whole-brain voxel-wise correlations of infra-slow (<0.1 Hz) BOLD signal fluctuations extracted during a resting period and depicts them as maps of brain connectivity. rsfMRI has been used extensively to reveal patterns of fc across and between large-scale neural networks (Damoiseaux et al., 2006). These patterns of correlations are believed to reflect an underlying dynamic but intrinsic neural architecture (Honey et al., 2009; Keller et al., 2011) driven by direct (e.g., mono-synaptic) and/or indirect (poly-synaptic) anatomical connectivity (Biswal

et al., 2010). Many proposed applications have capitalized on the inherent advantages of rsfMRI. For instance, rsfMRI has been to shown to identify intact language networks in the absence of verbal responses (Shimony et al., 2009). In MTLE patients, rsfMRI has revealed disrupted fc across regions commonly involved in the greater epilepsy network, primarily on the ipsilateral side to the seizure focus. Interestingly, increased fc within contralateral regions was also observed suggestive a possible cross-hemisphere compensatory mechanism (Bettus et al., 2009). As a follow-up investigation, the same group of investigators reported that fc increases observed contralateral to MTL pathology lead to high degree of specificity (>91%) for identification of the hemisphere that houses the seizure focus (Bettus et al., 2010).

More recent developments in rsfMRI methodology have begun to focus on patterns of connectivity specific to the local cortical environment (Zang et al., 2004). That is, measures of local connectivity mapping correlations restricted to a finite set of voxels within a ROI. One such method that has recently gained some popularity is Regional Homogeneity (ReHo), a technique that calculates a non-parametric cross-correlation coefficient between the timeseries of a center voxel with a local cluster of voxels of pre-defined sized (Zang et al., 2004; Zhong et al., 2011). To date, reports applying ReHo for seizure focus localization have not been published. A few studies have contrasted ReHo in epilepsy patients relative to control volunteers, observing for example significantly higher thalamic ReHo in a group of generalized tonic-clonic epilepsy patients, values that were negatively correlated with epilepsy duration (Zhong et al., 2011). The anatomical assumptions underlying local fc are built upon patterns of cortico-cortical connectivity. Variability across local cortical neighborhoods or "small-world networks" are therefore assumed to reflect weighted differences of connectivity across neighboring neuronal units (He et al., 2007; Bullmore and Sporns, 2009) leading to the concept of scale-free network properties inherent to the brain's innate architecture (Barabási and Albert, 1999).

While epileptogenic mechanisms and the underlying etiologies are widely variable in patients with focal, treatment-resistant epilepsy, it is well established that the aberrant nature of prolonged epileptic discharges lead to significant neuroanatomical alterations particularly within the epileptogenic zone (Thom, 2004). Animal models and neuropathological reports of resected human epileptogenic tissue have revealed that prolonged seizure activity results in (among many other well-established biochemical and pathological effects) significant neuronal injury and necrosis within the seizure network, particularly within neocortical pyramidal cells (Sankar et al., 1998; Chen and Wasterlain, 2006). Further, a wealth of animal studies has concluded that persistent seizures activity can lead to significant dendritic damage including alterations in spin morphology and an overall down regulation of dendritic spines (Multani et al., 1994; Wong and Guo, in press).

Our overall aim is to examine the capacity of rsfMRI local connectivity to serve as a useful adjunct in the pre-operative evaluation process of seizure focus localization. Based on the extent literature, we hypothesized that local fc in and around the seizure focus in patients with non-lesion, focal (NLF) epilepsy would be significantly lower relative to (1) controls, (2) the corresponding region within the contralateral hemisphere, and (3) ipsilateral ROIs outside of the epileptogenic zone. We choose a two-step analysis approach. First we examined fc in and around the seizure focus. To accomplish this we calculated whole-brain ReHo and averaged across different ROIs.We then tailored a more traditional fc approach to specifically contrast local fc at the seizure focus to the corresponding region within the contralateral hemisphere, an analysis we referred to as IntraRegional Connectivity (IRC).

#### **MATERIALS AND METHODS**

#### **SUBJECTS**

Four NLF epilepsy individuals (two female, mean age: 37.75 – **Table 1**) with unknown pathology (MRI negative) were scanned prior to epilepsy surgery at the University of Washington (UW). Scans from NLF patients were acquired on two different scanners (three on a clinical and one on a research magnet; both Philips 3T Achieva) using identical eight-channel SENSE head coils. **Table 1** details the biographical information for each NLF subject. Note: color-coding within **Table 1** is kept consistent throughout to denote results specific to each individual NLF subject. In order to minimize variance within the NLF data sets due to use of different scanners, we downloaded functional and anatomical data sets from 16 age- and gender-matched controls (**Table A1** in Appendix) from a multisite rsfMRI repository, the 1000 connectomes database<sup>1</sup> . Of the 16 controls, one quarter were specifically matched to one NLF subject. That is, four gender-matched controls with an age range of ±1 year were selected with specific reference to each NLF subject.

#### **IMAGING**

#### **MRI acquisition**

At each UW scan session (NLF subjects), the scanning protocol included a Magnetization prepared rapid gradient echo (MPRAGE) high-resolution T1 sequence (repetition time (TR)/echo time (TE)/flip angle: 6.5 ms/3 ms/8˚; matrix size of 256 × 256 and with 170 sagittally collected slices and a slice thickness of 1 mm) and a 6-min resting state, echo planar fMRI sequence (rsfMRI, TR/TE/FA: 2000/21/90˚). The clinical scan sequence consisted of 38 axially oriented slices and a matrix size 64 × 64, while the research scan sequence consisted of 41 axially oriented slices and a matrix size 80 × 80. For all subjects, five "dummy" volumes which were collected to stabilize T1 equilibration effects were excluded from analyses. Scan parameters for the 1000 connectomes control subjects varied according to acquisition site (see **Table A1** in Appendix for details).

<sup>1</sup>http://fcon\_1000.projects.nitrc.org

#### **Seizure focus identification**

After scanning, each NLF epilepsy subject underwent a craniotomy and long-term ECoG monitoring for epileptiform discharges. Ictal onset was defined clinically from video-ECoG and identification of concordant fast spiking, low voltage activity extending from the sub-dural montage. **Figure 1** (left column) shows the ECoG montage for the four NLF subjects. Electrodes highlighted in red denote the electrodes in the ictal onset zone. After ECoG monitoring, subjects underwent surgical resection of epileptic tissue. The red transparent areas (**Figure 1**, left column) reveal the approximate location of the resected tissue as outlined by post-op surgical notes. The location of the seizure focus was defined as the region containing an overlap between ECoG recorded ictal onset activity contained within the resection zone.

#### **ANALYSIS**

#### **Pre-processing**

At the individual level, standard rsfMRI pre-processing was conducted using FEAT (FMRI Expert Analysis Tool) Version 5.98, part of FSL (FMRIB's Software Library)<sup>2</sup> to remove nonneuronal sources of variance. These included skull stripping using BET, motion correction (realignment to the center volume) with FSL MCFLIRT, spatial smoothing using a 6 mm fullwidth half-maximum (FWHM) Gaussian kernel, grand-mean intensity normalization, and linear drift removal. Identified volumes exceeding 0.5 mm of motion in any direction or plane were eliminated (scrubbed) from further processing. Additionally, ventricular CSF signal was extracted, averaged, and removed from the overall whole-brain time-series. Each 4D data set was entered into a regression analysis, treating the movement parameters and CSF signal as nuisance variables. Finally, to limit the effect of physiological noise on fc, the overall timeseries was temporally low-passed filtered removing frequencies above 0.1 Hz.

#### **Regions of interests**

Our aim was to compare across cortical regions containing the seizure focus and control regions at the individual level. Thus, we parcellated each individual subject's brain into established, known ROIs using the MNI Harvard–Oxford (HO) probability atlas (included as part of the FSL anatomical toolkit; **Figure 2A**). Each of the 48 HO cortical ROIs (employing the 25% threshold criteria) were selected, degraded by an additional 25% to prevent overlap after warping into native space and then co-registered

<sup>2</sup>www.fmrib.ox.ac.uk/fsl



into native fMRI space through a three-step registration process using FSL FLIRT. First, the native high-resolution MPRAGE was

subjects high-res T1 MPRAGE scan (generated with FREESURFER

overlaps with the electrode falling within the seizure focus for each NLF subject.

registered into native fMRI space using a rigid-body transform. Second, the MNI 2-mm standard brain was registered onto the warped MRPRAGE using an affine transformation. The generated transformation matrices from standard-to-warped MPRAGE

**FIGURE 2 | Regions of interest and local connectivity**. **(A)** Shows the 48 thresholded HO ROIs overlaid in the MNI 152 brain. Using the structural detail inherent to the high-res T1 scans, all ROIs were warped into native fMRI space for each subject. Whole brain, normalized ReHo, and IRC values were then extracted and averaged from each ROI. **(B)** Reveals an example ReHo map form one epilepsy subject. Note the cross-hairs pinpoint a qualitative decrease in ReHo in and around the seizure focus within the right hemisphere, an effect that is absent from the left. **(C)** Plots raw normalized ReHo values across the 48 HO ROIs for the same NLF subject (green bars) and the mean of the four age- and gender-matched control subjects (white bars). Epilepsy and control values are sorted from lowest to highest for the NLF subject. In this NLF subject, the ROI that contains the epileptogenic zone (ROI 38) has one of the five lowest mean normalized ReHo values of all ROIs.

were then applied to all HO ROIs. Finally, for each patient the HO co-registered ROI that contained the electrode overlaying the seizure focus was identified and selected for statistical analysis (**Figure 2**, right-most column, cross-hairs, the ROI corresponding to the seizure focus is listed in the bottom right hand corner of the black box).

#### **FC ANALYSES**

#### **Regional homogeneity**

Each 4D pre-processed data set was then passed through ReHo analysis using the REST toolbox in MATLAB<sup>3</sup> . For each voxel, a mean correlation coefficient was computed using Kendal's coefficient of concordance (KCC-ReHo), relative to the time-series from the surrounding 27 voxel neighbors. Voxel-wise ReHo values were normalized by dividing by the global mean KCC-ReHo value (Mankinen et al., 2011). Greater ReHo values denote increased local connectivity (**Figure 2B**).

#### **IRC**

To specifically contrast fc between the HO ROI contralateral to the seizure focus, we adapted traditional fc methods by computing pair-wise correlation coefficients between all possible voxel pairs within each HO ROI. Coefficients within an ROI were then transformed into *z*-scores and a mean value of absolute *z*-scores was estimated. This score was then transformed back into an average correlation coefficient yielding a mean value of intra-nodal or local fc. Finally, a ratio of local fc in the left hemisphere ROI to the connectivity in the right hemisphere ROI was calculated. If the ratio is close to 1, the brain's fc is more symmetric, vice versa. The ratios were subsequently converted into log scale resulting in degrees of asymmetry (i.e., the larger the value, the more left ROI is locally connected in comparison to the right ROI).

#### **Statistics**

We are specifically interested in whether local fc within the HO ROI containing the seizure focus is lower relative to the same ROI in controls and non-seizure focus ROIs within each epilepsy subject (thus serving as his/her own control). Because of the small patient population presenting with refractory, non-lesion epilepsy combined with interest in comparing across different fc analyses, we used a non-parametric ranking metric to evaluate differences at the group level. For each subject, HO ROIs were ranked from lowest to highest with respect to the normalized ReHo values ipsilateral to the hemisphere housing the seizure focus (for an example ranking see **Figure 2C**). IRC ROIs were sorted according to the degree of left-to-right (or right-to-left depending on which hemisphere housed the focus) asymmetry. The two rankings were then averaged. Thus, stemming from our local connectivity analysis approaches, we generated three sets of rankings of 48 values of local connectivity for each subject. The rank value of each focus ROI for each of the three rankings were entered into an independent sample Wilcoxon Rank Sum test (two-sided, alpha level of 0.05), contrasting the rank value of that ROI for the four NLF against the 16 matched controls.

Additionally,we reasoned the translational value of rsfMRI fc as a pre-operative evaluation tool would come at the individual level, contrasting local fc values across brain regions for a given surgical candidate. To characterize the ranking valuesfor each NLF epilepsy subject, we took a parametric approach calculating the mean and standard deviation of the ranks from across all controls for each of the four seizure focus ROIs. For each of these four distributions, a

corresponding *z*-score and *p* value was estimated testing the null hypothesis that the local fc rank for a given NLF epilepsy subject was no different than the controls rank values.

Finally, the mean ReHo values from each ROI was standardized to a −1 to 1 distribution in order to average across the quantitative estimate of ReHo with the left-right IRC ratios (**Table 2**). For each HO ROI, a mean standardized ReHo and IRC ratios were averaged, ranked, and compared to the mean rank values.

#### **RESULTS**

The HO ROI containing the seizure focus for each epilepsy subject ordered in the bottom 10% for all within-subject rankings except for the IRC ranking for NLF4 (red text in **Table 2**). For example, the ReHo ranking for participant NLF1 was 2 indicating the HO ROI housing the seizure focus had the second lowest mean, normalized ReHo with respect to all ipsilateral ROIs. Further, the IRC ranking for this subject was 3, indicating that this ROI showed the third lowest local fc ranking when mean local fc was directly contrasted with its contralateral counterpart. The one exception was the IRC ranking for subject NLF4, indicated that the local fc showed a greater degree of contralateral connectivity relative to the seizure focus. **Figure 3** plots the rank value for the three ranking distributions revealing the raw values for each non-lesional, focal epilepsy patient as the colored bar.

#### **GROUP-LEVEL CONTRASTS**

To determine whether local fc in the seizure focus ROI was lower in the epilepsy group, we compared the rank value of the seizure focus ROI between NLF and controls across our three sets of rankings (ReHo, IRC, and mean rank). Both the ReHo (*p* = 0.0156, Wilcoxon Rank Sum test) and the mean rank (*p* = 0.0421, Wilcoxon Rank Sum test) were significantly lower averaged across the NLF subjects (**Figure 3**, color bars) relative to controls (**Figure 3**, mean value shown in gray bars) but not the IRC fc method (*p* = 0.0184). It should be noted that the unusual contralateral connectivity effect seen with in NLF4 subject likely contributed to the null statistical effect for the IRC method at the group level.

#### **INDIVIDUAL-LEVEL CONTRASTS**

To piece out ranking effects at the individual level, we calculated *z*-score statistics from the mean and SD across ranks values from the controls. For each seizure focus ROI across each of the three local fc rankings, we were able to reject the null hypothesis for only NLF1 subject (*p* = 0.0424) under the IRC rankings. Further, when the ReHo and IRC rankings were averaged together, both subjects NLF1 (*p* = 0.0409) and NLF3 (*p* = 0.0427) showed significantly lower rankings relative to controls.

We also directly contrasted the mean rankings (i.e., the average between ReHo and IRC) for each individual NLF subject with a mean value of the raw local fc estimations. For each ROI, quantitative local fc values were an average metric calculated from the normalized ReHo and the IRC ratio scores. The mean local fc value paralleled the average ranking for all four NLF epilepsy subjects. The red text items in **Table 2** reveal the ranking and raw local fc values for each of the seizure focus ROIs. As can be seen, across both the mean rankings and combined local fc estimates, the ROI

<sup>3</sup>http://www.restfmri.net/forum/index.php

#### **Table 2 | ReHo, IRC and mean ranking values for all HO ROIs.**


(Continued)

#### **Table 2 | Continued**


Raw values for local fc measurements with the corresponding rankings are shown across all ROIs. Red, bold text represents the seizure focus ROI for each subject.

housing the seizure focus was lower in value relative to either of the constituent values alone for three out of the four NLF patients. For example, in the patient NLF 2, the seizure focus ROI was the second lowest ranked and second lowest combined computed local fc value, but ranked third and fifth when individually sorting IRC and ReHo values, respectively.

### **DISCUSSION**

Here we observe that rsfMRI local fc shows some potential as a pre-operative mapping tool for seizure focus identification in individuals with NLF epilepsy. We examined two different methods of fc estimation, averaged across both methods and contrasted local fc at the site of the seizure focus between epilepsy individuals, normal controls, and within-subject ROIs. At a group level, we observed a decrease in both the ReHo ranking and the combined rankfor the ROI harboring thefocus compared to a matched group of control subjects. This suggests that in our cohort of local epileptics there was a marked decrease in one measure of local fc (ReHo) in the area around the seizure focus. Thus, at the group level, the disease process associated with epilepsy appears to alter local fc around the focus, a hypothesis that is consistent with the pathological effects typically seen in the epileptogenic zone (Thom,2004; Wong and Guo, in press). The real clinical value however of local fc to epilepsy surgery is an accurate estimation of the location of the epileptogenic zone at the individual level that is concordant with other modalities of investigation. This is particularly important for patients with focal epilepsy with a normal MRI where macroscopic structural abnormalities are not available as an initial guide for surgical planning.

#### **A WITHIN-SUBJECT METHOD FOR IDENTIFYING THE SEIZURE FOCUS**

We therefore examined whether rsfMRI could clarify the location of the seizure focus in NLF epilepsy at an individual level. Based on neuropathological reports and animal studies of focal epilepsy reporting significant neuronal necrosis at the seizure focus, we hypothesized that values of local connectivity would be abnormal in and around the seizure focus (Thom, 2004; Wong and Guo, in press). Based upon an extensive imaging literature showing compensatory effects within the contralateral hemisphere,we extended this hypothesis to a specific decrease in local fc within the ipsilateral relative to the contralateral cortical region (Won et al., 1999; Morgan et al., 2012). The approach providing the greatest potential for revealing our predicted effects was averaging across both ReHo and IRC and contrasting across all ROIs from a single subject (**Table 2**). This procedure revealed that for three out of the four NLF subjects the seizure focus ROI was either the lowest (NLF3) or second lowest (NLF1 and NLF2) ranked ROI (see **Figure 3**, colored bars, **Table 2**; for a specific discussion on the IRC ranking of NLF4, see below). That is, the predictive capacity of local fc rsfMRI in focal, non-lesion epilepsy is improved when combining a method that specifically computes local fc within the region around the seizure focus (ReHo) with an analysis that contrasts local fc with specific reference to the corresponding contralateral hemisphere (IRC).

We argue using the epilepsy patient as his or her own control while combining these two local fc approaches provides the most promise as a translational tool. rsfMRI provides wholebrain coverage. Thus, contrasting ReHo and IRC values across the brain is readily available when employing standard clinical rsfMRI sequences. Normative, population values for patterns of local fc have not been established, more importantly are not readily available in the clinical setting and will likely need to be developed for specific MR systems and imaging sequences.When combined with established physiological and anatomical and functional imaging abnormalities commonly associated with the seizure focus as well as including the possible compensatory effects observed in the contralateral hemisphere (Bettus et al., 2009), it is not surprising that factoring in both of these methodologies would improve the overall ability to identify the epileptogenic focus ROI.

This combined approach may also have an important role in patients early in the course of the epileptogenic process where potentially surgically remediable lesions can be identified at an incipient stage before neuroanatomical changes are observed on conventional MRI. Several studies have shown that surgical interventions early in the course of pharmacoresistant epilepsy leads to better quality of life and outcomes (Engel et al., 2012).

#### **THE LARGER SEIZURE NETWORK**

For NLF1, the HO ROI 34 (corresponding to anterior division of the parahippocampal gyrus) ranked lower after combining both ReHo and the IRC methods than the seizure focus ROI (corresponding to the anterior extent of the temporal fusiform – **Figure 1**). Portions of HO ROI 34 were resected in this patient. Thus under established criteria, the parahippocampal gyrus would be included as part of the epileptogenic zone (c.f. Laufs,2012). This region shares significant inter-connectivity with cortex throughout the medial temporal lobe, including the perirhinal and entorhinal cortices as well as with the hippocampus proper (Burwell, 2000). Accordingly, the parahippocampal cortex is heavily involved in recall and/or numerous memory-related processes (Eichenbaum et al., 2007). Intrinsic connectivity studies using rsfMRI have revealed significant fc with numerous neocortical association cortices including the posterior regions of the default mode network as well as inter-connectivity spread throughout the lateral temporal lobe (Ranganath and Ritchey, 2012) and extensively with the anterior extent of the inferior temporal lobe (Kahn et al., 2008). Not surprisingly, the parahippocampal gyrus is a key fixture in the larger network underlying MTL epilepsy and seizure propagation (McIntyre and Gilby, 2008). The widespread pattern of connectivity extending from the parahippocampal region throughout the temporal lobe provides an architecture that would easily promote temporal lobe seizure propagation. With specific reference to NLF1, the seizure focus is located in a densely connected adjacent portion of the anterior, inferior temporal lobe (the temporal lobe fusiform). Thus, the observation that these two regions show the lowest local fc estimates likely signifies that rsfMRI is revealing a broader epileptogenic zone or epilepsy network in this subject.

For NLF epilepsy subject 2, only the insular cortex ROI ranked lower in local fc relative to the seizure focus ROI (located within the posterior temporal fusiform). The insula is generally considered a multimodal integration site that shares a high level of connectivity with frontal and temporal cortex. A recent seed-based rsfMRI report noted significant fc between two different points along the anterior-posterior insular plane and the posterior fusiform (Taylor et al., 2009). Both ictal and IEDs originating from the insula have been reported in MTL epilepsy (Isnard et al., 2000). In this same report, it was observed that two patients with significant insular discharges continued to have seizures after temporal lobectomy. Moreover, lesions in the insula have been shown to develop into intractable epilepsy where resection of the lesion and the surrounding insular tissue yields seizure freedom (Roper et al., 1993). Based on these and similar reports, insula-based epilepsy has become more routinely recognized over the past few decades (Nguyen et al., 2009).

The converging notion from the current NLF epilepsy patients one through three is that alterations in local fc may identify the epileptogenic zone as well as the larger epilepsy network (Stufflebeam et al., 2011). The concept of widespread epilepsy networks has been identified using both imaging with MRS (Pan et al., 2012), SPECT (Sequeira et al., 2013), FDG-PET (Mauguière and Ryvlin, 2004), and electrophysiological studies (Muldoon et al., 2013). The observation that ROIs ranking lower in local fc relative to the seizure focus likely share rich patterns of connectivity with the seizure focus may be exposing a more widespread pathological consequence of the seizure propagation. Building upon the hypothesis that discrepancies in local fc are linked to local neuronal insults such as necrosis (or apoptosis), alterations in dendritic morphology, and potential compensation within the contralateral hemisphere, the currently applied techniques may be revealing the downstream consequences of seizure propagation across the entire epilepsy network.

#### **METHODOLOGICAL CONSIDERATIONS AND LIMITATIONS**

We choose to focus specifically on refractory, non-lesion epilepsy patients because of the added importance that functional-based modalities (i.e., electrophysiological and imaging based procedures) provide in the pre-surgical localization of the seizure focus. The number of patients presenting with NLF epilepsy that are candidates for surgery are however relatively small (<10% of all new cases per year; Duncan, 2010). Despite this limitation, the current results should be taken with a degree of caution due the small sample size. As a follow-up, future studies will clearly need to conduct similar analyses with larger samples. It is however likely that estimates of local fc may aid in the identification of the epileptogenic focus among patients presenting with various focal pathologies (i.e., cortical dysplasia, AVM, brain tumors etc.). Taken together with the lateralized fc differences throughout the medial temporal lobe previously reported in MTLE patients (Bettus et al., 2009), local fc would likely contribute to the pre-surgical evaluation even in the presence of an identified insult.

The current results would benefit from a more precise delineation of the epileptogenic zone. Other groups have identified the epileptogenic zone using a variety of additional techniques (c.f. Jayakar et al., 2008; Duncan, 2010). We were not able to use a more sophisticated means of defining the epileptogenic zone other than a description from post-op surgical notes of the extent and boundaries of the resected region. By choosing to parcellate the brain into ROIs using a well-established, probabilistic atlas combined with a sorting method based on mean local fc values, we ensured a completely unbiased process of identifying patterns of reduced local fc across subjects while maintaining relatively high anatomical specificity. One unfortunate and likely consequence of this procedure is a smearing of voxel types within an ROI. More specifically, it is unlikely that the ROI corresponding to the seizure focus in any given NLF patient contains voxels that would be exclusively labeled as falling in or exclusively out of the epileptogenic zone. Thus, it is likely that the mean values for each ROI in and around the epileptogenic zone are underestimated, and the true local fc value associated with the epileptogenic zone is likely lower. One possible solution for consideration in future studies is to contrast pre and post-resection MRI scans. This would generate a voxel mask of the resected tissue and by extension the extended epileptogenic zone. Furthermore, the current results would indeed benefit from the addition of simultaneously acquired EEG. Confirmation of the IED-related activity during rsfMRI acquisition would provide the ability to confirm the boundaries of epileptogenic zone. Provided the presence of IEDs during functional scanning, it may be feasible to select out specific periods of "IED-free" rsfMRI activity in order to determine whether the presence of IEDs are negatively (or positively) impacting local lc correlation coefficients. However, we reason that rsfMRI provides a simple yet powerful means of examining the underlying physiology of the epileptogenic zone that is also feasible in the clinical context (c.f. Fox and Greicius, 2010). Future studies will clearly need to address the influence of IEDs (as well as ictal discharges) on the rsfMRI BOLD activity and local fc estimates. Furthermore, future studies will need to address the concordance between rsfMRI local fc estimates in NLF epilepsy and more commonly used modalities such as FDG-PET. However, if local fc does indeed reflect the accurate location of the seizure focus and thereby supplementing more traditional evaluative modalities, then the need of simultaneous EEG would prove relatively superfluous.

NLF 4 did this not show the same pattern of IRC within the seizure focus ROI (located within the left middle temporal gyrus) as was observed in other NLF 3 patients. Although the raw and ranked ReHo values were within the bottom 10% of all sorted ROIs, the pattern of local fc under the IRC calculation was significantly greater within the ipsilateral hemisphere. The mechanism contributing to this effect is unknown. Results from the WADA test as well as clinical fMRI scans using various language screens concluded that language dominance was localized to the left hemisphere for this patient. It is possible that patterns of contralateral connectivity are not as vast within the middle, temporal lobe relative to noted contralateral compensatory effects stemming from medial temporal lobe (Bettus et al., 2009). It is also conceivable that scalars of local fc are greater in regions throughout the language dominant hemisphere relative to the contralateral counterparts. It is clear that future work will need to address baseline differences in local fc across both the temporal lobe as well as whole brain.

#### **CONCLUSION**

We present evidence suggesting local fc measurements from rsfMRI provide an accurate estimate of the location of the epileptogenic region in non-lesional, focal epilepsy. Structurally identified lesions are typically considered a reliable guide as a first pass for identifying the approximate location of the epileptogenic zone. Because the long-term benefits of epilepsy surgery are significant for individuals presenting with normal anatomical MRIs (Jayakar et al., 2008), accurate localization is a critical pre-operative function. In the absence of identified lesions, clinicians must rely more heavily on alterative methods to identify epileptogenic zones. Here we provide the first evidence that rsfMRI local fc may provide additional, confirmatory information about the location of the epileptogenic focus in refractory NLF epilepsy. These techniques may also identify the broader epilepsy network and identify comorbid neuropsychological dysfunction due to involvement of other functional networks.

#### **REFERENCES**


S., et al. (2012). Early surgical therapy for drug-resistant temporal lobe epilepsy: a randomized trial. *JAMA* 307, 922–930.


#### **ACKNOWLEDGMENTS**

This research was supported by NIH/NIMH grant 5K01 MH086118-03 (K. E. Weaver) and NIH/NINDS 5ROI NS065186- 03 & The Dreuding foundation (J. G. Ojemann).

lobe epilepsy. *Brain Res.* 1373, 221–229.


T., Joseph, J. E., et al. (2013). Perfusion network shift during seizures in medial temporal lobe epilepsy. *PLoS ONE* 8:e53204. doi:10.1371/journal.pone.0053204


tonic–clonic seizures. *Epilepsy Res.* 97, 83–91.

Zijlmans, M., Huiskamp, G., Hersevoort, M., Seppenwoolde, J.-H., van Huffelen, A. C., and Leijten, F. S. S. (2007). EEG-fMRI in the preoperative work-up for epilepsy surgery. *Brain* 130, 2343–2353.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 20 February 2013; paper pending published: 14 March 2013; accepted: 17 April 2013; published online: 01 May 2013.*

*Citation: Weaver KE, Chaovalitwongse WA, Novotny EJ, Poliakov A, Grabowski TG and Ojemann JG (2013) Local functional connectivity as a pre-surgical tool for seizure focus identification in nonlesion, focal epilepsy. Front. Neurol. 4:43. doi: 10.3389/fneur.2013.00043*

*This article was submitted to Frontiers in Epilepsy, a specialty of Frontiers in Neurology.*

*Copyright © 2013 Weaver, Chaovalitwongse, Novotny, Poliakov, Grabowski and Ojemann. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*

#### **APPENDIX**

**Table A1 | Control subject demographic information and scanning parameters.**


## Computer-aided diagnosis and localization of lateralized temporal lobe epilepsy using interictal FDG-PET

#### **Wesley T. Kerr 1,2\*, Stefan T. Nguyen<sup>3</sup> , AndrewY. Cho<sup>2</sup> , Edward P. Lau<sup>2</sup> , Daniel H. Silverman<sup>3</sup> , Pamela K. Douglas <sup>2</sup> , Navya M. Reddy <sup>3</sup> , Ariana Anderson<sup>2</sup> , Jennifer Bramen<sup>2</sup> , Noriko Salamon<sup>4</sup> , John M. Stern<sup>4</sup> and Mark S. Cohen2,5**

<sup>1</sup> Department of Biomathematics, David Geffen School of Medicine, University of California Los Angeles, Los Angeles, CA, USA

<sup>2</sup> Laboratory of Integrative Neuroimaging Technology, Department of Psychiatry, Neuropsychiatric Institute, University of California Los Angeles, Los Angeles, CA, USA

<sup>3</sup> Ahmanson Translational Imaging Division, Department of Molecular and Medical Pharmacology, David Geffen School of Medicine, University of California Los Angeles, Los Angeles, CA, USA

<sup>4</sup> Department of Neurology, Seizure Disorder Center, University of California Los Angeles, Los Angeles, CA, USA

<sup>5</sup> Laboratory of Integrative Neuroimaging Technology, Departments of Psychiatry, Neurology, Radiology, Biomedical Physics, Psychology and Bioengineering, University of California Los Angeles, Los Angeles, CA, USA

#### **Edited by:**

Mark Holmes, University of Washington, USA

#### **Reviewed by:**

Don Tucker, Electrical Geodesics, Inc. and the University of Oregon, USA Mark Holmes, University of Washington, USA

#### **\*Correspondence:**

Wesley T. Kerr, Laboratory of Neuroimaging Technology, Department of Biomathematics, University of California Los Angeles, 760 Westwood Plaza, Suit B8-169, Los Angeles, CA 90095, USA. e-mail: wesleytk@ucla.edu

Interictal FDG-PET (iPET) is a core tool for localizing the epileptogenic focus, potentially before structural MRI, that does not require rare and transient epileptiform discharges or seizures on EEG. The visual interpretation of iPET is challenging and requires years of epilepsy-specific expertise. We have developed an automated computer-aided diagnostic (CAD) tool that has the potential to work both independent of and synergistically with expert analysis. Our tool operates on distributed metabolic changes across the whole brain measured by iPET to both diagnose and lateralize temporal lobe epilepsy (TLE).When diagnosing left TLE (LTLE) or right TLE (RTLE) vs. non-epileptic seizures (NES), our accuracy in reproducing the results of the gold standard long term video-EEG monitoring was 82% [95% confidence interval (CI) 69–90%] or 88% (95% CI 76–94%), respectively. The classifier that both diagnosed and lateralized the disease had overall accuracy of 76% (95% CI 66–84%), where 89% (95% CI 77–96%) of patients correctly identified with epilepsy were correctly lateralized. When identifying LTLE, our CAD tool utilized metabolic changes across the entire brain. By contrast, only temporal regions and the right frontal lobe cortex, were needed to identify RTLE accurately, a finding consistent with clinical observations and indicative of a potential pathophysiological difference between RTLE and LTLE. The goal of CADs is to complement – not replace – expert analysis. In our dataset, the accuracy of manual analysis (MA) of iPET (∼80%) was similar to CAD. The square correlation between our CAD tool and MA, however, was only 30%, indicating that our CAD tool does not recreate MA. The addition of clinical information to our CAD, however, did not substantively change performance. These results suggest that automated analysis might provide clinically valuable information to focus treatment more effectively.

**Keywords: epilepsy, computer-aided diagnosis, mutual information, temporal lobe epilepsy, PET, fluoro-deoxyglucose positron emission tomography, machine learning**

#### **INTRODUCTION**

It is difficult to differentiate between patients with epilepsy (PWE), and those with non-epileptic seizures (NES). The clinical assessment relies on the report of untrained witnesses or the patients themselves. A non-epileptic seizure is defined as the presence of external seizure symptoms and/or signs with no electrographic features characteristic of epilepsy. Long term video-EEG monitoring has shown consistently that roughly one third of patients diagnosed with "medication refractory epilepsy" in fact suffer from NES (Kerr et al., 2012a). Because they don't suffer from epilepsy, these patients with NES (PWN) are not treated effectively with anti-epileptic drugs (AEDs). For the majority of PWN, the NES are a manifestation of dissociative or conversion disorder in which their psychological challenges manifest themselves

physically (Marchetti et al., 2008, 2009). A minority of PWN suffers from organic, non-epileptic maladies that can be confused with seizure disorder including, but not limited to, dementia and cardiovascular disease (Sahaya et al., 2011). The gold standard for the differential diagnosis and pre-surgical assessment of epilepsy includes 72 or more hours of video-EEG monitoring (Cragar et al., 2002; LaFrance and Devinsky, 2004). However, 10% of patients admitted for this extensive assessment leave with inconclusive results (Kerr et al., 2012a). Considering that one sixth of PWE are diagnosed with medication refractory epilepsy (Privitera, 2011), improved methods to effectively identify PWN who do not benefit from AEDs effectively could reduce the morbidity and both the financial and social cost of treating epilepsy.

Improved diagnostic tools could also help PWE. The difficulty in ruling out non-epileptic etiologies speaks to the challenge of adequately localizing and characterizing each patient's epileptic etiology. The major seizure type discriminations are focal vs. generalized; partial vs. complex; and lesional vs. nonlesional. Each of these key discriminations leads patients down a different treatment path. When medication or other novel treatments like the vagus nerve stimulator fails, as they frequently do, the patient is left to consider resective neurosurgery. Recent reports have shown that surgery is most effective earlier in the course of disease (Engel et al., 2012). Improved diagnostic tools could more quickly and effectively diagnose patients with epileptic seizures and therefore speed the progression toward considering the surgical option.

Ultimately, our goal is to establish a general, automated computer-aided diagnostic (CAD) tool that effectively combines clinical information, manual interpretation of EEG and imaging technologies as well as automated analysis of interictal FDG-PET (iPET), EEG, structural MRI (sMRI), and diffusion MRI for all subtypes of epilepsy and NES. To accomplish this, we first must develop effective CAD tools that harness the information from each modality for a limited set of epileptic localizations. We have begun already to address automated analysis of interictal EEG for a wide variety of epilepsy subtypes (Kerr et al., 2012a). Others have described effective CAD tools that diagnose and lateralize temporal lobe epilepsy (TLE) using structural and diffusion MRI (Farid et al., 2012; Focke et al., 2012; Keihaninejad et al., 2012).

The clinical, metabolic, and structural differences between left and right TLE can be subtle. Some theories suggest that TLE is inherently a bilateral disease. Potentially, due to the strong functional link between the hippocampi, the only clinical difference is that in the aura of patient with left TLE (LTLE) more frequently includes language dysfunction. Over time, patients with LTLE more commonly develop verbal memory deficits, compared to non-verbal memory deficits in right TLE (RTLE) (Delaney et al., 1980; Kim et al., 2003). This functional connection between the hippocampi may also lead some patients to be falsely-lateralized using scalp EEG because a small seizure onset zone (SOZ) in one hippocampus can induce larger scale ictal activity in the contralateral hippocampus with very little time delay. This can lead neurologists to falsely conclude that the SOZ is either bilateral or in the contralateral hippocampus. Structural and metabolic imaging can reduce these errors by demonstrating that that one temporal lobe is asymmetrically affected, as shown by the previous CAD tools that lateralize TLE (Farid et al., 2012; Focke et al., 2012; Keihaninejad et al., 2012). Studies of the functional connectivity of these epileptic networks, however, conclude that there are very few,if any, differences between the two lateralizations (Zhang et al., 2010; Liao et al., 2011; Morgan et al., 2011, 2012; Pittau et al., 2012; McCormick et al., 2013). Recently, Pereira et al. (2010) suggested that more patterns of functional connectivity change in LTLE compared to RTLE. However, after patients suffer from intractable seizures for 10 or more years, the intrahemispheric hippocampal connectivity linearly increases with the duration of disease, suggesting that over time lateralized disease may become bilateral disease (Morgan et al., 2011). Because patients with bilateral hippocampal disease are no longer considered surgical candidates, improved methods to distinguish left and right TLE early in the course of disease are needed.

In this manuscript, we discuss the development of an automated CAD tool to diagnose, and lateralize, TLE using iPET. We also begin to address how to combine our CAD tool with manual analysis (MA) and incorporate it into clinical practice. Using a mutual information-based feature selection technique, we examine how our methods reveal more about the distributed metabolic abnormalities that are associated with the different anatomical locations of the epileptogenic focus.

The realistic goal of CAD tools is to complement, not to replace, expert analysis. Therefore, we focus on how clinical information and expert analysis can work synergistically with our automated technology. To summarize the major clinical differences, patients with NES are characteristically females in the third decade of life with psychiatric co-morbidities (Sahaya et al., 2011). PWE, however, also have significant psychiatric co-morbidities including potentially reduced financial and social independence due to the suspension of their driver's and, frequently, professional license. Particularly in adult onset epilepsy, age-associated changes in metabolism may confound the interpretation of iPET, possibly leading to an increased diagnostic uncertainty. It is well established that 80–90% of medication refractory epilepsy is "PET positive" (Salamon et al., 2008; Lee and Salamon, 2009).The rate of PET positivity in NES has not been studied extensively, therefore the true positive predictive value of iPET is unclear. Although these differences in clinical presentation are salient, their quantitative effect on diagnostic probabilities is unknown. Therefore, we also examined how simple clinical information and expert manual interpretation can be incorporated into our quantitative CAD tool.

The standard of care for the pre-surgical assessment for epilepsy is the manual correlation of iPET with numerous other diagnostic modalities. The goal of this assessment is to simultaneously verify the diagnosis of epilepsy, characterize the seizure etiology, and identify the location and extent of the SOZ. Expert radiologists and neurologists can detect metabolic asymmetries indicative of the epileptogenic focus or foci (Person et al., 2010). The exact threshold at which asymmetric metabolism is attributed to pathologic change or seen as a variant of normal is part of the art of neuroradiology (Benbadis et al., 2000; Reuber et al., 2002). Once non-epileptic etiologies have been ruled out, our previous work demonstrated that the quantitative degree of metabolic asymmetry is correlated with surgical outcome (Lin et al., 2007). Surgical outcome is improved further when iPET is co-registered to sMRI because of improved characterization of the focus or foci (Chandra et al., 2006; Rastogi et al., 2008; Salamon et al., 2008; Lee and Salamon, 2009). These hypometabolic lesions are thought to be secondary to increased inhibitory neuron cell death, gliosis, and abnormal functional connectivity resulting in altered functional metabolism.

The size of the hypometabolic lesion tends to be larger than the SOZ, potentially due to functional changes in nearby tissue secondary to the presence of the epileptogenic lesion (Juhasz et al., 1999; Matheja et al., 2001; Henry and Roman, 2011). Such reports are major limitations to the wide implementation of iPET in epilepsy practices (Barrington et al., 1998; So et al., 2000; Henry et al., 2011). In addition to the limitation of counting statistics, that forces the quantitative radioactivity intensity of iPET to be less certain in hypometabolic lesions (Kerr and Lau, 2012), the biological hypothesis is that the epileptogenic abnormality induces metabolic abnormality at the SOZ and also at closely associated and/or functionally connected regions (Henry et al., 1990, 1993; Sperling et al., 1990; Sadzot et al., 1992; Arnold et al., 1996; Dlugos et al., 1999; Bouilleret et al., 2002; Rusu et al., 2005; Nelissen et al., 2006; Takaya et al., 2006; Lee et al., 2009). The epileptogenic lesion commonly is larger and more diffuse in left TLE then right TLE, potentially because of the high degree of functional connectivity between specialized foci within the left temporal lobe associated with language and other functions (Toga and Thompson, 2003; Barrick et al., 2005; Iturria-Medina et al., 2011; Haneef et al., 2012; Kucyi et al., 2012). These insights parallel the trend in dementia that atrophy starts focally then spreads more quickly to functionally connected regions (Zhou et al., 2012). The limited sensitivity of iPET unaligned with sMRI to characterize extratemporal lesions may be partly due to the insufficient description of the local functional network of each extratemporal focus and thereby reduced detection of a characteristic pattern of metabolic abnormalities associated with each focus. In general, an improved insight into the clinical interpretation and value of metabolic abnormalities outside the SOZ is needed. To overcome this limitation, the iPET analysis is used in combination with other diagnostic modalities determine which tissue to resect.

Clinical description, EEG, MRI, and FDG-PET each describe separate facets of the pathophysiological etiology, and therefore all play critical roles in the diagnosis of epilepsy, and in the identification of the epileptogenic lesion (Struck et al., 2011). Each modality, however, also has unique limitations. EEG provides an in-depth description of the seizures and interictal epileptiform spikes. These seizures and spikes, however, are rare events: only 50% of PWE exhibit diagnostic interictal epileptiform spikes and/or seizure activity during the first outpatient scalp EEG (Gilbert et al., 2003). The characteristic signs of epilepsy in structural and diffusion MRI may not be measurable until years after the first seizure because these methods require the detection of atrophic tissue and/or subtle regions of cortical dysplasia (Swartz et al., 1992; Reutens et al., 1996; Van Paesschen et al., 1998; Liu et al., 2002; Jung da and Lee, 2010; Bernasconi et al., 2011; Schmidt and Pohlmann-Eden, 2011; Dabbs et al., 2012). MA uses the contralateral structure to assess if atrophy is present but a certain degree of asymmetry is expected (Farid et al., 2012; Keihaninejad et al., 2012). It takes years of specific experience in manually analyzing sMRIs from PWE to reliably discriminate between normal variation and pathologic changes. Once these relatively large-scale changes in neural structure have occurred, it is less likely that both invasive and noninvasive treatments will be effective (Engel et al., 2012). iPET can localize the epileptogenic lesion without observing rare events and, potentially, before changes are detectible on sMRI and/or diffusion tensor imaging (DTI) (Theodore et al., 1990; Ryvlin et al., 1991; Swartz et al., 1992; Gaillard et al., 1995; Debets et al., 1997; Knowlton et al., 1997, 2008; Blum et al., 1998; Drzezga et al., 1999; Benedek et al., 2004; Carne et al., 2004; Chandra et al., 2006; Yun et al., 2006; Uijl et al., 2007; Willmann et al., 2007; Rastogi et al., 2008; Salamon et al., 2008; Duncan, 2009; Lee and Salamon, 2009; Lerner et al., 2009; Liew et al., 2009; Brodbeck et al., 2010; Chinchure et al., 2010; Kim et al., 2011; Jupp et al., 2012). As discussed above, the presence of metabolic abnormalities outside the SOZ, however, complicates the effective localization of the SOZ using iPET alone (Henry and Roman, 2011). An improved description of these induced changes outside the SOZ may help spare healthy tissue from resective surgery. Given the recent report that resective neurosurgery for epilepsy is more effective earlier in disease (Engel et al., 2012); we believe that iPET may play a critical role in characterizing patients with unremarkable MRIs and inconclusive EEGs earlier in the course of their disease.

### **MATERIALS AND METHODS**

#### **PATIENT DATA**

All of the 105 patients that were included in our analysis were admitted to the University of California, Los Angeles (UCLA) Seizure Disorder Center's video-EEG Epilepsy Monitoring Unit (EMU) between 2005 and 2012. Each patient's diagnosis was based on a consensus panel review of their clinical history, physical and neurological exam, neuropsychiatric testing, video-EEG, iPET, ictal FDG-PET, structural and diffusion MRI, and/or CT scan. This multimodal assessment is the gold standard for epilepsy diagnosis and localization of the epileptic focus (Cragar et al., 2002; LaFrance and Devinsky, 2004). The patients included in this analysis were chosen because they had an FDG-PET after 2005; had no history of penetrative neurotrauma,including neurosurgery; were determined by consensus diagnosis to have a single, lateralized epileptogenic focus; and had no suspicion of mixed non-epileptic and epileptic seizure disorder. These patients were diagnosed either with LTLE (*n* = 39), right TLE (RTLE, *n* = 34), or NES (NES, *n* = 32). PET images were determined to be interictal by clinical findings and concurrent scalp EEG.

PET and MRI images were acquired according to the best clinical practices at the time of acquisition. PET/CT studies were acquired using a Siemens Biograph scanner. After a minimum fasting period of 6 h, patients received 0.14 mCi/kg of 18F-FDG-PET intravenously. During the ensuing 40 min uptake period with concomitant EEG monitoring to confirm interictal status, the patients waited in a quiet, dimply lit room with their eyes open. PET images were reconstructed with an iterative algorithm (OSEM: 2 iterations, 8 subsets). CT images were reconstructed using filtered back projection at 3.4 mm axial intervals to match the slice separation of the PET data, and used for attenuation correction.

#### **COMPUTER-AIDED DIAGNOSTIC TOOL TRAINING AND VALIDATION**

Automated analysis of the iPET records was performed in four stages. (1) First, each image was screened for gross structural and/or metabolic abnormalities by S.T.N., N.M.R., and/or W.T.K. (*n* = 21). These excluded subjects are not reflected in the sample sizes quoted above. (2) NeuroQ (Syntermed, GA, USA) was used to segment each brain into 47 regions of interest (ROIs) and then to calculate the average radioactivity in each ROI, normalized by the whole brain radioactivity (**Table A1** in Appendix). (3) The minimum redundancy-maximum relevancy (mRMR) toolbox for MATLAB (Mathworks, MA, USA) was used to generate a ranked list of the ROI metabolisms (features) within each training set that were maximally relevant to the diagnosis of epilepsy and minimally redundant with all higher ranked features (Ding and Peng, 2005; Peng et al., 2005). The representative number of features to exclude and quantal levels was selected based on our method discussed previously (Kerr et al., 2012a,b) (see below). In each of the training sets, the feature ranking was determined exclusive of the test patient's data. We expect the ranked lists to be similar, but not identical, across training sets. For purely illustrative purposes, the full dataset was used to create the ranked list in **Table 2**. (4) Weka was used to implement leave-one-out cross-validation of a cost-sensitive Multilayer Perceptron (MLP) that was weighted to maximize balanced accuracy, defined by the mean of sensitivity and specificity (Bouckaert et al., 2010). Using this method, we examined our ability to diagnose either LTLE or RTLE from NES and assessed our ability to diagnose and lateralize disease simultaneously. For the remainder of this manuscript, the latter tool that discriminates LTLE vs. RTLE vs. NES is called the trinary classifier. Similarly, the binary CAD tools are referred to by the laterality of epilepsy that is being detected. The comparison to NES is not stated, but can be assumed. We then compared our CAD tool's performance to the results of MA alone.

#### **MACHINE LEARNING ALGORITHMIC DETAILS**

The MLP was implemented with default parameters in Weka (Bouckaert et al., 2010). All input features were normalized to values between negative and positive 1. No limit was set on the number of hidden layers or nodes within each hidden layer. These parameters were optimized within each training set independently. The learning rate and momentum were set to 0.3 and 0.2, respectively. Five hundred epochs were used for training. During training, models with more than 20 consecutive errors were excluded. The trinary classifier was created by decomposing the three class problem into three 1-against-1 problems that were combined using majority voting. No three-way ties occurred during training or testing.

Balanced accuracy was optimized using a cost-sensitive classifier in which a false positive was given a cost of *n*+ and a false negative was given the cost of n−, where n+ and n− represent the number of PWE and NES in the full sample, respectively. In the trinary classifier, the cost was set as the sum of the number of patients in the other two diagnostic classes.

Cyclical leave-one-out cross-validation (CL1OCV) was used to assess the performance of the MLP. In this paradigm, all but one patient was used to determine the features selected and train the algorithm. The single remaining patient is tested using the model built upon the other patients. The identity of the test patient is permuted until all patients have been the test case once and only once.

To determine the number and identity of the input features, the mRMR algorithm requires the number of input features, *F*, and quantal levels,*Q*, be set *a priori.* For the calculation of mutual information, the features were smoothed into *Q* quantal bins akin to the bins in a histogram. Classification, however, utilizes unsmoothed features. The choice of input features smoothed into quantal levels was determined to be most representative of the performance of the algorithm across a wide variety of choices of *F* and *Q* (Kerr et al., 2012b). This choice was made by selecting a point within a region of *F-Q* parameter space that performed significantly better than the naïve classifier with 95% confidence based on random field theory correction where the spatial smoothness is estimated directly from the data (for more details, see Worsley et al., 1992, 2004; Chauvin et al., 2005). The naïve classifier classifies all test exemplars as the most common class in the training set. Under the CL1OCV procedure, these input features were determined independently for each of the training samples. The illustrated rank order of features was calculated based on the full dataset, and does not necessarily match the rank list of any individual training sample.

When clinical information was incorporated into the algorithm, the same methodology was applied as above, except that all exemplars with missing data were excluded from analysis. In these additional analyses, we did not re-sample the parameter space of *F* and *Q.* We simply used the selections determined in the previous analysis.

#### **MANUAL ANALYSIS OF PET AND MRI RECORDS**

Manual analyses of the iPET and sMRI records were performed based on the review of clinical records primarily written by Dr. Noriko Salamon. Dr. Salamon has 10 years of experience in the pre-surgical assessment of epilepsy using FDG-PET and MRI. All manual interpretation was conducted for the clinical assessment of each patient when it occurred, prior to the CAD tool development. Therefore, Dr. Salamon was blinded to the automated results. Due to the unclear relationship between structural and metabolic abnormalities, asymmetries, and epilepsy, all abnormal results were interpreted to be consistent with some form of epilepsy. Not all patients had sMRI (*n* = 6) and iPET (*n* = 1) reports available; therefore all analysis regarding MA of neuroimaging includes only patients with available records. These patients had raw iPET data available; they therefore were included in the automated analysis.

#### **COMBINATION OF CLINICAL INFORMATION WITH COMPUTER-AIDED DIAGNOSTIC INFORMATION**

To examine the combined power of clinical knowledge, MA and our automated analysis, we assessed the linear correlation of detecting epilepsy with CAD compared to MA, and also incorporated clinical information and MA into our algorithm in two ways. First, the clinical literature suggests that patients with NES are more likely to be female, begin having seizures in the third decade of life, have a decreased duration of disease and have increased seizure frequency (**Table 1**). Although we did not see a significant difference seizure frequency within our dataset, we included this features to better match clinical practice. These clinical features were then added to the input and leave-one-out cross-validation was repeated. Secondly, to explore how our computational methods can complement clinical wisdom, we included the results of MA of the iPET and sMRI as two additional input features and re-evaluated CAD performance. For the trinary classifier only, we split each of the features describing the iPET and sMRI MA to indicate if a left and/or right sided abnormality was reported.

To assess the applicability of our CAD as a *separate* modality that could be considered as part of the clinical assessment of epilepsy, we calculated the likelihood ratios (LRs) of each of the combinations of our CAD with MA of iPET and/or sMRI.


#### **Table 1 | Clinical information and results of manual analysis.**

This table reflects the clinical information known before the application of the CAD tool. All times are listed in years (y) unless otherwise specified by days (d), weeks (w), or months (m). Manual analysis of all patients' iPET and sMRI were not done, therefore we list the number with available manual results. \*, § , or ¶ indicate that the value for NES vs. LTLE, NES vs. RTLE, or LTLE vs. RTLE, respectively, is statistically significant from both the LTLE and RTLE groups with at least 95% confidence using a two-sample z-test of proportions or Mann–Whitney U test, where appropriate. No other differences are statistically significant (p > 0.10).

This was done only for the binary classifiers, because LRs have a clear formulation only for binary outcomes. The likelihood ratio is defined by the likelihood that a patient with a certain combination of diagnostic outcomes has epilepsy, divided by the likelihood that the same patient has NES. Intuitively, a likelihood ratio of two implies that the patient is twice as likely to have epilepsy. The 95% confidence intervals of chance were calculated using exact binomial intervals by considering the likelihood ratio of a classifier that diagnosed patients according to their prior likelihood alone, conditioned upon the assumption that the same total number of patients would have the diagnostic outcome of interest. For example, 39 of 71 patients had LTLE when we discriminated between LTLE and NES, therefore the median LR is 1.2. Thirty-five patients from the NES vs. LTLE group had negative MA of their iPET. Therefore, we use a binomial distribution with 35 trials and success probability of 39 over 71 to yield a 95% confidence interval of 0.94–3.38.

#### **RESULTS**

All of our results are compared to the gold standard diagnosis from the consensus panel. The clinical trial statistics of each of our automated diagnostic tool matched, but were not redundant with, expert MA of both interictal PET and sMRI (**Figure 1**). All intervals reflect 95% confidence intervals and all *p*-values correspond to differences from anaïve classifier. The binary CAD tool for RTLE had accuracy of 88% (69–90%), compared to the accuracy of MA of iPET [85% (72–92%)] and sMRI [77% (63–85)]. The binary tool for LTLE had accuracy of 83% (69–90%), compared to the accuracy of MA of iPET [79% (66–88%)] and sMRI [70% (56–81%)]. The pattern in sensitivities, specificities, and odds ratios all parallel this trend where our automated diagnostic tools are non-statistically superior to MA oriPET, which, in turn, are non-statistically superior to MA of sMRI (**Figure 1**). The accuracy

of our trinary CAD tool that simultaneously diagnoses epilepsy and lateralize disease was 76% (66–84%), where 89% (77–96%) of patients correctly identified with epilepsy were also lateralized correctly. MA to diagnose and lateralize was 78% (69–86%) accurate with 89% (76–94%) correctly lateralized using iPET and 71% (61–80%) accurate with 91% (78–97%) correctly lateralized using sMRI.

The rank order of the features used in our algorithm parallel the clinical observation that the epileptogenic networks in LTLE are broader than in RTLE. The LTLE vs. NES classifier achieved its performance by utilizing trends across almost the entire brain by including 42 of the 47 features in the final algorithm. In contrast, the RTLE vs. NES classifier only needed to measure the metabolism in six regions – bilateral temporal cortex and two associated regions of cortex – to achieve its impressive performance (**Table 2**). As expected, the trinary classifier utilized an intermediate number of features to achieve its accuracy (30 of 47). The rank list of these features matches the biological intuition based on knowledge about the potential connectivity of epileptogenic networks (**Table 2**).

We then considered how this CAD information could be used in combination with clinical information or expert analysis. The squared correlation of our CAD tool with manually interpreted iPET was 0.25 (0.09–0.43), 0.32 (0.17–0.54), and 0.34 (0.17–0.46) for the LTLE, RTLE, and trinary classifiers, respectively (**Figure 2**). The squared correlation of our tool with manually interpreted sMRI was 0.07 (0.001–0.23),0.21 (0.06–0.40), and 0.11 (0.02–0.25) for the LTLE, RTLE, and trinary classifiers respectively. For comparison, the squared correlation between manually interpreted iPET and sMRI was 0.17 (0.06–0.33).

When the same automated analysis was used to combine clinical findings with our iPET data, performance did not change significantly. After the four clinical factors were added to the input

figures indicate the accuracy, sensitivity and specificity of the LTLE **(A)**, RTLE **(B)** and trinary **(C)** classifiers. The performance of our CAD tools matched that of MA and was superior to just using gender alone. The error bars indicate standard error of the mean performance for each measure. The translucent region indicates the performance of a naïve classifier. \*Indicates significant differences from the naïve classifier with a confidence level of 95% or more.

of our tools, the accuracy changed to 79% (66–88%), 68% (56– 79%), and 64% (54–73%) for the LTLE, RTLE, and trinary classifiers, respectively (**Figure 3**). These accuracies do not substantively


This table illustrates the top six informative and non-redundant regions of interest (ROIs) that may contribute to each of the CAD tools, as determined by the minimum redundancy-maximum relevancy criteria (mRMR; Ding and Peng, 2005; Peng et al., 2005). The illustrated rank order of features was calculated based on the full dataset and does not necessarily match the rank list of any individual training sample. The leading L or R indicates left or right. The lowercase letters indicate inferior (i), lateral (l), median (m), anterior (a), and posterior (p).The lagging C signifies cortex. Note that the LTLE vs. NES and trinary classifiers include information from 42 and 30 ROIs, respectively. To better understand the benefit of mRMR, this list can be directly compared to the list of ROIs ranked by t-statistics in**Table A1** in Appendix.

some information is shared, the majority of information provided by our CAD tools is not captured by MA. The correlation between MA of iPET and sMRI is similar in magnitude to the correlation of CAD with MA, therefore the CAD could potentially be seen as similar to another informative modality. \*Indicates significant differences of the correlation from zero with a confidence level of 95% or more.

change when only sex and duration of disease were considered (results not shown). Adding the results of MA of both iPET and sMRI to our iPET data changed the accuracy to 82% (73–91%), 77% (67–88%), and 68% (59–77%) for the LTLE, RTLE, and trinary classifiers, respectively. When all information sources contribute to the algorithm, the accuracy changed to 77% (68–88%),

**FIGURE 3 | Automated combination of clinical information with automated analysis of iPET images.** The automated combination of clinical information and/or MA with our analysis produced no significant change in performance for the LTLE **(A)**, RTLE **(B)** or trinary **(C)** classifiers, relative to the CAD operating on automated values alone. The unshaded bars indicate the performance of similarly constructed CAD tools using clinical information or the results of MA alone. The shaded bars indicate the modified performance when information from NeuroQ is added. The horizontal line indicates the mean accuracy of each CAD tool without clinical information. The translucent region indicates the performance of a naïve classifier.

74% (64–85%), and 76% (68–84%) for the LTLE, RTLE, and trinary classifiers, respectively.

We combined the results of MA were combined with our CAD tool manually using LRs. After doing so, the likelihood was generally only significant if all considered modalities agreed. Viewed alone, MA and our CAD increased the likelihood of the predicted outcome between two and ninefold (*p* < 0.02; **Figure 4A**). When two analysis streams were combined, if both analyses agreed, the likelihood of the predicted outcome was increased between 8- and 27-fold (*p* < 3 × 10−<sup>4</sup> ; **Figures 4B,C**). If all three analyses agreed, the likelihood of the predicted outcome increased more than 15 fold (*p* < 1.3 × 10−<sup>5</sup> ; **Figure 4D**). However, in most cases, if there was any disagreement, the likelihood did not change significantly, most probably due to the small numbers of patients with each potential outcome. There are two key exceptions: (1) Given iPET results indicating NES over RTLE using either MA or CAD, the sMRI could be largely ignored (*p* < 1.1 × 10−<sup>2</sup> ). (2) If both MA and CAD of iPET agreed that a patient suffered from LTLE and not NES, the sMRI results could be similarly ignored (*p* < 3.3 × 10−<sup>2</sup> ).

#### **DISCUSSION**

These results demonstrate how our CAD tool has the potential for clinically application, while also confirming and elucidating the distributed effects of epilepsy on the entire brain. Our CAD tool's diagnostic performance of TLE matches, but is not redundant with, expert MA of iPET and sMRI. When considered in the context of recent reports of CAD tools for epilepsy based on sMRI and interictal EEG data (Farid et al., 2012; Focke et al., 2012; Keihaninejad et al., 2012; Kerr et al., 2012a), CAD is proving especially applicable to epilepsy. Further, if more work confirms the hypothesis that metabolic changes in iPET are observable before the structural changes in sMRI, our iPET tool may have better clinical utility than these existing sMRI tools. In contrast to MA, this and other CAD tools can be quickly and efficiently applied by minimally trained technicians, emergency physicians, and primary care providers as preliminary analysis of the iPET images (van Ginneken et al., 2011; Kerr et al., 2012c). The performance of MA can vary with experience and fatigue of the observer; automated tools are consistent over time. Upon further validation, these CAD results could also be incorporated into the consensus diagnoses with minimal cost if iPET already has been obtained.

#### **CLINICAL IMPACT**

Our CAD tools could provide valuable clinical information that may help readily identify which treatments may be effective in patients who present with uncharacterized, and/or medication refractory seizures (Kerr et al., 2012a,c). In particular, 15 of our 105 patients were admitted twice to achieve definitive characterization or localization of their seizures. The appropriate binary classifier correctly diagnosed 12 (80%) of these challenging patients. This valuable information might reduce the need for multiple video-EEG admissions. Additionally, 28% (9/32) of our PWN were admitted for improved characterization of their previously-diagnosed "epilepsy," and 16% (12/73) of our PWE were admitted for the differential diagnosis of epilepsy, indicating that non-epileptic etiologies were not ruled out sufficiently. The trinary CAD effectively diagnosed 67% (14) of these

particularly challenging patients. Despite this impressive performance, the ultimate goal of CAD, however, is to complement – not replace – MA.

#### **COMBINATION OF AUTOMATED ANALYSIS WITH CLINICAL WISDOM**

Our finding that performance almost uniformly, but nonstatistically, decreased when the automated algorithm incorporated clinical information indicates that automated analysis cannot and should not replace manual interpretation across information modalities. We suspect that this performance decreased due to ineffective modeling of the contribution of the clinical information and over-fitting. The statistical distribution of the clinical factors was very different from the metabolic data therefore the same model likely cannot effectively utilize both modalities. The efficient incorporation of multimodality information into machine learning is an active area of theoretical research, and wellvalidated methods are not yet available. Now that CAD tools using interictal EEG (Kerr et al., 2012a), sMRI (Farid et al., 2012; Focke et al., 2012; Keihaninejad et al., 2012), and iPET have been published, we believe it will be extremely exciting to assess how these various tools can be combined.

We expected that the best performance would be achieved when our CAD is used synergistically with MA. The low correlations between the CAD results and MA suggest that our CAD tool provides information that is not evident on visual inspection. These results emphasize that PET is not redundant with MRI (Henry et al., 1999). Physicians could learn to view CAD as analogous to another imaging modality that provides valuable, but not perfectly diagnostic, clinical insight. This synergistic application of computer-aided diagnosis after manual interpretation already has proven beneficial in the detection of lung nodules by the FDA and is an active area of translational research (Kerr et al., 2012c; Wang et al., 2012). The key differences between MA and automated analysis are the ability to entirely ignore certain pieces of data, and to rule that the results are inconclusive.

The results summarized above, and the LRs for each analysis stream individually, show that both MA and CAD are useful clinically. If the analysis streams agree, the diagnostic certainty increases substantially, but at a cost: as more analyses are added, more patients have inconclusive results because the analyses did not agree, and the LRs are not significant. Even though our sample size is large compared to other studies of this type, there were not enough patients in our dataset with each diagnostic outcome to explain the clinical implication of disagreeing analyses adequately. This matter of inconclusive results is a common challenge faced in clinical practice. Physicians struggle regularly with those types of

decisions.When MA of iPET and sMRI are combined, they need to agree to yield meaningful results. However, our analysis shows that in some specific cases, if both the MA and CAD of iPET agree, the sMRI is not needed. This parallels the finding we suggested above: iPET may be more clinically useful than sMRI to diagnose and lateralized epilepsy. The hypometabolic abnormality may be present earlier in disease (Theodore et al., 1990; Ryvlin et al., 1991; Swartz et al., 1992; Gaillard et al., 1995; Debets et al., 1997; Knowlton et al., 1997, 2008; Blum et al., 1998; Drzezga et al., 1999; Benedek et al., 2004; Carne et al., 2004; Chandra et al., 2006;Yun et al., 2006; Uijl et al., 2007; Willmann et al., 2007; Rastogi et al., 2008; Salamon et al., 2008; Duncan, 2009; Lee and Salamon, 2009; Lerner et al., 2009; Liew et al., 2009; Brodbeck et al., 2010;Chinchure et al., 2010;Kim et al., 2011;Jupp et al., 2012), and it may provide slightly more accurate disease characterization, as seen in our dataset. In settings where the PET scanner is not combined with the MRI scanner, and/or when the cost of imaging is a limiting factor (both common occurrences) the effective application of our CAD could result in substantial cost savings.

#### **PATHOPHYSIOLOGICAL INSIGHTS**

Our methods also reveal a potential difference in the pathophysiology of left vs. right TLE. This may help explain why CAD tools perform slightly better when diagnosing RTLE compared to LTLE (Farid et al., 2012; Focke et al., 2012; Keihaninejad et al., 2012). The finding that mostly bilateral temporal ROIs, the right inferior frontal cortex and left sensorimotor cortex provide non-redundant diagnostic information for RTLE is consistent with the clinical wisdom that the epileptogenic network in RTLE is more focal than in LTLE. The inclusion of temporal regions echoes the conventional wisdom that focal hypometabolism and asymmetry reflect characteristic changes due to epilepsy. This suggests that conservative resection of the temporal lobe may result in increased rates of seizure freedom in RTLE compared to LTLE due to complete resection of the SOZ. Further, seizures that originate in the left temporal lobe may secondarily generalize more frequently in LTLE. These differences have not yet been studied clinically.

The trends in the extratemporal regions included in the algorithms suggest that the primary lesion may induce metabolic changes in functionally or anatomically associated regions. This is substantiated further by the finding that almost all regions of the brain provide informative diagnostic information in LTLE. This in turn mirrors the increased stereotypic connectivity of the left temporal lobe. Even though the interconnectivity of the right hemisphere is higher than the left hemisphere, the left hemisphere has strong connections between specialized foci (Barrick et al., 2005; Iturria-Medina et al., 2011; Kucyi et al., 2012). We hypothesize that the SOZ may induce abnormal metabolism along these strong, stereotyped connections. This change cannot be attributed to language specifically in our dataset because we did not identify the laterality of language dominance in our patients. Compared to our *t*-statistics ranking, it may seem surprising that the metabolism of the midbrain was ranked first by mRMR for LTLE vs. NES. This rank may indicate a non-linear change in the metabolism within the dorsal midbrain anticonvulsant zone, which has itself been identified in animals to be part of the network that modulates seizure threshold (Shehab et al., 1995). The

exact relationship between epilepsy and midbrain metabolism is unclear, however. The lack of distributed atrophy in LTLE measured by sMRI suggests that these changes are not associated with distributed cell death or gliosis (Farid et al., 2012; Focke et al., 2012; Keihaninejad et al., 2012). Instead, we hypothesize that this change instead reflects abnormal metabolism in these regions due to altered neural connectivity and/or activity secondary to the epileptogenic lesion. This is supported by the finding that LTLE was associated with more changes in functional connectivity than RTLE was (Pereira et al., 2010). This also explains why we observed metabolic changes in the right thalamus in RTLE: recent work demonstrates that the connectivity of the right thalamus with the right hippocampus is reduced in RTLE (Morgan et al., 2012). The presence of such distributed changes also supports the finding that the size of the hypometabolic lesion visualized on PET may be larger than the SOZ (Juhasz et al., 1999; Matheja et al., 2001; Henry and Roman, 2011). It is particularly interesting to note that the extent of these distributed changes is underappreciated by *t*statistics comparing LTLE to NES. This indicates that there is a complex, likely non-linear, relationship between the metabolism of the hypometabolic lesion and its associated tissue that may be better understood by mutual information.

The inclusion of the contralateral hippocampus in both of the binary classifiers lends itself to multiple interpretations that are all supported by biologically sound hypotheses. Firstly, a salient feature of LTLE or RLTE could be asymmetric metabolism, as suggested clinically; therefore the metabolism of the contralateral hippocampus was compared to the observed metabolism in the ipsilateral hippocampus. Alternatively, the interhemispheric connectivity between the hippocampi is high, therefore under our hypothesis that changes in metabolism spread according to functional connections, the metabolism in the contralateral hippocampus may be one of the first induced changes due to the epileptic lesion. Lastly, if LTLE and RTLE are inherently bilateral diseases then the metabolism in the contralateral hippocampus may also be abnormal. This also provides an explanation for why LTLE and RTLE were not perfectly distinguished.

In addition to diagnosing epilepsy, our algorithm lateralized disease efficiently with an accuracy of approximately 90% when epilepsy was diagnosed correctly. This impressive accuracy could be clinically useful for pre-surgical planning, when used in combination with other clinical and radiological information. Although our current sample size is too small to fully assess this potential fully, our results suggest that similar methodology could be applied to a larger dataset with more diverse and specific SOZ localizations to yield an objective and reliable tool to assist in pre-surgical SOZ localization. Our data suggest that this approach likely would identify and utilize distributed metabolic findings associated with each epileptic lesion to improve performance. Instead of blurring the boundary of the SOZ by detecting affected tissue outside the SOZ, the improved understanding of these distributed effects may lead to more refined characterization of this clinically vital SOZ. However, the spatial resolutions of our outcome classes were insufficient to assess the utility of this method directly to identify candidate lesions for resective surgery.

While our lateralization accuracy is exciting, there is also a potential clinical interpretation of the patients who were falsely-lateralized. Functional connectivity between the temporal lobes is particularly strong. In a minority of patients, this connectivity allows epileptogenic activity to spread quickly from the SOZ to the contralateral temporal lobe on EEG, resulting in the appearance of either bilateral or falsely-lateralized disease. Similarly to the distributed networks discussed above, this high degree of functional connectivity also may induce metabolic abnormalities in the contralateral temporal lobe that may be indistinguishable from the primary lesion. This hypothesis can be tested by comparing these falsely-lateralized patients to patients with bilateral TLE. This comparison requires a detailed methodological treatment of non-mutually exclusive classes in machine learning and therefore lies outside the scope of the current manuscript.

To characterize these and other pathophysiological insights, most studies utilize healthy neurologically normal controls. In contrast, we prefer the use of PWN as our control group. In brief, when constructing a control group, one aims to match the patients in the pathologic group in all aspects other than the pathology. In contrast to neurologically normal controls when compared to PWE, PWN's have been exposed similarly to AEDs and other medications, have increased prevalence of TBI and some other risk factors for epilepsy (Sahaya et al., 2011), have regular and frequent meetings with health care providers, and have much more strict inclusion criteria. Lastly, and perhaps most importantly, physicians do not consider whether all of their patients have epilepsy; they assess only the patients with seizures. Therefore, in our opinion, the use of PWN as the control group is a benefit in of our study because it maximizes the clinical relevance of our results while simultaneously improving its statistical selectivity.

#### **LIMITATIONS AND FUTURE DIRECTIONS**

Because our retrospective dataset was collected as part of clinical care, our approach has a few important limitations. The accuracy of MA reported in our patients is worse than the rates quoted in previous literature (Rastogi et al., 2008; Salamon et al., 2008; Lee and Salamon, 2009). Given UCLA's status as a tertiary referral center, the decrease in manual accuracy likely indicates that our patients had more heterogeneous etiologies and/or were more complex and difficult to diagnose than other centers. This suggests that our CAD tool may perform better on other datasets. Our iPETs and MRIs were collected on varying cameras with varying resolutions. This demonstrates the flexibility of our automated analysis using NeuroQ. The efficacy of the MA of older and limited resolution data may not be comparable to that of more current and higher resolution data. After establishing the efficacy of our method, we plan to both validate our tool prospectively on data from other centers, and to incorporate multi-center data into our algorithm to further improve its performance. Additionally, we only discuss the combination of CAD results with independently derived MA. Future work will examine the efficacy of CAD tools informed by MA and vice versa.

Critics of our approach might claim that the significant gender and age difference of the patients with NES compared to PWE may lead to our CAD simply detecting the age and/or gender of the patients. While we do not expect this to be the case for RTLE, the utilization of language areas by the LTLE classifier might reflect differences in gender, and not epileptogenic pathology. However, the performance of our CAD was significantly higher than when clinical information was used directly, therefore the algorithm utilized more information than just clinical data to achieve its strong performance. These significant differences in clinical factors largely mirror the observed differences in clinic; therefore our dataset better matches the population for which our CAD tool would be applied. The only notable exception is the significant age difference between LTLE and RTLE, which was unexpected. Due to the naturalistic nature of our data collection scheme, we did not correct for this difference. However, we note similarly to the NES group, the use of age alone was significantly worse than our tool and the addition of age to the iPET data to control for its effect did not significantly change performance.

Another key caveat to the direct clinical application of our tool to clinical practice is the fact that epilepsy is an extremely heterogeneous disease. The generalization of our method to bilateral TLE, extratemporal foci and multifocal epilepsy will be critical before it can be incorporated into clinical practice. In particular, even though NES mimic all types of seizures, it is uncommon for TLE to be mistaken for NES. Instead, it is more common that NES appear to have a focus in frontal cortex (LaFrance and Benbadis, 2011). Therefore, the literature suggests that the highest impact CAD tool would discriminate between frontal lobe epilepsy and NES and another, separate tool could be used to lateralize TLE. Based on our results above (see section Clinical Impact), we believe that our TLE-specific tool may be clinically applicable. For the first publication demonstrating the applicability of computer-aided diagnosis based on iPET data, we chose to focus on the diagnosis and lateralization of TLE, based, based on prior findings that the sensitivity of iPET is highest for TLE. Our future work then can address generalizing our methods to the other epilepsies, including bilateral TLE and frontal lobe epilepsy.

#### **CONCLUSION**

Despite a few caveats, and upon further validation with data from other centers, our automated methods could provide unique information for the effective and efficient characterization of epilepsy, with the potential to decrease the fraction of patients with NES that are being treated (inappropriately) with AEDs, and to more quickly triage patients with medication refractory epilepsy toward surgical intervention. This may help achieve the ultimate goal: a global reduction in seizures (Engel et al., 2012).

#### **ACKNOWLEDGMENTS**

The authors thank the reviewers for their helpful comments and would like to give special thanks to Maunoo Lee, Cheri Geist, Theresa Harrison, Henry Huang, and Klaus-Robert Müller for data management of the PET records, organizational support, neuroscientific insight, helpful comments, and machine learning expertise, respectively. This work was supported by the University of California (UCLA)-California Institute of Technology Medical Scientist Training Program (NIH T32 GM08042), the Systems and Integrative Biology Training Program at UCLA (NIH T32 GM008185), NIH R33 DA026109 (to Mark S. Cohen), and the UCLA Department of Biomathematics.

#### **REFERENCES**


al. (2006). FDG-PET/MRI coregistration and diffusion-tensor imaging distinguish epileptogenic tubers and cortex in patients with tuberous sclerosis complex: a preliminary report. *Epilepsia* 47, 1543–1549.


tomography abnormalities of thalamic nuclei in temporal lobe epilepsy. *Neurology* 53, 2037–2045.


between interictal electrical source imaging and PET hypometabolism. *Conf. Proc. IEEE Eng. Med. Biol. Soc.* 2010, 3723–3726.


Effects of locally administered excitatory amino acids or bicuculline on maximal electroshock seizures. *Neuroscience* 65, 671–679.


the use of semiautomated volumetry software. *Radiology* 262, 320–326.


univariate and multivariate random field theory. *Neuroimage* 23, S189– S195.


brain functional connectome. *Neuron* 73, 1216–1227.

**Conflict of Interest Statement:** Daniel H. Silverman is a co-inventor of NeuroQ, which is licensed to the University of California by Syntermed. The other authors have no conflicts of interest to disclose.

*Received: 18 January 2013; accepted: 18 March 2013; published online: 03 April 2013.*

*Citation: Kerr WT, Nguyen ST, Cho AY, Lau EP, Silverman DH, Douglas PK, Reddy NM, Anderson A, Bramen J, Salamon N, Stern JM and Cohen MS (2013) Computer-aided* *diagnosis and localization of lateralized temporal lobe epilepsy using interictal FDG-PET. Front. Neurol. 4:31. doi: 10.3389/fneur.2013.00031*

*This article was submitted to Frontiers in Epilepsy, a specialty of Frontiers in Neurology.*

*Copyright © 2013 Kerr, Nguyen, Cho, Lau, Silverman, Douglas, Reddy, Anderson, Bramen, Salamon, Stern and Cohen. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*

#### **APPENDIX**

**Table A1 | This table illustrates the ranking of informative regions of interest (ROIs) based on the maximum magnitude of t-statistic across the three contrasts.**



Negative t-values in the first two columns indicate the hypometabolism in epilepsy. Negative t-values in the last column indicate hypometabolism in LTLE compared to RTLE. The leading L or R indicates left or right. The lowercase letters indicate inferior (i), lateral (l), median (m), anterior (a), and posterior (p). The lagging C signifies cortex. Bold indicates significant differences at greater than the 95% confidence level. The markings of \*, \*\*, and \*\*\* indicate significance at the 95, 99, and 99.9% confidence level, respectively, without multiple testing correction. No t-statistics remain significant at the 95% confidence level after Bonferroni correction.

(Continued)

## Focal peak activities in spread of interictal-ictal discharges in epilepsy with beamformer MEG: evidence for an epileptic network?

#### **Douglas F. Rose<sup>1</sup>\*, Hisako Fujiwara<sup>1</sup> , Katherine Holland-Bouley <sup>1</sup> , Hansel M. Greiner <sup>1</sup> ,Todd Arthur <sup>1</sup> and Francesco T. Mangano<sup>2</sup>**

<sup>1</sup> Division of Neurology, Department of Pediatrics, Cincinnati Children's Hospital Medical Center, Cincinnati, OH, USA

<sup>2</sup> Division of Neurosurgery, Department of Neurosurgery, Cincinnati Children's Hospital Medical Center, Cincinnati, OH, USA

#### **Edited by:**

Don Tucker, Electrical Geodesics, Inc. and the University of Oregon, USA

#### **Reviewed by:**

Marino M. Bianchin, Universidade Federal do Rio Grande do Sul, Brazil Giridhar Padmanabhan Kalamangalam, University of Texas Health Science Center, USA

#### **\*Correspondence:**

Douglas F. Rose, Division of Neurology, Department of Pediatrics, Cincinnati Children's Hospital Medical Center, 3333 Burnet Avenue, ML11006, Cincinnati, OH 45229, USA. e-mail: douglas.rose@cchmc.org

Non-invasive studies to predict regions of seizure onset are important for planning intracranial grid locations for invasive cortical recordings prior to resective surgery for patients with medically intractable epilepsy. The neurosurgeon needs to know both the seizure onset zone (SOZ) and the region of immediate cortical spread to determine the epileptogenic zone to be resected.The immediate zone of spread may be immediately adjacent, on a nearby gyrus, in a different lobe, and sometimes even in the contralateral cerebral hemisphere.We reviewed consecutive simultaneous EEG/MEG recordings on 162 children with medically intractable epilepsy. We analyzed the MEG signals in the bandwidth 20–70 Hz with a beamformer algorithm, synthetic aperture magnetometry, at a 2.5 mm voxel spacing throughout the brain (virtual sensor locations, VSLs) with the kurtosis statistic (g2) to determine presence of excess kurtosis (γ2) consistent with intermittent increased high frequency spikiness of the background. The MEG time series was reconstructed (virtual sensor signals) at each of these VSLs.The VS signals were further examined with a relative peak amplitude spike detection algorithm.The time of VS spike detection was compared to the simultaneous EEG and MEG sensor signals for presence of conventional epileptiform spike morphology in the latter signals.The time ofVS spike detection was compared across VSLs to determine earliest and last VSL to show a VS spike. Seven subjects showed delay in activation across VS locations detectable on visual examination. We compared the VS locations that showed earliest and later VS spikes with the locations on intracranial grid locations by electrocorticography (ECoG) that showed spikes and both onset and spread of seizures. We compared completeness of resection of VS locations to postoperative outcome. The VS locations for spike onset and spread were similar to locations for ictal onset and spread by ECoG.

**Keywords: magnetoencephalography, beamformer, children, adolescents, intracranial EEG, outcome, network, localization**

#### **INTRODUCTION**

Non-invasive studies to predict regions of seizure onset are important for the planning of intracranial grid locations for epilepsy. The surgeon needs to know both the seizure onset zone (SOZ) and the region of immediate cortical spread to determine the epileptogenic zone (EZ) to be resected (Engel, 1996; Wiebe et al., 2001; Luders et al., 2006). The immediate zone of spread may be immediately adjacent, on a nearby gyrus, in a different lobe, and sometimes even in the contralateral cerebral hemisphere. Among non-invasive recording modalities, magnetoencephalography (MEG) has some advantages because of its very good time resolution and good spatial resolution, as signals are not altered by the skull (Knowlton et al., 2009; Stefan et al., 2011). Prediction of intracerebral locations of sources for recorded extracranial MEG signals requires mathematical source localization algorithms. Multiple categories of algorithms have been developed including single and multiple equivalent current dipoles (ECD), dipole scans such

as multiple signal classification (MUSIC) (Mosher, 1999), distributed dipoles (current density) such as minimum norm estimate (MNE) (Tanaka et al., 2010), standardized low resolution brain electromagnetic tomography (sLORETA) (Pascual-Marqui, 2002), vector and scalar beamformer (Robinson et al., 2004; Sekihara et al., 2005), and multiple other algorithms each with particular strengths. Source localization algorithms in the beamformer category are spatial filters and have very good signal to noise resolution. The beamformer algorithm, as a spatial filter, can "tune" the MEG sensor array to enhance the signal arising in a specific location while diminishing brain signals arising from other locations.

The tuning process has been referred to as creating a "virtual sensor" (VS) at the chosen spatial location (VS location) and the resulting spatially filtered signal is referred to as a VS signal. The tuning can be repeated iteratively for multiple locations intracranially on the same recording of the MEG signal, such that the whole brain can be scanned sequentially to evaluate source contributions to the overall signal that was recorded at the magnetometer sensors. We chose to evaluate the capability of a particular scalar beamformer algorithm, synthetic aperture magnetometry (SAM), to detect onset and spread of interictal spikes, interictal bursts of spikes, and ictal discharges (Robinson et al., 2004).

#### **METHODS**

#### **SUBJECTS**

We initially reviewed consecutive simultaneous EEG/MEG recordings between January 2006 and December 2008 for 162 children and adolescents with medically intractable epilepsy who were admitted for non-invasive Phase I presurgical evaluation. All subjects were recorded during spontaneous sleep after sleep deprivation, conscious sedation with chloral hydrate, or general anesthesia with dexmedetomidine (**Table 1**). This study was performed under Cincinnati Children's Hospital Medical Center Internal Review Board Protocol 2008-0403 as a retrospective review.

#### **ELECTROMAGNETIC RECORDINGS**

All subjects had simultaneous 23-channel scalp EEG and 275 channel MEG recorded with a whole-head CTF 275-channel magnetometer (VSM MedTech Systems Inc., Coquitlam, BC, Canada). Fiducial markers were placed at the nasion and left and right preauricular points and placements were photographed digitally. The same fiducial locations were marked with MRI visible targets. The MRI scans were thin slice, 1 mm in the sagittal plane. Each subject had a series of simultaneous MEG/EEG recordings for 2 and 10 min durations for a minimum recording time of 40 min. All subjects' recordings were digitized at 300 or 600 Hz for the 10-min recording and additionally at 4000 Hz for the 2-min recordings.

#### **DATA ANALYSIS**

We first visually reviewed each EEG/MEG recording in the bandwidth 1–70 Hz for ictal and interictal discharges, and removed any sections of the recordings with non-cerebral artifact, particularly sections containing muscle artifact. We further analyzed the remaining sections of MEG signals in the bandwidth 20–70 Hz with a beamformer algorithm, SAM, at a 2.5 mm voxel spacing throughout the brain (VS locations, VSLs) at approximately 300,000 locations in the three dimensional grid. We applied the kurtosis statistic (*g* <sup>2</sup>) to the signals at each VS location to determine presence of excess kurtosis (γ2) consistent with intermittent increased high frequency spikiness of the background (Ukai et al., 2004; Kirsch et al., 2006; Ishii et al., 2008; Prendergast et al., 2013). The kurtosis statistic is a measure of the peakedness and thickness of the tails of a unimodal distribution (DeCarlo, 1997). A positive kurtosis indicates more "outliers" in the tails of the distribution than expected. For this application, the "outliers" are intermittent additional high frequency activity within a constrained 20–70 Hz bandwidth compared to the background brain noise distribution in that bandwidth (Robinson et al., 2004). The VS locations were ordered based on highest to lowest excess kurtosis. The first 20 VS locations with highest excess kurtosis were chosen for further examination. A short note of explanation: Because the signal at each VS location is first evaluated for kurtosis separately from all others, a VS location with very infrequent additional high frequency activity may register a higher kurtosis than locations with

more frequent, yet still intermittent, additional high frequency activity occurrences. Also, the VS location with the highest kurtosis may not have the most frequent VS spikes or the earliest occurring of the spike components. Therefore, once the VS locations demonstrating increased kurtosis are identified, each location must still be evaluated to assess the nature and frequency of occurrence of the increased high frequency activity. For the further analysis, the MEG time series (VS signal) was reconstructed for each of these VS locations. The continuous VS signals for each VS location were evaluated with a simple adaptive spike detection algorithm based on peak amplitude relative to prior signal amplitude. An additional note here: sometimes a VS location with increased kurtosis does not show a spike in the reconstructed signal. In that case the increased high frequency activity may be low amplitude, infrequent and occur over a longer time span than the short time span of a spike. Thus, although increased activity is present, because it is more spread out over time, the activity does not rise up above the background as a visible spike. The opposite effect can also occur if spikes occur very frequently, for instance 2/s throughout the recording, which has happened in some of our pediatric patients. Then the spikes are a common, consistent feature of the background, contribute to the central "peak" of the distribution, and do not increase the kurtosis at the tails. Also, since we measured the kurtosis in the bandwidth 20–70 Hz, spikes in the original MEG/EEG signal, whose spectral power was confined below 20 Hz (this can happen in very young children), would not be detected by this method. Thus the reconstructed VS signals must always be compared to the simultaneous original MEG/EEG signals. The EEG, MEG, and the VS signals were displayed together as channels, with the VS signal at each VS location displayed as a separate VS channel (**Figure 3A** shows simultaneous EEG and VS channels, but not the MEG sensor channels because of space limitations). The time of each VS spike detection in each VS location (channel) was compared to the simultaneous original signals recorded at the EEG electrodes and MEG sensors (EEG/MEG sensor signals) for the presence of conventional epileptiform spike morphology in the latter signals. Only the VS spikes that were associated with conventional epileptiform spikes in the EEG/MEG sensor signals were chosen for further analysis and evaluation.

For each VS location, the peak of the VS spike was used as the time mark to average the spikes at that location. A time window 250 ms before and after the VS spike peak was chosen to make a 500-ms total duration epoch for averaging. The epochs were visually reviewed before averaging, including the simultaneous original EEG and MEG sensor signals. Only epochs that contained just one epileptiform discharge in the EEG/MEG sensor signals were included in the average. The simultaneous EEG, MEG, and VS signals from the other VS locations that occurred during that epoch were averaged at the same time based on the same time mark (**Figure 1A**). The number of VS spikes that occurred in the VS signals at a VS location varied from 2 to about 200, and in general depended on the sparseness or abundance of epileptiform discharges for that particular subject.

For the purposes of this article, the averaged VS spike at the VS location, along with the simultaneous averaged EEG and MEG sensor signals and the VS signals from the other VS locations, was considered a VS spike "type." This definition of a "type" is based


not done. \*A small amplitude deflection was seen shortly after VS spike at earliest VS location although the most prominent VS spike occurred at later timing.

**Frontiers in Neurology** | Epilepsy May 2013 | Volume 4 | Article 56 |

**FIGURE 1 | Subject 1**. **(A)** EEG in bipolar AP montage and five virtual sensor channels showing average of 146 interictal spike epochs in bandwidth 3–70 Hz, digitized at 300 Hz. Topographic maps for MEG and EEG signals are top of head views centered at vertex Cz with anterior (nose) at top of map. For this averaged spike, epochs were aligned based on spike peak detected at VS location V3 in right frontal lobe in the 20–70 Hz bandwidth. **(B)** Inset is enlarged view virtual sensor (VS) channels showing timing of peaks. V5 peak occurs earliest followed by V4 (40 ms), V3 and V2 (50 ms), and V1 (63 ms). Note that although largest amplitude peaks at V1–V4 are delayed after V5 peak, each has a small peak at about the time of V5 peak. Onset of spike peak at V5 location precedes apparent onset at the other four VS locations. Although the averaged epochs shown were aligned based on the spike at V3,

on a signal occurring at an intracranial location following spatial filtering, although it does include the simultaneous EEG and MEG sensor waveforms and VS waveforms forms from the other VS locations to further characterize the spike "type." The definition is somewhat different from the conventional definition of an EEG or MEG spike type that is based on the pattern of simultaneous EEG electrode or MEG sensor signals without any spatial filtering.

Since more than one VS location for a subject might have spikes detected in the VS signal, the procedure described above was repeated for the VS spikes at each VS location. Thus, if a subject had five VS locations that had VS spikes, the procedure would be repeated five times. The result would be five sets of averaged VS spikes (five VS spike types) that included the simultaneous EEG/MEG sensor signals and also the simultaneous VS signals from the other locations.

After averaging, the EEG, MEG, and VS averaged spikes were examined for earliest onset of the EEG, MEG, or VS spike and also for time of the waveform peak (**Figure 1A**). If more than one of the VS locations in an averaged epoch showed VS spikes, the time of VS spike onset and spike peak was compared across each of the VS locations to determine which VS location showed the first spike and the timing difference for EEG, MEG, and VS spikes at the other VS locations (**Figure 1A**). The averaged epochs from each of the VS locations were compared to determine whether the timing differences for the EEG, MEG, and VS spikes at all VS locations remained consistent.

For subjects showing a consistent timing difference among spikes across VS locations, the unaveraged EEG, MEG, and VS signals were then reviewed again at the time in the original EEG/MEG recording at which the VS spikes occurred to be certain that the timing differences seen in the averages of the EEG/MEG and VS spikes were also present in the unaveraged EEG, MEG, and VS signals.

For the subjects who had interictal or ictal spike bursts in the conventional EEG/MEG signals instead of individual spikes, the simultaneous EEG/MEG sensor signals and VS signals were compared for earliest time of burst onset, but no averaging was done across bursts.

#### **RESULTS**

All subjects had one or more VS locations identified by the SAM + *g* <sup>2</sup> [SAM(*g* <sup>2</sup>)], evaluation for excess kurtosis. Although the VS locations for each subject had been detected because of excess kurtosis, the subsequent analysis with the adaptive spike detection the earlier onset of a spike at V5 was detected. **(C)** Subjects MRI scan showing VS locations of V1–V5. Green dot is location of V5, earliest of the VS peaks. **(D)** Subject's MRI scan showing one possible pattern of spread based on timing of peaks and shortest distances between VS locations. The arrows are only to show relative timing, since propagation and pathways cannot be discerned solely with spike peaks separated in time. **(E)** 3D reconstruction of subject's MRI and subdural electrode locations. **(F)** Intraoperative photograph of left brain surface. Green dots represent electrodes showing spikes at onset of seizures. Ictal onset was in left frontal lobe at LIPF2–4 (red frontopolar electrodes and purple anterior frontal lobe electrodes) with rapid spread to left anterior temporal lobe (blue electrodes). Yellow overlay represents resection plan.

algorithm did not always detect VS spikes in the VS signals at all the VS locations.

Each VS location/channel that had VS spikes detections that were used to provide the marker for averaging a set of epochs of simultaneous EEG, MEG, and VS signals did show a spike in that location in the bandwidth 20–70 Hz following averaging the epochs, as would be expected. Usually a VS spike was also seen at that VS location in a wider bandwidth filter, such as 3–70 Hz. The other simultaneous VS signals/channels in an averaged epoch also sometimes showed VS spikes.

The VS spikes in the averaged epochs often showed an onset and peak before the simultaneously averaged EEG and MEG sensor spikes. For those averaged epochs that showed VS spikes at multiple VS locations, the timing of onset and peak of the VS spikes by visual inspection appeared near simultaneous across the VS locations, differing by often only 5–20 ms of each other. For a subset of seven subjects the timing of the peaks of VS spikes was longer than 20 ms and more easily discernible by visual inspection.

We studied in greater detail the seven subjects who had the relatively longer time difference among VS spike activation across their VS locations as we felt the slightly greater time separation could be more clearly analyzed. These seven subjects ranged in age from 1 year, 10 months to 16 years, 1 month (mean 9.9 years). Epileptiform patterns in the raw EEG and MEG signals included interictal spikes (three subjects), primarily interictal bursts (two subjects), mixed spikes and bursts (one subject), and both interictal and ictal bursts (one subject). Etiology was probable cortical dysplasia on imaging (one subject), prior brain tumor (one subject), tuberous sclerosis (two subjects), probable encephalitis (one subject) or undetermined at the time of study (two subjects). Of the seven subjects, five had Phase II intracranial monitoring with strip and grid electrode recordings and surgical resection; one subject had corpus callosotomy, and one subject chose not to have any surgical procedure. All six subjects with surgical treatment had a minimum of 2 years follow-up. The seven subjects' ages, presumed etiologies, epileptiform discharge types and locations, MEG findings, electrocorticography (ECoG) findings, surgical treatments, outcomes, and pathological findings (when available) are summarized in **Table 1**.

#### **SUBJECT 1**

On video/EEG, the subject had frequent moderate voltage spikes and sharp waves in the left temporal, left frontal, and right frontal head regions. A clinical complex partial seizure was recorded. EEG findings suggested an ictal onset zone at the left frontotemporal region with symptomatic zone at the left temporal head region. On combined MEG/EEG, the subject had 187 spikes and sharp waves that occurred topographically in the left middle temporal head region. SAM(*g* <sup>2</sup>) detected 5 VS locations: left frontal lobe (73 VS spikes averaged), left frontopolar lobe (85 VS spikes averaged), left temporal lobe (107 VS spikes averaged), right frontal lobe (146 VS spikes averaged), right frontopolar lobe (26 VS spikes averaged). The VS location with earliest VS spike was left frontopolar with spread both to right frontal and to left temporal and then to left inferior frontal (**Figure 1**).

Following craniotomy, subdural electrodes were placed over the lateral cortical surface for left frontal and temporal lobes, left orbitofrontal, and left frontopolar. Ictal ECoG demonstrated two main foci: one in the inferior prefrontal and posterior frontal lobe. In addition, interictal ECoG showed very frequent epileptiform discharges at the left anterior temporal region, which was involved during the two clinical seizures. The subject had a left frontal corticectomy sparing the opercular portion of inferior frontal gyrus and a left anterior temporal lobectomy. At 2 years postoperative follow-up, the subject was seizure free.

#### **SUBJECT 2**

On video/EEG frequent spikes occurred in the right hemisphere at O2, P4, and C4. Nine ictal events suggested right parieto-occipital onset. On combined MEG/EEG, the subject showed frequent right parieto-occipital spikes and frequent bursts of polyspikes in the central head regions that were bilateral but more prominent over the right hemisphere. SAM(*g* <sup>2</sup>) analysis of the recording revealed 7 VS locations. Three that showed VS spikes included: right insula (99 VS spikes averaged), right superior posterior frontal (110 VS spikes averaged), and right superior mesial anterior parietal (145 VS spikes averaged). The VS location showing earliest spike activation was right mesial parietal with almost synchronous spread to right posterior frontal and right insula For bursts, the earliest activation was right insula with almost synchronous spread to right posterior frontal and right mesial parietal (**Figure 2**).

Following craniotomy, subdural electrodes were placed over the lateral cortical surface for right temporal and right frontoparietal with one superior frontal strip. Earliest seizure onset was recorded at the superior temporal electrodes with spread to inferior posterior parietal and later posterior superior frontal. The subject had a right temporal lobectomy and right parietal corticectomy, which on pathology showed diffuse cortical dysplasia. The subject continued to have seizures and required a right hemispherectomy, but has been seizure free at 2 years follow-up.

#### **SUBJECT 3**

On video/EEG, the subject had frequent multifocal bilateral discharges and multiple clusters of electroclinical seizures that were tonic,atonic,or myoclonic,but electrographically the seizures were similar regardless of semiology. Initial changes appeared to occur in the left posterior region slightly before diffuse changes occurred, suggesting that these could be rapidly generalizing partial onset seizures. MRI imaging had revealed a lesion in the left mesial occipitoparietal region felt to be consistent with cortical dysplasia. On combined MEG/EEG recording, the subject showed multifocal spikes bilaterally and frequent diffuse bursts. SAM(*g* <sup>2</sup>) identified

15 VS locations that showed spike activity during the bursts. These were not averaged, but each burst was individually examined. The earliest onset was in the VS location in the left occipitoparietal junction, just lateral to the MRI lesion. There was spread to VS locations nearby and somewhat later to the left posterior frontal lobe (**Figure 3**).

Following craniotomy, subdural electrodes were placed over the lateral cortical surfaces for left temporal, parietal, occipital, and posterior frontal lobe. ECoG showed numerous clinicalelectrographic seizures arising from the left occipital pole and posterior subtemporo-occipital region spreading both superiorly into the left parietal region and anteriorly into the left temporal region. The surgical resection was left parieto-occipital and posterior temporal. The pathology of the MRI lesion was focal cortical dysplasia. The subject has been seizure free for 2 years.

#### **SUBJECT 4**

This subject had a brain tumor resected in the right frontal lobe 6 years before presurgical epilepsy evaluation and the MEG study. On video/EEG, the subject had frequent high voltage spikes and sharp waves in the left frontal, right anterior temporal, and right frontal head regions. He had 27 events recorded in the EMU, 26 of which were electroclinical seizures: 17 were tonic or generalized tonic clonic seizures (8 localized left frontal head region, 3 to left hemisphere, 6 showed no clear lateralization), 7 were complex partial seizures (5 showed a right frontal localization, 2 did not lateralize), 2 were of atypical semiology and did not show any clear lateralization. On combined MEG/EEG, the subject had 98 spikes seen in the right frontal head region in the EEG recording. The raw MEG signal was obscured by magnetic noise from the subjects vagal nerve stimulator (VNS). SAM(*g* <sup>2</sup>) detected 11 VS locations all in the right hemisphere. Four VS locations were in right frontal lobe (5, 6, 7, and 8 VS spikes), two in right posterior frontal (1 and 7 VS spikes), two in right lateral parietal (4 and 7 VS spikes), and three in right mesial parietal (1 VS spike in each VS location). The spikes in the right frontal lobe preceded those in the right mesial parietal lobe by about 65 ms (**Figure 4**).

Following craniotomy, subdural electrodes were placed over the lateral cortical surfaces for right frontal and anterior parietal lobe, and strips over lateral and inferior temporal lobe. ECoG recorded five clinical seizures characterized by asymmetric tonic posturing and one myoclonic seizure. The seizure onset was localized to the right mid and anterior frontal region with rapid spread to right anterior and mid temporal regions. A right frontal resection was performed. Postoperatively the subject had recurrent seizures.

#### **SUBJECT 5**

On video/EEG, the subject had occasional spike and slow waves in the left frontocentral and left centroparietal head regions. Twelve tonic seizures were recorded: nine began in the right temporal head region, one in the right centroparietal region, one in the left centroparietal region, and one was poorly localized. On combined MEG/EEG, the subject had multifocal interictal spikes: left frontal, left frontocentral, left centroparietal, right centroparietal, and right frontocentral head regions. SAM(*g* <sup>2</sup>) detected multiple VS locations bilaterally with VS spikes in left lateral superior frontal lobe (21 VS spikes averaged), left lateral superior parietal

occurred first followed by right mesial parietal V5 (80 ms), then right mid frontal V7 (240 ms), and right superior mid frontal V3 and V1 (260 ms), and finally right posterior frontal V2 (300 ms). **(B)** Subject's MRI scan showing one possible pattern of spread based on the timing of the VS peaks.

indicate the ictal onset (1) and spread pattern based on ECoG (2, 3). **(E)** Intraoperative photograph of brain and subdural grids with ictal onset right mid temporal (dark purple), spread to inferior parietal (pink), then superior posterior frontal (dark blue).

locations in this bandwidth (**B)** The same epileptiform discharge but shown in bandwidth 20–70 Hz. Earliest onset at VS location V10 left parietal followed by left posterior frontal V9 (43 ms), and right parietal V4 (43 ms), then right posterior frontal lobe V1 (173 ms). **(C)** Subject's sagittal MRI scan showing VS

relative timing of peaks. Green circle indicates location of left parietal MRI lesion. **(E)** 3D reconstruction of subject's MRI scan. Purple dots indicate the location of spikes on subdural grid electrodes preceding ictal fast activity. Posterior to red line indicates resection plan.

lobe (anterior – 5 and posterior – 17 VS spikes averaged), right frontal lobe (10 VS spikes averaged), right parietal lobe (lateral anterior superior – 72 and lateral mid superior – 20 VS spikes averaged). Based on the order of occurrence of VS spikes in the right hemisphere, the anterior and posterior parietal lobe locations showed the earliest spike, with spread into the right posterior parietal lobe and right lateral midfrontal lobe VS spike. In addition, the order of occurrence of VS spikes suggested there was simultaneous spread from the right posterior parietal to left posterior parietal and then to left anterior parietal, left posterior frontal, and left lateral midfrontal lobe locations. A second pattern of spread (again based on VS spike order of occurrence, but also based on averaging left hemisphere interictal spikes) suggested that some VS spikes began in the left lateral midfrontal lobe and spread to left

showing interictal spike in unaveraged signals in bandwidth 3–70 Hz, digitized at 300 Hz. The earliest peaks occurred at V4 and V6 right frontal lobe (red arrows), just posterior to the resection site for tumor, followed 100 ms later by spike peaks at V9 and V10 right parietal lobe (blue arrows). In an expanded timebased view (not shown), V9 and V10 also show smaller amplitude spikes earlier, 20 ms after V4 and V6. In that view, timing of the spike at V3, right frontal lobe, is also 20 ms after V4, V6. **(B)** Subject's MRI

at V4, V6 (red circle), and subsequent spike peaks anteriorly at V3 and posteriorly at V9, V10. Red arrows indicate relative timing of peaks. **(C)** 3D reconstruction of subject's MRI scan with subdural grids over right frontal and temporal lobe. **(D)** Summary of electrodes active at or shortly after seizure onset with the most intense red color corresponding to the most active electrodes in superior, mid, and inferior frontal lobe posterior to resection site and also mid temporal.

posterior frontal lobe, then to left anterior parietal lobe locations. This left lateral midfrontal lobe VS location was both the last VS location showing spread from the right parietal spikes and the first VS location showing activation and spread of VS spikes for spikes originating in the left hemisphere (**Figure 5**).

Based on the multimodality presurgical studies in addition to MEG, a right craniotomy was performed with electrode grids placed over the right frontal, temporal, and parietal lobes. Five typical seizures were recorded with onset detected in the superior right frontal lobe on the right mid superior frontoparietal grid, frontal grid, and middle interhemispheric strip. A right posterior frontal lobe resection was performed with multiple subpial resections over the motor cortex. Initially postoperatively the subject was seizure free. However, within 1 week seizures recurred. On repeat video/EEG, the subject showed a different clinical pattern of flexing the right elbow and flexing the back with EEG onset in the left central head region. Combined MEG/EEG was repeated at that time. The subject had fewer epileptiform discharges per unit time compared to the first study. The most frequent detections of VS spikes were in the left lateral midfrontal lobe location described above that preoperatively showed late activation after right parietal VS spikes and earliest activation of VS spikes for left hemisphere epileptiform discharges. The family declined a second operative procedure, and a second Phase II study was therefore not indicated. At 2 years postoperative follow-up, the subject continued to have intermittent seizures.

#### **SUBJECT 6**

On video/EEG, the subject had frequent myoclonic and atonic seizures with diffuse electrographic changes that were not welllocalized. On combined MEG/EEG, the subject had very frequent bilaterally diffuse and synchronous spikes and polyspikes. SAM(*g* <sup>2</sup>) detected multiple VS locations in the parietal and occipital lobes and fewer in the temporal and posterior frontal lobes bilaterally. Two VS locations that showed prominent VS spikes were left (154 VS spikes averaged) and right (129 VS spikes averaged) posterior parietal lobes. The left parietal VS spike usually preceded the right parietal VS spike by about 25 ms (**Figure 6**). Based on all of the subject's tests, the subject was judged clinically not to be a candidate for focal resection, did not have a Phase II study, but did have a corpus callosotomy. At 2 years follow-up, the subject continued to have seizures, but not the drop attacks.

#### **SUBJECT 7**

On video/EEG, the subject had frequent spikes and sharp waves in the left anterior temporal region head region. Frequent bursts of high voltage diffuse spikes and polyspikes and slow waves occurred with a frontocentral predominance. The subject had brief tonic seizures lasting 3–5 s and multiple complex partial seizures, but with somewhat variable semiology. All clinical seizures were poorly localized or lateralized on scalp EEG. On combined MEG/EEG, the subject had frequent focal spikes bilaterally and multiple bursts of polyspikes seen on EEG, six of which were associated with a bilateral tonic "shudder" of the arms. The MEG sensor signals could not be directly interpreted because of the magnetic noise from the subject's VNS. The spatial filter characteristics of the beamformer algorithm excluded this distant noise source in the SAM VS. SAM(*g* <sup>2</sup>) detected VS locations in left frontal lobe (424 VS spikes averaged), right frontal lobe (472 VS spikes averaged), and right opercular (106 VS spikes averaged), right parietal (118 VS spikes averaged), and occipital lobes (122 VS spikes averaged). No consistent pattern of spread from one VS location to another (based on VS spike onset or peaks) was found for averaged VS spikes. For the six brief clinical seizures, sometimes the earliest VS spikes (in the VS signals associated with the EEG polyspike bursts) was in the right hemisphere and at other times in the left hemisphere in the 20–70 Hz bandwidth (**Figure 7**). The VS waveforms at locations of peak spikiness suggested that the location of earliest onset during the electroclinical seizures of polyspike bursts varied and could begin in either the right or left hemisphere. The subject chose not to proceed with surgery, and Phase II evaluation was therefore not done.

#### **PATTERNS OF VS SPIKE OCCURRENCE DURING EEG/MEG SENSOR INTERICTAL AND ICTAL DISCHARGES**

The spatial separation between VS locations showed that spread of activation varied from intralobar (five subjects) through intrahemispheric (seven subjects) to interhemispheric (five subjects), with six subjects showing more than one distance of spread. Timing for spread intralobar ranged from 20 to 63 ms, intrahemispheric was 39 to 300 ms, and interhemispheric was 26 to 173 ms (**Table 1**; **Figures 1**–**7**). Usually, the VS location that showed the earliest visually noticeable change from baseline also showed the earliest peak activity; however,sometimes twoVS locations showed approximately similar first change from baseline, but one of the two VS locations might demonstrate an earlier peak. For some subjects, the VS locations that showed delayed spike peaks sometimes also showed a small amplitude deflection shortly after the VS spike at the VS location showing earliest VS spike peak.

#### **COMPARISON OF PREOPERATIVE MEG BEAMFORMER PUTATIVE SOURCE TIMING AND INTRACRANIAL ECOG INTERICTAL/ICTAL RECORDINGS**

For comparison of preoperative MEG beamformer predictions and ECoG findings, two of the seven subjects did not proceed to Phase II evaluation with ECoG, so for those two subjects no comparison between MEG and ECoG can be made.

For subject 1, MEG beamformer showed left frontopolar onset with spread to left temporal lobe, both of which were detected on ECoG and in that order for ictal onsets. However, later spread to the left inferior frontal was not detected, which could be MEG beamformer error, or possibly the left frontal grid did not extend sufficiently inferiorly over frontal lobe. ECoG could neither validate nor invalidate occurrence or relative timing of VS spikes in the right frontal lobe, since no intracranial electrodes were placed over that hemisphere. The subject has been seizure free for 2 years without any resection in right frontal lobe.

For subject 2, MEG beamformer detected an onset in right temporofrontal operculum or insula and subsequent occurrence in right mesial parietal, right mid frontal, and superior frontal. ECoG, on the surface of the frontal and temporal lobes, detected an onset in right middle temporal gyrus, then inferior parietal, and finally anterior superior parietal or posterior superior frontal. Overall both preoperative MEG beamformer and ECoG showed origin

**FIGURE 5 | Subject 5**. **(A)** EEG in average reference montage and seven virtual sensor channels showing average of 72 interictal spike epochs in bandwidth 3–70 Hz. For this averaged spike, epochs were aligned based on spike peak detected at VS location 1 in right parietal lobe in the 20–70 Hz bandwidth. The earliest onset occurred in V1 about 33 ms before the first onset in EEG channels C4 and P4. The panel also shows the first peak in V1 preceded the peak in V3 right mesial parietal and V6, right posterior frontal. Averages of epochs based on the timing of the spikes at each of the other six VS locations (not shown) were compared to obtain a composite of spike peak timing differences among all VS locations. **(B)** Subject's MRI scan showing the composite of spike peak timing differences. Red arrows indicate relative timing of peaks. The average of epochs based on MEG spike detections in the 20–70 Hz bandwidth for VS location V7 showed an additional pattern of spike peak timing that tracked back from left posterior frontal to anterior parietal, then posterior parietal, but did not cross the midline to the right hemisphere.

near/in temporal lobe with subsequent spread dorsally to parietal and then superior parietofrontal. No intracranial recording, stereotaxic/stereotactic (SEEG), was done of the insula/opercula nor ECoG of superior mesial parietal, so intracranial recordings did not validate or invalidate either location as an onset site or terminal extension of ictal activity as suggested by MEG interictal VS fast activity burst analysis. The resection was limited to temporal and inferior parietal, but the subject continued postoperatively to have seizures. Subsequently the subject had right hemispherectomy and thereafter has been seizure free for 2 years.

For subject 3, MEG beamformer showed left parietal onset in sulci just lateral to the MRI lesion followed both left posterior frontal and right parietal to right posterior frontal spread. ECoGdetected spikes in the left posterior parietal and occipital surfaces and rapid diffuse spread into the left anterior parietal and posterior temporal lobes, but not into left posterior frontal lobe. No intracranial electrodes were placed over the right hemisphere, so spread and timing there could not be validated/invalidate by the intracranial recordings. The resection was limited to left occipitotemporoparietal lobes, but the subject was 2 years seizure free.

For subject 4, MEG beamformer showed a frontal onset posterior the prior resection site and both more anterior spread mesially in frontal lobe and posterior spread to deep in the parietal lobe. ECoG showed ictal onset posterior to the lesion in frontal lobe but did not detect the parietal spread shown in MEG and frontopolar surface intracranial electrodes were not placed, perhaps because of scarring from the prior tumor resection. After a prefrontal lobe resection, 4 cm from the frontal tip, the subject continued to have seizures.

For subject 5, MEG beamformer detected a right frontoparietal onset but also a separate left posterior frontal onset. The other modality tests done clinically preoperatively also pointed to right frontoparietal region. ECoG over right temporofrontoparietal regions showed superior right parietal onset. No ECoG was done over the left hemisphere, so a second left posterior frontal onset could not be validated/invalidated. After a right posterior frontal lobe resection with multiple subpial resections over the motor cortex, the first ictal clinical pattern disappeared, but new seizures recurred suggesting left hemisphere onset based on semiology. A postoperative MEG showed left posterior frontal and **(C)** 3D reconstruction of the subject's MRI scan showing subdural electrode locations over right frontal, parietal, and temporal lobes. Ictal onset was in superior anterior parietal at the yellow X. [The red regions here represent region of activation from subtraction ictal SPECT co-registered with MRI (SISCOM) but were not analyzed as part of the current study.] **(D)** Combined EEG/MEG study performed 2 months postoperatively because subject continued to have seizures. Bipolar AP EEG and 10 virtual sensor channels showing average of 3 interictal spike epochs in bandwidth 5–70 Hz, digitized at 4 KHz. The peaks at VS locations 1–2 and 8–10 all seem to occur at approximately the same time. Spikes were seen in the MEG signals, but not well seen in the EEG. **(E)** The subject's MRI scan showing the VS locations from the postoperative MEG. Highest excess kurtosis was at left posterior frontal near V1 and V2 (red region). This location was also near the location of the second spread pattern seen preoperatively (Panel B) in left posterior frontal (V7 location in that panel).

parietal VS spikes that were not well seen in simultaneous 10– 20 scalp EEG. No further intracranial EEG studies were planned as a second resection was not planned, so the persistent left posterior frontal MEG spikes could not be confirmed as the ictal onset location in the left hemisphere.

#### **DISCUSSION**

This study examined source localization for epileptiform discharges with a beamformer algorithm that spatially filtered the MEG signals by "tuning" the MEG sensor array to provide a single VS signal at a single VS location. Each 2 or 10 min MEG recording was iteratively examined offline at multiple intracranial locations for source activity. The source activity at each location was evaluated with the *g* <sup>2</sup> excess kurtosis statistic.

Virtual sensor signals at the VS locations with highest excess kurtosis were evaluated for transient spikes. When spikes were detected, epochs surrounding the spikes were averaged including the simultaneous EEG/MEG sensor spikes and simultaneous VS signals at up to the 20 VS sensor locations with the highest excess kurtosis.

In review of all 162 subjects, the MEG VS spikes often had their onset at, or sometimes just before, the onset of the EEG/MEG sensor spikes (see **Figures 2** and **5**) and sometimes were shorter in duration than the associated EEG/MEG sensor spikes (see **Figure 6**). The finding of an earlier VS onset may have occurred because the spatial filtering lowered the contribution of background noise from other sites in the brain, improved the signal to noise ratio, and allowed a smaller amplitude, earlier occurring spike signal to be detected above the lowered background noise level. The shorter duration spike may indicate that activity at that location had ceased and propagated to another site; alternatively, the finding may indicate local contiguous spread created a more extended source that fit less well the forward solution of a focal dipole source, and lead to a lower amplitude VS signal.

For most of our 162 subjects whose averaged epochs contained VS spikes at multiple VS locations, and the timing of onset and peak of the VS spikes across the VS locations was similar. Based on visual inspection, these VS spike features appeared to occur either synchronously or varied by just 5–20 ms in timing. For a subset of seven subjects the timing of VS spike onsets and peaks differed more, and differences were discernible with visual inspection. Although time for a signal to cross from one cerebral hemisphere

channels showing average of 154 interictal spike epochs in bandwidth 3–70 Hz. For this averaged spike, epochs were aligned based on spike peak detected at VS location V1 in left occipital lobe in the 20–70 Hz bandwidth. **(B)** Inset shows close-up of VS peaks at VS locations. Although the averaged epochs shown were aligned on the timing for VS location V1, the average shows the earliest VS spike onsets were at VS locations V2 and V5 left occipitotemporal, then 13 ms later spike peak at V1, left occipital, and 13–26 ms later two peaks in V3, right occipital. **(C)** Bipolar AP EEG and fifteen virtual sensor channels showing average of 129 interictal spike

based on spike peak detected at VS location V3 in right occipital lobe in the 20–70 Hz bandwidth. **(D)** Inset shows close-up of VS peaks at VS locations. The spike peak at V3, right posterior occipital occurred earliest followed by V4 (26 ms), right anterior occipital, and V6 (39 ms), right posterior temporal. **(E)** Subject's axial MRI scan showing right V3 and left V2 occipital earliest spike peak locations. **(F)** Subject's 3D MRI scan showing relative timing of peak locations starting in either left or right occipital lobe and spreading anteriorly. Red arrows indicate relative timing of peaks, pointing toward later peaks.


**FIGURE 7 | Subject 7**. **(A)** Bipolar AP EEG and 12 virtual sensor channels showing onset of generalized ictal onset in unaveraged signals filtered in bandwidth 20–70 Hz, digitized at 600 Hz. Red arrow indicates VS location V7 in right parietal. Ictal onset at V7 precedes by 120 ms the ictal onset in left parietal V6 (blue arrow). **(B)** Subject's MRI scan showing VS V7 location in right parietal lobe (at red arrow and

crosshairs). **(C)** Overview of ictal onset from right parietal lobe. **(D)** Second generalized ictal onset. Blue arrow indicates VS location V6 in left parietal Ictal onset at V6 precedes by 83 ms the ictal onset in right parietal V7 (red arrow). **(E)** Subject's MRI scan showing VS V6 location in left parietal lobe (at blue arrow and crosshairs). **(F)** Overview of ictal onset from left parietal lobe.

to the other has been quoted as 10 ms (Barth et al., 1982) for axonal transmission through the corpus callosum, the timing differences between component VS spikes of the EEG/MEG sensor spike were greater than 10 ms. The reason for the greater delay seen for our seven subjects is not certain. The spikes in the original simultaneous scalp EEG did not appear unusual, although the epileptiform discharge for subject 6 was actually a sharp wave by time duration. For these subjects the delay may relate to a time period required for one cortical region to activate another cortical region. Similar delays in propagation time of interictal spikes were noted previously in a subset of subjects in a combined scalp and intracranial EEG study; however, in those subjects also nothing strikingly unusual was noted in the morphology of the scalp EEG spikes and/or sharps waves (Alarcon et al., 1994). As noted in the results section, for some subjects, the VS locations that showed delayed spike peaks sometimes also showed a small amplitude deflection shortly after the VS spike at the VS location showing the earliest VS spike peak. Axonal transmission times might be appropriate from one intralobar, intrahemispheric, or interhemispheric location to another, but activation of independent synchronized activity at the secondary sites might require the longer time delays.

#### **COMPARISON OF PREOPERATIVE MEG BEAMFORMER PUTATIVE SOURCE TIMING AND INTRACRANIAL ECoG INTERICTAL/ICTAL RECORDINGS**

Although two of seven subjects did not proceed to Phase II intracranial recordings, the MEG beamformer and intracranial EEG (ECoG) did share similarities of at least onset locations for the remaining five subjects. Corroboration of patterns of spread, predicted by MEG, was seen to a certain extent in for subjects 1 and 2 in the ECoG.

Several differences of the two recording modalities became relevant for attempts to corroborate non-invasive MEG prediction of onset locations and spread patterns using intracranial EEG recordings. First is the physics limitation that MEG is more likely to record tangential sources in sulci, which ECoG may not detect well, while ECoG records radial sources on the surface of the brain, the gyral crowns, which MEG does not detect well (Hillebrand and Barnes, 2002). Thus some spread to deeper locations may not be detected by ECoG, while spread from one gyral crown to another, perhaps by *U*-fibers, may not be detected by MEG (e.g., **Figure 2**, subject 2 – the second ECoG-detected location in lateral inferior parietal cortex but not detected by MEG).

The limitations of ECoG to see deeper sources in sulci may in part be solved by combined SEEG and ECoG, or just SEEG, but considerable preoperative planning is required for either evaluation (Koessler et al., 2010; Kakisaka et al., 2012). Non-invasive preoperative MEG detection of multiple contributory sources relative to gyral crown components may improve with simultaneous higher resolution scalp EEG with more closely spaced electrode placement such as 64, 128, or 256 channel arrays (Yamazaki et al., 2012). Beamformer analyses of EEG scalp recordings have been reported, but success may be dependent on the accuracy of the forward model and knowledge of variations in skull impedance over the head (Steinstrater et al., 2010; Dang et al., 2011; Jon Mohamadi et al., 2012). Combined beamformer analysis with

MEG and EEG may provide improved timing details regarding activation of cortical sources with both radial and tangential current sources (Schoffelen and Gross, 2009; Shahbazi Avarvand et al., 2012).

Postoperative outcomes may suggest support for sources seen in MEG beamformer, but not detected with intracranial electrodes, when those MEG identified structures have not been resected. This difference may be the case for the insular and mesial parietal MEG sources for subject two that were not resected initially and seizures continued, but were resected with the subsequent hemispherectomy and seizures ceased. The difference could also be the case for a deep parietal source for subject four predicted by MEG that was not resected and seizures continued in that subject. Nonetheless, the persistence of seizures after MEG-predicted source locations that were not resected, or the disappearance of seizures after procedures like hemispherectomy, does not validate or invalidate precise location predictions of MEG.

Another limitation occurs when MEG detects bilateral sources, but subsequent invasive intracranial EEG is performed only for one hemisphere. The features that make MEG detection of interhemispheric spread clinically relevant are not yet delineated. For subjects 1 and 3, interhemispheric spread of interictal spikes without resection of the MEG-detected sources in the contralateral hemisphere were not relevant, as without their resection, the subjects still had a seizure free outcome 2 years postoperatively. By contrast, the detection of spread of interictal spikes to a site in the contralateral hemisphere in subject five, and detection of independent spikes that arose from that contralateral site and spread locally there, did presage the subject's postoperative seizures in that contralateral hemisphere. Numbers of intracranial electrodes that can be placed clinically safely are limited, so extensive bihemispheric ECoG or SEEG coverage will not likely be the solution. Better understanding of the clinical relevance of spread of spikes to a contralateral hemisphere may require comparing MEG studies with other non-invasive whole-head modalities such as the high density scalp EEG recordings or EEG-fMRI in addition to tracking surgical outcomes.

#### **MACRO NETWORKS IN CLINICAL EPILEPSY**

The concept of a single epileptic focus to be surgically removed was popular in the early years of epilepsy surgery. The concept may be very appropriate still for seizures limited to mesial temporal lobe. However, in pediatrics most medically intractable epilepsy is extratemporal in origin and may originate in more than one location, especially when cortical dysplasia is involved.

Networks in the propagation of seizures can occur at multiple levels (Lemieux et al., 2011). At a single cortical location, a local network over a few square centimeters of cortex may be required to sustain repetitive fast activity at the beginning of a seizure. Larger networks may involve a single or several contiguous gyri and may include both a SOZ and a surrounding region of rapid spread to give an EZ, which must be resected entirely to achieve seizure control (Engel, 1996; Wiebe et al., 2001; Luders et al., 2006). A still larger network may involve a SOZ with a region of rapid spread that may be located several gyri away or even in a different lobe or the contralateral hemisphere. For these larger networks there may be more than one cortical region with lower threshold for seizure onset and for ictal spread. Based on clinical experience with ECoG, when more than one cortical region is present with apparent low ictal threshold, the roles of SOZ and region of rapid spread may alternate between the two or more cortical regions, depending on the particular ictal event evolving at the time. For these latter cases, surgical treatment may involve resections or corticectomies at multiple sites, or sometimes the decision not to proceed with resective surgery.

The hope and expectation of clinicians treating pediatric medication resistant epilepsy is that intracranial electrodes can be placed during invasive monitoring to identify and delineate all the SOZs and all regions of rapid spread. Preoperative identification of multiple intracranial sources contributing to the interictal epileptiform discharge or ictal event is an important goal. The knowledge may help the neurosurgeon better plan the surgical approach to place intracranial electrodes (Knowlton et al., 2009).

Accomplishing this goal based on scalp EEG may be difficult using the conventional 10–20 electrode set and head surface topography of scalp voltage potentials. Efficacy may be improved with additional scalp electrodes, though multiple intracranial sources may be difficult to delineate with scalp recording topography, since EEG sees a composite of radial and tangential dipoles; source localization algorithms may be required (Ebersole, 1994; Ochi et al., 2001).

MEG studies usually require mathematical algorithms to predict intracranial source locations. The conventional algorithm is the single ECD, but in its simplest form it localizes a single source. When two or more sources are active simultaneously or the timing of their activation overlaps, the single ECD may predict a location between the two sources, which in some cases may mislocalize to white matter or ventricle. The moving ECD algorithm can identify that the average location of the composite of sources changes over time, but the algorithm is not configured to identify where individual contributing sources are located. Current source density algorithms may better account for multiple sources; however, in general the timing of activation of the different source may be somewhat difficult to ascertain based on the current density map (de Gooijer-van de Groep et al., 2012). All of the algorithms that do not use spatial filtering suffer from including fluctuations in background noise in the source localization calculations. Although keeping the subject quiet, relaxed, and still are important for reducing extracranial noise sources, intracranial brain noise fluctuations are ever present and are still part of the signal that is being localized (Ward et al., 1999).

On the other hand, every methodology has limitations. For clinical studies we also evaluate the traditional ECD model and both dipole scan algorithms (MUSIC) and current density algorithms (MNE, sLORETA). A limitation for beamformers is that two sources that changed together with identical relative amplitude time courses, but at different locations, could not be distinguished by a beamformer algorithm and would be mislocalized to a third location (Diwakar et al., 2011). The observation that beamformers do as well as they do, suggests that different brain regions most often do not have exactly identical relative amplitude time courses, particularly at the higher frequency bandwidths and if greater than 3 cm distant from one another (Belardinelli et al., 2012; Moiseev and Herdman, 2013).

In this study we first visually inspected the EEG/MEG signals for artifact, then used quantitative algorithms to (1) spatially filter the MEG signal, (2) measure excess kurtosis to identify candidate brain regions of spikiness, and (3) identify timing of the transient spikes in the regions that contributed to the spikiness of the signal. Then we returned to simple visual inspection of the VS signals to compare the timing of the original EEG/MEG interictal and ictal discharges with the VS spike timing at each of the VS locations. We found that multiple sources identified by this method were seen also in the ECoG in similar, but not exactly identical, cortical locations. At this time it is not certain whether the differences in locations between ECoG and beamformer were secondary to beamformer relative mislocalization, or whether the different cortical substrates measured (ECoG – gyral crowns; beamformer – gyral sulci), contributed to some of the differences.

Visual inspection alone was able to detect timing differences among the EEG/MEG spikes and the VS signals in 7 of 162 (4%) subjects, which is a low percentage of all subjects. However, visual inspection of signals should be done at each step of analysis/processing of the signals to obtain a Gestalt of what changes are occurring in the signals and for quality control to detect signal artifacts in what may become clinical tools. The range of differences in timing VS spikes for the remaining 155 subjects were in the range of 5–20 ms. The signal digitization rates we are using now (600–4 KHz) may allow us to use pursue quantitative comparisons after visual inspection, to evaluate better the timing differences in onset and peak across VS spikes and EEG/MEG sensor spikes.

There are several limitations to the current study. The sample size is a small percentage of all subjects we studied in the enrollment time window, as noted above. Only five of seven subjects had resective surgery, and for subject 5, a second ECoG recording was not done after the second MEG study. Only three of the five with resective surgery were seizure free; a fourth (subject 5) who did not have the contralateral source resected, continued to have seizures. Although that finding might be considered supportive, subject 5 did not have a second ECoG, so the contralateral source could not be confirmed as a second SOZ. For these reasons, the clinical significance of VS spike activation delay remains uncertain until a larger population can be examined and the degree of usefulness of the analysis can be better understood.

#### **CONCLUSION**

A spatial filter analysis of MEG signals such as a beamformer may be helpful during presurgical epilepsy evaluation to distinguish subregions of spiking during recorded interictal spikes and onsets of ictal discharges. In a subset of subjects the differences in timing between these onsets may be great enough to be apparent by visual inspection. For our subject subset the onset and timing differences roughly correlated with the SOZs and regions of spread determined by later ECoG. These subjects appeared to have a network of spread of the seizures within a cerebral lobe, across lobes and from one hemisphere to another. If the significance of VS spike activation delay can be better ascertained in a larger population, this kind of presurgical evaluation may be helpful for presurgical planning for surgical approach and intracranial electrode placement.

#### **REFERENCES**


Spatially filtered magnetoencephalography compared with electrocorticography to identify intrinsically epileptogenic focal cortical dysplasia. *Epilepsy Res.* 81, 228–232.


(2001). Dipole localization for identification of neuronal generators in independent neighboring interictal EEG spike foci. *Epilepsia* 42, 483–490.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 19 February 2013; accepted: 30 April 2013; published online: 14 May 2013.*

*Citation: Rose DF, Fujiwara H, Holland-Bouley K, Greiner HM, Arthur T and Mangano FT (2013) Focal peak activities in spread of interictal-ictal discharges in epilepsy with beamformer MEG: evidence for an epileptic network? Front. Neurol. 4:56. doi: 10.3389/fneur.2013.00056*

*This article was submitted to Frontiers in Epilepsy, a specialty of Frontiers in Neurology.*

*Copyright © 2013 Rose, Fujiwara, Holland-Bouley, Greiner, Arthur and Mangano. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*

## Dense array EEG source estimation in neocortical epilepsy

#### **MadokaYamazaki 1,2\*, Don M. Tucker 3,4, Marie Terrill <sup>4</sup> , Ayataka Fujimoto<sup>2</sup> and TakamichiYamamoto<sup>2</sup>**


#### **Edited by:**

Mark Holmes, University of Washington, USA

#### **Reviewed by:**

Mario A. Vanegas, Instituto Nacional de Neurologia y Neurocirugia, Mexico Erik K. St. Louis, Mayo Clinic and Foundation, USA

#### **\*Correspondence:**

Madoka Yamazaki, Department of Health Science, Daito Bunka University, 560 Iwadono, Higashimatsuyama, Saitama, 430-8558, Japan. e-mail: madokaymzk@gmail.com

**Rationale:** Dense array EEG (dEEG) evenly covers the whole head surface with over 100 channels contributing to more accurate electrical source imaging due to the higher spatial and temporal resolution. Several studies have shown the clinical utility of dEEG in presurgical clinical evaluation of epilepsy. However validation studies measuring the accuracy of dEEG source imaging are still needed. This can be achieved through simultaneously recording both scalp dEEG with intracranial electrodes (icEEG), which is considered as the true measure of cortical activity at the source. The purpose of this study is to evaluate the accuracy of 256-channel dEEG electrical source estimation for interictal spikes.

**Methods:** Four patients with medically refractory neocortical epilepsy, all surgical candidates, underwent subdural electrode implantation to determine ictal onset and define functional areas. One patient showed a lesion on the magnetic resonance imaging in the right parietal lobe. The patient underwent simultaneous recording of interictal spikes by both scalp 256-channelsvdEEG and icEEG. The dEEG was used to non-invasively estimate the source of the interictal spikes detected by the 256-channel dEEG array, which was then compared to the activity measured directly at the source by the icEEG.

**Results:** From the four patients, a total of 287 interictal spikes were measured with the icEEG. One hundred fifty-five of the 287 spikes (54%) were visually detected by the dEEG upon examination of the 256 channel head surface array. The spike amplitudes detected by the 256-channel dEEG correlated with icEEG spike amplitudes (p < 0.01). All spikes detected in dEEG were localized to the same lobe correctly.

**Conclusion:** Our study demonstrates that 256-channel dEEG can reliably detect interictal spikes and localize them with reasonable accuracy. Two hundred fifty-six-channel dEEG may be clinically useful in the presurgical workup for epilepsy and also reduce the need for invasive EEG evaluation.

**Keywords: dense array EEG, source estimation, neocortical epilepsy, interictal spike, intracranial EEG**

#### **INTRODUCTION**

For patients with intractable epilepsy, surgical therapy is an important treatment option. Although the propagation of seizures may involve complex cerebral networks, surgical resection of the cortical zone of seizure onset may be effective, if this zone can be identified and delineated from surrounding tissue. The first stage of diagnostic tests is typically non-invasive, with the goal being to characterize the epileptogenic zone as comprehensively as possible in order to guide the next stage of invasive testing with intracranial electrodes that confirms the decision on resective surgery.

Magnetoencephalography (MEG) and electroencephalography (EEG) have both been used as non-invasive methods in the electrophysiological evaluation of the epileptogenic zone. A number of studies has shown that multi-channel MEG is useful for dipole localization of interictal spike events, especially in neocortical epilepsy, where the sources may be more superficial (closer to the skull) compared to epilepsy with deeper sources, such as mesial temporal lobe epilepsy (Ricci et al., 1987; Gallen et al., 1995; Knowlton et al., 1997; Stefan et al., 2003; Huiskamp et al., 2010). MEG is inherently insensitive to deeper brain sources because magnetic signals fall off by the square of the distance from the source. Although EEG has better depth sensitivity than MEG, and it is sensitive to radial (gyral) as well as tangential (sulcal) sources, the electrical volume conduction of EEG is distorted by the resistive skull, requiring detailed computational modeling for an accurate inverse estimation (Ebersole, 2000; Vanrumste et al., 2000; Baumgartner, 2004; Plummer et al., 2007).

Although accurate electrical head modeling is now available for EEG, based on the patient's magnetic resonance imaging (MRI) and CT (for skull conductivity), the conventional 21 channel scalp EEG does not have high enough spatial resolution to allow for accurate source localization in most cases (Tucker et al., 2004; Holmes et al., 2005). Adequate source localization requires adequate spatial sampling in addition to sophisticated source analysis techniques. Freeman et al. (2003) reported that the spatial Nyquist of the human scalp EEG necessitates

<sup>1</sup> Department of Health Science, Daito Bunka University, Saitama, Japan

inter-sensor distances of <10 mm. Recently it has become possible to record dense array EEG (dEEG) of up to 256 channels, providing 20–25 mm interelectrode distances for most adult head sizes. Several studies have shown that interictal source localization with dEEG suggests its usefulness in neurosurgical planning in epilepsy (Lantz et al., 2003a,b; Michel et al., 2004a; Holmes, 2008). However validation studies measuring the accuracy of dEEG source imaging are still needed. This can be achieved through simultaneously recording the scalp dEEG with intracranial electrodes (icEEG), which provides a local measure of electrical activity at the gyral cortical (pial) surface. In previous studies we have investigated the sensitivity and accuracy of dEEG source analysis comparing 256-channel dEEG source analysis of interictal activity and simultaneously recorded icEEG in patients with mesial temporal lobe epilepsy (Yamazaki et al., 2012a,b). In the present study, we extend this approach to investigate the accuracy of source estimation with dEEG in patients with neocortical epilepsy using simultaneous dEEG and icEEG recording.

#### **MATERIALS AND METHODS**

#### **PATIENTS**

We studied four patients, each of whom had suffered from intractable localization related epilepsy for at least 2 years. Each patient was a surgical candidate and underwent a presurgical workup including conventional long-term EEG monitoring, MRI, 125-Iomazenil (IMZ)-single-photon emission tomography (SPECT), 18-F fluorodeoxyglucose (FDG)-Positron emission tomography (PET), and neuropsychological testing. Clinical information for these patients is summarized in **Table 1**. Case 1 had right parietal cortical dysplasia.

We received approval for this study from Seirei Hamamatsu General Hospital Human Subject Committee and informed consent was obtained from all patients.

#### **INTRACRANIAL EEG RECORDING**

Subdural strip and grid electrodes (Unique medical, Tokyo, Japan) were implanted for each patient in order to delineate the epileptogenic zone for cortical excision, and to separate it from functional areas. All electrode contacts were platinum, and interelectrode distance was 10 mm. The location and types of subdural electrodes used for each patient are shown in **Table 1**.

#### **DENSE ARRAY EEG RECORDING**

The dEEG was recorded with the 256-channel Geodesic Sensor Net (Electrical Geodesics, Inc., Eugene, OR, USA), providing 257 electrodes (including vertex reference) covering the face and neck as well as cranium, with 20–25 mm interelectrode distances. The coverage of the face and neck is important for measuring the downward projection of electrical potentials from basal brain regions.

#### **SIMULTANEOUS icEEG AND dEEG ACQUISITION**

We conducted simultaneous icEEG recording with NicoletOne (CareFusion, Middleton, WI, USA) and dEEG with Net Amp300 (Electrical Geodesics, Inc., Eugene, OR, USA) at 1 kHz sampling with bandpass filter 0.1 and 400 Hz for approximately 30–40 min. A digital pulse from the icEEG system was provided to the dEEG acquisition system for synchronization. Prior to the simultaneous icEEG-dEEG recording, each patient recovered from the intracranial placement for at least 3 days, allowing surgical wounds to heal to avoid infection risk during simultaneous scalp recording. There were no complications due to any of the simultaneous recordings.

#### **DATA ANALYSIS**

Frequent epileptiform discharges were selected during artifact free periods in the simultaneous dEEG and icEEG recordings. Interictal spikes seen in icEEG were visually identified by a certified

#### **Table 1 | Case summary.**


SPS, simple partial seizure; CPS, complex partial seizure; sGTC, secondary generalized tonic clonic seizure; \* , hypoperfusion area; \*\*, hypometabolism area; mT, mesial temporal lobe; laT, lateral temporal lobe.

EEG technologist (MY) and confirmed by a board-certified clinical epileptologist (AF). As expected, many of the icEEG spikes were not detectible by inspection of the dEEG. We calculated the spike detection rate of 256-channel dEEG (with the icEEG detection as the criterion), and we measured the average maximum amplitude for both the icEEG detected spikes and the smaller subset of dEEG detected spikes.

Electrical source localization was conducted at the rising phase for the dEEG detected spikes with a linear inverse method (LAURA) using the GeoSource 1.0 software package<sup>1</sup> within the space of a 3D head model derived from the Montreal Neurological Institute's average adult MRI. The LAURA constraint provides results very similar to the low-resolution electromagnetic tomography (LORETA) spatial Laplacian constraint (Pascual-Marqui et al., 2002), and this smoothing constraint has been shown to provide stable source estimation of interictal epileptiform events in neurosurgical planning for epilepsy (Lantz et al., 2003a,b; Michel et al., 2004a,b).

#### **RESULTS**

#### **SPIKE DETECTION RATE**

A total of 287 spikes were recorded in icEEG with four patients during recordings of approximately 30–40 min each. One hundred

<sup>1</sup>http://www.egi.com

fifty-five of these spikes (54%) were also clearly distinguishable from background activity in the simultaneously recorded dEEG. The dEEG detection rate for the spikes verified by icEEG for each patient is summarized in **Table 2**. Case 2 showed two independent patterns of interictal spikes, one originating from the left frontal region and the other from the left mesial temporal region. The detection rate of neocortical spikes was 56% and that for mesial temporal spikes in this case was 39%.

#### **icEEG AMPLITUDE OF DETECTABLE SPIKES**

The average maximal icEEG amplitude of neocortical interictal spikes that were also detectable by the scalp dEEG was 968µV (standard deviation of 139µV), significantly higher (*p* < 0.01) than that of dEEG undetectable spikes (757µV; standard deviation of 102µV). The average maximal icEEG amplitude of mesial temporal spikes that were also detectable by the scalp dEEG was 1112µV (standard deviation of 197µV), significantly higher (*p* < 0.01) than that of dEEG undetectable spikes (763µV; standard deviation of 172µV).

#### **DENSE ARRAY EEG SOURCE ESTIMATION**

**Figure 1** shows dEEG source estimation for the dEEG detected icEEG spikes in each patient, visualized by displaying the voxel with maximal source amplitude (as well as voxels with similar

#### **Table 2 | Spike detection rate.**



amplitudes) using a section from the MNI typical brain image<sup>2</sup> consistent with the MNI atlas used for the cortical dipole locations in the GeoSource software. All spikes detected in dEEG were localized to the same lobe correctly and 141 of 155 spikes (91%) were well localized, close to the position confirmed by subdural electrodes.

**Figure 2** shows more detail for a typical example of icEEG and dEEG source estimation. For this selected spike, the icEEG and conventional EEG recordings are shown in **Figure 2A**. The locations of the icEEG strips are shown in **Figure 2B**. The head surface (scalp, face, and neck) distribution of the dEEG potentials is shown in **Figure 2C**. From inspection, this dEEG head surface topography is consistent with a mostly radial source in the right posterior midline, with diffuse inversions over both sides of the face indicating the approximate radial source orientation. The peak of the non-invasive dEEG source estimation for this spike is shown in **Figure 2D** consistent with visual inspection of the surface array, and the post-operative MRI shows the resected region in **Figure 2E**.

#### **SURGICAL OUTCOME**

All patients underwent resective surgery based on the intracranial icEEG findings, functional mapping, and imaging studies. Pathological findings are case 1: FCD type IIB, case 2: FCD type IIB, case 3: FCD type IA, and case 4: no abnormality. Surgical outcome is Engel class Ia (case 3, 4), IId (case 1), and IIIa (case 2) within the post-operative follow-up period 23–35 months at this time (mean, 30 months).

#### **DISCUSSION**

Simultaneous recording of non-invasive dEEG with icEEG allowed a direct comparison of the sensitivity of dEEG to epileptiform events (spikes) that were confirmed by intracranial recordings. In these four patients with neocortical (extra temporal) epilepsy, the non-invasive dEEG localization of spikes predicted not only the intracranial location of the spikes, but also the cortical location of seizure onset (**Figure 1**).

Clearly there were many spikes visible in the icEEG that were not detected through visual inspection of the dEEG. For the neocortical spikes identified by icEEG in the present study, the dEEG spike detection rate was 56%. As would be expected, the dEEG detected spikes were consistently larger in amplitude than the dEEG undetected spikes in each patient. For the mesial temporal spikes examined in our previous studies with simultaneous icEEG and dEEG (Yamazaki et al., 2012a,b), the dEEG detection rates were somewhat lower, 45 and 42%, respectively. Also as expected, to be detectable by visual inspection of the dEEG, the mesial temporal lobe spikes in those studies were typically larger in amplitude (averaging 1236µV; Yamazaki et al., 2012b) compared with the dEEG detected neocortical spikes of the present study (averaging 968µV). For the one patient of the present study with mesial temporal spikes, the 39% dEEG detection rate was similar to the previous rates for detecting mesial temporal spikes.

A higher non-invasive detection sensitivity to neocortical than mesial temporal spikes has also been observed in those MEG studies where validation has also been obtained from simultaneous icEEG recording. Compared with dEEG, MEG is particularly insensitive to the deep sources presented by mesial temporal lobe spikes. Huiskamp et al., 2010 reported that whole head MEG detected only 28% of the icEEG spikes generated in the mesial temporal lobe, whereas it detected 70% of the icEEG detected spikes from the lateral temporal lobe and extra temporal regions. Similarly, Oishi et al., 2002 observed that MEG detected only 26% of the mesial temporal spikes, but 53% of lateral frontal spikes that were confirmed by simultaneous icEEG recording. Considering our previous results and the present study, 256-channel dEEG is roughly equivalent to whole head MEG in detecting neocortical sites but superior to MEG in detecting mesial temporal lobe spikes.

Magnetic fields are less distorted by the skull, cerebral fluid, and scalp than electrical fields (Nakasato et al., 1994; Ebersole, 1999; Minassian et al., 1999; Otsubo et al., 1999; Morioka et al., 2000). However, the accuracy of MEG for deep sources, such as in the mesial temporal lobe, appears limited because magnetic signals fall off by the square of the distance (Mikuni et al., 1997; Oishi et al., 2002; Rampp and Stefan, 2007). In addition, the 256 channel dEEG sensor net includes sensors on the face and neck that improve the characterization of the electrical fields from the temporal lobe that are directed inferiorly and anteriorly, where even whole head MEG has limited coverage. Finally, the inclusion of accurate head conductivity information, such as from the MRI and CT atlas used for finite difference conductivity modeling in GeoSource software in the present analyses, is now helping to deal with the ambiguities of electrical volume conduction that plague EEG source localization.

Presurgical evaluation for what is typically called neocortical epilepsy, located outside the mesial temporal lobe (and archicortical hippocampus), requires precise accuracy because the epileptogenic zones are often located adjacent to critical functional areas of cortex (Otsubo et al., 1997; Minassian et al., 1999; Chitoku et al., 2001). For the four patients of the present study, the noninvasive dEEG source localization of the patient's typical interictal spikes was predictive of the icEEG localization of not only the interictal spikes but also the seizure onset. Although the interictal events do not always localize to the same region as seizure onset, the utility of careful non-invasive localization of interictal events as a first stage in the presurgical workup is consistent with previous dEEG studies on electric source imaging in epilepsy (Michel et al., 2004b; Brodbeck et al., 2010). Nonetheless, noninvasive localization of seizure onset is an important confirmation for the hypotheses generated from localizing interictal events, and may be an important form of converging evidence in preparing for icEEG mapping of seizure onset (Holmes et al., 2008, 2010). Non-invasive 256-channel dEEG is now available for long-term monitoring for seizure onset; although more challenging than conventional EEG monitoring, it is more practical than attempting seizure monitoring with MEG.

There were several important technical limitations of this study that should be addressed in future research. We used standard 3D sensor positions, rather than the precise 3D locations of the 256-channels for each individual net application, such as can be

<sup>2</sup>www.bic.mni.mcgill.ca

#### **FIGURE 2 | A typical example of icEEG and dEEG source estimation (case 1)**. **(A)** The icEEG (upper) shows a right parietal spike. Interictal spike is

shown at electrodes # 1–5–7 and 12–16 which are located over the parietal region. The EEG (lower) is simultaneous recorded 256-channel dEEG with 19-channel of 10/20 display. **(B)** Placement of subdural electrodes and the location of the interictal spike. N Indicates the electrodes which show the

interictal spike. **(C)** Two hundred fifty-six-channel dEEG topographic plot of the corresponding spike. The view is looking down on top of the head with nose at the top. The 256-channnel topographic plot was instructive in localizing the spike to the right parietal lobe. **(D)** The source estimation by dEEG superimposed on a standard MRI. The interictal spike is localized to right parietal head region. **(E)** Post-operative MRI.

obtained with photogrammetry (Russell et al., 2005). Even with careful placement of the net, the sensor positions can vary in relation to the cortex from patient to patient depending on head shape and the tissues of the face and neck. In addition, we used the GeoSource atlas finite difference model of head conductivity, based on the MNI brain template (and skull CT fit to that template), rather than building an electrical conductivity model based on the individual's MRI and CT for the source estimation. Brodbeck et al., 2011 examined the sensitivity and specificity of dEEG source estimation comparing individual MRI and standard MRI template and showed the benefit of using individual MRI. For example, Case 1 in the present series had a focal cortical dysplasia in the right parietal region; rather than using the GeoSource atlas it would have been useful if we could compare

#### **REFERENCES**


Holmes, M. D., Brown, M., and Tucker, D. M. (2005). Dense array EEG and source analysis reveal spatiotemporal dynamics of epileptiform discharges. *Epilepsia* 46(Suppl. 8), 136.


the source localization of the spikes with the lesion location precisely by building the conductivity model from the patient's MR and CT.

#### **CONCLUSION**

Simultaneous recording of dEEG and icEEG allows direct examination of the accuracy of localizing cortical electrical sources from non-invasive dEEG recordings. In four patients with neocortical epilepsy examined with 256-channel dEEG, we observed good prediction of both the interictal spikes and seizure onset as verified by icEEG. Careful non-invasive analysis of the patient's typical interictal events with dEEG may be an important first step in the neurosurgical planning for resection of the seizure onset zone.

clinical yield and localization precision. *J. Clin. Neurophysiol.* 21, 71–83.


adolescent with right frontocentral epilepsy. *Epilepsia* 40, 608–613.


detection and localization. *Epilepsy Res.* 98, 166–173.

Yamazaki, M., Terrill, M., Fujimoto, A., Yamamoto, T., and Tucker, D. M. (2012b). Integrating dense array EEG in the presurgical evaluation of temporal lobe epilepsy.*ISRN Neurol.* 2012, 924081.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 17 February 2013; paper pending published: 09 March 2013; accepted: 15 April 2013; published online: 13 May 2013.*

*Citation: Yamazaki M, Tucker DM, Terrill M, Fujimoto A and Yamamoto T (2013) Dense array EEG source estimation in neocortical epilepsy. Front. Neurol. 4:42. doi: 10.3389/fneur.2013.00042*

*This article was submitted to Frontiers in Epilepsy, a specialty of Frontiers in Neurology.*

*Copyright © 2013 Yamazaki, Tucker, Terrill, Fujimoto and Yamamoto. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*

## Methods for examining electrophysiological coherence in epileptic networks

#### **Jasmine Song<sup>1</sup>\*, Don M. Tucker 1,2,Tara Gilbert <sup>1</sup> , Jidong Hou<sup>1</sup> , Chelsea Mattson<sup>1</sup> , Phan Luu1,2 and Mark D. Holmes <sup>3</sup>**

<sup>1</sup> Electrical Geodesics, Inc., Eugene, OR, USA

<sup>2</sup> Department of Psychology, University of Oregon, Eugene, OR, USA

<sup>3</sup> Department of Neurology, Regional Epilepsy Center, University of Washington, Seattle, WA, USA

#### **Edited by:**

Don Tucker, Electrical Geodesics, Inc., USA; University of Oregon, USA

#### **Reviewed by:**

Andreas Schulze-Bonhage, University Medical Center Freiburg, Germany Pierre-Pascal Lenck-Santini, Geisel School of Medicine at Dartmouth, USA

#### **\*Correspondence:**

Jasmine Song, Electrical Geodesics, Inc., 1600 Millrace Dr. Suite 200, Eugene, OR 97403, USA. e-mail: jsong@egi.com

Epilepsy may reflect a focal abnormality of cerebral tissue, but the generation of seizures typically involves propagation of abnormal activity through cerebral networks. We examined epileptiform discharges (spikes) with dense array electroencephalography (dEEG) in five patients to search for the possible engagement of pathological networks. Source analysis was conducted with individual electrical head models for each patient, including sensor position measurement for registration with MRI with geodesic photogrammetry; tissue segmentation and skull conductivity modeling with an atlas skull warped to each patient's MRI; cortical surface extraction and tessellation into 1 cm<sup>2</sup> equivalent dipole patches; inverse source estimation with either minimum norm or cortical surface Laplacian constraints; and spectral coherence computed among equivalent dipoles aggregated within Brodmann areas with 1 Hz resolution from 1 to 70 Hz. These analyses revealed characteristic source coherence patterns in each patient during the pre-spike, spike, and post-spike intervals. For one patient with both spikes and seizure onset localized to a single temporal lobe, we observed a cluster of apparently abnormal coherences over the involved temporal lobe. For the other patients, there were apparently characteristic coherence patterns associated with the discharges, and in some cases these appeared to reflect abnormal temporal lobe synchronization, but the coherence patterns for these patients were not easily related to an unequivocal epileptogenic zone. In contrast, simple localization of the site of onset of the spike discharge, and/or the site of onset of the seizure, with non-invasive 256 dEEG was useful in predicting the characteristic site of seizure onset for those cases that were verified by intracranial EEG and/or by surgical outcome.

**Keywords: epilepsy, spike, coherence, networks, dEEG, cortical surface Laplacian**

#### **1. INTRODUCTION**

There is increasing evidence in the recent literature that functional networks of the human brain can be identified through correlation analysis of fluctuations in hemodynamic functional magnetic resonance imaging (fMRI) measures (Fox et al., 2005; Fair et al., 2007, 2008;Van Dijk et al., 2010). To capture the millisecond dynamics of the electrophysiological abnormalities of epilepsy, it is important to apply a similar approach to electrophysiological network analysis of epileptiform discharges, including both spikes and seizures (Gotman et al., 2006; Hamandi et al., 2006; Laufs et al., 2006; Gotman, 2008). Seizures may originate in a focal site of pathological tissue, such as a malformation or lesion of the cortex, but the clinically significant seizure typically involves propagation through some functional networks of the brain.

In a recent study from our laboratory (Ramon et al., 2008; Ramon and Holmes, 2012), non-linear dynamic measures of local correlation among the 256 dense array electroencephalography (dEEG) channels showed abnormally high levels of synchronization over cortical regions that later prove to be seizure onset zones. Because these analyses were conducted with the head surface (scalp) dEEG, they are less precise than synchronization analysis performed with cortical source analysis, assuming, of course, that the cortical source analysis is indeed accurate. In the present study, we developed computational models of the head geometry and conductivity for each patient, including extraction of the cortical surface and tessellation with oriented source dipoles, in order to improve the electrical source analysis of the dEEG measures of epileptiform events. The goal was to apply the synchronization analysis (spectral coherence) to the waveforms of the cortical sources directly.

#### **2. MATERIALS AND METHODS**

The workflow of the data analysis procedure employed for the present study is presented in **Figure 1**.

#### **2.1. PATIENTS AND CLINICAL SETTING**

Patients in this study are five individuals with medically refractory epilepsy who were referred to the University of Washington Regional Epilepsy Center for evaluation and treatment. All were considered as potential candidates for epilepsy surgery. At the time of the evaluation, the subjects ranged in age from 11 to 28 years (mean age = 21 years). Four of the five were males. Duration of

**FIGURE 1 | EEG was continuously recorded**. Spikes were detected. The continuous EEG was segmented. The head volume (MR/CT) was segmented and registered to the electrodes. The lead field matrix was calculated. The cortical source waveforms were estimated. The source waveforms were averaged over Brodmann Areas. The BA waveforms were transformed to the frequency domain by fast Fourier transformation. The coherences of the BA waveforms were computed. If the high performance computing (HPC) power is available, the step of averaging over BAs is not required.

epilepsy at the time of the assessment varied from 2 to 16 years (mean = 7 years). Each individual had failed trials of least four standard antiseizure medications, either as monotherapy or in combination, at the time of the presurgical evaluation (range 4–7 drugs).

All patients underwent as part of the presurgical assessment a comprehensive history and clinical examination, routine awakesleep electroencephalogram (EEG), high-resolution magnetic resonance imaging (MRI), formal neuropsychological testing, and standard EEG-video monitoring. In cases 1, 3, and 4, after a thorough review of these initial studies, the consensus of opinion rendered by epileptologists at Epilepsy Center was that intracranial, subdural strip, or grid electrode EEG-video monitoring was indicated, as the initial non-invasive evaluation failed to provide adequate information regarding seizure onsets. For case 2, intracranial recording was not done because the parents decided against surgery. For case 5, the evidence for temporal lobe onset was sufficiently clear from non-invasive recording to proceed to temporal lobe resection. Prior to invasive monitoring, all subjects underwent 256 dEEG video monitoring to record electrographic activity with a greater degree of spatial resolution and enhance the non-invasive estimate of ictal onsets.

#### **2.2. INTERICTAL SPIKE IDENTIFICATION**

Dense array EEG data were collected during long term monitoring on the epilepsy unit. Identification of spikes in the continuous EEG file was first performed automatically using Persyst (Persyst, San Diego, CA, USA). Briefly, Persyst employs a neural network technique to identify spikes and a statistical clustering algorithm to group spikes according to spatial similarity. The identified spikes were then reviewed by an epileptologist to confirm that identified spikes were indeed spikes and that their grouping were appropriate (i.e., that spikes were topographically homogeneous within each group).

Persyst identified multiple spike types for all patients. However, based on review by the epileptologist, in three patients, only one spike type was confirmed. In the other two patients, the epileptologist confirmed two spike types (see **Table 1**). For these two patients, we analyzed the data for from both groups but only report results from the analysis of the spike type consistent with their icEEG or resected zone (see Discussion). Once identified in the continuous record, the EEG was segmented centered on the spike peak, with 1.5 s before and after the spike for further analysis (see below).

#### **2.3. REQUIREMENTS FOR ELECTRICAL SOURCE ANALYSIS**

The understanding of the neural sources of the EEG has been improved in recent years by a number of factors. A first step has been more adequate spatial sampling of the potentials at the head surface with dense sensor arrays (Tucker, 1993). The estimation of the cortical sources of the surface EEG requires accurate characterization of tissue geometries, the EEG sensors positions relative to the tissues, and electrical conductivity of the tissues. Together, the combined information makes up the *electrical head model*.

To create the electrical head model, the EEG sensors must first be located precisely in relation to the head surface. Derivation of sensor positions can be accomplished by using photogrammetric techniques (Russell et al., 2005). Next, the cortical sources must be characterized in the electrical model, typically as electrical dipoles. If the orientation of each patch of cortex is known, such as from cortical surface extraction from the MRI (Dale and Sereno, 1993), then oriented sources may be used. Otherwise, xyz "triple" dipole models must be used for each cortical source, with a considerable loss of constraint, and thus loss of precision for the inverse estimation. Next, a mathematical model of the volume conduction of cortical sources to the head surface is constructed (Malmivuo et al., 1997). When a high-resolution, volumetric MRI can be segmented accurately, a finite difference (or finite element) model can be constructed, specifying the conductivity of each voxel of brain, cerebral spinal fluid, skull, and scalp (Salman et al., 2005a; Turovets et al., 2006). Spherical or boundary element models cannot account for the complexity of the conductive compartments, particularly for inferior head regions that are critical for the propagation of discharges basal brain structures (basal temporal and orbital frontal) to the head surface (including the face and neck). Because the resistive skull must be specified precisely, fitting an volumetric x-ray computed tomography (CT) image to the head model is a critical step (Salman et al., 2005b).

Finally, given the complex *forward model* specified by these facts of tissue geometry and electrical conductivity, the estimation of cortical source activity is made through an *inverse model*. When the parameters of the electrical head model are specified loosely or inaccurately, such as with electrodes only assumed to be at typical positions, the cortical sources at unknown orientations and therefore modeled as dipole "xyz" triples, and the


#### **Table 1 | Clinical features of the 5 patients.**

L, left; R, right; \*, typical spike type; temp, temporal; fron, frontal; mid, midline; occi, occipital; infe, inferior.

head conductivity modeled loosely as concentric spheres, the constraints on the source estimation are so loose as to create large uncertainty bounds. Fundamentally, the inverse model is *ill-posed*, and source estimation is an approximation in the best case. The ill-posed inverse model may be stabilized through regularizing with various methods, including the minimum norm (Ou et al., 2008), 3D Laplacian smoothing (Pascual-Marqui et al., 1994), statistical standardization of the projection of each source to the sensors (Pascual-Marqui, 2002), ECD (equivalent current dipole) (Fischer et al., 2005), or beamforming (Gross et al., 2001). In previous research on localization of spike and seizure onset, reasonable results have been created with dEEG (128 and 256 channel) measurement, using a realistic head conductivity atlas (finite difference model from an atlas MRI) with cortical triples and various methods of regularization and smoothing of the inverse (Michel et al., 2004; Tucker et al., 2007, 2009; Guggisberg et al., 2008; Brodbeck et al., 2010; Holmes et al., 2010a,b; Holmes, 2011; Bouet et al., 2012; Yamazaki et al., 2012). These results suggest that more accurate solutions for the inverse problem can be achieved with a more accurate volume conductor head model.

In the present study, we improved the accuracy of the electrical head model for each patient by specifying cortical source dipoles whose orientation was known from extracting the gyral-sulcal cortical surface for each patient. The goal was to develop an accurate source analysis for major regions of the cortex (Brodmann areas) that could then support analysis of electrophysiological coherence that may be related to the pathological activity in cerebral networks engaged by the patient's epileptic discharges.

#### **2.4. CONSTRUCTING THE ELECTRICAL HEAD MODEL**

The electrical head model was constructed with the BrainK software (Li et al., 2006). For each patient and each EEG recording session, sensor positions of the 256-channel Geodesic Sensor Net were determined with multi-camera geodesic photogrammetry system (GPS) (Russell et al., 2005). The data were then registered with the patient's head model, derived from structural MRIs. The patient's volumetric MRI was segmented into scalp, skull, CSF, and brain gray and white matter. Accurate segmentation was optimized with a unique *relative thresholding* algorithm, which corrects for the inhomogeneity of MR images. Characterization of the detailed geometry of the skull was accomplished through non-linear warping of an atlas skull CT to the patient's MRI image. Analysis of the electrical source activity within the patient's MRI allowed inspection of the convergence of anatomical features with the epileptic discharges. The MRI for Patient 3, for example, showed a blurred gray-white boundary in the left dorsal frontal region, possibly suggesting a malformation, and this structural feature proved to be co-located with the patient's typical source-localized epileptiform discharges.

With the cortical ribbon segmented by its inner and outer table, the outer (gray to CSF) cortical surface was meshed, and then cast into a geodesic graph theoretical representation with the Chaco graph algorithms (Sandia National Laboratories). The graph representation allows flexible tessellation of the cortical surface, and flexible representation of mathematical constraints defined on that surface, such as the Cortical Surface Laplacian (CSL) defined below. For the present analysis, 1 cm<sup>2</sup> cortical patches (about 1200 per hemisphere) were created. For each patch, the average orientation within that patch was characterized by computing the surface normal for each triangle in that patch, then computing the vector sum across all patch triangles as the representation of the equivalent dipole orientation for that patch.

#### **2.5. INVERSE METHOD**

The use of oriented dipoles on the cortical surface provides an important constraint on the inverse source estimation process that appears to improve accuracy considerably. Much of the ambiguity in the relation of cortical electrical activity to the head surface measurements is created by the extensive gyral and sulcal folding of the individual patient's cortex. The requirement, of course, is that the position of the EEG sensors in relation to the oriented sources must be specified precisely; otherwise the misalignment will result in incorrect source attributions.

Given accurate geometric alignment, accuracy in the conductivity of the volume conduction model is then required, including the precise geometry of the skull, such as from CT. Given these several measurement constraints on the forward model, we have observed that minimal regularization constraints on the inverse estimation, such as the minimum norm, yield accurate results. Support for this accuracy has come from initial (unpublished) validation studies with dEEG mapping of individual sensory and motor potentials, with similar results as seen with fMRI for that individual. These conclusions on the decreased importance of regularization are consistent with previous reports using oriented cortical sources (Dale and Sereno, 1993; Ou et al., 2008; Knosche et al., 2013).

#### **FIGURE 2 | Patient 1's non-invasive dEEG localization of a typical**

**(left-frontal) spike**. **(A)** Topographic map of scalp potentials at peak of spike (which is used for spike classification). Orientation is top looking down with the nose at the top; red = positive, blue = negative, white = zero. The right inferior (face and neck) spatial gradient (described by the isopotential lines) suggests a right anterior temporal source, even though the spike (i.e., negative features) was located at left-frontal recording sites. **(B)** Waveforms of EEG from all channels (top row), source waveforms from all 2240 sources (second row), source waveform from dipole that shows maximal activity at spike peak (third row, blue trace) and source waveform from dipole that shows maximal activity at 50% of spike peak (third row magenta trace), and

the frequency spectra of the blue source waveforms (third row) shown in forth row. Note that the dipole location of the blue trace is in the left temporal lobe and the magenta trace in the third row is in the right anterior temporal lobe, consistent with the icEEG data and resection zone. In **(B)**, each column represents 1 s of data from a 3-s segment. Middle column represent data centered on the spike peak. **(C)** Electrical source localization with the patient's individual electrical head model. Left: scalp potential reconstructed from the source estimate (at 50% of spike peak). Right: partially inflated brain illustrating the dominant source associated with the data illustrated in left figure (location of dipole described by magenta trace in **(B)**, third row). Data was thresholded to only show top 5% activity.

Although reasonably accurate with the orientation constraints, the minimum norm remains a non-specific constraint, and is not an optimal fit to the physiological properties of electrophysiological activity of the cortex. To optimize the inverse constraints that are appropriate to our knowledge of localized cortical activity, we have developed an inverse regularization or smoothing constraint that makes use of the knowledge of adjacency of dipoles on the extracted 2-dimensional cortical surface. This is the *Cortical Surface Laplacian*, computed as the Laplacian operator for a given dipole source against the (typically 6) neighboring dipoles on the 2D cortical surface. There are both physical and physiological rationales for this constraint. Physically, only synchronous activity of adjacent, laminar cortical neurons (thought to be pyramidal cells and their apical dendrites) is capable of generating the far field observed by head surface electrodes (Nunez and Srinivasan, 2005). Physiologically, the columnar and hyper-columnar organization of the human cortex (Jones, 2009) makes it plausible that synchronous activation – perhaps even as large as a 1-cm<sup>2</sup> patch – might be reasonable.

The CSL constraint is then used to smooth or regularize the inverse estimation, to avoid inappropriate mathematical results that are possible with the instabilities of the matrix inverse. The inverse estimation then: (1) minimizes the difference between the forward oriented cortical source model and the observed surface potentials (the *data fidelity term*) and (2) minimizes the CSL as the *regularization term*. With Φ as the measured surface potentials, *K* as the lead field or forward volume conduction matrix, *J* as the cortical sources, and ε as the error term, the surface potentials can be modeled as the projection of the cortical source activity through the head volume model, plus error:

$$
\Delta = Kl + v.
$$

The inverse formulation minimizes the data fidelity term plus the regularization term:

$$\hat{J} = \arg\min\_{J} \left\{ \left\| \Phi - K \right\|^2 + \alpha \left\| W \right\|^2 \right\},$$

where *W* comprises a discrete spatial Laplacian operator defined by

$$\mathcal{W}\_{\vec{ij}} = \begin{cases} -1 & \text{if } i = j \\ \frac{1}{N\_i} & \text{if } i \neq j \text{ and voxel } j \text{ is a neighbor of voxel } i, \\ 0 & \text{otherwise} \end{cases}$$

where *N<sup>i</sup>* is the number of cortical surface dipole neighbors (patches or voxels) to dipole *i*.

#### **2.6. INVESTIGATING PATHOLOGIES OF CORTICAL COHERENCE**

The goal of the present study was to implement and then evaluate these several advances in individual head modeling for dEEG source waveforms to investigate cortical network patterns (source coherence) in relation to epileptiform discharges (spikes) in patients being examined for possible neurosurgical resection of the epileptic focus. The hypothesis was that abnormal patterns of cortical source coherence during interictal events may reveal abnormal patterns of electrophysiological synchronization associated with the seizure onset zone. For each of the five patients, we selected characteristic interictal spikes and then clustered these to insure that all spikes in the cluster had a common head surface topography in the 256 surface dEEG array. The cortical source waveforms were computed with the patient's individual head model with the CSL constraint for three intervals: 1 s before the spike; 1 s centered on the spike, and 1 s after the spike.

Coherence is the standardized cross-spectral density computed between two signals, *x* and *y*:

$$\begin{aligned} &X, Y: \text{Fast Fourier Transform of } x \text{ and } y\\ &S\_{\text{XX}}(f): \text{autospectrum of Xat frequency } f\\ &S\_{YY}(f): \text{autospectrum of Yat frequency } f\\ &S\_{XY}(f): \text{cross - spectrum of X and Yat frequency } f\\ &Coherence(f) = \frac{S\_{XY}(f)}{\sqrt{S\_{\text{XX}}(f) \cdot S\_{YY}(f)}}.\end{aligned}$$

Just as the power spectrum describes frequency spectrum of the variance of one signal, the cross-spectrum describes the frequency spectrum of the covariance between two signals. Because the covariance is affected by the amplitude of the component signals, the coherence measure is standardized, by dividing the cross-spectrum by the square root of the product of the power spectra. Coherence is thus the frequency domain analog of the Pearson correlation coefficient (which is covariance of two variables divided by the square root of the product of their respective variances).

With the 1-s epoch for each interval (pre-spike, spike, and post-spike), the frequency resolution was 1 Hz, and coherences for 1–70 Hz were examined. The frequency range was divided into five frequency bands; delta, theta, alpha, beta, and gamma. With ∼2400 patches of size 1 cm<sup>2</sup> for the cortical surface, and thus 2400 equivalent source dipoles, there are 2,878,800 comparisons [(*N* <sup>2</sup> − *N*)*/* 2]. With 5 coherence values for each comparison, the coherence matrix for each 1-s interval becomes somewhat large (∼360 million entries). For this preliminary study, we therefore grouped the dipoles within the BAs for each hemisphere, and averaged all source waveforms in each BA. For the grouping, the individual's brain was transformed to Talairach space, and the BAs for each cortical patch were determined with the Talairach Daemon (Lancaster et al., 2000). This resulted in about 40 BAs per hemisphere (depending on the anatomy of the patient's cortex and its alignment with the Talairach transform). This resulted in ∼3160 cross-signal coherences in each matrix. For an overview visualization of the coherence matrices, each was plotted as an 80 × 80 matrix, with the coherence entries scaled by a color palette (**Figure 3**). Each matrix shows coherence values in the form of a color palette, with green to blue as lower coherence and yellow to red as higher coherence. The variables on the axes are the approximate Brodmann areas (slightly different areas are reflected by these numbers for each patient, depending on the patient's cortical anatomy and alignment with the Talairach atlas). In each matrix, the intra-hemispheric coherences for the right hemisphere are in the upper left, and intra-hemispheric coherences for the left hemisphere in the lower right, and inter-hemisphere coherences in the lower left and upper right quadrants. In general, for

all patients, intra-hemispheric coherences are higher than interhemispheric coherences. For a more anatomical visualization, the cortical patch dipoles were plotted in their 3D positions in a glass brain (the vertical or axial projection was made transparent), and the highest 100 coherences from each interval were represented by a line between the coherent signals (e.g., **Figure 4**).

#### **3. RESULTS**

For each patient, we first summarize the standard clinical evaluation: (1) non-invasive spike and seizure localization with conventional (Ten-Twenty System) EEG; (2) non-invasive spike and seizure localization with 256 dEEG; (3) icEEG when that was done; (4) the decision whether the patient was a surgical candidate; and (5) if so, the surgical outcome (Engel class) at the present time in **Table 1**. Next, the coherence analyses for the interictal

events (spikes) are presented, for the pre-spike, spike, and postspike intervals. One typical spike cluster (comprising spikes with similar head surface topographies) is described for each patient; the set of spike clusters obtained for that patient and differences in coherence patterns for the other clusters is summarized in each case. The peak coherences are shown for all 3 intervals (pre-spike, spike, and post-spike). In each case, the differences across differing spike clusters for that patient, and the differences between pre-spike, spike, and post-spike intervals are summarized.

#### **3.1. PATIENT 1**

EEG: right temporal seizures, not well localized on the scalp. dEEG: right medial temporal localization for both spikes and seizures.

icEEG: right medial temporal localization for both spikes and seizures.

Surgery: right temporal resection.

Outcome: seizure free (Engel Class I) for 1 year.

Coherence analysis: cortical patch tessellation resulted in 2240 patch dipoles for Patient 1 which were classified into 82 Brodmann areas (41 per hemisphere). For a typical cluster of leftfrontal spikes (*N* = 6), the coherence matrices are shown in **Figure 3**.

For this patient, the spike clustering show a group with a left-frontal distribution. However, analysis of the spike at 50% maximum amplitude (as recommended in previous work by Lantz et al., 2003) show a dominant right-lateralized source in the temporal region (see **Figure 2**). An overview summary of Patient 1's spike coherences is given in **Figure 3**. The pre-spike interval coherences appear unremarkable, with the possible exception of certain low inter-hemispheric values (blue) for the gamma band. The spike-interval coherences show higher values overall, including in the inter-hemispheric quadrants for the lower (delta, theta, alpha) bands. The post-spike interval coherences return to somewhat lower values, with perhaps a specific decrease in the inter-hemispheric values for the lower (delta, theta, alpha) bands. The peak coherences for this same spike cluster for Patient 1 are shown in **Figure 4**. The patterns for peak coherences for Patient 1 are similar for pre-spike, spike, and post-spike intervals, with a strong delta cluster (and perhaps theta) over the right temporal region. There seems to be a similar pattern across bands in which there are fairly distributed peak coherences across the left hemisphere, but a focal cluster in the temporal area for the right hemisphere. Somewhat higher values of the distributed left hemisphere coherences are observed for the alpha band in the pre-spike interval. The right temporal cluster may be particularly strong for the delta band in the spike interval. The post-spike interval includes what appears to be stronger occipital coherences, including inter-hemispheric, particularly for the theta band.

#### **3.2. PATIENT 2**

EEG: interictal events (spikes) were right sided. Seizures seemed to be on the right, but the pattern was unclear.

dEEG: right temporal localization for both spikes and seizure onset in **Figure 5**.

icEEG: not done.

Surgery: none. Parents decided not to pursue surgical resection. Coherence analysis: cortical patch tessellation resulted in 2272 patch dipoles for Patient 2 for 83 Brodmann Areas. A cluster of 15 right frontal spikes was identified. In the dEEG analysis (both atlas and individual cortical source localization), the onset of these spikes was localized to the right temporal area, and there was rapid propagation to sources in the right frontal area (creating the right frontal topography for the spike peak).

Patient 2's peak coherences were generally similar for the prespike, spike, and post-spike intervals, as was observed for all patients in this series (**Figure 6**). For the delta, theta, and alpha bands particularly, there was a cluster of peak coherence over the left-frontal and temporal regions. For beta and gamma bands, there was a cluster of peak coherence over the right posterior regions.

#### **3.3. PATIENT 3**

EEG: seizures observed with apparent onset in the left-frontal midline.

dEEG: similar left-frontal midline localization for both spikes and early seizure activity, but evidence of early onset for both spikes and seizures in left temporal region in **Figure 7**.

icEEG: seizure activity observed over left-frontal midline cortex. Surgery: resection of left-frontal midline.

Outcome: after 1 year, seizures are improved (Engel Class II).

Coherence analysis: cortical patch tessellation resulted in 2260 patch dipoles for Patient 3 for 80 Brodmann Areas. Spike clustering showed the predominance of spikes showed a left-frontal topography (*N* = 76); these were examined with coherence analysis.

The peak coherence analysis for Patient 3 (**Figure 8**) showed fairly similar patterns across the three intervals. Coherences are somewhat higher for right hemisphere cortical sources than left hemisphere cortical sources. Somewhat higher coherences are seen for the beta band than for other frequencies, and for the gamma band for the right hemisphere intra-hemispheric values particularly. The pattern seems like a mirror image of that for Patient 1, with a distributed frontotemporal cluster of peak coherences over the right hemisphere, but a more focal pattern of peak coherences over the left temporal region.

#### **3.4. PATIENT 4**

EEG: seizures appeared to involve midline structures.

dEEG: seizures appeared to involve right central and midline structures. Interictal analysis of several spike clusters (right inferior, left inferior) showed the onset of each of these differing spike discharges appeared to be in the left and right medial temporal cortex as well as bilateral mediofrontal cortex (**Figure 9**). icEEG: seizures involved both left and right midline cortex. Surgery: none. Patient was judged not to be a surgery candidate. Coherence analysis: cortical patch tessellation resulted in 2247 patch dipoles for Patient 4 for 81 Brodmann Areas. Spike clustering showed the three differing peak topographies: midline frontal, left inferior, and right inferior. The right inferior spike cluster was selected for coherence analysis (*N* = 12).

The peak coherences of the spike interval for Patient 4 were possibly more focal than for the other intervals, implicating the left inferior region for the delta band (**Figure 10**). However, as is typical in these analyses, the peak coherence patterns were roughly similar for the three (pre-spike, spike, and post-spike) intervals. For Patient 4, there were suggestions

of a tight temporal lobe cluster, in this case on the left, which seemed to engage other left posterior sites for the delta.

#### **3.5. PATIENT 5**

EEG: seizures appeared to involve right frontotemporal areas, and spikes were right sided.

**FIGURE 7 | Patient 3's non-invasive dEEG localization of a typical (left-frontal) spike**. **(A)** Topographic map of scalp potentials at peak of spike (which is used for spike classification). Orientation is top looking down with the nose at the top; red = positive, blue = negative, white = zero. **(B)** Waveforms of EEG from all channels (top row), source waveforms from all 2260 sources (second row), source waveform from dipole that shows maximal activity at 50% of spike peak (third row, blue trace) and source waveform from dipole at a distance from maximal location (third row magenta trace), and the frequency spectra of the blue source waveforms (third row) shown in forth row. Note that the dipole location of the blue trace is in the left-frontal lobe. In **(B)**, each column represents 1 s of data from a 3-s segment. Middle column represent data centered on the spike peak. **(C)** Electrical source localization with the patient's individual electrical head model. Left: scalp potential reconstructed from the source estimate (at 50% of spike peak). Right: partially inflated brain illustrating the dominant source associated with the data illustrated in left figure (location of dipole described by blue trace in **(B)**, third row). Data was thresholded to only show top 5% activity.

**frequency bands**.

**FIGURE 9 | Patient 4's non-invasive dEEG localization of a typical (right inferior) spike**. **(A)**Topographic map of scalp potentials at peak of spike (which is used for spike classification). Orientation is top looking down with the nose at the top; red = positive, blue = negative, white = zero. **(B)** Waveforms of EEG from all channels (top row), source waveforms from all 2247 sources (second row), source waveform from dipole that shows maximal activity at 50% of spike peak (third row, blue trace) and source waveform from dipole at a distance from maximal location (third row magenta trace), and the

frequency spectra of the blue source waveforms (third row) shown in forth row. In **(B)**, each column represents 1 s of data from a 3-s segment. Middle column represent data centered on the spike peak. **(C)** Electrical source localization with the patient's individual electrical head model. Left: scalp potential reconstructed from the source estimate (at 50% of spike peak). Right: partially inflated brain illustrating the dominant source associated with the data illustrated in left figure (location of dipole described by blue trace in **(B)**, third row). Data was thresholded to only show top 5% activity.

**five frequency bands**.

dEEG: both spikes and seizures showed a right medial temporal lobe onset in **Figure 11**.

icEEG: not done.

Surgery: right temporal lobe resection. Patient is seizure free (Engel Class I) for about 1 year.

Coherence analysis: cortical patch tessellation resulted in 2244 patch dipoles for Patient 5 for 82 Brodmann Areas. Spike clustering showed only one topography: right inferior spikes (*N* = 56).

Patient 5's peak coherences (**Figure 12**) show a unique pattern of left-frontal and temporal coherences, together with high coherences between these regions with a source in the right frontal area. This pattern is observed for the lower frequency bands. Whereas features of this pattern can be seen for the beta and gamma bands, there are increased coherences over right hemisphere regions, particularly for the frontal and temporal regions for the gamma band.

#### **4. DISCUSSION**

The results of the coherence analyses in the present study must be seen as preliminary. We are only beginning to validate the oriented cortical source analysis constructed from individual head models. The present data are the first we have examined with oriented cortical source analysis of epileptic discharges and the present results are perhaps most useful in illustrating the approach to characterizing cortical network electrophysiology rather than providing definitive evidence on epileptic networks.

In a similar vein, DICS (dynamic imaging of coherent sources) has been used to localize the epileptic spikes using oscillatory features (Guggisberg et al., 2008; Bouet et al., 2012). The DICS method focuses on the source localizations of the spike-related high frequency activity (>20 Hz). The present study considered low as well as high frequencies, using coherence analysis to attempt to discern pathological networks contributing to the epileptic spikes.

In considering the accuracy of the oriented cortical source analysis with the CSL constraint, we observed good correspondence in the present CSL analyses of spike and seizure onset and the atlas-based model (dipole triples with the MNI average cortex and a generic head conductivity model) implemented in Geosource 2.0. Given the previous validation of the atlas (Geosource 2.0) analysis in relation to intracranial recordings (Yamazaki et al., 2012) and surgical outcome (Holmes et al., 2008, 2010a) the present results imply that the CSL inverse is at least roughly correct when applied to the oriented cortical sources constructed for the individual patient.

Given the previous findings of characteristic patterns of EEG synchronization over the seizure onset zone in the scalp (head surface) EEG (Ramon and Holmes, 2012), we organized the present study to examine cortical source coherence, rather than EEG electrode channel coherence. On first principles, we would expect the cortical source coherence – even with the summation over the large Brodmann areas – to reveal more about cortical function than the badly superposed gyral and sulcal fields that are summed at any given surface electrode.

Whether the coherence analysis of the source activity before, during, and after the spike will yield clinically meaningful information remains to be seen. It is noted here that, because this study represent the first systematic study of the approach we employ here, we also analyzed several aspects of the data but do not report findings in detail here. First, for those patients with confirmed multiple spike types (e.g., Patients 1 and 4), the coherence patterns were not substantially different between the spike types at the three time intervals. Second, we also examined random data from spikefree intervals, and the coherence patterns for each patient was similar to their spike-interval coherence patterns. Nevertheless, the overall results do point to interesting possibilities.

Patient 1 did show a promising result. He manifested a single seizure focus in the right temporal lobe that, even though the peak of the spike indicated a left-frontal focus, when resected, resulted in his being seizure free. His characteristic interictal spikes were associated with a tight pattern of right temporal lobe peak coherences that dominated over other coherences in the right hemisphere. Of course, an obvious interest is whether seizure onset is associated with a similar pattern of abnormal synchronization. We chose to focus on interictal events for this preliminary study because the dEEG data quality is high, and there are many spikes that provide for statistical stability of the results. Nonetheless, it is clearly important to extend these methods to analysis of cortical synchronization associated with seizure onset.

Patient 5 also showed spikes and seizures localized to the right temporal lobe, and a right temporal lobectomy eliminated seizures (at least for the 1-year follow up to date). Yet Patient 5 showed a dominance of peak coherences over*left* frontotemporal regions for most frequency bands. It was the case that the gamma band showed the strongest coherences for Patient 5, where a right-lateralization of frontal and temporal coherences was observed. Nonetheless, there was no simple association of temporal lobe epileptic pathology that could be concluded from these two cases which appeared to be straightforward single sided temporal lobe epilepsy from the dEEG evaluation, and which have been seizure-free following temporal lobectomy.

Patient 2's parents decided not to pursue surgical intervention. Yet the dEEG localization of both spikes and seizure onset pointed to a right temporal localization. With no indications of other epileptic foci, the dEEG localization has generally proven accurate for predicting both icEEG results and surgical outcome (Holmes et al., 2008, 2010b; Holmes, 2011). Nonetheless, Patient 2 showed a consistent cluster of left-frontal and temporal peak coherences, for the delta and theta bands, a pattern that seems discrepant with the more standard dEEG localization of spike and seizure source amplitudes.

Patient 3 was a difficult case, in that the primary findings pointed to the left superior medial frontal area, and there was a possible cortical anomaly in this area, and yet surgical resection of this area was only effective in improving, rather than eliminating, seizures. The largest discharges for both spikes and seizure onset were clearly localized to the left superior medial frontal area for Patient 3. There were subtle clues in the dEEG amplitude analysis that these discharges (for both spikes and the large discharge at seizure onset) were preceded by a source amplitude increase in the left temporal area. However, temporal lobe source activity is often

**FIGURE 11 | Patient 5's non-invasive dEEG localization of a typical (right inferior) spike**. **(A)**Topographic map of scalp potentials at peak of spike (which is used for spike classification). Orientation is top looking down with the nose at the top; red = positive, blue = negative, white = zero. **(B)** Waveforms of EEG from all channels (top row), source waveforms from all 2244 sources (second row), source waveform from dipole that shows maximal activity at 50% of spike peak (third row, blue trace) and source waveform from dipole at a distance from maximal location (third row magenta trace), and the

frequency spectra of the blue source waveforms (third row) shown in forth row. In **(B)**, each column represents 1 s of data from a 3-s segment. Middle column represent data centered on the spike peak. **(C)** Electrical source localization with the patient's individual electrical head model. Left: scalp potential reconstructed from the source estimate (at 50% of spike peak). Right: partially inflated brain illustrating the dominant source associated with the data illustrated in left figure (location of dipole described by blue trace in **(B)**, third row). Data was thresholded to only show top 5% activity.

associated with the early milliseconds of spike discharges in our experience, but it is still unclear whether this has clinical significance. It is worth noting that Patient 3 showed a tight left temporal coherence cluster, in contrast to more distributed coherences over the right hemisphere (**Figure 8**), in a mirror image of the pattern observed for Patient 1. Whether or not this restriction of hemispheric peak coherences to the temporal lobe could be taken as a signal of temporal lobe involvement in seizure onset is unclear without further studies.

Patient 4 is another example where early source amplitude (rather than coherence) changes are observed in the temporal lobe at the onset of spikes, but with uncertain clinical significance. In examining each of Patient 4's several spike types (left inferior, right inferior), a left medial temporal source seemed to be the first indicator of the discharge. However, seizures appeared to show a midline onset, not only in the conventional EEG but in the dEEG as well. Furthermore, bilateral icEEG recordings showed apparent seizure onsets from both left and right midlines. The peak coherence patterns during spikes for Patient 4 showed left temporal and posterior clusters for delta and theta bands.

Although such clues are intriguing, the challenge for a novel method like cortical source coherence analysis will be to obtain validation from a number of converging perspectives, including normal and functional studies as well as epilepsy studies. We think it is impressive that the delta band of coherence shows apparently meaningful results in these analyses, and to a lesser extent theta band. Whereas the transient features of the spike amplitude involve spectral components in the delta and theta frequency bands, there is typically little power in the gamma band in the head surface EEG. Furthermore, gamma frequency of the EEG is often contaminated by muscle artifact. Nonetheless, although it is clearly important to consider artifactual explanations for the apparent EEG gamma patterns in high frequency analyses, it appears at the present time that meaningful coherence results can be obtained in this high frequency range in human subjects. The observation of high frequency oscillations in intracranial EEG at seizure onset (Kobayashi et al., 2010; Jacobs et al., 2011) makes it clear that high frequency features are an important target for non-invasive EEG analysis, and future research should examine the high gamma band (above 70 Hz).

A clear limitation of the present study is that we only examined EEG coherence around spike events rather than around

#### **REFERENCES**


reconstruction: a linear approach. *J. Cogn. Neurosci.* 5, 162–176.


seizure onset. Although interictal events are often localizing for seizure onset, there are very likely dynamics of cortical synchronization that can only be studied in relation to seizure onset.

Another limitation of the present study is that we grouped the oriented cortical sources into Brodmann areas, primarily for computational convenience. A goal for research on oriented cortical sources must be to implement the high performance computing necessary to examine the full cortical source coherence matrix, applying not only statistical analysis and decomposition but also effective visualization methods. Recent advances in diffusion imaging tractography (Scherrer and Warfield, 2010) have made it possible to compute a full set of likely cerebral fiber tracts (*N* = 15 million), and to align them with both subcortical structures and the cortical patch tessellation implemented in the present study. Anatomically specific measures of cortical electrophysiology will continue to bring important tools to the study of epilepsy in the near future. The challenge will be to apply these tools to the study of epileptic networks in order to understand their clinical significance.

#### **CONCLUSION**

EEG coherence provides a measure of the covariance among signals that resolves the frequency features of that covariance. In the present research, coherence was computed among cortical source waveforms, localized to the oriented cortical surface with detailed head models for individual patients. The improved source analysis with anatomical constraints promises new insights into the network properties that may be altered in epilepsy. The first step in these analyses was localizing the amplitude onset of the spike with 256 dEEG, which proved useful in predicting the onset of the seizure in each case. The utility of dEEG spike localization was verified by dEEG and icEEG analyses of seizure onset, and was confirmed in several cases by success of surgical resection. Whether cortical source coherence analysis adds to the clinical utility of spike amplitude localization remains to be seen. Perhaps the most intriguing observation was the frequent pattern of strong coherence centered on temporal lobe structures in several patients. Understanding the clinical significance of this pattern will require further studies of seizure onset, as well as contrast analyses with normal individuals and with apparently normal EEG intervals in the epileptic patients.

magnetoencephalography results to obtain favourable outcomes in epilepsy surgery. *Brain* 128, 153–157.


fMRI: a multimodal tool for epilepsy research.*J. Magn. Reson. Imaging* 23, 906–920.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 05 February 2013; accepted: 30 April 2013; published online: 15 May 2013.*

*Citation: Song J, Tucker DM, Gilbert T, Hou J, Mattson C, Luu P and Holmes MD (2013) Methods for examining electrophysiological coherence in epileptic networks. Front. Neurol. 4:55. doi: 10.3389/fneur.2013.00055*

*This article was submitted to Frontiers in Epilepsy, a specialty of Frontiers in Neurology.*

*Copyright © 2013 Song , Tucker, Gilbert , Hou, Mattson, Luu and Holmes. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*

### Stochastic behavior of phase synchronization index and cross-frequency couplings in epileptogenic zones during interictal periods measured with scalp dEEG

#### **Ceon Ramon1,2\* and Mark D. Holmes <sup>3</sup>**

<sup>1</sup> Department of Electrical Engineering, University of Washington, Seattle, WA, USA

<sup>2</sup> Department of Bioengineering, Reykjavik University, Reykjavik, Iceland

<sup>3</sup> Department of Neurology, University of Washington, Seattle, WA, USA

#### **Edited by:**

Don Tucker, Electrical Geodesics, Inc., USA; University of Oregon, USA

#### **Reviewed by:**

Jose F. Tellez-Zenteno, University of Saskatchewan, Canada Silvia Kochen, University of Buenos Aires, Argentina

#### **\*Correspondence:**

Ceon Ramon, Department of Electrical Engineering, University of Washington, Campus Box 352500, Seattle, WA 98195, USA. e-mail: ceon@u.washington.edu

These results in part were presented as a poster at the 66th Annual Meeting of American Epilepsy Society, 30 November–04 December 2012, San Diego, CA.

The stochastic behavior of the phase synchronization index (SI) and cross-frequency couplings on different days during a hospital stay of three epileptic patients was studied for non-invasive localization of the epileptogenic areas from high density, 256-channel, scalp EEG (dEEG) recordings. The study was performed with short-duration (0–180 s), seizure-free, epileptiform-free, and spike-free interictal dEEG data on different days of three subjects. The seizure areas were localized with subdural recordings with an 8 × 8 macro-electrode grid array and strip electrodes.The study was performed in theta (3–7 Hz), alpha (7–12 Hz), beta (12–30 Hz), and low gamma (30–50 Hz) bands. A detrended fluctuation analysis was used to find the long range temporal correlations in the SI that reveals the stochastic behavior of the SI in a given time period. The phase synchronization was computed after taking Hilbert transform of the EEG data. Contour plots were constructed with 20 s time-frames using a montage of the layout of 256 electrode positions. It was found that the stochastic behavior of the SI was higher in epileptogenic areas and in nearby areas on different days for each subject.The low gamma band was found to be the best to localize the epileptic sites. Also, a stable higher pattern of SI emerged after 60–120 s in the epileptogenic areas. The cross-frequency couplings of SI in theta–gamma, beta–gamma, and alpha–gamma bands were decreased and spatial patterns were fragmented in epileptogenic areas. Combinations of an increase in the stochastic behavior of the SI and decrease in cross-frequency couplings are potential markers to assist in localizing epileptogenic areas. These findings suggest that it is possible to localize the epileptogenic areas non-invasively from a short-duration (∼180 s), seizure-free and spike-free interictal scalp dEEG recordings.

**Keywords: epilepsy localization, dEEG, phase synchronization, stochastic behavior of EEG, cross-frequency couplings**

#### **INTRODUCTION**

In previous studies (Ramon and Holmes, 2012, 2013), we have demonstrated that the stochastic behavior of the phase synchronization index (SI) derived from the high density (256-channel) scalp EEG data has a potential to localize the epileptic sites in human subjects. A 256-channel scalp EEG data is often referred as dEEG or hdEEG data. These previous studies were limited either to a single day (Ramon and Holmes, 2012) or for multiple days but in a single EEG band of low gamma frequencies (Ramon and Holmes, 2013). The present work is an extension and a comprehensive summary of the previous work on multiple days in four EEG bands, viz., theta (3–7 Hz), alpha (7–12 Hz), beta (12–30 Hz), and low gamma (30–50 Hz) and it confirms our previous results that the stochastic behavior of the SI in the low gamma band was best suited to localize the epileptic sites derived from a short-duration (∼180 s) interictal scalp dEEG data.

The EEG phase synchronization plays an important role in studying the network connectivity of different regions of the brain under normal and diseased states, including various forms of epileptic seizures (Varela et al., 2001; Baier et al., 2012; Lang et al., 2012; Palmigiano et al., 2012). Local and long rage connectivity can be studied with EEG phase synchronizations (Varela et al., 2001). An increase in the phase synchronization has been observed in epileptogenic zones (Schevon et al., 2007; Warren et al., 2010) and, also, the connectivity patterns are different in epileptogenic zones as compared with other cortical zones not involved in epileptic activities (Varotto et al., 2012). Phase synchronization also influences the genesis of epileptic activity. From invasive recordings it has been observed and suggested that there is an increase in the phase synchronization activity in the epileptogenic regions of the brain and these regions are functionally isolated from the surrounding regions of the brain (Warren et al., 2010). It has also been observed with magnetoencephalogram (MEG) recordings in epileptic patients that there were fluctuations in synchrony between neighboring cortical networks (Dominguez et al., 2005).

The phase synchronization also randomly fluctuates with time. The stochastic behavior of these random fluctuations can be studied with detrended fluctuation analysis (DFA) which is, often, used to study long range temporal correlation (LRTC) in a time series data, such as, EEG data (Peng et al., 1995; Hardstone et al., 2012). It has been shown by Linkenkaer-Hansen et al. (2005) that DFA exposes LRTCs that are characteristic of epileptogenic neocortical networks, the areas where epilepsy begins. We examined the stochastic behavior of the SI on different days derived from a shortduration (∼3 min) interictal dEEG data. Our results show that the stochastic behavior of the SI is higher in the vicinity of the epileptogenic zones and possibly maybe useful to localize the epileptic sites. Here, the stochastic behavior of the SI or LRTC of SI is used in a synonymous fashion.

The brain also exhibits oscillatory activity in various frequency bands. The cross-frequency coupling, where one band modulates the activity of a different band is a very powerful tool to study the oscillatory activity of the brain. It has been observed with subdural recordings that the power of the high gamma (50–100 Hz) band is phase locked to theta oscillations (Canolty et al., 2006). The crossfrequency related phase synchrony among different bands in the frequency range of 3–80 Hz has also been observed in human MEG data (Palva et al., 2005) and also in scalp EEG data (Friese et al., 2012). We also found that during object naming tasks measured with dEEG data there was some coupling in the theta–gamma band (Ramon et al., 2009). This was observed by plotting the difference in phase synchronization between two EEG bands over a time period of 3 s. A similar analysis is used here to examine stochastic behavior of cross-frequency couplings. We found that cross-frequency couplings decreased in theta–gamma, beta– gamma, and in alpha–gamma bands. Also, some complex spatial patterns were observed. These results can also be used as an additional marker to non-invasively localize the seizure onset areas from interictal scalp EEG data.

#### **MATERIALS AND METHODS EEG DATA OF PATIENTS**

Our procedures for data collection and analysis have been described previously (Ramon et al., 2008; Ramon and Holmes, 2012, 2013). Only a brief summary is given here. Epileptic seizure areas in patients were localized with intracranial subdural EEG (ECoG) recordings with 8 × 8 contact grid electrodes and also with strip electrodes. The electrodes on the grid had an exposed surface area defined by 2.3 mm diameter and with a center-to-center, inter-electrode separation of 1.0 cm (Johnson et al., 2012). The strips had the same size electrodes with the same inter-electrode separation. Prior to this during pre-surgical evaluations, high density 256-channel scalp EEG data was collected continuously for 7–12 days. The data was collected with an EEG system developed by Electrical Geodesics, Inc. (Eugene, OR, USA). The electrode caps were filled with a conducting gel with an effective diameter of ∼1.0 cm. For an adult head, from the center of one electrode to the other, the inter-electrode separation was ∼2.0 cm or less (Tucker, 1993). The data was collected with a sampling rate of 250 Hz, i.e., the time difference between each sample was 4 ms. We used data of three adult subjects. All were candidates for resection surgery and after surgery they were cured of the epileptic seizures. Subjects were not on any medication during dEEG and ECoG monitoring. All data were collected at the Regional Epilepsy Center, University of Washington under the authorized human subjects protocol.

For each subject, we selected data on three randomly selected days for analysis. For a given day, approximately, 10 min long, seizure-free and spike-free data from each patient during sleep was selected for analysis. The selected data sets were not in close proximity to seizures. Out of this, a continuous 3 min long data was randomly selected and imported into MATLAB for further analysis. The analysis was repeated for other two randomly selected days of the data.

#### **COMPUTATIONS OF PHASE SYNCHRONIZATION INDEX**

The raw EEG data was normalized to their common averaged reference and then filtered using a FIR bandpass filter in the appropriate EEG band. Excessively noisy channels were eliminated by replacing them with the averages of their neighbors. For each subject, there was only one noisy channel. The 60 Hz power line artifact was eliminated with a matching pursuit filter (Gratkowski et al., 2006).

The synchronization between a pair of channels was inferred from a statistical tendency to maintain a nearly constant phase difference over a given period of time even though the analytic phase of each channel may change markedly during that time-frame (Freeman and Rogers, 2002). The Hilbert transform was applied on the pairs of EEG traces with a stepping window which is long enough to encompass at least two cycles of the lowest frequency in a given band. For example for the low gamma (30–50 Hz) band, the lowest frequency will be 30 Hz and based on that we selected a window of 80 ms. The analysis was repeated by stepping the window at 8 ms intervals, i.e., two digitized points. Similarly, the size of the stepping window was selected for other EEG bands while the step size of 8 ms was kept the same. The window sizes selected for theta, alpha and beta bands were, 680 ms, 300 ms, and 170 ms, respectively. The SI was computed for each pair of EEG traces.

The phase of the analytic signal has a sawtooth pattern which is unwrapped to produce a cumulative linear phase of the signal. The phase difference between the two channels was computed by subtracting the phase of one channel from the other. This phase difference was then used to determine the SI. The mathematical techniques for computing synchronization indices are given in detail elsewhere (Tass et al., 1998; Freeman and Rogers, 2002). We computed synchronization indices based on Shanon entropy function (Tass et al., 1998). Phase locking, i.e., synchronization between the phases of two signals within a stepping window was given by Shanon entropy function, *e*(*t*), defined as:

$$e\left(t\right) = -\sum\_{i=1}^{N} p\_i \ln p\_i \tag{1}$$

where *p<sup>i</sup>* was the relative frequency of finding the phase difference modulus of 2π in the *i*th bin. The function *e*(*t*) varied between zero and its maximum value of *e*max = ln *N*. We used 100 bins (*N* = 100) for the phase difference in a given stepping window. This phase locking, *q(t),* was normalized and is represented as:

$$q\left(t\right) = \frac{e\_{\text{max}} - e(t)}{e\_{\text{max}}} \tag{2}$$

**FIGURE 1 | The analysis of phase synchronization between two EEG signals**. **(A)** A short segment of two EEG traces, **(B)** unwrapped phases, ϕ<sup>1</sup> and ϕ<sup>2</sup> after taking Hilbert transform of EEG signals, **(C)** time derivative,

dϕ/dt, of the phase difference, ϕ = ϕ1−ϕ2, of two EEG signals, and **(D)** phase synchronization index, q(t), between the two EEG signals. Notice that q(t) is high when dϕ/dt is nearly zero.

**FIGURE 2 | (Left) each electrode is paired with six nearby electrodes for computation of phase synchronization indices**. There are 21 possible combinations. Few of them are shown by arrows. (Middle) average of phase synchronization indices over 21 combinations of electrode pairs. The phase

synchronization index of the white noise is shown with a dashed line. (Right) montage layout of 256 electrodes above the head for plotting purposes. The nose is on the top. The ellipse encloses the area covered by electrodes above the eye level but excludes the electrodes on the forehead.

The *q*(*t*) has a value of zero for uniform distribution of phase differences and a value of one for a spike or delta distribution of phase differences between two signals. This *q(t)*is also called phase syncronization index (SI).

The **Figure 1** shows a small segment of EEG traces for two nearby channels. It also shows the unwrapped phases,ϕ<sup>1</sup> and ϕ2, of the two signals, the time derivative of the phase difference, *d*(ϕ)/*dt*, and the *q*(*t*). The phase difference is defined as: ϕ = ϕ1−ϕ2. There

are two regions marked with light-green color shade where *d*ϕ/*dt* is nearly zero and the *q*(*t*) is high. This shows that in the shaded regions, the two EEG signals are almost in phase synchronization and due to this value of *q*(*t*) is high.

A global SI was also computed for each electrode by pairing it with the nearby six electrodes. There were 21 combinations of electrode pairs for each given electrode. The *q*(*t*) was averaged over these electrode pairs for each given electrode. The averaged *q*(*t*) represents the local cortical connectivity approximately over a circular area defined by a 4.0 cm diameter which is based on inter-electrode separation of 2.0 cm in the geodesic net used for collecting dEEG data (Tucker, 1993). **Figure 2** shows the electrode pairs, the global SI and a montage layout of 256 electrode positions for plotting purposes. The nose is on the top, back of the neck is at the bottom, left of the subject is on the left side of the plot and right side of the subject is on the right side of the plot. The ellipse, roughly covers the electrodes above the eye level including occipital areas, but, excludes the electrodes on the forehead. The horizontal and vertical axes for plots are in normalized length units.

Only few possible electrode pairs, out of possible 21 pairs are shown by the double headed arrows in **Figure 2**. The averaged *q*(*t*) given in **Figure 2** is slightly less than the *q*(*t*) given in **Figure 1** for a single electrode pair. In general, there is always some spatial decoherence involved between two electrode pairs (Freeman and Rogers, 2002). Thus, it is anticipated that by averaging over 21 electrode pairs, the averaged SI will be less than the SI of an individual electrode pair. However, it is still above the synchronization level of the white noise (Tass et al., 1998). Following the procedures given in Tass et al. (1998), a Gaussian white noise signal was generated, filtered and the SI was computed. The mean value of the SI of the white noise was found to be 0.18 with 95% confidence

interval. The noise level is marked as a dashed line in the middle plot in **Figure 2**. The averaged *q*(*t*) in **Figure 2** shows sustained levels of phase synchrony which is interrupted by dips below the noise level. Hereafter, the averaged *q*(*t*), for simplicity, is called the phase syncronization index (SI).

#### **COMPUTATIONS OF STOCHASTIC BEHAVIOR**

We used DFA to compute the stochastic behavior of the SI. The cumulative sum of each channel was calculated. This sum was divided into windows of 10 s length, i.e., windows were of 10, 20, 30, . . ., 170 and 180 s length. Within each window, a linear fit was found and the cumulative sum was detrended. Next, the root-mean-squared (RMS) fluctuation of this detrended sum was calculated. The median fluctuation at each window size was taken. The log of this median fluctuation was plotted against the log of the window size, and a linear fit was found. The slope of this linear fit, denoted as, α, is the result of the DFA which can be expressed as (Hardstone et al., 2012):

$$F\left(L\right) \propto L^{\alpha} \tag{3}$$

where *F*(*L*) is a fluctuation function and *L* is the window length. This is also called LRTC (Peng et al., 1995) of a given time series data, such as, EEG or SI which is derived from EEG. For simplicity, we will refer LRTC of SI as the stochastic behavior of the SI. The stochastic behavior of the SI, i.e., α, was computed for each channel as explained above. Color intensity plots of α were constructed using a montage layout of 256 electrode positions given in **Figure 2**.

The exponent, α, in Eq. 3 is also related to the power law representation of random walk processes after a walk of length *L* (Hardstone et al., 2012). This exponent, α, represents many properties of a time series data. If 0 < α < 0.5, the slope is negative and the time series sequence is anticorrelated or negatively correlated. If α ' 0.5, the time series sequence is an uncorrelated white noise. If 0.5 < α < 1, the time series sequence has LRTCs. If α ' 1.0, the time series sequence is a 1/*f* pink noise, where *f* is the frequency of the signal. If α > 1.0, the time series sequence is non-stationary, random walk like and has strong correlations which are not of a power law form. A special case is for α ' 3/2 which represents that the time series sequence is a Brownian noise. In general, the EEG signal is non-stationary while the white noise is stationary.

#### **RESULTS**

#### **SPATIOTEMPORAL PLOT**

For the subject #1,the spatiotemporal plots of the stochastic behavior of the SI are given in **Figure 3**. A rectangle marks the location of the seizure area as mapped with invasive subdural grid and strip electrode recordings. The subdural grid was near to the midline covering the right frontal, parietal and temporal areas of the brain. Two strip electrodes were covering the medial frontal and medial parietal areas of the brain. This subject had seizure in central parietal area and frontoparietal midline areas. The scalp dEEG electrodes also showed seizure activities in the same area. These seizures were more toward the right side of the brain from the midline. In **Figure 3**, the rectangle is shown in the 180th second of the plot. A stable pattern of stochastic behavior begins to emerge from 60 s onward. Similar plots were also constructed for other two subjects which are not shown here. The color bar gives the LRTC values, i.e., α values. Within the seizure area, the values of α are in the range of 1–1.5 indicating a strong LRTC which are not of a power law form.

#### **SUBJECT #1**

The stochastic behavior of the SI and cross-frequency couplings for the subject #1 are given in **Figure 4** and **Figure 5**, respectively. For this subject, as shown in **Figure 3**, 60–80 s duration of dEEG data was sufficient to localize the epileptic site. The analysis was done for the dEEG data collected on the first, second, and the eleventh day of the hospital stay. The seizure area is marked with a rectangle that was determined with invasive subdural grid and strip electrode recordings. The stochastic behavior of the SI in the seizure area is higher as compared with nearby surrounding areas on all 3 days in the low gamma band. Refer to **Figure 4**. It is strongest on the first day and fragmented on other days, but, still with hot spots within and near to the seizure area. In the beta band, on the first day, the stochastic activity is very low in the

seizure area. On the second day, it is wide spread, including in the seizure area. On the eleventh day, there is some higher activity at the upper edge of the seizure area. The theta and alpha bands do not show any hot spots in the seizure area.

The stochastic behavior of the cross-frequency couplings is given in **Figure 5**. In general, it is fragmented and less as compared with surrounding areas in all the bands and on all of 3 days, i.e., first, second, and the eleventh day. It is more noticeable in the beta–gamma coupling. The α values for theta–gamma and alpha– gamma couplings are in the range of 0.5–1.2. In comparison, the α values for beta–gamma couplings are higher and are in the range of 1.2–1.6. In **Figure 4**, there are other hot spots outside the seizure in all of these plots which could be related to the spread of the seizure activity in a larger area or due to the some other processes in the brain, such as, spontaneous brain activity. This makes it difficult to localize the seizure area from the stochastic behavior of the SI alone in different EEG bands. However, combining with cross-frequency couplings, it becomes feasible to localize the seizure areas.

#### **SUBJECT #2**

This subject had a complex seizure pattern spreading from right mid temporal area to the right occipital area. The spread of the seizure activity was seen in the grid and strip electrode recordings. The stochastic behavior of the SI is given in **Figure 6** for the first, second, and the third day of the hospital stay. For this subject, 100 s long EEG data was needed to find a stable pattern in the stochastic activity of the SI. The gamma band exhibits hot spots in the seizure are and also in the vicinity of seizure area on all 3 days. Also, there are hot spots in the beta band on the first 2 days in the seizure area. On the third day, the hot spots are toward the lower edge of the rectangle, very close to the occipital area. As compared with low gamma and beta bands, the activity in alpha and theta bands is less intense. For the first 2 days, the alpha band exhibits some low intensity activity in the seizure area and on the third day it is missing. The theta band does not show any activity in the seizure area.

The stochastic behavior of cross-frequency couplings is given in **Figure 7** for the subject #2. For the beta–gamma coupling within the seizure area, the activity is fragmented and exhibits a dipolar and quadrupolar patterns on all 3 days. There is also significant activity outside the seizure area on the first day which could be due to the spreading of the seizure activity and/or due to spontaneous brain activity. On the first day, there is some low level stochastic activity in the alpha–gamma coupling which is reduced

significantly on the following 2 days. The theta–gamma coupling is absent in the seizure area on all 3 days. Within the seizure area the values of α are in the range of 1.0–2.0 for alpha–gamma and beta–gamma couplings. In contrast for theta–gamma couplings, the values of α are in the range of 0.5–1.0. These plots show that the stochastic behavior of the cross-frequency coupling in the beta–gamma band is highest as compared with theta–gamma and alpha–gamma couplings.

#### **SUBJECT #3**

This subject had seizures in the left front temporal area with slow bilateral spread which were observed in subdural grid electrodes, strip electrodes and also in scalp dEEG recordings. The subdural grid was placed on left frontal and temporal areas and electrode strips were placed on left medial temporal, lateral frontal and occipital areas. The stochastic behavior of the SI is given in **Figure 8**

for the first, second, and the third day of the hospital stay. For this subject, 100–120 s long EEG data was needed to find a stable pattern in the stochastic activity of the SI. The gamma band exhibits hot spots in the seizure area on all 3 days. On the first day, the hot spot is observed only at the lower edge of the seizure area rectangle. The hot spots in the vicinity of right occipital area are probably due to the spread of the seizure activity. These seizures related hot spots are much stronger on the second day and also there is widespread activity in the frontal and right temporal areas which was also observed in strip electrodes. The low gamma band activity is very strong in the seizure area marked by the rectangle and also outside, below the rectangle. The beta band exhibits hot spots in the seizure area on the first 2 days while on the third day it is very subdued. The alpha band activity within the seizure area is absent on the first and the third day while some low level activity is present on the second day. Patterns similar to the alpha band

activity are also present in the theta band activity. No theta band activity in the seizure area on the first and the third day. There is a widespread activity in the frontal area, including the seizure area on the second day in beta and low gamma bands. This could be related to bilateral spread as observed in invasive grid and strip electrode recordings.

The stochastic behavior of cross-frequency couplings for the subject #3 is given in **Figure 9**. The stochastic behavior of the beta–gamma coupling within the seizure area is fragmented and exhibits complex spatial patterns on the first and the third day. On the second day, the stochastic behavior of the beta–gamma coupling is spread in a large area on the left side including the seizure area. The stochastic behavior of the alpha–gamma coupling is absent in the seizure area on the first and the second day, but some activity in visible on the third day. The stochastic behavior of the theta–gamma coupling is absent in the seizure area on the first day, slightly visible on the second day and becomes stronger on the third day.

#### **DISCUSSIONS**

For the subject #1, the stochastic behavior of the SI in the low gamma band (**Figures 3** and **4**) seems to show some higher patterns of activity in and near to the seizure area on all 3 days. In addition, the stochastic behavior of the beta–gamma coupling is fragmented in the seizure area. Also, behaviors of other crossfrequency couplings are also weaker in the seizure area with respect to the nearby areas. For this subject a stable pattern of stochastic activity (**Figure 3**) begins to appear from 60 s onward and thus, 80 s long interictal dEEG data was enough to localize the epileptic site.

For the subject #2, both the beta and the low gamma bands had the higher stochastic activity of the SI in and near to the seizure area and the stochastic behavior of the coss-frequency couplings was fragmented and relatively decreased. For this subject, 100 s long interictal dEEG data was needed to find stable patterns of stochastic activity in the seizure area.

For the subject #3 also, both the beta and the low gamma bands had the higher stochastic activity of the SI in and near to the seizure area and the stochastic behavior of the coss-frequency couplings was fragmented and relatively decreased. For this subject, 100– 120 s long interictal dEEG data was needed to find stable patterns of stochastic activity in the seizure area. An interesting feature to note is that for subjects #2 and #3 on the same days when the beta band activity is higher, in comparison, the low gamma band activity is lower.

From these results and from our previous work (Ramon and Holmes, 2012, 2013) one can conclude that in the seizure areas

and also in the vicinity of seizure areas, stochastic behavior of the SI in beta and low gamma bands are higher and also the stochastic behavior of cross-frequency couplings has decreased and spatial patterns are fragmented. Our results also suggest that the seizure related stochastic activity is present on a continuous basis in the interictal scalp EEG data which, possibly, could be useful for non-invasive localization of the epileptic sites. These are our preliminary results and show a promise that these have a potential to localize epileptic sites from a short-duration (1– 3 min), seizure-free, high density (256-channel) scalp EEG data. Further studies with more subjects are needed to substantiate these findings.

The SI is a time series sequence and in the seizure area, the values of α are in the range of 0.8–2.0 for beta and low gamma bands. Thus, in these two bands, SI exhibits a LRTC or strong correlations and also exhibit some properties of Brownian noise within the seizure area. In contrast, theta and alpha bands have lower values of α in the range of 0.4–1.2. This would suggest that the stochastic behavior of the SI in theta and alpha bands is an uncorrelated white noise in some parts of seizure area and in other parts, it has LRTC or strong correlations. Based on this, one can further infer that beta and low gamma bands in combination with cross-frequency couplings will be a better choice for localization of the epileptic sites.

The SI for a given electrode is averaged with the SI of nearby six electrodes. This averaged SI represents the local interconnected activity of cortical neurons. The stochastic behavior of the SI showed in **Figures 3**–**9** is often spread in large areas including the seizure areas. It could possibly represent the local and global connectivity of cortical neurons in the seizure areas as well as normal parts of the brain in the vicinity of the seizure area. Epileptic and interconnected normal cortical neurons, both play an important role in the genesis and spread of the seizure activity. It is also possible that the stochastic behavior of the SI in interictal periods is spread in an area larger than the seizure area. These concepts need to be further examined with seizure modeling and patient studies.

#### **REFERENCES**


oscillations. *Front. Physiol.* 3:450. doi:10.3389/fphys.2012.00450


An increase in fluctuations in the order parameter related to brain synchronization has also been observed and can be used as precursors for identification of epileptic seizures (Velazquez et al., 2011). The order/disorder parameter can also be used as a measure of EEG phase transition states of the brain. In view of these previous studies, our results based on the stochastic fluctuations in the phase synchronization indices are plausible.

derived from short duration interictal scalp dEEG. *Brain Topogr.* 26, 1–8.


stereo-EEG study. *Neuroimage* 61, 591–598.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 15 February 2013; paper pending published: 17 March 2013; accepted: 30 April 2013; published online: 16 May 2013.*

*Citation: Ramon C and Holmes MD (2013) Stochastic behavior of phase synchronization index and cross-frequency couplings in epileptogenic zones during interictal periods measured with scalp dEEG. Front. Neurol. 4:57. doi: 10.3389/fneur.2013.00057*

*This article was submitted to Frontiers in Epilepsy, a specialty of Frontiers in Neurology.*

*Copyright © 2013 Ramon and Holmes. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*