Author: Julia Jurkowska

Source localization and signal reconstruction - case study for oddball data#

Introduction#

In this tutorial, we will learn how to localize sources from EEG data and reconstruct signals at those sources using MVPURE_py, an extension to MNE-Python. Source localization allows us to move beyond sensor-level analysis to estimate where in the brain the measured activity originates. Once sources are identified, we can reconstruct time series from vertices of interest for further analysis.

We will cover the following steps:

  1. Reading all necessary data for the sample_subject. You can download this dataset here.

  2. Computing data and noise covariance (R and N, respectively).

  3. Analysis of \(RN^{-1}\) eigenvalues to guide the number of sources to localize and select an appropriate optimization parameter.

  4. Localizing the specified number of sources.

  5. Reconstructing source signals for vertices of interest and plotting the results.

  6. Comparing the results with standard LCMV beamformer and sLORETA methods.

All steps will be repeated for two time frames: “sensory” (50-200 ms after stimuli) and “cognitive” (350-600 ms after stimuli).

By the end of this tutorial, you will understand the basic workflow of source localization and signal reconstruction using the MVPURE-py package.

[1]:

import mne
import os
import numpy as np

mne.viz.set_3d_backend('pyvistaqt')

from mvpure_py import localizer, viz, utils

Using pyvistaqt 3d backend.

We will use data the sample_subject dataset provided on Figshare. If you wish to start from the beginning, please complete tutorial Preprocessing data from oddball paradigm first.

[2]:

subject = "sample_subject"
subjects_dir = "subjects"
mne.set_config('SUBJECTS_DIR', subjects_dir)

# Reading mne.Epochs
epoched = mne.read_epochs(os.path.join(subjects_dir, subject, "_eeg", "_pre", f"{subject}_oddball-epo.fif"))
forward_path = os.path.join(subjects_dir, subject, "forward", f"{subject}_ico4-fwd.fif")
trans_path = os.path.join(subjects_dir, subject, "_eeg", "trans", f"{subject}-fit_trans.fif")

# We will be using only data for 'target' stimuli
target = epoched['target']
sel_epoched = target.copy()
sel_epoched = sel_epoched.set_eeg_reference('average', projection=True)
sel_epoched.apply_proj()
sel_evoked = sel_epoched.average()

Reading /Volumes/UMK/oddball/subjects/sample_subject/_eeg/_pre/sample_subject_oddball-epo.fif ...
    Found the data of interest:
        t =    -199.22 ...     800.78 ms
        0 CTF compensation matrices available
Not setting metadata
621 matching events found
No baseline correction applied
0 projection items activated
EEG channel type selected for re-referencing
Adding average EEG reference projection.
1 projection items deactivated
Average reference projection was added, but has not been applied yet. Use the apply_proj method to apply it.
Created an SSP operator (subspace dimension = 1)
1 projection items activated
SSP projectors applied...

To perform source localization, we need a forward model that links activity at source locations to the sensors (in this case EEG channels). Here, we load the forward solution and convert it to a fixed-orientation representation.

[3]:

# Reading mne.Forward
fwd_vector = mne.read_forward_solution(forward_path)

# Using fixed orientation in forward solution
fwd = mne.convert_forward_solution(
    fwd_vector,
    surf_ori=True,
    force_fixed=True,
    use_cps=True
)

# Leadfield matrix
leadfield = fwd["sol"]["data"]

# Source positions extracted from forward model
src = fwd["src"]

Reading forward solution from /Volumes/UMK/oddball/subjects/sample_subject/forward/sample_subject_ico4-fwd.fif...
    Reading a source space...
    [done]
    Reading a source space...
    [done]
    2 source spaces read
    Desired named matrix (kind = 3523 (FIFF_MNE_FORWARD_SOLUTION_GRAD)) not available
    Read EEG forward solution (5124 sources, 128 channels, free orientations)
    Source spaces transformed to the forward solution coordinate frame
    No patch info available. The standard source space normals will be employed in the rotation to the local surface coordinates....
    Changing to fixed-orientation forward solution with surface-based source orientations...
    [done]

“Sensory” processing#

We will start with analysing processes in “sensory” time window.

In an oddball paradigm, participants are presented with a sequence of frequent (standard) and infrequent (target) stimuli. The early neural responses to these target stimuli reflect sensory processing - the brain’s initial registration of the incoming stimulus before higher-level cognitive mechanisms are engaged. We assume that sensory processing for given oddball paradigm occurs within the 50-200 ms window after the stimuli. We will therefore compute the data covariance in this time range. To estimate the noise covariance, we use a baseline period -200ms to 0 ms, i.e., the interval before stimulus onset. This baseline is assumed to be free of stimulus-locked activity and provides reference for separating signal from noise.

[4]:

# Compute noise covariance
noise_cov = mne.compute_covariance(
    sel_epoched,
    tmin=-0.2,
    tmax=0,
    method="empirical"
)

# Compute data covariance for range corresponding to sensory processing
data_cov_sen = mne.compute_covariance(
    sel_epoched,
    tmin=0.05,
    tmax=0.2,
    method="empirical"
)

