In this tutorial, we will do source reconstruction of time-frequency data with a beamformer method known as Dynamic Imaging of Coherent Sources (DICS; Gross et al. 2001). DICS is based on reconstructing sources that show strong dependency (coherence) in the frequency domain.

We will start by finding a time-frequency area of interest that want to know the underlying sources. From time-domain data, we will calculate the cross-spectral density of the frequency (or frequencies) of interest and then use the cross-spectral density to find underlying the sources. Most of this is done “under the hood” by FieldTrip with the function ft_sourceanalysis.

First, we need to prepare the raw data, define a source model, and calculate the lead field.

Set up general paths

Note that if work from the path where all the downloadable parts are downloaded to, you don’t need to change the paths(leave them as is, but do evaluate the sections). The paths serve as an example of how you can set up your analysis structure, which is especially useful if you have more than one subject

Change these to appropriate paths for your operating system and setup. Note that if work from the path where all the downloadable parts are downloaded to, you don’t need to change the paths (leave them as is, but do evaluate the section).

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

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/';

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

%% subjects and dates

subjects_and_dates = ...

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

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

Load files necessary for the beamforming

We are loading four different files here:

  1. We are loading the TFRs to locate interesting targets for beamformer source reconstruction
  2. We are loading the head model for the MEG data
  3. We are loading the head model for the EEG data
  4. We are loading the preprocessed data (baseline_data) on which we will do the actual source reconstruction

If you saved the files from the preprocessing tutorial, head model tutorial, and TFR analysis tutorial you should load those files. If not, you can find the following files in the tutorial material folder:

%% go to relevant path and load data

disp 'Loading input data'
load combined_tfrs.mat
load headmodel_meg.mat
load headmodel_eeg.mat
load baseline_data.mat
load cleaned_downsampled_data.mat
disp Done

Set channels

Here you can set which channels you want to base the source reconstruction on. For now, we will use the gradiometers, but you can always go back and try to re-run the tutoral with the other sensor types. You only need to change to one of the three options specified

%% set channels

channels = 'meggrad';  % either 'megmag', 'meggrad' or 'eeg'

if strcmp(channels, 'eeg')
    headmodel = headmodel_eeg;
    sensor_type = 'eeg';
    headmodel = headmodel_meg;
    sensor_type = 'meg';

Identify interesting features

Plot all the channels in the data. Notice (at least) three features of interest. Please explore the plots!

  1. The so-called beta rebound from about 700 to 1000 msec at around 16 Hz
  2. The so-called mu desynchronization from about 300 to 800 msec at around 10 Hz
  3. The so-called high-beta desynchronization from about 200 to 500 msec at around 22 Hz
%% identify interesting features in the data

% pick an event
event_index = 1;

cfg = [];
cfg.layout = 'neuromag306cmb.lay'; %% this is combined gradiometers, you can also choose <'neuromag306mag.lay'> or <'neuromag306eeg1005_natmeg.lay'
cfg.baselinetype = 'relative'; %% absolute power is not terribly meaningful; therefore we use 'relative' to look at increases and decreases in power relative to the overall power level
cfg.baseline = [-Inf Inf]; %% from min to max
cfg.colorbar = 'yes'; %% show the interpretation of the colours
% cfg.zlim = [0.6 1.6]; %% play around with this parameter to familiarize yourself with these plots

ft_multiplotTFR(cfg, combined_tfrs{event_index});

Visualizing channels

%% channels that we'll plot throughout

colours = ones(306, 3);
tactile_channel = 'MEG0432';
occipital_channel = 'MEG2142';
tactile_channel_index = find(strcmp(baseline_data.label, tactile_channel));
occipital_channel_index = find(strcmp(baseline_data.label, occipital_channel));

colours(tactile_channel_index, :)   = [1 0 0]; %% make red
colours(tactile_channel_index + 1, :)   = [1 0 0]; %% make red
colours(occipital_channel_index, :) = [0 0 1]; %% make blue
colours(occipital_channel_index + 1, :) = [0 0 1]; %% make blue

figure('units', 'normalized', 'outerposition', [0 0 1 1]);
hold on
sensors = cleaned_downsampled_data.grad;
ft_plot_sens(sensors, 'facecolor', colours, 'facealpha', 0.8);

ft_plot_vol(ft_convert_units(headmodel_eeg, 'cm'));

view([-45 25])

Showing where the channels are (we will use these throughout)