Source analysis

At the time of Hans Berger's measurements of the EEG (1924), it was thought that the signal observed at a given electrode location was generated by the cortex just below that electrode. Now we know this is not correct: the activity of an underlying source can also be 'seen' by distant sensors, because currents flow in the tissues surrounding active neuronal populations. The spatial spread of electromagnetic fields is generally referred to as 'volume conduction' or 'field spread'. Because of volume conduction, one recording electrode may pick up the activity of multiple neuronal sources. Interpretation of scalp-recorded signals is not straightforward, as the field spread corrupts these signals.

Several approaches have been proposed to estimate the position of the neuronal sources starting from the electrode potentials measured on the scalp. It is worth emphasizing from the outset that the so-called 'EEG inverse problem' – the process of inferring the distribution of neuronal sources from the electric potential distribution on the scalp – is an ill-posed problem, which does not have a unique solution. As stated by Helmholtz in 1853, the electric potential distribution recorded on a surface can be consistent with many different distributions of generators within the volume enclosed by that surface. Fortunately, approximate solutions to the EEG inverse problem exist, when electrophysiological and neuroanatomical constraints are plugged into the laws of electrodynamics. 

There are two main classes of methods able to approximate the inverse solution: 3D distributed linear tomographies and adaptive spatial filters (a.k.a., beamformers). In what follows, we will illustrate how to apply a beamformer technique in the time domain to localize the evoked sources that generate the ERPs previously analyzed in our oddball paradigm (see the section Time-domain analysis (ERP)). The aim is to reconstruct the neuronal activations underlying the event-related potentials, when participants are presented with standard and deviant tones. In particular, we want to reconstruct and compare the source distribution of two conditions (standard vs. deviant), assuming that the underlying neural populations are active in both conditions, but to a different degree. In such situations, it is advisable to compute a common spatial filter based on the combined conditions. The reason is that, in general, the more data we use for constructing the spatial filter, the more reliable the estimate of the filter. Another advantage of computing common filters is that possible differences in source activity can be interpreted more directly in terms of power differences between conditions. 

To begin with, we will define some variables that we have already used in previous sections.

Return to the wiki
home page

%% Initializations
dataPath = 'C:\Users\CIMeC\Desktop\wikiDatasets\'; % It depends on where the
                                                   % wikiDatasets folder is
                                                   % stored
fileNameDescription = 'Oddball.set';
 
subjects = {'S01' 'S02' 'S03' ...
            'S04' 'S05'}; % the variable subjects contains the list of 
                          % subject-specific directories
nSubjs = numel(subjects);

↑ Go up

In the next step, we will define trials using trigger values, combining standard and deviant tones. Unfortunately, EEGLAB and FieldTrip have different default coordinate systems. When importing EEGLAB data in FieldTrip, the position of each electrode appears rotated. In the last part of the following snippet of code, we fix this problem.

%% Defining trials (or epochs) using trigger values
allData = cell(1,nSubjs);
for i = 1:nSubjs
    fileName = strcat(subjects{i},'Oddball.set');
    dataSet = fullfile(dataPath,subjects{i},fileName);
 
    % Defining trials using trigger values - Standard tones
    cfg = [];
    cfg.dataset = dataSet; % file path
    cfg.trialfun = 'ft_trialfun_general'; % the default trial function 
    cfg.trialdef.eventtype = 'Stimulus';
    cfg.trialdef.eventvalue = {'S 103','S 104','S 105','S 106','S 107',...
                               'S 108','S 109','S 113','S 118' ... % standard tones
                               'S 203','S 204','S 205','S 206','S 207',...
                               'S 208','S 209','S 213','S 218'}; % deviant tones
    cfg.trialdef.prestim = 0.1; % in seconds
    cfg.trialdef.poststim = 0.6; % in seconds
    cfg = ft_definetrial(cfg);
    
    % Removing trials containing discontinuities
    trialsToReject = zeros(length(cfg.trl),1);
    for eventId = 1:length(cfg.event)
        if isempty(cfg.event(eventId).type)
            tmpSample = cfg.event(eventId).sample;
            for trialInfoId = 1:length(cfg.trl)
                if (tmpSample - cfg.trl(trialInfoId,1)) * ...
                        (tmpSample - cfg.trl(trialInfoId,2)) < 0
                    trialsToReject(trialInfoId) = trialInfoId;
                end
            end
        end
    end
    trialsToReject = trialsToReject(trialsToReject ~= 0);
    cfg.trl(trialsToReject,:) = [];
    
    % Reading trials
    data = ft_preprocessing(cfg);
    
    % Adjusting electrodes' position
    templateElec = 'standard_1005.elc'; % 'standard_1005.elc'
    elec = ft_read_sens(which(templateElec));
    [~,idx2] = ismember(data.label,elec.label);
    elec.chanpos = elec.chanpos(idx2,:);
    elec.chantype = elec.chantype(idx2);
    elec.chanunit = elec.chanunit(idx2);
    elec.elecpos = elec.elecpos(idx2,:);
    elec.label = elec.label(idx2);
    data.elec = elec;
    data.cfg = []; % to save memory
    
    allData{i} = data;
    
    clear idx2 templateElec tmpSample trialInfoId trialsToReject cfg data ...
        dataSet fileName eventId elec i