# Subset signal for given time range
signal_sen = sel_evoked.copy().crop(
    tmin=0.05,
    tmax=0.2
)

    Created an SSP operator (subspace dimension = 1)
    Setting small EEG eigenvalues to zero (without PCA)
Reducing data rank from 128 -> 127
Estimating covariance using EMPIRICAL
Done.
Number of samples used : 4056
[done]
    Created an SSP operator (subspace dimension = 1)
    Setting small EEG eigenvalues to zero (without PCA)
Reducing data rank from 128 -> 127
Estimating covariance using EMPIRICAL
Done.
Number of samples used : 3042
[done]

\(RN^{-1}\) eigenvalues analysis#

Before localization, we examine the eigenvalue spectrum of the whitened data covariance (RN⁻¹). Large eigenvalues indicate signal-dominated subspace; small ones reflect noise.

The plot shows:

  • Eigenvalues vs. rank: Helps determine optimal number of sources

  • Signal/noise separation: Drop-off indicates transition between signal and noise eigenvalues

For technical details on eigenvalue-based rank selection, see the MVPURE paper.

[5]:

fig_RN = viz.plot_RN_eigenvalues(
    R=data_cov_sen.data,
    N=noise_cov.data,
    subject=subject,
    s=14
)
../_images/tutorials_source_localization_and_reconstruction_on_data_from_oddball_paradigm_11_0.png

Localize#

Based on the eigenvalue spectrum above, we will localize \(l_0=5\) sources. We will use function mvpure_py.localizer.localize, which performs the actual source localization. The main parameters are:

  • subject: the subject ID (here: "sample_subject")

  • subjects_dir: directory containing the subject folders.

  • localizer_to_use: the algorithm variant. Here we choose "mai_mvp". For details, see paper or the function documentation.

  • n_sources_to_localize: number of sources to localize.

  • R: data covariance matrix

  • N: noise covariance matrix

  • forward: the mne.Forward object for this subject

  • r: optimization rank parameter. We will use list of ranks, therefore finally localized sources will be vertices localized in at least one iteration.

[6]:

l_0 = 5

locs_sen = localizer.localize(
    subject=subject,
    subjects_dir=subjects_dir,
    localizer_to_use=["mai_mvp"],
    n_sources_to_localize=l_0,
    R=data_cov_sen.data,
    N=noise_cov.data,
    forward=fwd,
    r=list(range(1,  l_0+1))
)

Rank: 1
Calculating activity index for localizer: mai_mvp

iter 1/5: 100%|██████████| 5124/5124 [00:00<00:00, 411866.65it/s]
iter 2/5: 100%|██████████| 5124/5124 [00:02<00:00, 2470.25it/s]
iter 3/5: 100%|██████████| 5124/5124 [00:00<00:00, 10469.17it/s]
iter 4/5: 100%|██████████| 5124/5124 [00:00<00:00, 11825.36it/s]
iter 5/5: 100%|██████████| 5124/5124 [00:00<00:00, 11814.32it/s]


[Activity Index Result]
  Selected indices (index_max): [31, 2697, 2042, 2747, 2554]
  Index max values: [1.21021186 1.45630239 1.65413543 1.8810501  1.98695636]
  Rank parameter (r): 1

Rank: 2
Calculating activity index for localizer: mai_mvp

iter 1/5: 100%|██████████| 5124/5124 [00:00<00:00, 804778.64it/s]
iter 2/5: 100%|██████████| 5124/5124 [00:00<00:00, 23801.55it/s]
iter 3/5: 100%|██████████| 5124/5124 [00:00<00:00, 18281.06it/s]
iter 4/5: 100%|██████████| 5124/5124 [00:00<00:00, 12398.25it/s]
iter 5/5: 100%|██████████| 5124/5124 [00:00<00:00, 9120.26it/s]


[Activity Index Result]
  Selected indices (index_max): [31, 2597, 4765, 2029, 2779]
  Index max values: [1.21021186 2.22236888 2.50811308 2.7559683  2.98024038]
  Rank parameter (r): 2

Rank: 3
Calculating activity index for localizer: mai_mvp

iter 1/5: 100%|██████████| 5124/5124 [00:00<00:00, 1167451.45it/s]
iter 2/5: 100%|██████████| 5124/5124 [00:00<00:00, 25179.68it/s]
iter 3/5: 100%|██████████| 5124/5124 [00:00<00:00, 24026.24it/s]
iter 4/5: 100%|██████████| 5124/5124 [00:00<00:00, 13863.52it/s]
iter 5/5: 100%|██████████| 5124/5124 [00:00<00:00, 11024.90it/s]


[Activity Index Result]
  Selected indices (index_max): [31, 2597, 1837, 4836, 1622]
  Index max values: [1.21021186 2.22236888 3.19893305 3.49029219 3.76719894]
  Rank parameter (r): 3

Rank: 4
Calculating activity index for localizer: mai_mvp

