# NEW INSIGHTS ON BASIC AND CLINICAL ASPECTS OF EEG AND MEG CONNECTOME

EDITED BY : Ryouhei Ishii, Roberto D. Pascual-Marqui, Leonides Canuet, Jing Xiang and William C. Gaetz PUBLISHED IN : Frontiers in Human Neuroscience

#### Frontiers Copyright Statement

© Copyright 2007-2018 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-88945-531-7 DOI 10.3389/978-2-88945-531-7

#### 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<EMAIL_ADDRESS>
## NEW INSIGHTS ON BASIC AND CLINICAL ASPECTS OF EEG AND MEG CONNECTOME

Topic Editors:

Ryouhei Ishii, Osaka University, Ashiya Municipal Hospital, Japan Roberto D. Pascual-Marqui, Universität Zürich, Switzerland; Kansai Medical University, Japan

Leonides Canuet, Universidad de La Laguna, Spain

Jing Xiang, Cincinnati Children's Hospital Medical Center, United States William C. Gaetz, Children's Hospital of Philadelphia, United States

Summary of the main statistically significant results comparing the network properties between eyes open and closed conditions. During eyes closed, the posterior cingulate significantly sends mostly alpha oscillations to all other regions. During eyes open, the anterior cingulate significantly sends mostly theta-alpha oscillations to the dorsolateral pre-frontal cortices. PCC, posterior cingulate cortex; ACC, anterior cingulate cortex; LIPL, RIPL, left and right inferior parietal lobule; LDLPFC, RDLPFC, left and right dorsolateral pre-frontal cortex. Image: Pascual-Marqui et al, 2014.

Pascual-Marqui RD, Biscay RJ, Bosch-Bayard J, Lehmann D, Kochi K, Kinoshita T, Yamada N and Sadato N (2014) Assessing direct paths of intracortical causal information flow of oscillatory activity with the isolated effective coherence (iCoh). *Front. Hum. Neurosci.* 8:448. doi: 10.3389/fnhum.2014.00448.

Recent advances in the neuroimaging field areas allow us to visualize the aggregate of neural connections at the macroscopic level within the brain, the so-called "connectome". In order to promote the development of the neurophysiological investigation of connectome of brain oscillations, this eBook aims at bringing together contributions from researchers in basic and clinical neuroscience using EEG and MEG connectome analysis. The most important focal point will be to address the functional roles of connectome of brain oscillations in contributing to understandings of higher cognitive processes in normal subjects and pathophysiology of psychiatric diseases. This Research Topic presented novel methodologies and various applications of neurophysiological connectome analysis. As a result, these papers were cited more than 120 times in these four years in total and threw light and impact on new directions for investigating the connectome of human brain.

Citation: Ishii, R., Pascual-Marqui, R D., Canuet, L., Xiang, J., Gaetz, W. C., eds. (2018). New Insights on Basic and Clinical Aspects of EEG and MEG Connectome. Lausanne: Frontiers Media. doi: 10.3389/978-2-88945-531-7

# Table of Contents

#### *05 Editorial: New Insights on Basic and Clinical Aspects of EEG and MEG Connectome*

Ryouhei Ishii, Roberto D. Pascual-Marqui, Leonides Canuet, Jing Xiang and William C. Gaetz

#### SECTION 1

#### EXPANDING THE BOUNDARIES OF BRAIN CONNECTOME ANALYSIS

Yasunori Aoki, Ryouhei Ishii, Roberto D. Pascual-Marqui, Leonides Canuet, Shunichiro Ikeda, Masahiro Hata, Kaoru Imajo, Haruyasu Matsuzaki, Toshimitsu Musha, Takashi Asada, Masao Iwase and Masatoshi Takeda

*31 LORETA EEG Phase Reset of the Default Mode Network* Robert W. Thatcher, Duane M. North and Carl J. Biver

### SECTION 2

#### APPLICATION OF CONNECTIVITY ANALYSIS FOR BASIC COGNITIVE STUDIES

Ryouhei Ishii, Leonides Canuet, Tsutomu Ishihara, Yasunori Aoki, Shunichiro Ikeda, Masahiro Hata, Themistoklis Katsimichas, Atsuko Gunji, Hidetoshi Takahashi, Takayuki Nakahachi, Masao Iwase and Masatoshi Takeda

*75 Magnetoencephalography Detection of High-Frequency Oscillations in the Developing Brain*

Kimberly Leiken, Jing Xiang, Fawen Zhang, Jingping Shi, Lu Tang, Hongxing Liu and Xiaoshan Wang

#### SECTION 3

#### NEW RESEARCH PROPOSAL AND REVIEW

## Editorial: New Insights on Basic and Clinical Aspects of EEG and MEG Connectome

Ryouhei Ishii 1,2 \*, Roberto D. Pascual-Marqui 3,4, Leonides Canuet 5 , Jing Xiang6 and William C. Gaetz 7

*1 Department of Psychiatry, Osaka University Graduate School of Medicine, Suita, Japan, 2 Neuroscience Center, Ashiya Municipal Hospital, Ashiya, Japan, 3 The KEY Institute for Brain-Mind Research, University of Zurich, Zurich, Switzerland, 4 Neuropsychiatry, Kansai Medical University, Moriguchi, Japan, 5 Department of Cognitive, Social and Organizational Psychology, La Laguna University, Tenerife, Spain, 6 Cincinnati Children's Hospital Medical Center, Cincinnati, OH, United States, 7 Children's Hospital of Philadelphia, Philadelphia, PA, United States*

Keywords: EEG, MEG, connectome, LORETA, default mode network (DMN)

**Editorial on the Research Topic**

#### **New Insights on Basic and Clinical Aspects of EEG and MEG connectome**

Previous studies have shown that, in electroencephalography (EEG) and magnetoencephalography (MEG), cortical oscillations in specific frequency bands (i.e., delta, theta, alpha, beta, and gamma) are functionally related to cognitive processing and behavior, and that abnormal patterns correlate with pathophysiological processes of neuropsychiatric disorders. Although these results have provided useful insights for neural communication in the human brain, they appeared to be highly speculative based on the level of EEG and MEG analyses performed (simple averaging method, wavelet analysis on the sensor level, etc.). The development of new analysis methods for EEG and MEG signals has allowed researchers to depict the information code of brain networks, by precisely localizing sources of oscillatory activity related to brain functions or pathological processes. However, this kind of approach, characterizing brain activity purely in terms of anatomically segregated responses, is not sufficient to explain the pathophysiology of complex neuropsychiatric disorders or the mechanisms underlying cognitive functions.

Recent advances in the neuroimaging field areas allow us to visualize the aggregate of neural connections at the macroscopic level within the brain, the so-called "connectome." In order to promote the development of the neurophysiological investigation of connectome of brain oscillations, this e-book aims at bringing together contributions from researchers in basic and clinical neuroscience using EEG and MEG connectome analysis. The most important focal point will be to address the functional roles of connectome of brain oscillations in contributing to understandings of higher cognitive processes in normal subjects and pathophysiology of psychiatric diseases.

The included papers can be roughly divided into three groups as follows: the first group of papers both tried to expand the boundaries of brain connectome analysis by applying new statistic solvers on exact low resolution electromagnetic tomography (eLORETA) analysis of EEG. Pascual-Marqui et al. who founded LORETA himself, presented namely the "isolated effective coherence" (iCoh) obtained from eLORETA, which consists of estimating the partial coherence under a multivariate autoregressive model, followed by setting all irrelevant associations to zero, other than the particular directional association of interest. They elucidated the direct directed connection of alpha activity from the posterior cingulate to all other regions during eyes closed and theta-alpha activity from the anterior cingulate to other frontal regions during eyes open in human

#### Edited by:

*Vladimir Litvak, Institute of Neurology, University College London, United Kingdom*

> Reviewed by: *Hisato Sugata, Oita University, Japan*

\*Correspondence: *Ryouhei Ishii<EMAIL_ADDRESS>
Received: *04 May 2018* Accepted: *22 May 2018* Published: *05 June 2018*

#### Citation:

*Ishii R, Pascual-Marqui RD, Canuet L, Xiang J and Gaetz WC (2018) Editorial: New Insights on Basic and Clinical Aspects of EEG and MEG Connectome. Front. Hum. Neurosci. 12:232. doi: 10.3389/fnhum.2018.00232* brain. Aoki et al. applied eLORETA-ICA analysis to resting-state EEG data in 80 healthy subjects using five frequency bands (delta, theta, alpha, beta and gamma band) and found five restingstate-independent-networks in alpha, beta and gamma frequency bands. Thatcher et al. analyzed the phase shift and lock duration of 3-dimensional current sources in 14 Brodmann areas located in the DMN using LORETA in the delta frequency band by using the Hilbert transform between all pairs of Brodmann areas. By depicting an inverse relation of phase shift and lock durations in an exponential way with distance between Brodmann areas, they showed a tremendous view of anatomical hubs in the brain behave like a "shutter" that opens and closes at specific durations as nodes of brain network. These cutting-edge methodologies were often applied in the recent publications and cited almost 60 times in these four years in total.

The application of these connectivity analysis for basic cognitive studies was described in the following papers the focus of the second group. To examine the relationship between functional connectivity and performance of brainmachine interfaces (BMIs), Sugata et al. analyzed the correlation coefficient between performance of neural decoding and functional connectivity over the whole brain using. They suggest that use of the strength of functional connectivity between M1 and motor association areas has the potential to improve the performance of BMIs to perform real and imagined movements. Asakawa et al. evaluated the effect of different anxiety states on information processing as measured by EEG using emotional stimuli on a smartphone. They found that the propagation speed of the low anxiety group at the medial coronal for resting stimuli for all time segments was higher than those of high anxiety group and suggested that neural information processes concerning emotional stimuli differ based on current anxiety state. Ishii et al. used spatially filtered MEG and permutation analysis to precisely localize cortical generators of the magnetic counterpart of frontal midline theta rhythm (Fmθ) and found the dorsal anterior cingulate and adjacent medial prefrontal cortex as the generetors of Fmθ and gamma event-related synchronization in right parietal regions subserving basic numerical processing and number-based spatial attention. They suggested that these multiple oscillatory activities might interact each other to proceed these cognitive tasks in the brain. Leiken et al. reported high-frequency brain signals in the developing brain using time-frequency analysis along with beamforming methods on MEG data. They suggested that the developing brain generates high-frequency signals that can be detected with the non-invasive technique of MEG. These studies showed several type of connectivity among certain brain areas in multiple frequency bands such as theta, alpha, beta, gamma bands related to various types of cognitive process and behavior.

The third group contains the new research proposal from Ibanez et al. assessing the combination of basal resting state influences together with the ongoing activity during a task and its evoked neural response to characterize brain dynamics changes related to preclinical Alzheimer's disease and the review article from Miki and Kakigi introducing their own three studies that focused on facial movements.

In summary, this ebook presented novel methodologies and various applications of neurophysiological connectome analysis. As a result, these papers were cited more than 120 times in these 4 years in total and threw light and impact on new directions for investigating the connectome of human brain.

### AUTHOR CONTRIBUTIONS

RI and LC discussed about this research topic and wrote the manuscript. RP-M, JX, and WG gave significant advise and helped the editing the manuscript.

### ACKNOWLEDGMENTS

This study was partly supported by a Viera y Clavijo fellowship from Tenerife, Spain.

**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.

Copyright © 2018 Ishii, Pascual-Marqui, Canuet, Xiang and Gaetz. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

## Assessing direct paths of intracortical causal information flow of oscillatory activity with the isolated effective coherence (iCoh)

#### *Roberto D. Pascual-Marqui 1,2\*†, Rolando J. Biscay3†, Jorge Bosch-Bayard4†, Dietrich Lehmann1, Kieko Kochi 1, Toshihiko Kinoshita2, Naoto Yamada5 and Norihiro Sadato6*

*1 The KEY Institute for Brain-Mind Research, University of Zurich, Zurich, Switzerland*

*3 CIMFAV, Universidad de Valparaiso, Valparaiso, Chile*

*4 Neuroinformatic Department, Cuban Neuroscience Center, Havana, Cuba*

*5 Department of Psychiatry, Shiga University of Medical Science, Shiga, Japan*

*6 Division of Cerebral Integration, National Institute for Physiological Sciences, Okazaki, Japan*

#### *Edited by:*

*Ryouhei Ishii, Osaka University Graduate School of Medicine, Japan*

#### *Reviewed by:*

*Guido Nolte, Fraunhofer FIRST, Germany Gennady Knyazev, Academy of Medical Sciences, Russia*

#### *\*Correspondence:*

*Roberto D. Pascual-Marqui, The KEY Institute for Brain-Mind Research, University Hospital of Psychiatry, Lenggstrasse 31, Zurich CH-8032, Switzerland*

*e-mail<EMAIL_ADDRESS>
*†These authors have contributed equally to this work.*

