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

^{1}The KEY Institute for Brain-Mind Research, University of Zurich, Zurich, Switzerland^{2}Department of Neuropsychiatry, Kansai Medical University, Osaka, Japan^{3}CIMFAV, Universidad de Valparaiso, Valparaiso, Chile^{4}Neuroinformatic Department, Cuban Neuroscience Center, Havana, Cuba^{5}Department of Psychiatry, Shiga University of Medical Science, Shiga, Japan^{6}Division of Cerebral Integration, National Institute for Physiological Sciences, Okazaki, Japan

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.

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

**Figure 1. Direct causal directed connections between nodes, corresponding to the toy example defined by the multivariate autoregressive model in Equation 13**.

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

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:

1. The noise contribution ratio (NCR) of Akaike (1968). This work has apparently gone unnoticed by most researchers in this field. It has been extensively used and published under other names, in particular by Saito and Harashima (1981), Kaminski and Blinowska (1991), and Baccala et al. (1998). Akaike's NCR method discovers connections (direct and indirect, without distinction), their directionality, and their spectral characteristics.

2. The partial directed coherence (PDC) of Baccala and Sameshima (2001), which is a measure designed to quantify direct connections which are not confounded by indirect paths, their directionality and their spectral characteristics. This is a very widely used measure (cited 672 times at the time of this writing according to “Google-Scholar”), practically considered the “golden standard” when all network properties (and not just part of them) are of interest.

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

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

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

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

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

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

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

## Methods

### The Multivariate Autoregressive Model

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

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

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

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

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

The frequency domain representation is:

where *X* (ω) ∈ ℂ^{q×1}, *A* (ω) ∈ ℂ^{q×q}, ∈ (ω) ∈ ℂ^{q×1} are the discrete Fourier transforms, and where ω denotes discrete frequency.

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

with :

where the superscript “*” denotes matrix transpose and complex conjugate, the superscript “−1” denotes matrix inversion, **I** is the identity matrix, and **S**_{ε} ∈ ℝ^{q×q} is the noise covariance.

### The Isolated Effective Coherence (iCoh)

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

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

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

and :

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

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

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

which clearly satisfies:

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:

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

The gPDC is:

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 *X _{cki}*, 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:

where Φ_{t} ∈ ℝ^{NE ×1} denotes the instantaneous scalp electric potential at *N _{E}* electrodes,

**J**

_{t}∈ ℝ

^{(3NV)×1}is the instantaneous current density vector field at

*N*cortical voxels (consisting of three components at each voxel), and

_{V}**K**∈ ℝ

^{NE ×(3NV)}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:

where **T** ∈ ℝ^{(3NV)×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 three-dimensional 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 **Ĵ**_{t} in Equation 12) consists of three time series, one for each moment component of the current density vector (i.e., dipole). In practice, this can be reduced to a single time series, due to the fact that the current density vector is anatomically constrained to have an orientation orthogonal to the cortical surface (see e.g., Baillet et al., 2001). Under this assumption, the 3 × 3 covariance matrix for the current density vector at each voxel must have rank 1, with the dipole orientation given by its largest eigenvector (Mosher et al., 1992; Mosher and Leahy, 1998). This procedure is applied in this study for the estimation of single time series of electric neuronal activity at each voxel. Note that this estimator for the current density vector field orientation is a maximum variance estimator.

## Materials

### A Toy Example: Five Time Series

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

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:

(1) Node 1 in isolation oscillates at peak frequency 28 Hz.

(2) Node 2 in isolation oscillates at peak frequency 16 Hz.

(3) Nodes 3, 4, and 5 in isolation oscillate at peak frequency 28 Hz.

(4) Nodes 3, 4, and 5 are receiving identical information from node 2.

### 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:

(1) Superior frontal gyrus (left), BA 10: *X* = −25, *Y* = 65, *Z* = −5

(2) Middle Occipital Gyrus (right), BA 18: *X* = 20, *Y* = −100, *Z* = 5