iter 1/5: 100%|██████████| 5124/5124 [00:00<00:00, 1070139.61it/s]
iter 2/5: 100%|██████████| 5124/5124 [00:00<00:00, 24746.92it/s]
iter 3/5: 100%|██████████| 5124/5124 [00:00<00:00, 23731.77it/s]
iter 4/5: 100%|██████████| 5124/5124 [00:00<00:00, 18485.44it/s]
iter 5/5: 100%|██████████| 5124/5124 [00:00<00:00, 9089.94it/s]


[Activity Index Result]
  Selected indices (index_max): [31, 2597, 1837, 3557, 2042]
  Index max values: [1.21021186 2.22236888 3.19893305 4.14451814 4.49089945]
  Rank parameter (r): 4

Rank: 5
Calculating activity index for localizer: mai_mvp

iter 1/5: 100%|██████████| 5124/5124 [00:00<00:00, 1166374.35it/s]
iter 2/5: 100%|██████████| 5124/5124 [00:00<00:00, 19519.41it/s]
iter 3/5: 100%|██████████| 5124/5124 [00:00<00:00, 21474.82it/s]
iter 4/5: 100%|██████████| 5124/5124 [00:00<00:00, 16581.93it/s]
iter 5/5: 100%|██████████| 5124/5124 [00:00<00:00, 12651.29it/s]


[Activity Index Result]
  Selected indices (index_max): [31, 2597, 1837, 3557, 1966]
  Index max values: [1.21021186 2.22236888 3.19893305 4.14451814 5.04837914]
  Rank parameter (r): 5



[7]:

# Optional visualization

# locs_sen.plot_localized_sources()

Signal reconstruction#

With sources localized, we reconstruct their time courses using LCMV beamforming:

  1. Restrict forward model to localized vertices only

  2. Build LCMV filters using unit-noise-gain normalization

  3. Apply filters to evoked data

[8]:

new_fwd = utils.subset_forward(
    old_fwd=fwd,
    localized=locs_sen,
    hemi="both"
)

[9]:

lcmv = mne.beamformer.make_lcmv(
    signal_sen.info,
    new_fwd,
    data_cov_sen,
    reg=0.05,
    noise_cov=noise_cov,
    pick_ori=None,  # not needed with fixed orientation forward
    weight_norm="nai",
    rank=None
)
stc_lcmv = mne.beamformer.apply_lcmv(signal_sen, lcmv)

Computing rank from covariance with rank=None
    Using tolerance 4.4e-13 (2.2e-16 eps * 128 dim * 15  max singular value)
    Estimated rank (eeg): 86
    EEG: rank 86 computed from 128 data channels with 1 projector
Computing rank from covariance with rank=None
    Using tolerance 3.6e-13 (2.2e-16 eps * 128 dim * 13  max singular value)
    Estimated rank (eeg): 86
    EEG: rank 86 computed from 128 data channels with 1 projector
Making LCMV beamformer with rank {'eeg': 86}
Computing inverse operator with 128 channels.
    128 out of 128 channels remain after picking
Selected 128 channels
Whitening the forward solution.
    Created an SSP operator (subspace dimension = 1)
Computing rank from covariance with rank={'eeg': 86}
    Setting small EEG eigenvalues to zero (without PCA)
Creating the source covariance matrix
Adjusting source covariance matrix.
Computing beamformer filters for 14 sources
Filter computation complete

Finally, let’s plot the localized sources with their reconstructed activity:

[10]:

viz.plot_sources_with_activity(
    subject=subject,
    stc=stc_lcmv
)

Using control points [0.92516798 1.13740896 1.58862482]

[10]:

<mne.viz._brain._brain.Brain at 0x128da7380>

Comparison with standard LCMV beamformer and sLORETA#

We can also compare the results of our MVPURE-based localization and reconstruction with standard LCMV beamformer and sLORETA methods.

Standard LCMV beamformer#

[11]:

standard_lcmv_sen = mne.beamformer.make_lcmv(
    signal_sen.info,
    fwd,
    data_cov_sen,
    reg=0.05,
    noise_cov=noise_cov,
    pick_ori=None,
    weight_norm="nai",
    rank=None
)
stc_standard_lcmv_sen = mne.beamformer.apply_lcmv(signal_sen, standard_lcmv_sen)

# Optional plot
stc_standard_lcmv_sen.plot(hemi="both")

Computing rank from covariance with rank=None
    Using tolerance 4.4e-13 (2.2e-16 eps * 128 dim * 15  max singular value)
    Estimated rank (eeg): 86
    EEG: rank 86 computed from 128 data channels with 1 projector
Computing rank from covariance with rank=None
    Using tolerance 3.6e-13 (2.2e-16 eps * 128 dim * 13  max singular value)
    Estimated rank (eeg): 86
    EEG: rank 86 computed from 128 data channels with 1 projector
Making LCMV beamformer with rank {'eeg': 86}
Computing inverse operator with 128 channels.
    128 out of 128 channels remain after picking
Selected 128 channels
Whitening the forward solution.
    Created an SSP operator (subspace dimension = 1)
Computing rank from covariance with rank={'eeg': 86}
    Setting small EEG eigenvalues to zero (without PCA)
Creating the source covariance matrix
Adjusting source covariance matrix.
Computing beamformer filters for 5124 sources
Filter computation complete
Using control points [0.44643495 0.51855953 1.26190296]

