Connectivity

Contents

Read raw data, rest, eyes open.

cd somsens
fileName='hb_c,rfhp0.1Hz';
cfg=[];
cfg.dataset=fileName;
cfg.trialdef.eventtype='TRIGGER';
cfg.trialdef.prestim=0;
cfg.trialdef.poststim=120;
cfg.trialdef.offset=0;
cfg.trialfun='BIUtrialfun';
cfg.trialdef.eventvalue= 90;
cfg.channel='MEG';
cfg=ft_definetrial(cfg);
cfg.padding=0.5;
cfg.bpfilter='yes';
cfg.bpfreq=[1 90];
cfg.demean='yes';
cfg.continuous='yes';
cfg.channel='MEG';
eyesOpen=ft_preprocessing(cfg);
evaluating trialfunction 'BIUtrialfun'
reading header...
applying digital weights in the gradiometer balancing matrix
undoing the Supine balancing
reading events from file...
applying digital weights in the gradiometer balancing matrix
undoing the Supine balancing
found 894 events
created 1 trials
the call to "ft_definetrial" took 3 seconds and an estimated 0 MB
applying digital weights in the gradiometer balancing matrix
undoing the Supine balancing
processing channel { 'A22' 'A2' 'A104' 'A241' 'A138' 'A214' 'A71' 'A26' 'A93' 'A39' 'A125' 
reading and preprocessing
reading and preprocessing trial 1 from 1
the call to "ft_preprocessing" took 5 seconds and an estimated 461 MB

View the data

see the onset of alpha for A113, A114 and A115

cfgb=[];
cfgb.layout='4D248.lay';
cfgb.continuous='yes';
cfgb.event.type='';
cfgb.event.sample=1;
cfgb.blocksize=10;
cfgb.viewmode='vertical';
cfgb.ylim=[-1e-12  1e-12];
cfgb.channel={'A112','A113','A114','A115','A116'};
comppic=ft_databrowser(cfgb,eyesOpen);
the input is raw data with 248 channels and 1 trials
detected   0 visual artifacts
redrawing with viewmode vertical
fetching data... done
fetching artifacts... done
preprocessing data... done
the call to "ft_prepare_layout" took 0 seconds and an estimated 462 MB
plotting artifacts...
plotting events...
plotting data...
done
the call to "ft_databrowser" took 7 seconds and an estimated 462 MB

Warning! cleaning data with PCA introduced false synchrony to dipoles

do not clean artifacts with component analysis before running connectivity analysis

Read somatosensory evoked responses.

reject bad trials

cfg=[];
cfg.dataset=fileName;
cfg.trialdef.eventtype='TRIGGER';
cfg.trialdef.prestim=0.1;
cfg.trialdef.poststim=0.3;
cfg.trialdef.offset=-0.1;
cfg.trialfun='BIUtrialfun';
cfg.trialdef.eventvalue= [104,102]; %left index finger
cfg1=ft_definetrial(cfg);
cfg1.demean='yes';
cfg1.baselinewindow=[-0.1 0];
cfg1.bpfilter='yes';
cfg1.bpfreq=[3 30];
cfg1.channel='MEG';
cfg1.feedback='no';
somsens=ft_preprocessing(cfg1);
cfg=[];
cfg.method='summary';
cfg.channel='MEG';
cfg.alim=1e-12;
somsensCln=ft_rejectvisual(cfg, somsens);
evaluating trialfunction 'BIUtrialfun'
reading header...
applying digital weights in the gradiometer balancing matrix
undoing the Supine balancing
reading events from file...
applying digital weights in the gradiometer balancing matrix
undoing the Supine balancing
found 894 events
created 200 trials
the call to "ft_definetrial" took 3 seconds and an estimated 0 MB
applying digital weights in the gradiometer balancing matrix
undoing the Supine balancing
processing channel { 'A22' 'A2' 'A104' 'A241' 'A138' 'A214' 'A71' 'A26' 'A93' 'A39' 'A125' 'A20'
the call to "ft_preprocessing" took 6 seconds and an estimated 116 MB
the input is raw data with 248 channels and 200 trials
showing a summary of the data for all channels and trials
computing metric [--------------------------------------------------------/]
196 trials marked as GOOD, 4 trials marked as BAD
248 channels marked as GOOD, 0 channels marked as BAD
the following trials were removed: 6, 146, 157, 181
the call to "ft_rejectvisual" took 53 seconds and an estimated 1 MB

Averaging left and right responses separately

cfg=[];
cfg.feedback='no';
cfg.trials=find(somsensCln.trialinfo==102);
rightInd=ft_timelockanalysis(cfg,somsensCln);
cfg.trials=find(somsensCln.trialinfo==104);
leftInd=ft_timelockanalysis(cfg,somsensCln);
the input is raw data with 248 channels and 196 trials
selecting 99 trials
selecting 99 trials
the call to "ft_timelockanalysis" took 0 seconds and an estimated 0 MB
the input is raw data with 248 channels and 196 trials
selecting 97 trials
selecting 97 trials
the call to "ft_timelockanalysis" took 0 seconds and an estimated 0 MB

Plot left finger response

select which channels to use for dipole fit

cfg=[];
cfg.layout='4D248.lay';
cfg.interactive='yes';
cfg.zlim=[-2e-13 2e-13];
cfg.xlim=[0.036 0.036];
figure;
ft_topoplotER(cfg,leftInd);
reading layout from file 4D248.lay
the call to "ft_prepare_layout" took 0 seconds and an estimated 0 MB
the call to "ft_topoplotTFR" took 0 seconds and an estimated 0 MB
the call to "ft_topoplotER" took 0 seconds and an estimated 0 MB

Plot right finger response

here too you have to select channels

figure;
ft_topoplotER(cfg,rightInd);
reading layout from file 4D248.lay
the call to "ft_prepare_layout" took 0 seconds and an estimated 0 MB
the call to "ft_topoplotTFR" took 0 seconds and an estimated 0 MB
the call to "ft_topoplotER" took 0 seconds and an estimated 0 MB

Head model

if exist('headmodel.mat','file')
    load headmodel
else
    [vol,grid,mesh,M1]=headmodel_BIU([],[],[],[],'localspheres');
    save headmodel vol grid mesh M1
end

Dipole fit

Here you use the channels selected at " Plot right finger response "

channelRH = {'A1', 'A2', 'A3', 'A4', 'A12', 'A13', 'A14', 'A15', 'A16', 'A17', 'A28', 'A29', 'A30', 'A31', 'A32', 'A33', 'A34', 'A35', 'A50', 'A51', 'A52', 'A53', 'A54', 'A55', 'A56', 'A57', 'A58', 'A78', 'A79', 'A80', 'A81', 'A82', 'A83', 'A84', 'A85', 'A86', 'A111', 'A112', 'A113', 'A114', 'A115', 'A116', 'A117', 'A144', 'A145', 'A146', 'A147', 'A148', 'A171', 'A172', 'A173', 'A174', 'A175', 'A193', 'A194', 'A210', 'A211'};
%channelRH = {'A11', 'A26', 'A27', 'A40', 'A41', 'A42', 'A43', 'A44', 'A45', 'A46' ... };
channelLH = {'A1', 'A2', 'A3', 'A7', 'A8', 'A9', 'A10', 'A11', 'A12', 'A21', 'A22', 'A23', 'A24', 'A25', 'A26', 'A27', 'A40', 'A41', 'A42', 'A43', 'A44', 'A45', 'A46', 'A65', 'A66', 'A67', 'A68', 'A69', 'A70', 'A71', 'A94', 'A95', 'A96', 'A97', 'A98', 'A99', 'A127', 'A128', 'A129', 'A155', 'A156', 'A157'};
t=0.0355;
cfg5 = [];
cfg5.latency = [t t];
cfg5.numdipoles = 1;
cfg5.vol=vol;
cfg5.feedback = 'textbar';
cfg5.gridsearch='yes';
cfg5.grad=ft_convert_units(leftInd.grad,'mm');
cfg5.grid=grid;
cfg5.channel=channelRH;
RHdip = ft_dipolefitting(cfg5, leftInd);
cfg5.channel=channelLH;
LHdip = ft_dipolefitting(cfg5, rightInd);
% plot dipoles
hs=ft_read_headshape('hs_file');
hs=ft_convert_units(hs,'mm');
hsx=hs.pnt(:,1);hsy=hs.pnt(:,2);hsz=hs.pnt(:,3);
figure;plot3(hsx,hsy,hsz,'r.');hold on;
ft_plot_dipole(LHdip.dip.pos,LHdip.dip.mom,'units','mm','color','b')
ft_plot_dipole(RHdip.dip.pos,RHdip.dip.mom,'units','mm','color','g')
converting units from 'm' to 'mm'
the input is timelock data with 248 channels and 407 timebins
using headmodel specified in the configuration
Your data and configuration allow for multiple sensor definitions.
Warning: using gradiometers specified in the configuration\n 
selected 57 channels
selected 1 topographies
creating dipole grid based on user specified dipole positions
3480 dipoles inside, 795 dipoles outside brain
the call to "ft_prepare_sourcemodel" took 0 seconds and an estimated 0 MB
Warning: The input units are unknown for points and S/unknown for conductivity 
scanning grid [------------------------------------------------------------]
found minimum after scanning on grid point [8.74663 -25.071 79.9769]
                                                        First-order 
 Iteration  Func-count       f(x)        Step-size       optimality
     0           4        0.0724747                        0.0021
     1          20        0.0689463            820        0.00155  
     2          28        0.0672847       0.408057        0.00155  
     3          32        0.0666039              1       0.000863  
     4          36        0.0662763              1        3.6e-05  
     5          40        0.0662757              1       7.09e-06  
     6          44        0.0662756              1        7.1e-06  
     7          48        0.0662754              1       1.74e-06  
     8          52        0.0662754              1       3.06e-07  
     9          56        0.0662754              1       7.75e-09  
    10          72        0.0662754      0.0354727       3.36e-09  

Local minimum possible.

fminunc stopped because the size of the current step is less than
the default value of the step size tolerance.



found minimum after non-linear optimization on [3.60683 -26.2694 78.6596]
the call to "ft_dipolefitting" took 8 seconds and an estimated 0 MB
the input is timelock data with 248 channels and 407 timebins
using headmodel specified in the configuration
Your data and configuration allow for multiple sensor definitions.
Warning: using gradiometers specified in the configuration\n 
selected 42 channels
selected 1 topographies
creating dipole grid based on user specified dipole positions
3480 dipoles inside, 795 dipoles outside brain
the call to "ft_prepare_sourcemodel" took 0 seconds and an estimated 0 MB
scanning grid [------------------------------------------------------------]
found minimum after scanning on grid point [9.45348 34.1831 87.5663]
                                                        First-order 
 Iteration  Func-count       f(x)        Step-size       optimality
     0           4         0.022875                       0.00155
     1          20        0.0212886            820        0.00152  
     2          24        0.0200368              1       0.000904  
     3          28        0.0195517              1       0.000273  
     4          32         0.019477              1       0.000289  
     5          36        0.0193226              1       0.000144  
     6          40        0.0193142              1          4e-05  
     7          44        0.0193134              1        6.9e-07  
     8          48        0.0193134              1       8.01e-09  
     9          52        0.0193134              1       4.05e-10  

Local minimum found.

Optimization completed because the size of the gradient is less than
the selected value of the function tolerance.



found minimum after non-linear optimization on [7.47257 37.7573 89.221]
the call to "ft_dipolefitting" took 7 seconds and an estimated 0 MB
converting units from 'm' to 'mm'

Plot dipole timecourse

cfg=[];cfg.channel=channelLH;
dataLH=ft_preprocessing(cfg,rightInd);

cfg.channel=channelRH;
dataRH=ft_preprocessing(cfg,leftInd);
timecourseLH=LHdip.dip.pot'*dataLH.avg;
timecourseRH=RHdip.dip.pot'*dataRH.avg;
figure;
plot(dataLH.time,timecourseLH,'b');
hold on;
plot(dataLH.time,timecourseRH,'g');
legend('LHdipole','RHdipole');
Warning: the trial definition in the configuration is inconsistent with the actual data 
Warning: reconstructing sampleinfo by assuming that the trials are consecutive segments of a continuous
recording 
preprocessing
preprocessing trial 1 from 1
the call to "ft_preprocessing" took 0 seconds and an estimated 0 MB
preprocessing
preprocessing trial 1 from 1
the call to "ft_preprocessing" took 0 seconds and an estimated 0 MB

Now see the somatosensory cortex activation in the rest state (Mu).

choose channels for data display

cfg=[];cfg.channel=channelLH;
dataLH=ft_preprocessing(cfg,eyesOpen);
cfg.channel=channelRH;
dataRH=ft_preprocessing(cfg,eyesOpen);
% time window with no MOG, 65 to 75s
smp1=nearest(dataLH.time{1,1},65);
smp2=nearest(dataLH.time{1,1},75);
timecourseLH=LHdip.dip.pot'*dataLH.trial{1,1}(:,smp1:smp2);
timecourseRH=RHdip.dip.pot'*dataRH.trial{1,1}(:,smp1:smp2);
timecourseLH=timecourseLH/max(abs(timecourseLH));
timecourseRH=timecourseRH/max(abs(timecourseRH));
figure1=figure;
set(figure1,'position',[100,100,1500,600])
plot(eyesOpen.time{1,1}(:,smp1:smp2),timecourseLH,'b');
hold on;
plot(eyesOpen.time{1,1}(:,smp1:smp2),timecourseRH+1.5,'g');
preprocessing
preprocessing trial 1 from 1
the call to "ft_preprocessing" took 0 seconds and an estimated 0 MB
preprocessing
preprocessing trial 1 from 1
the call to "ft_preprocessing" took 0 seconds and an estimated 0 MB

Make Hilbert envelope

vs=eyesOpen;
vs.trial{1,1}=[timecourseLH;timecourseRH];
vs.label={'Ldipole','Rdipole'};

cfg=[];
cfg.bpfilter='yes';
cfg.bpfreq=[7 13];
vsAlpha=ft_preprocessing(cfg,vs);
cfg.hilbert='yes';
vsHilb=ft_preprocessing(cfg,vs);
% plot the hilbert envelope and the data for the RH dipole
t1=1;t2=4; %sec
samp=round(1017.25.*[t1 t2])

plot(vsHilb.time{1,1}(1,samp(1):samp(2)),vs.trial{1,1}(1,samp(1):samp(2)),'b');
hold on;
plot(vsHilb.time{1,1}(1,samp(1):samp(2)),vsHilb.trial{1,1}(1,samp(1):samp(2)),'k');
plot(vsHilb.time{1,1}(1,samp(1):samp(2)),vsAlpha.trial{1,1}(1,samp(1):samp(2)),'m');
title('RH')
legend('raw','envelope','alpha')
figure; % now for the LH
plot(vsHilb.time{1,1}(1,samp(1):samp(2)),vs.trial{1,1}(2,samp(1):samp(2)),'b');
hold on;
plot(vsHilb.time{1,1}(1,samp(1):samp(2)),vsHilb.trial{1,1}(2,samp(1):samp(2)),'k');
plot(vsHilb.time{1,1}(1,samp(1):samp(2)),vsAlpha.trial{1,1}(2,samp(1):samp(2)),'m');
title('LH')
legend('raw','envelope','alpha')
preprocessing
preprocessing trial 1 from 1
the call to "ft_preprocessing" took 0 seconds and an estimated 0 MB
preprocessing
preprocessing trial 1 from 1
the call to "ft_preprocessing" took 0 seconds and an estimated 0 MB

samp =

        1017        4069

check correlations (don't expect much)

corAlpha=corrcoef(vsAlpha.trial{1,1}');corAlpha=corAlpha(1,2)
corHilb=corrcoef(vsHilb.trial{1,1}');corHilb=corHilb(1,2)
cor=corrcoef(vs.trial{1,1}');cor=cor(1,2)
corAbs=corrcoef(abs(vs.trial{1,1}'));corAbs=corAbs(1,2)
corAlpha =

   -0.1877


corHilb =

    0.0388


cor =

   -0.0507


corAbs =

    0.0353

let's check cross correlation for absolute MEG data

win=3000;
[r,lags]=xcorr(vs.trial{1,1}','coeff');
sZero=find(lags==0);
figure;
plot(lags((sZero-win):(sZero+win)),r((sZero-win):(sZero+win),[1 2 4]));
legend('LH - LH','LH - RH','RH - RH')

now for alpha hilbert envelopes

[r,lags]=xcorr(vsHilb.trial{1,1}','coeff');
sZero=find(lags==0);
figure;
plot(lags((sZero-win):(sZero+win)),r((sZero-win):(sZero+win),[1 2 4]));
legend('LH - LH','LH - RH','RH - RH')

Coherence

cfg           = [];
cfg.method    = 'mtmfft';
cfg.output    = 'fourier';
cfg.tapsmofrq = 4;
cfg.foi=1:50;
freq          = ft_freqanalysis(cfg, vs);
cfg           = [];
cfg.method    = 'coh';
cohLR          = ft_connectivityanalysis(cfg, freq);
figure;
plot(cohLR.freq,squeeze(cohLR.cohspctrm(1,2,:)))
xlabel('Frequency (Hz)')
ylabel('Coherence')
the input is raw data with 2 channels and 1 trials
processing trials
processing trial 1/1 nfft: 10174 samples, datalength: 10174 samples, 79 tapers
the call to "ft_freqanalysis" took 1 seconds and an estimated 0 MB
selection fourierspctrm along dimension 2
the call to "ft_connectivityanalysis" took 0 seconds and an estimated 0 MB

Another approach to left-right coherence

dataShort=eyesOpen;
dataShort.trial{1,1}=eyesOpen.trial{1,1}(:,smp1:smp2);
dataShort.time{1,1}=eyesOpen.time{1,1}(:,smp1:smp2);
%compute coh (not via fieldtrip) for every LR pair
temp={};
temp.label=dataShort.label;
temp.dimord='chan_freq';
temp.freq=(1:50);
temp.grad=dataShort.grad;
temp.powspctrm=ones(248,50);
coh=temp;
for cmbi=1:113
    chiL=find(strcmp(LRpairs{cmbi,1},dataShort.label));
    dataL=dataShort.trial{1,1}(chiL,:);
    chiR=find(strcmp(LRpairs{cmbi,2},dataShort.label));
    dataR=dataShort.trial{1,1}(chiR,:);
    [Cxy,F] = mscohere(dataL,dataR,[],[],[],dataShort.fsample);
    for freqi=1:50
        fi(freqi)=nearest(F,freqi); %frequency index
        coh.powspctrm(chiL,freqi)=Cxy(fi(freqi));
        coh.powspctrm(chiR,freqi)=Cxy(fi(freqi));
    end
end
cfg=[];
cfg.xlim=[11 11];
cfg.layout='4D248.lay';
figure;
ft_topoplotER(cfg,coh)
title('11Hz');
figure;
cfg.xlim=[3 3];
ft_topoplotER(cfg,coh)
title('3Hz')
reading layout from file 4D248.lay
the call to "ft_prepare_layout" took 0 seconds and an estimated 0 MB
the call to "ft_topoplotTFR" took 0 seconds and an estimated 0 MB
the call to "ft_topoplotER" took 0 seconds and an estimated 0 MB

ans = 

                 xlim: [11 11]
               layout: [1x1 struct]
          trackconfig: 'off'
          checkconfig: 'loose'
            checksize: 100000
         showcallinfo: 'yes'
                 ylim: 'maxmin'
                 zlim: 'maxmin'
                style: 'both'
            gridscale: 67
         interplimits: 'head'
        interpolation: 'v4'
           contournum: 6
             colorbar: 'no'
              shading: 'flat'
              comment: '06-May-2013
freq=[11 11]
powspctrm=[0.0116 1]'
           commentpos: 'leftbottom'
             fontsize: 8
             baseline: 'no'
               trials: 'all'
          interactive: 'no'
              hotkeys: 'no'
             renderer: []
               marker: 'on'
         markersymbol: 'o'
          markercolor: [0 0 0]
           markersize: 2
       markerfontsize: 8
            highlight: {'off'}
     highlightchannel: {{1x1 cell}}
      highlightsymbol: {'*'}
       highlightcolor: {[0 0 0]}
        highlightsize: {[6]}
    highlightfontsize: {[8]}
          labeloffset: 0.0050
        maskparameter: []
            component: []
       directionality: []
              channel: {248x1 cell}
            parameter: 'powspctrm'
             callinfo: [1x1 struct]
              version: [1x1 struct]
             previous: []

reading layout from file 4D248.lay
the call to "ft_prepare_layout" took 0 seconds and an estimated 0 MB
the call to "ft_topoplotTFR" took 0 seconds and an estimated 0 MB
the call to "ft_topoplotER" took 0 seconds and an estimated 0 MB

ans = 

                 xlim: [3 3]
               layout: [1x1 struct]
          trackconfig: 'off'
          checkconfig: 'loose'
            checksize: 100000
         showcallinfo: 'yes'
                 ylim: 'maxmin'
                 zlim: 'maxmin'
                style: 'both'
            gridscale: 67
         interplimits: 'head'
        interpolation: 'v4'
           contournum: 6
             colorbar: 'no'
              shading: 'flat'
              comment: '06-May-2013
freq=[3 3]
powspctrm=[0.00493 1]'
           commentpos: 'leftbottom'
             fontsize: 8
             baseline: 'no'
               trials: 'all'
          interactive: 'no'
              hotkeys: 'no'
             renderer: []
               marker: 'on'
         markersymbol: 'o'
          markercolor: [0 0 0]
           markersize: 2
       markerfontsize: 8
            highlight: {'off'}
     highlightchannel: {{1x1 cell}}
      highlightsymbol: {'*'}
       highlightcolor: {[0 0 0]}
        highlightsize: {[6]}
    highlightfontsize: {[8]}
          labeloffset: 0.0050
        maskparameter: []
            component: []
       directionality: []
              channel: {248x1 cell}
            parameter: 'powspctrm'
             callinfo: [1x1 struct]
              version: [1x1 struct]
             previous: []

Compute coherence for pairs of voxels with SAMwts

this requires running SAMcov and SAMwts first just watch the result in html version

% foi=11;
% [SAMHeader, ~, ActWgts]=readWeights('SAM/Raw,1-70Hz,Global.wts');
% cohlr=zeros(length(ActWgts),1);
% for Xi=-12:0.5:12
%     for Yi=-9:0.5:-0.5
%         for Zi=-2:0.5:15
%             [indR,~]=voxIndex([Xi,Yi,Zi],100.*[...
%                 SAMHeader.XStart SAMHeader.XEnd ...
%                 SAMHeader.YStart SAMHeader.YEnd ...
%                 SAMHeader.ZStart SAMHeader.ZEnd],...
%                 100.*SAMHeader.StepSize,0);
%             [indL,~]=voxIndex([Xi,-Yi,Zi],100.*[...
%                 SAMHeader.XStart SAMHeader.XEnd ...
%                 SAMHeader.YStart SAMHeader.YEnd ...
%                 SAMHeader.ZStart SAMHeader.ZEnd],...
%                 100.*SAMHeader.StepSize,0);
%             if ~isempty(find(ActWgts(indR,:))) && ~isempty(find(ActWgts(indL,:)))
%                 vsL=ActWgts(indL,:)*dataShort.trial{1,1};
%                 vsR=ActWgts(indR,:)*dataShort.trial{1,1};
%                 [Cxy,F] = mscohere(vsL,vsR,[],[],[],dataShort.fsample);
%                 fi=nearest(F,foi);
%                 cohlr(indR,1)=Cxy(fi);
%                 cohlr(indL,1)=cohlr(indR);
%                 display(['X Y Z = ',num2str([Xi,Yi,Zi])])
%             end
%         end
%     end
% end
% % plant ones in the midline
% Yi=0
% for Xi=-12:0.5:12
%     for Zi=-2:0.5:15
%         [indM,~]=voxIndex([Xi,Yi,Zi],100.*[...
%             SAMHeader.XStart SAMHeader.XEnd ...
%             SAMHeader.YStart SAMHeader.YEnd ...
%             SAMHeader.ZStart SAMHeader.ZEnd],...
%             100.*SAMHeader.StepSize,0);
%         cohlr(indM)=1;
%     end
% end
% % eval(['coh',conds{i},'=cohlr;'])
%
% cfg=[];
% cfg.step=5;
% cfg.boxSize=[-120 120 -90 90 -20 150];
% cfg.prefix='CohLR11Hz';
% VS2Brik(cfg,cohlr)