Functional connectivity is of central importance in understanding brain function. For this purpose, multiple time series of electric cortical activity can be used for assessing the properties of a network: the strength, directionality, and spectral characteristics (i.e., which oscillations are preferentially transmitted) of the connections. The partial directed coherence (PDC) of Baccala and Sameshima (2001) is a widely used method for this problem. The three aims of this study are: (1) To show that the PDC can misrepresent the frequency response under plausible realistic conditions, thus defeating the main purpose for which the measure was developed; (2) To provide a solution to this problem, namely the "isolated effective coherence" (iCoh), which consists of estimating the partial coherence under a multivariate autoregressive model, followed by setting all irrelevant associations to zero, other than the particular directional association of interest; and (3) To show that adequate iCoh estimators can be obtained from non-invasively computed cortical signals based on exact low resolution electromagnetic tomography (eLORETA) applied to scalp EEG recordings. To illustrate the severity of the problem with the PDC, and the solution achieved by the iCoh, three examples are given, based on: (1) Simulated time series with known dynamics; (2) Simulated cortical sources with known dynamics, used for generating EEG recordings, which are then used for estimating (with eLORETA) the source signals for the final connectivity assessment; and (3) EEG recordings in rats. Lastly, real human recordings are analyzed, where the iCoh between six cortical regions of interest are calculated and compared under eyes open and closed conditions, using 61-channel EEG recordings from 109 subjects. During eyes closed, the posterior cingulate sends alpha activity to all other regions. During eyes open, the anterior cingulate sends theta-alpha activity to other frontal regions.

**Keywords: causal intracortical connectivity, LORETA, isolated effective coherence, resting state electriphysiological connectivity, alpha oscillation connectivity**

#### **INTRODUCTION**

The type of problem that we are interested in can best be understood with an informal hypothetical example.

Consider time series of local electric potential differences measured at five sites (i.e., nodes) on the cortex (electrocorticogram, ECoG). Before connecting the five nodes, each one in isolation has its distinct activity. For instance, node 1 oscillates at 28 Hz, node 2 at 16 Hz, and nodes 3, 4, and 5 at 23 Hz. In the next construction step, some causal direct and directional connections with measurable time lags are established: node 1 sends to node 2; and node 2 sends identically to nodes 3, 4, and 5. The resulting connectivity graph is shown in **Figure 1**.

Instantaneous connections, as considered, e.g., by Faes et al. (2013), are not considered in this hypothetical example, i.e., ephaptic conduction is assumed to be absent, see e.g., Weiss et al. (2013).

Note the distinction between "direct" and "indirect" connection paths. Examples in this hypothetical network are: (A) The direct connection path from node 1 to node 2; (B) The indirect connection path from node 1 to node 3 mediated by node 2.

Time series measurements from this hypothetical network can be generated by means of a multivariate autoregressive model, as will be shown in a quantitatively precise manner below (Equations 1, 13).

*2 Department of Neuropsychiatry, Kansai Medical University, Osaka, Japan*

Given only the time series measurements, the problem of interest here is to recover the detailed properties of the network consisting of all the activity properties at each node, and the nature of the direct causal connections, i.e., their strength, direction, and spectral characteristics of the oscillations that are being transmitted.

This study is limited to this type of problem.

Moreover, this rather narrow and simple problem is, to this date, of great interest, as can be seen, for instance, in a recent publication by Plomp et al. (2014), which provides a brief review of methods and offers some benchmark data which will later be used in this present study.

Upon reviewing the history and state of the art in frequency domain causal connectivity studies, there are at least two noteworthy contributions to this field:

Recently, the PDC has been critically studied by Schelter et al. (2009). They pointed out that the normalization used in the PDC, i.e., the denominator in the PDC formula (see details below), contains influences from the sender node of interest to all receiver nodes, and as a consequence, the PDC decreases in the presence of many nodes, even if the relationship between a sender and receiver of particular interest remains unchanged. The solution to this problem was given in the form of a renormalization of the PDC, using the statistical variance of the strength of the connection.

In this present study, rather the aiming at a re-normalization of the PDC, such as that successfully achieved by Schelter et al. (2009), we reformulate the problem from scratch, estimating the partial coherence under a multivariate autoregressive model, followed by setting all irrelevant associations to zero, other than the particular directional association of interest. This procedure is akin to Pearl's "surgical intervention" for studying causality (Pearl, 2000). This approach gives the isolated effective coherence (iCoh) (Pascual-Marqui et al., 2014).

In the original Baccala and Sameshima paper (Baccala and Sameshima, 2001), a number of simple toy examples were designed to illustrate the superiority of the PDC as compared to other competing methods. Following the same style, we here provide a new simple toy example, which compellingly shows how the PDC can give incorrect information about the strength of a connection, and incorrect information on its spectral characteristics. And we show how the iCoh solves this problem.

To further illustrate the shortcomings of PDC as pointed out by Schelter et al. (2009), both PDC and iCoh are compared below in the analysis of publicly available benchmark data, consisting of somatosensory responses in rats, measured on 15 skull electrodes. The two methods produce dramatically different results, with much reduced PDC values in some cases, consistent with the observation made by Schelter et al. (2009), as discussed above.

In order to test the new iCoh measure as compared to the PDC under more adverse realistic conditions, the following simulation was performed. Five time series with known dynamics were generated, and used as the electrical activity assigned to 5 cortical sites. EEG recordings were computed by solving the forward equations (see e.g., Fuchs et al., 2002; Gomez-Herrero et al., 2008) at 19 scalp electrodes from these sources. The EEG was then given to an inverse solver, namely eLORETA [exact low resolution electromagnetic tomography (Pascual-Marqui, 2007, 2009; Pascual-Marqui et al., 2011)], producing estimated cortical signals which were then used to compute iCoh and PDC. As predicted and as will be shown below, iCoh recovers adequately all the information about the network, whereas PDC also does so, but reporting false results.

In a final example, real human recordings are analyzed, where the iCoh between six cortical regions of interest (ROIs) are calculated and compared under eyes open and closed conditions, using 61-channel EEG recordings from 109 subjects (EEGs from public data base, see Goldberger et al., 2000; Schalk et al., 2004). The ROIs consist of the anterior and posterior cingulate cortices, the inferior parietal lobules, and the dorsolateral pre-frontal cortices. Statistical comparisons for every pair of ROIs, and for every discrete frequency were based on non-parametric randomization of the maximum-statistic (see e.g., Nichols and Holmes, 2002), thus ensuring correction for multiple testing. During eyes closed, the posterior cingulate significantly sends alpha activity to all other regions. During eyes open, the anterior cingulate significantly sends theta-alpha activity to the dorsolateral pre-frontal cortices.

For the sake of reproducible research, the software code implementing the methods discussed here (using lazarus free-pascal "www.lazarus.freepascal.org"), including test data as a simple text file, are freely available at: https://sites.google.com/site/ pascualmarqui/home/icoh-isolated-effective-coherence.

#### **METHODS**

#### **THE MULTIVARIATE AUTOREGRESSIVE MODEL**

As described above, the definition of the iCoh is based on formulating a multivariate autoregressive model, and calculating the corresponding partial coherences after setting all irrelevant connections to zero. All technical details can be found in Pascual-Marqui et al. (2014). For the sake of completeness, a brief presentation is included here. General background and notation on multivariate autoregressive models, frequency domain causality, and spectral density matrices can be found, for instance, in Akaike (1968) and Yamashita et al. (2005).

A stable, stationary multivariate autoregressive model of order *p* ≥ 1, for *q* ≥ 2 time series **X**(*t*) ∈ R*q* × 1, is written as:

$$\mathbf{X}(t) = \sum\_{k=1}^{p} \mathbf{A}(k)\mathbf{X}(t-k) + \mathbf{e}(t) \tag{1}$$

where **A**(*k*) ∈ R*q*×*q* are the autoregressive coefficients, **ε**(*t*) ∈ R*q*×1 is the innovations (noise) vector, and *t* denotes discrete time.

In general, the autoregressive coefficients [**A**(*k*)]*ij*, i.e., the element *i*, *j* of the matrices **A**(*k*), quantify the direct causal influence for *j* → *i*. This corresponds to Granger causality (Granger, 1969; Lütkepohl, 2007; Valdes-Sosa et al., 2011).

Given data sampled in discrete time, and given an autoregressive order *p* ≥ 1, the autoregressive coefficients and the innovation covariance matrix can be estimated by any number of methods, one of which is least squares (see e.g., Akaike, 1968). The model order *p* can be estimated by means of Akaike's information criterion AIC (Akaike, 1974).

The frequency domain representation is:

$$\mathbf{X}(\omega) = \mathbf{A}(\omega)\mathbf{X}(\omega) + \mathbf{e}(\omega) \tag{2}$$

where *X* (ω) ∈ C*q*×1, *A* (ω) ∈ C*q*×*q*,ε(ω) ∈ C*q*×1 are the discrete Fourier transforms, and where ω denotes discrete frequency.

From Equation 2, the Hermitian covariance, i.e., the spectral density matrix, is:

$$\mathbf{S}\_{\mathbf{x}}(\omega) = \left(\check{\mathbf{A}}(\omega)\right)^{-1} \mathbf{S}\_{\boldsymbol{\ell}}\left(\check{\mathbf{A}}^{\*}(\boldsymbol{\alpha})\right)^{-1} \tag{3}$$

with :

$$
\dot{\mathbf{A}}(\omega) = \mathbf{I} - \mathbf{A}(\omega) \tag{4}
$$

where the superscript "∗" denotes matrix transpose and complex conjugate, the superscript "−1" denotes matrix inversion, **I** is the identity matrix, and **S**ε ∈ R*q*×*q* is the noise covariance.

#### **THE ISOLATED EFFECTIVE COHERENCE (ICoH)**

From the spectral density matrix (Equation 3), the partial coherences (see e.g., Brillinger, 2001) between any pair of nodes *i*, *j*  can be calculated. The significance of the partial coherence in a very general setting can be found in Radhakrishna Rao (1981). In simple terms, the partial coherence is a measure of association between two complex valued random variables after removing the effect of other measured variables.

The full general equation for the partial coherence as a function of all the autoregressive coefficients contains information on all possible connection paths. Technical details can be found in Equations 8, 9 within Pascual-Marqui et al. (2014). However, in order to "isolate" the direct and directional parts of a connection, all other possible paths must be severed. This is a procedure commonly used in causality analysis, metaphorically known as performing a "surgical intervention" (see e.g., Pearl, 2000).

For this reason, the isolated effective coherence (iCoh) for *j* → *i* is defined under the condition that the only non-zero association between the time series is due to **A**˘ (ω) *ij* = 0. This requires that all other possible associations be set to zero, i.e.,:

$$[\mathbf{A}(\omega)]\_{kl} \equiv \mathbf{0} \quad , \text{ for all } (k, l) \text{ such that}$$

$$(k, l) \neq \begin{pmatrix} i, \ j \end{pmatrix} \text{ and } k \neq l \tag{5}$$

and :

[**S**ε]*kl* ≡ 0 , for all (*k*, *l*) such that *k* = *l* (6)

Note that the diagonal elements of **S**ε and **A** (ω) remain unmodified, since they do not "associate" different nodes.

Emphasis must be placed on the fact that this procedure is meaningful only if the new system with a single association remains stable and stationary.

When the constraints in Equations 5, 6 are applied to the general partial coherence, we obtain the isolated effective coherence (iCoh). In particular, iCoh for *j* → *i* is defined as the squared modulus of the partial coherence between *i* and *j* under the constraints given by Equations 5, 6:

$$\kappa\_{i \gets j}(\boldsymbol{\omega}) = \frac{\left[\mathbf{S}\_{\boldsymbol{\varepsilon}}\right]\_{\boldsymbol{i}\boldsymbol{i}}^{-1} \left| \left[\check{\mathbf{A}}\left(\boldsymbol{\omega}\right)\right]\_{\boldsymbol{i}\boldsymbol{j}} \right|^{2}}{\left[\mathbf{S}\_{\boldsymbol{\varepsilon}}\right]\_{\boldsymbol{i}\boldsymbol{i}}^{-1} \left| \left[\check{\mathbf{A}}\left(\boldsymbol{\omega}\right)\right]\_{\boldsymbol{i}\boldsymbol{j}} \right|^{2} + \left[\mathbf{S}\_{\boldsymbol{\varepsilon}}\right]\_{\boldsymbol{j}\boldsymbol{i}}^{-1} \left| \left[\check{\mathbf{A}}\left(\boldsymbol{\omega}\right)\right]\_{\boldsymbol{j}\boldsymbol{i}} \right|^{2}} \tag{7}$$

which clearly satisfies:

$$0 \le \kappa\_{i \gets j}(\alpha) \le 1 \tag{8}$$

The detailed, step by step derivations are shown in Pascual-Marqui et al. (2014).

The iCoh can be described as the answer to the following question:

"Given a dynamic linear system characterized by its autoregressive parameters, what would be the equation for the partial coherence if all connections are severed, except for the single one of interest?"

Note that the algorithm for computing the iCoh requires:

(1) The estimation of the full, joint, multivariate autoregressive model (Equation 1). This step is performed only once.

(2) For any given pair of nodes and any direction such as *j* → *i*, compute Equation 7 using the parameters from step (1).

#### **THE PARTIAL DIRECTED COHERENCE (PDC) AND THE GENERALIZED PARTIAL DIRECTED COHERENCE (gPDC)**

These definitions are replicated here for the sake of completeness. The PDC is:

$$\left| \bar{\boldsymbol{\pi}}\_{\dot{\boldsymbol{\eta}}} (\boldsymbol{\omega}) \right|^{2} = \frac{\left| \left[ \check{\mathbf{A}} (\boldsymbol{\omega}) \right]\_{\dot{\boldsymbol{\eta}}} \right|^{2}}{\left[ \check{\mathbf{A}}^{\*} (\boldsymbol{\omega}) \check{\mathbf{A}} (\boldsymbol{\omega}) \right]\_{\dot{\boldsymbol{\eta}}}} = \frac{\left| \left[ \check{\mathbf{A}} (\boldsymbol{\omega}) \right]\_{\dot{\boldsymbol{\eta}}} \right|^{2}}{\sum\limits\_{k=1}^{q} \left| \left[ \check{\mathbf{A}} (\boldsymbol{\omega}) \right]\_{k\boldsymbol{\eta}} \right|^{2}} \tag{9}$$