[11]:

<mne.viz._brain._brain.Brain at 0x1368ff390>

sLORETA#

To use sLORETA, we need to compute the inverse operator first:

[12]:

inv_op_sen = mne.minimum_norm.make_inverse_operator(
    info=signal_sen.info,
    forward=fwd,
    noise_cov=noise_cov,
    loose=0.0,  # forward with fixed orientation
    depth=None
)

Computing inverse operator with 128 channels.
    128 out of 128 channels remain after picking
Selected 128 channels
Whitening the forward solution.
    Created an SSP operator (subspace dimension = 1)
Computing rank from covariance with rank=None
    Using tolerance 3.6e-13 (2.2e-16 eps * 128 dim * 13  max singular value)
    Estimated rank (eeg): 86
    EEG: rank 86 computed from 128 data channels with 1 projector
    Setting small EEG eigenvalues to zero (without PCA)
Creating the source covariance matrix
Adjusting source covariance matrix.
Computing SVD of whitened and weighted lead field matrix.
    largest singular value = 3.46898
    scaling factor to adjust the trace = 1.46719e+20 (nchan = 128 nzero = 42)

To obtain a data-driven regularization parameter, the SNR was estimated as the ratio of the RMS signal amplitude during the response window to the baseline standard deviation. The inverse regularization parameter was then computed as lambda2 = 1 / snr**2, following the convention used in MNE-Pyhton.

[14]:

baseline = sel_evoked.copy().crop(tmin=sel_epoched.tmin, tmax=0)
active = sel_evoked.copy().crop(tmin=0.05, tmax=0.2)

noise_std = np.std(baseline.data)
signal_rms = np.sqrt(np.mean(active.data**2))

snr_est = np.round(signal_rms / noise_std, 2)
print(f"SNR: {snr_est}")

SNR: 2.55

[15]:

# Compute sLORETA source estimate
stc_sloreta_sen = mne.minimum_norm.apply_inverse(
    signal_sen,
    inv_op_sen,
    lambda2=1.0 / snr_est ** 2,
    method="sLORETA"
)

# Optional plot
stc_sloreta_sen.plot(hemi="both")

Preparing the inverse operator for use...
    Scaled noise and source covariance from nave = 1 to nave = 78
    Created the regularized inverter
    Created an SSP operator (subspace dimension = 1)
    Created the whitener using a noise covariance matrix with rank 86 (42 small eigenvalues omitted)
    Computing noise-normalization factors (sLORETA)...
[done]
Applying inverse operator to "target"...
    Picked 128 channels from the data
    Computing inverse...
    Eigenleads need to be weighted ...
    Computing residual...
    Explained  81.1% variance
    sLORETA...
[done]
Using control points [2.08870182 2.46848157 5.19074348]

[15]:

<mne.viz._brain._brain.Brain at 0x1368c2990>

“Cognitive” task#

After examining the early sensory responses, we now turn to the later cognitive stage of processing in the oddball paradigm. In EEG, target stimuli typically evoke a P300 component — a positive deflection peaking around 300–600 ms after stimulus onset. This response is thought to reflect higher-level cognitive processes, such as attention allocation and stimulus evaluation, in contrast to the earlier sensory responses.

For this dataset, we will therefore define the cognitive time window as 350–600 ms. The pipeline remains the same as before:

  • Compute noise covariance (always from −200 to 0 ms).

  • Compute data covariance in the cognitive window (350–600 ms).

  • Subset the evoked signal to this time range.

[16]:

# Compute data covariance for range corresponding to sensory processing
data_cov_task = mne.compute_covariance(
    sel_epoched,
    tmin=0.35,
    tmax=0.6,
    method="empirical"
)
# There's no need to compute `noise_covariance` again as it is the same time interval

sel_evoked = sel_epoched.average()

# Subset signal for given time range
signal_task = sel_evoked.copy().crop(
    tmin=0.35,
    tmax=0.6
)

    Created an SSP operator (subspace dimension = 1)
    Setting small EEG eigenvalues to zero (without PCA)
Reducing data rank from 128 -> 127
Estimating covariance using EMPIRICAL
Done.
Number of samples used : 5070
[done]

From here, we can repeat the same steps as in the sensory section:

  • analyze eigenvalues of \(RN^{-1}\),

  • localize sources,

  • reconstruct signals with LCMV-NAI beamformer,

  • and finally visualize the results.

[17]:

fig_RN_task = viz.plot_RN_eigenvalues(
    R=data_cov_task.data,
    N=noise_cov.data,
    subject=subject,
    s=14
)
../_images/tutorials_source_localization_and_reconstruction_on_data_from_oddball_paradigm_31_0.png 
[18]:

l_0 = 10

# Localize
locs_task = localizer.localize(
    subject=subject,
    subjects_dir=subjects_dir,
    localizer_to_use=["mai_mvp"],
    n_sources_to_localize=l_0,
    R=data_cov_task.data,
    N=noise_cov.data,
    forward=fwd,
    r=list(range(1, l_0 + 1)),
)

Rank: 1
Calculating activity index for localizer: mai_mvp

