Head model and forward computation

The aim of this tutorial is to be a getting started for forward computation.

For more extensive details and presentation of the general concepts for forward modeling. See The forward solution.

import mne
from mne.datasets import sample
data_path = sample.data_path()

# the raw file containing the channel location + types
raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
# The paths to freesurfer reconstructions
subjects_dir = data_path + '/subjects'
subject = 'sample'

Computing the forward operator

To compute a forward operator we need:

  • a -trans.fif file that contains the coregistration info.
  • a source space
  • the BEM surfaces

Compute and visualize BEM surfaces

The BEM surfaces are the triangulations of the interfaces between different tissues needed for forward computation. These surfaces are for example the inner skull surface, the outer skull surface and the outer skill surface.

Computing the BEM surfaces requires FreeSurfer and makes use of either of the two following command line tools:

Here we’ll assume it’s already computed. It takes a few minutes per subject.

For EEG we use 3 layers (inner skull, outer skull, and skin) while for MEG 1 layer (inner skull) is enough.

Let’s look at these surfaces. The function mne.viz.plot_bem() assumes that you have the the bem folder of your subject FreeSurfer reconstruction the necessary files.

mne.viz.plot_bem(subject=subject, subjects_dir=subjects_dir,
                 brain_surfaces='white', orientation='coronal')
../_images/sphx_glr_plot_forward_001.png

Out:

Using surface: /home/ubuntu/mne_data/MNE-sample-data/subjects/sample/bem/inner_skull.surf
Using surface: /home/ubuntu/mne_data/MNE-sample-data/subjects/sample/bem/outer_skull.surf
Using surface: /home/ubuntu/mne_data/MNE-sample-data/subjects/sample/bem/outer_skin.surf
Triangle file: Converted tri file nvert = 2562 ntri = 5120
Triangle file: Converted tri file nvert = 2562 ntri = 5120
Triangle file: Converted tri file nvert = 2562 ntri = 5120
Triangle file: created by gramfort on Thu Sep 15 19:44:31 2011 nvert = 155407 ntri = 310810
Triangle file: created by gramfort on Thu Sep 15 22:04:14 2011 nvert = 156866 ntri = 313728

Visualization the coregistration

The coregistration is operation that allows to position the head and the sensors in a common coordinate system. In the MNE software the transformation to align the head and the sensors in stored in a so-called trans file. It is a FIF file that ends with -trans.fif. It can be obtained with mne_analyze (Unix tools), mne.gui.coregistration (in Python) or mrilab if you’re using a Neuromag system.

For the Python version see func:mne.gui.coregistration

Here we assume the coregistration is done, so we just visually check the alignment with the following code.

# The transformation file obtained by coregistration
trans = data_path + '/MEG/sample/sample_audvis_raw-trans.fif'

info = mne.io.read_info(raw_fname)
mne.viz.plot_trans(info, trans, subject=subject, dig=True,
                   meg_sensors=True, subjects_dir=subjects_dir)
../_images/sphx_glr_plot_forward_002.png

Out:

Read a total of 3 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
Using surface from /home/ubuntu/mne_data/MNE-sample-data/subjects/sample/bem/sample-5120-5120-5120-bem.fif
Getting helmet for system 306m

Compute Source Space

The source space defines the position of the candidate source locations. The following code compute such a cortical source space with an OCT-6 resolution.

See Setting up the source space for details on source space definition and spacing parameter.

src = mne.setup_source_space(subject, spacing='oct6',
                             subjects_dir=subjects_dir,
                             add_dist=False, overwrite=True)
print(src)

Out:

Setting up the source space with the following parameters:

SUBJECTS_DIR = /home/ubuntu/mne_data/MNE-sample-data/subjects
Subject      = sample
Surface      = white
Octahedron subdivision grade 6

>>> 1. Creating the source space file /home/ubuntu/mne_data/MNE-sample-data/subjects/sample/bem/sample-oct-6-src.fif...
Doing the octahedral vertex picking...
Loading /home/ubuntu/mne_data/MNE-sample-data/subjects/sample/surf/lh.white...
Triangle file: created by gramfort on Thu Sep 15 19:44:31 2011 nvert = 155407 ntri = 310810
    Triangle neighbors and vertex normals...