end

↑ Go up

From a mathematical standpoint, one of the most critical steps in beamformer techniques is the computation of the inverse of the covariance matrix calculated from the epochs of interest. Sometimes, data are rank deficient for different reasons (e.g., artifact rejection procedures based on independent component analysis). Mathematical inversion of rank deficient covariance matrices will likely introduce distortions in the source reconstruction. Prewhitening the data can help mitigate the negative effects of rank deficiency. Explaining how prewhitening works is beyond the scope of the present tutorial. We will just illustrate how to apply a spatial whitening procedure, using the activity from the baseline. To do this, we will employ the FieldTrip function ft_denoise_prewhiten.

It could be the case that ft_denoise_prewhiten throws the following error:

Error using ft_scalingfactor (line 181)
cannot convert V to unknown

Error in ft_scalingfactor (line 94)
factor = cellfun(@ft_scalingfactor, old(), new());

Error in ft_apply_montage (line 338)
scale = ft_scalingfactor(input.chanunit, montage.chanunitold);

Error in ft_denoise_prewhiten (line 190)
dataout.elec = ft_apply_montage(datain.elec,prewhiten,'balancename', 'prewhiten');

In that case, you need to edit the ft_chanunit.m file at lines 86 and 89 as follows:

Line 86: isfield(input, 'grad') should be replaced with isfield(input, 'elec')

Line 89: isfield(input, 'grad') should be replaced with isfield(input, 'opto')

Here is the snippet of code performing the prewhitening procedure.

%% Spatial whitening of the data, using the activity from the baseline
baselineWindow = [-0.1 0];
 
allDataWh = cell(1,nSubjs);
kappa = cell(1,nSubjs);
for i = 1:nSubjs
% Selecting the baseline
    cfg = [];
    cfg.latency = baselineWindow;
    baseline = ft_selectdata(cfg,allData{i});
    
    % Computing covariance
    cfg = [];
    cfg.covariance = 'yes';
    baseline = ft_timelockanalysis(cfg,baseline);
    
    % Detecting the location of the first large drop in the singular value spectrum
    [~,s,~] = svd(baseline.cov);
    d = -diff(log10(diag(s))); d = d./std(d);
    kappa{i} = find(d>5,1,'first');
    
    cfg = [];
    cfg.kappa = kappa{i};
    allDataWh{i} = ft_denoise_prewhiten(cfg,allData{i},baseline); 
    allDataWh{i}.cfg = []; % to save memory
    
    clear baseline s d kappa cfg
end
clear i

↑ Go up

Now we need to compute the covariance matrix between all pairs of channels in the whole time window.

%% Computing the covariance matrix of the prewhitened data
allErpDataWh = cell(1,nSubjs);
 
cfg = [];
cfg.preproc.demean = 'yes';
cfg.preproc.baselinewindow = baselineWindow;
cfg.covariance = 'yes';
for i = 1:nSubjs
    allErpDataWh{i} = ft_timelockanalysis(cfg,allDataWh{i}); 
    allErpDataWh{i}.cfg = []; % to save memory
end
clear cfg i

↑ Go up

In order to solve the EEG inverse problem, we first need to solve the EEG forward problem, which refers to estimating the electric potentials at the electrodes starting from a known electrical source. The forward solution (a.k.a. forward model) requires three geometric objects:

• A geometric description of the sensors (e.g., electrode positions)
• A volume conduction model of the head (a.k.a. head model)
• A source model (a set of dipole locations).

