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. DICS is an extension of an LCMV beamformer. The DICS beamformer is one of many different source reconstruction methods for estimating the origin of MEG (and EEG) signals. Some of the steps in the procedure below are common to other procedures, and some are not.

Common for MEG source reconstruction is that we extract a feature from our MEG data and use the head model (we created in the previous tutorial) and a source model to estimate the source of the signal. The source activity is not measured in the brain, but reconstructed “as if” measured from locations in the brain based upon a set of assumptions.

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

In the procedure, we define the source model specifying how we assume the sources of our signal to be distributed in the brain. We then calculate the lead field which specifies how activity from a point in the source model reaches the sensors. The leadfield is then used to “invert” the actual measured MEG signals to project them to the source space.

Set up general paths

Note that if you work from the path where all the downloadable parts are downloaded to, you do not 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

%% CLEAN AND SET PATHS

clear variables % clear workspace
restoredefaultpath


% CHANGE THE FieldTrip PATH ACCORDING TO YOUR SYSTEM
addpath /home/lau/matlab/fieldtrip-20180108_workshops_MEG_Nord/
ft_defaults

% CHANGE "data_path" ACCORDING TO YOUR SYSTEM PATH
data_path = '/home/share/workshops_MEG_Nord/mikkel_and_lau/data/';
meg_path = fullfile(data_path, 'sub-01/ses-meg/meg');
mri_path = fullfile(data_path, 'sub-01/ses-mri/anat');

Load files necessary for the beamforming

We are loading three different files here:

  1. The TFRs to locate interesting targets for beamformer source reconstruction
  2. The head model
  3. The preprocessed 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(fullfile(meg_path, 'combined_tfr.mat'));
load(fullfile(mri_path, 'headmodel.mat'));
load(fullfile(meg_path, 'cleaned_data.mat'));
cd(meg_path)
disp Done

Identify interesting features

In this tutorial, we will find the sources underlying the beta “rebound”: the relative increase in beta band (13-25 Hz) power that usually follows sensory-motor processes. Other features in the data

Plot all the channels in the data to find the beta rebound.

%% identify interesting features in the data

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 0]; %% 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
cfg.xlim            = [-1, 1];

ft_multiplotTFR(cfg, combined_tfr);

We crop the data to find the period of interest (0.460 to 0.660 s) and frequencies covering the beta rebound. We also define a baseline period of similar duration. We will compare the period of interest against this since the power estimates that the beamformer results in are not informative in themselves.

Feel free to change the time- and frequency-windows of interest if you like, or re-run the turorial with other features in the data.

%% Specify time-frequency windows of interest
rebound_toi  = [0.460 0.660];   % time/s
baseline_toi = [-0.500 -0.300]; % time/s
beta_foi     = [14 20];         % freq/Hz

Now use ft_redefinetrial to chop data into the time-intervals of interest.

%% time window of interest Beta Rebound
cfg = [];
cfg.toilim = rebound_toi;   
toi_rebound = ft_redefinetrial(cfg, cleaned_data);
    
cfg.toilim = baseline_toi;
toi_baseline = ft_redefinetrial(cfg, cleaned_data);

In this step, we make a combination of the time-window of interest and the baseline. We use the combined data to create a “common filter” for both the baseline and rebound periods.

%% combined data
cfg = [];
toi_combined = ft_appenddata(cfg, toi_rebound, toi_baseline);

Fourier analyses

Now we do a Fourier decompositions of the time courses in the cropped data for each of the times of interest. Here we use ft_freqanalysis similar to how we earlier calculated PSD and TFR. However, rather than focusing on the entire spectral range we confine ft_freqanalysis to only the narrow band we defined above. Also, we will not get the spectral power of the signal but keep the complex-valued Fourier representation. The complex-valued Fourier representation is used to calculate the cross-spectral density later. For this, we also need to keep the Fourier representation for the individual trials rather than averaging across trials.

Note that we compute three different Fourier decompositions: 1. One for the beta rebound time of interest. 2. One for the baseline time of interest. 3. One for the beta rebound and baseline combined.

%% fourier, decomposition Beta rebound
cfg = [];
cfg.method     = 'mtmfft';
cfg.output     = 'fourier';
cfg.taper      = 'hanning';
cfg.channel    = 'meggrad';
cfg.foilim     = beta_foi;
cfg.keeptrials = 'yes';
cfg.pad        = 'nextpow2';

% Run for each time segment and combined
fourier_rebound  = ft_freqanalysis(cfg, toi_rebound);
fourier_baseline = ft_freqanalysis(cfg, toi_baseline);
fourier_combined = ft_freqanalysis(cfg, toi_combined);

