EEG Source Localization with PhysioEEG

This vignette demonstrates EEG source localization methods available in PhysioEEG. Source localization (also called source imaging) addresses the inverse problem: estimating brain source activity from scalp-recorded potentials. The EEG inverse problem is ill-posed (infinitely many source configurations can produce the same scalp pattern), so additional constraints are needed to obtain a unique solution.

PhysioEEG provides:

Forward modeling with spherical head and simplified BEM models
Distributed source imaging via eLORETA, sLORETA, and MNE
Beamforming via LCMV and DICS spatial filters
Source-space frequency analysis for band-specific source maps
Visualization of source estimates on a 2D brain projection

Setup

library(PhysioEEG)

# Create 19-channel EEG data (10-20 system)
pe <- make_eeg(n_time = 2500, n_channels = 19, sr = 500)

# Assign standard electrode positions (required for source localization)
pe <- eegMontage(pe, system = "10-20")

# Verify positions are set
head(SummarizedExperiment::colData(pe))

For best results, preprocess the data before source localization:

# Bandpass filter to remove drift and high-frequency noise
pe <- eegFilter(pe, lowcut = 1, highcut = 40)

# Re-reference to average (recommended for source imaging)
pe <- eegRereference(pe, ref_type = "average", assay_name = "filtered")

Forward Model

The forward model describes how neural source activity maps to scalp electrode potentials. It is represented by a leadfield matrix of dimensions (n_electrodes x n_sources * 3), where each source has three orientation columns (x, y, z dipole directions).

Spherical Head Model

The simplest approach models the head as an infinite homogeneous conducting medium with conductivity 0.33 S/m:

# Build forward model with 100 source locations
fm_sphere <- eegForwardModel(pe, method = "spherical", n_sources = 100,
                              assay_name = "rereferenced")

# Inspect the leadfield
dim(fm_sphere$leadfield)
# [1]  19 300  (19 electrodes x 100 sources x 3 orientations)

# Source positions are distributed inside a sphere of radius 0.8
head(fm_sphere$source_positions)

# Electrode positions
print(fm_sphere$electrode_positions)

Simplified BEM Model

The Berg & Scherg (1994) 3-shell approximation accounts for the different conductivities of brain, skull, and scalp using virtual dipoles at scaled eccentricities:

fm_bem <- eegForwardModel(pe, method = "bem_simplified", n_sources = 100,
                           assay_name = "rereferenced")
dim(fm_bem$leadfield)
# Same dimensions, but different values due to 3-shell correction

Leadfield Interpretation

The leadfield encodes the sensitivity pattern of each electrode to each source:

# For source 1, orientation x: how much does each electrode "see"?
source1_x <- fm_sphere$leadfield[, 1]
names(source1_x) <- fm_sphere$electrode_positions$label
print(round(source1_x, 4))

# Electrodes closest to the source have the largest values

Source Estimation

Source estimation uses the forward model and data to reconstruct brain activity. PhysioEEG supports three distributed source imaging methods:

MNE (Minimum Norm Estimate)

MNE minimizes the L2 norm of the source distribution. It provides the minimum-energy solution but tends to favor superficial sources:

pe_mne <- eegSourceEstimate(pe, fm_sphere, method = "mne",
                             lambda = 0.05,
                             assay_name = "rereferenced")

# Source estimates are stored as a time x (n_sources*3) matrix
src_data <- SummarizedExperiment::assay(pe_mne, "source")
dim(src_data)
# [1] 2500  300  (2500 time points x 100 sources x 3 orientations)

sLORETA (Standardized Low-Resolution Brain Electromagnetic Tomography)

sLORETA normalizes MNE by the resolution matrix diagonal, providing zero-localization error for point sources:

pe_sloreta <- eegSourceEstimate(pe, fm_sphere, method = "sloreta",
                                 lambda = 0.05,
                                 assay_name = "rereferenced")

eLORETA (Exact Low-Resolution Electromagnetic Tomography)

eLORETA uses iterative weight estimation to achieve exact localization for point sources, even with noise:

pe_eloreta <- eegSourceEstimate(pe, fm_sphere, method = "eloreta",
                                 lambda = 0.05,
                                 assay_name = "rereferenced")

