[FieldTrip] Computing source-level ROI EEG connectivity
Wanze Xie
xiew1202 at gmail.com
Fri May 8 22:24:20 CEST 2020
Hi Jan-Mathijs,
I have two short questions regarding your suggestions to Michael: 1) any
reason for why doing svd for concatenated trials would work better than
doing svd per trial?
2) If we do svd with concatenated trials do we need to re-segment the data
into epochs/trials? I think ft_freqanalysis and ft_connectivity would work
way more efficiently if they are applied to data in epochs rather than one
long trial, although I am not 100% assured.
Look forward to hearing from you.
Best,
Wanze
Schoffelen, J.M. (Jan Mathijs) <jan.schoffelen at donders.ru.nl>于2020年5月8日
周五上午10:35写道：
> Hi Michael,
>
> The WFU atlas is a volumetrically defined atlas, correct? So, in order for
> the parcellation to work, the atlas needs to be mapped onto the cortically
> constrained sheet of dipoles that you used for source reconstruction.
> This is not entirely trivial, but doable. If you obtained your atlas
> already as represented on a cortically constrained set of locations, then
> you need to ensure that these locations topologically map to the dipole
> locations of your source model.
>
> Best wishes,
> Jan-Mathijs
>
>
>
>
>
> On 8 May 2020, at 16:15, Michael Glassen - Biomedical Engineer <
> MGlassen at kesslerfoundation.org> wrote:
>
> Thank you everyone for the help!
>
> Just one final question, I got the Brodmann area atlas and the hemisphere
> atlas from this toolbox SPM “WFU-PickAtlas”.
>
> Does anyone know if there is a Brodmann area atlas that is known to work
> with fieldtrip?
>
> Best,
> Michael Glassen
> ------------------------------
> *From:* fieldtrip <fieldtrip-bounces at science.ru.nl> on behalf of
> Schoffelen, J.M. (Jan Mathijs) <jan.schoffelen at donders.ru.nl>
> *Sent:* Friday, May 8, 2020 3:17:56 AM
> *To:* FieldTrip discussion list
> *Subject:* Re: [FieldTrip] Computing source-level ROI EEG connectivity
>
> Hi Michael,
>
> Adding to Wanze’s and Xavier’s replies:
>
>
> 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?
>
>
> If ‘rawtrial’ works, than this would be fine. I was a bit worried (and I’d
> need to look in the code myself in some more detail) because historically
> the ‘rawtrial’-functionality was implemented for the beamformer algorithms,
> and not for the MNE. As such, it has never been tested thoroughly for MNE,
> so although it seems to produce an output for you, I am not 100% sure that
> it’s correct. Also, anecdotally, the ‘rawtrial’/’singletrial’ functionality
> for the beamformers is shaky in itself, and we still have it on our
> to-do-list to make the code more internally consistent and correct.
>
>
> *-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.
>
>
> -With regard to thte parcellation: understood if memory limitations
> prevent you from using it. Yet, the error message you get when using the
> trial-averaged data suggests an issue with your atlas. Where did you get it
> from? It looks as if only 2611 out of the 8196 cortical locations have a
> label. Is that according to your expectations?
> -If indeed you are to manually parcellate (with the mask), you can do one
> of the following:
> * average x/y/z across dipole locations, followed by an svd (I’d
> recommend to do this svd on the trial-concatenated data, rather than on a
> per-trial basis)
> * (to mimick the ‘eig’ method): concatenate the time series across
> dipole locations (and trials), followed by a single svd on the whole matrix.
>
>
>
>
> *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)
>
>
> Irrespective of using granger or the disrecommended dtf/pdc, all
> techniques that require a spectral factorization need padding and
> smoothing. In your case I’d recommend a padding of 4 (cfg.pad = 4), and a
> cfg.tapsmofrq somewhere in the range of 3 and 6 (but the latter is a bit of
> ‘wet finger work’, as we sould say in Dutch -> google for ’natte
> vingerwerk’ if you want to learn more).
>
> Best wishes,
> Jan-Mathijs
>
>
> *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
>
>
>
>
> _______________________________________________
> 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/20200508/4767a65b/attachment.htm>
More information about the fieldtrip
mailing list