Connectivity
Contents
- Read raw data, rest, eyes open.
- View the data
- Warning! cleaning data with PCA introduced false synchrony to dipoles
- Read somatosensory evoked responses.
- Averaging left and right responses separately
- Plot left finger response
- Plot right finger response
- Head model
- Dipole fit
- Plot dipole timecourse
- Now see the somatosensory cortex activation in the rest state (Mu).
- Make Hilbert envelope
- check correlations (don't expect much)
- let's check cross correlation for absolute MEG data
- now for alpha hilbert envelopes
- Coherence
- Another approach to left-right coherence
- Compute coherence for pairs of voxels with SAMwts
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)
Published with MATLAB® 7.11