Loading geometry from /home/ubuntu/mne_data/MNE-sample-data/subjects/sample/surf/lh.sphere...
Triangle file: created by gramfort on Thu Sep 15 20:42:21 2011 nvert = 155407 ntri = 310810
Mapping lh sample -> oct (6) ...
Setting up the triangulation for the decimated surface...
loaded lh.white 4098/155407 selected to source space (oct = 6)

Loading /home/ubuntu/mne_data/MNE-sample-data/subjects/sample/surf/rh.white...
Triangle file: created by gramfort on Thu Sep 15 22:04:14 2011 nvert = 156866 ntri = 313728
    Triangle neighbors and vertex normals...
Loading geometry from /home/ubuntu/mne_data/MNE-sample-data/subjects/sample/surf/rh.sphere...
Triangle file: created by gramfort on Thu Sep 15 23:55:44 2011 nvert = 156866 ntri = 313728
Mapping rh sample -> oct (6) ...
Setting up the triangulation for the decimated surface...
loaded rh.white 4098/156866 selected to source space (oct = 6)

    Write a source space...
    [done]
    Write a source space...
    [done]
    2 source spaces written
Wrote /home/ubuntu/mne_data/MNE-sample-data/subjects/sample/bem/sample-oct-6-src.fif
You are now one step closer to computing the gain matrix
<SourceSpaces: [<surface (lh), n_vertices=155407, n_used=4098, coordinate_frame=MRI (surface RAS)>, <surface (rh), n_vertices=156866, n_used=4098, coordinate_frame=MRI (surface RAS)>]>

src contains two parts, one for the left hemisphere (4098 locations) and one for the right hemisphere (4098 locations). Sources can be visualized on top of the BEM surfaces.

mne.viz.plot_bem(subject=subject, subjects_dir=subjects_dir,
                 brain_surfaces='white', src=src, orientation='coronal')
../_images/sphx_glr_plot_forward_003.png

Out:

Using surface: /home/ubuntu/mne_data/MNE-sample-data/subjects/sample/bem/inner_skull.surf
Using surface: /home/ubuntu/mne_data/MNE-sample-data/subjects/sample/bem/outer_skull.surf
Using surface: /home/ubuntu/mne_data/MNE-sample-data/subjects/sample/bem/outer_skin.surf
Triangle file: Converted tri file nvert = 2562 ntri = 5120
Triangle file: Converted tri file nvert = 2562 ntri = 5120
Triangle file: Converted tri file nvert = 2562 ntri = 5120
Triangle file: created by gramfort on Thu Sep 15 19:44:31 2011 nvert = 155407 ntri = 310810
Triangle file: created by gramfort on Thu Sep 15 22:04:14 2011 nvert = 156866 ntri = 313728

However, only sources that lie in the plotted MRI slices are shown. Let’s write a few lines of mayavi to see all sources.

import numpy as np  # noqa
from mayavi import mlab  # noqa
from surfer import Brain  # noqa

brain = Brain('sample', 'lh', 'inflated', subjects_dir=subjects_dir)
surf = brain._geo

vertidx = np.where(src[0]['inuse'])[0]

mlab.points3d(surf.x[vertidx], surf.y[vertidx],
              surf.z[vertidx], color=(1, 1, 0), scale_factor=1.5)
../_images/sphx_glr_plot_forward_004.png

Compute forward solution

We can now compute the forward solution. To reduce computation we’ll just compute a single layer BEM (just inner skull) that can then be used for MEG (not EEG).

We specify if we want a one-layer or a three-layer BEM using the conductivity parameter.

The BEM solution requires a BEM model which describes the geometry of the head the conductivities of the different tissues.

conductivity = (0.3,)  # for single layer
# conductivity = (0.3, 0.006, 0.3)  # for three layers
model = mne.make_bem_model(subject='sample', ico=4,
                           conductivity=conductivity,
                           subjects_dir=subjects_dir)
