Misusing LORETA

 

This software computes LORETA from scalp electric potential differences (time domain EEG/ERP) or from EEG cross-spectra (frequency domain). One particular very incorrect usage is to cheat LORETA with the input. Examples of misuse:

1. Inputting scalp electric potential spectral powers will not output LORETA (current density) spectral powers.

2. Inputting scalp electric potential square roots of spectral powers will not output LORETA (current density) square roots of spectral powers.

3. Inputting scalp z-transformed-maps will not output LORETA (current density) z-transformed-values.

The three previous invalid inputs to LORETA violate the mathematics and the physics underlying all computations. Furthermore, they violate any correct usage of statistical analysis.

Some more technical details can be found in:

1. For time domain computations: Pascual-Marqui RD: Review of methods for solving the EEG inverse problem. International Journal of Bioelectromagnetism 1: 75-86, 1999.

2. For frequency domain computations: Frei E, Gamma A, Pascual-Marqui R, Lehmann D, Hell D, Vollenweider FX: Localization of MDMA-induced brain activity in healthy volunteers using low resolution brain electromagnetic tomography (LORETA). Human Brain Mapping 14: 152-165, 2001. See text and equations on pages 154-155 therein.

 

Proof is given in what follows:

 

Introduction

 

The forward problem of EEG is formulated as:

                                                                                                      Eq. 1

 

In Eq. 1,  is a vector containing scalp electric potentials measured at  cephalic electrodes, with respect to a common, arbitrary reference electrode located anywhere on the body.

 

Note: The “discrete” formulation in Eq. 1 is derived from the integral equation formulation. For example, the set of measured scalp electrode potential differences are collected into a vector.

 

The primary (impressed) current density  is defined as:

                                                                                     Eq. 2

where  for . At the lth voxel,  contains the three unknown dipole moments.

 

The cases considered here correspond to .

 

The superscript “T” denotes transpose.

 

The lead field  has the following structure:

                                                                          Eq. 3

with , for , and for . Note that , where  is the scalp electric potential at the ith electrode, due to a unit strength X-oriented dipole at the lth voxel;  is the scalp electric potential at the ith electrode, due to a unit strength Y- oriented dipole at the lth voxel; and  is the scalp electric potential at the ith electrode, due to a unit strength Z- oriented dipole at the lth voxel.

 

In Eq. 1, c is an arbitrary constant which embodies the fact that the electric potential is determined up to an arbitrary constant; and  is a vector of ones. The parameter c allows the use of any reference for the lead field and the measurements.

 

The functional of interest here is:

                                                                            Eq. 4

where  is a regularization parameter. This functional is to be minimized with respect to J and c, for given K, F, and a. In Eq. 4, the matrix B implements a discrete 3D Laplacian (for smoothness), and W normalizes the columns of the lead field. This is explained in detail in {Pascual-Marqui RD. Review of methods for solving the EEG inverse problem. International Journal of Bioelectromagnetism 1999, 1: 75-86.}.

 

The explicit solution to this minimization problem is:

                                                                                                            Eq. 5

where:

                                         Eq. 6

                                                                                                Eq. 7

with  denoting the centering matrix;  the identity matrix; and  is a vector of ones.

 

For any matrix ,  denotes its Moore-Penrose pseudoinverse.

 

The centering matrix H in Eq. 7 is the average reference operator.

 

Note: The general form of Eq. 5 is valid for any discrete, distributed, linear inverse. The equation for T (e.g. Eq. 6) will change with the particular pseudo-inverse.

 

The time domain inverse solution (LORETA)

 

This corresponds to Eq. 5, which will now be written explicitly as:

                                                                                                           Eq. 8

where the subscript t signifies an instantaneous measurement at time instant t.

 

For a time series, with , Eq. 8 is then generalized to:

                                                                                                        Eq. 9

where , and .

 

Note that the columns of  and  correspond to all discrete time instants.

 

The frequency domain inverse solution (LORETA)

 

The discrete Fourier transform operator is denoted as the matrix . This holds, regardless of the “fast” algorithms used in practice. Fast algorithms implement multiplication with F in a fast way, nothing more. This operator gives the DFT:

                                                                                              Eq. 10

or equivalently:

                                                                                                       Eq. 11

with:

                                                                                     Eq. 12

                                                                                        Eq. 13

 

Note that the columns of  and  correspond to all discrete frequencies.

 

For a single discrete frequency w, with , Eq. 11 is:

                                                                                                         Eq. 14

with , and .

 

Note: Although there are  discrete frequencies, there is redundancy in the fact that about half the discrete frequencies are conjugates of the other half. The exact relations can be found in classical FFT books. This is of no interest here.

 

VERY IMPORTANT Note: Up to here, all transforms are linear. However, we are not yet dealing with spectra!

 

The computation of spectral density for electric neuronal activity (LORETA spectral density)

 

Let , for , denote the DFT at frequency w, for the ith EEG epoch. For example, an EEG epoch may contain 256 time frames (discrete time instants), sampled at 128 Hz, corresponding to 2 seconds worth of EEG per epoch. We typically use about N=20 such epochs.

 

From Eq. 14, the cross-spectral matrices are constructed as follows.

                                                                                                  Eq. 15

                                                            Eq. 16

where:

                                                                                             Eq. 17

and where the superscript “*” denotes transpose and complex conjugate. This is equivalent to:

                                                                                                   Eq. 18

where:

                                                                                        Eq. 19

                                                                                        Eq. 20

 

In Eqs. 19 and 20,  and  denote the cross-spectral matrices for the current density and the scalp electric potential differences, respectively. These are Hermitian matrices.

 

The concept of cross-spectral matrices can be found in {Brillinger DR (1981): Time series: data analysis and theory. New York: McGraw-Hill}.

 

Note: The spectra at each voxel, i.e. the spectra of the “electric neuronal activity” time series estimated via LORETA at each voxel in cortical grey matter, corresponds to the diagonal elements of the matrix , obtained via Eq. 18. This means that at frequency w, we are concerned with:

                                                                               Eq. 21

where the “diag” operator takes a Hermitian matrix and returns a real diagonal matrix with elements corresponding to the diagonal elements of the original Hermitian matrix.

 

VERY IMPORTANT Note: These transforms are not linear. Spectra (and cross-spectra) are second order moment statistics!

 

VERY IMPORTANT Note: There is no way to obtain the LORETA spectra with only the spectra of the electric potentials.

 

The following inequality holds:

                                                                           Eq. 22

This inequality holds in general, since the “diag” operator in Eq. 21 does not commute with matrix multiplication. If you don’t believe me, one example suffices to prove that the “equality” does not hold. This was confirmed empirically with randomly generated matrices.

 

The following inequality holds:

                                                                   Eq. 23

Based on the inequality in Eq. 22, this one is even worse, i.e., it’s much more of an “inequality”.

 

Let:

                                                                                                  Eq. 24

where “vdiag” is an operator such that, given a square Hermitian matrix M, V is a vector containing the diagonal elements of M.

 

Using the notation in Eq. 24, note that the correct LORETA spectra can be written as:

                                                                            Eq. 25

 

The following inequality holds:

                                               Eq. 26

This inequality holds in general. If you don’t believe me, one example suffices to prove that the “equality” does not hold. This was confirmed empirically with randomly generated matrices.