## Abstract

Scalp electric potentials (electroencephalogram; EEG) are contingent to the impressed current density unleashed by cortical pyramidal neurons undergoing post-synaptic processes. EEG neuroimaging consists of estimating the cortical current density from scalp recordings. We report a solution to this inverse problem that attains exact localization: exact low-resolution brain electromagnetic tomography (eLORETA). This non-invasive method yields high time-resolution intracranial signals that can be used for assessing functional dynamic connectivity in the brain, quantified by coherence and phase synchronization. However, these measures are non-physiologically high because of volume conduction and low spatial resolution. We present a new method to solve this problem by decomposing them into instantaneous and lagged components, with the lagged part having almost pure physiological origin.

## 1. Introduction

We aim to use the electroencephalogram (EEG), consisting of high time resolution measurements of scalp electric potential differences, for studying brain function. These measurements are employed for computing the current density on the cortex. The new data now consist of time series of current density, typically sampled at rates of 100–1000 samples per second, and virtually recorded from all over the cortex, at typically 5 mm spatial resolution. This massive amount of data contains essential information on brain function. This poses the problem of how to process and extract information.

Traditionally, topographic scalp maps have been used for studying brain function. Compelling studies that demonstrate the functional brain state information in scalp maps can be found in Lehmann [1], Saletu *et al*. [2] and John *et al*. [3]. However, their use for localization inference is questionable, since cortical electric neuronal activity does not map radially onto the scalp. We emphasize that this problem also applies to the magnetoencephalogram (MEG). This justifies the need for an inverse solution that correctly localizes cortical activity from scalp potentials.

The aims of this study are: (i) to present an inverse solution to the EEG problem that computes cortical current density (i.e. electric neuronal activity), with optimal localization properties, and (ii) to present methods for properly estimating dynamic functional connectivity in the brain based on the estimated current density signals.

## 2. Electroencephalograms and intracranial current density

### (a) General formulation

The intracranial sources of scalp electric potential differences are thought to be cortical pyramidal neurons, which undergo post-synaptic potentials that create an active, impressed current density. The dipolar nature (current density vector) of these sources is documented in Mitzdorf [4] and Martin [5].

The corresponding relation between the current density vector field and the scalp potentials can be expressed as [6]
2.1
where **Φ**_{c}∈R^{NE×1} denotes the potentials at *N*_{E} scalp electrodes, all measured with respect to the same reference electrode; **J**∈R^{NV×1} is the current density at *N*_{V} cortical voxels; *c* is a scalar determined by the reference electrode (which can be arbitrary); **1**∈R^{NE×1} is a vector of ones; and **K**_{c}∈R^{NE×NV} is the lead field (determined by the geometry and conductivity profile of the head). The subscript ‘c’ in **Φ**_{c} and **K**_{c} emphasizes a fact of nature: potentials are determined up to an arbitrary additive constant (which in this case is related to the choice for the reference electrode). Equation (2.1) expresses a deterministic law of physics. Further on, this will be embedded in a stochastic setting, by considering both measurement and biological noise.

Technical details on the definition and computation of the lead field can be found in Sarvas [7] and Fuchs *et al*. [8]. In previous simple simulation studies, such as those reported in Pascual-Marqui [9], a spherical head model was used, for which analytical expressions exist for the lead field [10]. In this study, the lead field was computed numerically for a realistic head shape, using the boundary-element method [8].

This forward equation corresponds to an instantaneous discrete sampling of the measurement space (scalp electrodes) and the solution space (cortical voxels). For simplicity, we assume that the orientation of the current density vector is known (the general case with full unknown current density vector can be found in Pascual-Marqui [6]).

The inverse problem of interest is the estimation of the unknown electric neuronal activity **J**, given the lead field **K**_{c} and the scalp potential measurements **Φ**_{c}. The nuisance parameter *c*, related to the arbitrary choice of the reference electrode, must be accounted for.

### (b) The final solution to the so-called ‘reference electrode problem’