It is important to point out that these geometric objects must be represented in the same metric units.

To increase the accuracy of source reconstruction, it is advisable to record the exact position of all electrodes, as well as to construct the head and source models based on an individual subject’s MRI. Unfortunately, these pieces of information are often missing in EEG studies. The only way to proceed is then to use standard templates. In fact, we have already assigned standard positions to our electrodes during the preprocessing phase. Now we will load the standard Boundary Element Method (bem) head model provided by FieldTrip.

%% Loading standard head model
templateHeadModel = 'standard_bem.mat';
load(which(templateHeadModel)); % loads the variable 'vol'
headModel = vol;
headModel = ft_convert_units(headModel,allErpDataWh{1}.elec.unit);
 
figure;
ft_plot_mesh(headModel.bnd(3),'facecolor','r','edgecolor','none');
ft_plot_mesh(headModel.bnd(2),'facecolor','g','edgecolor','none');
ft_plot_mesh(headModel.bnd(1),'facecolor','b','edgecolor','none');
view([45 0 0]); alpha 0.3; rotate3d
 
clear vol templateHeadModel

↑ Go up

The source model we are going to use contains a set of dipole locations defined on the cortical sheet of the averaged 'Colin27' MRI template. Prior to loading the source model, you need to download it from here (Subject00 Sourcemodel 4k).

%% Loading standard source model
load('Subject00_sourcemodel_4k.mat');
sourceModel = sourcemodel;
sourceModel = ft_convert_units(sourceModel,allErpDataWh{1}.elec.unit);
figure; ft_plot_mesh(sourceModel,'facecolor','brain','edgecolor','none')
camlight('infinite'); lighting gouraud; rotate3d
 
clear sourcemodel

↑ Go up

It is always useful to check that head model, source model and sensor array are aligned correctly.

%% Plot the final source model together with the head model and the sensor array
figure;
ft_plot_mesh(headModel.bnd(3),'facecolor','r','edgecolor','none');
ft_plot_mesh(headModel.bnd(2),'facecolor','g','edgecolor','none');
ft_plot_mesh(headModel.bnd(1),'facecolor','b','edgecolor','none');
ft_plot_mesh(sourceModel,'facecolor','brain','edgecolor','none'); camlight;
ft_plot_sens(allErpDataWh{1}.elec);
view([45 0 0]); alpha 0.3; rotate3d

↑ Go up

We are now ready to compute the forward solution (i.e., the lead field matrix) for the whitened data.

%% Computing the leadfield (forward solution) for the whitened data
leadFieldWh = cell(1,nSubjs);
for i = 1:nSubjs
    cfg = [];
    cfg.headmodel = headModel;
    cfg.sourcemodel = sourceModel;
    % input of the whitened data ensures the correct sensor definition to be used
    leadFieldWh{i} = ft_prepare_leadfield(cfg,allErpDataWh{i});
    leadFieldWh{i}.cfg = []; % to save memory
end
clear cfg i

↑ Go up

Once the forward problem is solved, it is possible to reconstruct the electric source distribution for each participant.

%% Computing source reconstruction of the prewhitened data
allSrcLcmv = cell(1,nSubjs);
for i = 1:nSubjs
    [~,s,~] = svd(allErpDataWh{i}.cov);
    d = -diff(log10(diag(s))); d = d./std(d);
    kappa = find(d>5,1,'first');
    
    cfg = [];
    cfg.method = 'lcmv'; % Linearly Constrained Minimum Variance 
                         % (van Veen et al., 1997)
    cfg.lcmv.kappa = kappa; % parameter for the regularization
    cfg.lcmv.keepfilter = 'yes'; % to keep the spatial filters in the output
    cfg.lcmv.fixedori = 'yes'; % for each location, 
                               % we assume a fixed dipole orientation
    cfg.lcmv.weightnorm = 'unitnoisegain';
    cfg.headmodel = headModel;
    cfg.sourcemodel = leadFieldWh{i};
    allSrcLcmv{i} = ft_sourceanalysis(cfg,allErpDataWh{i}); 
    allSrcLcmv{i}.cfg = []; % to save memory
    clear s d kappa
end
clear cfg i

↑ Go up

The most interesting variable is 'mom', contained in the 'avg' structure of each participant. The 'mom' field represents the so-called 'virtual channels', namely the time courses of the ERPs at the source level. Another important variable is 'filter' that contains the beamformer spatial filter. We will use the spatial filters later on to reconstruct the condition-specific virtual channels.

