[FieldTrip] Non-parametric granger causality stripe artifact

[Canedo Riedel, Walter]

Dear community,

My name is Walter Cañedo Riedel and I am working at Kelschlab (Mainz University). I am currently analyzing LFP data; there are several channels and I am trying to obtain their Granger causality spectra (by pairs of channels) with the non-parametric approach starting with a wavelet transform. My LFP data is at a resolution of 1000Hz.

I would appreciate if you help me understand how to run the non-parametric Granger causality correctly.  The GC spectra that I obtain have all the following artifact: every second frequency is exalted with respect to its immediate neighbor frequencies, forming a zebra-stripe pattern (image: GC from a channel to another).

I made a simple simulation to look into it. Channel 1 and channel 2 are modulated by different frequencies plus noise. Ch1 from time t-3 to t-1 influences Ch2 linearly at time t.
If Ch1 is modulated by a frequency below 10 Hz, then I get the stripe artifact for its Granger causality to Ch2. (Image: GC has an artifact...)  The script for the simulation is found below after my signature.

Things that I've tried to no avail:

  *   Doing the analysis up to the Nyquist limit (by down-sampling the data to 200Hz and analyzing frequencies 1:100).
  *   High-pass filtering my data (I still get stripes even after filtering out frequencies below 15Hz).
  *   Trying out 'mtmconvol' with hanning or with dpss instead of 'wavelet' for ft_freqanalysis cfg.method.
  *   Trying cfg.output = 'fourier' or 'powandcsd'.
  *   Varying the size of the width for the wavelet or of the respective window for mtmconvol.


[cid:b088c18a-8599-4446-8130-f8bf549fe488]

[cid:7683ed51-e509-48ee-bf96-835fc2867626]


All the best,
Walter Cañedo Riedel

Central Institute of Mental Health (ZI)
AG Kelsch

+49 (0)621 1703 6203



%% Simulation
clear
lfpchannels = NaN(1200000, 2); % Two channels (CH1 & CH2) with million two hundred thousand datapoints
frequency = 1000; % I have a 1000Hz resolution

SaverGC = [];
SaverPW = [];
for fr = 1:10
    lfpchannels(:, :) = rand(size(lfpchannels)); % Initial noise in the channels

    % Modulate CH1 at (fr)Hz
    dt = 1/frequency;
    StopTime = size(lfpchannels,1)/frequency;
    t = (0:dt:StopTime-dt)';
    F = fr; % Sine wave frequency (Hz)
    Modulator1 = sin(2*pi*F*t) / 10;
    lfpchannels(1:numel(Modulator1), 1) = lfpchannels(1:numel(Modulator1), 1) + Modulator1;
    % Modulate CH2 at 38Hz
    F = 38; % Sine wave frequency (Hz)
    Modulator = sin(2*pi*F*t) / 10;
    lfpchannels(1:numel(Modulator), 2) = lfpchannels(1:numel(Modulator), 2) + Modulator;

    % Random weights for CH1(t-3:t-1) influencing CH2(t).
    FRA_s = 3;
    FRA = rand(FRA_s, 1);
    % Make CH1 influence CH2
    for p = FRA_s+1:size(lfpchannels,1)
        lfpchannels(p, 2) = lfpchannels(p, 2) + mean(FRA.*lfpchannels(p-FRA_s : p-1, 1));
    end
    % z score channels to have same variance
    lfpchannels = (lfpchannels - mean(lfpchannels,1)) ./ std(lfpchannels, 0, 1);

    %%% Make data structure with trials
    TrialMarker = 1:numel(-1.999:0.001:14)/2:(numel(-1.999:0.001:14)/2)*150 + 1;
    data            = struct;
    data.label      = {'1';'2'};
    data.fsample    = 1000;
    data.trial      = {};
    data.time       = {};
    for tr = 1:50
        data.trial{tr} = lfpchannels(TrialMarker(tr):TrialMarker(tr+2)-1, :)';
        data.time{tr} = -1.999:0.001:14;
    end
    data.sampleinfo = [1:50;101:150];
    data.trialinfo = ones(1, 50);

    %%% Run wavelet transform
    cfg        = [];
    cfg.method = 'wavelet';
    cfg.output = 'powandcsd';
    cfg.foi    = 1:100;
    cfg.width  = 7;
    cfg.toi    = (3);
    cfg.pad    = 20;
    freq = ft_freqanalysis(cfg, data);

    cfg           = [];
    cfg.method    = 'granger';
    cfg.granger.sfmethod = 'bivariate';
    Granger = ft_connectivityanalysis(cfg, freq);
    SaverGC(fr, :) = Granger.grangerspctrm(1, :);
    SaverPW(fr, :) = freq.powspctrm(1, :);
end

figure
tiledlayout(1, 2)
nexttile
imagesc(SaverPW)
xlabel("Power spectrum (Hz)")
ylabel("Modulating signal (Hz)")
nexttile
imagesc(SaverGC)
xlabel("GC spectrum (Hz)")
sgtitle("GC has an artifact for the channel modulated by a signal of less than 10Hz")


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20230718/7da13c6e/attachment-0001.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: pastedImage.png
Type: image/png
Size: 127433 bytes
Desc: pastedImage.png
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20230718/7da13c6e/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: pastedImage.png
Type: image/png
Size: 86480 bytes
Desc: pastedImage.png
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20230718/7da13c6e/attachment-0003.png>