We first solve the so-called ‘reference electrode problem’, which corresponds to solving
2.2
where ∥•∥ denotes the Euclidean norm. This means that *c* must satisfy, as best as possible, the forward equation (2.1). The squared Euclidean distance can be generalized to a Mahalanobis distance using the covariance matrix of the measurement noise. Solving the problem in equation (2.2) and plugging the solution back into equation (2.1) gives
2.3
where **H** is
2.4
where the superscript ‘T’ denotes the vector–matrix transpose; and with **I**∈R^{NE×NE} being the identity matrix. The essential property of **H** is
2.5
The matrix **H** is known in statistics as the centring matrix (e.g. [11]), which subtracts the average value.

This has a profound implication, demonstrating that inverse solutions are reference independent as they must be obtained from the reference-independent forward equation (2.3), which is invariant to any change of reference.

Henceforth, the reference-independent forward equation (2.3) will be written as
2.6
with
2.7
and
2.8
The effect of the operator **H** is known in the EEG literature as the average reference [6,12].

An important consequence of these initial derivations is the proof that there is no such thing as an ‘ideal reference’, if the aim is localization. Obviously, an inverse solution will change with different samplings (i.e. with the number and the locations of scalp electrodes); but for a given sampling, the solution to a properly posed inverse problem that explicitly models the reference electrode (as in equations (2.1) and (2.2)) will not depend on the choice of the reference.

### (c) The inverse solution

We seek to solve the linear system of equations (2.6) for the unknown current density **J**. From the outset, we clarify that this is a multi-variate problem, analogous, for example, to X-ray tomography [13], which is not solved as a collection of independent univariate least squares for each voxel separately (repeated single dipole fitting).

Typically, the number of electrodes is much smaller than the number of voxels, *N*_{E}≪*N*_{V}, and the problem is underdetermined, with infinitely many solutions. In a general setting, Helmholtz [14] demonstrated that many different current density distributions exist in three-dimensional space that are consistent with the electric potential distribution on a surface enclosing the volume.

Similar to the setting in equation (2.2), we look for particular linear solutions of the form
2.9
which is a regularized, weighted minimum norm problem, where *α*≥0 is the Tikhonov regularization parameter [15]; and **W**∈R^{NV×NV} is a symmetric positive definite weight matrix (see Pascual-Marqui [16] for a generalization to non-negative definite matrices). It is important to note that the problem in equation (2.9) has at least two different interpretations. On the one hand, it is a conventional form studied in mathematical functional analysis [15]. On the other hand, it can be derived from a Bayesian formulation of the inverse problem [17]. An excellent general review of other functionals for solving the inverse problem can be found in Valdes-Sosa *et al*. [18].

The general solution to the problem in equation (2.9) is 2.10 where the superscript ‘+’ denotes the Moore–Penrose pseudoinverse [19].

Setting **W**=**I** (identity matrix) gives the classical non-weighted minimum norm solution [20]
2.11
It was shown in Pascual-Marqui [9], both empirically and theoretically, that the non-weighted minimum norm has very bad localization properties, misplacing deep sources to the surface. The reason is because this solution is a harmonic function that attains its extreme values only at the boundary of the solution space [21]. This basic physics property of the minimum norm invalidates its use for localization, regardless of the fact that it is the simplest solution (see [22]).

In the method known as low-resolution electromagnetic tomography (LORETA) [23], the matrix **W** implements the squared three-dimensional spatial Laplacian operator. In its simplest discrete implementation, the Laplacian compares the current density at one voxel with that of its closest neighbours, as explained in detail in Pascual-Marqui *et al*. [23] and Pascual-Marqui [9]. Minimization of its squared value forces the current density at each voxel to be as similar as possible to that of its neighbours, thus forcing spatial smoothness. This condition is a macroscopic implementation of what occurs at the cellular level: non-negligible EEG is only possible if neighbouring pyramidal neurons are very highly synchronized (e.g. [24]). It should be noted that a three-dimensional Laplacian does not respect smoothness along the cortical surface with its many foldings (e.g. smoothing would occur across opposite walls of a sulcus). A more accurate implementation would consider the Laplacian along the cortical manifold, as in Grova *et al*. [25], where the Laplacian was restricted to the cortical mesh.