which corresponds to Baccala and Sameshima (2001), Equation 18 therein.

The gPDC is:

$$\left| \pi\_{ij}^{\rm w}(\omega) \right|^{2} = \frac{\left[ \mathbf{S}\_{\boldsymbol{\varepsilon}} \right]\_{\vec{\boldsymbol{\mu}}}^{-1} \left| \left[ \check{\mathbf{A}}(\omega) \right]\_{\vec{\boldsymbol{y}}} \right|^{2}}{\sum\limits\_{k=1}^{q} \left[ \mathbf{S}\_{\boldsymbol{\varepsilon}} \right]\_{\vec{\boldsymbol{k}}}^{-1} \left| \left[ \check{\mathbf{A}}(\omega) \right]\_{\vec{\boldsymbol{k}}} \right|^{2}} \tag{10}$$

which corresponds to Baccalá et al. (2007), Equation 11 therein.

Note that these measures are not proper partial coherences. The squared modulus of a proper coherence or partial coherence has a value between zero and one, and they do not have a column or row sum value of 1.

#### **STATISTICS**

In some instances, results will be presented simply as the estimated values of connectivities, without performing an actual statistical test. This type of result is akin to showing the effect size.

In other cases, where specified, statistical tests are carried out based on the method of non-parametric randomization of the maximum-statistic, which has the advantage of correcting for multiple testing, and of not relying on Gaussianity (Blair and Karniski, 1994; Karniski et al., 1994; Nichols and Holmes, 2002; Nichols, 2012).

A brief description of the multivariate non-parametric randomization method follows. Technical details are not included here because they can found in the specialized literature, see e.g., Nichols and Holmes (2002) and the cited literature therein. Consider an example where the data is represented as *Xcki*, consisting of *i* = 1 ... *R* variables, measured on *k* = 1 ... *N* subjects, under two conditions *c* = 1 and *c* = 2. The variables may correspond to cortical spectral power at each voxel and each frequency, or to direct and directed connection strength between each pair of regions of interest and each frequency.

In this example, the aim is the discovery of the variables that are significantly different between the two conditions. For this purpose, the simple variable-by-variable *t*-statistic can be used as a statistical measure of "distance" between the two conditions. Other choices of statistics are equally valid. From the set of "*R*" *t*-statistics (one for each variable), the absolute maximum is chosen. Then its empirical probability distribution is estimated by repeatedly randomizing the conditions "*c*," and recalculation the maximum-*t*'s under the null hypothesis. This empirical probability gives the threshold with correction for multiple ("*R*" tests) testing, as explained in Nichols and Holmes (2002). The correction is exact (in the sense of Fisher's exact test) for a large number of randomizations, regardless of the original probability distribution of the variables.

#### **EEG: FORWARD AND INVERSE PROBLEMS**

The equation of electrodynamics that links current density in the brain to scalp electric potential differences is known as the "forward" equation of EEG. This forward problem, which has a well-defined solution, is typically solved with numerical methods.

Simulated EEG is easily created by placing sources of time varying electric neuronal activity at any number of cortical sites, and calculating the electric potential differences on scalp electrodes, by means of the forward equation.

Formally, the forward equation of EEG in discrete form at time instant "*t*" can be written as:

$$\Phi\_t = \mathbf{K} \mathbf{J}\_t \tag{11}$$

where *t* ∈ R*NE* × 1 denotes the instantaneous scalp electric potential at *NE* electrodes, **J***t* ∈ R(3*NV* )×1 is the instantaneous current density vector field at *NV* cortical voxels (consisting of three components at each voxel), and **K** ∈ R*NE* × (3*NV* ) is the lead field.

The inverse problem, which consists of estimating the cortical activity (current density vector field) from measured scalp EEG, is known to have no unique solution (see e.g., Helmholtz, 1853; Pascual-Marqui, 2009). This is the reason for the existence of many different inverse solutions found in the literature. In this study, the method known as exact low resolution electromagnetic tomography (eLORETA; Pascual-Marqui, 2007, 2009; Pascual-Marqui et al., 2011) is used for estimating sources in both simulated EEG and for real human EEG measurements. The eLORETA solution has the following generic form:

$$\hat{\mathbf{J}}\_t = \mathbf{T} \Phi\_t \tag{12}$$

where **T** ∈ R(3*NV* )×*NE* is the eLORETA pseudoinverse (Pascual-Marqui, 2007; Pascual-Marqui et al., 2011).

In the current implementation of eLORETA, computations are made in a realistic head model (Fuchs et al., 2002), using the MNI152 template (Mazziotta et al., 2001), with the threedimensional solution space restricted to cortical gray matter, as determined by the probabilistic Talairach atlas (Lancaster et al., 2000). The standard electrode positions on the MNI152 scalp were taken from Jurcak et al. (2007). A total of 6239 cortical gray matter voxels at 5 mm spatial resolution constitute the solution space.

The estimated time varying electric neuronal activity at each cortical voxel (given by ˆ **J***t* in Equation 12) consists of three time series, one for each moment component of the current density vector (i.e., dipole). In practice, this can be reduced to a single time series, due to the fact that the current density vector is anatomically constrained to have an orientation orthogonal to the cortical surface (see e.g., Baillet et al., 2001). Under this assumption, the 3 × 3 covariance matrix for the current density vector at each voxel must have rank 1, with the dipole orientation given by its largest eigenvector (Mosher et al., 1992; Mosher and Leahy, 1998). This procedure is applied in this study for the estimation of single time series of electric neuronal activity at each voxel. Note that this estimator for the current density vector field orientation is a maximum variance estimator.

#### **MATERIALS**

#### **A TOY EXAMPLE: FIVE TIME SERIES**

Simulated recordings from five time series were generated from the following stable, stationary multivariate autoregressive model of order 2:

$$\mathbf{A}(1) = \begin{pmatrix} 1.5 & -0.25 & 0 & 0 & 0 \\ -0.2 & 1.8 & 0 & 0 & 0 \\ 0 & 0.9 & 1.65 & 0 & 0 \\ 0 & 0.9 & 0 & 1.65 & 0 \\ 0 & 0.9 & 0 & 0 & 1.65 \end{pmatrix};$$

$$\mathbf{A}(2) = \begin{pmatrix} -0.95 & 0 & 0 & 0 & 0 \\ 0 & -0.96 & 0 & 0 & 0 \\ 0 & -0.8 & -0.95 & 0 & 0 \\ 0 & -0.8 & 0 & -0.95 & 0 \\ 0 & -0.8 & 0 & 0 & -0.95 \end{pmatrix};$$

$$\mathbf{S}\_{\varepsilon} = \mathbf{I} \tag{13}$$

The direct causal directed connections between nodes are illustrated as arrows in **Figure 1**.

Assuming a Gaussian distribution for the noise (zero mean, unit variance, as shown in Equation 13), 25600 time samples were generated (after discarding the first 1000 time samples) and used for all estimation procedures.

Assuming that the times series are sampled at 256 Hz, the main spectral properties of this network, by construction, are the following:

#### **SIMULATED EEG**

In a different setting, the five time series generated in the previous subsection were used as the time varying electric neuronal activities at the following cortical locations:

where (X, Y, Z) denotes the MNI coordinates in millimeters, and BA denotes Brodmann area.

**Figure 2** illustrates the five cortical locations.

Use was made of the forward equations previously explained for generating EEG recordings at 19 scalp electrodes, corresponding to the 10/20 electrode placement system. In this generation process, two relatively large sources of noise were added:

#### **RAT EEG**

EEG recorded at 2 kHz sampling rate from 15 electrodes placed directly onto the skull of rats, during a somatosensory experiment, are publicly available from Plomp et al. (2014). A single recording from their repository, the data from the file named "RN060915A2\_STIMD," was taken for analysis. This corresponds to the average evoked response for one particular animal, with unilateral whisker stimulation. The recording starts at −60 ms relative to stimulus onset, and has a total duration of 180 ms.

#### **HUMAN EEG RECORDINGS**

Real human EEG recordings under eyes open and closed conditions, using 64-channel EEG recordings from 109 subjects, are publicly available from Goldberger et al. (2000), Schalk et al. (2004). Each recording (218 in total) consists of 1 min. EEG, sampled at 160 Hz. Three electrodes (T9, T10, and Iz) were discarded for analysis, because they were spatial outliers relative to the other 61 electrodes that cover the scalp in an approximate uniformly distributed manner.

#### **RESULTS**

#### **A TOY EXAMPLE: FIVE TIME SERIES**

**Figure 3** shows the iCoh (Equation 7) and the gPDC (Equation 10) calculated for the network in **Figure 1**. In both cases, the same estimated multivariate autoregressive model of order *p* = 3 was used. The results were essentially identical for autoregressive order *p* = 2.

Of importance to note in **Figure 3**: the two methods give very different results with respect to node #2 as sender (column 2).

#### **SIMULATED EEG**

The simulated EEG time series for 19 scalp electrodes, using as generators the five cortical locations described in the materials section (**Figure 2**), with time dynamics from the previous example, were analyzed with eLORETA.  However, connectivity computations are presented for the estimated electrical activities at the same cortical sites as in **Figure 2**. **Figure 4** shows the estimated iCoh and gPDC. In both cases, the same estimated multivariate autoregressive model of order *p* = 3 was used. The results were essentially identical for autoregressive order *p* = 2.

Of importance to note in **Figure 4**:

#### **RAT EEG**

The average somatosensory evoked response for one rat was analyzed with a multivariate autoregressive model of order *p* = 8, based on the model order determined in the original publication. Although this data is clearly not stationary, it was analyzed as such in the original publication (Plomp et al., 2014), using a recursive least squares (RLS) algorithm with a forgetting factor, in order to implement a time varying version of the autoregressive model.

#### **FIGURE 3 | Estimated connectivity properties for the network in Figure 1.** Isolated effective coherence (iCoh) shown in RED, and the

generalized partial directed coherence (gPDC) shown in BLUE. Overlap of both curves is shown in BLACK. Vertical axis: 0 to 1. Frequency axis: 1 to

127 Hz. Columns are senders, rows are receivers. Coherence peak in column 1 occurs at 28 Hz. Coherence peak for iCoh in column 2 occurs at 16 Hz. Coherence peak for gPDC in column 2, row 1 occurs at 1 Hz; and Coherence peak for gPDC in column 2, rows 3, 4, and 5 occur at 23 Hz.

The only particular and differentiating feature in this current study is that there is no forgetting factor.

But regardless of these considerations, the sole purpose of this rat data analysis is to show the extreme differences in network properties estimated by iCoh and gPDC, as shown in **Figure 5**.

Of importance to note in **Figure 5**: the extremely low connectivity values produced by gPDC as compared to iCoh, and the number of missing spectral peaks in gPDC as compared to iCoh.

#### **HUMAN EEG RECORDINGS**

EEGs recorded from 109 subjects under eyes open and eyes closed conditions were analyzed. Resting state, awake, eyes closed EEG is characterized by the presence of alpha rhythm, as compared to the eyes open condition.

In a first analysis step, the spectral density of electric neuronal activity throughout the cortex was calculated with eLORETA at all 6239 cortical voxels. The technical details on calculating cortical activity spectra can be found in Frei et al. (2001). A voxel by voxel, frequency by frequency comparison between eyes open and closed conditions was performed. **Figure 6** shows the three main statistically significant results. Eyes open is characterized by significantly stronger activity in frontal cortical regions oscillating at 3 Hz and in the beta band 23–28 Hz. Eyes closed is characterized by significantly stronger activity in occipital cortical regions oscillating at 10 Hz.

In a second analysis step, time series of electric neuronal activity were estimated with eLORETA at 6239 cortical voxels, from which six cortical regions of interest were used for the connectivity analyses. This procedure was applied to the EEGs recorded in 109 subjects, under eyes open and eyes closed conditions. The

iCoh was estimated for the six time series, in the 218 recordings, using an autoregressive order *p* = 7, which corresponds to the median order for all EEG recordings based on Akaike's AIC (see subsection "The multivariate autoregressive model"). A statistical comparison between eyes open and eyes closed conditions was carried out, for each frequency, for each pair of regions of interest, and for each direction of connection. The significant differences at probability 0.05 with correction for multiple testing, are shown in **Figure 7**.

**Figure 8** summarizes the main statistically significant results. During eyes closed, the posterior cingulate significantly sends activity to all other regions. During eyes open, the anterior cingulate significantly sends activity to the dorsolateral pre-frontal cortices.

**FIGURE 5 | Estimated connectivity properties for rat EEG recordings from 15 skull electrodes.** Isolated effective coherence (iCoh) shown in RED, and the generalized partial directed coherence (gPDC) shown in BLUE. Vertical axis: 0 to 1. Frequency axis: 7.8 to 250 Hz. Columns are senders, rows are receivers.

#### **DISCUSSION**

The results corresponding to the toy example with five time series demonstrate very clearly a major problem with the partial directed coherence (PDC; Baccala and Sameshima, 2001) as well as with its generalized version (gPDC; Baccalá et al., 2007). By construction, node 2 sends identical information to nodes 3, 4, and 5. And yet the gPDC gives very different results. Moreover, the frequency responses of the gPDC are incorrect, with one missing spectral peak and other peaks at incorrect frequencies.