The following lines of code are aimed at visualizing the activity of virtual channels. Remember that so far we have analyzed data combing standard and deviant tones.

%% Visualisation of virtual channel data
cfg = [];
cfg.colormap = 'jet';
cfg.parameter = 'mom';
cfg.data_labels = {'Subj 1'};
ft_sourceplot_interactive(cfg,allSrcLcmv{1}); % source activity of 
                                              % the first participant
clear cfg

↑ Go up

Two figures pop out: one figure shows how the source activity, averaged across all dipoles, evolve over time; the other figure shows the reconstructed topographical map on the cortical surface. You can vary the latency at which the topographical map is shown, by clicking on the figure with the time course of the source activity. You can also inspect the virtual ERP at a given location, by pressing the Shift key while clicking on that location on the cortical surface.

It is convenient to compute the absolute value of the dipole moments to make the topographic maps less patchy.

%% Computing the absolute value of the dipole moment
for i = 1:nSubjs
    for j = 1:length(allSrcLcmv{i}.avg.mom)
        allSrcLcmv{i}.avg.mom{j} = abs(allSrcLcmv{i}.avg.mom{j});
    end
end
clear i j

↑ Go up

We will use again the FieldTrip function ft_sourceplot_interactive to see how the topographic maps appear now. You can also inspect two or more participants simultaneously.

%% Simultaneous visualisation of two participants
cfg = [];
cfg.colormap = 'jet';
cfg.parameter = 'mom';
cfg.data_labels = {'Subj 1','Subj 4'};
ft_sourceplot_interactive(cfg,allSrcLcmv{1},allSrcLcmv{4});
clear cfg

↑ Go up

Using the spatial filters estimated from all combined trials, we are now able to split the different conditions and perform the source analysis separately for each condition.

%% Computing source reconstruction of the whitened data separately per condition
% Computing covariance matrix for each subject and condition
allErpDataWhStd = cell(1,nSubjs);
allErpDataWhDev = cell(1,nSubjs);
cfg = [];
cfg.preproc.demean = 'yes';
cfg.preproc.baselinewindow = baselineWindow;
cfg.covariance = 'yes';
for i = 1:nSubjs
    % Standard tones
    cfg.trials = find(ismember(allDataWh{i}.trialinfo(:,1),...
        [103 104 105 106 107 108 109 113 118]));
    allErpDataWhStd{i} = ft_timelockanalysis(cfg,allDataWh{i}); 
    allErpDataWhStd{i}.cfg = []; % to save memory
    
    % Deviant tones
    cfg.trials = find(ismember(allDataWh{i}.trialinfo(:,1),...
        [203 204 205 206 207 208 209 213 218]));
    allErpDataWhDev{i} = ft_timelockanalysis(cfg,allDataWh{i}); 
    allErpDataWhDev{i}.cfg = []; % to save memory
end
 
% Computing source reconstruction for each subject and condition
allSrcLcmvStd = cell(1,nSubjs);
allSrcLcmvDev = cell(1,nSubjs);
for i = 1:nSubjs
    cfg = [];
    cfg.method = 'lcmv';
    cfg.lcmv.keepfilter = 'yes';
    % With precomputed spatial filters, options for fixed orientation or 
    % weight normalisation are not applied
    %cfg.lcmv.fixedori = 'yes';
    %cfg.lcmv.weightnorm = 'unitnoisegain';
    cfg.headmodel = headModel;
    cfg.sourcemodel = leadFieldWh{i};
    cfg.sourcemodel.filter = allSrcLcmv{i}.avg.filter;
    cfg.sourcemodel.filterdimord = allSrcLcmv{i}.avg.filterdimord;
    
    % Standard tones
    [~,s,~] = svd(allErpDataWhStd{i}.cov);
    d = -diff(log10(diag(s))); d = d./std(d);
    kappa = find(d>5,1,'first');
    cfg.lcmv.kappa = kappa;
    allSrcLcmvStd{i} = ft_sourceanalysis(cfg,allErpDataWhStd{i});
    allSrcLcmvStd{i}.cfg = []; % to save memory
    clear s d kappa
    
    % Deviant tones
    [~,s,~] = svd(allErpDataWhDev{i}.cov);
    d = -diff(log10(diag(s))); d = d./std(d);
    kappa = find(d>5,1,'first');
    cfg.lcmv.kappa = kappa;
    allSrcLcmvDev{i} = ft_sourceanalysis(cfg,allErpDataWhDev{i});
    allSrcLcmvDev{i}.cfg = []; % to save memory
    clear s d kappa
    
    % Computing the absolute value of the dipole moment
    for j = 1:length(allSrcLcmvStd{i}.avg.mom)
        allSrcLcmvStd{i}.avg.mom{j} = abs(allSrcLcmvStd{i}.avg.mom{j});
    end
    for j = 1:length(allSrcLcmvDev{i}.avg.mom)
        allSrcLcmvDev{i}.avg.mom{j} = abs(allSrcLcmvDev{i}.avg.mom{j});
    end
