Tags: tutorial timelock source meg headmodel mri plotting meg-language

Source reconstruction of event-related fields using minimum-norm estimation


In this tutorial you can find information about how to do source reconstruction using minimum-norm estimation, to reconstruct the event-related fields (MEG) of a single subject. We will be working on the dataset described in the Preprocessing - Segmenting and reading trial-based EEG and MEG data and Event-related averaging and MEG planar gradient tutorials, and we will use also the anatomical images that belong to the same subject. We will repeat code to select the trials and preprocess the data as described in the Event-related averaging and MEG planar gradient tutorial. We assume that preprocessing and event-related averaging is already clear for the reader. To preprocess the anatomical data, we will use two other software packages (FreeSurfer and MNE Suite).

This tutorial will not show how to do group-averaging and statistics on the source-level. It will also not describe how to do source-localization of oscillatory activation. You can check the Localizing oscillatory sources using beamformer techniques tutorial if you are interested in this latest.


In the Event-related averaging and MEG planar gradient tutorial time-locked averages of event-related fields of three conditions have been computed and the Cluster-based permutation tests on event-related fields tutorial showed that there was a significant difference among two conditions. The topographical distribution of the ERFs belonging to each conditions and ERFs belonging to those differences have been plotted. The aim of this tutorial is to calculate a distributed representation of the underlying neuronal activity that resulted in the brain activity observed at the sensor level.

To calculate distributed neuronal activation we will use the minimum-norm estimation. This approach is favored for analyzing evoked responses and for tracking the wide-spread activation over time. It is a distributed inverse solution that discretizes the source space into locations on the cortical surface or in the brain volume using a large number of equivalent current dipoles. It estimates the amplitude of all modeled source locations simultaneously and recovers a source distribution with minimum overall energy that produces data consistent with the measurement (Ou, W., Hämäläinen, M., Golland, P., 2008, A Distributed Spatio-temporal EEG/MEG Inverse Solver. Jensen, O., Hesse, C., 2010, Estimating distributed representation of evoked responses and oscillatory brain activity, In: MEG: An Introduction to Methods, ed. by Hansen, P., Kringelbach, M., Salmelin, R., doi:10.1093/acprof:oso/9780195307238.001.0001). The reference for the implemented method is Dale et al. (2000).


Figure 1 shows a schematic of the steps needed for the calculation of the minimum-norm estimate. It shows that the computation of the inverse solution is based on the outputs of two independent processing steps: the processing of the anatomical images that leads to a forward model and the processing of the MEG data. To create a useable source model, additional software is needed, for example FreeSurfer (for the creation of a model of the cortical sheet), and MNE Suite or HCP workbench (to get a minimally distorted low-resolution version of the cortical sheet).

Figure 1. A schematic overview of the steps needed for the calculation of the minimum-norm estimate

The forward model requires three geometric object

  • A volume conduction model of the head, also known as headmodel.
  • A sourcemodel, we advocate a minimally distorted low-resolution description of the cortical sheet.
  • A geometric description of the sensor-array (electrode positions + referencing information for EEG, coil positions/orientation and balancing information for MEG).

The sourcemodel and headmodel are ideally generated from a subject-specific MRI image. The description of the sensor-array typically is represented in the data (MEG), or needs to be constructed, for example with a Polhemus device (EEG). The construction of the head- and sourcemodels that are needed for the remainder of this tutorial is described in the following tutorial

Once we have the headmodel and sourcemodel, we perform the following step

Processing of functional data

The following will use the MEG data belonging to Subject01 which can be obtained from our download server. For both preprocessing and averaging, we will follow the steps that have been written in the Event-related averaging and MEG planar gradient tutorial. We will use trials belonging to two conditions (FC and FIC) and we will calculate their difference.

Preprocessing of MEG data

We will now read and preprocess the data. If you would like to continue directly with the already preprocessed data, you can download dataFIC_LP.mat and dataFC_LP.mat. Load the data into MATLAB with the command load and skip to Averaging and noise-covariance estimation.

Otherwise run the following code:

The ft_definetrial and ft_preprocessing functions require the original MEG dataset, which is available here

Do the trial definition for the all conditions together:

cfg                         = [];
cfg.dataset                 = 'Subject01.ds';
cfg.trialfun                = 'ft_trialfun_general'; % this is the default
cfg.trialdef.eventtype      = 'backpanel trigger';
cfg.trialdef.eventvalue     = [3 5 9]; % the values of the stimulus trigger for the three conditions
% 3 = fully incongruent (FIC), 5 = initially congruent (IC), 9 = fully congruent (FC)
cfg.trialdef.prestim        = 1; % in seconds
cfg.trialdef.poststim       = 2; % in seconds

cfg = ft_definetrial(cfg);


% remove the trials that have artifacts from the trl
cfg.trl([2, 5, 6, 8, 9, 10, 12, 39, 43, 46, 49, 52, 58, 84, 102, 107, 114, 115, 116, 119, 121, 123, 126, 127, 128, 133, 137, 143, 144, 147, 149, 158, 181, 229, 230, 233, 241, 243, 245, 250, 254, 260],:) = [];

% preprocess the data
cfg.channel    = {'MEG', '-MLP31', '-MLO12'};        % read all MEG channels except MLP31 and MLO12
cfg.demean     = 'yes';
cfg.baselinewindow  = [-0.2 0];
cfg.lpfilter   = 'yes';                              % apply lowpass filter
cfg.lpfreq     = 35;                                 % lowpass at 35 Hz.

data_all = ft_preprocessing(cfg);

These data have been cleaned from artifacts by removing several trials and two sensors; see the visual artifact rejection tutorial.