iter 1/10: 100%|██████████| 5124/5124 [00:00<00:00, 888964.83it/s]
iter 2/10: 100%|██████████| 5124/5124 [00:00<00:00, 8378.00it/s]
iter 3/10: 100%|██████████| 5124/5124 [00:00<00:00, 6862.69it/s]
iter 4/10: 100%|██████████| 5124/5124 [00:00<00:00, 7467.71it/s]
iter 5/10: 100%|██████████| 5124/5124 [00:00<00:00, 7979.12it/s]
iter 6/10: 100%|██████████| 5124/5124 [00:00<00:00, 11787.54it/s]
iter 7/10: 100%|██████████| 5124/5124 [00:00<00:00, 14687.92it/s]
iter 8/10: 100%|██████████| 5124/5124 [00:00<00:00, 12704.20it/s]
iter 9/10: 100%|██████████| 5124/5124 [00:00<00:00, 11896.28it/s]
iter 10/10: 100%|██████████| 5124/5124 [00:00<00:00, 12578.46it/s]


[Activity Index Result]
  Selected indices (index_max): [2160, 4755, 419, 1240, 996, 1546, 3978, 3790, 3861, 441]
  Index max values: [1.77758988 2.25869353 2.55953156 2.85677162 2.98744644 3.06891462
 3.13975703 3.20066062 3.3391399  3.48255084]
  Rank parameter (r): 1

Rank: 2
Calculating activity index for localizer: mai_mvp

iter 1/10: 100%|██████████| 5124/5124 [00:00<00:00, 993648.05it/s]
iter 2/10: 100%|██████████| 5124/5124 [00:00<00:00, 27620.99it/s]
iter 3/10: 100%|██████████| 5124/5124 [00:00<00:00, 21019.72it/s]
iter 4/10: 100%|██████████| 5124/5124 [00:00<00:00, 13141.39it/s]
iter 5/10: 100%|██████████| 5124/5124 [00:00<00:00, 10400.24it/s]
iter 6/10: 100%|██████████| 5124/5124 [00:00<00:00, 10018.55it/s]
iter 7/10: 100%|██████████| 5124/5124 [00:00<00:00, 14864.97it/s]
iter 8/10: 100%|██████████| 5124/5124 [00:01<00:00, 5113.66it/s]
iter 9/10: 100%|██████████| 5124/5124 [00:00<00:00, 6746.69it/s]
iter 10/10: 100%|██████████| 5124/5124 [00:00<00:00, 8054.88it/s]


[Activity Index Result]
  Selected indices (index_max): [2160, 2300, 4897, 4755, 419, 1808, 2461, 1240, 2248, 2834]
  Index max values: [1.77758988 3.2754443  4.04566418 4.46424172 4.79131535 5.08992702
 5.34796421 5.62534082 5.86618604 6.15678146]
  Rank parameter (r): 2

Rank: 3
Calculating activity index for localizer: mai_mvp

iter 1/10: 100%|██████████| 5124/5124 [00:00<00:00, 990853.56it/s]
iter 2/10: 100%|██████████| 5124/5124 [00:00<00:00, 15807.20it/s]
iter 3/10: 100%|██████████| 5124/5124 [00:00<00:00, 21354.96it/s]
iter 4/10: 100%|██████████| 5124/5124 [00:00<00:00, 16599.34it/s]
iter 5/10: 100%|██████████| 5124/5124 [00:00<00:00, 9891.17it/s]
iter 6/10: 100%|██████████| 5124/5124 [00:00<00:00, 12161.66it/s]
iter 7/10: 100%|██████████| 5124/5124 [00:00<00:00, 15177.84it/s]
iter 8/10: 100%|██████████| 5124/5124 [00:00<00:00, 13158.00it/s]
iter 9/10: 100%|██████████| 5124/5124 [00:00<00:00, 11424.85it/s]
iter 10/10: 100%|██████████| 5124/5124 [00:00<00:00, 12800.40it/s]


[Activity Index Result]
  Selected indices (index_max): [2160, 2300, 4755, 1368, 4897, 2224, 1832, 1843, 419, 769]
  Index max values: [1.77758988 3.2754443  4.85038569 5.52222705 6.10691151 6.45017419
 6.99421624 7.34994298 7.64392056 7.89816082]
  Rank parameter (r): 3

Rank: 4
Calculating activity index for localizer: mai_mvp

iter 1/10: 100%|██████████| 5124/5124 [00:00<00:00, 1302522.04it/s]
iter 2/10: 100%|██████████| 5124/5124 [00:00<00:00, 27756.58it/s]
iter 3/10: 100%|██████████| 5124/5124 [00:00<00:00, 23749.40it/s]
iter 4/10: 100%|██████████| 5124/5124 [00:00<00:00, 18777.50it/s]
iter 5/10: 100%|██████████| 5124/5124 [00:00<00:00, 11719.40it/s]
iter 6/10: 100%|██████████| 5124/5124 [00:00<00:00, 10106.85it/s]
iter 7/10: 100%|██████████| 5124/5124 [00:00<00:00, 14545.87it/s]
iter 8/10: 100%|██████████| 5124/5124 [00:00<00:00, 12244.84it/s]
iter 9/10: 100%|██████████| 5124/5124 [00:00<00:00, 13265.75it/s]
iter 10/10: 100%|██████████| 5124/5124 [00:00<00:00, 11298.09it/s]