end
clear cfg i j

↑ Go up

In the next snippet of code, we will first calculate the grand average source activity for standard and deviant tones, and then we will plot the results.

%% Computing grand average source activity for each condition
cfg = [];
cfg.parameter = 'mom';
gAvgSrcLcmvStd = ft_sourcegrandaverage(cfg,allSrcLcmvStd{:});
gAvgSrcLcmvStd.tri = allSrcLcmvStd{1}.tri;
gAvgSrcLcmvStd.cfg = []; % to save memory
gAvgSrcLcmvDev = ft_sourcegrandaverage(cfg,allSrcLcmvDev{:});
gAvgSrcLcmvDev.tri = allSrcLcmvDev{1}.tri;
gAvgSrcLcmvDev.cfg = []; % to save memory
clear cfg
 
%% Plotting grand average source activity
% Defining grand average mismatch negativity (MMN)
cfg = [];
cfg.operation = 'subtract';
cfg.parameter = 'mom';
gAvgSrcLcmvMMN = ft_math(cfg,gAvgSrcLcmvDev,gAvgSrcLcmvStd);
 
% Plotting grand average source activity
cfg = [];
cfg.parameter = 'mom';
cfg.has_diff = true;
cfg.data_labels = {'Std','Dev','MMN'};
ft_sourceplot_interactive(cfg,gAvgSrcLcmvStd,gAvgSrcLcmvDev,gAvgSrcLcmvMMN);
 
clear cfg

↑ Go up

So far, we have reconstructed the source activity for standard and deviant tones, using a source model that contains a large number of dipole positions located on the cortical surface of a standard brain. As far as we know about the EEG spatial resolution, this number of dipoles (8004, to be exact) is too large. It would not be appropriate to assume that each dipole location could indeed represent an independent neural source. It would be more appropriate to reduce the data dimensionality instead. To do this, we will adopt a parcellation scheme that has been proposed recently (Glasser et al., 2016). You need to download the brain atlas from here (Atlas MMP1.0 4k), before executing the following code.

%% Loading atlas
load('atlas_MMP1.0_4k.mat')
atlasMMP1 = atlas;
atlasMMP1 = ft_convert_units(atlasMMP1,sourceModel.unit);
clear atlas
 
%% Applying a parcellation scheme to reduce data dimensionality
atlasMMP1.pos = sourceModel.pos;
 
allSrcLcmvPrclStd = cell(1,nSubjs);
allSrcLcmvPrclDev = cell(1,nSubjs);
for i = 1:nSubjs
    cfg = [];
    cfg.parcellation = 'parcellation';
    cfg.method = 'eig';
    cfg.parameter = 'mom';
    allSrcLcmvPrclStd{i} = ft_sourceparcellate(cfg,allSrcLcmvStd{i},atlasMMP1);
    allSrcLcmvPrclStd{i}.mom = abs(allSrcLcmvPrclStd{i}.mom);
    allSrcLcmvPrclStd{i}.avg = allSrcLcmvPrclStd{i}.mom;
    allSrcLcmvPrclStd{i}.dimord = allSrcLcmvPrclStd{i}.momdimord;
    allSrcLcmvPrclStd{i} = removefields(allSrcLcmvPrclStd{i},{'mom','momdimord'});
    allSrcLcmvPrclStd{i}.cfg = []; % to save memory
    allSrcLcmvPrclDev{i} = ft_sourceparcellate(cfg,allSrcLcmvDev{i},atlasMMP1);
    allSrcLcmvPrclDev{i}.mom = abs(allSrcLcmvPrclDev{i}.mom);
    allSrcLcmvPrclDev{i}.avg = allSrcLcmvPrclDev{i}.mom;
    allSrcLcmvPrclDev{i}.dimord = allSrcLcmvPrclDev{i}.momdimord;
    allSrcLcmvPrclDev{i} = removefields(allSrcLcmvPrclDev{i},{'mom','momdimord'});
    allSrcLcmvPrclDev{i}.cfg = []; % to save memory