We then average across frequencies. We are not interested in doing source reconstructions of each frequency separately in this tutorial.

%% Average over frequncies
fourier_rebound.fourierspctrm = mean(fourier_rebound.fourierspctrm, 3);
fourier_rebound.freq = mean(fourier_rebound.freq);

fourier_baseline.fourierspctrm = mean(fourier_baseline.fourierspctrm, 3);
fourier_baseline.freq = mean(fourier_baseline.freq);

fourier_combined.fourierspctrm = mean(fourier_combined.fourierspctrm, 3);
fourier_combined.freq = mean(fourier_combined.freq);

Create a leadfield and grid

Create a grid around the brain and estimate the leadfield for each of the grid points in the brain. The leadfield specifies for any given source (a grid point inside the brain) is projected picked up by the MEG sensors. One might say that it says: “For a given source, if it is active, how would the different sensors see it.” The leadfiled is also sometimes called the forward model.

%% leadfield beamformer

cfg = [];
cfg.grad            = toi_baseline.grad;   % magnetometer and gradiometer specification
cfg.headmodel       = headmodel;            % headmodel used
cfg.channel         = 'meggrad';
cfg.senstype        = 'meg';

% This is where we define the source model
cfg.grid.resolution = 1;                    % resolution of grid...
cfg.grid.unit       = 'cm';                 % ...in units of cm

leadfield = ft_prepare_leadfield(cfg);

cfg.channel = 'megmag';
leadfield_mag = ft_prepare_leadfield(cfg);

Take a look at the source model by plotting it toghether with the volume conductor model:

figure('units', 'normalized', 'outerposition', [0 0 1 1]);
hold on
ft_plot_mesh(ft_convert_units(leadfield, 'mm'));
ft_plot_vol(headmodel,'facealpha',0);
view([-45 20])

What does the source model look like?

How many sources are there?

Try to explore the headmodel structure to find out. The plot command below might be of help:

lf = ft_convert_units(leadfield, 'mm');
lf.pos = lf.pos(lf.inside,:);
ft_plot_mesh(lf,'vertexcolor','r')
# ft_plot_vol(headmodel,'facealpha',.1);     # Comment/uncomment to toggle headmodel

The red points are where we assume the sources in the brain are. It can be difficult to see, so try rotating the figure.

Now that we have the Fourier transformed data, the source space, and the lead field, we are finally able to do the source analysis using the DICS beamformer.
First, however, we will take a closer look at the lead field

Plotting lead fields (code to generate not shown)

  1. The centre of the blue circle is the grid point of the source (the size is just for visibility)
  2. The hotter the red color, the more strongly the given sensor would see it
  3. Each of the sources (~2000) modelled have a lead field associated with it, which we could plot as below.
  4. Thus, the lead fields link the sources to the sensors, and do the source reconstruction we have to invert this link

Source analysis using DICS beamformer

The beamformer is a spatially adaptive filter which allows us to estimate the amount of activity at any given location in the brain. The inverse filter is based on minimizing the source power (or variance) at a given location, subject to unit-gain constraint. Unit-gain constraint means that, if a source had the power of amplitude one and was projected to the sensors by the lead field, the inverse filter applied to the sensors should then reconstruct power of amplitude one at that location.

In this step, ft_sourceanalysis will take all the previous ingredients and do the source reconstruction. We specify that we want to use DICS with cfg.method = ‘dics’. When we specify method = ‘dics’ FieldTrip automatically knows that it has to calculate the cross-spectral density from the complex Fourier representation.

Also, we give parameters that tell how the DICS beamformer will regularise noise (cfg.lambda), and ask to keep the spatial filters that are estimated (cfg.dics.keepfilter = ‘yes’).

We start by calculating a spatial filter that is shared between the rebound and the baseline based on the combined Fourier analysis. The combined filter is then used on the beta rebound and baseline data separately.

%% Get common filter
    
cfg = [];
cfg.method              = 'dics';                   % Dynamic Imaging of Coherent Sources
cfg.frequency           = fourier_rebound.freq;     % the frequency from the fourier analysis (as defined above to be 16 Hz)
cfg.grid                = leadfield;                % Our grid and the leadfield
cfg.headmodel           = headmodel;                % our headmodel (tells us how the magnetic field/electrical potential is propagated)
cfg.dics.projectnoise   = 'yes';                    % estimate noise
cfg.dics.lambda         = '10%';                    % how to regularise
cfg.dics.keepfilter     = 'yes';                    % keep the spatial filter in the output
cfg.dics.realfilter     = 'yes';                    % retain the real values
cfg.channel             = 'meggrad';
cfg.senstype            = 'meg';
cfg.grad                = toi_baseline.grad;