bem = mne.make_bem_solution(model)

Out:

Creating the BEM geometry...
Going from 4th to 4th subdivision of an icosahedron (n_tri: 5120 -> 5120)
inner skull CM is   0.67 -10.01  44.26 mm
Surfaces passed the basic topology checks.
Complete.

Approximation method : Linear collocation

Homogeneous model surface loaded.
Computing the linear collocation solution...
    Matrix coefficients...
        inner_skull (2562) -> inner_skull (2562) ...
    Inverting the coefficient matrix...
Solution ready.
BEM geometry computations complete.

Note that the BEM does not involve any use of the trans file. The BEM only depends on the head geometry and conductivities. It is therefore independent from the MEG data and the head position.

Let’s now compute the forward operator, commonly referred to as the gain or leadfield matrix.

See mne.make_forward_solution() for details on parameters meaning.

fwd = mne.make_forward_solution(raw_fname, trans=trans, src=src, bem=bem,
                                fname=None, meg=True, eeg=False,
                                mindist=5.0, n_jobs=2)
print(fwd)

Out:

Source space                 : <SourceSpaces: [<surface (lh), n_vertices=155407, n_used=4098, coordinate_frame=MRI (surface RAS)>, <surface (rh), n_vertices=156866, n_used=4098, coordinate_frame=MRI (surface RAS)>]>
MRI -> head transform source : /home/ubuntu/mne_data/MNE-sample-data/MEG/sample/sample_audvis_raw-trans.fif
Measurement data             : sample_audvis_raw.fif
BEM model                    : dict
Accurate field computations
Do computations in head coordinates
Free source orientations
Destination for the solution : None

Read 2 source spaces a total of 8196 active source locations

Coordinate transformation: MRI (surface RAS) -> head
     0.999310  0.009985 -0.035787      -3.17 mm
     0.012759  0.812405  0.582954       6.86 mm
     0.034894 -0.583008  0.811716      28.88 mm
     0.000000  0.000000  0.000000       1.00

Read 306 MEG channels from info
72 coil definitions read
Coordinate transformation: MEG device -> head
     0.991420 -0.039936 -0.124467      -6.13 mm
     0.060661  0.984012  0.167456       0.06 mm
     0.115790 -0.173570  0.977991      64.74 mm
     0.000000  0.000000  0.000000       1.00
MEG coil definitions created in head coordinates.
Source spaces are now in head coordinates.

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

Source spaces are in head coordinates.
Checking that the sources are inside the bounding surface and at least    5.0 mm away (will take a few...)
2 source space points omitted because they are outside the inner skull surface.
364 source space points omitted because of the    5.0-mm distance limit.
1 source space point omitted because it is outside the inner skull surface.
331 source space point omitted because of the    5.0-mm distance limit.
Thank you for waiting.

Setting up compensation data...
    No compensation set. Nothing more to do.

Composing the field computation matrix...
Computing MEG at 7498 source locations (free orientations)...

Finished.
<Forward | MEG channels: 306 | EEG channels: 0 | Source space: Surface with 7498 vertices | Source orientation: Free>

We can explore the content of fwd to access the numpy array that contains the gain matrix.

leadfield = fwd['sol']['data']
print("Leadfield size : %d sensors x %d dipoles" % leadfield.shape)

Out:

Leadfield size : 306 sensors x 22494 dipoles

To save to disk a forward solution you can use mne.write_forward_solution() and to read it back from disk mne.read_forward_solution(). Don’t forget that FIF files containing forward solution should end with -fwd.fif.

Exercise

By looking at sphx_glr_auto_examples_forward_plot_read_forward.py plot the sensitivity maps for EEG and compare it with the MEG, can you justify the claims that:

  • MEG is not sensitive to radial sources
  • EEG is more sensitive to deep sources

How will the MEG sensitivity maps and histograms change if you use a free instead if a fixed/surface oriented orientation?

Try this changing the mode parameter in mne.sensitivity_map() accordingly. Why don’t we see any dipoles on the gyri?

Total running time of the script: ( 1 minutes 1.310 seconds)

Generated by Sphinx-Gallery