end
 
% Defining relevant contrast
allSrcLcmvPrclMMN = cell(1,nSubjs);
 
cfg = [];
cfg.operation = 'subtract';
cfg.parameter = 'avg';
for i = 1:nSubjs
    allSrcLcmvPrclMMN{i} = ft_math(cfg,allSrcLcmvPrclDev{i},allSrcLcmvPrclStd{i});
    allSrcLcmvPrclMMN{i}.avg = abs(allSrcLcmvPrclMMN{i}.avg);
end
clear cfg i

↑ Go up

Now we can compute the grand average of source activity parcellated for each condition.

%% Computing grand average of parcellated source activity
cfg = [];
cfg.parameter = 'avg';
gAvgSrcLcmvPrclStd = ft_timelockgrandaverage(cfg,allSrcLcmvPrclStd{:}); 
gAvgSrcLcmvPrclStd.cfg = []; % to save memory
gAvgSrcLcmvPrclDev = ft_timelockgrandaverage(cfg,allSrcLcmvPrclDev{:}); 
gAvgSrcLcmvPrclDev.cfg = []; % to save memory
gAvgSrcLcmvPrclMMN = ft_timelockgrandaverage(cfg,allSrcLcmvPrclMMN{:}); 
gAvgSrcLcmvPrclMMN.cfg = []; % to save memory
clear cfg

↑ Go up

It is worth observing that we can statistically analyze parcellated data in the same way as we have analyzed ERP data. To begin with, let us visualize the time courses of the new virtual channels. Firstly, we need to define a layout structure for the 362 parcels.

%% Preparing a new layout for virtual channels
cfg = [];
cfg.rows = 20;
cfg.columns = 20;
cfg.layout = 'ordered';
layoutPrcl = ft_prepare_layout(cfg,gAvgSrcLcmvPrclStd);
clear cfg
figure; ft_plot_layout(layoutPrcl,'interpreter','none','fontsize',8)

↑ Go up

Now let us plot the time courses of the parcellated virtual channels.

%% Plotting parcellated virtual channels
cfg = [];
cfg.layout = layoutPrcl;
cfg.linecolor = 'brg';
cfg.linestyle = {'-','-',':'};
cfg.showcomment = 'no';
cfg.showscale = 'no';
cfg.xlim = [-0.1 0.6];
cfg.ylim = [0 6];
ft_multiplotER(cfg,gAvgSrcLcmvPrclStd,gAvgSrcLcmvPrclDev,gAvgSrcLcmvPrclMMN);
clear cfg
 
%% Plotting single channel
chan = 9;
 
cfg = [];
cfg.channel = chan;
cfg.linecolor = 'brg';
cfg.linestyle = {'-','-',':'};
cfg.xlim = [-0.1 0.6];
cfg.ylim = [0 6];
ft_singleplotER(cfg,gAvgSrcLcmvPrclStd,gAvgSrcLcmvPrclDev,gAvgSrcLcmvPrclMMN);
legend({'Std','Dev','MMN'})
clear cfg chan
 
%% Plotting the brain area corresponding to a single parcel
chan = 9;
 
tmpData = allSrcLcmvPrclStd{1};
color = ones(length(tmpData.brainordinate.parcellation),3) * 0.7; % light grey
color(tmpData.brainordinate.parcellation==chan,1) = 1; % red
color(tmpData.brainordinate.parcellation==chan,2) = 0;
color(tmpData.brainordinate.parcellation==chan,3) = 0;
figure; ft_plot_mesh(tmpData.brainordinate,'vertexcolor',color)
title(['Channel ' num2str(chan) ': ' ... 
    strrep(atlasMMP1.parcellationlabel{chan},'_','\_')])
% set display properties
view([0 90]); camlight('infinite'); lighting gouraud; 
material dull; shading interp; axis vis3d equal off; rotate3d;
clear tmpData color chan

↑ Go up

To conclude, we will perform the statistical analysis at the group level and we will plot the results.

