Set up general paths

Change these to appropriate paths for your operating system and setup

%% paths

addpath /home/lau/matlab/fieldtrip-20170613_workshop_natmeg/
ft_defaults

raw_meg_path = '/archive/20067_workshop_source_reconstruction/MEG/';
meg_path = '/home/share/workshop_source_reconstruction/data/MEG/';
mri_path = '/home/share/workshop_source_reconstruction/data/MRI/';

Set up subject specific paths

Make subject and recording specific paths (the cell array “subjects_and_dates” can be expanded)

%% subjects and dates

subjects_and_dates = ...
                    {
                        'NatMEG_0177/170424/'
                    };

output_path = fullfile(meg_path, subjects_and_dates{1});
cd(output_path)

events = [1 2 4 8 16]; %% right little, right ring, right middle, right index, right thum

Import source space

Minimum-norm estimate (MNE) use a cortical surface as source model. In the tutorial documentation “Prepare source space for MNE” (http://natmeg.se/onewebmedia/prepare_MNE_sourcespace.html) you can read how to export MRI to Freesurfer to extract the cortical surface and setup a source model with MNE-C. This creates a set of point equally sampled across the cortical surface that is used for source modeling.

You will find the file “Sub02-oct-6-src.fif” in the tutorial datafiles.

We read the MNE-C source space into MATLAB by using ft_read_headshape.

sourcemodel = ft_read_headshape('Sub02-oct-6-src.fif', 'format', 'mne_source');

Note that the units are in meters. You can use ft_convert_units to make the units milimeters instead:

sourcemodel = ft_convert_units(sourcemodel, 'mm');

You can plot the source model with the function ft_plot_mesh:

figure; ft_plot_mesh(sourcemodel, 'edgecolor', 'k'); camlight

Looks nice!

Each point on the grid will represent a dipole in the MNE source reconstruction.

Load data

load headmodel_mne.mat
load timelockeds.mat

Now you should have the three main ingredients for doing MNE source reconstruction with MEG:

  1. Processed evoked MEG data.
  2. A volume model of the brain.
  3. A source space model of the cortical surface.

Take a look at the headmodel and sourcemodel toghether:

figure; hold on
ft_plot_vol(headmodel_mne, 'facealpha', 0.5, 'edgecolor', 'none');
ft_plot_mesh(sourcemodel, 'edgecolor', 'k'); camlight
  • What is going on?

Aligning source space and volume model

We have aligned the MRI that we used to create the volume model to the Neuromag coordinate system used by the MEG scanner. However, in order to create the cortical surface in Freesurfer, we had to convert the MRI to the spm coordinate system and exported it to Freesurfer (see http://natmeg.se/onewebmedia/prepare_MNE_sourcespace.html). When we imported the source space back into MATLAB, it kept this coordinate system. Our head model and source space are now in two different coordinate systems.

To get the cortical surface back into the head, you need to transform the points in source space to the coordinate system of the head model. Load the transformation matrix of the resliced MRI in the spm coordinate system. We will use this to create a transformation from the spm coordinate system to the neuromag coordinate system. If you have run the prepare MNE source space tutorial, you must read the transformation files you created yourself.

% from voxel to spm-coordsys.
load transform_vox2spm.mat

% from voxel to neuromag-coordsys.
transform_vox2neuromag = mri_segmented_mne.transform;

% Get transformation matrix from spm to neuromag
T = transform_vox2neuromag/transform_vox2spm_rs;

Now transform the source space using the transformation matrix T. Note, that if you run the following line several times, it will apply the transformation each time, transforming from the current position to a new position, which will be wrong. If in doubt, save the source model before proceeding.

sourcemodel = ft_transform_geometry(T, sourcemodel);

Take a look at the headmodel and sourcespace again:

figure; hold on
ft_plot_vol(headmodel_mne, 'facealpha', 0.5, 'edgecolor', 'none');
ft_plot_mesh(sourcemodel, 'edgecolor', 'k'); camlight

Remember to save the source model.

save('sourcemodel.mat','sourcemodel')
disp('Done');

MNE source reconstruction on MEG data

Now we can create the leadfield and do the source reconstruction. Load the relevant files (if you do not already have then loaded):

load timelockeds.mat
load sourcemodel.mat
load headmodel_mne.mat

For this tutorial we will use the stimulation of the indexfinger – corresponding the the 4th condition:

data_meg = timelockeds{4};

Calculate the leadfield using ft_prepare_leadfield with the MEG data, the sourcemodel, and the appropriate head model. Here we use the gradiometers (meggrad).


cfg = [];
cfg.grad                = data_meg.grad;              % sensor positions
cfg.channel             = 'meggrad';                  % the used channels
cfg.senstype            = 'meg';
cfg.grid.pos            = sourcemodel.pos;            % source points
cfg.grid.inside         = 1:size(sourcemodel.pos,1);  % all source points are inside of the brain
cfg.headmodel           = headmodel_mne;          % volume conduction model

leadfield_mne = ft_prepare_leadfield(cfg,data_meg);
plot3(leadfield_mne.pos(:,1),leadfield_mne.pos(:,2),leadfield_mne.pos(:,3),'0')
plot3(leadfield.pos(:,1),leadfield.pos(:,2),leadfield.pos(:,3),'0')

Finally, we do the source reconstruction using ft_sourceanalysis. Specify cfg.method = ‘mne’ to use MNE. MNE also require additional parameters indicating how to scale the noise covariance. Note that the noise covariance matrix already is in the “data” structure. This was calculated already in the preprocessing steps when the evoked fields were calculated.

cfg                     = [];
cfg.method              = 'mne';
cfg.channel             = 'meggrad';
cfg.senstype            = 'meg';
cfg.grid                = leadfield_mne;
cfg.headmodel           = headmodel_mne;
cfg.mne.prewhiten       = 'yes';
cfg.mne.lambda          = 3;
cfg.mne.scalesourcecov  = 'yes';
source_mne  = ft_sourceanalysis(cfg,data_meg);

save source_mne.mat source_mne;
disp('Done')

Visualize

Plot the MNE source reconstruction on the grid. The source_mne structure contains values for all grid points and all timepoints. To visualize the source reconstruction on the grid, we need to decide what time points to plot.

Go back to the evoked data to find some interesting times to plot using ft_multiplotER. Since we used the gradiometers, we will start by plotting the fields from the combined gradiometers. For example, let us compare the first peak around 50 ms and the larger peak from 140-150 ms. Use the drag and click tool on the multiplot to plot topographical plots.

data_cmb = ft_combineplanar([],data_meg); %Combine gradiometers

cfg = [];
cfg.layout = 'neuromag306cmb.lay';
figure; ft_multiplotER(cfg, data_cmb);

Select the time-windows of interest and plot on the sourcespace:

source_mne.tri = sourcemodel.tri; %ft_sourceplot need the triangulation information, add to source reconstruction.

cfg = [];
cfg.method          = 'surface';
cfg.funparameter    = 'pow';
cfg.funcolormap     = 'jet';
cfg.latency         = .060;
cfg.colorbar        = 'no';

ft_sourceplot(cfg, source_mne)

This image shows a single time point for all estimated sources. Each source (i.e. each grid-point on the cortical surface) has its own time series. You can in principle treat each source as its own “channel”. E.g. you can view the activation of all sources over time, analogous to viewing the activation in MEG or EEG channels. Use MATLAB to plot source activation across time (here only plotting every third time-series to save memory)

plot(source_mne.time,source_mne.avg.pow(1:3:end,:));

Now try to use ft_sourceplot to plot the source space topography of other time points.

You can visualize the average over time over with ft_sourceplot by changing the cfg.latency to an interval [start stop] and giving an additional cfg argument cfg.avgovertime = ‘yes’.

cfg = [];
cfg.method          = 'surface';
cfg.funparameter    = 'pow';
cfg.funcolormap     = 'jet';
cfg.latency         = [.350 .40]; %350 ms to 400 ms
cfg.avgovertime     = 'yes';
cfg.colorbar        = 'no';

ft_sourceplot(cfg, source_mne)

To get an easy overview of both the temporal and spatial features of the source reconstruction, you can use ft_sourcemovie to make a movie of the source activation. Note that the triangulation field must be present in the source_mne structure for ft_sourcemovie to work. You can play around with the scaling to get a nicer looking movie.

cfg = [];
cfg.funparameter = 'pow';
ft_sourcemovie(cfg,source_mne);