This type of problem was already pointed out by Schelter et al. (2009), and here we show a compelling example of how severe it can be.

In contrast, the isolated effective coherence (iCoh) introduced in this study recovers and reports correctly all the network properties for the toy example.

The results corresponding to the evoked response recordings from an animal experiment (Plomp et al., 2014) demonstrate that gPDC and iCoh can give very different results with real experimental data. Because of the complex biological nature of the data, the ground truth is not unambiguously known. Regardless, this demonstrates that when there are many nodes (15 in this case), the gPDC can give very low connectivity values and almost featureless spectral properties as compared to iCoh.

By its very nature and definition, the gPDC does not report information on a particular direct directed connection. Instead, it is affected by many other connections, in such a way that it can report incorrect values and spectral properties of the "senderreceiver" pair of interest. The iCoh solves this problem, by its very nature and definition (Pascual-Marqui et al., 2014).

In the neurosciences, one of the most interesting applications of a method such as the iCoh is the elucidation of effective cortical connections based on measurements of electric neuronal activity. However, these are extremely invasive measurements. In order to solve this problem non-invasively, one possible approach is to use scalp EEG measurements, and to estimate with an inverse solution the electric neuronal activity at any number of cortical locations. It is then very important to prove that the estimated time series are of sufficient quality to calculate iCoh reliably.

This was the aim of the experiment with simulated EEG. Cortical signals were used for computing EEG, which was then analyzed with the eLORETA inverse solution (which has no prior information about the locations or about the dynamics of the actual sources). Despite the low spatial resolution of eLORETA, and despite the use of only 19 scalp electrodes, iCoh was estimated very reliably. This is partly due to the fact that a measure such as iCoh separates rather well instantaneous and lagged connections, especially if the instantaneous connections are mediated by the noise covariances, which explicitly do not affect iCoh (see Equations 6, 7). However, the low spatial resolution of eLORETA will mix the autoregressive coefficients, as is shown in Gomez-Herrero et al. (2008). Both Gomez-Herrero et al. (2008) and Faes et al. (2013) propose solutions to this problem, which can be applied also to eLORETA signals.

Finally, an eLORETA-iCoh study was performed on real human EEG which is available from a public repository (Goldberger et al., 2000; Schalk et al., 2004). The research aim here was to search for differences in brain function between two resting states, namely eyes open and eyes closed. Two aspects of brain function were explored:

1. The cortical location of the generators of different oscillatory activity (functional localization).

2. The network properties among a group of six very important cortical sites (functional "effective" connectivity).

This type of study is of interest in understanding the resting state of the brain. In particular, in understanding the functional role of the alpha rhythm (Knyazev et al., 2011; Klimesch, 2012; Sigala et al., 2014), and in understanding the functional changes during the eyes open condition (Jao et al., 2013). The results show that eyes closed alpha activity was localized to occipital areas, while delta and beta activities were located in frontal cortical regions. With respect to the network properties, the iCoh analysis demonstrated that the posterior cingulate cortex is a major sender of mainly alpha oscillations to all other regions. Interestingly, during eyes open, this function is turned off, and the anterior cingulate activates as a sender of mainly theta-alpha oscillations to the dorsolateral pre-frontal cortices.

#### **OUTLOOK AND LIMITATIONS**

The iCoh method can be extended to other conditions, different from the particular ones considered here. For instance,

**FIGURE 7 | t-Statistics comparing eyes open minus eyes closed iCoh for 109 subjects, in six regions of interest: anterior cingulate, posterior cingulate, left and right inferior parietal, and left and right dorsolateral pre-frontal cortices.** Frequency axis: 1 to 30 Hz. Corrected *p* = 0.05 was at t-threshold = 4.3, with vertical axis: −7 to +7. Blue

color denotes eyes closed significantly larger, red color denotes eyes open significantly larger. The three numbers indicate the frequencies (Hz) for the significant results: start, end, and the most significant oscillation indicated with a superscript "∗." Columns are senders (prefix "s"), rows are receivers (prefix "r").

instantaneous connections such as those considered in the generalized autoregressive model of Faes et al. (2013) can be directly applied to iCoh. Moreover, since iCoh solely depends on the estimated autoregressive coefficients and the noise variances, these parameters can be estimated under non-stationary, time-varying conditions, as for example in Plomp et al. (2014).

However, if the actual dynamics are non-linear or simply do not follow a linear autoregressive model, then iCoh might be invalid. Non-linear causality measures have been reviewed in Marinazzo et al. (2011), where a novel method is proposed: "kernel Granger causality." Another method is the phase slope index (Nolte et al., 2008), which is of a more nonparametric nature, not relying on the parametric form of the linear autoregression. However, these two alternative methods do not distinguish the direct or indirect nature of the connections.

Interestingly, a very recent book entitled "Directed Information Measures in Neuroscience" (Wibral et al., 2014) barely deals with methods that reveal all properties of a neural network, namely the spectral content of information flow, the direct or indirect nature of the connections, and the actual direction. One exception is a single chapter that refers to the PDC of Baccala and Sameshima (2001), and to Geweke's method (Geweke, 1984) (which is based on the predictive approach of Granger).

In our present study, the method of Geweke was not studied, and certainly deserves more attention in future research. However, we note that Geweke's method has been criticized elsewhere (Chen et al., 2006) because it often produces negative connectivity values that render it meaningless.

It is important to emphasize that the EEG simulation example presented here is very limited, and only constitutes a "proof of principle," since the cortical signals used for analysis were close to the actual locations of the sources. The effect of the choice of the number of regions of interest and of their locations relative to the actual unknown active network requires further study.

One common problem in all models that depend on fitting a multivariate autoregressive model is the curse of dimensionality: for a large number of nodes and for a high autoregressive order, the number of parameters to be estimated can be too large to produce reliable estimators. One possible solution is the estimation of sparse multivariate autoregressions as developed by Valdes-Sosa et al. (2005). Alternatively, stable high dimensional autoregressive models can be successfully estimated under spatiotemporal constraints, such as those considered by Jiménez et al. (1995).

The eLORETA method was used in this study. Other inverse solutions can be used. The only requirement is that the selected method needs to be capable of correct estimation of the neuronal current density. This was the reason for choosing eLORETA, because it is an improvement over the previous related tomographies known as LORETA (Pascual-Marqui et al., 1994) and sLORETA (Pascual-Marqui, 2002), which have received considerable and substantial validation (Pascual-Marqui et al., 2011). We note that all these techniques can equally be applied to MEG measurements as well.

**conditions.** During eyes closed, the posterior cingulate significantly sends mostly alpha oscillations to all other regions. During eyes open, the anterior cingulate significantly sends mostly theta-alpha oscillations to the dorsolateral pre-frontal cortices. PCC, posterior cingulate cortex; ACC, anterior cingulate cortex; LIPL, RIPL, left and right inferior parietal lobule; LDLPFC, RDLPFC, left and right dorsolateral pre-frontal cortex.

There is a severe limitation in the use and interpretation of all connectivity measures (including iCoh) if they are applied to scalp EEG signals. In this case, the results should never be interpreted as representing cortical connections. The reason is that cortical activity does not project radially onto the scalp (see e.g., Lehmann et al., 2006, 2012). This problem applies to all EEG analyses that naively map scalp measurements and features onto the underlying cortex, which in general produce incorrect results.

In conclusion, iCoh is most certainly not intended as the general solution to the problem of identifying network properties. It is a very simple and particular measure for correctly assessing direct connections that causally transmit oscillatory information between nodes, under the assumption of a linear autoregressive model. It is distinct from the PDC method of Baccala and Sameshima (2001); Baccalá et al. (2007), which can produce incorrect results.

#### **ACKNOWLEDGMENTS**

Rolando J. Biscay was partly supported by the CONICYT ACT 1112 and CONACYT 131771 projects. Part of this work was performed while Roberto D. Pascual-Marqui was at the Department of Psychiatry in Shiga University of Medical Science.

#### **REFERENCES**

functional brain mapping. *Hum. Brain Mapp.* 10, 120–31. doi: 10.1002/1097- 0193(200007)10:3%3C120::AID-HBM30%3E3.0.CO;2-8

**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: 08 May 2014; accepted: 03 June 2014; published online: 20 June 2014. Citation: Pascual-Marqui RD, Biscay RJ, Bosch-Bayard J, Lehmann D, Kochi K, Kinoshita T, Yamada N and Sadato N (2014) Assessing direct paths of intracortical causal information flow of oscillatory activity with the isolated effective coherence (iCoh). Front. Hum. Neurosci. 8:448. doi: 10.3389/fnhum.2014.00448 This article was submitted to the journal Frontiers in Human Neuroscience.*

*Copyright © 2014 Pascual-Marqui, Biscay, Bosch-Bayard, Lehmann, Kochi, Kinoshita, Yamada and Sadato. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with*

*accepted academic practice. No use, distribution or reproduction is permitted which*

*does not comply with these terms.*

## Detection of EEG-resting state independent networks by eLORETA-ICA method

#### **Yasunori Aoki 1 , Ryouhei Ishii 1\*, Roberto D. Pascual-Marqui 2,3 , Leonides Canuet 4 , Shunichiro Ikeda1 , Masahiro Hata1 , Kaoru Imajo5 , Haruyasu Matsuzaki 6 , Toshimitsu Musha6 , Takashi Asada7 , Masao Iwase1 and Masatoshi Takeda1**

#### **Edited by:**

Jing Xiang, Cincinnati Children's Hospital Medical Center, USA

**Reviewed by:** Fawen Zhang, University of Cincinnat, USA Mitsuru Kikuchi, Kanazawa University, Japan

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

Ryouhei Ishii, Department of Psychiatry, Osaka University Graduate School of Medicine, D3 2-2 Yamada-oka, Suita, Osaka 565-0871, Japan e-mail<EMAIL_ADDRESS>Recent functional magnetic resonance imaging (fMRI) studies have shown that functional networks can be extracted even from resting state data, the so called "Resting State independent Networks" (RS-independent-Ns) by applying independent component analysis (ICA). However, compared to fMRI, electroencephalography (EEG) and magnetoencephalography (MEG) have much higher temporal resolution and provide a direct estimation of cortical activity. To date, MEG studies have applied ICA for separate frequency bands only, disregarding cross-frequency couplings. In this study, we aimed to detect EEG-RS-independent-Ns and their interactions in all frequency bands. We applied exact low resolution brain electromagnetic tomography-ICA (eLORETA-ICA) to restingstate EEG data in 80 healthy subjects using five frequency bands (delta, theta, alpha, beta and gamma band) and found five RS-independent-Ns in alpha, beta and gamma frequency bands. Next, taking into account previous neuroimaging findings, five RS-independent-Ns were identified: (1) the visual network in alpha frequency band, (2) dual-process of visual perception network, characterized by a negative correlation between the right ventral visual pathway (VVP) in alpha and beta frequency bands and left posterior dorsal visual pathway (DVP) in alpha frequency band, (3) self-referential processing network, characterized by a negative correlation between the medial prefrontal cortex (mPFC) in beta frequency band and right temporoparietal junction (TPJ) in alpha frequency band, (4) dual-process of memory perception network, functionally related to a negative correlation between the left VVP and the precuneus in alpha frequency band; and (5) sensorimotor network in beta and gamma frequency bands. We selected eLORETA-ICA which has many advantages over the other network visualization methods and overall findings indicate that eLORETA-ICA with EEG data can identify five RS-independent-Ns in their intrinsic frequency bands, and correct correlations within RS-independent-Ns.

**Keywords: eLORETA-ICA, LORETA, resting state network, independent component analysis, ICA, EEG**

#### **INTRODUCTION**

The brain intrinsically interacts between distant regions, building cortical networks during motor and cognitive tasks. Interestingly, one network enhances its activity in no-task resting state. In particular, the so called default mode network (DMN) is known to be active during resting and attenuates during task performance. However, recent findings suggest that the DMN is also involved in internally focused processes such as self-referential thoughts, envisioning one's future and autobiographical memory retrieval (Raichle et al., 2001; Buckner et al., 2008). Furthermore, it has been reported that several cortical networks cooperate with each other positively or negatively during performance of complex cognitive tasks (Spreng and Schacter, 2012). These functional networks have been investigated by lesional and anatomical studies and during functional tasks with functional magnetic resonance imaging (fMRI), which measures regional cerebral blood flow (rCBF) changes. However, one mathematical method called independent component analysis (ICA) have received growing attention (Bell and Sejnowski, 1995; Hyvärinen and Oja, 2000). ICA is a mathematical decomposing method which separates mixture of signals like electroencephalography (EEG), magnetoencephalography (MEG) and fMRI data into a set of statistical independent components (ICs) that are artifact signals and physiological signals. In addition, it should be noted that using ICA these task positive or negative functional networks can be extracted from "resting state" fMRI data and MEG data (Beckmann et al., 2005; Allen et al., 2011; Brookes et al., 2011). These led to the concept of "Resting State independent Network" (RS-independent-N). Also, there are some other methods used for the discovery of interactions in the brain which are seedbased correlation analyses. These analyses has extracted Resting State correlated Networks (RS-correlated-Ns) from resting state fMRI data or MEG data (Biswal et al., 1995; Vincent et al., 2008; Brookes et al., 2011; Raichle, 2011; Hipp et al., 2012). In this way, ICA and seed-based correlation analyses with fMRI data has identified several RS-independent-Ns and RScorrelated-Ns, including the basal ganglia network, auditory network, sensorimotor network, visual network, DMN, ventral and dorsal visual pathway (VVP and DVP), and the frontal network (Biswal et al., 1995; Allen et al., 2011; Joel et al., 2011; Raichle, 2011; Meyer et al., 2013). However, correlation analysis has a problem of an implicit assumption of Gaussianity of the signal where fMRI signals are approximately Gaussian (Hlinka et al., 2011) but EEG and MEG signals are non-Gaussian (Stam, 2005). Thus, RS-correlated-Ns derived from correlation analysis of EEG and MEG data are not independent with each other in a precise sense because of non- Gaussianity of EEG and MEG data (Hyvärinen and Oja, 2000; Stam, 2005). In addition, correlation analyses emphasize the special role of some preselected brain region. However, unlike the seed-based methods, ICA is appropriate for the discovery of distributed networks, giving equal importance to all brain voxels (Joel et al., 2011). Furthermore, ICA can remove artifacts such as electromyogram or base line shift by separating out artifact components (Custo et al., 2014).