cfg = [];
cfg.trials = data_all.trialinfo == 3;
dataFIC_LP = ft_redefinetrial(cfg, data_all);

cfg = [];
cfg.trials = data_all.trialinfo == 9;
dataFC_LP = ft_redefinetrial(cfg, data_all);

Subsequently you can save the data to disk.

  save dataFIC_LP dataFIC_LP
  save dataFC_LP dataFC_LP

Averaging and noise-covariance estimation

The function ft_timelockanalysis makes averages of all the trials in a data structure and also estimates the noise-covariance. For a correct noise-covariance estimation it is important that you used the cfg.demean = ‘yes’ option when the function ft_preprocessing was applied.

The trials belonging to one condition will now be averaged with the onset of the stimulus time aligned to the zero-time point (the onset of the last word in the sentence). This is done with the function ft_timelockanalysis. The input to this procedure is the dataFC_LP structure generated by ft_preprocessing. At the same time, we need to compute the noise-covariance matrix, therefore cfg.covariance = ‘yes’ has to be specified as well as the time window where the noise-covariance will be estimated. Here, we use the baseline where there is no signal of interest yet.

  load dataFC_LP
  load dataFIC_LP

  cfg = [];
  cfg.covariance = 'yes';
  cfg.covariancewindow = [-inf 0]; %it will calculate the covariance matrix
                                   % on the timepoints that are
                                   % before the zero-time point in the trials
  tlckFC = ft_timelockanalysis(cfg, dataFC_LP);
  tlckFIC = ft_timelockanalysis(cfg, dataFIC_LP);
  save tlck tlckFC tlckFIC;

Forward solution

The source space, the volume conduction model and the position of the sensors are necessary inputs for creating the leadfield (forward solution) with the ft_prepare_leadfield function. The sensor positions are contained in the grad field of the averaged data. However, the grad field contains the positions of all channels, therefore, the used channels have to be also specified. The source model file can be obtained here and the head model file here.

load tlck
load Subject01_sourcemodel_15684
load Subject01_headmodel

cfg         = [];
cfg.grad    = tlckFC.grad;   % sensor information
cfg.channel = tlckFC.label;  % the used channels
cfg.grid    = sourcemodel;   % source points
cfg.headmodel = headmodel;   % volume conduction model
cfg.singleshell.batchsize = 5000; % speeds up the computation
leadfield   = ft_prepare_leadfield(cfg);

save Subject01_leadfield leadfield;

Inverse solution

The ft_sourceanalysis function calculates the inverse solution. The method used (minimum-norm estimation) has to be specified with the cfg.method option. The averaged functional data, the forward solution (the output of the ft_prepare_leadfield function), the volume conduction model (in this case, the output of the ft_prepare_headmodel function) and the noise-covariance matrix (the cov field of the output of the ft_timelockanalysis function) have to be provided.

The lambda value is a scaling factor that is responsible for scaling the noise-covariance matrix. If it is zero the noise-covariance estimation will be not taken into account during the computation of the inverse solution. Noise-covariance is estimated in each trial separately and then averaged, while the functional data (of which we calculate the source-analysis) is simply averaged across all the trials. Therefore, the higher the number of trials the lower the noise is in the averaged, functional data, but the number trials is not reducing the noise in the noise-covariance estimation. This is the reason while it is useful to use a scaling factor for the noise-covariance matrix if we want to estimate more realistically the amount of noise.

You do not have to specify of the noise-covariance matrix separately, because it is in the tlckFC.cov and in the tlckFIC.cov fields, and ft_sourceanalysis will take it into account automatically.

load tlck;
load Subject01_leadfield;
load Subject01_headmodel;

cfg               = [];
cfg.method        = 'mne';
cfg.grid          = leadfield;
cfg.headmodel     = headmodel;
cfg.mne.prewhiten = 'yes';
cfg.mne.lambda    = 3;
cfg.mne.scalesourcecov = 'yes';
sourceFC          = ft_sourceanalysis(cfg,tlckFC);
sourceFIC         = ft_sourceanalysis(cfg, tlckFIC);

save source sourceFC sourceFIC;


You can plot the inverse solution onto the source-space at a specific time-point with the ft_plot_mesh function.

load source;

m=sourceFIC.avg.pow(:,450); % plotting the result at the 450th time-point that is
                         % 500 ms after the zero time-point
ft_plot_mesh(sourceFIC, 'vertexcolor', m);
view([180 0]); h = light; set(h, 'position', [0 1 0.2]); lighting gouraud; material dull

Figure 6. The result of the source reconstruction of the FIC condition plotted onto the source-space at 500 ms after the 0 time-point

But we would like to know where the difference between the conditions can be localized. Therefore, we calculate the difference of the two conditions, and we use ft_sourcemovie to visualize the results.

cfg = [];
cfg.projectmom = 'yes';
sdFC  = ft_sourcedescriptives(cfg,sourceFC);
sdFIC = ft_sourcedescriptives(cfg, sourceFIC);

sdDIFF         = sdFIC;
sdDIFF.avg.pow = sdFIC.avg.pow - sdFC.avg.pow;

save sd sdDIFF;

cfg = [];
cfg.funparameter = 'pow';

Figure 7. One frame from the movie that shows the differences of the two source reconstructions

Summary and further readings

In this tutorial we showed how to do MNE source reconstruction method on a single subject data. We compared the averaged ERF in two conditions and we reconstructed the sources and we calculated the difference of the two source reconstruction. We showed also how you can visualize the results.

Functions and tutorial pages that show how to average, and how to analyze statistically source reconstructions across subjects or how to compare those to a template brain are still under development.


Example script