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

## Introduction

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.

## Background

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).

## Procedure

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

- Creating a volume conduction model of the head for source reconstruction of MEG data
- Creating a sourcemodel for source reconstruction of MEG or EEG data

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

- compute the forward solution using
**ft_prepare_leadfield**; - preprocess the MEG data using
**ft_definetrial**and**ft_preprocessing**; - compute the average over trials and estimate the noise-covariance using
**ft_timelockanalysis**; - compute the inverse solution using
**ft_sourceanalysis**and**ft_sourcedescriptives**; - visualize the results with
**ft_plot_mesh**and**ft_sourcemovie**.

## 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);
```

### Cleaning

```
% 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;
```

## Visualization

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';
ft_sourcemovie(cfg,sdDIFF);
```

*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.

FAQ

- Where is the anterior commissure?
- Can I do combined EEG and MEG source reconstruction?
- Can I restrict the source reconstruction to the grey matter?
- How are the different head and MRI coordinate systems defined?
- How are the Left and Right Pre-Auricular (LPA and RPA) points defined?
- How can I check whether the grid that I have is aligned to the segmented volume and to the sensor gradiometer?
- How can I determine the anatomical label of a source or electrode?
- How can I fine-tune my BEM volume conduction model?
- How can I visualize the different geometrical objects that are needed for forward and inverse computations?
- How to coregister an anatomical MRI with the gradiometer or electrode positions?
- Is it good or bad to have dipole locations outside of the brain for which the source reconstruction is computed?
- Is it important to have accurate measurements of electrode locations for EEG source reconstruction?
- What is the conductivity of the brain, CSF, skull and skin tissue?
- What kind of volume conduction models of the head (head models) are implemented?
- Where can I find the dipoli command-line executable?
- Why is there a rim around the brain for which the source reconstruction is not computed?
- Why should I use an average reference for EEG source reconstruction?

Example script

- Combined EEG and MEG source reconstruction
- Common filters in beamforming
- Compute forward simulated data using ft_dipolesimulation
- Compute forward simulated data and apply a beamformer scan
- Compute forward simulated data and apply a dipole fit
- Compute forward simulated data with the low-level ft_compute_leadfield
- Check the quality of the anatomical coregistration
- Localizing the sources underlying the difference in event-related fields
- Fit a dipole to the tactile ERF after mechanical stimulation
- Align EEG electrode positions to BEM headmodel
- How to create a head model if you do not have an individual MRI
- How to import data from MNE-Python and FreeSurfer
- Make leadfields using different headmodels
- Plotting the result of source reconstruction on a cortical mesh
- Source statistics
- Create MNI-aligned grids in individual head-space
- Symmetric dipole pairs for beamforming
- Testing BEM created EEG lead fields
- Use your own forward leadfield model in an inverse beamformer computation