Unlike fMRI, which measures hemodynamic changes that occur in response to cortical activity, neurophysiological techniques, such as EEG and MEG measure cortical electrical/magnetic activity directly and noninvasively with a high temporal resolution (1–2 ms) (Canuet et al., 2011). Thus, EEG has been widely used in clinical practice to support clinical diagnosis and management of neuropsychiatric diseases such as epilepsy, disturbance of consciousness and dementia, and also in neuroscience to investigate cortical electrical activities and functions (Ishii et al., 1999; Canuet et al., 2011; Kurimoto et al., 2012; Aoki et al., 2013a,b).

Recent findings of EEG and MEG analyses indicate that electromagnetic oscillatory activity of the functional networks varies its frequency from lower sensory areas to higher-order control areas. For instance, intra-cortical investigations using depth electrodes with syllable auditory task reported that cortical electrical activity of auditory area changed from evoked activity (phase-locked to the stimulus) to induced activity (non-phaselocked to the stimulus) and also its frequency changed from theta and low gamma to beta and high gamma, as activity shifted from primary auditory cortices to associative auditory cortex (Morillon et al., 2012). Another MEG study using a visuospatial attentional task found that the cortical electrical activity of the DVP changed from alpha evoked activity to beta induced activity as it shifted from early visual areas to prefrontal control areas (Siegel et al., 2008). And recent fMRI and MEG studies using decomposing methods have repeatedly shown that these functional networks can also be seen during resting state with changing its power of activity (Smith et al., 2009; Grady et al., 2010; Brookes et al., 2011). From these accumulating evidences, we can assume that RS-independent-Ns are associated with several frequency bands of electromagnetic activity depending on the function subserved by the different cortical regions. In support of this notion, a simultaneous fMRI and EEG study showed that blood oxygenation level dependent (BOLD) signals of RS-independent-Ns correlated with EEG waveforms in several frequency bands (Mantini et al., 2007). In addition, Jonmohamadi et al. (2014) and Mantini et al. (2011) showed that ICA decomposition of EEG and MEG data becomes more correct in localization and more robust to artifacts when applied after source reconstruction. Taken together, in order to visualize RS-independent-Ns across several frequency bands, we consider appropriate to apply ICA to cortical electrical activity reconstructed from EEG or MEG data, analyzing all frequency bands. To our knowledge, there is one previous EEG-RS-independent-N study. However, ICA was applied to scalp recorded EEG data in the time domain, followed by a second step using a sLORETA source reconstruction on the ICA-scalp topographies; in contrast, we apply ICA directly to the reconstructed cortical electrical activity by eLORETA in the frequency domain. And the results of cortical electrical distributions of ICs were rather different from known functional networks (Chen et al., 2013). Also there is a few previous MEG-RS-independent-N studies. In their studies, ICA was applied to cortical electrical activity reconstructed from MEG data, however, in separate frequency bands, disregarding possible cross-frequency coupling. Furthermore, sample sizes of these studies were small (Brookes et al., 2011, 2012; Luckhoo et al., 2012).

Also, ICA of EEG data has been widely used for various purposes, such as artifact rejection by separating out artifact components (Custo et al., 2014) and examination of the EEG resting states (infra-slow EEG fluctuations and EEG microstates). For instance, Hiltunen et al. (2014), found correlations between the filtered ICA time series (using ultra-low frequencies) of the EEG with BOLD time series in specific fMRI RS-independent-Ns. And Yuan et al. (2012), performing ICA on EEG microstates to decompose into ICs (independent microstates), found that each fMRI RS-independent-N was characterized by one to a combination of several independent microstates.

Exact low resolution brain electromagnetic tomography (eLORETA) is a linear inverse solution method that can reconstruct cortical electrical activity with correct localization from the scalp EEG data (Pascual-Marqui et al., 2011; Aoki et al., 2013a). The implementation of ICA in the eLORETA software with EEG data allows for decomposition of cortical electrical activity which is non-Gaussian into ICs in different frequency bands (Pascual-Marqui and Biscay-Lirio, 2011). Other decomposing methods (e.g., principal component analysis or correlation analysis) with EEG data cannot strictly to do so (Bell and Sejnowski, 1997; Hyvärinen and Oja, 2000; Mantini et al., 2011). Furthermore, electromagnetic tomography-ICA (eLORETA-ICA) uses all frequency information of EEG data in analysis. In this study, we selected eLORETA-ICA which has many advantages over the other network visualization methods as we explained above and applied it to EEG data to obtain complete set of EEG-RS-independent-Ns across several frequency bands for the first time.

### **METHODS**

#### **SUBJECTS**

We recruited 306 healthy elderly subjects who had no history of neurological or psychiatric disorders. Elderly subjects over 60 years old underwent clinical tests to ensure that memory and other cognitive functions were within normal limits (MMSE > 24, CDR = 0). From the participants, 146 subjects were healthy volunteers, and the remaining 160 subjects were ascertained from an epidemiological study among inhabitants in Tone, Ibaraki, Japan. This study was approved by the Ethics Committee of Osaka University Hospital and followed the Declaration of Helsinki. Written informed consent was obtained from the subjects.

#### **EEG RECORDING AND DATA ACQUISITION**

Subjects underwent EEG recordings in a resting state, eyes closed condition for about 5 min. Subjects were instructed to keep their eyes closed while staying awake during the recordings. Spontaneous cortical electrical activity was recorded with a 19-channel EEG system (EEG-1000/EEG-1200, Nihon Kohden, Inc., Tokyo, Japan), filtered through a 0.53–120 Hz band-pass filter, and sampled at 500 Hz. EEG was recorded with the electrodes positioned according to the International 10–20 system (i.e., Fp1, Fp2, F3, F4, C3, C4, P3, P4, O1, O2, F7, F8, T3, T4, T5, T6, Fz, Cz, Pz) using a linked ears reference. Electrode impedances were kept below 5 kΩ. For each subject, 120-s artifact-free, resting-awake segments were manually selected by visual inspection using Neuroworkbench software (Nihon Kohden, Inc., Tokyo, Japan).

#### **EEG-SOURCE RECONSTRUCTION METHOD**

We used eLORETA (exact low resolution brain electromagnetic tomography) to compute the cortical electrical distribution from the scalp electrical potentials measured at the electrode sites (Pascual-Marqui et al., 2011). The eLORETA method is a weighted minimum norm inverse solution, where the weights are unique and endow the inverse solution with the property of exact localization for any point source in the brain. Thus, due to the principles of linearity and superposition, any arbitrary distribution will be correctly localized, albeit with low spatial resolution. In the current eLORETA version, the solution space consists of 6239 cortical gray matter voxels at 5 mm spatial resolution, in a realistic head model (Fuchs et al., 2002), using the MNI152 template (Mazziotta et al., 2001). The LORETA method has been previously used and validated with real human data during diverse sensory stimulation and in neuropsychiatric patients (Dierks et al., 2000; Vitacco et al., 2002; Pascual-Marqui et al., 2011; Aoki et al., 2013a). A further property of eLORETA is that it has correct localization even in the presence of structured noise (Pascual-Marqui et al., 2011). In this sense, eLORETA is an improvement over previously related versions of LORETA (Pascual-Marqui et al., 1994) and sLORETA (Pascual-Marqui, 2002). eLORETA images of spectral density were computed for five frequency bands: delta (2–4 Hz), theta (4–8 Hz), alpha (8– 13 Hz), beta (13–30 Hz), and gamma (30–60 Hz) (Canuet et al., 2012).

#### **FUNCTIONAL ICA**

In most of the resting state network (RSN) literature, ICA is the method most widely used for the discovery of sets of regions that work together as networks. There are numerous different processing strategies that are being used in the RS-independent-N literature, as reviewed by Calhoun (Calhoun et al., 2009). For instance, in typical fMRI group studies for the discovery of RS-independent-Ns, the time series images for each subject are first heavily pre-processed (see Calhoun et al., 2009 for details), and then all subjects' time series images are concatenated. This produces a matrix, where one dimension consists of "space" (i.e., the brain voxels), and the other dimension consists of time. Finally, an ICA algorithm is applied to this matrix, which will produce a set of spatial components (i.e., images), where each "component image" consists of a so-called "network". In order to interpret a network image, one must threshold its values appropriately, displaying the brain regions that have highest loadings. This post-processing is achieved by z-transforming the component network image, and using an empirical threshold, as in for example (McKeown et al., 1998; Calhoun et al., 2004; Kelly et al., 2010; Agcaoglu et al., 2014). In this way, each network image will display areas whose activities are tightly linked (i.e., they work together as a network).

In contrast to relatively slow hemodynamic images, high time resolution images of electrical neuronal activity can be computed using eLORETA applied to EEG recordings. In an implicit manner, these images contain an additional dimension of frequency. Whereas fMRI images have their spectrum concentrated below 0.1 Hz, EEG contain a wealth of differential functional information in the classical range from 2 to 60 Hz. In order to take into account this additional dimension of information, the classical ICA as applied in fMRI was generalized. All the technical details can be found in Pascual-Marqui and Biscay-Lirio (2011).

For the sake of completeness, a brief description follows. The EEG recording of each subject is first transformed to the frequency domain, using the discrete Fourier transform. This will produce a set of cross-spectral EEG matrices, for each frequency of interest, such as those described above. This information is then used for calculating the spectral density for each cortical voxel and for each frequency band, using the methodology described in detail in Frei et al. (2001). With this initial procedure, each subject contributes five eLORETA images of cortical spectral density (one for each frequency band: delta, theta, alpha, beta, and gamma). From the point of view of mathematics, these data correspond to a "function" of space (cortical voxel) and frequency. In the next step, the data from each subject is concatenated, thus producing a matrix where one dimension corresponds to the different subjects, and the other dimension corresponds jointly to space-frequency. This approach is common in a relatively young field of statistics known as functional data analysis (Ramsay and Silverman, 2005). When ICs analysis is applied to this matrix, a more general form of networks are discovered, and the method

is described as functional ICA, given its origin in the field of functional data analysis. Each functional network consists of a set of five images, one for each frequency, because space and frequency and all their possible interactions are now jointly expressed. In contrast to a classical fMRI network image which corresponds to brain regions that "work" together over time, an EEG-eLORETA based functional network corresponds to brain regions and frequencies that "work" together across a population of subjects. This allows not only for the discovery of regions that work together, but also for the discovery of cross-frequency coupling.

In this paper, the number of ICs (networks) is estimated by sphericity test (Bartlett, 1954). In the eLORETA-ICA algorithm, ICs were obtained by maximizing the independence between components which was measured by non-Gaussianity. In particular, non-Gaussianity was calculated by fourth-order cumulant (Cardoso, 1989; Cichocki and Amari, 2002). Then, ICs were ranked according to total EEG power and color coded with a z-score threshold of 3.0, in complete analogy to the methods used in practice in fMRI-ICA networks (as explained in detail above). In the color–coded maps, red and blue colors represent power increase and decrease with increasing IC coefficient which indicates activity of IC, respectively.

#### **RESULTS**

Artifact-free 120-s epochs were obtained in 80 out of 306 healthy subjects. The age distribution of the 80 healthy subjects (57 men and 23 women) was as follows: 18–29 years (25 men and 2 women), 30–49 years (15 men and 4 women), 50–69 years (14 men and 11 women) and 70–87 years (3 men and 6 women) (44 ± 20 (mean ± standard deviation)). The median of MMSE scores over 60 years old was 30 (interquartile range; 29–30). It can be seen an overall male predominance, which may reflect a bias of our healthy volunteers, and the female predominance in the 70–87 years group, which may reflect a delay of age-related cognitive decline in female. The number of ICs estimated by the sphericity test was 12.0. Subsequently, we applied eLORETA-ICA as the number of components varied from 11 to 13. Then, 11 ICs were most consistent with physiological assumption that is topography and frequency of some known networks and artifacts such as electromyogram is at frontal or temporal cortex in gamma frequency band, therefore we selected 11 as the number of components. Next, we identified, based on spatial distributions of power and frequency ranges, IC4, IC5, IC6, IC9 and IC10 as RS-independent-Ns (**Figures 1**–**5**); IC1, IC2, IC3, IC7, IC8 and IC11 as artifacts of frontal and temporal electromyogram or frontal and occipital baseline shifts (**Figure 6**).

