[FieldTrip] Granger Causality to Compare Trial Conditions

[Harrison Fisher]

Hello,

I am a student researcher at Bowdoin College, working on applying granger causality to look at information flow during theta oscillations in an episodic memory retrieval task. I have several questions about the workflow I am using to execute GC analysis and plot the results. I have 30+ subjects, each with 2 sessions, and 16 conditions for the data trials. I ultimately want to look at the difference between connectivity in one condition compared to another. 

I’ve cleaned the data for each condition separately and used a multi taper for frequency transform (‘mtmfft’). Then I localized with the DICS beam former and created a difference data set for the 2 conditions of interest in order to compute a grand average of theta power to plot on an MRI head model to select regions of interest with high theta power. 

With those regions of interest, I created two virtual channels from the time locked dataset of all conditions concatenated and performed a mtmconvol analysis to get a wavelet (fr_frequencyanalysis, mtmconvol) for each of the virtual channels. Then I computed granger over the wavelet (ft_connectivityanalysis) and use ft_singleplotTFR to visualize. 

My main question is at what point do I need to separate out by condition in order to do a comparison of granger causality on the difference between two conditions? I am working off of a previous student’s scripts, and it appears that she was running the beamformer creation of the virtual channels on a trial by trial basis (see script excerpt below). So does this mean that the construction of the virtual channels will be unaffected if I just select a smaller subset of the trials corresponding to the individual conditions? Will this affect the subsequent wavelet and granger analyses? Or is there a way to pull the conditions out of the wavelet or granger matrices created for the entire appended dataset of all the conditions? 

Thanks for the help!

Harrison Fisher
Bowdoin College 
Class of 2017

time  = timelock.All %select the appended dataset

        % apply LCMV spatial filter to location 01 of interest
        cfg = [];
        cfg.method = 'lcmv';
        cfg.lcmv.keepfilter = 'yes';
        cfg.headmodel = vol;
        cfg.elec = egi;
        cfg.grid.pos = roi1;
        source01 = ft_sourceanalysis(cfg, time);

        % construct 3-D virtual channel at location 01
        beamformer01 = source01.avg.filter{1};

        chansel = ft_channelselection('all' , data.label); % find the names
        chansel = match_str(data.label, chansel); % find the indices

        chan01_3D = [];
        chan01_3D.label = {'x', 'y', 'z'};
        chan01_3D.time = data.time;
        for i=1:length(data.trial)
            chan01_3D.trial{i} = beamformer01 * data.trial{i}(chansel,:);
        end

        % construct a single virtual channel in the maximum power orientation
        timeseries = cat(2, chan01_3D.trial{:});
        [u, s, v] = svd(timeseries, 'econ'); 
        timeseriesmaxproj = u(:, 1)' * timeseries; 
        
        chan01 = [];
        chan01.label = {'source01'};
        chan01.time = data.time;
        for i = 1:length(data.trial)
            chan01.trial{i} = u(:, 1)' * beamformer01 * data.trial{i}(chansel, :);
        end