It was empirically shown in simulation studies [9], under ideal noiseless conditions, using 148 electrodes and 818 uniformly distributed voxels that LORETA achieves an average localization error of only one voxel unit, uniformly at all depths, when compared with the non-weighted minimum norm that has maximal error for deep sources.

It has been thought that better localization can be achieved with depth-weighted minimum norm, by giving larger weights to deep voxels. A recent attempt in this direction can be found in Lin *et al*. [26], where some improvement is reported. However, the weights only achieve a modest reduction in localization error compared with the minimum norm, as reported in Pascual-Marqui [16].

### (d) Standardized current density tomographies

In this approach, after the current density is estimated, a second post-processing step is applied, consisting of the statistical standardization of the current density values. In this sense, localization inference is not based on the estimated current density directly, but on its standardized value, which is unitless.

Typically, these methods use the minimum norm solution as the current density estimator in the first step, which by itself has a large localization error. The statistically standardized current density is defined as
2.12
where [**J**]_{i} is the estimated current density at the *i*th voxel (e.g. from equation (2.11)); **S**_{J}∈R^{NV×NV} is the covariance matrix for the current density; and [**S**_{J}]_{ii} is its *i*th diagonal element corresponding to the variance at the *i*th voxel. Note that localization inference is based on the images of unitless values *z*_{i}, and not on the current density [**J**]_{i}.

There is no unique methodology for selecting the standard deviation [**S**_{J}]_{ii} of the current density at each voxel. In the Bayesian formulation used by Dale *et al*. [27], it is assumed that the standard deviation for the current density at each voxel is due exclusively to measurement noise, i.e. only measurement noise contributes to the total EEG covariance. Their method is known as dynamic statistical parametric mapping (dSPM).

A very different derivation for the standard deviation is given by Pascual-Marqui [28], using a functional analysis formulation. This method is known as standardized low-resolution electromagnetic tomography (sLORETA). Unlike the Dale *et al*. [27] method, it is assumed in sLORETA that the total EEG covariance receives contributions from noise in the scalp measurements and from neuronal generator noise.

A detailed comparison of the linear imaging methods of dSPM, sLORETA and the non-weighted minimum norm solution was performed using a realistic head model under ideal (low measurement noise) conditions in Pascual-Marqui [28]. All possible single-point test sources (Dirac deltas) were used for the generation of the scalp potentials, which were then given to each imaging method. Localization errors were defined as the distance between the location of the absolute maximum value (|*z*_{i}| or |[**J**]_{i}|) and the actual location of the test source. The mean location errors for the minimum norm method and for the Dale *et al*. [27] method were 38 and 34 mm, respectively. This indicates that the standardization of the Dale *et al*. method produces only a slight improvement over the minimum norm solution. The sLORETA method showed zero localization error under ideal (no-noise) conditions. These first empirical results demonstrated the zero-error property of sLORETA.

Soon after, theoretical proof for the zero-error property of sLORETA under ideal conditions was independently provided by Sekihara *et al*. [29] and Greenblatt *et al*. [30]. It was later shown theoretically that sLORETA has no localization bias, even in the presence of both measurement and biological noise [31].

### (e) Exact low-resolution electromagnetic tomography

Historically, soon after the publication of the first tomography in this field in 1984, namely the non-weighted minimum norm solution [20], many attempts were made to improve the mislocalization of deep sources. A long sought solution has been to find an appropriate weight matrix **W** in equation (2.10), such that the distributed linear inverse solution has zero localization error when tested with point sources anywhere in the brain, under ideal (no-noise) conditions.

The weights for exact localization with zero error under ideal (no-noise) conditions are obtained from the following nonlinear system of equations:
2.13
where *w*_{i}, for *i*=1,…,*N*_{V}, are the elements of the diagonal weight matrix **W**, and **K**_{i}∈R^{NE×1} denotes the *i*th column of the lead field matrix **K**. Note that the weights depend nonlinearly on the lead field columns, but the inverse solution remains linear (equation (2.10)). Equation (2.13) corresponds to the algorithm for solving the weights (which do not depend on the measurements).

