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.
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');We are loading three different files here:
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 DoneIn 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/HzNow 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);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 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])