# 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: researchtopics@frontiersin.org

## 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 <sup>5</sup> , Jing Xiang<sup>6</sup> and William C. Gaetz <sup>7</sup>

*<sup>1</sup> Department of Psychiatry, Osaka University Graduate School of Medicine, Suita, Japan, <sup>2</sup> Neuroscience Center, Ashiya Municipal Hospital, Ashiya, Japan, <sup>3</sup> The KEY Institute for Brain-Mind Research, University of Zurich, Zurich, Switzerland, <sup>4</sup> Neuropsychiatry, Kansai Medical University, Moriguchi, Japan, <sup>5</sup> Department of Cognitive, Social and Organizational Psychology, La Laguna University, Tenerife, Spain, <sup>6</sup> Cincinnati Children's Hospital Medical Center, Cincinnati, OH, United States, <sup>7</sup> 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 ishii@psy.med.osaka-u.ac.jp*

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*

*<sup>1</sup> The KEY Institute for Brain-Mind Research, University of Zurich, Zurich, Switzerland*

*<sup>3</sup> CIMFAV, Universidad de Valparaiso, Valparaiso, Chile*

*<sup>4</sup> Neuroinformatic Department, Cuban Neuroscience Center, Havana, Cuba*

*<sup>5</sup> Department of Psychiatry, Shiga University of Medical Science, Shiga, Japan*

*<sup>6</sup> 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: pascualmarqui@key.uzh.ch*

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

*<sup>2</sup> 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 *<sup>p</sup>* <sup>≥</sup> 1, for *<sup>q</sup>* <sup>≥</sup> 2 time series **<sup>X</sup>**(*t*) <sup>∈</sup> <sup>R</sup>*<sup>q</sup>* <sup>×</sup> 1, is written as:

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