(3) Post-central Gyrus (left), BA 3: *X* = −50, *Y* = −20, *Z* = 60

(4) Middle Temporal Gyrus (left), BA 21: *X* = −65, *Y* = −15, *Z* = −15

(5) Middle Temporal Gyrus (right), BA 21: *X* = 70, *Y* = −20, *Z* = −10

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

Figure 2 illustrates the five cortical locations.

**Figure 2. Schematic representation of the anatomical locations of five cortical point sources used for generating EEG**.

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:

(1) Biological noise, where independent and identically distributed uniform [0, 0.2] random values were assigned to the current density at each cortical voxel and at each time sample.

(2) Measurement noise, where independent and identically distributed uniform [0, 0.2] random values were multiplied by the potential at each electrode and time sample, and added to potential.

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

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

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.

**Figure 4. Estimated connectivity properties for simulated EEG signals (see Figure 2).** Isolated effective coherence (iCoh) shown in RED, and the generalized partial directed coherence (gPDC) shown in BLUE. 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.

Of importance to note in Figure 4:

(1) Ideally, Figure 4 should be identical to Figure 3. This is the case to a very good approximation from a qualitative point of view, despite the use of estimated signals using eLORETA, from as few as 19 electrodes, and corrupted with relative high levels of biological and measurement noise.

(2) The two methods give very different results with respect to node #2 as sender (column 2).

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

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

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.

**Figure 6. Comparison of electric neuronal activity (eLORETA) between eyes open and closed conditions.** A Log-F-ratio statistic with correction for multiple testing was used, with corrected *p* = 0.05 at LogF = 0.91. Eyes open is characterized by significantly stronger activity in frontal cortical regions oscillating at 3 Hz **(A)** and in the beta band 23–28 Hz **(C)**. Eyes closed is characterized by significantly stronger activity in occipital cortical regions oscillating at 10 Hz **(B)**.

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 regions of interest are:

(1) Anterior Cingulate, BA 32: *X* = 0, *Y* = 45, *Z* = 10

(2) Posterior Cingulate, Precuneus, BAs 23, 31: *X* = 0, *Y* = −50, *Z* = 30

(3) Inferior Parietal Lobule (left), BA 40: *X* = −45, *Y* = −45, *Z* = 50

(4) Inferior Parietal Lobule (right), BA 40: *X* = 45, *Y* = −45, *Z* = 50

(5) Dorsolateral Pre-frontal (left), BA 10: *X* = −40, *Y* = 40, *Z* = 25

(6) Dorsolateral Pre-frontal (right), BA 10: *X* = 40, *Y* = 40, *Z* = 25

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 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”).

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

## 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 “sender-receiver” 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, 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 non-parametric 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 spatio-temporal 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.

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.

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

## 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

Akaike, H. (1968). On the use of a linear model for the identification of feedback systems. *Ann. Inst. Stat. Math*. 20, 425–439. doi: 10.1007/BF02911655

Akaike, H. (1974). A new look at the statistical model identification. *Autom. Control IEEE Trans*. 19, 716–723. doi: 10.1109/TAC.1974.1100705

Baccala, L., Sameshima, K., Ballester, G., Do Valle, A., and Timo-Iaria, C. (1998). Studying the interaction between brain structures via directed coherence and Granger causality. *Appl. Signal Process*. 5, 40. doi: 10.1007/s005290050005

Baccala, L. A., and Sameshima, K. (2001). Partial directed coherence: a new concept in neural structure determination. *Biol. Cybern*. 84, 463–74. doi: 10.1007/PL00007990

Baccalá, L. A., Sameshima, K., and Takahashi, D. (2007). “Generalized partial directed coherence,” in *15th International Conference on Digital Signal Processing IEEE* (Wales). 163–166. doi: 10.1109/ICDSP.2007.4288544

Baillet, S., Mosher, J. C., and Leahy, R. M. (2001). Electromagnetic brain mapping. *Signal Process. Mag. IEEE* 18, 14–30. doi: 10.1109/79.962275