When identifying the different ICs derived from our analyses, we found that IC4 corresponded to the occipital visual network in alpha frequency band (**Figure 1**). IC5 consisted of the right VVP, corresponding to the right occipitotemporal cortex and the right ventral prefrontal cortex (vPFC), and the left posterior DVP. The right VVP linked right occipitotemporal cortex in alpha frequency band to the right vPFC in beta frequency band. The left posterior DVP, comprised the ipsilateral posterior occipito-parietal cortex, caudal intraparietal sulcus (cIPS) and posterior end of middle temporal gyrus (MT+) in alpha frequency band, which correlated negatively with the areas of the right VVP (**Figure 2**). IC6 was formed by the medial PFC (mPFC) in beta frequency band and the right temporoparietal junction (TPJ) in alpha frequency band, which showed negative correlation (**Figure 3**). IC9 comprised the precuneus in alpha frequency band and the left VVP in alpha frequency band, which also showed negative correlation (**Figure 4**). IC10 comprised the medial postcentral regions (Brodmann area 5 and 7 (BA 5–7)) in beta frequency band and the pre supplementary motor area (pre-SMA) in gamma frequency band, which showed positive correlation (**Figure 5**).

#### **DISCUSSION**

In this study, using eLORETA-ICA, we could identify five RSindependent-Ns corresponding to (1) the occipital visual network in alpha frequency band (IC4), (2) the right VVP in alpha and

beta frequency bands and left posterior DVP in alpha frequency band (IC5), (3) the mPFC in beta frequency band and right TPJ in alpha frequency band (IC6), (4) the precuneus and left VVP in alpha frequency band (IC9); and (5) the medial postcentral regions in beta frequency band and the pre-SMA in gamma frequency band (IC10).

#### **INDEPENDENT COMPONENT 4**

IC4 was found at the occipital cortex in alpha frequency band (**Figure 1**). It is well known that the occipital cortex is involved in visual perception processing. Consistent with this fact and with our result, previous neurophysiological studies found that visual processing related activity in the occipital regions occurred in the alpha frequency band. In particular, alpha oscillation in the occipital regions is enhanced during no expectation of visual stimulus and is reduced during expectation and presentation of visual stimulus (Klimesch et al., 1998).

#### **INDEPENDENT COMPONENT 5**

IC5 was found at the right occipitotemporal cortex in alpha frequency band and at the right vPFC in beta frequency band with left posterior occipito-parietal cortex, cIPS and MT+ in alpha frequency band (**Figure 2**). Left IC5 regions (the left posterior occipito-parietal cortex, cIPS and MT+) corresponds to left posterior DVP and right IC5 regions (the right occipitotemporal cortex, TPJ, parahippocampal gyrus, fusiform gyrus and vPFC) corresponds to right VVP. DVP is a functional network involved in automatic visual guidance of spatial movements. Within this network cIPS and MT+ is linked to action-relevant features of objects such as shape and orientation from visual information

processed in the occipital lobe. Right VVP is a visual recognition network where visual information that has flowed from the occipital lobe is compared to visual/spatial memory in right temporal cortex then identified in right temporal cortex or right vPFC (Fairhall and Ishai, 2007; Kravitz et al., 2011, 2013; Milner, 2012). Taking into account these findings, IC5 corresponds to a network that activity of the right VVP correlated negatively with left posterior DVP activity. Previous accumulating studies revealed that function of DVP is"online" "unconsciously occurred (automatic)" visual perception of spatial components to guide spatial movements, while function of VVP is "offline" "conscious" visual perception and recognition of feature components (Kravitz et al., 2011, 2013; Harvey and Rossit, 2012; Milner, 2012). Therefore, we can assume IC5 as dualprocess of visual perception: the left posterior DVP for automatic visual perception to guide spatial movements and right VVP for detailed perception and recognition of visual input. Our result of negative correlation between right VVP and left posterior DVP is consistent with dual-process of visual perception. In addition, our result of emergence of VVP only on the right side also fit with the fact that right dominant engagement of VVP in visuospatial search and recognition (Corbetta et al., 2005). This negative correlation was also seen in visuospatial neglect patients, who injured right VVP area, enhanced left posterior DVP activity (not whole DVP) at acute stage and attenuated its activity with clinical recovery (Corbetta et al., 2005; He et al., 2007; Rossit et al., 2012).

#### **INDEPENDENT COMPONENT 6**

IC6 was found at the mPFC in beta frequency band and right TPJ in alpha frequency band (**Figure 3**). Medial PFC is anterior hub of the DMN and right TPJ is a hub of the right VAN (Corbetta

and Shulman, 2002; Buckner et al., 2008). Connectivity analysis of resting fMRI data has showed that mPFC has maximal positive connectivity with right posterior TPJ (Mars et al., 2012; Kubit and Jack, 2013). Taking into account these findings, IC6 corresponds to a network that activity of anterior hub of DMN (mPFC) positively correlated with that of right VVP. The DMN enhance its activity in autobiographical memory retrieval (Cabeza et al., 2004). However autobiographical memory retrieval involves both self-referential processing and memory retrieval process. So, Kim (2012), by subtracting fMRI activity in laboratory-based memory retrieval from autobiographical memory retrieval, found that self-referential processing was related to mPFC, right parahippocampal cortex and posterior cingulate cortex (PCC). So, we can speculate IC6 as self-referential processing. In support to this notion, there is a case report of a patient with loss of the sense of self-ownership who also showed hypometabolism in the right inferior temporal cortex as well as in the right parietooccipital junction and precentral cortex (Zahn et al., 2008).

#### **INDEPENDENT COMPONENT 9**

IC9 was found at the precuneus and left VVP in alpha frequency band (**Figure 4**). The precuneus is dominantly related to familiarity of the memory (Yonelinas et al., 2005) and left VVP is memory recognition area whose activation reflects retrieval and identification of memory (Cabeza, 2008; Ravizza et al., 2011; Angel et al., 2013). IC9 showed the precuneus was negatively correlated with left VVP in alpha frequency band. Consistent with our result, EEG study using sLORETA showed the same correlation between decreasing alpha power in the precuneus and increasing alpha power in the left temporal cortex with WM load during WM retention period in some healthy subjects (Michels et al., 2008). Dual-process models of memory recognition have been proposed by many researchers which suggest memory has two separate systems: familiarity of the memory (sense of knowing) and recollection (Yonelinas, 2002). In memory retrieval, the precuneus engages in familiarity, while left VVP regions (left TPJ, parahippocampal cortex and

hippocampal formation) engage in episodic memory retrieval (Yonelinas et al., 2005; Sestieri et al., 2011). Familiarity is a working memory which is a sense of knowing temporarily occurred (several tens of seconds) after encoding. That is, familiarity is "unconsciously occurred (automatic)" "online" "sensory component" of short-term memory to be manipulated in multiple cognitive processes (working memory). On the other hand, episodic memory retrieval is a "conscious" "offline" "detailed" perception and recognition of long-term episodic memory (Baddeley and Hitch, 1974; Huijbers et al., 2010, 2012). Therefore, we can conclude that familiarity and episodic memory have properties of the DVP and the VVP, respectively (please refer to the discussion of IC5). In fact, the precuneus showed strong coherence with DVP by fMRI connectivity analysis (Huijbers et al., 2012). Taken together, we can speculate that IC9 reflects dual-process of memory perception: the precuneus for automatic sensory component of the memory to guide multiple cognitive processes in memory domain and left VVP for detailed perception and recognition of episodic memory. Our results elucidated that similarity of perception and recognition between vision (IC5) and memory (IC9). Lesion studies also presented a case of neglect in memory domain analogous to visuospatial neglect: patients who had bilateral TPJ lesions showed a deficit in detailed memory retrieval in free recall (subserved by the left VVP), although they can access to these memories when guided by probe questions (function subserved by the precuneus; Berryhill et al., 2007; Cabeza, 2008).

#### **INDEPENDENT COMPONENT 10**

IC10 was found at the medial postcentral regions (BA5–7) in beta frequency band, and at the pre-SMA in gamma frequency band (**Figure 5**). Beta activity in medial sensory regions is known as Rolandic beta rhythm, which is typically observed in resting state and suppressed by voluntary movements (Pfurtscheller, 1981). This beta oscillation is thought as an idling rhythm of sensory regions (Ritter et al., 2009). From our result, gamma oscillation in the pre-SMA can also be assumed as idling rhythm of pre-SMA. In support of this notion, the gamma oscillation in the pre-SMA was suppressed by voluntary movements (Hosaka et al., 2014). Taking into account these findings, we identified IC10 as sensorimotor network.

Overall, topographies of alpha and beta frequency bands is consistent with their roles: alpha oscillation for inhibition of the visual pathway (Snyder and Foxe, 2010; Capotosto et al., 2012; Capilla et al., 2014), beta oscillation in PFC for higher cognitive functions such as evaluation and prediction (Arnal et al., 2011; Hanslmayr et al., 2011; Buschman et al., 2012; Aoki et al., 2013b; Kawasaki and Yamaguchi, 2013) and beta oscillation in sensorimotor area for motor control (Engel and Fries, 2010).

This is the first study presenting ICs using eLORETA-ICA with resting state EEG data, and more importantly, which highlight the differences in some aspects from the previous RS-independent-Ns using ICA with resting state fMRI data. First, eLORETA-ICA of EEG data presented right and left VVP separately, strikingly different from ICA results of fMRI data showing VVP bilaterally. However, de Pasquale et al. (2010) using correlation analysis showed that MEG has greater correlations between intrahemispheric nodes than inter-hemispheric nodes in RSNs. They elucidated that this difference stemmed from the difference of temporal resolution: EEG and MEG have much higher temporal resolution (1–2 ms) of cortical activity than fMRI, which has 2 s temporal resolution. These findings indicate that only EEG and MEG, which have millisecond temporal resolution, combined with ICA can detect the correct ICs of cortical activity. Furthermore, our result of right and left separation of VVP is consistent with previous findings that left lateralized activation of VVP during episodic memory retrieval and right lateralized activation of VVP during visual target detection (Corbetta et al., 2005; Daselaar et al., 2006; Angel et al., 2013). Second, our results were restricted to cortical areas whereas RSNs derived from fMRI data included deep brain structures such as basal ganglia, hippocampus and cingulate cortex. This caused from the fact that EEG cannot detect electrical activity of the deep brain because electrical potential drastically attenuated by conduction from deep brain to the surface of the head. Therefore, for instance, we cannot determine the PCC is involved in IC6 or IC9, although controversy exists whether the PCC should be involved in selfreferential processing or episodic memory retrieval (Kim, 2012; Angel et al., 2013).

Although the fact is known that cortical electrical activity reconstructed from EEG data using sLORETA showed several topographic distributions somewhat similar to RS-independent-Ns for a short period (microstate; Musso et al., 2010), no one could extracted independent sets of cortical electrical activity (EEG-RS-independent-Ns). And there are some other decomposing methods such as principal component analysis and correlation analysis, they cannot decompose cortical electrical activity into ICs in a precise sense because cortical electrical activity is non-Gaussian (Bell and Sejnowski, 1997; Hyvärinen and Oja, 2000; Stam, 2005; Mantini et al., 2011). Therefore, we selected eLORETA-ICA to detect EEG-RS-independent-Ns.

Our results should be interpreted with caution based on the following limitations. First, relative small number of electrodes (19 electrodes) and realistic head model in eLORETA may affect the source localization results. However, the good localization property of the LORETA tomography was validated in several studies as we mentioned in the Methods section and our source localization results of eLORETA-ICA are consistent with neuroimaging findings of RSNs. Second, low spatial resolution of eLORETA, which blur the cortical sources, may cause non-detection of the low-electrical-activity cortical sources. Thus, subsequent ICA may have missed some low activity RS-independent-Ns. Third, our present study has made use of the hypothesis that healthy subjects have common RS-independent-Ns which are consistent throughout a very wide age range, thus aging-related changes are restricted to activities of RS-independent-Ns (IC coefficients). However, we confirmed that occipital basic oscillations of all subjects were in the alpha frequency band by visual inspection and almost all the RS-independent-N results did not change even excluding 9 subjects aged 70 years or more from eLORETA-ICA. In addition, our source localization results of eLORETA-ICA are consistent with neuroimaging findings of RSNs. Fourth, we supposed correspondences between RS-independent-Ns and functional networks. However, these correspondences should be confirmed by comparing the activities of RS-independent-Ns with cognitive test scores in the future study.

#### **CONCLUSION**

We selected eLORETA-ICA which has many advantages over the other network visualization methods and overall findings indicate that eLORETA-ICA with EEG data can identify five RS-independent-Ns with their intrinsic oscillatory activities, as well as functional correlations within these networks, while conventional methods used to examine RSNs such as fMRI with functional tasks or fMRI with ICA have not been shown to do so. Moreover, once RS-independent-Ns are determined by eLORETA-ICA, this method can accurately identify activity of each RS-independent-N from EEG data as it removes EEG artifacts by separating artifact components. Therefore, eLORETA-ICA with EEG data may represent a useful and powerful tool to assess activities of RS-independent-Ns, which correspond to specific functions, in patients with neuropsychiatric disease such as dementia and depression.

#### **ACKNOWLEDGMENTS**

This research is partially supported by the Center of Innovation Program from Japan Science and Technology Agency, JST.

#### **REFERENCES**

**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: 24 July 2014; accepted: 12 January 2015; published online: 10 February 2015*. *Citation: Aoki Y, Ishii R, Pascual-Marqui RD, Canuet L, Ikeda S, Hata M, Imajo K, Matsuzaki H, Musha T, Asada T, Iwase M and Takeda M (2015) Detection of EEGresting state independent networks by eLORETA-ICA method. Front. Hum. Neurosci. 9:31. doi: 10.3389/fnhum.2015.00031*

