Skip to content

Binder

Computation of the lead-field matrix

The lead-field matrix is a mathematical representation that describes the relationship between the electrical activity generated by neural sources in the brain and the resulting electrical measurements observed on the scalp's surface during EEG recordings.

The focus in this Notebook is on constructing a lead-field matrix that corresponds to the Automated Anatomical Labeling 2 (AAL2) atlas, a widely used brain parcellation scheme that divides the brain into regions based on anatomical landmarks. This atlas provides a comprehensive map of brain regions, allowing to associate EEG measurements with specific anatomical structures.

We utilized the standard 1020 EEG electrode configuration, which is a well-established system for scalp EEG recordings, involves employing a specific arrangement of electrode placement points on the scalp. By aligning the Automated Anatomical Labeling 2 (AAL2) atlas with the 1020 EEG electrode configuration, the primary aim of this Notebook is to establish a direct link between distinct brain regions defined by the AAL2 atlas and the specific EEG recording channels.

# Authors: Mohammad Orabe <orabe.mhd@gmail.com>
#          Zixuan liu <zixuan.liu@campus.tu-berlin.de> 
# 
# change to the root directory of the project
import os

if os.getcwd().split("/")[-1] == "examples":
    os.chdir("..")

# This will reload all imports as soon as the code changes
%load_ext autoreload
%autoreload 2
from importlib import import_module
import subprocess
import sys

named_libs = [
    ("matplotlib.pyplot", "plt"),
    ("nibabel", "nib"),
    ("mne", "mne"),
    ("neurolib", "neurolib"),
]

for name, short in named_libs:
    try:
        globals()[short] = import_module(name)
    except ImportError:
        subprocess.check_call([sys.executable, "-m", "pip", "install", name])
        globals()[short] = import_module(name)
from neurolib.utils.leadfield import LeadfieldGenerator
subject = "fsaverage"

atlas_xml_path = os.path.join(
    "examples", "data", "AAL2_atlas_data", "AAL2.xml"
)  # Load the NIfTI file of the AAL2 atlas

atlas_nii_path = os.path.join(
    "examples", "data", "AAL2_atlas_data", "AAL2.nii"
)  # Load the XML file of the AAL2 atlas.

Download and load the data

Download the template MRI "fsaverage", user-specific data can replace the dataset.

object = LeadfieldGenerator(subject)
object.load_data(subjects_dir=None, subject=subject)
object.load_transformation_file(trans_path=None, subject=subject)
0 files missing from root.txt in /home/ben/mne_data/MNE-fsaverage-data
0 files missing from bem.txt in /home/ben/mne_data/MNE-fsaverage-data/fsaverage
Load template data 'fsaverage'
Load default transformation file of 'fsaverage'

Set conductivity of tissues and compute BEM

Setting conductivity based on the type of the tissues, calculate the BEM of the data, further visualize BEM surfaces.

bem, plot_bem_kwargs_default = object.build_BEM(
    conductivity=(0.3, 0.006, 0.3),
    visualization=True,
    brain_surfaces="white",
    orientation="coronal",
    slices=[50, 100, 150, 200],
)
Creating the BEM geometry...
Going from 5th to 4th subdivision of an icosahedron (n_tri: 20480 -> 5120)
Going from 5th to 4th subdivision of an icosahedron (n_tri: 20480 -> 5120)
Going from 5th to 4th subdivision of an icosahedron (n_tri: 20480 -> 5120)
outer skin  CM is  -0.21 -19.38  -0.23 mm
outer skull CM is  -0.19 -19.34  -0.49 mm
inner skull CM is  -0.53 -21.10   6.21 mm
Checking that surface outer skull is inside surface outer skin  ...
Checking that surface inner skull is inside surface outer skull ...
Checking distance between outer skin  and outer skull surfaces...
Minimum distance between the outer skin  and outer skull surfaces is approximately    1.6 mm
Checking distance between outer skull and inner skull surfaces...
Minimum distance between the outer skull and inner skull surfaces is approximately    5.4 mm
Surfaces passed the basic topology checks.
Complete.