cfg.ylim = [0 6];
figure;
ft_multiplotER(cfg,gAvgSrcLcmvPrclStd,gAvgSrcLcmvPrclDev);
graphChild = get(gca,'children');
tmp = cell(length(graphChild),1);
for i = 1:length(graphChild)
    tmp{i} = graphChild(i).Type;
end
[~,idx] = sort(tmp);
set(gca,'children',graphChild(idx))
clear i tmp idx graphChild
 
%% Highlighting the significant differences on the cortical surface
color = nan(size(allSrcLcmvPrclStd{1}.brainordinate.pos,1),1);
for i = 1:numel(allSrcLcmvPrclStd{1}.label)
  color(allSrcLcmvPrclStd{1}.brainordinate.parcellation==i) = ...
      min(statVc.prob(i,:),[],2);
end
figure; ft_plot_mesh(allSrcLcmvPrclStd{1}.brainordinate,'vertexcolor',color)
view(0,90); camlight('infinite'); lighting gouraud; 
material dull; shading interp; axis vis3d equal off; rotate3d;
colormap('hot'); c = colorbar; c.Title.String = 'p-value'; caxis([0 0.1]);
clear i c color cfg

↑ Go up

We will conclude the first part of this section by saving the relevant variables to disk.

%% Saving relevant variables
save('workspaceSrcLcmv','allSrcLcmvStd','allSrcLcmvDev',...
    'allSrcLcmvPrclStd','allSrcLcmvPrclDev');

↑ Go up

Part 2

In the second part of this section, we will analyze the neural activity in the resting state, using a different approach. In particular, we will adopt the eLORETA method (Pascual-Marqui, 2007) to reconstruct the source activity in the frequency domain, and more precisely, in the alpha band. Let us start by reading and segmenting the continuous data.

%% Reading the continuous data and segmenting into 2-second epochs
cfg = [];
cfg.dataset = 'C:\Users\CIMeC\Desktop\wikiDatasets\S01RestPreprocessed.set';
cfg.continuous = 'yes';
data = ft_preprocessing(cfg);
 
cfg = [];
cfg.length = 2;
cfg.overlap = 0.5;
data = ft_redefinetrial(cfg,data);
 
% Removing the DC-component
cfg = [];
cfg.demean = 'yes';
data = ft_preprocessing(cfg,data);
 
% Making the time field constant across all trials
for i = 2:numel(data.time)
    data.time{i} = data.time{1};
end
 
% Adjusting electrodes' position
templateElec = 'standard_1005.elc'; % 'standard_1005.elc'
elec = ft_read_sens(which(templateElec));
[~,idx2] = ismember(data.label,elec.label);
elec.chanpos = elec.chanpos(idx2,:);
elec.chantype = elec.chantype(idx2);
elec.chanunit = elec.chanunit(idx2);
elec.elecpos = elec.elecpos(idx2,:);
elec.label = elec.label(idx2);
data.elec = elec;
 
clear cfg i templateElec elec idx2

↑ Go up

As a first step, we will compute the power spectrum of the data and we will plot the results. As you can see from the topography image, alpha-band activity seems to be mostly localized to the parietal cortex.

%% Computing the power spectrum
cfg = [];
cfg.output = 'pow';
cfg.foilim = [0 40];
cfg.method = 'mtmfft';
cfg.taper = 'hanning';
cfg.keeptrials = 'no';
dataPow = ft_freqanalysis(cfg,data);
 
clear cfg
 
%% Plotting the topography map and the mean power spectrum across parietal channels
cfg = [];
cfg.layout = 'easycapM1.mat';
cfg.xlim = [8 12]; % mean value in the alpha band
figure; ft_topoplotER(cfg,dataPow);
 
cfg = [];
cfg.channel = {'P1','Pz','P2'};
ft_singleplotER(cfg,dataPow);
 
clear cfg

↑ Go up

To estimate the source distribution using eLORETA, it is useful to compute the Fourier spectrum.

%% Computing the Fourier spectrum in the alpha band
cfg = [];
cfg.output = 'fourier';
cfg.foilim = [8 12];
cfg.method = 'mtmfft';
cfg.taper = 'hanning';
cfg.keeptrials = 'yes';
freq = ft_freqanalysis(cfg,data);
 
clear cfg

↑ Go up

As explained previously, we need the source model, the volume conduction model and the position of the electrodes to compute the forward solution (the lead field).