Initialize the diagonal weight matrix

**W**with elements*w*_{i}=1, for*i*=1,…,*N*_{V}.Compute 2.14

Holding

**C**fixed, for*i*=1,…,*N*_{V}, compute new weights 2.15Using the new weights, go to step 2 until convergence (i.e. stop when the change in the weight matrix is sufficiently small).

Theoretical proof follows for the exact localization property of the linear inverse solution. From equation (2.10), the current density at the *i*th voxel is
2.16
Substituting equation (2.13) into equation (2.16) gives
2.17
with **C** given by equation (2.14). The squared current density ([**J**]_{i})^{2} can be shown to attain its maximum value when the actual scalp potential is generated by a point source at the target. For instance, a point source of strength *a* at the *j*th voxel gives
2.18
and
2.19
The partial derivative of the squared current density (equation (2.19)) with respect to **K**_{i} gives
2.20
which attains zero value when **K**_{i}=**K**_{j}.

This proves that the linear tomography in equation (2.10), with weights given by equation (2.13), solves the multi-variate inverse problem and, at the same time, achieves exact localization to test point sources under ideal (no-noise) conditions.

The linear tomography defined by equation (2.10) with weights given by equation (2.13) is denoted as exact low-resolution electromagnetic tomography (eLORETA) [6,31].

### (f) Exact low-resolution electromagnetic tomography validation

eLORETA was tested under computer-controlled conditions, using a realistic head model, with 71 electrodes and 7002 cortical voxels. In this case, non-ideal conditions were used, with noise in measurements (signal to noise ratio SNR=10), which will necessarily produce non-zero localization errors. Table 1 shows a number of performance measures for several weighted linear solutions; eLORETA outperforms all other linear solutions. In addition, simulations were carried out under ideal noiseless conditions, with eLORETA attaining zero localization error.

eLORETA was validated with real human EEG recordings obtained under diverse stimulation conditions, in order to verify the correct localization of activity in sensory cortices. Figure 1 shows correct localization of visual, auditory and somato-sensory cortices.

## 3. Intracranial coherence and phase synchronization

### (a) Basic definitions

Dynamic functional connectivity between two brain regions will be quantified here as the ‘similarity’ between time-varying signals recorded at the two regions [33]. The generic term ‘similarity’ allows for a variety of measures that have been evaluated and used in the analysis of time series of electric neuronal activity. Quian-Quiroga *et al*. [34] and Dauwels *et al*. [35] provide excellent reviews.

Coherence provides a measure of linear similarity between signals in the frequency domain (e.g. [36]). Phase synchronization corresponds to a very particular form of nonlinear similarity. Both measures can be extended to time-varying Fourier transforms (or any complex valued wavelet transform).

Let *x*_{jt} and *y*_{jt} denote two stationary time series, with *j*=1,…,*N*_{R}, denoting the *j*th segment or epoch. Define the real-valued vector
3.1
Its Fourier transform at frequency *ω* is denoted as
3.2
It will be assumed that *x*_{ω} and *y*_{ω} have zero mean. The Hermitian covariance, which is proportional to the cross-spectral matrix, is
3.3
where the superscript asterisk denotes transpose and complex conjugate.

A classic expression for the complex-valued coherence is 3.4 Phase synchronization is equivalent to the coherence, but based on normalized (unit modulus) Fourier transforms, which gives it its nonlinear character. Formally, define the vector 3.5 with and denoting the normalized Fourier transforms, which are complex numbers of the unit circle. This means that , and the complex valued phase synchronization is 3.6

### (b) Problems owing to volume conduction and low spatial resolution