beamformer_combined = ft_sourceanalysis(cfg, fourier_combined);

Now we will apply the common filter to use on the beta rebound data and baseline data.

%% Copy filter to cfg
cfg.grid.filter = beamformer_combined.avg.filter;

beamformer_rebound = ft_sourceanalysis(cfg, fourier_rebound);
beamformer_baseline = ft_sourceanalysis(cfg, fourier_baseline);

Plot sensitivity profiles of two channels on the brain

Before we look at the results, let us plot the filter for a MEG channel above sensory-motor areas and a MEG channel above occipital areas to see how they “pick up” the sources. Note the centre of head bias.
This plot is the opp

“Tactile” channel

“Occipital” channel

Plot beamformer source reconstructions

We now plot the source reconstruction. We plot the power of the sources on top of the resliced structural MRI from the previous tutorial.

Before plotting the result of the source reconstruction, make an educated guess about where you will see active sources?

%% plot non-contrasted

load(fullfile(mri_path,'mri_resliced_cm.mat'));

cfg = [];
cfg.downsample = 2;
cfg.parameter = 'pow';

beam_inter = ft_sourceinterpolate(cfg, beamformer_rebound, mri_resliced_cm);

cfg = [];
cfg.funparameter = 'pow';
ft_sourceplot(cfg, beam_inter);
set(gcf,'units','normalized','outerposition',[0 0 1 1])

Do this look like what we should expect?

The center of the head bias means that the beamformer source reconstructions in themselves are not readily interpretable. The center of the head bias is a result of the way the weights of the spatial filters works.

Contrast the rebound and the baseline

Here we take the difference between the source reconstructions and divide it by the baseline power. This will give us the relative increase in power in the beta rebound time-window.

%% Make contrasts

beamformer_contrast = beamformer_rebound; % just make a copy
beamformer_contrast.avg.pow = (beamformer_rebound.avg.pow - beamformer_baseline.avg.pow) ./ beamformer_baseline.avg.pow;

Interpolate pow data onto the MRI of the participant to plot the results:

%% contrasts interpolated onto MRI
    
cfg = [];
cfg.downsample = 2;
cfg.parameter = 'pow';

beamformer_contrast_interpolated = ft_sourceinterpolate(cfg, beamformer_contrast, mri_resliced_cm);

Now plot the contrasts.

%% plot source analyses

max_pow = max(beamformer_contrast_interpolated.pow);
min_pow = min(beamformer_contrast_interpolated.pow);

cfg = [];
cfg.method          = 'ortho';
cfg.funparameter    = 'pow';
cfg.maskparameter   = cfg.funparameter;
cfg.funcolorlim     = [0.3 max_pow];
cfg.opacitylim      = [0.3 max_pow];

ft_sourceplot(cfg, beamformer_contrast_interpolated);
set(gcf,'units','normalized','outerposition',[0 0 1 1])

The peak of the contrast roughly corresponds to the sensory-motor areas, as we would have expected. Notice that there is a spread beyond the sensory-motor cortex. There are also patches spread throughout the brain which might suggest bilateral activation of the sensory cortex as well as frontal and parietal areas. It is, however, only with some care that we can make inferences regarding location. What we regard as “active” is dependent on the threshold we apply. Try, e.g., to change the color scale of the plots above and plot again:

cfg.funcolorlim     = [0.3 max_pow];
cfg.opacitylim      = [0.3 max_pow];

This limits the “active” area to only the sensory-motor cortex.

Different thresholds can be used to explore different aspects of the data. In the end, we might want to set a threshold based on a statistical threshold derived from multiple subjects. How to do this is covered in another workshop.

Summary

This tutorial concludes the single subject analysis from raw data to source reconstruction. The pipeline you have worked through is, as mentioned several times, one example of a way to go from raw MEG data to source reconstruction. The type of beamformer source reconstruction we presented here is useful for finding the source of frequency-specific oscillatory activity, here the origin of the beta rebound after a sensory stimulation. Other signals and other experimental procedures might call for other processing pipelines.

When you report how you did your pipeline, it is essential to make it clear what steps you did. Specify how you pre-processed your data (filters, epoch length, etc.), how you did the TFR decomposition, how you defined the source space and volume conduction/head model, and what type of source reconstruction method you used to get the source estimates. Though different method exists, most of these steps come into play in one way or another. Even though this tutorial only has presented one type of pipeline, to consider the different steps you went through is relevant for all analysis. Making it clear what you did helps reproducibility and transparency.