Introduction to source reconstruction

Problem outline

One of the holy grails of EEG analysis would be to infer where in the brain the signals are generated that constitute the EEG. It is well established that the EEG does not reflect the activity of single neurons, but rather the combined electrical activity of spatially extended patches of simultaneously active cortical tissue with a minimal area in the order of several square centimeters. To go from these patches of simultaneous active cortex to EEG signals, we will use the terminology and mathematical notation described in (Grech et al., 2008)

Assumptions

In this module, we will make a few key assumptions to develop elegant mathematical models for the forward and inverse problems in EEG. These assumptions are:

Forward problem

To transform the electrical signals recorded on the EEG into the source space (i.e., to estimate where in the brain they originated), we first need to understand how brain activity generates signals measurable at the scalp. This is known as the forward problem in EEG. The first step is to define how neural sources are represented. The best current approximation models each source as an electric dipole — a pair of equal and opposite charges separated by a very small distance. Because of this structure, a dipole has both a magnitude (strength) and a direction, and can therefore be represented as a vector.

This dipole model is based on the properties of pyramidal neurons, which are long, oriented cells in the cerebral cortex. When many of these neurons are active and aligned, their combined electrical fields produce measurable EEG signals. Thus, the dipole serves as a simplified but biologically grounded model of how neuronal populations contribute to scalp potentials.

In the diagram below, adapted from (Grech et al., 2008), a simplification of the forward problem is shown.
Schematic three-layer head model (scalp, skull, brain) with a single dipole at position r_dip_i and one measurement electrode at r. Dipole orientation decomposed into x and z axes.
Diagram adapted from Grech et al., 2008. A three-layer approximation of the head (scalp, skull, brain) is shown with a single dipole located at position \( \mathbf{r}_{\text{dip}_i} \) from the origin. In this two-dimensional plane, any dipole orientation at that location can be decomposed into two orthogonal axes, \( \mathbf{e}_{\text{x}} \) and \( \mathbf{e}_{\text{z}} \). For simplicity, one measurement electrode \( \mathbf{m}_{\text{r}} \) is drawn at radius \( \mathbf{r} \).


The source space can be represented in one of two main ways. In a volumetric model, dipoles oriented along the x, y, and z axes are placed at discrete points throughout the whole brain volume. Alternatively, sources can be restricted to the cortical surface, where dipoles are oriented perpendicular (normal) to the cortex. This approach reduces the number of parameters to estimate and reflects the organization of cortical pyramidal neurons.

By solving the Poisson equation—using the tissues’ electrical conductivities, geometry, source locations, and EEG electrode positions as inputs—we can estimate how each dipole influences each electrode. Two examples of dipoles and their influence on the EEG electrodes are given below

Normal dipole
The electric field and influence on EEG electrodes for a normal oriented dipole
Tangential dipole
The electric field and influence on EEG electrodes for a tangential oriented dipole


According to the principle of superposition, the EEG signal measured at electrode \( m \) at time \( t \) can be expressed as a linear combination of dipole contributions, each weighted by a gain factor \( g \). All gain factors can be collected into an \( M \times N \) matrix, where \( M \) is the number of electrodes and \( N \) is the number of sources. Typically, \( M \ll N \), making the system underdetermined; in practice, \( M \) ranges from roughly 20 to 256 electrodes, whereas \( N \) is on the order of several thousand potential sources. $$ \mathbf{M} = \begin{bmatrix} m(\mathbf{r}_1) \\ \vdots \\ m(\mathbf{r}_N) \end{bmatrix} = \begin{bmatrix} g(\mathbf{r}_1, \mathbf{r}_{\text{dip}_1}) & \cdots & g(\mathbf{r}_1, \mathbf{r}_{\text{dip}_p}) \\ \vdots & \ddots & \vdots \\ g(\mathbf{r}_N, \mathbf{r}_{\text{dip}_1}) & \cdots & g(\mathbf{r}_N, \mathbf{r}_{\text{dip}_p}) \end{bmatrix} \begin{bmatrix} d_1 \mathbf{e}_1 \\ \vdots \\ d_p \mathbf{e}_p \end{bmatrix} $$

