[FieldTrip] Computing source-level ROI EEG connectivity
Wanze Xie
xiew1202 at gmail.com
Fri May 8 02:40:38 CEST 2020
Hi Michael,
In a recent paper we published
<https://bmcmedicine.biomedcentral.com/articles/10.1186/s12916-019-1431-5>,
I basically did what Jan-Mathijs recommended in his email. I generated the
eLORETA filter with the lead-field matrix and head models. I then used the
filter only to perform the cortical source reconstruction to get the
moments in 3 orientations for each dipole. The next step was to use the SVD
function (like PCA) to reduce the moments into one dimensional signal. I
copied and pasted the part of the codes I used for these steps in the end
of this email. You can also find the entire program I used here:
https://github.com/happytudouni/SourceSpaceFunctionalConnectivityAnalysis
The templates we used were children's MRI templates, so I did the
parcellation outside of Fieldtrip. You can find that part of the codes in
the program as well.
I hope this helps and good luck.
Wanze
sourcedata = [];
ntrial = size(EEG_FT.trial,1);
ndipole = find(~cellfun(@isempty,efilter)); %This array equals
find(grid.inside);
%ndipole = find(grid.inside);
ntime = size(EEG_FT.time,2);
sourcedata.trial = NaN(ntrial,length(ndipole),ntime);
tic;
for i=1:ntrial
fprintf('%.0f out of %.0f \n',i,ntrial); %Inform the process (n out of the
total trial number);
for j = 1:length(ndipole)
moments = efilter{ndipole(j)} *squeeze(EEG_FT.trial(i,:,:));
[u, s, v] = svd(moments, 'econ'); %Find the dimension
explaining/representing the most variance;
m = u(:,1)'*moments; %this m equals to or has a linear relationship with
v(:,1)', e.g., all the elements are 3 times bigger than the elements in
v(:,1)';
%Or simply do: %m = v(:,1)';
sourcedata.trial(i,j,:) = m;
end
end
analysistime = toc;
disp(['The Source analysis took ' num2str(analysistime) 's']);
Michael Glassen - Biomedical Engineer <MGlassen at kesslerfoundation.org>
于2020年5月7日周四 上午11:48写道：
> Hi Jan,
>
>
>
> Thank you so much for your help, I just had a few follow up questions.
>
>
>
> *-I would check whether the individual trials after source reconstruction
> yield meaningful values in the ‘mom’-field. Particularly, are the different
> trials different? Do the time series look meaningful to begin with (both in
> terms of signal magnitude and ’shape’)? The reason I ask, is that I am not
> fully sure whether the ‘rawtrial’ option is behaving according to your
> expectations.*
>
> *-Another way for getting the mom would be to use cfg.keepfilter in
> ft_sourceanalysis, and left-multiply each dipole location’s ’spatial
> filter’ (in sourceEEG.avg.filter) with a dataEEG.trial matrix.*
>
>
>
> I went through and plotted a couple different time series
> from the mom field from different trials and the same dipole, they are not
> equal and the shape looks ok to me, the magnitudes range from 10^-7 to
> 10^-4, is that normal? Source was calculated with everything in meters if
> that matters. Can you explain a little more the other method for
> calculating mom? The only option I have to add is keepfilter to
> sourceanalysis? And then I multiply each of these filters against the trial
> data matrix and it will leave me with what rawtrial is giving me but more
> reliably? Can I then still use the parcellate option on the output here? Or
> is this not necessary as rawtrial seems to be working?
>
>
>
> *-In order to get a ‘parcellated’ representation of your source level
> data, I recommend to use ft_sourceparcellate, which should achieve the
> parcellation you performed by hand (with the binary masks).*
>
>
>
> I tried using the ft_sourceparcellate function but I was
> running into a memory issue when ran on the rawtrial source, and when I
> remove that option and run it on the default trial-averaged sources I get a
> different error, which is pasted below. mri2 is the broddman area atlas
> interpolated onto my source model. Also, is there any way to parcellate
> twice with the function like I am doing by hand? The broddman area atlas I
> have is not divided into left and right hemispheres so I have another atlas
> I am using that is divided into hemispheres.
>
>
>
> *mri2 = *
>
>
>
> * struct with fields:*
>
>
>
> * pos: [8196×3 double]*
>
> * tri: [16384×3 double]*
>
> * unit: 'm'*
>
> * tissue: [8196×1 double]*
>
> * tissuelabel: {1×70 cell}*
>
> * cfg: [1×1 struct]*
>
>
>
> * cfg = *
>
>
>
> * struct with fields:*
>
>
>
> * method: 'mean'*
>
> * parcellation: 'tissue'*
>
> * parameter: 'mom'*
>
>
>
> *ft_sourceparcellate(cfg,sourceEEG,mri2)*
>
> *there is no "Dentate" in "tissue"*
>
> *there is no "Hypothalamus" in "tissue"*
>
> *there is no "Substania Nigra" in "tissue"*
>
> *there is no "Caudate Tail" in "tissue"*
>
> *there is no "Putamen" in "tissue"*
>
> *there is no "Medial Geniculum Body" in "tissue"*
>
> *there is no "Lateral Globus Pallidus" in "tissue"*
>
> *there is no "Subthalamic Nucleus" in "tissue"*
>
> *there is no "Medial Globus Pallidus" in "tissue"*
>
> *there is no "Lateral Geniculum Body" in "tissue"*
>
> *there is no "Anterior Commissure" in "tissue"*
>
> *there is no "Ventral Posterior Lateral Nucleus" in "tissue"*
>
> *there is no "Ventral Posterior Medial Nucleus" in "tissue"*
>
> *there is no "Ventral Lateral Nucleus" in "tissue"*
>
> *there is no "Ventral Anterior Nucleus" in "tissue"*
>
> *there is no "Caudate Body" in "tissue"*
>
> *there is no "Lateral Posterior Nucleus" in "tissue"*
>
> *there is no "Lateral Dorsal Nucleus" in "tissue"*
>
> *there is no "Midline Nucleus" in "tissue"*
>
> *there are in total 8196 positions, 8190 positions are inside the brain,
> 2611 positions have a label*
>
> *2609 of the positions inside the brain have a label*
>
> *2609 of the labeled positions are inside the brain*
>
> *5581 of the positions inside the brain do not have a label*
>
> *creating 70 parcels for parameter mom by taking the mean*
>
> *computing parcellation*
>
> *Index exceeds the number of array elements (0).*
>
>
>
> *Error in ft_sourceparcellate>cellmean1 (line 488)*
>
> *y = x{1};*
>
>
>
> *Error in ft_sourceparcellate (line 226)*
>
> * tmp(j,:,:) =
> cellmean1(dat(seg==j));*
>
>
>
> -If I cannot get source parcellate to work, does my method
> of creating a mask and pulling out the source powers seem sound? I see that
> I should use mom instead of pow, so would I then just take all the voxels
> in an ROI and average the x direction together, the y direction together
> and the z together? Then find the direction with the largest power and use
> that as my 1by x time series power(using svd as described in
> documentation)? This is what I’m assuming sourceparcellate does if I select
> ‘eig’ as the method.
>
>
>
> *ALSO IMPORTANT: if you go for a granger-like measure (granger/dtf/pdc),
> in the spectral analysis you should probably use a decent amount of zero
> padding (cfg.pad) and sufficient tapsmofrq in order to get spectral data
> that stands a chance to get a stable factorization.*
>
>
>
>
>
> -Finally, if I choose to go for granger, what kind of
> cfg.pad and tapsmofrq would you recommend? I am mostly interested in the
> beta band(12-30Hz) and my trial lengths are 357 samples(1.3 seconds)
>
>
>
> Thank you again so much for all your help!
>
>
>
> Best,
>
> Michael Glassen
>
>
>
>
>
>
>
>
>
>
>
> -
>
>
>
>
>
>
>
> *From:* fieldtrip <fieldtrip-bounces at science.ru.nl> *On Behalf Of *Schoffelen,
> J.M. (Jan Mathijs)
> *Sent:* Thursday, May 7, 2020 3:17 AM
> *To:* FieldTrip discussion list <fieldtrip at science.ru.nl>
> *Subject:* Re: [FieldTrip] Computing source-level ROI EEG connectivity
>
>
>
> Dear Michael,
>
>
>
> It looks that you already got quite far! I agree that the DTF values do
> not make much sense.
>
>
>
> a few pointers:
>
>
>
> -I would use some regularisation for the MNE reconstruction.
>
> -I would check whether the individual trials after source reconstruction
> yield meaningful values in the ‘mom’-field. Particularly, are the different
> trials different? Do the time series look meaningful to begin with (both in
> terms of signal magnitude and ’shape’)? The reason I ask, is that I am not
> fully sure whether the ‘rawtrial’ option is behaving according to your
> expectations.
>
> -Another way for getting the mom would be to use cfg.keepfilter in
> ft_sourceanalysis, and left-multiply each dipole location’s ’spatial
> filter’ (in sourceEEG.avg.filter) with a dataEEG.trial matrix.
>
> -IMPORTANT: don’t use the ‘pow’ for the subsequent spectral analysis and
> connectivity estimation. You should work with the mom.
>
> -In order to get a ‘parcellated’ representation of your source level data,
> I recommend to use ft_sourceparcellate, which should achieve the
> parcellation you performed by hand (with the binary masks).
>
> -The parcellation operation can be an average of the mom across the
> dipoles that belong to a parcel, but an alternative might be the the 1st
> principal component. The advantage of the latter would be that the
> per-parcel signal will be reduced to a single time series (rather then 3
> x/y/z orientation time courses), which is convenient (and required) for the
> interpretation of the subsequent computations.
>
> -You may want to think a bit about the depth weighting, since the MNE
> favors superficial locations in terms of amplitude, and this will also
> percolate to the way the parcellated time series are computed.
>
> -IMPORTANT: using ‘dtf’ without any additional specification of parameters
> is probably not what you want here. The reason is that under the hood a
> multivariate factorization of the spectral matrix is performed (at least:
> an attempt is made), which will miserably fail in your case. The reason for
> this failure (thinking a bit more about this, I am quite surprised that it
> worked to begin with) is the fact that you ask for a factorization of a
> size NparcelxNparcel(xfrequency bins) matrix, where the effective rank of
> the matrix is less than the number of channels in your EEG array. Even if
> Nparcel is less than the number of EEG electrodes, this factorization will
> be very poorly estimated and nonsensical. I’d recommend using ‘granger’ as
> a method, in combination with cfg.granger.sfmethod = ‘bivariate’. This
> results in a parcel-pairwise factorization (which is hopefully stable most
> of the time). I am not a fan of dtf/pdc anyway given that the measures are
> claimed to provide estimates unproblematic estimates that are interpretable
> without any caveats.
>
> -ALSO IMPORTANT: if you go for a granger-like measure (granger/dtf/pdc),
> in the spectral analysis you should probably use a decent amount of zero
> padding (cfg.pad) and sufficient tapsmofrq in order to get spectral data
> that stands a chance to get a stable factorization.
>
>
>
> So, no definitive answers unfortunately, but perhaps some food for thought
> that might get you closer to a good result.
>
>
>
> Good luck and keep up the good work,
>
>
>
> Jan-Mathijs
>
>
>
>
>
>
>
> J.M.Schoffelen, MD PhD
> Associate PI, VIDI-fellow - PI, language in interaction
> Telephone: +31-24-3614793
>
> Physical location: room 00.028
>
> Donders Centre for Cognitive Neuroimaging, Nijmegen, The Netherlands
>
>
>
>
>
> On 6 May 2020, at 22:34, Michael Glassen - Biomedical Engineer <
> MGlassen at kesslerfoundation.org> wrote:
>
>
>
> Hi all,
>
>
>
> I am having issues using fieldtrip for source reconstruction and
> connectivity analysis with preprocessed EEG data. I believe the issue is
> occurring when I try to project my sources to specific ROIs before getting
> the connectivity. I have outlined my pipeline below:
>
>
>
> -Because I do not have subject MRIs, I am using a template headmodel and
> sourcemodel that come included with fieldtrip
>
> - headmodel = standard_bem.mat
>
> -sourcemodel = cortex_8196.surf.gii
>
>
>
> -After importing EEG data, I project the EEG electrodes to the scalp using
> the code below:
>
> if strcmp(headmodel.type, 'bemcp')
>
> scalp_index = 3;
>
> elseif strcmp(headmodel.type, 'dipoli')
>
> scalp_index = 1;
>
> end
>
> cfg = [];
>
> cfg.method = 'project';
>
> cfg.headshape = headmodel.bnd(scalp_index);
>
>
>
> eegData.elec = ft_electroderealign(cfg, eegData.elec);
>
>
>
>
>
> -I then calculate covariance on a trial basis and calculate the leadfield
> matrix:
>
>
>
> cfg = [];
>
> cfg.covariance = 'yes';
>
> cfg.covariancewindow = [-inf 0];
>
> cfg.keeptrials = 'yes';
>
> tlckEeg = ft_timelockanalysis(cfg, eegData);
>
> cfg = [];
>
> cfg.elec = tlckEeg.elec; % sensor information
>
> cfg.channel = tlckEeg.label; % the used channels
>
> cfg.grid = sourcemodel; % source points
>
> cfg.headmodel = headmodel; % volume conduction model
>
> cfg.singleshell.batchsize = 5000; % speeds up the computation
>
> leadfield = ft_prepare_leadfield(cfg);
>
>
>
>
>
> -Now with the leadfield calculated, I calculate sources on a trial basis
> here:
>
>
>
> cfg = [];
>
> cfg.method = 'mne';
>
> cfg.sourcemodel = leadfield;
>
> cfg.headmodel = headmodel;
>
> cfg.mne.prewhiten = 'yes';
>
> cfg.mne.lambda = 0;
>
> cfg.mne.scalesourcecov = 'yes';
>
> cfg.rawtrial = 'yes';
>
> sourceEEG = ft_sourceanalysis(cfg,tlckEeg);
>
> sourceEEG =
>
>
>
> struct with fields:
>
>
>
> time: [1×357 double]
>
> inside: [8196×1 logical]
>
> pos: [8196×3 double]
>
> tri: [16384×3 double]
>
> method: 'rawtrial'
>
> trial: [1×400 struct]
>
> df: 400
>
> cfg: [1×1 struct]
>
>
>
>
>
> sourceEEG.trial =
>
>
>
> 1×400 struct array with fields:
>
>
>
> mom
>
> pow
>
> noisecov
>
>
>
> size(sourceEEG.trial(1).pow)
>
>
>
> ans =
>
>
>
> 8196 357
>
>
>
> -So at this point I have 400 trials, and for each of these trials I have a
> time course of source power at 8196 voxels and 357 time points for each of
> those voxels.
>
> -I then import two atlases, one that includes broddman areas and one that
> includes left and right hemispheres. I interpolate these atlases onto my
> source model and get the following:
>
> broddman =
>
> struct with fields:
>
>
>
> pos: [8196×3 double]
>
> tri: [16384×3 double]
>
> unit: 'm'
>
> tissue: [8196×1 double] (*each of
> these is a value 1-70, corresponding to which atlas entry the voxel is
> located in)*
>
> tissuelabel: {1×70 cell}*(ROI
> labels for the atlas, includes broddman areas not separated by hemisphere)*
>
> cfg: [1×1 struct]
>
> hemis2 =
>
> struct with fields:
>
> pos: [8196×3 double]
>
> tri: [16384×3 double]
>
> unit: 'm'
>
> tissue: [8196×1 double] (*each of
> these is a value 1-7, corresponding to which atlas entry the voxel is
> located in)*
>
>
>
> tissuelabel: {1×7 cell} *(ROI
> Labels for hemisphere atlas, contains 3 Left areas, 3 Right Areas and one
> inter-hemispheric)*
>
> cfg: [1×1 struct]
>
>
>
>
>
> Now I create a binary mask for each left and right broddman area. For each
> ROI in the broddman area atlas(1-70), I find all voxels that are located
> within. This gives me a list of voxels for each broddman area. I then take
> each of these lists and find which are in the left hemisphere and which are
> in the right to get left and right broddman areas.
>
> After these steps I am left with 140 binary masks, (size 1x8196), which
> equal 1 if the voxel is contained in the specific ROI and 0 if it is not.
> Using these masks, I extract the time course power for each ROI by taking
> the mean power of every voxel in the roi at each time point:
>
> roiData(i,i3,i2) =
> mean(sourceEEG.trial(i2).pow(Masks{i},i3));
>
>
>
> where i = roi number,
>
>
>
> i2 = trial number,
>
>
>
> i3 = time point
>
>
>
> In order to get this data into fieldtrip format, I put ROI data into an
> EEGLAB structure along with 6 emg channels that were recorded at the same
> time and use eeglab2fieldtrip function to have a fieldtrip structure
> containing the broddman area sources as channels. (contained in *header*)
>
>
>
> After this I perform the final step of computing dtf connectivity between
> each of the source rois and the emg channels with the code below:
>
>
>
> cfg = [];
>
> cfg.method = 'mtmfft';
>
> cfg.taper = 'dpss';
>
> cfg.output = 'fourier';
>
> cfg.tapsmofrq = 2;
>
> freq = ft_freqanalysis(cfg, header);
>
> cfg = [];
>
> cfg.method = 'dtf';
>
> dtf = ft_connectivityanalysis(cfg, freq);
>
>
>
>
>
> The values I get in DTF are extremely low, 10^-20. I am not sure why but I
> believe it has to do with the way I am projecting source activations onto
> specific ROIs. Any help or advice on how to change my pipeline is
> appreciated.
>
>
>
>
>
> Best,
>
> Michael Glassen
>
>
>
>
>
>
> _______________________________________________
> fieldtrip mailing list
> https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
> https://doi.org/10.1371/journal.pcbi.1002202
>
>
>
>
>
>
> _______________________________________________
> fieldtrip mailing list
> https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
> https://doi.org/10.1371/journal.pcbi.1002202
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20200507/9da8f273/attachment.htm>
More information about the fieldtrip
mailing list