Comparing Methods

The three methods differ in their spatial resolution, depth bias, and noise sensitivity:

# Compute source power (RMS across orientations) for each method
.rms_power <- function(pe_src, n_sources) {
  src <- SummarizedExperiment::assay(pe_src, "source")
  power <- numeric(n_sources)
  for (j in seq_len(n_sources)) {
    cols <- ((j - 1) * 3 + 1):(j * 3)
    power[j] <- mean(rowSums(src[, cols]^2))
  }
  power
}

mne_power <- .rms_power(pe_mne, 100)
sloreta_power <- .rms_power(pe_sloreta, 100)
eloreta_power <- .rms_power(pe_eloreta, 100)

# Find peak source for each method
cat("MNE peak source:", which.max(mne_power), "\n")
cat("sLORETA peak source:", which.max(sloreta_power), "\n")
cat("eLORETA peak source:", which.max(eloreta_power), "\n")

Beamformer

Beamformers are spatial filters that estimate the power at each source location by constructing a filter that passes activity from the target location while suppressing contributions from elsewhere.

LCMV Beamformer

The Linearly Constrained Minimum Variance (LCMV) beamformer operates in the time domain using the data covariance matrix:

pe_lcmv <- eegBeamformer(pe, fm_sphere, method = "lcmv",
                          assay_name = "rereferenced")

# Source power is stored as a n_sources x 1 matrix
bf_power <- SummarizedExperiment::assay(pe_lcmv, "beamformer")
dim(bf_power)
# [1] 100   1

# Find the source with maximum power
peak_source <- which.max(bf_power[, 1])
cat("LCMV peak source:", peak_source, "\n")
cat("Peak power:", bf_power[peak_source, 1], "\n")

DICS Beamformer

Dynamic Imaging of Coherent Sources (DICS) operates in the frequency domain, using the cross-spectral density matrix in a specified frequency band:

# Alpha-band (8-13 Hz) source localization
pe_dics <- eegBeamformer(pe, fm_sphere, method = "dics",
                          freq_range = c(8, 13),
                          assay_name = "rereferenced")

dics_power <- SummarizedExperiment::assay(pe_dics, "beamformer")
peak_dics <- which.max(dics_power[, 1])
cat("DICS alpha peak source:", peak_dics, "\n")

# Compare with beta-band source
pe_dics_beta <- eegBeamformer(pe, fm_sphere, method = "dics",
                               freq_range = c(13, 30),
                               assay_name = "rereferenced")
dics_beta <- SummarizedExperiment::assay(pe_dics_beta, "beamformer")
cat("DICS beta peak source:", which.max(dics_beta[, 1]), "\n")

Source Power

eegSourcePower() computes frequency-band source maps from previously estimated source time series. It requires prior source estimation via eegSourceEstimate():

# Compute band power at each source location
bp <- eegSourcePower(pe_mne, source_assay = "source")
head(bp)
# source_id, band, power

# Reshape to wide format for comparison
bp_wide <- reshape(bp, direction = "wide", idvar = "source_id",
                   timevar = "band", v.names = "power")
head(bp_wide)

Custom Frequency Bands

# Narrow-band source analysis
custom_bands <- list(
  mu     = c(8, 12),
  beta   = c(15, 25),
  gamma1 = c(30, 45)
)

bp_custom <- eegSourcePower(pe_mne, bands = custom_bands,
                             source_assay = "source")
head(bp_custom)

Visualization

PhysioEEG provides eegPlotSource() for visualizing source localization results on a 2D brain projection.

Scatter Plot

Source locations are shown as points, sized and colored by amplitude. Only sources above a percentile threshold are displayed:

# Create source data for visualization
src_positions <- fm_sphere$source_positions
src_amplitudes <- bf_power[, 1]

src_df <- data.frame(
  x = src_positions$x,
  y = src_positions$y,
  amplitude = src_amplitudes
)

# Plot with 80th percentile threshold (show top 20% sources)
eegPlotSource(pe, source_data = src_df, method = "scatter",
              threshold_pct = 80)

Flat Map

The flat map mode uses interpolation to create a continuous image of source activity:

eegPlotSource(pe, source_data = src_df, method = "flatmap",
              threshold_pct = 70)

Comparing Methods Visually

# LCMV result
lcmv_df <- data.frame(
  x = src_positions$x,
  y = src_positions$y,
  amplitude = bf_power[, 1]
)
p_lcmv <- eegPlotSource(pe, source_data = lcmv_df, method = "scatter",
                          threshold_pct = 80)

# DICS alpha result
dics_df <- data.frame(
  x = src_positions$x,
  y = src_positions$y,
  amplitude = dics_power[, 1]
)
p_dics <- eegPlotSource(pe, source_data = dics_df, method = "scatter",
                          threshold_pct = 80)

# Display side by side (requires patchwork or gridExtra)
print(p_lcmv)
print(p_dics)

Method Comparison

Distributed Source Imaging vs. Beamforming

Feature	Distributed (MNE/LORETA)	Beamformer (LCMV/DICS)
Output	Full time course per source	Power per source
Multiple sources	Well-suited	Can suppress correlated sources
Depth bias	MNE biased to surface	Less depth bias
Regularization	Lambda parameter	Data covariance
Frequency specificity	Post-hoc via eegSourcePower	DICS: built-in

Choosing a Method

MNE: Good starting point, interpretable, but surface-biased.
sLORETA: Best for localization accuracy of focal sources; zero localization error for point sources in noiseless data.
eLORETA: Exact localization even with noise, but more computationally expensive due to iterative estimation.
LCMV: Effective for broadband source power estimation; robust to noise.
DICS: Best for frequency-specific source analysis (e.g., alpha generators, beta oscillations in motor cortex).

Best Practices

Electrode Density

Source localization accuracy improves with more electrodes. The 19-channel 10-20 system provides limited spatial resolution. For research applications, 64 or more channels are recommended:

# For demonstration with 19 channels
pe_19 <- make_eeg(n_time = 2500, n_channels = 19, sr = 500)
pe_19 <- eegMontage(pe_19, system = "10-20")
fm_19 <- eegForwardModel(pe_19, method = "spherical", n_sources = 50)
# Limited spatial resolution with 19 channels

SNR Considerations

Higher signal-to-noise ratio (SNR) improves source estimation. Strategies to improve SNR include:

Averaging across epochs (ERP paradigms)
Narrowband filtering before source estimation
Artifact rejection and ICA cleaning
Increasing regularization (lambda) for noisy data

# Higher regularization for noisy data
pe_noisy <- eegSourceEstimate(pe, fm_sphere, method = "sloreta",
                               lambda = 0.2,
                               assay_name = "rereferenced")

# Lower regularization for clean data
pe_clean <- eegSourceEstimate(pe, fm_sphere, method = "sloreta",
                               lambda = 0.01,
                               assay_name = "rereferenced")

Regularization

The regularization parameter lambda controls the trade-off between data fit and solution smoothness:

Small lambda (e.g., 0.01): Better data fit, but sensitive to noise and may produce spurious peaks.
Large lambda (e.g., 0.5): Smoother solution, more robust to noise, but reduced spatial resolution.
Typical range: 0.01 to 0.1 for well-preprocessed data.

# Sweep regularization values
lambdas <- c(0.001, 0.01, 0.05, 0.1, 0.5)
for (lam in lambdas) {
  pe_src <- eegSourceEstimate(pe, fm_sphere, method = "mne",
                               lambda = lam,
                               assay_name = "rereferenced")
  src <- SummarizedExperiment::assay(pe_src, "source")
  max_power <- max(rowSums(src^2))
  cat(sprintf("lambda = %.3f -> max source power = %.2e\n", lam, max_power))
}

Preprocessing Recommendations

Filter: Bandpass 1-40 Hz (or narrower for specific questions)
Re-reference: Average reference is standard for source imaging
Artifact rejection: Remove or correct eye blinks, muscle, and cardiac artifacts
Electrode positions: Use digitized positions when available; standard montage positions are an approximation
Forward model: Use BEM for more realistic head modeling when individual MRIs are available (not implemented in this package; consider MNE-Python for advanced forward modeling)