%% Loading standard head model
templateHeadModel = 'standard_bem.mat';
load(which(templateHeadModel)); % loads the variable 'vol'
headModel = vol;
headModel = ft_convert_units(headModel,freq.elec.unit);
 
figure;
ft_plot_mesh(headModel.bnd(3),'facecolor','r','edgecolor','none');
ft_plot_mesh(headModel.bnd(2),'facecolor','g','edgecolor','none');
ft_plot_mesh(headModel.bnd(1),'facecolor','b','edgecolor','none');
view([45 0 0]); alpha 0.3; rotate3d
 
clear vol templateHeadModel
 
%% Loading standard source model
load('Subject00_sourcemodel_4k.mat');
sourceModel = sourcemodel;
sourceModel = ft_convert_units(sourceModel,freq.elec.unit);
figure; ft_plot_mesh(sourceModel,'facecolor','brain','edgecolor','none')
camlight('infinite'); lighting gouraud; rotate3d
 
clear sourcemodel
 
%% Plot the final source model together with the head model and the sensor array
figure;
ft_plot_mesh(headModel.bnd(3),'facecolor','r','edgecolor','none');
ft_plot_mesh(headModel.bnd(2),'facecolor','g','edgecolor','none');
ft_plot_mesh(headModel.bnd(1),'facecolor','b','edgecolor','none');
ft_plot_mesh(sourceModel,'facecolor','brain','edgecolor','none'); camlight;
ft_plot_sens(freq.elec);
view([45 0 0]); alpha 0.3; rotate3d
 
%% Computing the leadfield (forward solution) for the whitened data
cfg = [];
cfg.headmodel = headModel;
cfg.sourcemodel = sourceModel;
leadField = ft_prepare_leadfield(cfg,freq);
leadField.cfg = []; % to save memory
 
clear cfg

↑ Go up

Now we can compute and plot the source reconstruction.

%% Computing source reconstruction using eLORETA
cfg = [];
cfg.method = 'eloreta'; % eLORETA method
cfg.eloreta.lambda = 0.05; % parameter for the regularization
cfg.eloreta.keepfilter = 'yes'; % to keep the spatial filters in the output
cfg.headmodel = headModel;
cfg.sourcemodel = leadField;
srcELoreta = ft_sourceanalysis(cfg,freq); 
srcELoreta.cfg = []; % to save memory
srcELoreta.avg.momdimord = '{pos}_ori_rpttap_freq';
 
% Plotting
cfg = [];
cfg.method = 'surface';
cfg.funparameter = 'pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim = [0 150];
cfg.funcolormap = 'jet';
cfg.colorbar = 'no';
tmpSrc = srcELoreta;
tmpSrc2 = tmpSrc;
for i = 1:length(tmpSrc.freq)
    tmpSrc2.avg.pow = tmpSrc.avg.pow(:,i);
    ft_sourceplot(cfg,tmpSrc2);
    title([num2str(tmpSrc.freq(i)) ' Hz'])
    view([0 90]); rotate3d;
    light('style','infinite','position',[0 -200 200]);
end
tmpSrc2.avg.pow = mean(tmpSrc.avg.pow,2);
ft_sourceplot(cfg,tmpSrc2);
title('Average across frequencies')
view([0 90]); rotate3d;
light('style','infinite','position',[0 -200 200]);
 
clear cfg i tmpSrc tmpSrc2

↑ Go up

Do not forget to save the relevant variables. It is worthnoting that the variable 'srcELoreta' is larger than 2 GB and to save it, we should append the '–v7.3' flag in the Matlab 'save' command. However, because saving large files to v7.3 is not optimal and causes performance issues, we will save all the variables that allow recomputing the source activity.

%% Saving relevant variables to recompute srcELoreta
save('workspaceSrcELoreta','data','freq','headModel','leadField')

↑ Go up

References

• Van Veen BD, van Drongelen W, Yuchtman M, Suzuki A (1997). Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Trans Biomed Eng, 44(9), 867-80.
• Glasser MF, Coalson TS, Robinson EC, Hacker CD, Harwell J et al. (2016). A multi-modal parcellation of human cerebral cortex. Nature, 536, 171-178.
• Pascual-Marqui RD (2007). Discrete, 3D distributed, linear imaging methods of electric neuronal activity. Part 1: exact, zero error localization. arXiv:0710.3341.