[Activity Index Result]
  Selected indices (index_max): [2160, 2300, 4755, 1368, 4897, 1986, 419, 1843, 2554, 4147]
  Index max values: [1.77758988 3.2754443  4.85038569 6.4063041  7.08410596 7.67033714
 8.09155456 8.48531314 8.94506573 9.48490676]
  Rank parameter (r): 4

Rank: 5
Calculating activity index for localizer: mai_mvp

iter 1/10: 100%|██████████| 5124/5124 [00:00<00:00, 1180403.89it/s]
iter 2/10: 100%|██████████| 5124/5124 [00:00<00:00, 29186.60it/s]
iter 3/10: 100%|██████████| 5124/5124 [00:00<00:00, 26246.38it/s]
iter 4/10: 100%|██████████| 5124/5124 [00:00<00:00, 16430.39it/s]
iter 5/10: 100%|██████████| 5124/5124 [00:00<00:00, 11812.64it/s]
iter 6/10: 100%|██████████| 5124/5124 [00:00<00:00, 12403.11it/s]
iter 7/10: 100%|██████████| 5124/5124 [00:00<00:00, 15545.31it/s]
iter 8/10: 100%|██████████| 5124/5124 [00:00<00:00, 10017.54it/s]
iter 9/10: 100%|██████████| 5124/5124 [00:00<00:00, 12525.49it/s]
iter 10/10: 100%|██████████| 5124/5124 [00:00<00:00, 10419.19it/s]


[Activity Index Result]
  Selected indices (index_max): [2160, 2300, 4755, 1368, 1689, 4897, 1967, 2037, 1121, 4587]
  Index max values: [ 1.77758988  3.2754443   4.85038569  6.4063041   7.90005161  8.56096489
  9.15113909  9.62022658 10.09576615 10.54620627]
  Rank parameter (r): 5

Rank: 6
Calculating activity index for localizer: mai_mvp

iter 1/10: 100%|██████████| 5124/5124 [00:00<00:00, 1424040.13it/s]
iter 2/10: 100%|██████████| 5124/5124 [00:00<00:00, 27496.21it/s]
iter 3/10: 100%|██████████| 5124/5124 [00:00<00:00, 26382.74it/s]
iter 4/10: 100%|██████████| 5124/5124 [00:00<00:00, 19414.53it/s]
iter 5/10: 100%|██████████| 5124/5124 [00:00<00:00, 10836.60it/s]
iter 6/10: 100%|██████████| 5124/5124 [00:00<00:00, 15466.34it/s]
iter 7/10: 100%|██████████| 5124/5124 [00:00<00:00, 15021.11it/s]
iter 8/10: 100%|██████████| 5124/5124 [00:00<00:00, 12913.21it/s]
iter 9/10: 100%|██████████| 5124/5124 [00:00<00:00, 11685.79it/s]
iter 10/10: 100%|██████████| 5124/5124 [00:00<00:00, 12603.99it/s]


[Activity Index Result]
  Selected indices (index_max): [2160, 2300, 4755, 1368, 1689, 1835, 3842, 1875, 4794, 5079]
  Index max values: [ 1.77758988  3.2754443   4.85038569  6.4063041   7.90005161  9.38758629
 10.06701259 10.68034168 11.27555351 11.83302435]
  Rank parameter (r): 6

Rank: 7
Calculating activity index for localizer: mai_mvp

iter 1/10: 100%|██████████| 5124/5124 [00:00<00:00, 1168784.73it/s]
iter 2/10: 100%|██████████| 5124/5124 [00:00<00:00, 18726.63it/s]
iter 3/10: 100%|██████████| 5124/5124 [00:00<00:00, 22185.37it/s]
iter 4/10: 100%|██████████| 5124/5124 [00:00<00:00, 17633.03it/s]
iter 5/10: 100%|██████████| 5124/5124 [00:00<00:00, 13379.93it/s]
iter 6/10: 100%|██████████| 5124/5124 [00:00<00:00, 11536.61it/s]
iter 7/10: 100%|██████████| 5124/5124 [00:00<00:00, 17132.76it/s]
iter 8/10: 100%|██████████| 5124/5124 [00:00<00:00, 14138.33it/s]
iter 9/10: 100%|██████████| 5124/5124 [00:00<00:00, 13097.38it/s]
iter 10/10: 100%|██████████| 5124/5124 [00:00<00:00, 10019.89it/s]


[Activity Index Result]
  Selected indices (index_max): [2160, 2300, 4755, 1368, 1689, 1835, 1875, 3842, 4973, 4569]
  Index max values: [ 1.77758988  3.2754443   4.85038569  6.4063041   7.90005161  9.38758629
 10.84390766 11.57196752 12.18799488 12.76289158]
  Rank parameter (r): 7

Rank: 8
Calculating activity index for localizer: mai_mvp

