Head Model User Guide

What is a Head Model and what is it used for?

In order to estimate how neural currents produce electro-cortical potential differences and magnetic fields at external sensors used in neuroimaging recordings like electroencephalography (EEG) or magnetoencephalography (MEG), head models are created by computing a forward solution known as a "gain" or "leadfield" matrix to explain this transformation while accounting for the brain anatomy and different head tissues (white/grey matter, cerebrospinal fluid (CSF), skull, bone, skin, etc.).

Once a head model is created then various source localization methods (aka "inverse solutions") can be used to estimate the current density for each source point (vertex) specified in the anatomy of the head model. These individual sources can then be grouped together (anatomically or functionally) and segmented (parcellated) into specific regions of interest (ROIs) such that the segmented sources and label names are stored together in a dictionary called a "brain atlas". These atlases can be custom made to study specific regions of interest or there are standardized, pre-built atlases that use common ROIs for consistent comparisons across different studies (e.g. Desikan–Killiany, Destrieux) (see https://surfer.nmr.mgh.harvard.edu/fswiki/CorticalParcellation).

Head models are computed using one of many different forward solution methods (e.g. boundary element method (BEM), finite element method (FEM), etc.) that are implemented in a number of available software applications (e.g. Brainstorm, MNE, EEGLAB headModel toolbox, OpenMEEG, etc.) available to do so. Because each software package has its own methods for creating head models and these methods are already well documented, this guide is not meant to explain the head model creation process. NeuroPype currently does not have the ability to create a head model. Instead, this guide will assume you have already created your own head model from an external application and show you how to properly format it and prepare it to be used within NeuroPype.

What is contained in a Head Model?

A head model is essentially the result of calculating how the source space (brain areas) and external sensors (EEG channels or MEG sensors) are related. A basic head model should contain a channel montage, which is information about the EEG/MEG channels/sensors including the label names and coordinates of their positions. It will also contain the anatomical information (meshes) related to the brain model (e.g. cortex surface, CSF, skull, scalp). And finally, it will contain the computed model that is a transformation matrix of brain sources to external sensors called the gain/leadfield matrix (known as the "forward model"). Optionally, the head model can store or link to one or many brain atlases that provide segmented groups of sources and labels for each group. These pieces of information are stored and named in various ways depending on the software used to create the model. It doesn't matter which software is used to create the head model as long as the data can be exported and stored in the specific structure and format that NeuroPype expects.

NeuroPype Head Model Format

In NeuroPype, head models must be stored as MATLAB v7.3 files. Typically, we use Brainstorm (which is based in MATLAB) to create our head models so we save the head models out with a ".hm" extension to signify "head model", but in reality these files are exactly the same as MATLAB ".mat" files. Because of this, NeuroPype .hm files can be easily imported into MATLAB using the simple "load" function. Example:

HM = load('<path_to_cpe>/resources/headmodels/Colin-339ch-4495v-scalar-4shell-1.1.hm', '-mat');

A complete head model is a struct variable (in this example 'HM') that contains 6 fields of data (with multiple subfields) named precisely (fieldname spelling and casing is important) as follows (field order is not important):

1. HM.leadfield

HM.leadfield.matrix: double [Nsensors x Nsources]

This is the gain/leadfield matrix (aka the "forward model") and the most important piece of information. Typically, this is computed as an unconstrained model, which are 3 orthogonal orientations (x,y,z) at each grid point (p1, p2, etc). The information relative to each pair sensor <-> grid source point is stored as successive columns of the matrix and ordered as: [p1_x, p1_y, p1_z, p2_x, p2_y, p2_z ...]. However, NeuroPype uses the "orientation-constrained" version of this model where the orientation of each dipole is fixed and normal to the cortex surface. Brainstorm has a function to automatically do this conversion (copied from https://neuroimage.usc.edu/brainstorm/Tutorials/HeadModel):

To convert this unconstrained leadfield matrix to that of an orientation-constrained model, where the orientation of each dipole is fixed and normal to the cortex surface:

1. Export the head model file to the "HeadModel" structure: Right-click > File > Export to Matlab. Name the workspace variable "HeadModel".

2. At the Matlab prompt:

      `Gain_constrained = bst_gain_orient(HeadModel.Gain, HeadModel.GridOrient);`

  The dimension of the output matrix size is three times smaller (now only one source orientation at each location): [Nsensors x Nsources]

2. HM.sensors * HM.sensors.labels: cell array [1 x Nsensors]

Each cell stores the channel/sensor labels (as string)

HM.sensors.coordinates: double [Nsensors x 3]

Data matrix of channels/sensors by X, Y, and Z coordinates (in millimeters). X should correspond to "right" direction, Y to "front" direction, and Z to "up" direction.

3. HM.laplacian

HM.laplacian.matrix: double [Nsources by Nsources]

Because the cortex mesh might be irregular, a laplacian interpolation should be calculated (ref: Oostendorp, Oosterom & Huiskamp (1989), Interpolation on a triangulated 3D surface. Journal of Computational Physics, 80: 331-343). This can be obtained using an open source Matlab function here: https://www.mathworks.com/matlabcentral/fileexchange/1875-mesh-laplacian

4. HM.meshes: 1 instance per tissue (e.g. Brain, CSF, Skull, Scalp)

HM.meshes.vertices: double [Nvertices x 3]

Data matrix of mesh vertices by X, Y, and Z coordinates (in millimeters).

HM.meshes.faces: double [Nfaces x 3]

Data matrix of triangles constituting the surface mesh. Each value is a pointer to an indexed vertex number (int).

NOTE: This is an index value and different softwares use different indexing systems (e.g. 0 or 1 base). Because our head models are created in Brainstorm/MATLAB, it is recommended to use MATLAB's 1-based indexing system. For instance, in MATLAB the min(HM.atlases.labeling) should be 1. In Python, numpy.min(HM.atlases.labeling) would be 0. NeuroPype by default assumes a 1-based indexing was used in the head model and converts this field to 0-based indexing after importing for downstream processing since NueroPype is based in Python. Because of this, it is best to just leave index values in 1-base.

HM.meshes.name: string

Name of each mesh. Currently, NeuroPype only accepts "Brain", "CSF". "Skull", and "Scalp".

5. HM.atlases: 1 instance per atlas (e.g. Desikan–Killiany, Destrieux). This can store an unlimited number of atlases.

HM.atlases.labeling: double [Nsources]

This stores a single value (int) for each source/vertex that is an index corresponding to which ROI label the source falls in.

NOTE: This is an index value and different softwares use different indexing systems (e.g. 0 or 1 base). Because our head models are created in Brainstorm/MATLAB, it is recommended to use MATLAB's 1-based indexing system. For instance, in MATLAB the min(HM.atlases.labeling) should be 1. In Python, numpy.min(HM.atlases.labeling) would be 0. NeuroPype by default assumes a 1-based indexing was used in the head model and converts this field to 0-based indexing after importing for downstream processing since NeuroPype is based in Python. Because of this, it is best to just leave index values in 1-base.
HM.atlases.labels: cell array [Nrois]

Each cell stores the region of interest (ROI) labels (as string)

HM.atlases.name: string

Name of the atlas

6. HM.meta

This stores the meta information regarding the coordinate system being used. For compatibility, we recommend using the following system:

X should correspond to "right" direction, Y to "front" direction, and Z to "up" direction. All values are in millimeters and in the MNI coordinate system. This is currently the default system that NeuroPype supports.

HM.meta.coordinates: struct
- HM.meta.coordinates.unit: string
Coordinate unit type (e.g. "millimeters")
- HM.meta.coordinates.x: string
Direction of positive x-axis ['left', 'right', 'front', 'up']
- HM.meta.coordinates.y: string
Direction of positive y-axis ['left', 'right', 'front', 'up']
- HM.meta.coordinates.z: string
Direction of positive z-axis ['left', 'right', 'front', 'up']

Tutorial Example for Extracting a Head Model from Brainstorm

We have chosen to use Brainstorm for calculating our head models supplied in NeuroPype. We have based these head models on the Colin 27 Average Brain (see http://www.bic.mni.mcgill.ca/ServicesAtlases/Colin27). We have also used a custom EEG channel montage of 339 channels across many different EEG caps and naming systems. The following example shows how to export all the necessary data from the Brainstorm GUI after the head model has been calculated.

1. Extract and constrain the Gain matrix

Once you have successfully computed the Gain Matrix in Brainstorm, if you right click on it and goto "File" -> "View file contents" you can inspect all the variables that it stores (shown below):

image alt text

Note the "GridLoc" variable indicates there are 4496 vertices (sources). However, the size of the "Gain" variable is showing number of channels (339) by three times the number of vertices (4496 * 3 = 13488). Neuropype expects the "constrained" version of the gain matrix which should be [Nsensors x Nsources]. To do obtain this, first export the Gain Matrix structure to Matlab by right clicking on "File" -> "Export to Matlab" and then give your variable a name like "Gain_Matrix". This variable will now appear in your Matlab workspace. In the Matlab command window, type:

Gain = bst_gain_orient(Gain_Matrix.Gain, Gain_Matrix.GridOrient);

You will now have a new variable, "Gain", that is the correct number of channels (339) by number of vertices (4496). Note the gain matrix units from Brainstorm. "All internal units in Brainstorm, with possibly minor exceptions in a few cases for user clarity, are in MKS, so the forward model is designed to transform current dipoles expressed as Amps-meter (A-m) into Volts, and therefore the units of the gain matrix are (for EEG): V / (A-m)." (source: https://neuroimage.usc.edu/forums/t/eeg-units/1499). Because NeuroPype expects units in millimeters, this should be accounted for by dividing by 1000.

Place this matrix into the final head model struct variable, HM:

HM.leadfield.matrix = Gain / 1000;

2. Extract sensor information

Similarly to the Gain matrix, first export the channel/sensor object by right clicking on the channel montage, select "File" -> "Export to Matlab" and give it a variable name like "EGI". First add the channel label names to the HM.sensors.labels field:

HM.sensors.labels = {EGI.Channel.Name};

Next, add the channel coordinates:

HM.sensors.coordinates = [EGI.Channel.Loc]';

Note: NeuroPype reads from the "meta.coordinates" field in the head model to account for units and orientations. It is important to keep track of the format of the sensors based on the channel montage you used to compute your head model. In this step, it may be necessary to re-arrange x, y, and z coordinate columns, apply a sign change, and/or scale values in order to be consistent with the meta information stored. We recommend using Neuropype's default expectations: X should correspond to "right" direction, Y to "front" direction, and Z to "up" direction. All values are in millimeters and in the MNI coordinate system.

For instance, if your montage used meters for location values, you should multiply by 1000 to obtain millimeters. If your X-axis is pointing to the left direction (positive values point left), you should multiply by -1 so X-axis is now pointing "right".

3. Extract meshes

Here you will need to display the anatomical tab that shows each mesh layer used in your head model. Export each one separately to MATLAB using the same procedure as above and copy the "Vertices" and "Faces" fields to HM.meshes struct.

image alt text

For brain mesh (exported to MATLAB as "brain" variable):

HM.meshes(1).name = 'Brain';

HM.meshes(1).vertices = brain.Vertices;

HM.meshes(1).faces = brain.Faces;

Repeat this process for each of the anatomical layers that were used for computing the gain matrix.

4. Compute Laplacian interpolation of cortex and store in head model

Once you have exported/stored the "brain" or "cortex" mesh, you will need to run the MATLAB function "mesh_laplacian" (link to download this function provided in previous section).

[lap, edge] = mesh_laplacian(HM.meshes(1).vertices, HM.meshes(1).faces);

HM.laplacian.matrix = lap;

5. Store the brain atlases

For NeuroPype head models you only need to store the label names and the index value of the corresponding label name of each ROI (called "Scouts" in Brainstorm) that each source vertex resides in. NeuroPype head models can also store unlimited number of atlases in the head model (create a new instance for each).

Once you have imported a given atlas into Brainstorm, export it to a MATLAB variable (e.g. "atlas") the same way as the other variables. You can use the following template script to correctly merge it into the HM struct head model:

HM.atlases = struct('name','','labels','','labeling','');
HM.atlases.name = atlas.Atlas(1).Name;
HM.atlases.labels = {atlas.Atlas(1).Scouts.Label};
HM.atlases.labeling = zeros(size(atlas.Vertices,1),1);
for iiScout = 1:numel(atlas.Atlas(1).Scouts)
  HM.atlases.labeling(atlas.Atlas(1).Scouts(iiScout).Vertices) = iiScout;
end

6. Store the meta data

HM.meta.coordinates.unit = 'millimeters';

HM.meta.coordinates.x = 'right';

HM.meta.coordinates.y = 'front';

HM.meta.coordinates.z = 'up';

7. Save the head model

save(<OutputFile>,'-struct','HM','-mat','-v7.3');

Your head model is now compatible with NeuroPype. When you use a NeuroPype node that requires a head model (e.g. sLORETA), simply provide the filepath to the head model file you created. If there are any issues with your head model, NeuroPype will provide error messages explaining what might be wrong. If you can't solve your issues, note the error messages and send them along with your head model file to support for review.