*This article was submitted to the journal Frontiers in Human Neuroscience*.

*Copyright © 2015 Aoki, Ishii, Pascual-Marqui, Canuet, Ikeda, Hata, Imajo, Matsuzaki, Musha, Asada, Iwase and Takeda. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution and reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms*.

### LORETA EEG phase reset of the default mode network

#### *Robert W. Thatcher\*, Duane M. North and Carl J. Biver*

*EEG and NeuroImaging Laboratory, Applied Neuroscience Research Institute, Seminole, FL, USA*

#### *Edited by:*

*Jing Xiang, Cincinnati Children's Hospital Medical Center, USA*

#### *Reviewed by:*

*Arun Bokde, Trinity College Dublin, Ireland Xiaoshan Wang, Nanjing Brain Hospital, China*

#### *\*Correspondence:*

*Robert W. Thatcher, EEG and NeuroImaging Laboratory, Applied Neuroscience Research Institute, 7985 113th Street, Suite 210, Seminole, FL 33772, USA e-mail<EMAIL_ADDRESS>
**Objectives:** The purpose of this study was to explore phase reset of 3-dimensional current sources in Brodmann areas located in the human default mode network (DMN) using Low Resolution Electromagnetic Tomography (LORETA) of the human electroencephalogram (EEG).

**Methods:** The EEG was recorded from 19 scalp locations from 70 healthy normal subjects ranging in age from 13 to 20 years. A time point by time point computation of LORETA current sources were computed for 14 Brodmann areas comprising the DMN in the delta frequency band. The Hilbert transform of the LORETA time series was used to compute the instantaneous phase differences between all pairs of Brodmann areas. Phase shift and lock durations were calculated based on the 1st and 2nd derivatives of the time series of phase differences.

**Results:** Phase shift duration exhibited three discrete modes at approximately: (1) 25 ms, (2) 50 ms, and (3) 65 ms. Phase lock duration present primarily at: (1) 300–350 ms and (2) 350–450 ms. Phase shift and lock durations were inversely related and exhibited an exponential change with distance between Brodmann areas.

**Conclusions:** The results are explained by local neural packing density of network hubs and an exponential decrease in connections with distance from a hub. The results are consistent with a discrete temporal model of brain function where anatomical hubs behave like a "shutter" that opens and closes at specific durations as nodes of a network giving rise to temporarily phase locked clusters of neurons for specific durations.

**Keywords: LORETA, EEG phase reset, phase lock, phase shift, chaos, stability, self-organized criticality**

#### **INTRODUCTION**

When one is at rest and not engaged in a task and absorbed in a ruminating self-narrative about the past and future then it is during these reflective moments that the default mode network (DMN) is activated and the attention network is anti-correlated or reciprocally deactivated (Raichle et al., 2001; Raichle, 2010). The insula appears to function as a switch that is correlated with phase shifting of the attention and default networks activation vs. suppression (Bressler and Menon, 2010). Petersen and Posner (2012) review the functional MRI (fMRI) studies of the attention network and the DMN in attention deficit disorders that are characterized by the intrusion of the self-narrative in academic situations resulting in poor grades. The reciprocal relationship between the DMN related to an ongoing internal self-narrative and the attention network focused on the external world is an important dynamic, however, fMRI has a limited temporal resolution and is unable to resolve millisecond periods of phase lock and phase shift of neurons located in network nodes and functional connections that comprise the DMN.

The EEG has slightly less spatial resolution than fMRI, but adequate spatial resolution to measure the average current density of Brodmann areas at 2 cm3 to about 3 cm3 volumes in the millisecond time domain (Pascual-Marqui, 1999; Vitacco et al., 2002; Mulert et al., 2004; Grech et al., 2011).

There are changes in the synaptic synchrony of millions of neurons connected at varying time delays and frequencies in the DMN. The "DMN" is constituted primarily by the cingulate gyrus, hippocampus, medial frontal lobes, temporal lobes and parietal lobes with approximately five times the number of synaptic connections than any other cortical network (Buckner et al., 2008; Hagmann et al., 2008). Activation of the DMN significantly increases demand on blood glucose and oxygen as well as changes in the synchrony of synaptic potentials on the dendrites and cell bodies of cortical pyramidal neurons as measured in the human EEG using 3-dimensional electrical neuroimaging methods (Pascual-Marqui et al., 1994; Pascual-Marqui, 1999; Michel et al., 2009) also referred to as EEG Tomography (tEEG) (Cannon et al., 2009; Thatcher, 2011; Thatcher et al., 2011) or Brain Electromagnetic Tomography (BET) (Valdés-Sosa et al., 1992; Bosch-Bayard et al., 2001; Hernandez-Gonzalez et al., 2011).

Hughes and Crunelli (2007) and Buzsaki (2006) review how action potentials occur when neurons are in-phase with respect to the local field potentials (LFPs) and how action potentials are suppressed when neurons are shifted anti-phase with respect to the LFP where phase shift is a high speed switch for large collections of neurons to functionally synchronize with sub-sets of neurons in different network nodes in the millisecond time domain. The human EEG is the summation of LFPs arising from pyramidal neuron synapses. The process of phase shift and phase lock of the EEG between different scalp locations has been shown to be correlated with a range of functional conditions. For example, measures of EEG phase reset (PR) have been correlated to various frequency bands during cognitive tasks (Tesche and Karhu, 2000; Kirschfeld, 2005; Kahana, 2006), working memory (John, 1968; Damasio, 1989; Tallon-Baudry et al., 2001; Rizzuto et al., 2003), sensory-motor interactions (Vaadia et al., 1995; Roelfsema et al., 1997), hippocampal long-term potentiation (McCartney et al., 2004), brain development (Thatcher et al., 2008b), autism (Thatcher et al., 2009b), and consciousness (Varela et al., 2001; John, 2002, 2005; Cosmelli et al., 2004).

These studies indicate that measures of phase shift duration and phase lock duration between groups of neurons measured at the scalp EEG must also be capable of being measured at the level of the 3-dimensional sources of the EEG using inverse methods. The initiation of in-phase to anti-phase dynamics of neurons is modeled by long distant excitatory inputs on dendrites where phase shift duration is dependent on the density of neurons in local loops and phase lock duration is determined by the reaction to the short duration excitatory inputs producing long duration inhibitory synaptic potentials (van Drongelen et al., 2004; Ko and Ermentrout, 2007; Tiesinga and Sejnowski, 2010; Li and Zhou, 2011). These studies suggest that phase shift duration is directly proportional to the temporal compactness or density of activated synapses, because, bursting of in-phase action potentials results in an average synaptic potential shift in the frequency of neurons and consequently phase shift duration is expected to be inversely related to neural density, that is, the higher the packing density then the shorter the phase shift duration. This is like local neighbors quickly communicating whereas long distance cousins take more time and effort to synchronize (i.e., long phase shift durations).

A series of EEG PR studies by Lehmann et al. (2006) and Thatcher et al. (2008a, 2009a,b) are consistent with physiological models of EEG PR and have added to earlier studies by Varela (1995), Breakspear (2002, 2004), Freeman (2003), and Freeman et al. (2003) by measuring discontinuities of electrical potentials and current sources of the two main physiological processes that underlie PR, namely, phase shift followed by phase lock. Lehmann et al. (2006) and Thatcher et al. (2007) demonstrated temporal discontinuities of EEG current sources of about 40–250 ms. Thatcher et al. (2009a) studied the development of scalp electrode distance and PR times by measuring phase shift durations (range of about 30–70 ms) and phase lock durations (100–800 ms) from birth to 16 years of age where short distance inter-electrode pairing (6 cm) exhibited shorter phase shift duration and longer phase lock duration than longer distance inter-electrode parings (18–24 cm). Furthermore, it has been shown that phase shift duration is positively related to intelligence while phase lock duration is negatively related to intelligence measured by WISC-R I.Q. test (Thatcher et al., 2008a). The findings of Thatcher et al. (2008a, 2009a) and Lehmann et al. (2006) are consistent with the hypothesis that phase shift is a process involved in the recruitment of available neurons at a given moment of time and phase lock duration is the binding or synchrony of groups of neurons that simultaneously mediate different functions in different brain regions (i.e., sustained commitment of neurons). It is also possible that phase lock reflects the inhibition of billions of "irrelevant" neurons that are excluded or restricted resulting in the "protection" of a small subset of neural loops that mediate a multidimensional sub-network. That is, the large spatial inhibition isolates a spatially smaller subset of synchronized neurons that are masked or invisible to the scalp recorded EEG. The present study is a further exploration of the scalp surface EEG studies of phase shift and phase lock duration by applying 4-dimensional neuroimaging (tEEG) of current sources using Low Resolution Electromagnetic Tomography (LORETA) using the time series of the center voxel of each of 88 Brodmann areas that comprise the DMN (Pascual-Marqui et al., 1994; Pascual-Marqui, 1999; Lehmann et al., 2006; Canuet et al., 2011; Langer et al., 2011). The present study is designed to explore the nature of sudden phase shifts followed by phase lock in small 3-dimensional volumes of EEG current density located in the center of Brodmann areas that constitute the DMN (Buckner et al., 2008). Because of the large number of possible network combinations and frequencies, we limited this study to the delta frequency band (1–4 Hz) and the DMN. Analyses of different frequency bands and locations show similar basic time domain measures with evidence of spatial-frequency "preferences." These analyses will be published in the future.

#### **METHODS**

#### **SUBJECTS**

A total of 70 subjects ranging in age from 13.01 to 19.98 years (males = 41) were included in this study. The subjects in the study were recruited using newspaper advertisements in rural and urban Maryland (Thatcher et al., 1987, 2003, 2007). The inclusion/exclusion criteria were no history of neurological disorders such as epilepsy, head injuries and reported normal development and successful school performance. None of the subjects had taken medication of any kind at least 24 h before testing. All of the subjects were within the normal range of intelligence as measured by the WISC-R and were performing at grade level in reading, spelling and arithmetic as measured by the WRAT and none were classified as learning disabled nor were any of the school aged children in special education classes. All subjects were given an eight-item "laterality" test consisting of three tasks to determine eye dominance, two tasks to determine foot dominance, and three tasks to determine hand dominance. Scores ranged from −8 (representing strong sinistral preference or left handedness), to +8 (representing strong dextral preference or right handedness). Dextral dominant children were defined as having a laterality score of ≥2 and sinistral dominant children were defined as having a laterality score of ≤ −2. Only 9% of the subjects had laterality scores ≤ −2 and 87% of the subjects had laterality scores ≥2 and thus the majority of subjects were right side dominant.

#### **EEG RECORDING**

The EEG was recorded from 19 scalp locations based on the International 10/20 system of electrode placement, using linked ears as a reference. Eye movement electrodes were applied to the inner and outer canthus to monitor artifact and all EEG records were visually inspected and manually edited to remove any visible artifact. Two 5 min of EEG was recorded in the eyes closed and in the eyes open condition. The order of recording for the eyes open followed by closed conditions and vice versa was counter-balanced across subjects. Each EEG record was plotted and visually examined and split-half reliability and test– retest reliability measures of the artifacted data were computed using the Neuroguide software program (NeuroGuide, v2.6.9). The amplifier bandwidths were nominally 1.0–30 Hz, the outputs being 3 db down at these frequencies and the EEG was digitized at 100 Hz. Analyses were performed on 58 s to 2 min 17 s segments of EEG. Split-half reliability tests were conducted on the edited EEG segments and only records with *>*90% reliability were entered into the spectral analyses. Phase shift and lock duration were computed only on contiguous EEG segments.

#### **LORETA TIME DOMAIN COMPUTATION**