iter 1/10: 100%|██████████| 5124/5124 [00:00<00:00, 1447441.66it/s]
iter 2/10: 100%|██████████| 5124/5124 [00:00<00:00, 30162.82it/s]
iter 3/10: 100%|██████████| 5124/5124 [00:00<00:00, 24891.32it/s]
iter 4/10: 100%|██████████| 5124/5124 [00:00<00:00, 20275.80it/s]
iter 5/10: 100%|██████████| 5124/5124 [00:00<00:00, 8525.79it/s]
iter 6/10: 100%|██████████| 5124/5124 [00:00<00:00, 15371.78it/s]
iter 7/10: 100%|██████████| 5124/5124 [00:00<00:00, 17809.91it/s]
iter 8/10: 100%|██████████| 5124/5124 [00:00<00:00, 17884.64it/s]
iter 9/10: 100%|██████████| 5124/5124 [00:00<00:00, 13870.51it/s]
iter 10/10: 100%|██████████| 5124/5124 [00:00<00:00, 12671.69it/s]


[Activity Index Result]
  Selected indices (index_max): [2160, 2300, 4755, 1368, 1689, 1835, 1875, 4405, 4701, 4231]
  Index max values: [ 1.77758988  3.2754443   4.85038569  6.4063041   7.90005161  9.38758629
 10.84390766 12.2470246  12.9417525  13.68085168]
  Rank parameter (r): 8

Rank: 9
Calculating activity index for localizer: mai_mvp

iter 1/10: 100%|██████████| 5124/5124 [00:00<00:00, 1189155.85it/s]
iter 2/10: 100%|██████████| 5124/5124 [00:00<00:00, 25884.03it/s]
iter 3/10: 100%|██████████| 5124/5124 [00:00<00:00, 23616.04it/s]
iter 4/10: 100%|██████████| 5124/5124 [00:00<00:00, 19388.75it/s]
iter 5/10: 100%|██████████| 5124/5124 [00:00<00:00, 13321.03it/s]
iter 6/10: 100%|██████████| 5124/5124 [00:00<00:00, 12150.82it/s]
iter 7/10: 100%|██████████| 5124/5124 [00:00<00:00, 15276.32it/s]
iter 8/10: 100%|██████████| 5124/5124 [00:00<00:00, 15621.20it/s]
iter 9/10: 100%|██████████| 5124/5124 [00:00<00:00, 10204.14it/s]
iter 10/10: 100%|██████████| 5124/5124 [00:00<00:00, 10814.94it/s]


[Activity Index Result]
  Selected indices (index_max): [2160, 2300, 4755, 1368, 1689, 1835, 1875, 4405, 600, 4701]
  Index max values: [ 1.77758988  3.2754443   4.85038569  6.4063041   7.90005161  9.38758629
 10.84390766 12.2470246  13.62102771 14.3351125 ]
  Rank parameter (r): 9

Rank: 10
Calculating activity index for localizer: mai_mvp

iter 1/10: 100%|██████████| 5124/5124 [00:00<00:00, 1307196.26it/s]
iter 2/10: 100%|██████████| 5124/5124 [00:00<00:00, 27305.15it/s]
iter 3/10: 100%|██████████| 5124/5124 [00:00<00:00, 22590.20it/s]
iter 4/10: 100%|██████████| 5124/5124 [00:00<00:00, 16476.62it/s]
iter 5/10: 100%|██████████| 5124/5124 [00:00<00:00, 9448.84it/s]
iter 6/10: 100%|██████████| 5124/5124 [00:00<00:00, 13026.82it/s]
iter 7/10: 100%|██████████| 5124/5124 [00:00<00:00, 16204.68it/s]
iter 8/10: 100%|██████████| 5124/5124 [00:00<00:00, 15458.13it/s]
iter 9/10: 100%|██████████| 5124/5124 [00:00<00:00, 13535.00it/s]
iter 10/10: 100%|██████████| 5124/5124 [00:00<00:00, 12961.10it/s]


[Activity Index Result]
  Selected indices (index_max): [2160, 2300, 4755, 1368, 1689, 1835, 1875, 4405, 600, 2431]
  Index max values: [ 1.77758988  3.2754443   4.85038569  6.4063041   7.90005161  9.38758629
 10.84390766 12.2470246  13.62102771 14.93465609]
  Rank parameter (r): 10



[19]:

# Optional visualization of localized sources
# locs_task.plot_localized_sources(hemi="both", scale_factor=0.7, scale_mapping=False)

[20]:

new_fwd_task = utils.subset_forward(
    old_fwd=fwd,
    localized=locs_task,
    hemi="both"
)
lcmv_task = mne.beamformer.make_lcmv(
    signal_task.info,
    new_fwd_task,
    data_cov_task,
    reg=0.05,
    noise_cov=noise_cov,
    pick_ori=None,  # not needed with fixed orientation forward
    weight_norm="nai",
    rank=None
)
stc_lcmv_task = mne.beamformer.apply_lcmv(signal_task, lcmv_task)

Computing rank from covariance with rank=None
    Using tolerance 5.2e-13 (2.2e-16 eps * 128 dim * 18  max singular value)
    Estimated rank (eeg): 86
    EEG: rank 86 computed from 128 data channels with 1 projector