The moduli of these measures (i.e. |*r*_{ω}| and |*ϕ*_{ω}|) are in the range zero (no similarity) to unity (perfect similarity). When these measures are computed for invasive intracranial recordings (i.e. for time series of local electric potential differences), they validly correspond to connectivity. However, for scalp EEG signals, it is invalid to assume that these measures establish connectivity between electrode sites, since electric neuronal activity does not project radially to the scalp. These connectivity measures can be validly applied to eLORETA signals, but must be corrected for bias towards higher values owing to the low spatial resolution of the method, which makes signals appear extremely similar at zero lag.

This problem has been dealt with previously in the literature, although the solutions were aimed at correcting for the volume conduction effect in scalp signals. Nolte *et al*. [37] proposed the use of the imaginary part of the coherence (denoted as as an index of true physiological connectivity, and Stam *et al*. [38] proposed an estimator for the lagged part of phase synchronization.

### (c) The instantaneous and lagged frequency components of the total connectivity

Here, we appropriately decompose the total connectivity into instantaneous and lagged contributions in such a way that the lagged component is minimally affected by the low spatial resolution artefact, thus containing almost pure physiological information. We follow the seminal work of Geweke [39], in which total connectivity is additively expressed in terms of instantaneous and Granger-causal components 3.7 In particular, we use the observation by Parzen [40], which noted that these measures correspond to likelihood tests for different hypotheses on dependence.

Let
3.8
denote the coherence matrix, where Diag(⋅) is the diagonal matrix of the argu- ment. The total connectivity measure (including instantaneous and lagged compo- nents) is proportional to the log-likelihood statistic for testing *H*_{0}:**R**_{ZZω}=**I**, which is
3.9
where Det(⋅) denotes the determinant of the argument.

Next, we define the instantaneous (zero-lag) connectivity at discrete frequency *ω*. The direct, straightforward definition is based on the time domain zero-lag covariance of the filtered time series. For this purpose, consider the time series **Z**_{jt} (equation (3.1)), and let denote its filtered version to the single discrete frequency *ω*. The time domain zero-lag covariance is
3.10
The instantaneous connectivity corresponds to the log-likelihood statistic for testing that the correlation matrix of **A**_{ZZ} is the identity matrix. However, based on Parseval's theorem, the time-domain covariance **A**_{ZZ} is proportional to the real part of the frequency-domain Hermitian covariance (equation (3.3)). Formally, the instantaneous (zero-lag) connectivity measure is proportional to the log-likelihood statistic for testing *H*_{0}:**R**^{re}_{ZZω}=**I**, which is
3.11
where the superscript ‘re’ denotes the real part of the complex-valued matrix or number.

Finally, from equation (3.7), lagged (non-instantaneous) connectivity is
3.12
which can be expressed explicitly as
3.13
Note that these measures of connectivity ‘*F*’ take values in the range zero (no similarity) to infinity (perfect similarity). They can be transformed to the [0,1] interval, as a squared coherence or synchronization measure, by applying the transformation [41]
3.14
Note however, that these squared coherences do not satisfy the additive property in equation (3.7).

Thus, the total squared coherence is the classical definition 3.15 where the superscript ‘im’ denotes the imaginary part. The instantaneous squared coherence is related to the real part of the coherence 3.16 and the lagged squared coherence is 3.17 These measures can be applied to the normalized Fourier transforms, as defined above (equations (3.5) and (3.6)), giving the same decomposition for the phase synchronization, where the lagged component is pure physiological, and affected minimally by low spatial resolution, which affects the instantaneous component.

This methodology has been generalized to include measures of similarity between two or more multi-variate time series [42]. Furthermore, the definitions extend directly to frequency bands, by using the pooled Hermitian covariance (equation (3.3) or based on phase information from equation (3.5)) over the group of discrete frequencies of choice.

We note that this approach can be applied to single discrete frequencies. For such a case, the autoregressive model degenerates and cannot be used, since the autoregression is deterministic and not stochastic for time series that are purely sinusoidal in nature. However, when the signals have a broad-band spectrum, the direct use of autoregressive type models can be very informative (e.g. [35,43,44]). A further limitation of the methods presented here is that the actual causality cannot be resolved when the signals are filtered to single frequencies, since the concept of lag for pure sinusoidal signals is ambiguous. This means the lagged connectivity measure cannot be unambiguously decomposed into two causal components, as is the case with autoregressive modelling of more complex signals (e.g. [45]).

### (d) A comparison of non-instantaneous connectivity measures

We compare in a simple setting the new ‘lagged connectivity’ measure introduced here ( in equation (3.17)), with the measure proposed by Nolte *et al*. [37] for true brain interaction, given by the imaginary part of the coherence
3.18
Specifically, we generated data as
and
3.19
where the time series *c*_{jt}, *z*_{jt}, *d*_{jt} and *e*_{jt}, were independent and identically distri- buted uniform random variables, with *c*,*z*∼*U*[−1,+1] and *d*,*e*∼*U*[−0.1,+0.1]. The simulations consisted of using different values for *γ* in the range 0.2–4.8, and for each value, generating 500 epochs of 1 s duration each, sampled at 256 Hz. We computed the lagged connectivity (this work) and the imaginary part of the coherence [37], and show the results at 8 Hz frequency. Ideally, any change in the instantaneous component (determined by *γ*) should have little effect on the lagged or imaginary connectivities. Figure 2 shows that as the instantaneous component increases, both measures decrease, with the imaginary coherence tending to zero very quickly, while the lagged connectivity *ρ*^{2}_{Lag} remains non-zero. The relative decrease is also much higher for the imaginary coherence.

The important point to note is that when there exists a lagged connection, the imaginary part of the coherence [37] fails to detect it by tending to zero if the instantaneous component is large. This is not the case for the lagged coherence in equation (3.17), which asymptotically tends to a non-zero value, detecting the presence of a physiological lagged connection.

### (e) An example: disconnections in schizophrenia

In a previous study [46], the neuronal generators of oscillatory activity were compared between a group of patients with schizophrenia and a healthy control group. Specific patterns of cortical locations at specific EEG frequencies were observed to be different between the groups.

Using the same material [46], connectivity matrices were computed for the subjects in each group, for instantaneous and lagged connectivity measures (equations (3.11) and (3.13)). The intracranial signals were now computed with eLORETA, at 19 cortical sites, located under the electrodes of the 10/20 system. A statistical comparison of the connectivity matrices between the two groups was carried out, using non-parametric randomization techniques with correction for multiple testing.

No significant differences were found for the instantaneous connections. However, for the physiological lagged measure, a significant generalized, widespread disconnection was observed in schizophrenia, particularly notable in the theta and lower alpha bands. Figure 3 shows the significantly disconnected areas (as blue wires) in schizophrenia, for the lower alpha band (8.5–10 Hz). These results are in agreement with a number of other studies that report abnormal connectivity in schizophrenia (e.g. [47]).

## 4. Summary and outlook

We presented a new method for calculating cortical electric neuronal activity distributions from EEG measurements. This can be used for functional localization, as in classical neuroimaging; but more importantly, it provides non-invasive intracranial recordings for the assessment of dynamic functional connectivity. We also presented a method for assessing connectivity between pairs of brain regions that is minimally affected by volume conduction and low spatial resolution, thus revealing pure physiological connectivity.

Using visual, auditory and somatosensory evoked responses obtained from different laboratories, eLORETA is shown to correctly localize function in the primary and secondary sensory cortices. This evidence validates the eLORETA method in terms of functional localization. Unfortunately, experimental data do not exist that can be considered as a gold standard for testing connectivity. Nevertheless, our connectivity analysis showing that schizophrenia is characterized by a highly disconnected brain is in general agreement with the literature.

These techniques provide information on brain function, in terms of localization and interactions. While there has been much development in statistics for the analysis of functional localization, there is still much need of tools for summarizing and interpreting connectivity data. Two approaches in the literature seem very promising: (i) graph-theoretical methods [48] and (ii) independent components and singular value decomposition methods [33,49]. Our next goal consists of adapting, extending and applying these connectivity analysis techniques to high time-resolution intracranial signals.

## Footnotes

One contribution of 11 to a Theme Issue ‘The complexity of sleep’.

- This journal is © 2011 The Royal Society