Brain Electric Model

Initialize NeuroAnalyzer 
using NeuroAnalyzer
using LinearAlgebra
NeuroAnalyzer.Random.seed!(1)

The brain electric model is the physical and mathematical framework describing how neural activity generates the electric potentials measured on the scalp (EEG) or intracranially (ECoG, sEEG), and the magnetic fields measured outside the head (MEG). It is the basis of forward modelling (predicting the measured signals from known sources) and of source localization (the inverse problem).

Concept Description
Neural sources currents generated by synchronously active cortical pyramidal neurons, modelled as current dipoles
Volume conduction spread of the currents and potentials through the tissues of the head
Head model the geometry and conductivities of the head tissues (spheres, or realistic models from MRI)
Forward problem calculating the scalp potentials (or magnetic fields) for given sources and a head model
Inverse problem estimating the sources from the measured scalp potentials (or fields)

See also Introduction to EEG for the physiological origin of the EEG.

Physics

Neural sources

The synchronous postsynaptic activity of the apical dendrites of pyramidal neurons is modelled as a current dipole with:

  • location \(\mathbf{r}_0\)
  • orientation (the direction of the intracellular current, perpendicular to the cortical surface)
  • moment \(\mathbf{q}\) (current × distance, in A·m; typical evoked responses: 10–100 nAm, i.e. the synchronous activity of about \(10^4\)–\(10^5\) pyramidal neurons)

Volume conduction

The head is a volume conductor. At the frequencies of brain activity, the quasi-static approximation of Maxwell’s equations applies (capacitive and inductive effects are negligible), and the potential \(\phi\) satisfies Poisson’s equation:

\[ \nabla \cdot \left( \sigma(\mathbf{r}) \nabla \phi(\mathbf{r}) \right) = \nabla \cdot \mathbf{J}^p(\mathbf{r}) \]

where \(\sigma(\mathbf{r})\) is the conductivity of the tissue at \(\mathbf{r}\) and \(\mathbf{J}^p\) is the primary (neural) current density; the divergence of the primary current acts as the current source density. The total current is the sum of the primary current and the volume (return) currents \(-\sigma \nabla \phi\). On the scalp surface, no current flows into the air (zero normal current).

Conductivities of the head tissues (approximate values; they vary between studies and subjects):

Tissue Conductivity (\(\sigma\)) Notes
Gray matter 0.33 S/m
White matter 0.14 S/m anisotropic: higher along the fibers
CSF 1.79 S/m the most conductive tissue: shunts currents and spreads the potentials
Skull 0.004–0.02 S/m 15–80× lower than the brain and scalp (the ratio of 80 of older models is now considered too high)
Scalp 0.33–0.43 S/m

The low conductivity of the skull attenuates the scalp potentials and smears them over a large area, which limits the spatial resolution of the EEG. MEG is much less affected, because the magnetic permeability of all tissues is that of free space.

Potential of a dipole in an infinite homogeneous medium (the simplest model):

\[ \phi(\mathbf{r}) = \frac{1}{4\pi\sigma} \frac{\mathbf{q} \cdot (\mathbf{r} - \mathbf{r}_0)}{|\mathbf{r} - \mathbf{r}_0|^3} \]

The potential is positive in the direction of the dipole moment, negative in the opposite direction, and decreases with the square of the distance. In a bounded head, the insulating boundary (the scalp-air interface) increases the scalp potentials, and the layered conductivities change their distribution; these effects are described by the head models below.

Forward Problem

The forward problem is calculating the potentials (or magnetic fields) at the sensors for given sources and a head model. It has a unique solution.

  1. Head model: spherical (analytical solutions) or realistic, from the individual MRI (numerical solutions).
  2. Source model: one or a few dipoles, or a distributed model with thousands of dipoles on a grid or on the cortical surface.
  3. Lead field: the potentials at the sensors produced by unit dipoles at each source location and orientation; for \(N\) sensors and \(M\) source components it is an \(N \times M\) matrix \(\mathbf{L}\), and the measured data are a linear combination of its columns:

\[ \mathbf{m} = \mathbf{L} \mathbf{j} + \mathbf{n} \]

where \(\mathbf{m}\) is the vector of measured potentials, \(\mathbf{j}\) the source strengths and \(\mathbf{n}\) the noise. The lead field depends on the reference: average-referenced data require an average-referenced lead field.

Example: dipoles in a spherical head

Scalp potentials of a dipole in a homogeneous conducting sphere (radius 9 cm, conductivity 0.33 S/m), calculated for the 47 electrodes of the 10-10 system. NeuroAnalyzer uses this model (the Brody series of the surface potential) in dipole_fit(); its internal function _leadfield() returns the lead field of a dipole at a given position.

