[FieldTrip] Help with beamformer localizing tremor-coherence peaks

[Alexandra Kathrin Steina]

Dear Fieldtrip community,

I am working with MEG and LFP data from patients with Parkinson’s 
disease/essential tremor.
I want to compare how coherence between DBS-targets and cortex changes 
between tremor and tremor-free epochs at tremor frequency.
For this, I’ve calculated coherence between LFPs and MEG and used a 
dics-beamformer to localize coherent sources.
I have encountered something that I find weird/confusing when localizing 
the tremor-coherence peaks (LFP-MEG coherence at tremor frequency) of 
most of the subjects.
I've uploaded images of one example on WeTransfer under the following 
link: https://we.tl/t-B6dw2hF4aT On sensor-level there is a clear peak 
visible at tremor frequency which is absent during no-tremor (figure 1). 
 From what I see on sensor-level I would expect to also see much higher 
activity on source-level and a clear localization. However, when I apply 
dics-beamformer the maximal coherence value for the tremor-epochs is 
much lower, than for the tremor-free epochs (figure 2).
In addition, the localization seems a bit strange to me, especially when 
compared with the sensor image.
I’ve also applied a common filter, because I thought that it might solve 
the problem, but it did not.
I’ve used the exact same pipeline and exact same script to localize 
coherence between LFP and MEG at eg. high-beta or alpha frequencies and 
it gives me “good” results (see image 3, data from a different subject 
which did not have tremor).
I have used ‘fieldtrip-20201214‘ to generate the images, but have also 
tried fieldtrip versions from 2018 to the most recent one, all yielding 
similar results.

I would be grateful for advice on why this is happening!

Best,
Alex

Code-snippet I've used:

     refchannel = 'LFP-right-2B-3B';
     freqs = 1:35;

     cfg=[];
     cfg.method='mtmfft';
     cfg.output='powandcsd';
     cfg.taper= 'hanning';
     cfg.foi = freqs;
     cfg.keeptrials = 'no';
     pow_csd_tremor =ft_freqanalysis(cfg,data_seg_tremor);
     pow_csd_notremor =ft_freqanalysis(cfg,data_seg_notremor);

     cfg = [];
     cfg.method = 'coh';
     coh_tremor = ft_connectivityanalysis(cfg,pow_csd_tremor);
     coh_notremor = ft_connectivityanalysis(cfg,pow_csd_notremor);

     cfg = [];
     cfg.interactive = 'yes';
     cfg.parameter = 'cohspctrm';
     cfg.refchannel = refchannel;
     cfg.layout = 'neuromag306planar.lay';
     ft_multiplotER(cfg,coh_tremor,coh_notremor);

     %% dics beamformer

     common_filter_dics = 1;

     for f = 1:length(freqs)

         dat_appended = 
ft_appenddata([],data_seg_tremor,data_seg_notremor);

         cfg = [];
         cfg.method = 'mtmfft';
         cfg.output = 'powandcsd';
         cfg.foi = freqs(f);
         cfg.taper = 'hanning';
         cfg.keeptrials = 'yes';
         freq_tremor = ft_freqanalysis(cfg,data_seg_tremor);
         freq_notremor = ft_freqanalysis(cfg,data_seg_notremor);
         freq_appended = ft_freqanalysis(cfg,dat_appended);

         cfg=[];
         cfg.method='dics';
         cfg.dics.lambda='5%';
         cfg.dics.keepcsd = 'yes';
         cfg.dics.keepfilter = 'yes';
         cfg.dics.reafilter = 'yes';
         cfg.frequency = freqs(f);
         cfg.channel = 'MEG';
         cfg.headmodel=hdm;
         cfg.sourcemodel=grid;
         cfg.reducerank=2;
         cfg.refchan = refchannel;
         source_appended_forfilter = ft_sourceanalysis(cfg, 
freq_appended);


         cfg=[];
         cfg.method='dics';
         cfg.dics.lambda='5%';
         cfg.dics.keepcsd = 'yes';
         cfg.headmodel=hdm;
         cfg.sourcemodel=grid;
         cfg.frequency = freqs(f);
         cfg.channel = 'MEG';
         cfg.reducerank=2;
         cfg.refchan = refchannel;
         if common_filter_dics == 1
             cfg.sourcemodel.filter = 
source_appended_forfilter.avg.filter;
             cfg.sourcemodel.label = source_appended_forfilter.avg.label;
         end
         source_tremor_dics{f}= ft_sourceanalysis(cfg, freq_tremor);
         source_tremor_dics{f}.pos = template_grid.pos;
         source_tremor_dics{f}.dim = template_grid.dim;
         source_notremor_dics{f}= ft_sourceanalysis(cfg, freq_notremor);
         source_notremor_dics{f}.pos = template_grid.pos;
         source_notremor_dics{f}.dim = template_grid.dim;

     end

     %% plot all results

     freq_to_plot = 4;

     cfg = [];
     cfg.parameter = 'coh';
     cfg.interpmethod        = 'cubic';
     it_tremor_dics = 
ft_sourceinterpolate(cfg,source_tremor_dics{freq_to_plot},t1);
     it_notremor_dics = 
ft_sourceinterpolate(cfg,source_notremor_dics{freq_to_plot},t1);

     cfg =  [];
     cfg.method = 'ortho';
     cfg.funparameter = 'coh';
     cfg.maskparameter = 'coh';
     cfg.interactive = 'yes';
     figure;
     cfg.funcolorlim = [0 max(it_tremor_dics.coh)];
     ft_sourceplot(cfg,it_tremor_dics);
     cfg.funcolorlim = [0 max(it_notremor_dics.coh)];
     ft_sourceplot(cfg,it_notremor_dics);