Or in compact notation, with added gaussian noise from the amplifying system \( n \sim \mathcal{N}(0,\sigma^{2}) \) $$M = GD+n$$

The minimum norm inverse solution

Since the forward problem is non-square due to the vast amount of sources compared to the limited amount of measuring electrods, the matrix \( G \) cannot be simply inverted and applied to the EEG signal \( M \) to find the sources \( D \) . Consequently, there is an infinite amount of solutions to this problem. Choosing an inverse method is thus choosing a solution with certain properties. The most simple solution is choosing a solution where the sources have minimal power with regards to a certain norm. This family of estimates is usually referred to as the Minimum norm estimates (MNE) To find the minimum norm solution with regularization, we solve the following minimization problem:

\[ \min_D \| M - G D \|_2^2 + \lambda \| D \|_2^2 \]

Here, \( \lambda > 0 \) is a regularization parameter that balances the fitting of the residual and the norm of the sources.

Taking the derivative with respect to \( D \) and set it to zero:

\[ -2 G^\top (M - G D) + 2 \lambda D = 0 \]

Rearranging:

\[ G^\top G D + \lambda D = G^\top M \quad \Rightarrow \quad (G^\top G + \lambda I) D = G^\top M \]

Solve for \( D \):

\[ \boxed{ D_\text{min-norm} = (G^\top G + \lambda I)^{-1} G^\top M = G^\top (G G^\top + \lambda I)^{-1}M } \]

These are called the regularized minimum norm solutions, where the second form is obtained by applying the Woodbury matrix identity. Since \(G \) is an \( M \times N \) matrix, the first form requires the inversion of an \( N \times N \) matrix and the second form only a \( M \times M \) matrix, obviously, for computational purposes, the second form is preferred. Inherently, this solution is biased towards shallower sources, since sources closer to the measurement electrodes require less amplitude to give a similar signal with respect to deeper sources. One way to solve this is by changing the minimization problem to include a weighting matrix \( W \) which is diagonal with entries equal to the norms of the lead fields. Since sources with smaller lead field norms are typically deeper, this weighting effectively compensates for depth bias.

\[ \min_D \| M - G D \|_2^2 + \lambda \| WD \|_2^2 \]

This leads to the solution for the weighted minimum norm where again, the Woodbury matrix identity has been used to find the computationally more efficient form.

\[ \boxed{ D_\text{weighted-min-norm} = (G^\top G + \lambda W^\top W)^{-1} G^\top M = W^{-1}G^\top (G (W^\top W)^{-1}G^\top + \lambda I)^{-1}M } \]


Parametric methods

The minimum norm methods belong to the family of distributed source models, which are built on the assumption that the brain activity can be represented by many small sources. These models estimate the contribution of all potential sources simultaneously. In contrast, parametric methods attempt to explain the measured voltages by a limited amount of equivalent current dipoles.

The most intuitive of these methods is the dipole fitting method, which identifies the dipole (or dipoles) located that best fits the observed scalp potential topography. This means fitting six parameters per dipole, 3 for the spatial location and three for its moment (orientation and strength). This is a non-linear minimization of the cost function defined by :

\[ || M - G({r_j,r_{dip_i}})D || \]


Another approach is to focus on the sources rather than the measurements. A beamformer applies spatial filtering to estimate the activity of a source at a specific location. It computes a weighted linear combination of sensor signals that maximizes sensitivity to the target source while minimizing sensitivity to other brain regions. The weights are derived from the lead field and the covariance of the measured data, making the method a data-driven approach to source localization.

Other types of inverse solutions

he solutions described in the previous section are among the most fundamental, because their mathematical foundations are elegant and relatively straightforward to describe. Countless other algorithms have been developed, each constraining the solution to the problem in different ways. These include methods that provide exact localization of the source in a noiseless situation (sLORETA, eLORETA) and methods that introduce sparsity in the source space (FOCUSS, MxNE).