Connectivity analysisBrain functioning is based on the dynamic interactions between functional specialized areas. Many connectivity measures have been proposed to quantify the interactions occurring in the brain, and each metric has its own advantages and disadvantages. An important caveat to consider is that connectivity estimates between sensor recordings – namely between signals picked up by the electrodes – are difficult to interpret because of the volume conduction problem. Therefore, we strongly advise to study neuronal interactions on the level of the reconstructed sources. Another important caveat is that, even when working at the source level, different problems other than the volume conduction may affect the interpretation of the estimated brain interactions (see Bastos and Schoffelen, 2016). In the second part of the previous section (Source analysis), we have computed the source activity from signals recorded during resting state. We will recompute the variable ‘srcELoreta’ as a starting point to estimate source connectivity. |
Return to the wiki
|
%% Loading relevant variables to recompute srcELoreta
load('workspaceSrcELoreta.mat')
% Recomputing 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';
srcELoreta.avg.oridimord = '{pos}_ori_freq';
clear cfg
|
↑ Go up |
|
Notice that the field ‘srcELoreta.avg.mom’ contains three channels at each of the 8004 positions on the cortical surface. The three channels represent the x-, the y- and the z-component of the Fourier transform of the equivalent current dipole source. Computing connectivity between triplets of channels makes the interpretation of results much more difficult. Therefore, we will project the frequency representations of the source activity along the dipole orientation that explains most variance. This projection can be performed using the singular value decomposition that allow us determining the largest eigenvector. |
|
%% Projecting along the strongest dipole direction
srcELoretaProj = srcELoreta;
inside = srcELoreta.inside;
tmpMom = cell(1,length(srcELoreta.avg.mom));
for i = 1:length(srcELoreta.avg.mom)
if inside(i)
for j = 1:length(srcELoreta.freq)
% Using singular value decomposition
[u,s,~] = svd(real(srcELoreta.avg.mom{i}(:,:,j)),'econ');
tmpMom{i}(1,:,j) = u(:,1).' * srcELoreta.avg.filter{i} * ...
freq.fourierspctrm(:,:,j).';
clear u s
end
end
end
srcELoretaProj.avg.mom = tmpMom;
clear tmpMom i j inside
|
↑ Go up |
|
We will use the FieldTrip function ‘ft_connectivityanalysis’ for the actual computation of connectivity. In particular, we will estimate the imaginary part of the coherency (Nolte et al., 2004), by specifying cfg.method = ‘coh’ and cfg.complex = ‘absimag’. |
|
%% Connectivity analysis
% Computing connectivity for each frequency using the imaginary part of coherency
tmpConn = cell(1,length(srcELoretaProj.freq));
inside = leadField.inside;
cfg = [];
cfg.method = 'coh';
cfg.complex = 'absimag';
for i = 1:length(srcELoretaProj.freq)
tmpMom = srcELoretaProj.avg.mom;
for j = 1:length(tmpMom)
if inside(j)
tmpMom{j} = tmpMom{j}(:,:,i);
end
end
tmpSrc = srcELoretaProj;
tmpSrc.avg.mom = tmpMom;
tmpConn{i} = ft_connectivityanalysis(cfg,tmpSrc);
end
% Collapsing the cell array tmpConn in a single structure
srcConn = tmpConn{1};
tmp = nan(size(srcConn.cohspctrm,1),size(srcConn.cohspctrm,2), ...
length(srcELoretaProj.freq));
for i = 1:length(srcELoretaProj.freq)
tmp(:,:,i) = tmpConn{i}.cohspctrm;
end
srcConn.cohspctrm = tmp;
% Computing connectivity collapsing all frequencies
tmpSrc = srcELoretaProj;
for j = 1:length(tmpSrc.avg.mom)
if inside(j)
tmpSrc.avg.mom{j} = mean(tmpSrc.avg.mom{j},3);
end
end
srcConnAvg = ft_connectivityanalysis(cfg,tmpSrc);
srcConnAvg.freq = mean(srcConnAvg.freq);
srcConnAvg.dimord = 'pos_pos';
clear i j cfg tmp tmpMom tmpSrc tmpConn inside
|
↑ Go up |
|
The field ‘srcConn.cohspctrm’ contains the connectivity matrix estimated for each frequency. Each connectivity matrix contains about 64 million (= 8004 x 8004) elements and, because of the low spatial resolution of EEG, it is convenient to reduce the matrix dimensions using a parcellation approach (see also the section ‘Source analysis’). |
|
%% Loading atlas and applying a parcellation scheme to reduce data dimensionality
load('atlas_MMP1.0_4k.mat')
atlasMMP1 = atlas;
atlasMMP1 = ft_convert_units(atlasMMP1,leadField.unit);
clear atlas
atlasMMP1.pos = leadField.pos;
cfg = [];
cfg.parcellation = 'parcellation';
cfg.method = 'mean';
cfg.parameter = 'cohspctrm';
srcConnPrcl = ft_sourceparcellate(cfg,srcConn,atlasMMP1);
srcConnAvgPrcl = ft_sourceparcellate(cfg,srcConnAvg,atlasMMP1);
clear cfg
|
↑ Go up |
|
Now we can plot the results. |
|
%% Plotting the connectivity matrix between all pairs of parcels
for i = 1:length(srcConnPrcl.freq)+1
if i <= length(srcConnPrcl.freq)
figure; imagesc(srcConnPrcl.cohspctrm(:,:,i));
title([num2str(srcConnPrcl.freq(i)) ' Hz']);
else
figure; imagesc(srcConnAvgPrcl.cohspctrm);
title('Average across frequencies')
end
colorbar; caxis([0 0.25]);
end
clear i
% Plotting the connectivity between a pair of parcels as a function of frequency
ch1 = 10;
ch2 = 222;
figure; plot(srcConnPrcl.freq,squeeze(srcConnPrcl.cohspctrm(ch1,ch2,:)))
title(sprintf('Connectivity between %s and %s',srcConnPrcl.label{ch1},...
srcConnPrcl.label{ch2}),'Interpreter','none');
xlabel('freq (Hz)')
ylabel('Im(coherence)')
clear ch1 ch2
|
↑ Go up |
|
The imaginary part of coherency is a symmetric measure that does not take into account the direction of information flow between any pair of signals. To study directionality in connectivity, FieldTrip provides methods based on Granger causality. The interested reader can straightforwardly compute Granger causality by specifying cfg.method = ‘granger’. Before concluding this section, we will save the relevant variables. |
|
%% Saving relevant variables to disk
save('workspaceSrcConn','srcConnPrcl','srcConnAvgPrcl')
|
↑ Go up |
References
|
↑ Go up |