Blair, R. C., and Karniski, W. (1994). *Distribution-Free Statistical Analyses of Surface and Volumetric Maps, Functional Neuroimaging: Technical Foundations*, San Diego, CA: Academic Press.

Brillinger, D. R. (2001). *Time Series: Data Analysis and Theory*. Philadelphia, PA: Siam. doi: 10.1137/1.9780898719246

Chen, Y., Bressler, S. L., and Ding, M. (2006). Frequency decomposition of conditional Granger causality and application to multivariate neural field potential data. *J. Neurosci. Methods* 150, 228–237. doi: 10.1016/j.jneumeth.2005.06.011

Faes, L., Erla, S., Porta, A., and Nollo, G. (2013). A framework for assessing frequency domain causality in physiological time series with instantaneous effects. *Philos. Trans. A Math. Phys. Eng. Sci*. 371:20110618. doi: 10.1098/rsta.2011.0618

Frei, E., Gamma, A., Pascual-Marqui, R., Lehmann, D., Hell, D., and Vollenweider, F. X. (2001). Localization of MDMA−induced brain activity in healthy volunteers using low resolution brain electromagnetic tomography (LORETA). *Hum. Brain Mapp*. 14, 152–165. doi: 10.1002/hbm.1049

Fuchs, M., Kastner, J., Wagner, M., Hawes, S., and Ebersole, J. S. (2002). A standardized boundary element method volume conductor model. *Clin. Neurophysiol*. 113, 702–712. doi: 10.1016/S1388-2457(02)00030-5

Geweke, J. F. (1984). Measures of conditional linear dependence and feedback between time series. *J. Am. Stat. Assoc*. 79, 907–915. doi: 10.1080/01621459.1984.10477110

Goldberger, A. L., Amaral, L. A., Glass, L., Hausdorff, J. M., Ivanov, P. C., Mark, R. G., et al. (2000). PhysioBank, PhysioToolkit, and PhysioNet: components of a new research resource for complex physiologic signals. *Circulation* 101, E215–E220. doi: 10.1161/01.CIR.101.23.e215

Gomez-Herrero, G., Atienza, M., Egiazarian, K., and Cantero, J. L. (2008). Measuring directional coupling between EEG sources. *Neuroimage* 43, 497–508. doi: 10.1016/j.neuroimage.2008.07.032

Granger, C. W. (1969). Investigating causal relations by econometric models and cross-spectral methods. *Econometrica* 37, 424–438. doi: 10.2307/1912791

Helmholtz, H. (1853). Ueber einige Gesetze der Vertheilung elektrischer Ströme in körperlichen Leitern mit Anwendung auf die thierisch-elektrischen Versuche. *Ann. Phys*. 165, 211–233. doi: 10.1002/andp.18531650603

Jao, T., Vertes, P. E., Alexander-Bloch, A. F., Tang, I. N., Yu, Y. C., Chen, J. H., et al. (2013). Volitional eyes opening perturbs brain dynamics and functional connectivity regardless of light input. *Neuroimage* 69, 21–34. doi: 10.1016/j.neuroimage.2012.12.007

Jiménez, J. C., Biscay, R., and Montoto, O. (1995). Modeling the electroencephalogram by means of spatial spline smoothing and temporal autoregression *Biol. Cybern*. 72, 249–259. doi: 10.1007/BF00201488

Jurcak, V., Tsuzuki, D., and Dan, I. (2007). 10/20, 10/10, and 10/5 systems revisited: their validity as relative head-surface-based positioning systems. *Neuroimage* 34, 1600–1611. doi: 10.1016/j.neuroimage.2006.09.024

Kaminski, M., and Blinowska, K. (1991). A new method of the description of the information flow in the brain structures. *Biol. Cybern*. 65, 203–210. doi: 10.1007/BF00198091

Karniski, W., Blair, R. C., and Snider, A. D. (1994). An exact statistical method for comparing topographic maps, with any number of subjects and electrodes. *Brain Topogr*. 6, 203–210. doi: 10.1007/BF01187710