where **<sup>A</sup>**(*k*) <sup>∈</sup> <sup>R</sup>*q*×*<sup>q</sup>* are the autoregressive coefficients, **<sup>ε</sup>**(*t*) <sup>∈</sup> R*q*×<sup>1</sup> 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 *<sup>X</sup>* (ω) <sup>∈</sup> <sup>C</sup>*q*×1, *<sup>A</sup>* (ω) <sup>∈</sup> <sup>C</sup>*q*×*q*,ε(ω) <sup>∈</sup> <sup>C</sup>*q*×<sup>1</sup> 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 **<sup>S</sup>**<sup>ε</sup> <sup>∈</sup> <sup>R</sup>*q*×*<sup>q</sup>* 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**<sup>ε</sup> 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 *<sup>t</sup>* <sup>∈</sup> <sup>R</sup>*NE* <sup>×</sup> <sup>1</sup> denotes the instantaneous scalp electric potential at *NE* electrodes, **<sup>J</sup>***<sup>t</sup>* <sup>∈</sup> <sup>R</sup>(3*NV* )×<sup>1</sup> is the instantaneous current density vector field at *NV* cortical voxels (consisting of three components at each voxel), and **<sup>K</sup>** <sup>∈</sup> <sup>R</sup>*NE* <sup>×</sup> (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 **<sup>T</sup>** <sup>∈</sup> <sup>R</sup>(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***<sup>t</sup>* 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. We emphasize that this EEG was corrupted with relatively large amounts of additive biological and measurement noise. eLORETA was computed at all 6239 cortical voxels. 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 <sup>1</sup> , Ryouhei Ishii <sup>1</sup>\*, Roberto D. Pascual-Marqui 2,3 , Leonides Canuet <sup>4</sup> , Shunichiro Ikeda<sup>1</sup> , Masahiro Hata<sup>1</sup> , Kaoru Imajo<sup>5</sup> , Haruyasu Matsuzaki <sup>6</sup> , Toshimitsu Musha<sup>6</sup> , Takashi Asada<sup>7</sup> , Masao Iwase<sup>1</sup> and Masatoshi Takeda<sup>1</sup>**


#### **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: ishii@psy.med.osaka-u.ac.jp 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: rwthatcher2@yahoo.com*

**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 cm<sup>3</sup> to about 3 cm<sup>3</sup> 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*<sup>2</sup> + *y*<sup>2</sup> + *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** *<sup>R</sup>* **<sup>=</sup> -** *x***<sup>2</sup> +** *y***<sup>2</sup> +** *z***<sup>2</sup> 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*<sup>3</sup> 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*<sup>1</sup> will result in shifts in the starting point of the function.

**LORETA time series directions and the resultant vector in the lower right where** *<sup>R</sup>* **<sup>=</sup> -** *x***<sup>2</sup> +** *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 <sup>×</sup> <sup>10</sup><sup>10</sup> cm/s (Feynman et al., 1963; Malmivuo and Plonsey, 1995). For the distance of the order of centimeters the delay is approximately 3.3 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>s</sup> 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**

This study extends the investigation of the spatial and temporal properties of the scalp surface EEG PR to the current sources of the EEG. There was an approximate seven year age range of

**Table 3 | Coefficients of equation (1) (***b***1,** *b***2, and** *b***3) for the resultant vector for phase lock and phase shift durations in the left and right hemispheres used in Figure 5.**


**Table 4 | The results of the MANOVA for eyes closed vs. eyes open conditions and hemisphere for phase lock and phase shift durations.**


subjects in this study, however, no significant correlations with age were present. A new finding is discrete temporal modes of phase shift and phase lock duration that are unique for different Brodmann area associations of the DMN. A second finding is significantly greater complexity and higher dimensionality of phase shift duration in comparison to phase lock duration, i.e., a temporal dimension reduction from phase shift to phase lock duration. A third finding is an inverse relationship between spatial distance between Brodmann areas and PR metrics where short distances are exponentially related to long phase lock durations and short phase shift duration, whereas, long distances are exponentially related to long phase shift durations and short phase lock duration. Overall, the findings are consistent with previous studies demonstrating that temporal "frames" or "chunks" are discrete and different for different Brodmann areas of the human brain. The results can be explained by a local vs. distant connection model which is unique for different Brodmann areas and thereby accounting for discrete temporal quanta while also exhibiting an exponential decrease in the number of connections with distance that corresponds to exponential changes in PR as a function of distance. Furthermore, the data indicate that discrete peaks in PR duration reflect a different or unique ratio of long vs. short distance connection densities within different Brodmann areas.

Studies by Ko and Ermentrout (2007) and Tiesinga and Sejnowski (2010) use mathematical models that best explain the onset of a phase shift to be initiated by cortico-cortical long distant excitatory dendritic synapses. The mathematical models and experiments fit the physiological facts that local distance connections are exponentially more numerous and exhibit temporal compactness while distant connections are less numerous and exhibit temporal dispersion. Therefore, the ratio of local to distant connections determines the duration of phase shift and this ratio varies as an exponential function of distance between a given Brodmann area. This conclusion is consistent with studies showing anti-phase shifts are related to neural packing density (Li and Zhou, 2011) and with the surface EEG PR where increased local packing density was hypothesized to explain the difference between PR in the posterior-to-anterior direction (e.g., O1-P3) vs. the anterior-to-posterior directions (e.g., Fp1-F3) because of

increased packing density in the posterior cortex as compared to frontal cortex (Thatcher et al., 2009a). The results of the present study are also consistent with the hypothesis that phase shift is related to recruitment of neurons and phase lock duration is directly related to synchrony or binding of neurons as an exponential function of packing density and inter-node (Brodmann area) distance.

#### **LIMITATIONS OF THIS STUDY**

One limitation is the sample rate used where 100 Hz sample rates limit analyses to 10 ms resolution and therefore phase shift durations at shorter durations may be present but were missed. However, the time resolutions appears to be due more to physiology than the speed of computer measurement because we have measured phase shift and lock duration at sample rates of 512 Hz (1.95 ms resolution) and find the same minimum phase shift durations of about 20 ms as reported in this study. Mathematically the limiting relationship is between distance and time as described by Equation (1) that can be evaluated using the coefficients in **Table 3**. If *d* = 0 then duration = infinity. If *d* = infinity then duration = 0. The values between zero and infinity are a good fit of how packing density of neurons as a function of distance are related to PR duration. Higher sample rates should be used to further test the temporal synchrony limits as reported in this study. Another limitation are analyses only of the delta frequency band. Because of the large volume of data and limitations of publication space it was necessary to limit the analyses to a single frequency band, i.e., the delta frequency band (1–4 Hz). Different relations are present at different frequency bands as well as cross-frequency band coupling and these relations have been measured and will be the subject of future publications. Another limitation is the use of a center frequency of a narrow band that is necessary with the Hilbert transform. The Gabor transform is independent of a center frequency and provides optimal time frequency resolution (Witte and Schack, 2003) and higher temporal-spatial resolution and can yield more detailed time-frequency information than Complex Demodulation. While cross-frequency phase shift and phase lock measurements is important, however, because of page limitations this topic will be also presented in future studies.

#### **DISCRETE DURATIONS AND THE DEFAULT MODE NETWORK**

Shallace (1964), Efron (1967, 1970a,b), Allport (1968), and others (Sanford, 1971; Varela, 1995; Varela et al., 2001) have shown a minimum perceptual frame from approximately 40 ms for auditory stimuli to 140 ms for visual stimuli that are durations that temporally distinguish events as being successive in time where duration is defined at T1 – T2 = 0, i.e., simultaneity where there is no perceived time difference between two distinct events. These studies as well as others show that learning-dependent changes in neural networks is not a continuous process but rather a discontinuous sequencing of narrow time windows (Thatcher and John, 1977; John, 2005; Lehmann et al., 2006; Thatcher et al., 2007, 2008b, 2009a). Thatcher et al. (2009a)indicated a linkage between spontaneous and ongoing perceptual frames and event related desynchronization (ERD) by considering phase shift duration and phase lock duration as elemental "atoms" that underlie the duration of perceptual frames and ERD. For example, in Thatcher et al. (2009a) the mode of scalp surface EEG PR was temporally bounded with a minimal phase shift duration of about 45 ms and a maximum phase shift duration of about 70 ms. Phase shift was followed by phase lock that was temporally bounded from about 150 ms to about 800 ms with the most frequent phase locking intervals between 200 and 450 ms (Thatcher et al., 2009a). The findings in the present study are consistent with the earlier surface EEG analyses and indicate that 3-dimensional current source phase shift and phase lock are ongoing spontaneous processes that occur between network nodes and that discrete phase shift durations (i.e., discrete quanta of time) operate at higher temporal-spatial precision than the surface EEG. Network nodes are defined as clusters of neurons connected to other clusters (nodes) and the present study shows that a fundamental property of nodes or clusters is to operate like temporal "shutters" that open and close at specific durations. Each Brodmann area maintains a different ratio of local vs. distant connections but nonetheless follows the general rule of an exponential decrease in local connections as a function of distance. This anatomical fact may explain the findings of discrete phase shift and lock duration due to unique local vs. distant excitatory connection densities in different Brodmann areas (Sporns, 2011).

**Figure 7** illustrates an hypothesis to explain the findings in this study by fitting the data to a single exponential model based on the ratio of local and distant excitatory dendritic synapses. **Table 3** provides the coefficients of Equation (1) to allow one to experiment with different Talariach distances, e.g., enter the *x*, *y*, and *z* Euclidean distance between Hagmann et al. (2008) Modules in millimeters into Equation (1) and then calculate the predicted LORETA phase shift and phase lock duration.

The EEG is the summation of LFPs therefore phase shifts are necessarily related to connection density. However, the arrival of distant synaptic action potentials exhibits temporal dispersion whereas local excitatory connections are higher in number and temporally compact. Therefore, it is hypothesized that the ratio of distant to local connections varies as a function of distance from any Brodmann area and as a consequence longer phase shift duration occurs as the ratio shifts toward distant excitatory inputs. It is hypothesized that phase lock is inversely proportional to phase shift duration based on the spatial-temporal GABA connections and delays between excitatory neurotransmitter EPSPs (and excitatory neuromodulators) in local and long distance loops. It is also hypothesized that the physiological differences in the genesis of phase shift vs. phase lock is related to the "Gap" of time between distant and local EPSP excitation that produce the phase shift followed by long duration inhibitory synaptic potentials that contribute to phase lock duration. The collective resonance of the rebound from local inhibition and the arrival of long distant EPSPs contribute to the onset of the PR process.

Whether or not the phase shift and phase lock are time correlated to a task is irrelevant since "self-organized criticality" is an ongoing background emergent process that on the average produces an approximately 20–80 ms period of phase shift or "uncertainty" or approximate duration of "chaos" followed on the average by a 200–800 ms period of phase locking

**FIGURE 7 | (A)** (Top) is an example of discrete durations or "temporal quanta" of phase shift duration in different Brodmann areas of the DMN. The eyes vs. eyes closed shift to longer durations are due to increased functional connectivity with increased input that is like a "shutter" whose duration is proportional to the number of recruited neurons. The discontinuities are due to different packing densities in different Brodmann areas. The bottom plot **(B)** is the evaluation of Equation (1) using the mean LORETA phase shift and lock duration for both eyes closed and open conditions for the 91 Brodmann area combinations: *<sup>T</sup>* <sup>=</sup> *<sup>b</sup>*<sup>1</sup> <sup>+</sup> *<sup>e</sup>b*2+*(b*3/ *<sup>d</sup>)* The link between EEG phase shift and neural packing density is the physiological observation of action potential bursts when in-phase to LFPs vs. suppression of action potentials when ant-phase to LFPs Hughes and Crunelli (2007). EEG is the summation of LFPs therefore phase shifts are necessarily related to neural packing density, i.e., the higher the packing density than the longer the phase shift as a property of summation. It is hypothesized that the increased shift duration

between eyes closed vs. eyes open is due to increased arousal and increased depolarization resulting in increased functional connectivity (increased neural resource). Phase lock is inversely proportional to phase shift duration based on the spatial-temporal GABA connections and delays between excitatory neurotransmitter EPSPs in local and long distance loops. The physiological differences in the genesis of phase shift vs. phase lock is related to the "Gap" of time which is a transition time between local EPSP excitation and re-enterant long distance EPSPs that produce the phase shift followed by long duration inhibitory synaptic potentials that contribute to phaselock duration. The rebound from inhibition and arrival of EPSPs starts the phase reset process. Predicted phase shift and lock durations can be evaluated by using the coefficients in **Table 3**. For example, the distance between Module1 and Module2 from Hagmann et al. (2008) and Thatcher et al. (2011, **Table 3**) is 43.4 mm which that predicts a phase shift duration = 56 ms and phase lock duration = 300 ms based on Equation (1).

or "stability." This process represent a continuous sequence of meta-stable states (Rabinovich et al., 2012). The loop network background process includes "blank" periods when large assemblies of neurons are in a PR mode (i.e., phase shift and phase locking) and otherwise not available to participate in loops which represents an approximate average 135 ms "Gap" or period of time between when neurons recruited by a phase shift are followed by phase lock onset. This gap is a statistical region with an average and range between the end of a phase shift and the beginning of a subsequent phase lock. This appears to be a fundamental "unconscious" transition time that is near to the flicker fusion frequency and may be why TV viewers rely upon instant replay to confirm a touch down in football or a replay in baseball, etc. (i.e., human time frames are too long and require higher speed instant replay to determine the reality of an event under question). The "gap" interval of about 135 ms is spatially distributed and temporally averaged across billions of loops in the brain and results in a continuous or smooth appearance of reality.

Another type of "refractory" period or "gap" is when phase locked neurons are unavailable for allocation by a different cluster of neurons at a different moment of time. Long distance phase locking of local clusters of neurons can result in a reduction in the amplitude of the surface EEG because phase locking occurs over long distances and thus reduces the size of the number of synchronized local cluster of neurons by spatial differentiation. This is consistent with studies of schizophrenia that show hyperconnectivity in local frontal and parietal regions associated with increased local current density (Canuet et al., 2011). The hypothesis of linking PR during the background spontaneous EEG provides a new definition of the term "desynchronization" used to describe ERD and the waxing and waning of the spontaneous EEG. That is desynchronization is actually "spatially differentiated phase reset" or "micro bonding" of local clusters of neurons connected across long distances to other local clusters for brief periods of time (Thatcher et al., 2009a).

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

The eyes closed vs. eyes open conditions were measured at different times with a pause between recordings and thus are within subject replications. Both conditions exhibited significant fits to the exponential Equation (1) and there were no significant differences between eyes open and eyes closed phase shift duration (see **Table 4**). The eyes closed condition did exhibit a small but significantly longer phase lock duration which may reflect the relatively small differences in neural excitability between these two states. The lack of a large difference between resting eyes closed and eyes open states is consistent with fMRI studies of the DMN that is evident in the resting state (Raichle, 2010). This is because the resting state is characterized by quiet repose with either eyes closed or eyes open and with or without visual fixation. During the resting state subjects typically experience an ongoing state of conscious awareness filled with "stimulus-independent thoughts" (SITs; Antrobus, 1968) or more commonly, day dreaming or mind wandering (Mason et al., 2007). Internal thoughts, self-awareness and rumination commonly dominate the resting state as opposed to active task engagement (Sridharan et al., 2009; Raichle, 2010). The right hemisphere exhibited longer phase lock durations than the left hemisphere (see **Table 4**) which is consistent with the findings of lower packing density in the right hemisphere in comparison to the left, especially in the planum temporale (Gur et al., 1980; Buchell et al., 2004).

#### **ACKNOWLEDGMENTS**

We are indebted to Drs. Rebecca McAlaster, David Cantor and Michael Lester and Ms. Sheila Ignasius and Ms. Diane Pruit for their involvement in the recruitment, EEG testing and evaluation of subjects and Rebecca Walker and Richard Curtin for database management. An IRB approved consent form was signed by the parents and guardians of subjects 18 years and younger and by the subjects older than 18 years in this study. This research was supported by grants from USDA HRD-0200, USDA CSRS 801- 12-09C, and NIH Grant RR-08079-09.

#### **REFERENCES**


Velleman, P. F. (1997). Data Desk, V. 6. Ithaca, NY: Data Description, Inc.


**Conflict of Interest Statement:** The authors are on the board of the company Applied Neuroscience, Inc. However, no company activity was involved in this research study and the data were collected prior to the creation of the company.

*Received: 10 April 2014; accepted: 30 June 2014; published online: 23 July 2014. Citation: Thatcher RW, North DM and Biver CJ (2014) LORETA EEG phase reset of the default mode network. Front. Hum. Neurosci. 8:529. doi: 10.3389/fnhum. 2014.00529*

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

*Copyright © 2014 Thatcher, North and Biver. 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.*

### Alpha band functional connectivity correlates with the performance of brain–machine interfaces to decode real and imagined movements

#### *Hisato Sugata1, Masayuki Hirata1,2 \*,Takufumi Yanagisawa1,2,3 , Morris Shayne1, Kojiro Matsushita1, Tetsu Goto1,2 , Shiro Yorifuji <sup>2</sup> and Toshiki Yoshimine1*

<sup>1</sup> Department of Neurosurgery, Osaka University Medical School, Suita, Japan

<sup>2</sup> Division of Functional Diagnostic Science, Graduate School of Medicine, Osaka University, Suita, Japan

<sup>3</sup> ATR Computational Neuroscience Laboratories, Kyoto, Japan

#### *Edited by:*

Ryouhei Ishii, Osaka University Graduate School of Medicine, Japan

#### *Reviewed by:*

William C. Gaetz, The Children's Hospital of Philadelphia, USA Leonides Canuet, Center for Biomedical Technology, Spain

#### *\*Correspondence:*

Masayuki Hirata, Department of Neurosurgery, Osaka University Medical School, 2-2 E6 Yamadaoka, Suita, Osaka 565-0871, Japan e-mail: mhirata@nsurg.med. osaka-u.ac.jp

Brain signals recorded from the primary motor cortex (M1) are known to serve a significant role in coding the information brain–machine interfaces (BMIs) need to perform real and imagined movements, and also to form several functional networks with motor association areas. However, whether functional networks between M1 and other brain regions, such as these motor association areas, are related to the performance of BMIs is unclear. To examine the relationship between functional connectivity and performance of BMIs, we analyzed the correlation coefficient between performance of neural decoding and functional connectivity over the whole brain using magnetoencephalography. Ten healthy participants were instructed to execute or imagine three simple right upper limb movements. To decode the movement type, we extracted 40 virtual channels in the left M1 via the beam forming approach, and used them as a decoding feature. In addition, seed-based functional connectivities of activities in the alpha band during real and imagined movements were calculated using imaginary coherence. Seed voxels were set as the same virtual channels in M1. After calculating the imaginary coherence in individuals, the correlation coefficient between decoding accuracy and strength of imaginary coherence was calculated over the whole brain. The significant correlations were distributed mainly to motor association areas for both real and imagined movements. These regions largely overlapped with brain regions that had significant connectivity to M1. Our results 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.

**Keywords: brain–machine interfaces, functional connectivity, alpha band, real movement, imagined movement, magnetoencephalography, primary motor area, motor association area**

#### **INTRODUCTION**

The brain signals recorded from the primary motor cortex (M1) are known to serve a significant role in providing the information necessary for brain–machine interfaces (BMIs). This technology is expected to offer patients who have lost control of voluntary movements, including those with amyotrophic lateral sclerosis (ALS) and spinal cord injury, greater independence, and a higher quality of life by enabling them to control external devices to communicate with others and to manipulate their environment at will (Wolpaw et al., 2002; Birbaumer, 2006; Hirata et al., 2012; Hochberg et al., 2012; Collinger et al., 2013). Recently, many studies reported the importance of M1 signals in providing the information necessary for BMIs using various types of signal platforms to execute real and imagined movements; for example, electroencephalography (EEG; Bradberry et al., 2010; Shindo et al., 2011), magnetoencephalography (MEG; Mellinger et al., 2007; Buch et al., 2008; Waldert et al., 2008; Wang et al., 2010; Sugata et al., 2012a), and electrocorticography (ECoG; Leuthardt et al., 2004; Schalk et al., 2007; Yanagisawa et al., 2011, 2012a).

In electrophysiological studies, particular ranges of neural oscillations, which are usually classified into alpha (8–13 Hz), beta (14–25 Hz) and gamma (30–90 Hz), were shown to be associated with motor control (Lopes da Silva, 2013), and their applications to BMIs have been investigated (Wolpaw and McFarland, 2004; Birbaumer, 2006; Yanagisawa et al., 2011). Rhythmic activity in the alpha range observed over the region of the Rolandic fissure is typically not in the form of a sinusoidal curve (Pfurtscheller and Neuper, 1997) and variably referred to as mu rhythm (Gastaut, 1952). It can be observed along with beta band activity during movement (Pfurtscheller and Aranibar, 1977; Cheyne, 2013) and tactile stimulation (Gaetz and Cheyne, 2006). Since modeling of non-sinusoidal waveforms requires the use of higher frequency harmonic components in addition to a fundamental frequency, beta rhythm activity associated with mu rhythm might result from the non-sinusoidal nature of mu rhythms, rather than an independent physiological processes (Jurgens et al., 1995). In addition, the gamma band has been shown to correlate with the firing activities of neurons representing neural information (Ray et al., 2008; Quian Quiroga and Panzeri, 2009). These frequency bands

compose task-specific spatial connectivity patterns in movement related neural networks such as those involving the M1, premotor cortex (PMC), and supplementary motor area (SMA; Herz et al., 2012). Among these frequency bands, recently, functional connectivity within the range of alpha band between the sensorimotor area and motor association area was shown to be relevant to post-stroke recovery potential (Westlake et al., 2012). In this study, neural oscillations of the alpha band were used to calculate functional connectivity because of the higher signal-to-noise ratio compared to oscillations in other frequency bands (e.g., theta, beta) and because they play an important role in controlling cortical excitability. In addition, alpha oscillations are also relevant in the controlling of motor execution through the modulation of gamma-band activity (Yanagisawa et al., 2012b). This alpha oscillation in the sensorimotor cortex, i.e., mu rhythm, has been observed in relation to not only motor execution (Salmelin and Hari, 1994; Leocani et al., 2001) but also motor preparation (Pfurtscheller et al., 1997; Pineda, 2005) and motor imagery (Pfurtscheller et al., 2006; Llanos et al., 2013) as well as beta oscillations, and is considered a mechanism for improving information processing during these tasks (Basar et al., 2001; Palva and Palva, 2007; Sabate et al., 2012). Furthermore, functional connectivities within the range of alpha band activity are suggested to be related to physical and mental fitness (Douw et al., 2014). Such neurophysiological aspects have also been proposed as useful markers of impaired brain states, such as schizophrenia (Hinkley et al., 2011), Alzheimer's disease (Canuet et al., 2012), and multiple sclerosis (Cover et al., 2006). However, in the field of BMIs, there have been few studies focusing on the relationship between functional connectivity within the range of alpha band activities and the performance of BMIs. Based on the above findings, we hypothesized that alpha band activity is a key component to revealing the relationship between the functional connectivity of M1 and performance of BMIs in decoding real and imagined movements. We further hypothesized that brain regions possessing strong alpha functional connectivity with M1 contribute to the performance of BMIs.

The aim of this study was to clarify the relationship between alpha functional connectivity and the performance of BMIs. For this purpose, we used MEG to examine the relationship between the performance of neural decoding, which has been also termed "decoding accuracy" in several studies (Waldert et al., 2008; Bradberry et al., 2010; Yanagisawa et al., 2011), and functional connectivity of activities within the alpha band (8–13 Hz). MEG has several advantages for analyzing functional connectivity compared with EEG and fMRI. MEG has a higher spatial resolution than EEG, and can record a direct correlate of neural activity with high temporal resolution compared with fMRI. We extracted 40 virtual channels in the left M1 using a beam forming approach and used them as a decoding feature. In addition, we calculated seed-based functional connectivity over the whole brain using alpha band activity. Seed voxels corresponded to the same locations as the 40 virtual channels set in the left M1. We then computed the task-related functional connectivity instead of that during the resting state because previous studies using MEG (Bardouille and Boe, 2012) and fMRI (Newton et al., 2007; Treserras et al., 2009) showed that functional

connectivity during motor tasks is greater than that in the resting state. After calculating the task-related functional connectivity, the correlation coefficients between decoding accuracy and strength of functional connectivity were calculated over the whole brain.

#### **MATERIALS AND METHODS PARTICIPANTS**

Ten healthy volunteers participated in this study (five males and five females; mean age 29.8 ± 13.2 years). All participants were confirmed to be right-handed using the Edinburgh Handedness Inventory (Oldfield, 1971; all participants had a score of 100), had no history of neurological or psychiatric diseases, and had normal vision. The protocol of this study was approved by the ethics committee of Osaka University Hospital and all participants provided informed, written consent.

#### **TASKS**

The experimental paradigm is shown in **Figure 1A**. We prepared two tasks: a real movement task and an imagined movement task. We have previously shown the contribution of M1 signals in classifying movement types using these motor tasks based on ECoG (Yanagisawa et al., 2009) and MEG (Sugata et al., 2012b). An epoch started with a 4-s rest phase, and a black fixation cross (+) was presented to fix the participant's eyes on the screen. Then, a Japanese word representing one of the three right upper limb movements (grasping, pinching, or elbow flexion) was presented for 1 s to instruct the participant which movement to perform or imagine after the appearance of the execution cue. Two timing cues, "> <" and "> <," were then sequentially presented for 1 s each to enable the participants to prepare the execution of the real or imagined movements. In the real movement task, the participants were instructed to perform the instructed movement presented on the display immediately after the appearance of the execution cue (×). In the imagined movement task, the participants were instructed to imagine performing the movement immediately after the appearance of the execution cue. Each of the three types of movements was performed 60 times during the real movement trials, and the movement in any given epoch was selected randomly. Then the imagined movement trials were conducted in the same manner.

#### **MEG MEASUREMENTS**

Neuromagnetic activity was recorded in a magnetically shielded room using a 160-channel whole-head MEG system equipped with coaxial type gradiometers (MEG vision NEO;Yokogawa Electric Corporation, Kanazawa, Japan). The participant lay on a bed in the supine position with their head centered. The head position was measured before and after recording using five coils placed on the face (the external meatus of each ear and three points on the forehead). Visual stimuli were displayed on a projection screen positioned 325 mm from the participant's eyes using a visual presentation system (Presentation; Neurobehavioral Systems, Albany, CA, USA) and a liquid crystal projector (LVP-HC6800; Mitsubishi Electric, Tokyo, Japan). Data were sampled at a rate of 1000 Hz with an online low-pass filter at 200 Hz. To reduce contamination from muscle activity and eye movements, we instructed the participants to rest their elbows on a cushion to

Participants performed a real movement task and an imagined movement task following the same task paradigm. Each trial consisted of four phases: a rest phase, an instruction phase, a preparation phase, and an execution phase. In the rest phase, participants fixed their eyes on a black fixation cross "+" presented for 4 s. A Japanese word representing one of three movements was then presented for 1 s during the instruction phase. Then, two timing cues, "> <" and "> <," were presented during the preparation phase to enable the participants to prepare the execution of real or imagined movements. Finally, the

avoid shoulder movements, and to watch the center of the display without ocular movements and blinking. In addition, to monitor unwanted muscular artifacts, electromyograms (EMG) were simultaneously recorded with electrodes on the *flexor pollicis brevis*, *flexor digitorum superficialis*, and *biceps brachii* muscles during the tasks.

movement presented during the instruction phase. Each of the three movements was performed 60 times. **(B)** Analysis procedure. The beam forming approach was used to extract 40 virtual channels from the left M1, and decoding accuracy was calculated using these channels. Seed-based functional connectivity of activities within the alpha band between M1 virtual channels and target voxels over the rest of the whole brain was calculated using imaginary coherence (IC) in individual participants. Then, the correlation coefficient between decoding accuracy and IC was calculated over the participants.

After data acquisition, a 60-Hz notch filter was applied to eliminate the AC line noise, and eye blink artifacts were rejected applying the signal-space projection (SSP), one of the approaches implanted in Brainstorm<sup>1</sup> to reject external disturbances (Tadel

<sup>1</sup>http://neuroimage.usc.edu/brainstorm

et al., 2011). In addition, to align the MEG data with individual MRI data, three-dimensional data of facial surfaces obtained by laser scanning were superimposed on the anatomical facial surface provided by the individual MRI data with an anatomical accuracy <1 mm.

#### **VIRTUAL CHANNELS AND PREPROCESSING**

To extract M1 signals from the MEG sensor, we used an adaptive, spatial filtering beamforming technique (Sekihara et al., 2002). This approach is used to estimate the temporal course of neural activity at a particular site in the brain marked by an imaging voxel, such as that derived from MRI. The output of such a spatial filter is termed a *virtual channel* or virtual sensor (Robinson and Vrba, 1999). The beamformer is constructed to project signals exclusively from the targeted voxels, while removing residual noise to suppress signals from other parts of the brain. Thus, virtual channels provide data concerning neural activity at target voxels with a considerably higher signal-to-noise ratio than that of raw MEG data (Robinson and Vrba, 1999).

The target location of the virtual channels for the present study was the left M1 gyrus. Forty virtual channels were selected in the M1 with an approximately 2.5 mm inter-sensor spacing using the Montreal Neurological Institute (MNI) coordinates (**Figure 1B**). Then, the virtual channel location coordinates on individual MRIs were extracted utilizing MNI coordinates and warping parameters calculated by Statistical Parametric Mapping 8 (SPM8; Wellcome Department of Imaging Neuroscience, London, UK) using an MRI-T1 template and individual MRI-T1 images. A tomographic reconstruction of the data was created by generating a singlesphere head model based on the head shape obtained from the structural MRIs of each individual participant.

Presentation of the execution cue was defined as the onset of real and imagined movements (0 ms), and all time windows were relative to this time. Epochs from –4000 to 500 ms were analyzed. The baseline was set from –4000 to –3500 ms, during the resting phase. Data from each epoch were normalized by subtracting the mean and then dividing by the SD of the baseline values.

#### **FUNCTIONAL CONNECTIVITY ANALYSIS**

The MEG sensor data were reconstructed in source space with the same beamformer approach described above with 5-mm voxel spacing over the whole brain. The frequency component of the alpha band was chosen to calculate source-space, and seed-based functional connectivity. The functional connectivity at 0–500 ms was calculated with imaginary coherence (IC), one of the connectivity analysis approaches that can reduce overestimation biases in EEG/MEG data generated from common references, cross-talk, and volume conduction (Nolte et al., 2004; Guggisberg et al., 2008; Hinkley et al., 2011). IC rules out real parts of coherence containing similarities with zero time lag, and uses imaginary parts of coherence which contains similarities with a certain time lag, because phase similarities with zero time delay among time series are likely to be caused by crosstalk or volume conduction. Using this method, we can evaluate the "true" interactions between brain areas occurring with a certain time lag. Seed voxels were set at the 40 virtual channels in the left M1 at the same locations described above, and the targets were

set as voxels over the remaining whole brain (i.e., except the left M1). The connectivity at each voxel was estimated by averaging across all its Fisher's Z-transformed connections (Nolte et al., 2004; Guggisberg et al., 2008; Hinkley et al., 2011). All ICs calculated from 40 seed voxels were averaged and used as the strength of functional connectivity between M1 and the target voxel.

Group statistical maps were generated to reveal the brain regions with significant ICs during real and imagined movements. The statistical significance of IC across participants was tested with SPM8. The functional images were normalized using the MNI template in SPM8. A one-sample *t*-test at the voxel level was performed using a *t*-statistic incorporating variance smoothing with an 8-mm Gaussian kernel. Voxels with differences at *p* < 0.01 (familywise error rate, FWER) were considered statistically significant, and were superimposed on the template of the inflated cortical surface brain extracted by FreeSurfer2.

#### **DECODING ANALYSES**

Several studies reported that the amplitudes of brain waveforms yield higher performances for BMIs than their power spectrums, such as alpha, beta, and gamma bands (Schalk et al., 2007; Waldert et al., 2008; Yanagisawa et al., 2009). In addition, although the high gamma band activity of ECoG signals is also known to provide high BMI performance (Leuthardt et al., 2004; Yanagisawa et al., 2011), it is difficult for MEG to record high gamma band activity and to obtain high BMI performances. With this in mind, we chose a low frequency component to decode the movement types. The normalized amplitude of the signal recorded at each M1 virtual channel from 0 to 500 ms was resampled over an average 100-ms time window, sliding by 50 ms (9 time points) and then used as a decoding feature. In our preliminary analysis, we also examined other features based on the power spectra of the 40 virtual channels (theta; 4–8 Hz, alpha; 8–13 Hz, beta; 13–25 Hz, low-gamma; 25–50 Hz), but such features did not outperform the normalized amplitudes. Thus, we focused on the decoding results obtained from the low frequency component of the normalized amplitudes.

To examine decoding accuracy, we used a support vector machine (SVM) operating on MATLAB 2013a software (Math-Works, Natick, MA, USA), which was extended to discriminate multiple movements (Kamitani and Tong, 2005; **Figure 1B**). Decoding accuracy was evaluated using 10-fold cross-validation. Each dataset was divided into 10 parts; the classifiers were determined from 90% of the dataset (training set) and tested on the remaining 10% so that the testing dataset was independent from the training dataset for each time point. This procedure was then repeated 10 times. The averaged decoding accuracy over all runs was used as a measure of decoder performance. The binomial test was used to confirm that the decoding performance significantly exceeded chance levels.

#### **CORRELATION BETWEEN IC AND DECODING ACCURACY**

To examine whether functional connectivity is associated with decoding accuracy, we calculated the correlation coefficient

<sup>2</sup>http://freesurfer.net/

between the IC and the decoding accuracy among the ten participants using the Spearman's rank correlation test over the whole brain (**Figure 1B**). All correlation analyses were corrected for multiple comparisons with a false discovery rate (FDR). Voxels with differences at *p* < 0.05 were considered statistically significant, and were superimposed on the template of the inflated cortical surface brain extracted by FreeSurfer.

#### **RESULTS**

#### **FUNCTIONAL CONNECTIVITY DURING REAL AND IMAGINED MOVEMENTS**

During real movements, statistically significant ICs were observed in the bilateral superior and middle frontal gyri, including the SMA and PMC, in the left parietal lobe and the temporal lobe, and in the right sensorimotor area (**Figure 2** left and **Table 1**). During imagined movements, statistically significant ICs were localized only in the left hemisphere, including the left inferior and superior parietal lobules (IPL, SPL), the superior and middle frontal gyri, and the postcentral gyrus (**Figure 2** right and **Table 2**).

#### **DECODING ACCURACY FOR REAL AND IMAGINED MOVEMENTS**

For the real movements, the averaged decoding accuracy among all participants was 67.1 ± 12.5% (mean ± SD), which was significantly higher than chance level (binomial test, *p* < 0.05; **Figure 3**). For the imagined movements, decoding accuracy was also significantly higher than chance level (48.7 ± 8.7%; binomial test, *p* < 0.05), although it was lower than that for the real movements.

#### **CORRELATION OF IC AND DECODING ACCURACY FOR REAL AND IMAGINED MOVEMENTS**

After calculating the ICs in individuals, we examined the correlation coefficient between strength of IC and decoding accuracy among all participants during real and imagined movements. **Figure 4** depicts the distribution of significant correlations between IC and decoding accuracy over the whole brain during real movements (*p* < 0.05, FDR-corrected). Significant correlations were localized mainly to the left PMC, postcentral gyrus, and right sensorimotor area (**Figure 4** upper panel and **Table 3**). On the other hand, significant correlations between IC and decoding accuracy for imagined movements were more widely distributed than those of the real movements (**Figure 4** lower panel). In particular, large clusters were observed in the left IPL and SPL and the right inferior frontal gyrus (IFG). Other significant correlations were observed in the left prefrontal cortex (including dorsolateral prefrontal cortex; DLPFC) and right sensorimotor area (**Figure 4** lower panel and **Table 4**).

**Figure 5** depicts the overlay map of the distribution of significant correlations between strength of IC and decoding accuracy and significant IC during real and imagined movements. The significant correlations were mainly distributed in or around the brain regions that exhibited significant IC during real and imagined movements. No overlap between correlations and IC was observed in the right hemisphere during imagined movements.

#### **DISCUSSION**

To explore the contribution of functional connectivity to the performance of BMIs, we examined the relationship between neural decoding and alpha band IC with the left M1 during real and imagined movements. The brain regions with significant functional connectivity with M1 during both real and imagined movements were distributed mainly in motor association areas, including the SMA, PMC, and parietal area. In addition, the significant correlations between decoding accuracy and IC strength were distributed in or around the brain regions with significant IC. These results indicate that functional connectivity

**FIGURE 2 | Cortical connectivity maps during real and imagined movements.** Group (N = 10 participants) results of anatomically constrained imaginary coherence (IC) visualized on the inflated cortical surface during real and imagined movements. The brain regions with significant IC with the left M1 are represented in blue (p < 0.01, FWER-corrected).


**Table 1 | Brain regions showing significant IC with the left M1 during real movements.**

IC, imaginary coherence; M1, primary motor cortex; MNI, Montreal Neurological Institute template.

**Table 2 | Brain regions showing significant IC with the left M1 during imagined movements.**


IC, imaginary coherence; M1, primary motor cortex; MNI, Montreal Neurological Institute template.

of alpha band activity between M1 and the motor association area is involved in neural decoding of real and imagined movements.

#### **FUNCTIONAL CONNECTIVITY IN THE ALPHA BAND DURING REAL AND IMAGINED MOVEMENTS**

Recent studies have suggested that fluctuations in alpha band oscillations facilitate processing in task-relevant cortical regions or suppress processing distracting input in task-irrelevant regions

to improve task performance (Palva and Palva, 2007; Mazaheri and Jensen, 2010; Gould et al., 2011; Haegens et al., 2011a,b). By holding and releasing high gamma activity during a movement task, we previously demonstrated a functional role of alpha band activity in movement that modulates motor representation in the sensorimotor cortex (Yanagisawa et al., 2012b). Because the main body of alpha band activity recorded over the sensorimotor area is thought to be due to the mu rhythm (Salmelin and Hari, 1994; Leocani et al., 2001), fluctuations in the mu rhythm may play a significant role in controlling the cortical excitability in M1. In fact, it has been demonstrated that data processing improves when the phase of the mu rhythm is modified, and data processing is inhibited when its phase is unlocked (Sabate et al., 2012). Furthermore, the power of the mu rhythm in the sensorimotor area was recently demonstrated to play an important role in cortico-cortical connectivity (Ronnqvist et al., 2013). In the present study, we calculated seed-based functional connectivity from alpha band activity over the whole brain during real and imagined movements. Seed voxels were set as 40 virtual channels in M1. The significant functional connectivity was distributed to the motor association area over frontal and parietal areas during real and imagined movements. The results were mostly concordant with previous studies using fMRI (Solodkin et al., 2004; Gao et al., 2011). Given that these motor association areas have been shown to play important roles in movement planning (Nachev et al., 2008; Andersen and Cui, 2009), movement preparation (Desmurget et al., 2009), and movement intention (Desmurget et al., 2009) by collaborating with M1, the functional connectivity observed in this study may represent the cortical networks of the mu rhythm related to controls of M1 activity during real and imagined movements.

On the other hand, during imagined movements, no significant functional connectivity was observed between the SMA and M1, whereas previous studies reported functional connectivity between these two regions during both imagined and real movements (Solodkin et al., 2004; Chen et al., 2009; Gao et al., 2011). These studies mainly used complex or sequential motor tasks to calculate functional connectivity, while we used three

**FIGURE 4 | Spatial distributions of significant correlations between decoding accuracy and connectivity during real and imagined movements.** Correlation coefficients between decoding accuracy and

imaginary coherence with the left M1 during real and imagined movements were calculated over the whole brain. Brain regions with significant correlations are represented in orange (p < 0.05, FDR-corrected).

#### **Table 3 | Coordinates of significant correlation coefficients between decoding accuracy and strength of IC during real movements.**


simple right upper limb movements for both real and imagined movements. Because the SMA has been shown to be more sensitive to complex and sequential actions than to simple ones (Nachev et al., 2008), it may be reasonable to expect that we would not observe significant functional connectivity between the SMA and M1 during imagined movements. Furthermore, several studies suggested that the network involved in real movements has a positive influence from SMA on M1, and during imagined movements, the SMA exerts a suppressive influence on M1 (Solodkin et al., 2004; Kasess et al., 2008). These results indicate that the functional connectivity between the SMA and M1 has different characteristics for information processing of real and imagined movements.

#### **RELATIONSHIP BETWEEN DECODING ACCURACY AND FUNCTIONAL CONNECTIVITY**

As described above, the importance of M1 signals for decoding movement types or movement directions has been previously demonstrated using amplitude of waveforms or low frequency components (Waldert et al., 2008; Yanagisawa et al., 2009),



sensorimotor rhythm (Mellinger et al., 2007), and high gamma power (Yanagisawa et al., 2012a). An MEG study by Waldert et al. (2008) showed that the low frequency components of neuromagnetic signals are more important for obtaining high decoding accuracy than either alpha-beta or gamma power. We also showed that the fluctuations in the amplitude of low frequency MEG signals carry enough information about hand and arm movements to decode movement kinematics (Sugata et al., 2012b). In the present study, we obtained significantly high decoding accuracy for both real and imagined movements using smoothed M1 signals from 40 virtual channels. Given that the amplitudes of waveforms or

**coefficient results for real and imagined movements.** Most of the significant correlations (orange) were located in or around brain regions (white dotted circles) with significant IC (blue). For imagined movements, there was no co-localization between significant correlations and

there. The lower panels indicate magnified figures with brain regions shown using white dotted lines and with corresponding letters on the upper panels. Black dotted lines in the lower panel indicate the central sulcus.

the low frequency components of the signals have higher signalto-noise ratios than the high frequency components, the decoding feature used in this study may be suited for classifying unilateral upper limb movements with single trial MEG signals, even though MEG is inferior to invasive cortical recordings with respect to the sensitivity in weak signals in the high frequency band.

Significant correlations between strength of functional connectivity and decoding accuracy during real movements were observed in motor association areas, such as the left postcentral gyrus and PMC, and the right sensorimotor area. These regions largely overlapped with or were located close to the brain regions with significant IC. Several previous studies reported that the activities of such motor association areas modulate M1 activity by integrating sensory-motor information and transforming the sensory information into motor representation (Luppino and Rizzolatti, 2000; Solodkin et al., 2004; Hoshi and Tanji, 2007; Kantak et al., 2012; Xu et al., 2014). In addition, activity of the PMC during real movement was shown to resemble the activity of M1 neurons, suggesting that the PMC is directly relevant to motor execution (Lee and van Donkelaar, 2006). Thus, the representation of motor information in M1 during real movement may depend

on the intensity of sensory-motor integration and activity of the PMC. Furthermore, regarding the correlation in the right sensorimotor area, it is possible that interhemispheric communication between the bilateral M1s, which contribute to controlling the cortico-spinal output from M1 (Avanzino et al., 2007; Lee et al., 2007), was relevant to the representation of motor information in M1.

Significant correlations between functional connectivity and decoding accuracy were observed at the left prefrontal cortex (including DLPFC), IPL, SPL, right IFG, and sensorimotor area during imagined movements. In the left hemisphere, significant correlations around the DLPFC and parietal area overlapped with the brain regions with significant functional connectivity observed in the group analysis. Previous studies showed that the DLPFC and parietal areas were more activated during imagined movement than during real movement (Vry et al., 2012) and play an important role in working memory (Smith and Jonides, 1999; Baddeley, 2003). Motor imagery, which involves simulating movement through the manipulation of visual and kinesthetic information, is a cognitive process that requires working memory (Munzert et al., 2009). Thus, these regions are thought to

play a significant role in generating clear motor imagery using working memory. Actually, deactivation of DLPFC was shown to decrease information processing during imagined movement in Parkinson's disease (Jahanshahi et al., 1995; Samuel et al., 1997, 2001). Furthermore, previous studies reported that the parietal area is associated with accuracy of the imagined movement (Hanakawa et al., 2003) and that lesions to this area reduce motor imagery abilities (Jackson et al., 2001; Lotze and Halsband, 2006; Mulder, 2007). Considering that our results showed that these brain regions are functionally connected with M1 and exhibit significant correlations between decoding accuracy and strength of functional connectivity, it is possible that the DLPFC and parietal area modulate the M1 activity related to the representation of motor information by interacting with the M1 during imagined movement.

In the right hemisphere, significant correlation between decoding accuracy and strength of functional connectivity was observed in the IFG during imagined movement but there was no spatial overlap with significant functional connectivity. The right IFG, but not the left, is relevant to the suppression of movement (Aron et al., 2004; Coxon et al., 2009), and impairment in this region resulted in a loss of suppression of movements in the inhibitory control task (Aron et al., 2004). In addition, a recent study suggested that activation of the right IFG during imagined movement may be relevant to an active inhibition process in the prevention of actual movement (Fleming et al., 2010). Furthermore, the right IFG was reported to be more activated in good imagers than in poor imagers for imagined sequential finger movements (Guillot et al., 2008). More recently, Gaetz et al. (2013) successfully demonstrated that gamma-band activity from the right IFG is observed for tasks involving response interference. Because our results showed that participants with strong functional connectivity between left M1 and right IFG exhibited high decoding accuracy, such inhibition processes in the IFG may work to generate the clear imagery of the movement necessary to decode the imagined movements. Nevertheless, no significant functional connectivity was observed between M1 and IFC in the group analysis, suggesting that connectivity between the two regions is not necessarily required for imagined movement and that the strategy for imagined movement may vary between individuals.

The significant correlation between decoding accuracy and strength of functional connectivity during imagined movement was also distributed to the right sensorimotor area, although there was no spatial overlap with significant functional connectivity. Previous studies reported that the sensorimotor area was activated on not only the contralateral side, but also the ipsilateral side during imagined movement (Kim et al., 1993; Porro et al., 2000). In addition, interhemispheric communication between the bilateral M1s was recently shown during imagined movement as well as real movement, suggesting that bilateral interactions of M1 play a crucial role in the modulation of the motor system during imagined movement (Liang et al., 2014). On the basis of these findings, our results suggest that the strength of functional connectivity between bilateral sensorimotor areas observed in this study may contribute to the

modulation of interhemispheric communication between the two regions, and that the subjects with strong connectivity may create more vivid motor imagery related to motor information than the subjects with low connectivity. However, as there was no significant functional connectivity between the bilateral sensorimotor areas in the group analysis, that connectivity may be not an essential component of imagined movement, but rather may be an additive one to generate the M1 activity similar to real movement.

Several studies reported anatomical and functional connectivities between the SMA and M1 (Muakkassa and Strick, 1979; Solodkin et al., 2004; Matsumoto et al., 2007; Nachev et al., 2008). This region was shown to play an important role in movement preparation and movement intention (Lau et al., 2004; Nachev et al., 2008). However, our results showed that there is no significant correlation between decoding accuracy and functional connectivity over the SMA for either real or imagined movements. As described above, the SMA is more sensitive to complex and sequential actions, while the task used in this study was a simple, unilateral, upper limb movement, so that the strength of functional connectivity between M1 and SMA may not have contributed to the decoding performance. It is possible that we can observe the contribution of functional connectivity between M1 and SMA to decoding performance if a more complex task is used. Further studies are needed to clarify this speculation.

#### **IMPLICATIONS FOR CLINICAL NON-INVASIVE BMIs**

To date, many researchers have tried to apply BMIs to patients with severe motor dysfunction using invasive methods, such as ECoG and local field potentials. When we put these invasive BMIs into clinical use, it is indispensable to perform a pre-operative, non-invasive evaluation to determine whether an invasive BMI would be effective. In addition, considering that BMIs are likely to be practically applied to patients with severe motor dysfunction, it is important to improve the performance of decoding accuracy for imagined movements. Recently, several studies reported predictors for the performance of BMIs using sensorimotor rhythm (Blankertz et al., 2010; Hammer et al., 2012), near-infrared spectroscopy activity (Fazli et al., 2012), high theta and low alpha powers (Ahn et al., 2013b), and gamma band activity (Ahn et al., 2013a) during real and imagined movements. In the present study, we focused on the relationship between alpha band functional connectivity and the decoding accuracy of real and imagined movements. As a result, significant correlations between these two aspects were mainly obtained in motor association areas, such as the PMC, sensorimotor areas, and the parietal area. This result suggests that we may be able to predict and improve decoding accuracy by evaluating and enhancing the functional connectivity between M1 and these brain regions, perhaps using neurofeedback methods as previously reported (Shibata et al., 2011; Koush et al., 2013). In particular, for imagined movement, because the strength of functional connectivity observed in this study may be relevant to generation of a vivid imagined movement for decoding movement types, enhancing these networks are important for improving the performance of imagery-based BMIs. Although we used brain signals extracted from the M1 gyrus in the present study, M1 signals from the central sulcus provides a rich

source of information representing movement types (Yanagisawa et al.,2009) and have a better signal-to-noise ratio for MEG recordings. They may contribute more to performance of imagery-based BMIs. Also, since the activities of sub-cortical regions as well as the cerebellum are also associated with the generation of voluntary movement (Shibasaki and Hallett, 2006), further investigation of whether functional connectivity between M1 and sub-cortical regions or the cerebellum can be detected and if they are related to the performance of BMIs using this method is thought to be necessary. Additional investigation may lead to the establishment a method for pre-operative evaluation or for the application of the present findings to clinical tools such as neurorehabilitation.

#### **CONCLUSION**

In this study, we examined the relationship between decoding accuracy and alpha functional connectivity during real and imagined movements. The significant correlations between decoding accuracy and the strength of alpha functional connectivity were mainly distributed to motor association areas. Our results indicate that alpha functional connectivity between M1 and the motor association area is important for the improved neural decoding of real and imagined movements. Further investigation may lead to the establishment of a method for pre-operative evaluation or for the application of the present findings to clinical tools such as neurorehabilitation.

#### **ACKNOWLEDGMENTS**

This work was supported by a grant for "Development of BMI Technologies for Clinical Application"from the Strategic Research Program for Brain Sciences by the Ministry of Education, Culture, Sports, Science, and Technology of Japan, a Health Labour Sciences Research Grant (23100101) from the Ministry of Health, Labour, andWelfare of Japan, and KAKENHI (23390347, 26282165) grants funded by the Japan Society for the Promotion of Science (JSPS).

#### **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: 22 April 2014; accepted: 23 July 2014; published online: 08 August 2014. Citation: Sugata H, Hirata M, Yanagisawa T, Shayne M, Matsushita K, Goto T, Yorifuji S and Yoshimine T (2014) Alpha band functional connectivity correlates with the performance of brain–machine interfaces to decode real and imagined movements. Front. Hum. Neurosci. 8:620. doi: 10.3389/fnhum.2014.00620*

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

*Copyright © 2014 Sugata, Hirata, Yanagisawa, Shayne, Matsushita, Goto, Yorifuji and Yoshimine. 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.*

### Comparison of EEG propagation speeds under emotional stimuli on smartphone between the different anxiety states

#### *Tetsuya Asakawa1 \*, Ayumi Muramatsu2, Takuto Hayashi 3, Tatsuya Urata4, Masato Taya5 and Yuko Mizuno-Matsumoto2*

*<sup>1</sup> Department of Physiology and Biological Information, Dokkyo Medical University, Tochigi, Japan*

*<sup>2</sup> Graduate School of Applied Informatics, University of Hyogo, Hyogo, Japan*

*<sup>3</sup> Department of Medical Engineering, Aino University, Osaka, Japan*

*<sup>4</sup> Faculty of The Physical Education, Osaka University of Health and Sports Science, Osaka, Japan*

*<sup>5</sup> KDDI R & D Laboratories Inc., Saitama, Japan*

#### *Edited by:*

*Leonides Canuet, Center for Biomedical Technology, Spain*

#### *Reviewed by:*

*Jing Xiang, Cincinnati Children's Hospital Medical Center, USA Leonides Canuet, Center for Biomedical Technology, Spain Yasunori Aoki, Osaka University Graduate School of Medicine, Japan*

#### *\*Correspondence:*

*Tetsuya Asakawa, Department of Physiology and Biological Information, Dokkyo Medical University, 880 Kitakobayashi, Mibu-machi, Shimotsuga-gun, Tochigi 321-0293, Japan e-mail: ici34774@nifty.com*

The current study evaluated the effect of different anxiety states on information processing as measured by an electroencephalography (EEG) using emotional stimuli on a smartphone. Twenty-three healthy subjects were assessed for their anxiety states using The State Trait Anxiety Inventory (STAI) and divided into two groups: low anxiety (I, II) or high anxiety (III and IV, V). An EEG was performed while the participant was presented with emotionally laden audiovisual stimuli (resting, pleasant, and unpleasant sessions) and emotionally laden sentence stimuli (pleasant sentence, unpleasant sentence sessions) and EEG data was analyzed using propagation speed analysis. 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. The low anxiety group propagation speeds at the medial sagittal for unpleasant stimuli in the 0–30 and 60–150 s time frames were higher than those of high anxiety group. The propagation speeds at 150 s for all stimuli in the low anxiety group were significantly higher than the correspondent propagation speeds of the high anxiety group. These events suggest that neural information processes concerning emotional stimuli differ based on current anxiety state.

**Keywords: EEG, smartphone, emotional stimuli, anxiety states, propagation speed**

#### **INTRODUCTION**

The Ministry of Internal Affairs and Communications of Japan has reported that the market penetration of mobile communication terminals such as cellular phones has increased from 0.3% in 1989 to 98.0% in 2012 (Ministry of Health, Labor and Welfare, 2012) and it continues to grow. Recently, smartphones have become popular as mobile communication terminals, acting not only as telephones but as information and communication terminals. The use of smartphones has been reported to produce new relationships over age, distance, and status (Ministry of Health, Labor and Welfare, 2012). Use of the smartphone has spread rapidly in recent years. Out of the world population, the penetration rate of smartphones has reached 70%. For many, the smartphone has become a necessity in private work and in daily life (ITmedia, 2008; Nikkei BPnet, 2008). On the other hand, smartphones are thought to be a cause of physical and mental stress, and overdependence or overuse of smartphones has been reported to cause depression, sleep disorders, and other symptoms of stress (Thomée et al., 2011). Stress induced by smartphone use has not been clarified and may be categorized as a new type. Biological responses to emotional stress may be different from responses to emotional stress and may cause an as yet undetermined condition (Krause et al., 1998).

The current study is based on emotional stress stimulation using audiovisual stimuli which have been used in previous studies (Hayashi et al., 2009; Mizuno-Matsumoto et al., 2011). Researchers examined the effect of emotional stimuli on a specific emotional stress reaction within a living organism using emotional stimuli on a smartphone. Although it is possible to display images, video, voice, and words using audio-visual stimuli to stimulate an emotional stress response, the primary feature of smartphones is logging in to e-mail and Web pages that communicate with friends and family, such as Facebook. Therefore, using the smartphone differs from conventional emotional stimuli because others may provide stimuli in displayed text and images.

Part of the current study included performing propagation speed analysis between neighboring electrodes. Results were also evaluated using power spectral analysis in EEG (Kamo et al., 2013). There are two possible methods which can be used for analysis between electrodes: coherence analysis and phase analysis. Propagation speed analysis is the primary analysis used in the current study which evolved from phase analysis. It is considered that high propagation speeds correlate with fast information processing. In previous studies, information propagation analysis has not been used for the evaluation of EEGs (electroencephalography) using emotional stimuli and differentiating between differences in anxiety state. The effect of anxiety states on information processing as measured using propagation speed analysis has not been previously studied.

The purpose of the current study is to evaluate the relationship between different anxiety states and neural information processing using emotional stimuli on a smartphone.

#### **MATERIALS AND METHODS**

#### **PARTICIPANTS**

The subjects included 23 healthy adults (9 males and 14 females) aged 22–45 years old (mean ages 25.5 ± 1.2 years). Written informed consent was obtained from all subjects before the start of the experiments. They had no history of neurological or psychiatric illness. The Ethical Committee of the University of Hyogo approved this investigation and informed consent was obtained according to the Declaration of Helsinki.

#### **PSYCHOLOGICAL TEST**

Prior to the start of experiments, psychological tests were conducted to assess the anxiety conditions of the subjects. The State Trait Anxiety Inventory (STAI), a psychological inventory based on a 4-point Likert scale and consisting of 40 questions on a self-report basis was used to assess subjects' anxiety levels (Julian, 2011). The STAI measures two types of anxiety: state anxiety, or anxiety about an event, and trait anxiety, or anxiety level as a personal characteristic. Higher scores are positively correlated with higher levels of anxiety (Hidano et al., 2000; Ueno, 2003).

In this study, all participants were assessed for the presence of anxiety states using the trait anxiety measurement of the STAI. They were categorized into two categories: psychosomatically healthy, comprising subjects with scores of I and II, and psychosomatically ill, comprising subjects with scores of III, IV, and V. According to the results of the STAI, subjects were divided into two groups. Subjects with scores I and II were classified as the "low anxiety group" and those with scores III, IV, and V as the "high anxiety group." Then the emotional stimuli-related cerebral activities of the two groups were compared.

#### **EEG EXPERIMENT**

After a brief instruction about the physiological test and the attachment of electrodes, subjects were placed in a shielded room (band cut filter of 0.5–500 MHz, attenuation of 40 dB) in a sitting position.

As an objective physiological evaluation, an EEG was performed for each subject. Electrooculogram (EOG) was recorded for the removal of eye-movement artifacts. The Ag-AgCl EEG electrodes were placed on 16 sites (Fp1, Fp2, F3, F4, C3, C4, P3, P4, O1, O2, F7, F8, T3, T4, T5, and T6) according to the International 10–20 System, which is recognized as a standard in cranial landmarks. The linked earlobes (Aav) were adopted as a reference. During recording, impedances of all electrodes were kept below 10 k-. The EEG data was digitized at a sampling rate of 500 Hz. **Figure 2** shows a subject's position during the sessions.

#### **EXPERIMENT PROTOCOL**

Each experiment consisted of a total of five sessions presenting emotional audio-visual and sentence stimuli using a smartphone. We performed Desire (HTC Co. Ltd.) of smartphone (see **Figure 1A**). Before use, the SIM card of the smartphone used in this study was unplugged so as not to cause artifacts.

Two different sessions using audio-visual stimuli was used: pleasant and unpleasant. Similarly, sentence stimuli were separated into pleasant and unpleasant sessions. The fifth condition was a rest session using a neutral audiovisual stimulus. The subjects were to follow the instructions that were displayed on the smartphone. As shown in **Figure 1B**, the screen showed random numbers during the audio-visual stimuli; the subjects selected the maximum or minimum numbers. Emotional stimuli sessions provided audio-visual stimuli of four types every 10 s for 40 s during the stimulation. **Figure 1C** shows an example of e-mail on a smartphone. At the same time, audio stimuli, sound effects, were presented.

Participants viewed relaxing pictures such as landscapes in the resting session, funny pictures such as animals in the pleasant session, and terror pictures such as horror movies in the unpleasant session. In the emotional sentence stimuli sessions, participants viewed funny sentences found in e-mail as pleasant sentence stimuli, and anxiety-provoking sentences from e-mail as the unpleasant sentence stimuli.

**Figure 2** shows that the procedure for a session consisted of three steps. The first step was "Control," in which the subject kept a quiescent state with eyes closed for 180 s; this phase was defined as a no-load state (Control). After the no-load state, the second step was "Task," in which the subject watched stimuli for 40 s. The third step of a session was "Recalling," in which the subject recalled the contents of the emotionally stimuli with eyes closed for 180 s. The second and third steps were performed three times consecutively. Thereafter, this series was performed in two replicates. This protocol was applied to each stimulus as 1 session, and a total of 5 sessions were performed which included resting, pleasant, and unpleasant stimuli. EEG measurements were performed continuously for 5 sessions. EEG results were analyzed in the "Recalling" phases. This study was designed based on experimental designs of previous research (Hayashi et al., 2009; Mizuno-Matsumoto et al., 2010; Miyagawa et al., 2013). **Figure 3** shows a subject's position during the sessions.

#### **ANALYSIS**

#### *Propagation speed analysis between neighboring regions*

Artifacts were removed from the EEG data with an IIR band-pass filter of 2–50 Hz; large artifacts were manually removed from the analysis periods.

Propagation speed analysis between the brain regions was conducted to investigate the functional relationship between regions in the EEG data. Our in-house programs in MATLAB ver.7.7 were used for analysis.

The analysis epoch was 150 out of 180 s for recall. The number of points for analysis was 2048 and the frequency resolution was 0.488 Hz. The propagation speed for 4.1 s per 1 analysis period was calculated in 23 combinations between neighboring brain regions.

Furthermore, the 150 s per task were divided into 5 time windows (30 s per window), and propagation speeds for each time window were calculated between each pair of regions. We performed propagation speed in the α (8.0–14.0 Hz) band. Propagation speeds were calculated for each of the 5 kinds of tasks: resting, pleasant, unpleasant, pleasant sentence, and unpleasant sentence stimuli.

#### *The definition of propagation analysis*

The definition of propagation analysis is as below (Saiwaki et al., 1997; Mizuno-Matsumoto et al., 2000a). The Propagation analysis is based on the coherence and phase analyses. This expression used on cross-correlation *Cxy*(*uptau*), cross-spectral **Sxy**(**ω**), real part *Kxy*(**ω**) and imaginary part *Qxy*(**ω**). The cross-correlation *Cxy*(**τ**) of two time-series signals x(t) and y(t) was defined by the following equation. In previous research, propagation speed has been used for analysis of abnormal EEG epileptic activity. It was able to capture propagation of abnormal EEG activity to localize brain lesions (Mizuno-Matsumoto et al., 2000b). Overall, propagation speed analysis is basically used to measure the propagation of postsynaptic potentials between electrodes and their abnormalities.

$$C\_{\mathbf{xy}}(\mathbf{r}) = \lim\_{T \to \infty} \frac{1}{2T} \int\_{-T}^{T} \mathbf{x}\left(t\right)\mathbf{y}\left(t+\mathbf{r}\right)dt$$

Here, τ was showed time interval. We performed the FFT of crosscorrelation function, and required the cross-spectral *Sxy*(**ω**) was

defined by the following equation.

$$\mathcal{S}\_{\text{xy}}\left(\omega\right) = \int\_{-\infty}^{\infty} \mathcal{C}\_{\text{xy}}\left(\mathbf{r}\right) e^{-j\omega \mathbf{r}} d\mathbf{r}$$

ω was showed angular frequency. *Sxy* (**ω**)divided into real part *Kxy*(**ω**) and imaginary part *Qxy* (**ω**). The phase by angular frequency ω between two signals was defined by the following equation.

$$\theta\_{\mathbf{x}\mathbf{y}} = \tan^{-1}\left\{\frac{\mathbf{Q}\_{\mathbf{x}\mathbf{y}}\left(\omega\right)}{\mathbf{K}\_{\mathbf{x}\mathbf{y}}\left(\omega\right)}\right\} \left(-\pi < \theta\_{\mathbf{x}\mathbf{y}} < \pi\right),$$

Furthermore, time lag of frequency component ω was defined by the following equation.

$$\mathfrak{r}\left(\omega\right) = \frac{\theta\_{\mathfrak{xy}}\left(\omega\right)}{\omega}$$

The propagation speed *vp*,*<sup>q</sup>* by distance *rp*,*<sup>q</sup>* between two electrodes p and q was defined by the following equation.

$$\nu\_{p,q} \left[ m/s \right] = \frac{r\_{p,q} \left[ mm \right]}{\mathfrak{r} \left( \omega \right) \left[ ms \right]}$$

#### *Statistical analysis*

We divided the analysis of propagation speeds between neighboring regions into four regions: medial sagittal, lateral sagittal, lateral coronal, and medial coronal (see **Figure 4**). Statistical analysis was used to measure the propagation speed between neighboring regions of each subject. The values between each region for each stimulus, time course, and group were compared using a Two-Way repeated-measures analysis of variance (ANOVA) with the Bonferroni correction for multiple comparisons. IBM SPSS Statistics ver. 21 was used for statistical analysis.

#### **RESULTS**

#### **PSYCHOLOGICAL TESTS**

Using the STAI, 7 subjects had a score of I and 6 subjects had a score of II. Four subjects were categorized as score III, 5 subjects scored as IV, and 1 subject had score V. Therefore, we divided the subjects into two groups split at score II. Thus, 13 subjects were categorized as Low anxiety group, and 10 subjects were categorized as High anxiety group, respectively (see **Table 1**).

#### **Table 1 | The number of the subjects and STAI.**


*Thirteen subjects were categorized as Low anxiety group, and 10 subjects were categorized as High anxiety group, respectively.*

#### **PROPAGATION SPEED**

**Figures 5**, **6** show a graphical representation of propagation speeds between each region in each stimulus of time series in Low anxiety group and High anxiety group, respectively. The thickness and color of the lines displayed represent the propagation speed. Bold red lines show a high level (greater than or equal to 0.5 [m/s] and less than 25 [m/s]); green moderate-thickness lines show an intermediate level (greater than or equal to 0.05 [m/s] and less than 0.5 [m/s]), and blue thin lines represent a low level (more than 0 [m/s] and less than 0.05 [m/s]).

The results show that propagation speeds in the interhemisphere (Fp1-Fp2, F3-F4, C3-C4, P3-P4, O1-O2) of the low anxiety group were fast propagation speeds from 0 to 150 s in resting stimuli; from 0 to 90 and 120–150 s in pleasant stimuli; or 30–150 s in unpleasant stimuli' 0–150 s in pleasant sentence stimuli; and 0–150 s in unpleasant sentence stimuli. The same goes for both hemispheres of directions. Notably, all stimuli were propagated in a sagittal direction in the low anxiety group. The propagation speeds in the inter-hemisphere of the high anxiety group from 0 to 30 s in unpleasant and unpleasant sentence stimuli were faster than 30–150 s in those.

#### **PROPAGATION SPEEDS STATICALLY ANALYSIS**

**Figure 7** shows the results of a Two-Way factorial ANOVA for the propagation speed in all time intervals. The vertical axis shows propagation speed and the horizontal axis shows each region examined in each task.

In comparing the two groups, the propagation speeds in the medial coronal region for the low anxiety group were significantly higher than that in the high anxiety group (*p* < 0.05) during resting and unpleasant sentence stimuli of 0–30 s. The high anxiety group showed significantly higher propagation speeds in the medial coronal region for resting and unpleasant stimuli of 120– 150 s than that in Low anxiety group (*p* < 0.05). Apart from the 0–30 s interval for unpleasant stimuli and the 60–90 s interval for unpleasant stimuli, the propagation speeds for resting, pleasant, and unpleasant stimuli in all time intervals for the medial sagittal region were significantly higher in the low anxiety group than that in High anxiety group (*p* < 0.05).

In comparison between regions, apart from the 0–30 s interval for resting stimuli and the 90–120 s interval for unpleasant stimuli in the high anxiety group between lateral coronal and medial coronal (*p* < 0.05), the propagation speed of the medial sagittal and medial coronal regions was significantly higher than the lateral sagittal and lateral coronal regions in both groups in all time intervals.

**time series.** Drawing propagation speeds between each region in each stimulus of time series in Low anxiety group. And the thickness and color of level; green moderate-thickness lines show an intermediate level, and blue thin lines represent a low level.

**time series.** Drawing propagation speeds between each region in each stimulus of time series in High anxiety group. And the thickness and color of lines displayed represent the coherence value. Bold red lines show a high level; green moderate-thickness lines show an intermediate level, and blue thin lines represent a low level.

**FIGURE 7 | Regional propagation speed in each time.** MS, medial sagittal; LS, lateral sagittal; LC, lateral coronal; MC, medial coronal. Two-Way factorial ANOVA for the propagation speed in all time. The vertical axis shows propagation speed and the horizontal axis shows each region examined in each task.

#### **DISCUSSION**

The present study was designed to evaluate information processing through EEG propagation speed analysis in the α (8.0–14.0 Hz) band for emotional stress stimuli presented on a smartphone. The difference in neural information processing based on the anxiety state of the subject was then extracted.

#### **PROPAGATION ANALYSIS IN THE PRESENT STUDY**

This study calculated propagation speeds between neighboring electrodes by performing an EEG propagation speed analysis using emotional stimuli presented using a smartphone. Previous studies have shown that analysis between electrodes was coherence analysis, phase-lock analysis, and phase synchronization analysis, and methodology and clinical cases is almost (Mizuno-Matsumoto et al., 2002, 2005; Cohen et al., 2009; Thatcher, 2012). Additionally, processing information related to word and object encodings has been reported (Stein et al., 1999; Weiss and Rappelsberger, 2000). Propagation speed analysis related to emotion stimuli has not been used in previous studies. The propagation speed analysis was developed based on delay analysis. Govindan et al. (2006) reported that by using propagation speed analysis, researchers could understand the nature of information flow, and ascribe the direction of information flow between centers based on delay analysis. Thus, it can be used to establish the connection and direction between the two signals using delay analysis (Govindan et al., 2005).

#### **EMOTION-RELATED PROPAGATION SPEED CHANGES IN ALPHA ACTIVITY**

Propagation speeds in α frequency were evaluated in relation to the difference between anxiety states of subjects.

In generally, α band activity occurs in an arousal state with resting eye close (Niedermeyer and Silva, 2005). And Knyazev (2007) was reported that α frequency was related to the visual identification function for emotion. However, it was noted that it increased by mental activity, internal activity, and short-term memory such as mental calculation, image, and maintenance of working memory (Kostyunina and Kulikov, 1996; Palva and Palva, 2007). Further, Jensen et al. (2002) were reported that α band activity increased in association with the maintenance of simple memory. In addition, α band activity is reduced by intense mental activity, including anxiety, vigilance, and attention to the stimulus (Avram et al., 2010). Neurotic subjects manifest a low appearance ratio of α frequency, which was associated with anxiety (Markand, 1990). We previously reported that α band activity in the non-stress group was significantly increased during the presentation of unpleasant emotional stimuli, indicating that brain functional activity influences stress resistance (Hayashi et al., 2006).

The results show that propagation speeds for resting and unpleasant sentence stimuli in the inter-hemisphere of the low anxiety group were fast at 0–150 s for both stimuli. However, those of the high anxiety group remained unchanged and showed a low value of propagation speeds.

#### **EMOTION-RELATED PROPAGATION SPEED CHANGES OF REGIONAL PARTS**

We divided the analysis of propagation speeds under emotional stimuli between neighboring regions into four regions in this study. The results of a Two-Way factorial ANOVA in **Figure 7** shows that more significant differences were noted in medial sagittal and medial coronal. In addition, propagation speeds were high for inter-hemisphere (Fp1-Fp2, F3-F4, C3-C4, P3-P4, O1- O2) and both hemispheres of directions (Fp1-F3, Fp2 -F4, F3-C3, F4-C4, C3-P3, P3-P4, P3-O1, P4-O2) in **Figures 5**, **6**. And the direction of propagation speeds at alpha frequency in low anxiety group under unpleasant and unpleasant sentence were from left hemisphere to right hemisphere, that in high anxiety group, by contrast, were from right to left.

The brain inputs emotional visual stimuli from the external word, and performs simple visual processing in the occipital area (Martin et al., 2002). Then, Vision Society of Japan (2001) reported that information is propagated through the cerebral area from the occipital area to the frontal area. Researchers have considered that the emergence of alpha frequency in the occipital region was controlled by a thalamic pacemaker. The thalamic pacemaker is mediated by sagittal nerve fibers that communicate between the frontal and occipital regions in the cerebral hemispheres, and suggested that the α frequency in the occipital region was propagated to the frontal area. Therefore, Okuma (1999) suggested that α frequency emerged from the frontal area. Previous studies have shown a relationship with α frequency and emotional stimuli (Asakawa et al., 2012; Kamo et al., 2013). Prior research indicates that high anxiety subjects manifested reduced information processing at inter-hemispheric region using threatening images tasks (Compton et al., 2008). In addition, Weissman and Banich (1999) hypothesized the callosum as a type of selective filter that can adaptively control information flow between the hemispheres.

The current results of propagation speed may indicate that subjects with high anxiety were not as responsive compared to those with low anxiety. We have suggested that non-properly information processing in the brain was performed by the anxiety state. Additionally, it seems that subjects with high anxiety felt more disgust and anxiety in response to unpleasant stimuli.

These events suggest that information processing in the brain for emotional stimuli differ based on the anxiety state of the subject. In the future, smartphone manufacturers and providers should be educated about the need to develop software and hardware optimal for minimizing anxiety responses.

#### **ACKNOWLEDGMENTS**

This work is partially supported by a Grant-in-Aid from the Ministry of Education, Culture, Sports, Science and Technology, JP (23500527).

#### **REFERENCES**


Okuma, T. (1999). *Clinical Electroencephalography*. Tokyo: Igaku-shoin.


**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: 25 July 2014; accepted: 25 November 2014; published online: 10 December 2014.*

*Citation: Asakawa T, Muramatsu A, Hayashi T, Urata T, Taya M and Mizuno-Matsumoto Y (2014) Comparison of EEG propagation speeds under emotional stimuli on smartphone between the different anxiety states. Front. Hum. Neurosci. 8:1006. doi: 10.3389/fnhum.2014.01006*

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

*Copyright © 2014 Asakawa, Muramatsu, Hayashi, Urata, Taya and Mizuno-Matsumoto. 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.*

### Frontal midline theta rhythm and gamma power changes during focused attention on mental calculation: an MEG beamformer analysis

#### *Ryouhei Ishii <sup>1</sup> \*, Leonides Canuet 2, Tsutomu Ishihara1, Yasunori Aoki 1, Shunichiro Ikeda1,3, Masahiro Hata1, Themistoklis Katsimichas 1, Atsuko Gunji 4, Hidetoshi Takahashi 5, Takayuki Nakahachi 5, Masao Iwase1 and Masatoshi Takeda1*

*<sup>1</sup> Department of Psychiatry, Osaka University Graduate School of Medicine, Suita, Japan*

*<sup>2</sup> Department of Cognitive and Computational Neuroscience, Centre for Biomedical Technology, Complutense University of Madrid, UPM, Madrid, Spain*

*<sup>3</sup> Osaka Psychiatric Medical Center, Hirakata, Japan*

*<sup>4</sup> Faculty of Education and Human Sciences, Course of School Education, Yokohama National University, Yokohama, Japan*

*<sup>5</sup> Department of Child and Adolescent Mental Health, National Center for Neurology and Psychiatry, National Institute of Mental Health, Kodaira, Japan*

#### *Edited by:*

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

#### *Reviewed by:*

*Fernando Maestú, Complutense University of Madrid, Spain William C. Gaetz, The Children's Hospital of Philadelphia, USA*

#### *\*Correspondence:*

*Ryouhei Ishii, Department of Psychiatry, Osaka University Graduate School of Medicine, D3 2-2, Yamada-oka, Suita, Osaka 565-0871, Japan e-mail: ishii@psy.med.osaka-u.ac.jp* Frontal midline theta rhythm (Fmθ) appears widely distributed over medial prefrontal areas in EEG recordings, indicating focused attention. Although mental calculation is often used as an attention-demanding task, little has been reported on calculation-related activation in Fmθ experiments. In this study we used spatially filtered MEG and permutation analysis to precisely localize cortical generators of the magnetic counterpart of Fmθ, as well as other sources of oscillatory activity associated with mental calculation processing (i.e., arithmetic subtraction). Our results confirmed and extended earlier EEG/MEG studies indicating that Fmθ during mental calculation is generated in the dorsal anterior cingulate and adjacent medial prefrontal cortex. Mental subtraction was also associated with gamma event-related synchronization, as an index of activation, in right parietal regions subserving basic numerical processing and number-based spatial attention. Gamma eventrelated desynchronization appeared in the right lateral prefrontal cortex, likely representing a mechanism to interrupt neural activity that can interfere with the ongoing cognitive task.

**Keywords: frontal midline theta, focused attention, arithmetic calculation, gamma band, magnetoencephalography (MEG), synthetic aperture magnetometry (SAM), beamformer, spatial filtering**

#### **INTRODUCTION**

Frontal midline theta rhythm (Fmθ) is a distinct train of focal 5–7 Hz theta waves which appears over medial frontal areas in the EEG of normal subjects when performing a broad range of cognitive tasks demanding mental concentration (Ishihara and Yoshi, 1972; Iramina et al., 1996; Sasaki et al., 1996; Ishii et al., 1999). Thus, this brain activity is considered to reflect focused attentional processing. Reports of enhanced Fmθ in the pre-shot phase of rifle shooting (Doppelmayr et al., 2008), during car driving (Laukka et al., 1995) and in meditation states (Aftanas and Golocheikine, 2001) provide support to this notion.

Since the first report of Fmθ made by Ishihara and Yoshi (1972), Fmθ has been investigated in a number of neurophysiological and neuroimaging studies. Earlier studies of Fmθ using scalp EEG reported widespread distribution of this activity in midfrontal sites, but an accurate identification of its cortical generators within the medial frontal cortex was lacking (Mizuki et al., 1980; Laukka et al., 1995; Iramina et al., 1996). This is mainly due to the low spatial resolution of EEG. To overcome this problem, a few studies looking at the anatomical correlates of Fmθ used fMRI with simultaneous EEG recording. Despite the fact that fMRI scanner noise may affect mental concentration in some individuals when engaged in mental reasoning tasks (Pripfl et al., 2006), these fMRI studies clearly visualized Fmθ activity localized to the anteromedial frontal cortex (Gevins et al., 1997; Mizuhara et al., 2004; Sammer et al., 2007).

In an attempt to visualize the magnetic counterpart of the EEG-recorded Fmθ activity and clarify its cortical sources, we previously used MEG in four normal subjects and found that a large area over the bilateral medial prefrontal cortex generated Fmθ during continuous mental calculation (Ishii et al., 1999). Thus, our findings provided further support for a specific role of the prefrontal cortex in focused attentional processing. Findings from another MEG study using different types of attention-demanding tasks, including mental calculation, suggested that Fmθ reflects alternative activation of the prefrontal and anterior cingulate cortex (ACC) in the human brain (Asada et al., 1999). Overall, neuroimaging MEG research taking advantage of the excellent temporal resolution and higher spatial resolution of this technique compared to EEG, has shed some light on the generators of Fmθ, and supported the concept that Fmθ reflects activation of neural networks involved in allocation of attention related to various types of cognitive stimuli.

To generate Fmθ and understand the neural correlates of attentional processing, arithmetic calculation has often been used as the attention demanding task (Mizuki et al., 1980; Iramina et al., 1996; Sasaki et al., 1996; Asada et al., 1999; Ishii et al., 1999). However, apart from attention-induced neural activity, little has been reported on activation patterns related to calculation itself in Fmθ experiments. Indeed, there have been few attempts to identify whether source-power changes in frequency bands other than theta also emerge when focusing attention on mental calculation. Recent neurophysiological studies looking at EEG event-related responses in mental addition and subtraction using a calculation strategy approach have focused on theta and alpha oscillatory activity (De Smedt et al., 2009; Grabner and De Smedt, 2011). Evidence from neuropsychological and neuroimaging studies indicate that several cortical areas across hemispheres are implicated in arithmetic processing (Menon et al., 2000; Gruber et al., 2001; Dehaene et al., 2003, 2004; Kong et al., 2005; Fehr et al., 2007; Ischebeck et al., 2009). For instance, multiplication operations, which require retrieval of arithmetic facts stored in rote verbal memory (verbal number manipulation) mainly induce activation of the left angular gyrus (Gruber et al., 2001; Dehaene et al., 2003; Ischebeck et al., 2009). This area is also implicated in complex arithmetic operations (Menon et al., 2000; Dehaene et al., 2003; Grabner et al., 2007) along with other regions such as the left inferior temporal gyrus (Gruber et al., 2001; Kong et al., 2005) and the inferior and medial parietal cortex (Chochon et al., 1999; Kong et al., 2005). In contrast, addition and subtraction, which require genuine numerical calculation (quantity representations), have often been reported associated with activation of the parietal cortex (Chochon et al., 1999; Menon et al., 2000; Dehaene et al., 2003, 2004; Fehr et al., 2007). This indicates that the large variation in cortical sources of arithmetic-induced activation across studies may be due to the existence of sharp disassociations between arithmetic operations (i.e., addition, subtraction, multiplication, and division) and calculation complexity.

Although substantial progress has been made toward characterizing the anatomical correlates of arithmetic operations, the underlying neural activity (i.e., oscillations in different frequency bands) has been largely unexplored. Findings from EEG studies on arithmetic processing suggested that engagement in simple mental calculation may be associated particularly with oscillatory activity power changes in the gamma band (Micheloyannis et al., 2002). There is increasing evidence that gamma oscillations are also involved in a variety of cognitive processes including visuospatial focused attention (Kaiser and Lutzenberger, 2003), visual perception, learning and memory (Kaiser and Lutzenberger, 2005). However, little has been reported on the possible implication of gamma oscillations in calculation-related attention or in arithmetic operations. In fact, previous MEG investigations based on mental calculation paradigms often used single dipole analysis to localize specifically theta activity sources (Iramina et al., 1996; Sasaki et al., 1996; Asada et al., 1999). Because cognitive processing is functionally related to serial and parallel activation of multiple brain regions (Ishii et al., 2009) as well as to cortical oscillations in different frequency bands (Pfurtscheller and Lopes da Silva, 1999), applying MEG-dipole models which identify center of gravity rather than the volume of activation (Herdman et al., 2003), and focusing exclusively in theta oscillations, might not be sufficient to visualize the extended network of sources related to focused attention and mental calculation in the human brain. Hence, the application of methods which can detect cognitive task-induced oscillatory response and localize the underlying cortical sources may help elucidate the role of cortical oscillations in mental arithmetic processing.

Considerable insight into the dynamics of oscillatory activity across the cortex is provided by beamformer, owing to its action as a spatially selective filter to MEG signals. This allows estimation of the oscillatory activity coming from a given location in the brain (Hillebrand et al., 2005). Furthermore, by applying beamformer in both the active and control time windows (e.g., during and prior to stimulation), task-related power changes in brain electric or magnetic activity can be assessed, as well (Vrba and Robinson, 2001; Brookes et al., 2007). Synthetic aperture magnetometry (SAM), a spatially filtering technique based on non-linear constrained minimum-variance beamformer, permits unambiguous three-dimensional mapping of cortical power changes within specific frequency bands during task performance. The accuracy of this map, however, relies on the correctness of the beamformer assumptions for the given data set (Robinson and Vrba, 1998; Ishii et al., 1999, 2009; Hillebrand et al., 2005). Using this method along with permutation tests for statistical group analysis of MEG data, we could accurately identify neural sources and underlying oscillatory activity power changes that were functionally engaged in auditory attention and memory updating process (Ishii et al., 2009), cortical organization of sensorimotor areas (Ishii et al., 2002), singing and vocalization (Gunji et al., 2007), and perceptual information processing (Herdman et al., 2003; Doesburg et al., 2013). Thus, SAM beamformer and permutation analysis have proven to be useful methods to visualize sources of cognitive task-induced oscillatory activity in the brain.

The purpose of this study was to use spatially filtered MEG by SAM technique, and permutation analysis to precisely localize cortical generators of the magnetic counterpart of Fmθ as well as other sources of oscillatory activity associated with mental calculation, particularly with arithmetic subtraction, as it requires genuine numerical calculation or quantity representation.

#### **METHODS**

#### **SUBJECTS**

Eleven healthy volunteers, who were mainly researchers at Osaka university (six males, aged 27–36 years, mean age 32 years), participated in this study. The subjects had no specific education or background that could facilitate mental calculation. All subjects were right-handed, as assessed by the Edinburgh Handedness Inventory (Oldfield, 1971). Informed consent was obtained from all subjects prior to the experiments. The study was performed in accordance with the Declaration of Helsinki, and approved by the Ethics Committee of the Osaka University Hospital.

#### **EXPERIMENTAL DESIGN**

The experiment consisted of two conditions: (1) eye-closed resting state (control interval) and (2) mental arithmetic state (active interval). During the mental arithmetic state, the subjects were asked to serially subtract 7 from 1000, as fast as possible, with their eyes closed, thereby generating Fmθ. This paradigm is also assumed to elicit activation of cortical areas subserving genuine numerical calculation or quantity manipulation, as subtraction problems involve more quantity representation compared to division operations, multiplication tables and small exact addition facts, that can be stored in rote verbal memory (Menon et al., 2000; Dehaene et al., 2003). When enhanced, rhythmic theta oscillations lasting for at least 10 s. were visually identified in the MEG recordings, a beeping sound was given to indicate the end of the arithmetic state and the beginning of the resting state for 10 s. (**Figure 1**). Thus, resting and mental calculation states were alternately recorded on each subject for a total of 8 trials, each one of 20-s duration. The task was carried out purely mentally to avoid movement-related artifacts. Subject's performance of mental arithmetic was not controlled during the MEG recording. However, prior to the experiments, like in previous functional neuroimaging studies using covert mental arithmetic (Rueckert et al., 1996; Kawashima et al., 2004; Mizuhara et al., 2004), practice sessions (serial subtractions) were performed. In these sessions, the subjects provided answers orally, and the accuracy of the calculations was checked. This served as an index of each subjects' arithmetic proficiency. All study participants showed an excellent performance in the practice trials.

#### **MEG DATA ACQUISITION**

MEG recordings were performed on all subjects in a magnetically shielded room using a helmet-shaped whole-head array of 64-channel SQUID sensor (NeuroSQUID Model 100, CTF Systems Inc.). Each of the 64 primary sensors used a first-order axial gradiometer flux transformer. Ambient magnetic noise was reduced further by synthesizing third-order gradiometer response in firmware using the reference SQUID sensor array (Vrba and Robinson, 2001). This was especially effective at reducing lowfrequency noise. The data were recorded with the subject sitting on a comfortable chair with the head positioned in the helmetshaped Dewar. A head position indicator with three small coils, placed at the nasion and bilateral preauricular points, was fixed on the scalp. MEG signals were digitized at a sample rate of 250 Hz, and filtered using a 60 Hz notch filter and 100 Hz low pass filter. The resulting data were recorded on disk and analyzed offline.

#### **ANATOMICAL MRI**

To convert the sources of MEG oscillatory activities into subjects' brain images, magnetic resonance imaging (MRI) scans were obtained for all subjects using a 1.0-T MRI system (Magneton Impact, SIEMENS Inc., Germany) or a 1.5-T Siemens Magnetom

during focused attention on mental calculation (active state) followed by a 10-s non-arithmetic period (control state).

Vision plus system (Siemens, Erlangen, Germany). MRI data consisted of T1-weighted axial anatomical images with an inplane resolution of 256 × 192 and 124 sagittal slices (1.4 mm thickness). Anatomical landmarks (i.e., nasion and bilateral preauricular points) were used to create an MEG head-based three-dimensional coordinate system.

#### **TIME FREQUENCY ANALYSIS**

Brain Electrical Source Analysis (BESA 5.0, MEGIS Software GmbH, Grafelfing, Germany) software was used to visualize timefrequency representations for MEG sensors in each subject. The BESA beamformer applies complex demodulation to transform time-domain MEG data into time-frequency data (Hoechstetter et al., 2004). This provides information on the envelope amplitude and the phase of a specified frequency band as a function of time. The complex demodulation consisted of a multiplication of the time-domain signal by a complex periodic potential function with a frequency equal to the frequency analyzed, followed by a low-pass finite impulse response (FIR) filter of Gaussian shape. This is equivalent to a wavelet transformation with constant wavelet width across frequencies. In the resulting complex signal, its magnitude corresponds to half the envelope amplitude of the filtered frequency band and its phase to the compound phase at that frequency. We analyzed the frequency range of 4–60 Hz in 2-Hz steps with a time sampling rate of 25 ms-steps. To obtain power values, the time-series MEG data were squared and averaged across all single trials under the respective conditions. From this time-frequency transformation, event-related synchronization (ERS) and event-related desynchronization (ERD) measures are obtained. ERS was denoted as an increase in power of oscillatory activity in a given frequency band during the mental arithmetic state (active interval) relative to the mean power during the resting state (control interval). The opposite phenomenon, suppression of rhythmic brain activity during the mental task compared to a rebound after the task, was denoted as ERD.

#### **WITHIN-SUBJECT SOURCE LOCALIZATION ANALYSIS (SAM ANALYSIS)**

The spatial distributions of ERS or ERD in different frequency bands functionally related to focused attention on mental calculation were estimated from the unaveraged MEG measurements using Synthetic Aperture Magnetometry (SAM) analysis (Robinson and Vrba, 1998; Ishii et al., 1999, 2009). The data were first subjected to the following bandpass filters: 4–8 Hz (theta), 8–15 Hz (alpha), 15–30 Hz (beta), and 30–60 Hz (gamma). Then, SAM was used to generate a 16 × 12 × 12 cm volumetric image of root-mean squared (RMS) source activity from the filtered MEG signals, with a 2.5 mm voxel resolution. As an adaptive beamformer, SAM applies a spatial filter specific for each brain voxel, to suppress the interference of unwanted signals from other locations including the environmental noise, thus estimating source power with high spatial resolution (Robinson and Vrba, 1998). The spatial filter at a given location is a linear projection operator defined by a set of coefficients, with one coefficient for each sensor, which is determined by minimizing the source power under a constraint of unity gain at the location of interest. Finally, the Student's t statistic was computed, on a voxel-by-voxel basis, as the difference between the estimated source power for the active (8 epochs of 10-s Fmθ during mental calculation) and the control (8 epochs of 10-s non-Fmθ during resting following mental calculation) states, divided by their ensemble standard error that included both instrumental (SQUID sensor) noise and experimental variance. The resulting functional image represents a Student's t statistic parametric map, which was then fused with the corresponding MRI, relating brain anatomy to function. The optimum orientation for the spatial filter at location of interest is determined by rotation of the source orientation in the tangential plane; the orientation at which the pseudo-Z statistic is maximized is then used as the optimum orientation for the source strength estimate at each location. Details on SAM procedures are described elsewhere (Robinson and Vrba, 1998; Ishii et al., 2009).

#### **STATISTIC GROUP ANALYSIS (PERMUTATION TESTS)**

For statistical analyses of the group data, the distribution of each individual's SAM image was transformed into a common anatomical space, the SPM T1 template space (Barnes and Hillebrand, 2003). First, the SAM volumes of each subject were co-registered with his/her three-dimensional anatomical MRI based on fiducial positions measured during the MEG acquisition. Transformation parameters that map the subject's MRI to the template space were then determined using SPM99 software (Wellcome Department of Cognitive Neurology, London, UK). The spatial normalized subject's data were subsequently obtained by applying the above transformation to the SAM volumes. A non-parametric permutation technique was applied to the normalized SAM results (SAM-permutation statistics) to determine voxels with significant values by comparing the grand mean pseudo *t-*value of a voxel and the distribution of permuted pseudo *t-*values. This distribution was computed by randomly rearranging the active and control conditions and averaging the newly calculated pseudo *t-*values. The omnibus null hypothesis of "no activation" anywhere in the brain was rejected if at least one *t-*value was above the critical threshold for α *<* 0.05 determined by 1024 permutations, thus correcting for multiple testing. Voxels with pseudo *t-*values above this critical 0.05 threshold were deemed regions of activation, and the corresponding voxels were then overlaid on a normalized structural MRI. For details of this procedure see a study by Chau et al. (2004) and our previous report (Ishii et al., 2013).

#### **RESULTS**

Focusing attention on mental calculation, specifically on serial arithmetic subtraction, resulted in significant source-power changes in theta and gamma frequency bands over different cortical regions. There were no significant power changes in any of the other frequency bands. **Table 1** summarizes the cortical distribution of task-related activation, as indicated by SAM-permutation analysis.

#### **THETA POWER CHANGES: FRONTAL MIDLINE THETA (Fm**θ**) ACTIVITY**

The visual inspection of the MEG recordings revealed increased rhythmic theta activity at around 5–7 Hz when the subjects were engaged in mental calculation compared to the resting condition. This rhythmic theta activity appeared over the frontal regions bilaterally (**Figure 2**). **Figure 3** shows the results of SAM analysis **Table 1 | Cortical regions showing significant task-related activation or deactivation in different frequency bands.**


*ERS, Event-related synchronization; ERD, Event-related desynchronization; BAs, Brodmann areas.*

of the subjects (*n* = 8) with prominent theta waves in medial frontal areas during the arithmetic task compared to the resting state. In addition to midfrontal theta oscillations seen in all subjects, theta ERS was also seen in other cortical areas in some subjects, with inter-subject variability. These areas included the anterior temporal, orbitofrontal, and dorsolateral prefrontal cortex (**Figure 3**). The permutation test results indicated a significant increase in theta activity, or theta ERS, in the medial prefrontal cortex (BAs 8 and 9) and the adjacent dorsal part of the ACC (BA 32) during periods of focused attention on mental calculation (**Figure 4**). **Table 1** summarizes the cortical distribution and values of task-induced activation or deactivation, as indicated by SAM-permutation analysis.

#### **POWER CHANGES IN GAMMA FREQUENCY BAND**

Focusing attention on mental calculation, specifically on serial arithmetic subtraction, resulted in significant source-power changes in gamma frequency bands over different cortical regions. There were no significant power changes in any of the other frequency bands. Time-frequency analysis for MEG channels in individual subjects showed a similar pattern of enhancement and suppression of oscillatory activity in the gamma band over parietal and frontal areas, respectively, predominantly in the right hemisphere. The statistical group analysis provided by SAM-permutation tests revealed that gamma activity (30–60 Hz) exhibited both significant ERS and ERD during the mental calculation periods. Gamma ERS was observed in the right intraparietal sulcus (IPS) and the adjacent posterosuperior and inferior parietal lobules, whereas the ERD was observed over the inferior frontal gyrus (BA 44) in the same hemisphere (**Figure 5**). Calculation-related gamma ERS and ERD values are provided in **Table 1**.

#### **DISCUSSION**

Using MEG and SAM-permutation analysis during continuous mental calculation, we clearly identified pronounced theta ERS, representing Fmθ, distributed over bilateral medial prefrontal regions and the dorsal area of the ACC (**Figure 4**). A striking finding was the identification of significant gamma power

**FIGURE 2 | MEG waveforms during the active and control conditions.** The location of channels showing frontal theta enhancement are indicated in red color on the MEG sensor map.

changes, in particular gamma ERS in the right IPS and adjacent cortex, and gamma ERD in the inferior frontal cortex that appeared concomitantly with Fmθ (**Figure 5**). This clearly shows that focusing attention on mental calculation results not only in Fmθ generation but also in the activation of neural networks involving the parietal and lateral prefrontal cortex, likely associated with the arithmetic processing of the task, with power changes in the gamma band representing the underlying neural activity.

#### **FRONTAL MIDLINE THETA SOURCES**

Despite the body of information obtained from earlier EEG and MEG studies, to date, a precise delineation of brain structures involved in Fmθ generation has been difficult. This is due in part to a lower spatial resolution of EEG compared to MEG, the use of dipole modeling which can only detect center of gravity of activated regions rather than cortical volume of activation (Herdman et al., 2003), and the limitation of the lack of group statistical analysis in previous MEG-SAM studies (Ishii et al., 1999). In the present study, SAM-permutation analysis showed theta ERS specifically over the dorsal part of the ACC and adjacent

**FIGURE 4 | SAM-permutation images of source power changes (event-related synchronization) in the theta (4–8 Hz) band.** Responses were calculated for the mental calculation (active state) vs. non-arithmetic condition (control state). The color bar represents pseudo-*t*-values. L, Left; R, Right; A, Anterior; P, Posterior.

medial prefrontal cortex bilaterally, when subjects were engaged in continuous mental calculation (**Figure 4**).

The ACC encompasses numerous specialized subdivisions that subserve a vast array of cognitive, emotional, executive, nociceptive and visuospatial functions (Bush et al., 2000; Womelsdorf et al., 2010; Ovaysikia et al., 2011). Interestingly, we noted that within the ACC, significant increase in theta power was observed particularly over the dorsal area (BA 32), which corresponds to the cognitive subdivision of the ACC. Among other functions such as control of motivation, error detection and working memory processing, the dorsal part of the ACC is implicated in modulation of attention, which is why this area is regarded as part of a distributed attentional network (Bush et al., 2000). Conflict resolution is another function traditionally associated with the ACC. In the context of arithmetic subtraction, this ACC function might be necessary to inhibit potential error responses during the task. Thus, our results are consistent with and extend the findings of previous investigations proposing the dorsal part or cognitive subdivision of the ACC and adjacent medial prefrontal cortex as the generators of Fmθ associated with focused attention and other cognitive functions during mental calculation (Ishihara and Yoshi, 1972; Sasaki et al., 1996; Asada et al., 1999; Ishii et al., 1999; Enriquez-Geppert et al., 2013).

It is noteworthy that enhanced theta activity is thought to reflect working memory processes, as well. For instance, frontal theta ERS has been associated with working memory load and attention demands in several neurophysiological studies (Gevins et al., 1997; Kahana et al., 2001; Jensen and Tesche, 2002; Onton et al., 2005). The simple subtraction task used in this study requires working memory processing, in particular the maintenance and manipulation of information, which is closely Anterior; P, Posterior.

associated with medial prefrontal cortex (i.e., Brodmann areas 8 and 9) function. Taking this into account, we believe that, in addition to focused attention-related activation, the theta ERS observed in our study might also reflect working memory process associated with mental calculation.

In our study, three subjects out of eleven didn't show any prominent theta activity during the calculation task, even though they carried out the task like other subjects. There are several previous studies suggesting some influential parameters which could affect the appearance rate of Fmθ among individuals, such as anxiety level and personality traits (Inanaga, 1998). Although we didn't check the anxiety level and personality traits in this study, some possible behavioral, structural and genetic factors which might be associated with prominent theta activity can be investigated as our future application.

#### **PARIETAL CORTEX INVOLVEMENT IN NUMBER PROCESSING**

Little has been reported on gamma band activity and arithmetic processing in previous Fmθ experiments using mental calculation as attention-demanding task. Indeed, several EEG/MEG studies on Fmθ induced by arithmetic tasks focused exclusively on source localization of theta oscillations and did not analyze fast frequencies, including the gamma band (Iramina et al., 1996; Sasaki et al., 1996; Mizuhara et al., 2004; Missonnier et al., 2006; Doppelmayr et al., 2008). Consistent with our findings of pronounced gamma ERS in the right IPS (**Figure 5**), this cortical area has been a major site of activation in neuroimaging studies on number processing. Based on its systematic activation whenever numbers are manipulated, regardless of number notation, the IPS has been regarded as a potential substrate for quantity or numbers representation common to all arithmetic tasks (Rueckert et al., 1996; Simon et al., 2002; Dehaene et al., 2003, 2004; Kong et al., 2005). Of note, subtraction operation, which was used for mental arithmetic in this study, usually elicits greater IPS activation compared to other arithmetic operations, in particular multiplication and division. This is partly explained by the fact that multiplication tables, and even small addition facts, can be stored in rote verbal memory, hence placing minimal requirements on quantity manipulation, whereas subtraction problems are not learned by rote and therefore require genuine quantity manipulation (Simon et al., 2002; Dehaene et al., 2003).

We also noted significant gamma ERS in cortical areas adjacent to the IPS, namely the inferior and postero-superior parietal (PSP) lobules, which have also been implicated in numerical operations, including subtraction of two or more digits and counting (Dehaene et al., 2003). However, unlike the IPS, the PSP region is not specific to the number domain. Rather, it plays a central role in a variety of visuospatial tasks including orienting of attention, regardless of whether working memory process is involved or not (Mitchell and Cusack, 2008; Olson and Berryhill, 2009). This suggests that PSP cortex activation mainly reflects the processing of attended items. In this context it is interesting that the so called internal "number line," a quasispatial representation on which numbers are organized by their proximity, can be likened to the core semantic representation of numerical quantity. Thus, it is conceivable that the same process of covert attention that operates to select locations in space can also be engaged when attending to specific quantities on the "number line" during mental calculation (Dehaene et al., 2003). This number-based spatial attention hypothesis supports the involvement of the PSP lobule not only in visuospatial processing, as previously reported, but also in non-visual mental arithmetic tasks, as suggested by our results.

#### **GAMMA EVENT-RELATED SYNCHRONIZATION (ERS) IN COGNITIVE PROCESSING**

It is noteworthy that activation in areas subserving basic arithmetic processing and number-based spatial attention in the parietal lobe manifested specifically as an increase in power in the gamma band (**Figure 5**). This supports previous observations of sustained EEG gamma oscillations during high-level mental activities, such as reading, learning, emotion and arithmetic subtraction (Fitzgibbon et al., 2004; Luo et al., 2013). Our finding is consistent with evidence demonstrating that activation of cortical regions induced by cognitive processes generally translates into synchronization of rhythmic neural activity at frequencies above 40 Hz, the so-called gamma synchronization (Lachaux et al., 2008). Induced gamma response has also been reported to reflect activation of task modality-dependent networks or stimulusrelated sensory/cognitive function in primary or association areas subserving the specific stimulus information processing (Jensen et al., 2007). Furthermore, synchronized gamma activity is considered to be involved in object representation, including internally driven representations (Bertrand and Tallon-Baudry, 2000) and in specific modalities of attention (Kaiser and Lutzenberger, 2003, 2005; Jensen et al., 2007). Taken together, our findings and those of earlier neuroimaging studies suggest that gamma band synchronization underlies number-related cognitive processing.

A separate line of research indicates that theta phase can modulate gamma power in certain brain regions (Scheffer-Teixeira et al., 2012). Kaplan et al. using MEG found that theta-gamma phase-amplitude coupling between medial prefrontal areas and medial temporal areas was linked to memory retrieval (Kaplan et al., 2014). Based on these findings, it would be interesting to explore in future studies whether theta and gamma phaseamplitude coupling between the medial prefrontal and parietal cortex, rather than activation of isolated cortical regions, might play a role in arithmetic processing.

#### **GAMMA EVENT-RELATED DESYNCHRONIZATION (ERD) IN COGNITIVE PROCESSING**

In contrast to a gamma ERS phenomenon in the parietal cortex, as an index of stimulus-related local activation, the significant power changes in the inferior frontal gyrus consisted of gamma ERD (**Figure 5**). This speaks in favor of a functional disassociation between prefrontal and parietal cortices during arithmetic processing, as proposed by an fMRI study by Menon et al. (2000). Although there is no clear explanation for the existence of simultaneous gamma ERS and ERD during focused attention on mental calculation, this finding is in line with recent evidence indicating that performing attentiondemanding cognitive tasks require not only activation of specific cortical regions but also deactivation of other regions that can interfere with the ongoing cognitive task, either in low-level sensory areas or high-level structures, such as the prefrontal cortex (Lachaux et al., 2008).

In this context, the gamma ERD found in the inferior frontal cortex may indicate that this region was not directly involved in the numerical processing. Rather, this region may play a more supportive role, for instance in managing parallel processes that might interfere with fast continuous subtractions, such as working memory-related processes or the speed of the calculation without seriously compromising arithmetic performance(Menon et al., 2000). Reports indicating that inhibitory control (inhibition of irrelevant information with working memory demands) is associated with a distributed network, involving the right dorsolateral prefrontal cortex (Garavan et al., 1999; MacDonald et al., 2000), ACC (Rubia et al., 2003), and the inferior parietal cortex (MacDonald et al., 2000; Garavan et al., 2002) provide further support to this argument.

In this study, we noted a right hemisphere laterality of the gamma source-power changes. Interestingly, previous fMRI and neuropsychological studies also found a predominant right hemisphere activation of frontal and parietal regions in association with mental subtraction (Fehr et al., 2007). Further, there is evidence that engagement in simple arithmetic, in particular subtraction operation, activates a neural network predominantly in the right hemisphere. This network is thought to serve as a common basis to which more regions in the left hemisphere are recruited for more difficult problems or different arithmetic operations (Kong et al., 2005).

Our findings should be interpreted with caution based on the limitation of the lack of behavioral performance data. This was due to the use of a mental arithmetic task to avoid movementrelated artifacts during the MEG recordings. Previous neuroimaging studies using fMRI or EEG have also used this paradigm (Rueckert et al., 1996; Kawashima et al., 2004; Mizuhara et al., 2004). However, prior to experiments, the subjects in this study performed practice trials, during which answers were given orally and the accuracy of the calculations was checked. This allowed us to measure each subjects' arithmetic proficiency and to ensure that they were capable to perform well on the task. Despite the use of a purely mental task, we found cortical activation, as indicated by specific oscillatory power changes, in a frontal-parietal network thought to be involved in focused attention (Ishihara and Yoshi, 1972; Ishii et al., 1999) and numeric/quantity representations (Dehaene et al., 2003, 2004), concomitantly with Fmθ. This strongly suggests that the subjects were actually performing the arithmetic task during the recordings.

#### **LATERALITY OF SOURCE-POWER CHANGES**

Consistent with results of earlier studies, we found that Fmθ appeared bilaterally over the medial prefrontal and ACC cortex. This confirms the involvement of both hemispheres in focused attentional processing. However, the gamma sourcepower changes in number processing areas and the prefrontal cortex, showed a unilateral distribution in the right hemisphere. This is unlikely to be related to the limitation of SAM beamformer of suppressing sources highly correlated in time across hemispheres to reduce noise in the signal (Brookes et al., 2007, 2008) because this effect generally applies to evoked-related responses which are time-locked to the onset of a specific stimulus in each single trial. Cognitive stimulus-induced oscillatory responses, however, jitter in the long range of time window after the stimulus onset from one trial to another (Pfurtscheller and Lopes da Silva, 1999; Tallon-Baudry and Bertrand, 1999). Thus, SAM is appropriate to reconstruct those task-related source-power changes in oscillatory activity or induced response, as it has been demonstrated in several MEG studies (Herdman et al., 2003; Gunji et al., 2007; Ishii et al., 2009).

A possible explanation for the right hemisphere laterality of the gamma source-power changes in this study may be based on the existence of sharp disassociations between arithmetic operations and calculation complexity (Menon et al., 2000; Gruber et al., 2001; Dehaene et al., 2003, 2004; Kong et al., 2005; Fehr et al., 2007; Ischebeck et al., 2009). As mentioned above, multiplication operations, which require arithmetic fact retrieval and rote memory, induce activation mainly of the left angular gyrus (Gruber et al., 2001; Dehaene et al., 2003, 2004; Ischebeck et al., 2009). This area which is associated with verbal number manipulation, appears to be also involved in complex arithmetic operations (Menon et al., 2000; Dehaene et al., 2003, 2004; Fehr et al., 2007; Grabner et al., 2007) along with the left inferior temporal (Gruber et al., 2001; Kong et al., 2005) and inferior parietal (Kong et al., 2005) cortex. Taking into account that we used simple subtraction as basic arithmetic operation, which is related to genuine quantity manipulations (Menon et al., 2000; Dehaene et al., 2003, 2004), it is not surprising that all those areas in the left hemisphere showing operation-specific activation were not found significantly activated in this study. Previous fMRI and neuropsychological reports suggesting a predominant right hemisphere activation of frontal and parietal regions in association with mental subtraction (Fehr et al., 2007) and with Arithmetical Reasoning Test performance (Langdon and Warrington, 1997) provide further support to our results.

Lesional studies have also provided evidence of right parietal cortex involvement in arithmetic processing. Dehaene and Cohen (1997) reported two acalculic patients who had structural lesions in the left subcortical areas or the right parietal cortex. They noted that the left lesional case had impaired rote arithmetic facts processing with preserved knowledge of numerical quantities. However, in line with our findings, the patient with right inferior parietal lesion showed a specific impairment of quantitative numerical knowledge, which was particularly remarkable for subtraction tasks (Dehaene and Cohen, 1997). Additionally, recent reports of intraoperative cortical electrostimulation in patients with brain tumors have confirmed an anatomofunctional organization for arithmetic processing within the right parietal cortex (Della Puppa et al., 2013).

Consistent with our SAM-permutation results, Micheloyannis et al. (2002)study using linear and non-linear EEG measures indicated that right hemisphere activation during simple arithmetic is manifested as increased power in the gamma band, while left or bilateral theta and alpha responses appear to be associated with the calculation strategy applied (De Smedt et al., 2009; Grabner and De Smedt, 2011). Further, there is evidence suggesting that engagement in simple arithmetic, in particular subtraction operation, activates a neural network predominantly in the right hemisphere, which serves as a common basis to which more regions in the left hemisphere are recruited for more difficult problems or different arithmetic operations (Kong et al., 2005). Based on this argument, we can speculate that the significant gamma synchronization seen in our study represent activation and connectivity of different brain areas, engaged in simple subtraction in the right hemisphere (i.e., right parietal cortex), and contralateral cortical areas could have been recruited if the complexity of the task was manipulated or other arithmetic operations were used. Although not directly tested in the current study, we suggest that MEG connectome might be able to provide novel insights on the biological mechanisms of higher cognitive functions and pathophysiology of neuropsychiatric diseases.

#### **CONCLUSION**

The findings of this study should be interpreted in the context of brain activation associated particularly with focused attention on mental arithmetic subtraction. Using MEG and SAMpermutation analysis, our results confirm and extend those of previous EEG and MEG studies indicating that Fmθ is generated in medial prefrontal cortex and dorsal ACC during cognitive tasks requiring focused attention and working memory process. Moreover, our results suggest that right parietal cortical areas which subserve basic numerical processing and numberbased spatial attention, namely the IPS and the adjacent posterosuperior parietal lobule, respectively (Dehaene et al., 2003), are activated during mental subtraction operations. The main contribution of our work, however, is the identification of gamma ERS as the underlying neural activity of this parietal sources. In addition, gamma ERD occurred in the right lateral frontal cortex during performance of the mental task, likely representing a mechanism to interrupt transiently local neural communication in cortical regions not relevant to the ongoing cognitive task. Overall, our findings demonstrate the feasibility of using MEG and SAM-permutation analysis to determine cortical network of sources related to focused attention or arithmetic calculation, and their underlying neural activity.

#### **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: 02 February 2014; accepted: 21 May 2014; published online: 11 June 2014. Citation: Ishii R, Canuet L, Ishihara T, Aoki Y, Ikeda S, Hata M, Katsimichas T, Gunji A, Takahashi H, Nakahachi T, Iwase M and Takeda M (2014) Frontal midline theta rhythm and gamma power changes during focused attention on mental calculation: an MEG beamformer analysis. Front. Hum. Neurosci. 8:406. doi: 10.3389/fnhum. 2014.00406*

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

*Copyright © 2014 Ishii, Canuet, Ishihara, Aoki, Ikeda, Hata, Katsimichas, Gunji, Takahashi, Nakahachi, Iwase and Takeda. 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.*

## Magnetoencephalography detection of high-frequency oscillations in the developing brain

#### **Kimberly Leiken1,2, Jing Xiang1,2\*, Fawen Zhang<sup>3</sup> , Jingping Shi <sup>4</sup> , Lu Tang1,2,4, Hongxing Liu1,2,4 and XiaoshanWang<sup>4</sup>**

<sup>1</sup> Department of Pediatrics, Magnetoencephalography (MEG) Center, Cincinnati Children's Hospital Medical Center, Cincinnati, OH, USA

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

<sup>3</sup> Department of Communication Sciences and Disorders, University of Cincinnati, Cincinnati, OH, USA

<sup>4</sup> Department of Neurology, Nanjing Brain Hospital, Nanjing Medical University, Jiangsu, China

#### **Edited by:**

Ryouhei Ishii, Osaka University Graduate School of Medicine, Japan

#### **Reviewed by:**

Ryouhei Ishii, Osaka University Graduate School of Medicine, Japan Leonides Canuet, Center for Biomedical Technology (CTB), Spain

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

Jing Xiang, Division of Neurology, MEG Center, Cincinnati Children's Hospital Medical Center, 3333 Burnet Avenue, Cincinnati, OH 45229, USA e-mail: jing.xiang@cchmc.org

Increasing evidence from invasive intracranial recordings suggests that the matured brain generates both physiological and pathological high-frequency signals. The present study was designed to detect high-frequency brain signals in the developing brain using newly developed magnetoencephalography (MEG) methods. Twenty healthy children were studied with a high-sampling rate MEG system. Functional high-frequency brain signals were evoked by electrical stimulation applied to the index fingers. To determine if the highfrequency neuromagnetic signals are true brain responses in high-frequency range, we analyzed the MEG data using the conventional averaging as well as newly developed timefrequency analysis along with beamforming. The data of healthy children showed that very high-frequency brain signals (>1000 Hz) in the somatosensory cortex in the developing brain could be detected and localized using MEG. The amplitude of very high-frequency brain signals was significantly weaker than that of the low-frequency brain signals. Very high-frequency brain signals showed a much earlier latency than those of a low-frequency. Magnetic source imaging (MSI) revealed that a portion of the high-frequency signals was from the somatosensory cortex, another portion of the high-frequency signals was probably from the thalamus. Our results provide evidence that the developing brain generates high-frequency signals that can be detected with the non-invasive technique of MEG. MEG detection of high-frequency brain signals may open a new window for the study of developing brain function.

**Keywords: magnetoencephalography, high-frequency oscillations, somatosensory cortex, wavelet, beamformer, pediatrics**

#### **INTRODUCTION**

Increasing evidence indicates that the brain generates both very low- and very high-frequency brain signals (Bowyer et al., 2012; Zijlmans et al., 2012; Haegelen et al., 2013). Low-frequency brain signals are typically referred to as "direct current" or "infraslow activity" (Bowyer et al., 2012). High-frequency brain signals are also called "high-gammas" (Huo et al., 2010), "high-frequency oscillations" (HFOs), "ripples," or "fast ripples" (Worrell et al., 2008; Engel et al., 2009; Gotman, 2010; Haegelen et al., 2013). High-frequency brain signals are potential new biomarkers for the study of brain function (Ozaki and Hashimoto, 2011).

In comparison with conventional low-frequency brain signals (<70 Hz), high-frequency brain signals are typically very weak (low in amplitude) (Worrell and Gotman, 2011). Neural signals from the brain are the spatiotemporal summation of synchronous firing from at least 10,000–50,000 neurons (Murakami and Okada, 2006). Therefore, it has been postulated that the functional organization of brain activity is encoded in multi-frequency ranges of that neuronal activity (Lina et al., 2014). The detection and recording of such high-frequency brain signals requires a highsampling rate. Therefore, how to handle the large datasets digitized

at a high-sampling rate has become an important issue for those investigating these signals (Worrell et al., 2012). In particular, using non-invasive techniques to detect and localize very high-frequency brain signals is rapidly becoming a new technique in many areas of research and clinical work (Miao et al., 2014; Rampp et al., 2014).

The most recognized HFOs are identified in the somatosensory cortex (Kotecha et al., 2009b). Previous reports have shown that the early portion of this activity is generated by action potentials of thalamocortical fibers and the late somatosensory HFO burst results from these action potentials arriving at the somatosensory cortices; specifically Brodmann areas 3b and 1 (Ozaki and Hashimoto, 2005). Though the exact range of HFOs remains unclear, HFOs up to 2632 Hz are detectable from the human somatosensory cortex with invasive recordings (Sakura et al., 2009). Since these oscillations are very weak, previous studies have required more than 1000 trials to obtain reliable results from adults (Ozaki et al., 1998). However, modern magnetoencephalography (MEG) systems have a whole-cortex sensor array, which can capture the spatial information of brain activity from variety of angles (Barkley and Baumgartner, 2003). This mechanism allows HFOs to be detected and localized with a much smaller number of

trials. Time-frequency analysis and spatial filtering (Kotecha et al., 2009b), described in further detail below, can be used for HFO investigations with MEG.

The objective of the present study was to assess high-frequency brain signals in children using newly developed MEG methods. The present paper will detail the mathematic algorithms of our new MEG analysis technique. Our central hypothesis was that somatosensory HFOs can be non-invasively localized using wavelet-based source localization. In comparison to previous reports on high-frequency brain signals (Jacobs et al., 2012), the major innovation of the present study was the optimization of MEG approaches for localizing very high-frequency neuromagnetic signals (>1000 Hz) in the developing brain. We consider this study of paramount importance because the development of reliable high-frequency signal detectors will likely have a significant impact on future clinical applications of MEG, such as presurgical brain mapping in those suffering from epilepsy (Gloss et al., 2014).

#### **MATERIALS AND METHODS**

#### **SUBJECTS**

Twenty healthy children (age range: 6–17 years, mean and SD: 12.3 ± 2.7 years, 10 girls and 10 boys) were recruited for this study. Inclusion criteria were (1) healthy without a history of neurological disorders or brain injuries; (2) age-appropriate levels of hearing, vision, and hand movement. Exclusion criteria were (1) inability to remain still inside the MEG scanner; and (2) presence of any non-removable metal device, such as a cochlear implant, a pacemaker, or a neurostimulator containing electrical circuitry, generating magnetic signals, or having other metal that could produce visible magnetic noise in the MEG data; (3) visually identifiable magnetic noise in the subject's recording (amplitude of waveforms >6 Pt). Since head movement during MEG recording could affect the accuracy of source estimation, data from a subject with head movement larger than 5 mm would be excluded from analysis. Written consents, formally approved by the Institutional Review Board (IRB) at Cincinnati Children's Hospital Medical Center (CCHMC) or Nanjing Brain Hospital, were obtained from each participant prior to testing.

#### **STIMULUS**

Two Digitimer Constant Current Stimulator model DS7A electrical stimulation systems (Digitimer Ltd., Welwyn Garden City, England) were used to stimulate the participant's right and left index fingers independently via Digital Ring Electrodes (Oxford Instruments Medical, Hawthorne, NY, USA). One hundred stimuli were delivered to each finger. The duration of the electrical stimulus was 0.3 ms (Kotecha et al., 2009b). The interval between two stimuli was 1010–1030 ms as the electrical stimuli were randomized by varying interstimulus intervals (10 ~ 30 ms) using the software BrainX (Kotecha et al., 2009b). The stimuli were adjusted to an intensity that delivered stimulation large enough to be detected by the participant yet not cause pain, as determined by previous studies (Kotecha et al., 2009b).

#### **MEG RECORDINGS**

Magnetoencephalography signals were recorded in a magnetically shielded room using a whole head CTF 275-Channel MEG system (VSM MedTech Systems Inc., Coquitlam, BC, Canada) in the MEG Center at CCHMC and at Nanjing Brain Hospital. Before data acquisition commenced, three electromagnetic coils were attached to the nasion, left and right pre-auricular points of each subject. These three coils were subsequently activated at different frequencies for measuring each subject's head position relative to the MEG sensors. Each subject lay comfortably in the supine position with his or her arms resting on either side, during the entire procedure. For the study of somatosensory evoked magnetic fields (SEFs), 100 trials were recorded for each finger (200 trials for each subject). The duration of each trial recording was 1000 ms (400 ms prestimulation baseline and 600 ms post-stimulation time-window). The sampling rate of the MEG recording was 6000 Hz per channel. MEG data were recorded with a noise cancelation of third-order gradients and without online filtering. To identify system and environmental noise, we routinely recorded one MEG dataset without a subject immediately prior to the experiment.

#### **MAGNETIC RESONANCE IMAGING SCAN**

Three-dimensional magnetization-prepared rapid acquisition gradient echo (MP\_RAGE) sequences were obtained for all subjects with a 3-T scanner (Siemens Medical Solutions, Malvern, PA, USA). Three fiducial points were placed in identical locations to the positions of the three coils used in the MEG recordings, with the aid of markers and digital photographs, to allow for an accurate co-registration of the two data sets. Subsequently, all anatomical landmarks digitized in the MEG study were made identifiable in the MR images. Pediatric brain templates developed by the MEG Center at CCHMC were also used for data comparison and visualization (Kotecha et al., 2009a).

#### **WAVEFORM AND TIME-FREQUENCY ANALYSES**

We visually inspected our MEG data and marked any possible artifacts over 6 Pt before data analyses. MEG waveforms evoked by 100 finger stimuli were averaged together over each finger for analyzing neural response amplitude and latency. Since somatosensory evoked activation has mainly been found in the range of 10– 120 ms (Kotecha et al., 2009b), our data analyses focused on this latency range. To analyze high-frequency brain signals in waveforms, we performed conventional waveform filtering (Kotecha et al., 2009b) with band-pass filters of 1–20, 20–500, 500–1000, and 1000–2000 Hz. Since somatosensory evoked MEG data were digitized at a sampling rate of 6000 Hz, a band-pass filter of 2000– 3000 Hz was also applied to the data. In this study, Butterworth filters were used (the phase shift was 0; the slope was −24 dB/oct).

The analysis technique of Morlet continuous wavelet transform was employed to transform time-domain data to time-frequencydomain data. The Morlet wavelet was used because brain activity is non-stationary and the wavelet is better suited for non-stationary data (Ghuman et al., 2011). The Morlet wavelet is described by the following equation:

$$G(t,f) = C\_{\sigma} \pi^{-\frac{1}{4}} e^{-\frac{1}{2}t^2} (e^{i\sigma t} - \kappa\_{\sigma}) \tag{1}$$

In the above formula, *t* indicates time and *f* indicates frequency (or a scale in the wavelet mother function for a specific frequency). Each wavelet transform has its own sigma value. κ<sup>σ</sup> represents the admissibility and *C*σ represents a normalized constant. If signals appeared in the given time-sensitive (a small sigma value) and frequency-sensitive (a large sigma value) ranges, they would be enhanced.

#### **FORWARD AND INVERSE SOLUTIONS**

To detect MEG signals at source levels, a three-dimensional source grid (3D grid) was developed. In the 3D grid, each grid node represented a possible source. Differing from the conventional volumetric source imaging or distributed source map, each grid node consisted of multiple data items including the strength and frequency of the source activity. Similar to previous reports (Mosher and Leahy, 1998; Vrba and Robinson, 2001; De Gooijer-Van De Groep et al., 2013), the sources of activity were determined with following equations:

$$B = LQ + N\tag{2}$$

In Eq. 2, *B* represents the MEG data, *L* represents the lead field, *Q* represents the source strength, and *N* represents the noise. For a given MEG dataset, *B* is known and *L* can be computed for each node with a forward solution. The forward solution in this study was computed according to Sarvas' formula for outside hemispherical conductors in Cartesian coordinates (Sarvas, 1987).

The determination of source strength and orientation of *Q* has been a challenge as discussed in many previous reports (Mosher et al., 1999; Huang et al., 2004; Robinson, 2004; De Munck and Bijma, 2009; Ou et al., 2009). According to our tests, for a given MEG dataset in multiple frequency ranges within a limited time window, the positions of the sensor array and the 3D source grid were fixed; consequently, lead fields could be computed once and then used for both low- and high-frequency ranges. Under these assumptions, we propose using single value decomposition (SVD) to decompose the lead field as in the following:

$$L = USV^T\tag{3}$$

Where *U*∈*R* mxm is an orthogonal (unitary in the complex case) matrix. The columns of *U* are the left singular vectors of *L*. *V*∈*R* mxm is an orthogonal matrix. The columns of *V* are right singular vectors of *L*. *S* = diag (σ1, σ2,. . . , σ*p*) is an *M* × *N* diagonal matrix with *p* = min (*m*, *n*) and σ1, σ2,. . . , σ*<sup>p</sup>* are the singular values of *L*. *M* indicates the number of sensors and *N* indicates the number of source orientations. For a single source, σ is <=3. The Moore–Penrose pseudo-inverse of *L* is given by:

$$L^{+} = V\mathcal{S}^{+}U^{T} \tag{4}$$

Where *S* <sup>+</sup> is a diagonal formed by the multiplicative inverses of the non-zero singular values of *L*. The correlation between measured MEG data, *B*, and the lead field is defined by:

$$B = LQ = USV^T Q \tag{5}$$

$$Q = BL^{-1} \tag{6}$$

By replacing *L* −1 in Eq. 6 with *L* <sup>+</sup> in Eq. 4, the estimated moment *Q*E can be computed with a SVD back substitution:

$$\vec{Q} = BVS^{+}U^{T} \tag{7}$$

Of note, *L* <sup>+</sup> could be computed once and used for the analysis of data in all frequency ranges, which makes the computation of source strength and probability more efficient. The probability of source activity was assessed with the correlation *t* value that was computed for the measured MEG signal and the computed MEG signals with Eq. 7. The threshold for correlations was statistically established to be 0.6. The equation for computing the statistical *t* is:

$$t = \frac{r}{\text{spt}[(1-r^2)/(N-2)]}\tag{8}$$

Where *r* represents the correlation coefficients and *N* indicates the number of samples. The parameters used for establishing the threshold are (a) sample size; (b) alpha value (0.05); (c) two-tailed test; and (d) type of correlation coefficient: Pearson's correlation. Of note, each voxel in our magnetic source imaging (MSI) had multiple values, which included source strength and source probability.

#### **MAGNETIC SOURCE IMAGING**

To go beyond localizing magnetic sources with aforementioned algorithms, we also developed a new technique, accumulated source imaging (ASI), to localize high-frequency signals. ASI was based on the previous observation that an accumulated spectrogram of virtual sensors could reliably detect high-frequency signals (Xiang et al., 2004, 2009a,b;Xiang and Xiao, 2009). ASI was defined as the volumetric summation of source activity over a period of time. In the present study of somatosensory evoked signals, ASI was the summation of source activity of all electrical pulse trials on each finger. ASI can be described as the following equation:

$$\text{Asi}(r, s) = \sum\_{t=1}^{n} Q(r, t) \tag{9}$$

In Eq. 9, ASI represents accumulated source strength at location *r*; *s* indicates the time slice; *t* indicates time point of MEG data; *n* indicates the number of trials; and *Q* indicates the source activity at source *r* and at a specific trial *t*. From a computer program point of view, the use of computer memory and storage space by Eq. 9 is dependent on the *s* for a fixed source imaging configuration (e.g., spatial resolution and dimension). Even though *n* could be infinitely increasing, the requirements for computer memory and storage remain the same. Consequently, the approach automatically avoided possible "overflow" or "out-of-space" problems in a large number of trials. The basic principle of ASI was that high-frequency MEG signals from the brain activity were locked to a spatial location and certain frequency bands. By accumulating all the source data computed for each location and each frequency band from multi-trial recordings, noise in a random space and frequency would be minimized and the signal-to-noise ratio in source imaging would then be increased.

Since this method spatially accumulates the results of source data, it is different from previously employed methods, which compute a covariance matrix or kurtosis of sensor data for a long recording. Specifically, using a covariance matrix or kurtosis for source localization is based on the assumption that the source was stationary during the long recording. Our approach, on the other hand, did not make this assumption, thereby taking into account minimal head movement typically found in pediatric MEG recordings. Therefore, our approach had the capability to detect both stationary and non-stationary source activity. Since high-frequency signals are very weak, we implemented the algorithms with C/C++ in double precision (64 bits). Therefore, the combination of spatial accumulation and double precision computation are well-suited to detect high-frequency signals over the course of at least 100 trials.

#### **STATISTICS**

Statistical comparisons between different frequency bands were performed using a Student's *t*-test. The normality of MEG data was tested with the Shapiro–Wilk test. For MEG data, which not found to be normally distributed, Mann–Whitney tests were used for comparisons. Statistical analyses were performed using SPSS version 16.0 for Windows (SPSS Inc., Chicago, IL, USA). The threshold of statistical significance for differences was set at *p* < 0.05. For multiple comparisons, a Bonferroni multiple comparisons correction was applied. Therefore, for the comparisons of 4 frequency bands, the *p* value would need to decrease to 0.012 (0.05/4).

#### **RESULTS**

Somatosensory evoked magnetic fields were analyzed with both band-pass filtering and time-frequency transforms as described above. The amplitude and latency of neuromagnetic responses are summarized in **Table 1**. We noted that the number of neuromagnetic responses in 500–1000 Hz was less than that in the other four frequency ranges. **Figure 1** shows examples of SEFs in the frequency ranges of 1–20, 20–500, and 500–1000 Hz in a healthy subject. **Figure 2** shows examples of SEFs in 1000– 2000 and 2000–3000 Hz. The amplitude and latency of SEFs in all healthy subjects are shown in **Table 1**. The result of statistical group analysis revealed that the amplitude of high-frequency signals was significantly weaker (lower) than that of the lowfrequency signals (*p* < 0.001). In addition, low-frequency signals appeared in a later latency while high-frequency signals appeared in an earlier latency (*p* < 0.001). **Figure 3** shows examples of global spectrograms and spectral contour maps of MEG signals in 1–20, 20–500, and 500–1000 Hz for a healthy subject. **Figure 4** shows examples of global spectrograms and spectral contour maps of MEG signals in 1000–2000 and 2000–3000 Hz for a healthy subject. There was no statistically significant difference in amplitude and latency between the left- and right-finger stimulations. We also did not find significant differences in amplitude and latency between the two age groups of subjects (6–11 vs. 12–17 years).

To better understand the spatial distribution, the oscillation covariance was computed for very high-frequency signals in all sensors. The contour map of oscillation covariance showed that very high-frequency signals among sensors over the somatosensory cortex were well-correlated (0.71 ± 0.12, *M* ± SD), which indicated that the signals at the sensors around the somatosensory region were from the same source. **Figure 4** shows the examples of three-dimensional contour maps of the oscillation covariance. **Table 1 | Latency and amplitude of somatosensory elicited magnetic fields in multi-frequency ranges (mean** ± **SD)**.


<sup>a</sup>Number of subjects that showed a response (deflections).

\*Compared with 1–20 Hz, p < 0.01.

\*\*Compared with 1–20 Hz, p < 0.001.

#Compared with 1000–2000 Hz, p < 0.01.

We also noted that very high-frequency signals in the occipital and frontal sensors were also partially correlated (0.61 ± 0.14, *M* ± SD), which implied deep brain activation (DA). The correlation between the occipital and frontal regions was relatively weak as compared with that in the somatosensory regions (**Figure 4**, "DA" and "SEF"). Since the spectral power of vHFOs in the occipital and frontal regions was very high, but the sensors in these regions showed a weak correlation (0.13 ± 0.02), those very highfrequency signals were considered to be noise or artifacts. In sum, the data indicated that very high-frequency signals were from the somatosensory cortex and deep brain regions, and may have included some occipital artifacts.

Magnetic source imaging revealed MEG signals at 1–20 and 20– 500 Hz from the somatosensory cortex in 17 subjects (17/20, 85%), MEG signals at 500–1000 Hz from the somatosensory cortex in 6 subjects (6/20, 30%), and MEG signals at 1000–2000 and 2000– 3000 Hz from the somatosensory in 18 subjects (18/20, 90%). MEG signals within the range of 1000–2000 and 2000–3000 Hz were also localized to the deep brain area, likely the thalamus, in 18 subjects (18/20, 90%) and 14 subjects (14/20, 70%), respectively. The result of statistical group analysis revealed that the source strength of MEG signals at 1–20 Hz was significantly stronger than that of MEG signals at 2000–3000 Hz in the primary somatosensory cortex (SI) for both left and right stimulation (*p* < 0.001). The result of statistical group analysis did not reveal significantly different SI source coordinates for each of the frequency bands (1–20, 20–500, 500–1000, 1000–2000, and 2000–3000 Hz). However, the source coordinates of the DA in 1000–2000 Hz was significantly different from that of the SI activation (*p* < 0.001). **Figure 5** shows the source locations of MEG signals in all thefrequency rangesfor a healthy subject. The sources of very high-frequency signals around the occipital regions were localized to the posterior regions, which were outside of the brain and were determined to be muscle artifacts. **Figure 6** shows an example of one of these muscle artifacts in the posterior regions.

**FIGURE 1 | Averaged waveforms and contour maps from a healthy subject showing neuromagnetic activation in 1–1000 Hz evoked by left finger stimulation**. The three sets of waveforms are filtered by three band-pass filters of 1–20, 20–500, and 500–1000 Hz, respectively. The amplitude of neuromagnetic activation decreases with the increase of frequency range. The latency of neuromagnetic activation shortens with the increase of frequency range. The three sets of waveforms have the same time ("5.3 ms") and amplitude ("40.0 ft") scales (the bottom). Waveforms in the 20–500 Hz range and the 500–1000 Hz range reflect a full time scale, as well as a "zoomed in" time scale from 10 to 30 ms to show the detailed deflections. The color contour maps show the distributions of neuromagnetic activation. All the contour maps have the same color scales.

#### **DISCUSSION**

The present study has demonstrated that the developing brain generates high-frequency neuromagnetic signals. Though the cerebral mechanisms underlying neuromagnetic high-frequency signals remain unclear, we noted that stimulation-induced highfrequency neuromagnetic signals were much weaker (lower amplitude) as compared with corresponding low-frequency neuromagnetic signals. Previous literature shows that magnetic signals detected by MEG are typically from ~10,000 to 50,000 synchronously active neurons (Murakami and Okada, 2006). If we assume that the 10,000–50,000 synchronously active neurons comprise

a "dipolar source" that can generate a signal strong enough to be recorded by MEG, then high-frequency neuromagnetic signals at the sensor level may be from multiple "dipolar sources." In other words, the frequency signature of neuromagnetic signals is not equal to the frequency signature of a single neuron firing; instead, the frequency signature of neuromagnetic signals may reflect the spatiotemporal and spectral patterns of multiple sources.

The non-invasive detection of high-frequency brain signals is still a new area of investigation (Gotman, 2010), the reliability and accuracy of the new methods are of utmost concern. We have been careful to eliminate artifacts in the present data. To eliminate possible artifacts, we routinely conducted noise tests just

before each clinical and research MEG recording. In addition, the patient's head position was monitored with three coils and our method volumetrically scanned the sources of the entire brain. In particular, our method was able to localize muscle artifacts to the occipital region. If the high-frequency signals obtained in the present study were artifacts, they should also be localized to a random place or out of the brain. Furthermore, we confirmed these results with time-frequency data, oscillation covariance maps, and MSI. Therefore, it is unlikely that the results reported here are due to the measurement of artifacts. Building on the data

from multiple approaches, we conclude that the measurements we have obtained are true neuromagnetic signals in high-frequency ranges.

The finding of high-frequency signals in the somatosensory cortex in the developing brain is consistent with previous reports of HFOs found in the somatosensory cortex in adults (Urasaki et al., 2002; Waterstraat et al., 2012). Compared with previous reports, there are several unique features in the present study. First, previous reports on somatosensory evoked fields were typically recorded by stimulating left and right median nerves

(Hashimoto et al., 1999). According to our clinical practice, the finding of median nerves for children can be challenging due to time constraints. Thus, the present study used finger stimulation with 100 trials; a location and trial number more suitable for children. Second, the present study used both the conventional band-pass filters and the newer time-frequency analyses. We demonstrated that the latency and amplitude of neuromagnetic responses to finger stimulation changed with different filters. Third, to our knowledge, this is the first study showing very highfrequency evoked signals (>1000 Hz) from the somatosensory cortex in the developing brain using MEG. Though the present study showed that high-frequency signals could be identified in the somatosensory cortex in children and adolescents, there were no statistical differences between two age groups of subjects (6–11 vs. 12–17 years). Thus, we assume that the developmental changes that occur between these two age groups are minor and/or the number of subjects was not large enough for the difference to reach significance. However, the present results suggest that future studies with a large number of subjects can employ high-frequency MEG signals to study the development of brain function.

The precise frequency ranges of high-frequency neuromagnetic signals from the somatosensory cortex vary among subjects. The cerebral mechanisms of individual variations among subjects in terms of frequency ranges are not yet well understood. Gobbele et al. (2000) have found 600-Hz bursts in the SEFs. We also identified activation in 20–500 and 500–1000 Hz. The percentage of high-frequency bursts in 500–1000 Hz appeared to be lower than that of the other frequency ranges. Though the exact cerebral mechanism remains unclear, we postulate that the stage of brain development and the methodologies used (e.g., stimulation paradigms, data analysis methods) play a role in these findings. Since our results have demonstrated that the developing brain generates high-frequency somatosensory evoked signals, it would be very interesting to standardize MEG protocols for future highfrequency brain signal investigations in children, adolescents, and adults in the future.

Though the main activation was seen in the contralateral somatosensory area (SII) following finger stimulation, we noted deep source activation coming from the ipsilateral thalamus in healthy subjects (see **Figure 5**, for example). Since the somatosensory tracts are already decussated, the activation in the ipsilateral thalamus might be related to the interhemispherical interactivation of the somatosensory system. One possibility is interhemispherical inhibition. Building on previous reports that finger stimulation is associated with deactivation of the ipsilateral SI (Hlushchuk and Hari, 2006), we postulate that the activation in the ipsilateral thalamus in our data might be related to the ipsilateral SI deactivation. In addition, we consider that activation in the ipsilateral thalamus might be also involved in mediating the activation of the secondary SII ipsilateral to the side of stimulation (Stancak et al., 2002).

As pointed out by Benar et al. (2010), high-pass filtering of waveforms for detection of oscillatory activity should be performed with great care. We have therefore used both filters and time-frequency analyses to verify the HFO findings in the developing brain. Our results from the two approaches strongly suggest that HFOs in the somatosensory system are non-invasively detectable. Of note, MEG HFOs can potentially be viewed as new biomarkers of brain activity, and MEG detection of highfrequency signals may open a new avenue for the study of the brain.

One potential flaw of the present study is the number of stimuli (100). This number was much less that used in previous studies [typically 3,000–5000 stimuli, e.g., Ozaki et al. (1998)]. This shorter paradigm was selected due to the fact that pediatric

populations tend to have shorter attention spans than the adults in the previous studies. However, the benefit of the new method (Morelet wavelet, source localization, and accumulated spectrogram) is that it is assumed to require much fewer trials in order to detect high-frequency signals. We recognize the lower trial number as a possible limitation of the present exploratory study. The software and supplementary materials, which implemented the aforementioned algorithms, are freely available from the following website (https://sites.google.com/site/braincloudx/ home) for other researchers to test, reproduce, and improve upon.

In summary, the results have demonstrated that somatosensory high-frequency activation can be non-invasively detected with MEG and advanced signal processing methodology. The proposed method was further validated with previously established conventional methods. MEG detection of high-frequency brain activity may open a new avenue in the study of the human brain function in the future.

#### **ACKNOWLEDGMENTS**

value that does not have a specific unit.

We thank Dr. Douglas Rose, Dr. Nat Hemasilpin, and Ms. Hisako Fujiwara for helping with MEG recordings. We would like to thank Mr. Kendall O'Brien for his technical help in performing MRI scans. We also thank Dr. Yingying Wang for helping with data analysis and management. We appreciate Ms. Chelsea Benson for her comments and edits. *Funding:* This project was partially supported by Grant Number R21NS072817 from the National Institute of Neurological Disorders and Stroke (NINDS), and the National Institute of Health. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

#### **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: 21 July 2014; accepted: 13 November 2014; published online: 12 December 2014.*

*Citation: Leiken K, Xiang J, Zhang F, Shi J, Tang L, Liu H and Wang X (2014) Magnetoencephalography detection of high-frequency oscillations in the developing brain. Front. Hum. Neurosci. 8:969. doi: 10.3389/fnhum.2014.00969*

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

*Copyright © 2014 Leiken, Xiang , Zhang , Shi, Tang , Liu and Wang . 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) 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.*

### Mapping memory binding onto the connectome's temporal dynamics: toward a combined biomarker for Alzheimer's disease

#### *Agustin Ibanez 1,2,3,4\* and Mario A. Parra5,6*

*<sup>1</sup> Laboratory of Experimental Psychology and Neuroscience (LPEN), Institute of Cognitive Neurology (INECO), Favaloro University, Buenos Aires, Argentina*

*<sup>2</sup> Laboratory of Cognitive and Social Neuroscience, UDP-INECO Foundation Core on Neuroscience, Diego Portales University, Santiago, Chile*

*<sup>3</sup> National Scientific and Technical Research Council, Buenos Aires, Argentina*

*<sup>4</sup> Universidad Autónoma del Caribe, Barranquilla, Colombia*

*<sup>5</sup> Psychology, Human Cognitive Neuroscience and Centre for Cognitive Ageing and Cognitive Epidemiology, University of Edinburgh, Edinburgh, UK*

*<sup>6</sup> Alzheimer Scotland Dementia Research Centre, Scottish Dementia Clinical Research Network, Edinburgh, UK*

*\*Correspondence: aibanez@ineco.org.ar*

#### *Edited by:*

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

### *Reviewed by:*

*Masafumi Yoshimura, Kansai Medical University, Japan*

**Keywords: Alzheimer disease, connectome, EEG/MEG, biomarkers, early detection, mild cognitive impairment**

The classical amyloid cascade hypothesis of Alzheimer's disease (AD) has driven research and clinical practice for several decades. It states that the deposition of the amyloid-β peptide in the brain parenchyma initiates a sequence of events that ultimately lead to atrophy and AD dementia. This proposal stimulated the study of specific brain regions mapped along the neurodegeneration sequence (e.g., hippocampus) and their associated impaired functions (e.g., episodic memory). Although anticipated by Mesulam more than 20 years ago (e.g., Mesulam and Asuncion Moran, 1987), it was not until recently that this view has started to change, largely due to the disappointing results of trials relying on the beta-amyloid cascade hypothesis and its variants, e.g., the synaptic beta-amyloid hypothesis.

A new approach based on neuroplasticity of neural networks has shifted the attention from one region to the orchestration of several brain hubs. The connectivity account would play an important role in revealing the transneural spread of misfolded proteins through neural networks in neurodegenerative disease (Pievani et al., 2011; Ibanez and Manes, 2012), and specifically in AD (Raj et al., 2012). In line with this view, the metabolism hypothesis (MH) has been proposed, which suggests that changes in the default mode network (DMN, the ongoing low-frequency fluctuations during resting state between the anterior and posterior cingulate cortex as well as the precuneus) stimulate an activitydependent or metabolism-dependent cascade that promotes the development of the AD pathology (Buckner et al., 2005). Notably, hyperactive neurons are observed near amyloid plaques in animal models (Busche et al., 2008) and in humans, connectivity hubs overlap the anatomy of A-β deposition (Buckner et al., 2009). Abnormal DMN activity discriminates between Mild Cognitive Impairment (MCI), AD and controls (Rombouts et al., 2005; Petrella et al., 2011; Seo et al., 2013; Wang et al., 2013b), and predicts AD conversion (Pievani et al., 2011). Thus, default connectivity seems to be a promising approach to reveal novel mechanisms leading to AD.

However, the DMN seems to be affected by many other diseases (Sonuga-Barke and Castellanos, 2007; Whitfield-Gabrieli and Ford, 2012). Moreover, although the altered DMN might interrupt or affect the brain dynamic in AD patients, it actually reflects a resting state activity unlikely to explain, on its own, profiles of cognitive decline in AD. The combined analysis of brain connectivity associated with specific cognitive processes affected by AD early in its course with resting DMN is an unexplored area that could help overcome these limitations. This approach may reveal markers for the early or even preclinical detection of neurocognitive impairments in AD. A potential strategy would be to assess the neural connectivity associated with episodic memory tasks (Wang et al., 2013a). However, these tasks have only detected AD-related changes in its prodromal or clinical stages (Fields et al., 2011). A recently developed methodology, namely short-term memory binding (STMB, Parra et al., 2009, 2010, 2011), is intrinsically related to brain networks activation, and appears to be more promising for the preclinical detection of AD. Binding functions, as originally investigated in perception, require a large-scale network integration mechanism (Varela et al., 2001). In AD research, emerging evidence suggests that binding impairments occur at the short term memory level. STMB is a cognitive function responsible for retaining, on a temporary basis, intra-item features thus contributing to the formation of objects' identity. It has been recently assessed with a change detection tasks, which ask participants to judge whether the content of two consecutive arrays of shapes, colors or shape-color combinations is the same or different. STMB is impaired early in AD (Parra et al., 2009) and also in preclinical familial AD (Parra et al., 2010, 2011), preserved in healthy aging (Brockmole et al., 2008), and declines earlier in AD than in other dementias (Della Sala et al., 2012). Moreover, STMB seems to be indexing a pre-hippocampal phase of AD (Reiman et al., 2010) and recruits other regions (Parra et al., 2014). Thus, a combination of both neural (DMN) and cognitive (STMB) integration processes may contribute an early and specific marker of AD progression. A new research agenda linking resting brain dynamics (DMN) with an active task which (1) relies on neural network integration and (2) is highly specific and sensitive for AD, would represent a powerful approach to the early detection of AD. Importantly, the resting and active connectivity measures drawn from this novel approach can be tracked both with neuroimaging (fMRI) and electromagnetic methods (EEG and MEG).

STMB tasks are well suited to be tracked with EEG or MEG measurement (Luria and Vogel, 2011; Wilson et al., 2012). Typically, the temporal dynamic of the integrative functions assessed by this task would be in the order of milliseconds. These techniques can capture the evoked responses and their neural connections during the whole process of STMB. This task also provides several advantages for EEG/MEG procedures such as high number of trials, stimuli-evoked activity pre and post memory binding process, temporal sequencing, and categories with different levels of difficulty.

Notwithstanding the adequacy of the task's parameters for EEG/MEG recordings, the question of why techniques with low spatial resolution should be considered instead of neuroimaging methods stands out. Several factors support this selection. First, high-density electroencephalography (hd-EEG) and other electromagnetic techniques permit an easy, low-cost, non-invasive, and accessible approach for large-scale multisite studies around the world. Second, high density EEG/MEG technologies have provided an increased spatial resolution of fine temporal dynamics both at the analytical (mathematical methods) and technical (high number of channels, photometry methods, and individual MRI co-recordings) level, making them more suitable to investigate AD-related changes. Third, EEG/MEG techniques have proven useful for characterizing AD and also for detecting changes in preclinical familial AD and MCI (Jackson and Snyder, 2008; Stam, 2010). For example, source EEG functional network disruption in AD is associated with cognitive decline (Gianotti et al., 2007; Ishii et al., 2010) (see (Kurimoto et al., 2008; Hsiao et al., 2013) for similar results in MCI), APOE genotype (Canuet et al., 2012) and differentiates between other dementias (Babiloni et al., 2004). Moreover, loss of interregional synchronization between different functional brain regions also reflects cognitive decline in AD (De Haan et al., 2012). Furthermore, the EEG/MEG based connectivity analysis (EMCA) can also be used to track the effect of medication on AD (Babiloni et al., 2006; Gianotti et al., 2008).

There are other direct advantages of using EMCA. The theoretical frame of interdependence between spontaneous and evoked neuroelectric oscillations in terms of frequency and phase reset has been forwarded earlier (e.g., by Basar' group in the eighties). Nevertheless, over the last 10 years, a real increase in technical and mathematical sophistication in EMCA has produced new research possibilities with practical applications (Basar et al., 2013; Larson-Prior et al., 2013). Regarding brain global properties, current graph theoretical network studies of the brain have shown a self-organized small-world network characterized by a combination of focal connectivity and long-distance connections. Graph theory is one of the most powerful forms of connectivity analysis for AD (Pievani et al., 2011; Tijms et al., 2013), and it can be correctly implemented with EEG/MEG signals (e.g., Stam and Van Straaten, 2012; Barttfeld et al., 2013). Regarding time dynamics, different networks are orchestrated in our brain in time windows of milliseconds and the connectivity within and between them is not a static process. Our brain has rapid rhythms that allow for communication between different regions at several frequencies. The high time-resolution of intracranial signals from EEG sources can be quantified by coherence and phase synchronization, two methods that have proved informative in AD (Czigler et al., 2006; Knyazeva et al., 2012). Moreover, recent methods provide a better characterization of the physiological signal with better spatial location and provide solutions for classical problems of volume conduction (Pascual-Marqui et al., 2011). They also permit the comparison of spontaneous and stimulus-induced activations and the identification of commonalities between them (Lehmann et al., 2010). Moreover, oscillatory neuronal dynamics in the human brain using connectivity analysis of source estimated event-related synchronization at different frequencies is now an available method (Ishii et al., 2013). Transient momentary events (e.g., thoughts) in electromagnetic signals might be incorporated in temporal chunks of processing (10–100 ms) as quasi-stable brain states (Lehmann et al., 2006). In brief, EMCA can track the brain dynamic of rapid fluctuations (Barttfeld et al., 2014) and also of transient activity during very short periods, especially those supporting binding or transient integration.

Thus, assessing the combination of basal resting state influences together with the ongoing activity during a task and its evoked neural response may allow for investigation and characterization of preclinical AD-related changes in brain dynamics. This novel approach offers new possibilities to better understand the cognitive binding problem in the course of AD as well as the dynamics of cortical integration. This research proposal also requires tackling several methodological and empirical challenges. Although promising, current methods for combining connectivity measures during ongoing activity and evoked responses do not yet fully capture the single trial dynamics. The potential use of connectivity metrics as predictor of patients' clinical outcome is not well understood at present. No single study using EMCA to assess AD or MCI patients has combined active and resting recordings and the analysis of source connectivity using individual MRI. We believe this combined approach would disclose unknown AD mechanisms. The research growth in these domains of cognitive neuroscience will offer support to key strategies such as combining STMB and EMCA to provide a neurocognitive marker for the preclinical detection of AD. In support to this proposal, recent neuroimaging studies carried out in cases of preclinical familial AD have revealed a temporal proximity between the onset or appearance of STMB deficits and amyloid-β deposition. By the average age at which amyloid-β depositions reach a plateau (Fleisher et al., 2012), STMB deficits become detectable behaviorally (Parra et al., 2010). It is worth noticing that these mutation carriers would be otherwise completely asymptomatic. This evidence warrants investigation of the hypothesis of a link between connectivity problems as assessed by STMB and EMCA and neurodegenerative changes in AD. Such research would shed further light into the link, or lack thereof, between amyloid changes, cognition and AD.

#### **ACKNOWLEDGMENTS**

This study was supported by grants CONICYT/FONDECYT Regular (1130920 and 1140114), FONCyT-PICT 2012-0412, FONCyT-PICT 2012- 1309, CONICET and the INECO Foundation. Mario A. Parra work is supported by Alzheimer's Society, Grant # AS-R42303.

#### **REFERENCES**


posterior alpha activity is correlated with cognitive impairment in early Alzheimer's disease. *Psychogeriatrics* 10, 138–143. doi: 10.1111/j.1479- 8301.2010.00326.x


in psychopathology. *Annu. Rev. Clin. Psychol.* 8, 49–76. doi: 10.1146/annurev-clinpsy-032511- 143049

Wilson, K. E., Adamo, M., Barense, M. D., and Ferber, S. (2012). To bind or not to bind. *J. Vis.* 12, 14. doi: 10.1167/12.8.14

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

*Received: 21 January 2014; paper pending published: 19 March 2014; accepted: 01 April 2014; published online: 22 April 2014.*

*Citation: Ibanez A and Parra MA (2014) Mapping memory binding onto the connectome's temporal dynamics: toward a combined biomarker for Alzheimer's disease. Front. Hum. Neurosci. 8:237. doi: 10.3389/ fnhum.2014.00237*

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

*Copyright © 2014 Ibanez and Parra. 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) 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.*

### Magnetoencephalographic study on facial movements

#### **Kensaku Miki 1,2\* and Ryusuke Kakigi 1,2**

<sup>1</sup> Department of Integrative Physiology, National Institute for Physiological Sciences, Okazaki, Japan

<sup>2</sup> Department of Physiological Sciences, School of Life Science, The Graduate University for Advanced Studies (SOKENDAI), Hayama, Kanagawa, Japan

#### **Edited by:**

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

#### **Reviewed by:**

Ryouhei Ishii, Osaka University Graduate School of Medicine, Japan Gavin Perry, Cardiff University, UK

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

Kensaku Miki, Department of Integrative Physiology, National Institute for Physiological Sciences, 38 Nishigonaka Myodaiji, Okazaki 444-8585, Japan

e-mail: kensaku@nips.ac.jp

In this review, we introduced our three studies that focused on facial movements. In the first study, we examined the temporal characteristics of neural responses elicited by viewing mouth movements, and assessed differences between the responses to mouth opening and closing movements and an averting eyes condition. Our results showed that the occipitotemporal area, the human MT/V5 homologue, was active in the perception of both mouth and eye motions. Viewing mouth and eye movements did not elicit significantly different activity in the occipitotemporal area, which indicated that perception of the movement of facial parts may be processed in the same manner, and this is different from motion in general. In the second study, we investigated whether early activity in the occipitotemporal region evoked by eye movements was influenced by the facial contour and/or features such as the mouth. Our results revealed specific information processing for eye movements in the occipitotemporal region, and this activity was significantly influenced by whether movements appeared with the facial contour and/or features, in other words, whether the eyes moved, even if the movement itself was the same. In the third study, we examined the effects of inverting the facial contour (hair and chin) and features (eyes, nose, and mouth) on processing for static and dynamic face perception. Our results showed the following: (1) In static face perception, activity in the right fusiform area was affected more by the inversion of features while that in the left fusiform area was affected more by a disruption in the spatial relationship between the contour and features; and (2) In dynamic face perception, activity in the right occipitotemporal area was affected by the inversion of the facial contour.

**Keywords: MEG, facial movements, MT/V5, fusiform area, occipitotemporal area**

#### **INTRODUCTION**

The "Face" provides much important information in our daily lives. Many studies on face perception that used microelectrodes on monkeys and humans detected face-specific neurons in the temporal cortex, mainly in the superior temporal sulcus (STS), and convexity of the inferior temporal (IT) cortex. Face perception processes have been reported in psychological studies (e.g., Bruce and Young, 1986), and many studies have examined the mechanisms underlying human face perception in detail using neuroimaging methods, including electroencephalography (EEG) recorded from the scalp (event-related potentials, ERPs) and cortical surface (electrocorticography, ECoG), magnetoencephalography (MEG), functional magnetic resonance imaging (fMRI), and near infrared spectoroscopy (NIRS).

Two important factors for face perception are static face perception and facial movement perception. The fusiform gyrus in the IT area may be specific to static face perception. A negative ERP component, being maximum at approximately 170 ms, is evoked by face stimuli in the bilateral temporal area. This has been referred to as N170. N170 is known to be larger for the face than for other objects (for example, a car and chair), and this reflects face perception (e.g., Rossion and Jacques, 2008). In ECoG studies, a negative component (N200) was detected in the small regions in the fusiform and IT gyrus evoked by faces but not by other stimuli (Allison et al., 1994a,b). The ECoG can been used to investigate the temporal and spatial aspects of the mechanisms responsible for face in detail, but it is an invasive method. MEG also has high temporal and spatial resolutions, and it is a noninvasive method. Therefore, MEG is useful to investigate face perception in normal volunteers. In MEG studies, a component was found to be evoked by the face, M170, which corresponded to N170 in EEG studies, and its activity was estimated to be in the fusiform area (e.g., Watanabe et al., 1999; Halgren et al., 2000). In fMRI studies with a very high spatial resolution, the fusiform face area (FFA) was identified in the fusiform gyrus and is selectively involved in the perception of faces (Kanwisher et al., 1997).

Recognizing facial movements is very important in addition to recognizing the static face, for example, in social communication and non-language communication. In monkeys, a neuronal population in the anterior superior temporal polysensory area (STPa) responded selectively to the motion of animate objects, including bodies and faces (e.g., Oram and Perrett, 1994). In recent monkey studies using fMRI, a different pattern was observed in the anterior STS, which responded more to dynamic than static faces, but was not sensitive to dot motion (Furl et al., 2012).

Many human studies have been conducted on static face perception, whereas a smaller number of studies have investigated facial movement perception using non-invasive neuroimaging methods. In fMRI studies, the STS in addition to MT/V5, which is considered to play an important role in motion perception, was activated when facial movement was viewed (Puce et al., 1998; Schultz and Pilz, 2009). In ERP studies (Wheaton et al., 2001; Thompson et al., 2002; Puce and Perrett, 2003; Puce et al., 2003), the amplitude of N170 in the occipitotemporal area evoked by facial (mouth) movement with a linedrawn face was larger than that evoked by general motion with a spatially "scrambled" line-drawn face. In a recent eventrelated spectral perturbations (ERSPs) study, occipitotemporal beta and gamma activities differentiated between facial and non-facial motion (Rossi et al., 2014). In our previous MEG studies (Watanabe et al., 2001; Miki et al., 2004, 2007), the activities by eyes and mouth movements of the occipitotemporal area, corresponding to human MT/V5, was different from movements in general. Consistent with our findings, a recent MEG study reported that cortical responses to eye blinks were clearly differently than expected on the basis of simple physical characteristics (Mandel et al., 2014). Moreover, Ulloa et al. (2014) reported that the initial gaze change elicited a significantly larger M170 under the deviated than mutual attention scenario.

Based on previous studies, we determined (1) whether processing of the perception of facial movements was specific and different from motion in general; (2) what information within the face was important to the processing of facial movement perception if this processing was specific; and (3) whether the right hemisphere played a more important role in facial movement perception than the left. Since we previously studied brain activities evoked by viewing various types of human facial movements based on our hypothesis, we herein introduced three representative studies: (a) Magnetoencephalographic study of occipitotemporal activity elicited by viewing mouth movements (Miki et al., 2004); (b) Effects of face contour and features on early occipitotemporal activity when viewing eye movement (Miki et al., 2007); and (c) Effects of inverting contour and features on processing for static and dynamic face perception: an MEG study (Miki et al., 2011).

#### **MAGNETOENCEPHALOGRAPHIC STUDY OF OCCIPITOTEMPORAL ACTIVITY ELICITED BY VIEWING MOUTH MOVEMENTS**

Facial movements are useful for social communication in humans. For example, the direction of the eye gaze is used to assess the social attention of others, and moreover, it becomes markedly easier to understand speech when we can see the mouth movements of the speaker. In a previous MEG study (Watanabe et al., 2001), a specific region for the perception of eye movements was detected within the occipitotemporal area, corresponding to human MT/V5, and this was different from motion in general. We hypothesized that the perception of the movement of facial parts may also be processed in a similar manner, which is different from motion in general. Therefore, the main objectives of the first study were to examine the temporal characteristics of the brain activity elicited by viewing mouth movements (opening and closing), and compare them to those of eye aversion movements and motion in general.

Seventeen right-handed adults (4 females, 13 males: 24–43 (mean age 32.2) years) with normal or corrected visual acuity participated in this study. We used apparent motion, in which the first stimulus (S1) was replaced by a second stimulus (S2) with no inter-stimulus interval as follows (**Figure 1**):


A large clear component, 1M, was elicited by all conditions (M-OP, M-CL, EYES, and RADIAL) within 200 ms of the stimulus onset (**Figure 2**). Concerning the peak latency of 1M, the means and standard deviations were 159.8 ± 17.3, 161.9 ± 15.0, 161.2 ± 18.9, and 140.1 ± 18.0 ms for M-OP, M-CL, EYES and RADIAL in the right hemisphere, respectively, and 162.4 ± 11.6, 160.9 ± 9.8, 164.6 ± 14.2, and 138.4 ± 9.0 ms for M-OP, M-CL, EYES, and RADIAL in the left, respectively. The latency for RADIAL was significantly shorter than that for the facial motion conditions (*p* < 0.05). No significant differences were observed in 1M latency between M-OP, M-CL, or EYES.

We used a multi-dipole model, brain electric source analysis (BESA; Scherg and Buchner, 1993) (Neuroscan, McLean, VA) computation of theoretical source generators in a threelayer spherical head model, and estimated activity in the occipitotemporal area, the human MT/V5 area homologue, from 1M. The means and standard deviations of the dipole moment of the estimated dipoles from 1M was 7.9 ± 1.9, 7.8 ± 3.2, 10.0 ± 6.8, and 13.8 ± 4.9 nAm for M-OP, M-CL, EYES, and RADIAL in the right hemisphere, respectively, and 7.4 ± 2.8, 6.7 ± 3.0, 9.3 ± 4.3, and 13.6 ± 1.8 nAm for M-OP,

**superimposed display for all conditions in a representative subject (adopted from Miki et al., 2004).**

M-CL, EYES, and RADIAL in the left, respectively. No significant differences were observed in the dipole moment (strength) for M-OP, M-CL and EYES between either hemisphere. However, M-OP and M-CL were significantly smaller than RADIAL (*p* < 0.05) in the right hemisphere and M-OP, M-CL, and EYES were significantly smaller than RADIAL (*p* < 0.05) in the left.

The results of the first study indicated that the occipitotemporal (human MT/V5) area was active in the perception of both mouth and eye movement. Furthermore, viewing mouth and eye movements did not elicit significantly different activity in the occipitotemporal (human MT/V5) area, which suggested that the perception of the movement of facial parts may be processed in the same manner, and this is different from motion in general.

#### **EFFECTS OF FACE CONTOUR AND FEATURES ON EARLY OCCIPITOTEMPORAL ACTIVITY WHEN VIEWING EYE MOVEMENT**

The first study showed that the perception of the movement of facial parts may be processed in the same manner, and this is different from motion in general. However, the main factor(s) causing differences in recognizing facial versus general movement have yet to be elucidated in detail.

Many studies have investigated effect facial contour and features using a static face. A previous study reported that it took longer to recognize an eyes-only stimulus or only facial features (eyes, nose, and mouth) than a full-face stimulus with a contour (e.g., Watanabe et al., 1999), and the contour of the face is important in static face recognition. However, to the best of our knowledge, the effects of the facial contour and features on facial movement recognition have not yet been investigated. Therefore, the main objectives of the second study were to investigate the effects of facial contour and features on early occipitotemporal activity evoked by facial movement. We also used a schematic face because a simple schematic drawing with a circle for a contour, two dots for eyes, and a straight line for lips, was recognized as a face even though each individual component of the drawing by itself was not. Previous studies using a schematic face showed that N170 was evoked by schematic faces as well as photographs of real faces (Sagiv and Bentin, 2001; Latinus and Taylor, 2006).

Thirteen right-handed adults (6 females, 7 males: 24–46 (means; 33.6) years) with normal or corrected visual acuity participated in this study. We used apparent motion and presented the following four conditions (**Figure 3**):


Subjects described the simple movement of dots for D, whereas eye movement for CDL, CD, and DL, though movement modalities were the same for all conditions. In source modeling, we used a single equivalent current dipole (ECD) model (Hämäläinen et al., 1993) within 145–220 ms of the stimulus onset.

Clear MEG responses were elicited in all conditions (CDL, CD, DL, and D) at the sensors in the bilateral occipitotemporal area (**Figure 4**). The means and standard deviations of the peak latency of the estimated dipole was 179.3 ± 26.3, 183.0 ± 16.9, 180.9 ± 20.8, and 180.3 ± 23.7 ms for CDL, CD, DL, and D in the right hemisphere, respectively, and 180.2 ± 14.9, 180.5 ± 24.8, 174.0 ± 24.9, and 177.7 ± 20.3 ms for CDL, CD, DL, and D in the left, respectively. No significant differences were observed among any condition.

The means and standard deviations of the dipole moment was 14.4 ± 6.2, 11.2 ± 7.9, 9.6 ± 6.5, and 10.3 ± 5.3 nAm for CDL, CD, DL, and D in the right hemisphere, respectively, and 12.7 ± 6.7, 11.1 ± 6.1, 9.6 ± 5.4, and 8.9 ± 5.5 nAm for CDL, CD, DL, and D in the left, respectively. The moment was significantly larger for CDL than for CD (*p* < 0.05), DL (*p* < 0.01), and D(*p* < 0.01) in the right hemisphere, and for CDL than for DL and D (*p* < 0.01) in the left hemisphere.

Our results in the second study demonstrated specific information processing for eye movements, which was different from motion in general, and activity in the occipitotemporal (human MT/V5) area related to this processing was influenced by whether movements appeared with the contour and/or features of the face.

#### **EFFECTS OF INVERTING CONTOUR AND FEATURES ON PROCESSING FOR STATIC AND DYNAMIC FACE PERCEPTION: AN MEG STUDY**

The second study showed that the activity evoked in the occipitotemporal area by eye movements was influenced by the existence of the contour and/or features of the face. However, it remained unclear whether this activity was influenced by the orientation of the contour and/or features of the face.

In static face perception, the N170 component was longer and larger for an inverted face than for upright face (Bentin et al., 1996; Sagiv and Bentin, 2001; Itier and Taylor, 2004; Latinus and Taylor, 2006), which indicated that N170 was affected by inversion of the face, i.e., the face inversion effect. In addition, the latency of N170 was longer for scrambled features than for upright faces (George et al., 1996; Latinus and Taylor, 2006), which confirmed that N170 was affected by a disruption in the spatial relation between the facial contour and features.

Based on the findings of previous N170 studies on static face perception, we hypothesized that the perception of eye movements may mainly be affected by information on the contour and other facial features. Therefore, the main objectives of the third study were to investigate the effects of inverting the facial contour and features on the occipitotemporal (human MT/V5) area related to a dynamic face and what information within the face was important for processing dynamic face perception. We also investigated the effects of inverting the facial contour and features of the face on the fusiform area related static face perception to compare with the occipitotemporal area.

indicates the main response after the S2 onset (adopted from Miki et al.,

We recruited 10 right-handed adults (3 females and 7 males: 24–47 (means; 30.6) years) with normal or corrected visual acuity. We used apparent motion and presented the following three conditions (**Figure 5**):

2007).


The eyes were averted to the right of the viewer under all conditions.

As in the second study, we used a single ECD model (Hämäläinen et al., 1993) and estimated the dipole within 105– 200 ms of the S1 onset (static face) and 115–210 ms of the S2 onset (eye movements).

In static face perception (S1 onset), ECDs were estimated to lie in the fusiform area from MEG following S1 in all conditions.

The means and standard deviations of the peak latency in activity in the fusiform area was 133.6 ± 18.5, 148.4 ± 22.4, and 151.4 ± 24.1 ms for U&U, U&I, and I&I in the right hemisphere, respectively, and 143.2 ± 19.7, 162.2 ± 21.0, and 148.4 ± 13.6 ms for U&U, U&I, and I&I in the left, respectively. Latency was significantly longer for U&I (Upright contour and Inverted features) (*p* < 0.05) and I&I (Inverted contour and Inverted features) (*p* < 0.05) than for U&U in the right hemisphere, and also for U&I than for U&U (*p* < 0.01) and I&I (*p* < 0.05) in the left (**Figure 6**).

The means and standard deviations of the strength (the maximum of the dipole moment) of activity in the fusiform area was 28.6 ± 21.1, 35.5 ± 18.5, and 34.4 ± 18.5 nAm for U&U, U&I, and I&I in the right hemisphere, respectively, and 21.9 ± 13.5, 20.3 ± 11.8, and 21.9 ± 14.5 nAm for U&U, U&I, and I&I in the left, respectively. No significant differences were observed in the maximum of the dipole moment among the three conditions.

In dynamic face perception (S2 onset), ECDs were estimated to lie in the occipitotemporal area, the human MT/V5 area homologue, from MEG following S2 in all conditions.

The means and standard deviations of the peak latency of activity in the occipitotemporal area was 163.8 ± 31.3, 159.4 ± 22.3, and 159.7 ± 25.1 ms for U&U, U&I, and I&I in the right hemisphere, respectively, and 151.6 ± 22.1, 157.8 ± 24.4, and 151.7 ± 23.8 ms for U&U, U&I, and I&I in the left, respectively. No significant differences were observed in the peak latency among the three conditions.

The means and standard deviations of the strength (the maximum of the dipole moment) of activity in the occipitotemporal area was 11.3 ± 4.6, 11.5 ± 4.8, and 15.1 ± 5.9 nAm for U&U, U&I, and I&I in the right hemisphere, respectively, and 9.6 ± 3.7, 10.0 ± 5.6, and 8.8 ± 4.6 nAm for U&U, U&I, and I&I in the left,

respectively. The maximum of the dipole moment was larger for I&I than for U&U and U&I in the right hemisphere (*p* < 0.01), but not the left (**Figure 7**).

The results of the third study indicated the following: (a) considering the fusiform area related to static face perception, activity was affected more by the inversion of features in the right hemisphere while it was affected more by a disruption in the spatial relationship between the facial contour and features in the left hemisphere; and (b) considering the occipitotemporal (human MT/V5) area related to dynamic face perception, activity was affected by the inversion of the facial contour in the right, but not in the left hemisphere.

#### **SUMMARY AND CONCLUSION**

In our three studies, we focused on activity in the occipitotemporal area, the human MT/V5 homologue, related to facial parts movement. We summarized our results as followings: (1) viewing mouth and eye movements did not elicit significantly different activity in the occipitotemporal area; (2) neuronal activities in the occipitotemporal area evoked by facial (eye) movements were affected by whether the contour and/or features of the face were in the stimulus; (3) the activity of the right occipitotemporal area evoked by facial (eye) movements was affected by the inversion of the facial contour, and these results indicated the following: (1) processing of the perception on facial movements is specific and is different from motion in general; (2) the existence of the facial contour and face parts are important factors in the perception of facial movements; (3) the orientation of the contour and spatial relationship between the contour and facial parts are also important; and (4) the right occiptiotemporal area is more important in the perception of the facial movements than the left. Based on the results of the three experiments, it still remains unclear how the transmission of facial movement processing was modulated by facial form information.

Connectivity models that modeled communication between the ventral form and dorsal motion pathways were tested in a recent fMRI study related to perception of dynamic faces (Furl et al., 2014), and the findings obtained clearly showed that facial form information modulated the transmission of motion from V5 to the STS. Based on these findings, we hypothesized that information on the facial contour and parts, transmitting via FFA and/or OFA (occipital face area), may gate the transmission of information regarding facial motion via MT/V5, and that facial form and motion information may have been integrated in the STS. We consider these results and hypothesis to be useful for investigating the functional roles of human brain connectomes and also provide an insight into facial movement processes.

#### **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: 27 April 2014; accepted: 07 July 2014; published online: 29 July 2014*. *Citation: Miki K and Kakigi R (2014) Magnetoencephalographic study on facial movements. Front. Hum. Neurosci. 8:550. doi: 10.3389/fnhum.2014.00550 This article was submitted to the journal Frontiers in Human Neuroscience*.

*Copyright © 2014 Miki and Kakigi. 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*.