The dipole is placed 7 cm from the center of the head, under the electrode C3, with a moment of 20 nAm, either radial (perpendicular to the scalp) or tangential (parallel to the scalp, pointing anteriorly):

locs = import_locs("files/locs/standard-10-10-cap47.ced")
xyz = Matrix(locs[:, [:loc_x, :loc_y, :loc_z]])
u = xyz ./ sqrt.(sum(abs2, xyz; dims = 2))        # unit vectors to electrodes
R = 0.09                                          # head radius [m]
σ = 0.33                                          # conductivity [S/m]

c3 = u[findfirst(==("C3"), locs.label), :]
r0 = 0.07 .* c3                                   # dipole position [m]
q_rad = 20e-9 .* c3                               # radial moment [A·m]
q_tan = 20e-9 .* normalize(cross(c3, [0, 0, 1]))  # tangential moment

L = NeuroAnalyzer._leadfield(r0, u, R, σ)         # lead field [V / A·m]
v_rad = L * q_rad .* 1e6                          # scalp potentials [μV]
v_tan = L * q_tan .* 1e6
println("radial dipole:     ", round.(extrema(v_rad); digits = 1), " μV")
println("tangential dipole: ", round.(extrema(v_tan); digits = 1), " μV")
radial dipole:     (-0.7, 26.8) μV
tangential dipole: (-3.7, 3.7) μV

The radial dipole produces a single focal maximum directly above the source:

plot_topo(v_rad;
          locs = locs,
          cb_title = "[normalized]",
          title = "radial dipole")

The tangential dipole produces a bipolar pattern: a positive and a negative pole on both sides of the source, with zero potential directly above it, and lower amplitudes:

plot_topo(v_tan;
          locs = locs,
          cb_title = "[normalized]",
          title = "tangential dipole")

The maximum of the EEG is therefore not always above the source: for tangential (e.g. sulcal) sources, the source lies between the positive and the negative pole.

Inverse Problem

The inverse problem is estimating the sources (locations, orientations and strengths) from the measured data. It is ill-posed:

  1. Non-uniqueness: infinitely many source configurations produce the same scalp potentials (e.g. sources whose fields cancel, “silent sources”); additional assumptions (constraints) are needed to obtain a unique solution.
  2. Instability: small changes of the data (noise) can produce large changes of the solution, so the estimate must be regularized.
  3. Model errors: errors of the head model (geometry, conductivities, electrode positions) cause localization errors.

Dipole fitting

Assumption: the data are generated by a small number of focal sources (equivalent current dipoles).

  1. Initialization: the number of dipoles and their starting positions (e.g. from a grid search).
  2. Optimization: the dipole positions are found by non-linear optimization (e.g. Nelder-Mead, Levenberg-Marquardt) minimizing the residual variance; for a given position, the moments are found by linear least squares.
  3. Validation: the goodness of fit (the explained variance), the plausibility of the location and the confidence volume.

Suited for focal activity (early sensory responses, epileptic spikes); the number of dipoles must be known, and the fit can end in local minima.

dipole_fit() fits a single dipole in a homogeneous sphere to the simulated data (in the same model, the solution is exact):

true_pos = r0 ./ R                                 # in the units of locs
for (name, v) in (("radial", v_rad), ("tangential", v_tan))
    d = dipole_fit(reshape(v, :, 1), locs; loc_scale = R)
    err = 1000 * R * norm(vec(d.pos) .- true_pos)
    println(rpad(name, 11), "  error: ", round(err; digits = 1), " mm",
            "  goodness of fit: ", round(d.gof[1]; digits = 3))
end
radial       error: 0.0 mm  goodness of fit: 1.0
tangential   error: 0.0 mm  goodness of fit: 1.0

With noise (10% of the maximum amplitude of each map), the localization error increases; how much depends on the source (its depth and orientation), the number and positions of the electrodes and the noise level:

for (name, v) in (("radial", v_rad), ("tangential", v_tan))
    v_noisy = v .+ 0.1 * maximum(abs.(v)) .* randn(length(v))
    d = dipole_fit(reshape(v_noisy, :, 1), locs; loc_scale = R)
    err = 1000 * R * norm(vec(d.pos) .- true_pos)
    println(rpad(name, 11), "  error: ", round(err; digits = 1), " mm",
            "  goodness of fit: ", round(d.gof[1]; digits = 3))
end
radial       error: 23.8 mm  goodness of fit: 0.76
tangential   error: 4.6 mm  goodness of fit: 0.941

Distributed source models

Assumption: the activity is distributed over many sources (e.g. thousands of dipoles on the cortical surface); the number of sources is much larger than the number of sensors, so constraints are added:

  1. Minimum-norm estimate (MNE): the solution with the smallest overall source strength (L2 norm) that explains the data:

\[ \hat{\mathbf{j}} = \arg\min_{\mathbf{j}} \left( \|\mathbf{m} - \mathbf{L}\mathbf{j}\|^2 + \lambda \|\mathbf{j}\|^2 \right) = \mathbf{L}^T (\mathbf{L}\mathbf{L}^T + \lambda \mathbf{I})^{-1} \mathbf{m} \]

where \(\lambda\) is the regularization parameter, balancing the fit to the data and the stability of the solution (it depends on the noise level; the noise covariance is used instead of \(\mathbf{I}\) in practice). The two equivalent forms \((\mathbf{L}^T\mathbf{L} + \lambda \mathbf{I})^{-1} \mathbf{L}^T \mathbf{m}\) and \(\mathbf{L}^T (\mathbf{L}\mathbf{L}^T + \lambda \mathbf{I})^{-1} \mathbf{m}\) give the same result; the second inverts only an \(N \times N\) matrix (sensors). MNE favours superficial sources (they explain the data with smaller strengths). 2. Weighted MNE (wMNE): depth weighting reduces the bias towards superficial sources. 3. LORETA (low-resolution electromagnetic tomography): the smoothest solution (minimum norm of the spatial Laplacian of the sources); low spatial resolution, but less depth bias. 4. Noise-normalized estimates: dSPM and sLORETA divide the MNE by its estimated noise (or resolution) at each location, which gives statistical maps and better localization of the maximum (sLORETA: zero localization error for single sources without noise); eLORETA uses weights giving exact localization of single sources. 5. Bayesian methods (e.g. multiple sparse priors, Champagne) estimate the sources and the hyperparameters (e.g. source variances) from the data.

Distributed models need no assumption about the number of sources and capture extended activity, but the estimates are blurred and depend on the regularization.

Beamforming

Beamformers (e.g. LCMV in the time domain, DICS in the frequency domain) construct a spatial filter for each location that passes the activity from that location and suppresses the activity from elsewhere (minimizing the output variance under the constraint of unit gain at the location); the brain is scanned location by location.

They give a good spatial resolution for well-separated sources and are robust to noise, but they assume that the sources are uncorrelated: correlated sources (e.g. both auditory cortices) are suppressed. They need a good estimate of the data covariance (enough data).

Machine learning

Deep neural networks (e.g. convolutional networks) can be trained on simulated data (from forward models) to map scalp topographies to source distributions; they depend on the realism of the training simulations.

Head Models

Spherical models:

  1. Homogeneous sphere: one compartment; analytical solution (used above).
  2. Three-shell model: concentric spheres for the brain, skull and scalp with different conductivities; analytical solution (series of Legendre polynomials); the standard simple model for EEG.
  3. Four-shell model: adds the CSF.

Spherical models fit the upper part of the head reasonably well but not the temporal and frontal regions; for MEG, a single sphere is often sufficient (the conductivities do not affect the field outside a spherically symmetric conductor).

Realistic models (from the individual MRI, or a template):

  1. Boundary element method (BEM): the head is described by closed surfaces (e.g. brain, inner and outer skull, scalp) separating compartments of constant, isotropic conductivity; efficient, the standard realistic model for EEG and MEG.
  2. Finite element method (FEM): the head is divided into a volume mesh (e.g. tetrahedra), each element with its own conductivity, which allows any geometry, holes in the skull and anisotropic conductivity (white matter, skull layers); most accurate, highest computational cost.
  3. Finite difference method (FDM): a regular grid (e.g. the MRI voxels); simple, but geometry is represented only approximately.

Applications

Application Description
Source localization epileptic foci, eloquent cortex, generators of ERP components
Source-space analysis power and connectivity of brain regions instead of channels (less affected by volume conduction)
Simulations generating realistic test data with known sources, e.g. to evaluate analysis methods
Brain-computer interfaces source estimates as features for decoding (e.g. motor imagery)
Brain stimulation modelling the electric fields of transcranial electric (tDCS, tACS) and magnetic stimulation (TMS)

Forward vs. Inverse Problem

Aspect Forward problem Inverse problem
Definition calculate the scalp EEG/MEG from known sources estimate the sources from the EEG/MEG
Uniqueness unique solution infinitely many solutions; constraints needed
Input sources and head model measured data, head model (lead field)
Output potentials or fields at the sensors source locations, orientations and strengths
Methods analytical (spheres) and numerical (BEM, FEM, FDM) solutions dipole fitting, minimum-norm methods, beamformers, Bayesian and machine-learning methods
Main difficulties accuracy of geometry and conductivities, computational cost non-uniqueness, noise, model errors, choice of constraints and regularization