Klimesch, W. (2012). Alpha-band oscillations, attention, and controlled access to stored information. *Trends Cogn. Sci*. 16, 606–617. doi: 10.1016/j.tics.2012.10.007

Knyazev, G. G., Slobodskoj-Plusnin, J. Y., Bocharov, A. V., and Pylkova, L. V. (2011). The default mode network and EEG alpha oscillations: an independent component analysis. *Brain Res*. 1402, 67–79. doi: 10.1016/j.brainres.2011.05.052

Lancaster, J. L., Woldorff, M. G., Parsons, L. M., Liotti, M., Freitas, C. S., Rainey, L., et al. (2000). Automated Talairach atlas labels for 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

Lehmann, D., Faber, P. L., Gianotti, L. R., Kochi, K., and Pascual-Marqui, R. D. (2006). Coherence and phase locking in the scalp EEG and between LORETA model sources, and microstates as putative mechanisms of brain temporo-spatial functional organization. *J. Physiol. Paris* 99, 29–36. doi: 10.1016/j.jphysparis.2005.06.005

Lehmann, D., Faber, P. L., Tei, S., Pascual-Marqui, R. D., Milz, P., and Kochi, K. (2012). Reduced functional connectivity between cortical sources in five meditation traditions detected with lagged coherence using EEG tomography. *Neuroimage* 60, 1574–1586. doi: 10.1016/j.neuroimage.2012.01.042

Marinazzo, D., Liao, W., Chen, H., and Stramaglia, S. (2011). Nonlinear connectivity by Granger causality. *Neuroimage* 58, 330–338. doi: 10.1016/j.neuroimage.2010.01.099

Mazziotta, J., Toga, A., Evans, A., Fox, P., Lancaster, J., Zilles, K., et al. (2001). A probabilistic atlas and reference system for the human brain: international consortium for brain mapping (ICBM). *Philos. Trans. R. Soc. Lond. B Biol. Sci*. 356, 1293–1322. doi: 10.1098/rstb.2001.0915

Mosher, J. C., and Leahy, R. M. (1998). Recursive MUSIC: a framework for EEG and MEG source localization. *Biomed. Eng. IEEE Trans*. 45, 1342–1354. doi: 10.1109/10.725331

Mosher, J. C., Lewis, P. S., and Leahy, R. M. (1992). Multiple dipole modeling and localization from spatio-temporal MEG data. *Biomed. Eng. IEEE Trans*. 39, 541–557. doi: 10.1109/10.141192

Nichols, T. E. (2012). Multiple testing corrections, nonparametric methods, and random field theory. *Neuroimage* 62, 811–815. doi: 10.1016/j.neuroimage.2012.04.014

Nichols, T. E., and Holmes, A. P. (2002). Nonparametric permutation tests for functional neuroimaging: a primer with examples. *Hum. Brain Mapp*. 15, 1–25. doi: 10.1002/hbm.1058

Nolte, G., Ziehe, A., Nikulin, V. V., Schlogl, A., Kramer, N., Brismar, T., et al. (2008). Robustly estimating the flow direction of information in complex physical systems. *Phys. Rev. Lett*. 100:234101. doi: 10.1103/PhysRevLett.100.234101

Pascual-Marqui, R. (2007). Discrete, 3D distributed, linear imaging methods of electric neuronal activity. Part 1: exact, zero error localization. Avilable online at: http://arxiv.org/pdf/0710.3341

Pascual-Marqui, R., Biscay, R., Bosch-Bayard, J., Lehmann, D., Kochi, K., Yamada, N., et al. (2014). Isolated effective coherence (iCoh): causal information flow excluding indirect paths. arXiv preprint arXiv:1402.4887.

Pascual-Marqui, R. D. (2002). Standardized low-resolution brain electromagnetic tomography (sLORETA): technical details. *Methods Find Exp. Clin. Pharmacol*. 24(Suppl. D), 5–12.