Three-layer model surfaces loaded.
Computing the linear collocation solution...
    Matrix coefficients...
        outer skin  (2562) -> outer skin  (2562) ...
        outer skin  (2562) -> outer skull (2562) ...
        outer skin  (2562) -> inner skull (2562) ...
        outer skull (2562) -> outer skin  (2562) ...
        outer skull (2562) -> outer skull (2562) ...
        outer skull (2562) -> inner skull (2562) ...
        inner skull (2562) -> outer skin  (2562) ...
        inner skull (2562) -> outer skull (2562) ...
        inner skull (2562) -> inner skull (2562) ...
    Inverting the coefficient matrix...
IP approach required...
    Matrix coefficients (homog)...
        inner skull (2562) -> inner skull (2562) ...
    Inverting the coefficient matrix (homog)...
    Modify the original solution to incorporate IP approach...
        Combining...
        Scaling...
Solution ready.
BEM geometry computations complete.
Using surface: /home/ben/mne_data/MNE-fsaverage-data/fsaverage/bem/inner_skull.surf
Using surface: /home/ben/mne_data/MNE-fsaverage-data/fsaverage/bem/outer_skull.surf
Using surface: /home/ben/mne_data/MNE-fsaverage-data/fsaverage/bem/outer_skin.surf

/home/ben/anaconda3/envs/mne_py39/lib/python3.9/site-packages/mne/viz/utils.py:151: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.
  (fig or plt).show(**kwargs)

Generate surface source space

The method to achieve AAL atlas source space is to calculate the average dipole value of all dipoles in each annotation of the atlas, so at first place, the surface source space needed to be generated. The data "fsaverage" has a pre build-up surface source space, nevertheless, the source space can also be calculated.

src = object.generate_surface_source_space(
    plot_bem_kwargs_default, spacing="ico4", add_dist="patch", visualization=True
)
Using surface: /home/ben/mne_data/MNE-fsaverage-data/fsaverage/bem/inner_skull.surf
Using surface: /home/ben/mne_data/MNE-fsaverage-data/fsaverage/bem/outer_skull.surf
Using surface: /home/ben/mne_data/MNE-fsaverage-data/fsaverage/bem/outer_skin.surf
    Reading a source space...
    [done]
    Reading a source space...
    [done]
    2 source spaces read

/home/ben/anaconda3/envs/mne_py39/lib/python3.9/site-packages/mne/viz/utils.py:151: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.
  (fig or plt).show(**kwargs)

Load and coregistrate standard EEG configuration

The standard 1020 EEG electrode locations are already calculated in fsaverage's space (MNI space)

raw = object.EEG_coregistration(src, configuration="standard_1020", visualization=False)
Extracting EDF parameters from /home/ben/mne_data/MNE-eegbci-data/files/eegmmidb/1.0.0/S001/S001R06.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19999  =      0.000 ...   124.994 secs...
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.

Calculate the general forward solution

The overall forward solution between surface source model and standard EEG montage is computed

fwd = object.calculate_general_forward_solution(raw, src, bem, eeg=True, mindist=5.0)
Source space          : /home/ben/mne_data/MNE-fsaverage-data/fsaverage/bem/fsaverage-ico-5-src.fif
MRI -> head transform : /home/ben/mne_data/MNE-fsaverage-data/fsaverage/bem/fsaverage-trans.fif
Measurement data      : instance of Info
Conductor model   : instance of ConductorModel
Accurate field computations
Do computations in head coordinates
Free source orientations

Reading /home/ben/mne_data/MNE-fsaverage-data/fsaverage/bem/fsaverage-ico-5-src.fif...
Read 2 source spaces a total of 20484 active source locations

Coordinate transformation: MRI (surface RAS) -> head
     0.999994  0.003552  0.000202      -1.76 mm
    -0.003558  0.998389  0.056626      31.09 mm
    -0.000001 -0.056626  0.998395      39.60 mm
     0.000000  0.000000  0.000000       1.00

Read  64 EEG channels from info
Head coordinate coil definitions created.
Source spaces are now in head coordinates.

Employing the head->MRI coordinate transform with the BEM model.
BEM model instance of ConductorModel is now set up

Source spaces are in head coordinates.
Checking that the sources are inside the surface and at least    5.0 mm away (will take a few...)
Checking surface interior status for 10242 points...
    Found  2433/10242 points inside  an interior sphere of radius   47.7 mm
    Found     0/10242 points outside an exterior sphere of radius   98.3 mm
    Found     0/ 7809 points outside using surface Qhull
    Found     0/ 7809 points outside using solid angles
    Total 10242/10242 points inside the surface