Computing rank from covariance with rank=None
    Using tolerance 3.6e-13 (2.2e-16 eps * 128 dim * 13  max singular value)
    Estimated rank (eeg): 86
    EEG: rank 86 computed from 128 data channels with 1 projector
Making LCMV beamformer with rank {'eeg': 86}
Computing inverse operator with 128 channels.
    128 out of 128 channels remain after picking
Selected 128 channels
Whitening the forward solution.
    Created an SSP operator (subspace dimension = 1)
Computing rank from covariance with rank={'eeg': 86}
    Setting small EEG eigenvalues to zero (without PCA)
Creating the source covariance matrix
Adjusting source covariance matrix.
Computing beamformer filters for 41 sources
Filter computation complete

[21]:

# Plot localized sources with reconstructed activity
viz.plot_sources_with_activity(
    subject=subject,
    stc=stc_lcmv_task,
)

Using control points [0.59890943 0.65281948 1.00913456]

[21]:

<mne.viz._brain._brain.Brain at 0x1369a6490>

Comparison with standard LCMV beamformer and sLORETA#

Standard LCMV beamformer#

[22]:

standard_lcmv_task = mne.beamformer.make_lcmv(
    signal_task.info,
    fwd,
    data_cov_task,
    reg=0.05,
    noise_cov=noise_cov,
    pick_ori=None,
    weight_norm="nai",
    rank=None
)

stc_standard_lcmv_task = mne.beamformer.apply_lcmv(signal_task, standard_lcmv_task)

# Optional plot
# stc_standard_lcmv_task.plot(hemi="both")

Computing rank from covariance with rank=None
    Using tolerance 5.2e-13 (2.2e-16 eps * 128 dim * 18  max singular value)
    Estimated rank (eeg): 86
    EEG: rank 86 computed from 128 data channels with 1 projector
Computing rank from covariance with rank=None
    Using tolerance 3.6e-13 (2.2e-16 eps * 128 dim * 13  max singular value)
    Estimated rank (eeg): 86
    EEG: rank 86 computed from 128 data channels with 1 projector
Making LCMV beamformer with rank {'eeg': 86}
Computing inverse operator with 128 channels.
    128 out of 128 channels remain after picking
Selected 128 channels
Whitening the forward solution.
    Created an SSP operator (subspace dimension = 1)
Computing rank from covariance with rank={'eeg': 86}
    Setting small EEG eigenvalues to zero (without PCA)
Creating the source covariance matrix
Adjusting source covariance matrix.
Computing beamformer filters for 5124 sources
Filter computation complete
Using control points [0.59661989 0.65132542 1.08688352]

[22]:

<mne.viz._brain._brain.Brain at 0x13697a350>

sLORETA#

[23]:

inv_op_task = mne.minimum_norm.make_inverse_operator(
    info=signal_task.info,
    forward=fwd,
    noise_cov=noise_cov,
    loose=0.0,  # forward with fixed orientation
    depth=None
)

Computing inverse operator with 128 channels.
    128 out of 128 channels remain after picking
Selected 128 channels
Whitening the forward solution.
    Created an SSP operator (subspace dimension = 1)
Computing rank from covariance with rank=None
    Using tolerance 3.6e-13 (2.2e-16 eps * 128 dim * 13  max singular value)
    Estimated rank (eeg): 86
    EEG: rank 86 computed from 128 data channels with 1 projector
    Setting small EEG eigenvalues to zero (without PCA)
Creating the source covariance matrix
Adjusting source covariance matrix.
Computing SVD of whitened and weighted lead field matrix.
    largest singular value = 3.46898
    scaling factor to adjust the trace = 1.46719e+20 (nchan = 128 nzero = 42)

[24]:

baseline = sel_evoked.copy().crop(tmin=sel_epoched.tmin, tmax=0)
active = sel_evoked.copy().crop(tmin=0.35, tmax=0.6)

noise_std = np.std(baseline.data)
signal_rms = np.sqrt(np.mean(active.data**2))

snr_est = np.round(signal_rms / noise_std, 2)
print(f"SNR: {snr_est}")

SNR: 3.5

[26]:

stc_sloreta_task = mne.minimum_norm.apply_inverse(
    signal_task,
    inv_op_task,
    lambda2=1.0 / snr_est ** 2,
    method="sLORETA"
)

# Optional plot
stc_sloreta_task.plot(hemi="both")

Preparing the inverse operator for use...
    Scaled noise and source covariance from nave = 1 to nave = 78
    Created the regularized inverter
    Created an SSP operator (subspace dimension = 1)
    Created the whitener using a noise covariance matrix with rank 86 (42 small eigenvalues omitted)
    Computing noise-normalization factors (sLORETA)...
[done]
Applying inverse operator to "target"...
    Picked 128 channels from the data
    Computing inverse...
    Eigenleads need to be weighted ...
    Computing residual...
    Explained  86.3% variance
    sLORETA...
[done]
Using control points [1.98567814 2.18181795 3.69853747]

[26]:

<mne.viz._brain._brain.Brain at 0x1369a7b10>