Pascual-Marqui, R. D. (2009). “Theory of the EEG inverse problem,” in *Quantitative EEG Analysis: Methods and Clinical Applications, Chapter 5*, eds S. Tong and N. V Thakor (Boston, MA: Artech House). 121–140.

Pascual-Marqui, R. D., Lehmann, D., Koukkou, M., Kochi, K., Anderer, P., Saletu, B., et al. (2011). Assessing interactions in the brain with exact low-resolution electromagnetic tomography. *Philos. Trans. A Math. Phys. Eng. Sci*. 369, 3768–3784. doi: 10.1098/rsta.2011.0081

Pascual-Marqui, R. D., Michel, C. M., and Lehmann, D. (1994). Low resolution electromagnetic tomography: a new method for localizing electrical activity in the brain. *Int. J. Psychophysiol*. 18, 49–65. doi: 10.1016/0167-8760(84)90014-X

Plomp, G., Quairiaux, C., Michel, C. M., and Astolfi, L. (2014). The physiological plausibility of time-varying Granger-causal modeling: normalization and weighting by spectral power. *Neuroimage*. 97, 206–216. doi: 10.1016/j.neuroimage.2014.04.016

Radhakrishna Rao, C. (1981). A lemma on g-inverse of a matrix and computation of correlation coefficients in the singular case. *Commun. Stat.Theor. Methods* 10, 1–10. doi: 10.1080/03610928108828015

Saito, Y., and Harashima, H. (1981). “Tracking of information within multichannel {EEG} record causal analysis in {EEG},” in *Recent Advances in {EEG} and {EMG} Data Processing*, eds N. Yamaguchi and K. Fujisawa (Amsterdam: Elsevier). 133–146.

Schalk, G., McFarland, D. J., Hinterberger, T., Birbaumer, N., and Wolpaw, J. R. (2004). BCI2000: a general-purpose brain-computer interface (BCI) system. *IEEE Trans. Biomed. Eng*. 51, 1034–1043. doi: 10.1109/TBME.2004.827072

Schelter, B., Timmer, J., and Eichler, M. (2009). Assessing the strength of directed influences among neural signals using renormalized partial directed coherence. *J. Neurosci. Methods* 179, 121–130. doi: 10.1016/j.jneumeth.2009.01.006

Sigala, R., Haufe, S., Roy, D., Dinse, H. R., and Ritter, P. (2014). The role of alpha-rhythm states in perceptual learning: insights from experiments and computational models. *Front. Comput. Neurosci*. 8:36. doi: 10.3389/fncom.2014.00036

Valdes-Sosa, P. A., Roebroeck, A., Daunizeau, J., and Friston, K. (2011). Effective connectivity: influence, causality and biophysical modeling. *Neuroimage* 58, 339–361. doi: 10.1016/j.neuroimage.2011.03.058

Valdes-Sosa, P. A., Sanchez-Bornot, J. M., Lage-Castellanos, A., Vega-Hernandez, M., Bosch-Bayard, J., Melie-Garcia, L., et al. (2005). Estimating brain functional connectivity with sparse multivariate autoregression. *Philos. Trans. R. Soc. Lond. B Biol. Sci*. 360, 969–981. doi: 10.1098/rstb.2005.1654

Weiss, S. A., McKhann, G. Jr. Goodman, R., Emerson, R. G., Trevelyan, A., Bikson, M., et al. (2013). Field effects and ictal synchronization: insights from in homine observations. *Front. Hum. Neurosci*. 7:828. doi: 10.3389/fnhum.2013.00828

Wibral, M., Vicente, R., and Lizier, J. T. (2014). *Directed Information Measures in Neuroscience*. Heidelberg: Springer. doi: 10.1007/978-3-642-54474-3

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

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

Received: 08 May 2014; Accepted: 03 June 2014;

Published online: 20 June 2014.

Edited by:

Ryouhei Ishii, Osaka University Graduate School of Medicine, JapanReviewed by:

Guido Nolte, Fraunhofer FIRST, GermanyGennady Knyazev, Academy of Medical Sciences, Russia

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.

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