The standard procedures for the computation of LORETA were followed according to Pascual-Marqui et al. (1994, 2001), Pascual-Marqui (1999), Gomez and Thatcher (2001), Thatcher et al. (2005a,b). Numerous studies have demonstrated that 19 scalp electrodes are sufficient in number to measure intracranial sources including from the hippocampus (see the listing of 795 publications at: http://www*.*uzh*.*ch/keyinst/NewLORETA/ QuoteLORETA/PapersThatQuoteLORETA05*.*htm).

The Talairach Atlas coordinates of the Montreal Neurological Institute's MRI average of 305 brains (Pascual-Marqui, 1999; Lancaster et al., 2000) was used and the linkage to the standard anatomical 7 × 7 × 7 mm voxels each with a distinct Talairach Atlas Coordinate. Groups of voxels are also defined by the clear anatomical landmarks established by Brodmann in 1909 and referred to as Brodmann areas. The time series of current source vectors in the *x*, *y*, and *z* directions were computed at the center voxel of each of 14 Brodmann area that comprises the DMN. In addition to the x, y, and z time series from each voxel the resultant vector was computed as the square root of the sum of the squares for the *x*, *y*, and *z* source moments.

#### **HILBERT TRANSFORM AND COMPLEX DEMODULATION**

The Hilbert transform of the LORETA time series was computed using complex demodulation to compute instantaneous coherence and phase-differences between each pair of the Brodmann area time series with Talaraich atlas coordinates described in **Table 1** (Granger and Hatanka, 1964; Otnes and Enochson, 1978; Bloomfield, 2000). A total of 91 pairs of the LORETA Brodmann area time series were used to compute "instantaneous" phase differences. This method is an analytic linear shift-invariant transform that first multiplies a time series by the complex function of a sine and cosine at a specific center frequency (Center frequency = 2.5 Hz) followed by a low pass filter (6th order low-pass Butterworth, bandwidth = 1–4 Hz) which removes all but very low frequencies (shifts frequency to 0) and transforms the time series into instantaneous amplitude and phase and an "instantaneous" spectrum (Bloomfield, 2000). We place quotations around the term "instantaneous" to emphasize that there is always a trade-off between time resolution and frequency resolution. The broader the band width the higher the time resolution but the **Table 1 | List of the Brodmann areas of the Default Mode Network (DMN) based on Buckner et al. (2008).**

*The Key Institute's Talairach atlas coordinates in the x, y, and z directions from Lancaster et al. (2000) and Pascual-Marqui (2004). These coordinate values were used to calculate the Euclidean Distances as applied to Equation (1) and in Figure 5.*

lower the frequency resolution and vice versa. Mathematical details are in Thatcher et al. (2008a).

#### **COMPUTATION OF THE 1ST AND 2ND DERIVATIVES OF THE TIME SERIES OF PHASE DIFFERENCES**

The 1st and 2nd derivatives of the time series of instantaneous phase-differences was calculated for all pair wise combinations of DMN voxels in the *x*, *y*, and *z* direction in order to detect instantaneous advancements and reductions of phase-differences. The same mathematical procedures published for measuring phase shift and phase lock duration of the scalp surface EEG were used for the computation of PR of the LORETA time series (Thatcher et al., 2008a, 2009a,b).

PR is composed of two events: (1) a phase shift of a finite duration (SD) and (2) followed by an extended period of phase locking as measured by the phase lock duration (LD) and PR = SD + LD. Phase Shift duration (SD) is the interval of time from the onset of a phase shift to the termination of phase shift (5◦ threshold) where the termination is defined by two conditions: (1) a peak in the 1st derivative (i.e., 1st derivative changes sign from positive to zero to negative) and (2) a peak in the 2nd derivative or inflection on the declining side of the time series of first derivatives. The peak of the 2nd derivative marked the end of the phase shift period. Phase shift duration is the difference in time between phase shift onset and phase shift offset or *SD*(*t*) = *S*(*t*)onset − *S*(*t*)offset. Phase lock duration (LD) was defined as the interval of time between the end of a significant phase shift (i.e., peak of the 2nd derivative) and the beginning of a subsequent

significant phase shift, i.e., marked by the peak of the 2nd derivative and the presence of a peak in the 1st derivative or LD(*t*) = *S*(*t*)offset − *S*(*t*)onset. In summary, two measures of phase dynamics were computed: (1) Phase shift duration (ms) (SD) and, (2) Phase lock duration (ms) (LD). Given the range of epoch sizes from 58 s to 2 min the range of PRs per subject was 66 to 457. **Figure 1** shows an example of the computation of PR metrics in a single subject.

#### **DEFAULT MODE NETWORK (DMN)**

Complex Demodulation was used to compute the Hilbert transform of the current source density time series from 14 left and 14 right hemisphere Brodmann areas were selected based on the review of the human DMN by Buckner et al. (2008). The mathematical description of the equivalence of complex demodulation and the Hilbert transform is by Bloomfield (2000). **Table 1** shows a listing of the nearest Brodmann area fit to the center of the Brodmann areas comprising the DMN to the LORETA Talairach Atlas coordinates (Pascual-Marqui, 1999; Lancaster et al., 2000) and the linkage to the standard anatomical 7 × 7 × 7 mm voxels in the approximate center of each Brodmann areas. All pair wise combinations of the 14 DMN Brodmann areas produced 91 pairs from the left and 91 pairs from the right hemisphere. Only intra-hemisphere analyses of the delta frequency band (1–4 Hz) were included in this study. Different frequencies and cross-hemisphere combinations will be analyzed in a future study. **Table 1** shows the Brodmann areas and Talairach Atlas coordinates of the voxels in the present study which are used to calculate distances between Brodmann area center voxels including the computation of the Euclidean distance between Brodmann areas: *D* = - *x*2 + *y*2 + *z*2.

**Figure 2** shows the LORETA saggital, coronal and horizontal sections of the DMN voxels used in this study.

#### **RESULTS**

#### **VARIMAX FACTOR ANALYSIS**

A principal components analysis followed by a varimax rotation were computed for the 91 left and 91 right intra-hemisphere Brodmann area combinations that comprised the DMN for both phase shift and phase lock duration. A varimax rotation was used because it uses a min/max method to maximize loadings on a given component. **Table 2** shows that with an eigenvalue cutoff of 1.0 that phase shift involved more factors than did phase lock duration and that the same number of factors were involved for the left and right hemispheres. The results of the factor analysis also demonstrated well ordered orthogonality or independent clusters of phase shift and phase lock duration where the variable loadings were essentially the same in the *x*, *y*, *z* directions and for the resultant vector within each hemisphere.

#### **DISCRETE TEMPORAL FRAMES OR "QUANTA" OF PHASE SHIFT DURATION**

Examination of the variables in the factor analyses showed that there are three main "modes" or time frames of duration for

**voxels (Lancaster et al., 2000; Pascual-Marqui, 2004).**

**Table 2 | Number of factors that had Eigenvalues of 1.0 or greater for LORETA phase reset in 91 combinations of Brodmann areas for phase shift and lock in the 3-dimensions (***x, y, z***) and also the resultant vector** *R* **= -** *x***2 +** *y***2 +** *z***2 in the left and right hemispheres.**

phase shift duration. **Figure 3** shows phase shift duration in the x-axis and the percentage of subjects in the study showing specific phase shift durations on the y-axis. The *x*, *y*, and *z* directions all showed essentially the same phase shift duration modes. The three duration modes or "quanta" for phase shift using the resultant vector for both the eyes closed and eyes open conditions demonstrated that discrete durations are present and that there is no or minimal overlap between modes. It can be seen in **Figure 3** that the subjects clustered in three different and discrete time frames or modes. Mode 1 showed a peak phase shift duration of approximately 25 ms, Mode 2 was approximately 50 ms and Mode 3 was approximately 65 ms. There was little overlap between phase duration modes and all variables exhibited only one of the three modes which demonstrates discrete time frames of phase shift duration in specific groups of Brodmann areas that comprise the DMN.

#### **DISCRETE TEMPORAL FRAMES OR "QUANTA" OF PHASE LOCK DURATION**

Examination of the factor analyses showed that there are two main "modes" or time frames of duration for phase lock duration. The *x*, *y*, and *z* directions all showed essentially the same duration modes and **Figure 4** shows the two duration frames for phase lock using the resultant vector for both the eyes closed and eyes open conditions. The x-axis is phase lock duration in milliseconds and the y-axis are the percent of subjects (*N* = 70) in both the eyes closed and open conditions. It can be seen in **Figure 4** that the subjects clustered in two different and discrete time frames or modes. Mode 1 showed a peak phase lock duration of approximately 250 ms and Mode 2 was approximately 425 ms. There was a minor mode at approximately 800 ms which was much smaller than modes 1 and 2. There was little overlap between

phase duration modes and all variables exhibited only one of the two modes with no bi-modal distributions which demonstrates discrete time frames of phase lock duration in specific groups of Brodmann areas that comprise the DMN.

#### **SHORT vs. LONG DISTANCE CONNECTIONS AND LORETA PHASE RESET**

The finding of discrete time frames or modes was explored further by correlating the Euclidean distance of the Talarich atlas coordinates of the DMN (i.e., square root of the sum of squares of *x* + *y* + *z*) with phase shift and lock duration. **Figure 5** are the results of regression analyses in which statistically significant inverse relationships were present between phase shift vs. phase lock duration. That is, phase shift duration tended to be short when Brodmann areas were closer together and lengthened as the distance between Brodmann areas increased. In contrast, phase lock duration tended to be long when Brodmann areas were closer together and short when Brodmann areas were more distant.

These analyses showed that both a non-linear relationship involving "modes" or discrete time frames were present if a small number of Brodmann areas are compared as well as a relatively smooth exponential decline as a function of distance between Brodmann areas when all 91 Brodmann areas are examined. In order to further explore both the general and discrete aspect of the findings we evaluated each combination of Brodmann areas by fitting an exponential equation to the data points. Exponentials provided the best fit of the data, for example, the red line is the evaluation of the equation:

$$T = b\_1 + e^{b\_2 + \binom{b\_3}{d}}\tag{1}$$

where *T* = duration (ms), *d* = distance between Brodmann areas (millimeters) and *b*1, *b*2, and *b*3 are coefficients. The coefficients were derived using Data Desk statistical package (Velleman, 1997). The regression to Equation (1) was statistically significant but in opposite directions for phase shift vs. phase lock duration. Brodmann area pairs that exhibited short phase shift durations also exhibited long phase lock durations while Brodmann areas that exhibited long phase shift durations exhibited short phase lock durations as shown in **Figure 5**. The coefficients of the equation are in **Table 3**. The coefficients in **Table 3** allow one to test the hypotheses in this study by simply entering Talairach Euclidian distances between voxels of the brain into an Excel worksheet and then compute the predicted phase lock and phase shift durations. Shifts in the intercept *b*1 will result in shifts in the starting point of the function.

**LORETA time series directions and the resultant vector in the lower right where** *R* **= -** *x***2 +** *y***+***z***2.** The x-axis is phase lock duration in milliseconds and the y-axis is the percent of subjects that exhibited a given phase lock duration for different Brodmann area pairs. The solid line is the eyes closed condition and the dashed line is the eyes open condition. All of the subjects are

represented within each curve. For example, 100% of the subjects exhibited a phase shift duration between 250 and 500 ms for Brodmann areas 8 and 10 (upper left panel x-lock) and similarly for each Brodmann area pair. The finding of discrete phase lock durations with none or little overlap of data points under each phase lock duration curve was a dominant feature of phase lock duration and demonstrates discrete "temporal quanta."

These data show that nearest neighbors exhibit short phase shift and long phase lock durations while long distance relations exhibit long phase shift and short phase lock durations. This is analogous to every day life where one quickly engages and maintains communication with a local neighbor but it takes longer to engage and remain connected to a long distance neighbor.

#### **EYES OPEN vs. EYES CLOSED CONDITIONS AND LEFT vs. RIGHT HEMISPHERE**

The eyes open condition was an independent replication of the eyes closed condition because the eye conditions were recorded at different times and the order of recording was counter balanced across subjects. Multivariate analysis of variance (MANOVA) with Bonferroni correction was performed where eyes open and eyes closed conditions and left vs. right hemisphere as factors for all 91 Brodmann area combinations. The MANOVA was used because it controls for the intercorrelation between the dependent and between the independent variables. **Table 4** shows the results of these analyses which showed that eyes open vs. eyes closed conditions were not statistically significantly different in the *x*, *y*, *z* directions and resultant vector for phase shift duration. However, eyes open vs. closed conditions were statistically significant for phase lock duration in the x and z directions and for the resultant vector. In all cases, phase lock duration in the eyes closed condition was longer than in the eyes open condition (mean difference = 9.52 ms).

As seen in **Table 4**, there were no statistically significant differences between left vs. right hemisphere in phase shift duration, however, there were statistically significant differences in phase lock duration in the x and y directions. In all cases, phase lock duration was longer in the right hemisphere than in the left hemisphere (mean = 11.05 ms).

#### **VOLUME CONDUCTION TEST**

Finally, the role of volume conduction is important to evaluate. Volume conduction is the propagation of an electomagnetic field at the speed of light or 3 × 1010 cm/s (Feynman et al., 1963; Malmivuo and Plonsey, 1995). For the distance of the order of centimeters the delay is approximately 3.3 × 10−9 s which is a delay that is extremely small and therefore approximates zero phase difference at all points in a volume. Phase

The y-axis is phase shift duration **(Top)** and phase lock duration **(Bottom)** of the Brodmann areas described in **Table 1**. The left row are the left hemisphere Brodmann areas and the right row are the right hemisphere Brodmann areas. The red line is the fit of an exponential equation

correlation and *p* = statistical probability. Phase shift and phase lock are inversely related where Brodmann areas with short phase shift duration exhibit long phase lock durations while Brodmann areas with short phase lock durations exhibit long phase shift durations.

difference is measured by the "imaginary number" component of the cross-spectrum and PR is the 1st derivative of the "imaginary number" that by definition is not volume conduction if a phase difference is greater than zero. A test of the possible influence of volume conduction involved measuring the absolute phase difference between Brodmann areas and then determining if the phase difference was greater than zero and how phase difference changes with distance. If greater than zero then the phase lock measures cannot be explained by volume conduction. Further, because connection density decreases as a function of distance and conduction delays linearly accumulate with distance then if there is an increase in phase differences as a function of distance than this also cannot be explained by volume conduction. If absolute phase differences are significantly greater than zero and increase as a function of distance than this is consistent with physiological connections between connected parts of a network in which conduction velocity, synaptic rise times and synaptic delays produce accumulative time delays.

**Figure 6** shows the results of this test in which mean phase differences were not equal to zero and instead varied as a function of distance which is what is expected in a network of connections and cannot be explained by volume conduction. In fact, phase difference values were many times greater than zero and even the most close or nearby neighbor Brodmann areas 8 and 9 differed by 2–4◦ and therefore also cannot be explained solely by volume conduction. Further, phase lock duration is maximal at short distances and declines with distance that also cannot be explained by volume conduction because a sustained phase difference over time must be maintained which is consistent with the known physiology of local and distant loop iteration that can extend for many milliseconds (Beste and Dinse, 2013).

#### **DISCUSSION**