Interior check completed in 1602.6 ms
Checking surface interior status for 10242 points...
    Found  2241/10242 points inside  an interior sphere of radius   47.7 mm
    Found     0/10242 points outside an exterior sphere of radius   98.3 mm
    Found     0/ 8001 points outside using surface Qhull
    Found     0/ 8001 points outside using solid angles
    Total 10242/10242 points inside the surface
Interior check completed in 1642.3 ms

Setting up for EEG...
Computing EEG at 20484 source locations (free orientations)...

Finished.
The general forward solution: <Forward | MEG channels: 0 | EEG channels: 64 | Source space: Surface with 20484 vertices | Source orientation: Free>
=====================================================

Downsample the forward solution to achieve lead-field matrix

With the forward solution that being calculated above, compute the average dipole value of the dipoles in each AAL atlas to acquire the lead-field matrix.

This part is majorly the contribution from Dr Nikola Jajcay and Martin Krück

leadfield_downsampled, unique_labels = object.compute_downsampled_leadfield(
    fwd,
    atlas_nii_path=atlas_nii_path,
    atlas_xml_path=atlas_xml_path,
    atlas="aal2_cortical",
    cortex_parts="only_cortical_parts",
    path_to_save=None,
)
    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]

WARNING:root:Atlas doesn't start at 0, reindexing...

Lead-field matrix's size : 64 sensors x 80 dipoles
=====================================================
Downsampled lead-field matrix: [[ 46.0499115   -8.10420036  51.77446747 ...  -2.94830823 -23.09635735
  -19.40752983]
 [ 63.07761765  -3.93809652  81.55727386 ...  -4.3759985  -29.43257332
  -22.51317406]
 [ 45.61634445   4.58903265  83.16187286 ...  -5.45599174 -29.02531624
  -24.48542023]
 ...
 [-17.10632133 -12.98373508 -18.83823586 ... -17.25921631   6.94994116
    5.47431898]
 [-17.82714272 -12.14234829 -22.45002937 ... -18.51730156   3.97969389
   17.09500885]
 [-18.07602882 -14.23802567 -20.28145409 ... -16.23429108  19.68803215
   16.43236351]]
=====================================================

Investigation of the missing regions of AAL atlas

Can be implemented if the downsample lead-field matrix does possess not an anticipated size.

object.check_atlas_missing_regions(
    atlas_xml_path=atlas_xml_path, unique_labels=unique_labels
)
WARNING:root:Atlas doesn't start at 0, reindexing...

total region quantity: 120
missed region quantity:  40
=====================================================
missed region labels: [4101 4102 4111 4112 4201 4202 7001 7002 7011 7012 7021 7022 7101 7102
 9001 9002 9011 9012 9021 9022 9031 9032 9041 9042 9051 9052 9061 9062
 9071 9072 9081 9082 9100 9110 9120 9130 9140 9150 9160 9170]
=====================================================
missed region names: ['Hippocampus_L', 'Hippocampus_R', 'ParaHippocampal_L', 'ParaHippocampal_R', 'Amygdala_L', 'Amygdala_R', 'Caudate_L', 'Caudate_R', 'Putamen_L', 'Putamen_R', 'Pallidum_L', 'Pallidum_R', 'Thalamus_L', 'Thalamus_R', 'Cerebelum_Crus1_L', 'Cerebelum_Crus1_R', 'Cerebelum_Crus2_L', 'Cerebelum_Crus2_R', 'Cerebelum_3_L', 'Cerebelum_3_R', 'Cerebelum_4_5_L', 'Cerebelum_4_5_R', 'Cerebelum_6_L', 'Cerebelum_6_R', 'Cerebelum_7b_L', 'Cerebelum_7b_R', 'Cerebelum_8_L', 'Cerebelum_8_R', 'Cerebelum_9_L', 'Cerebelum_9_R', 'Cerebelum_10_L', 'Cerebelum_10_R', 'Vermis_1_2', 'Vermis_3', 'Vermis_4_5', 'Vermis_6', 'Vermis_7', 'Vermis_8', 'Vermis_9', 'Vermis_10']
=====================================================
missed region index: [ 41  42  43  44  45  46  75  76  77  78  79  80  81  82  95  96  97  98
  99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116
 117 118 119 120]
=====================================================