% ck_StatsAndPlots % Takes the pre-processed fitting result tables % Calculates stats and creates figures SaveFigs=false; % autosave generated figs CloseFigs=false; % autoclose generated figs %% ======================================================================== % INITIATE --------------------------------------------------------------- % ======================================================================== %% Paths ================================================================== fprintf('Setting up paths\n') BaseFld = pwd; cd ../../../ SHARED_ROOT_FLD = pwd; cd(BaseFld) ResFld = fullfile(SHARED_ROOT_FLD,'FitResults','MULTIMODAL','cv1'); % add toolboxes (colorbrewer,matplotlib,boundedline) addpath(genpath(fullfile(SHARED_ROOT_FLD,'Toolboxes'))); def_cmap = 'Spectral'; ResType = 'mean'; TT = ['Tables_' ResType]; % figure saving folder pngfld = fullfile(pwd,'fig_png'); [~,~] = mkdir(pngfld); svgfld = fullfile(pwd,'fig_svg'); [~,~] = mkdir(svgfld); %% Load =================================================================== fprintf(['Loading results table. '... 'Please be patient, this will take a while..\n']); tic; load(fullfile(ResFld,TT)); t_load=toc; fprintf(['Loading took ' num2str(t_load) ' s\n']); %% Correct for result type ================================================ fprintf('Generalizing variable names...\n'); fprintf('MRI...'); eval(['tMRI = tMRI_' ResType '; clear tMRI_' ResType ';']); fprintf('MUA...'); eval(['tMUA = tMUA_' ResType '; clear tMUA_' ResType ';']); fprintf('LFP...'); eval(['tLFP = tLFP_' ResType '; clear tLFP_' ResType ';']); fprintf('>> DONE\n'); %% Split MRI table by model =============================================== fprintf('Splitting MRI table by model\n') m = unique(tMRI.Model); for mi = 1:length(m) M = tMRI(strcmp(tMRI.Model,m{mi}),:); T(mi).mod = M; T(mi).name = m{mi}; modidx.(m{mi}) = mi; end %% Initiate some information ============================================== fprintf('Generating ROI, model, and subject labels\n') % proces ROI info rois = {... 'V1', [34];... % occipital / visual 'V2', [131];... % occipital / visual 'V3', [60,93];... % occipital / visual 'V3A' [123];... % occipital / visual 'V4', [20,39,75];... % mid-visual 'V6', [73];... % mid-visual 'V6A', [141,56];... 'MT', [95];... % mid-visual 'MST', [99];... % mid-visual 'TEO', [125];... % temporal 'TAa', [152];... % temporal 'Tpt', [97];... % temporal (temporal/parietal) 'TPO', [159];... % temporal 'FST', [53];... % temporal 'A1', [80];... % temporal 'ML' [25];... % temporal 'AL', [46];... % temporal 'PULV', [197,198];... % subcortical 'LGN', [200,201];... % subcortical 'STR', [175];... % subcortical 'MIP', [89];... 'PIP', [90];... 'LIP', [31,130];... % parietal 'VIP', [30];... % parietal '5', [134];... % parietal '7', [91,121];... % parietal 'SI', [50,137];... % Prim Somatosensory 'SII', [63];... % Sec Somatosensory 'F2', [153];... % premotor 'F4', [146];... % premotor 'F5', [129];... % premotor 'F7', [126];... % premotor '8', [32,51,57,148];... % frontal FEF 8A: 51, 148; 8B: 32, 57 % combine these 'CINp', [6,17,27];... % posterior cingulate 'CINa', [45,98];... % anterior cingulate 'OFC', [55,107];... % frontal 'INS', [18,87,128];... % frontal 'DLPFC',[47,76,127];... % frontal 'VMPFC',[5];... % frontal }; roi=rois(:,2); roilabels=rois(:,1); MRI_MODELS={... 'linhrf_cv1_mhrf','linhrf_cv1_dhrf';... 'linhrf_cv1_mhrf_neggain','linhrf_cv1_dhrf_neggain';... 'csshrf_cv1_mhrf','csshrf_cv1_dhrf';... 'doghrf_cv1_mhrf','doghrf_cv1_dhrf';... }; ephys_MOD={'linear_ephys_cv1','linear_ephys_cv1_neggain',... 'css_ephys_cv1','dog_ephys_cv1'}; MMS={... 'LIN_m','LIN_d';... 'LIN-N_m','LIN-N_d';... 'CSS_m','CSS_d';... 'DOG_m','DOG_d';... }; cSUBS = unique(tMRI.Monkey); SUBS = cellstr(cSUBS); %% ======================================================================== % FMRI PRF ANALYSIS ------------------------------------------------------ % ======================================================================== %% Fig 2C: # significant voxels per ROI (split by monkey and model) ======= RTHRES = 5; % R2 inclusion threshold ff = figure; set(ff,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap)); set(ff,'Position',[10 10 1600 1000]); paneln=1; for s = 1:length(SUBS) fprintf(['SUB: ' SUBS{s} '\n']) for rowmod=1:4 subplot(4,2,paneln); hold on; nvox=[]; for r=1:length(roi) nvox = [nvox; sum(... T(modidx.(MRI_MODELS{rowmod,1})).mod.R2 > RTHRES & ... strcmp(T(modidx.(MRI_MODELS{rowmod,1})).mod.Monkey, cSUBS(s)) & ... ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,... ck_GetROIidx(roilabels(r),rois) ) ) ... sum(strcmp(T(modidx.(MRI_MODELS{rowmod,1})).mod.Monkey, cSUBS(s)) & ... ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,... ck_GetROIidx(roilabels(r),rois) ) ) ]; end for xval=1:size(nvox,1) bar(xval,nvox(xval,1)); end title(['SUB ' SUBS{s} ', ' MMS{rowmod,1} ' nVox with R2 > ' ... num2str(RTHRES)],'interpreter','none'); xlabel('ROI'); ylabel('nVox','interpreter','none'); set(gca, 'Box','off','yscale','log'); set(gca,'xtick',1:length(nvox),'xticklabels',roilabels,'TickDir','out'); xtickangle(45) paneln=paneln+1; end end if SaveFigs saveas(ff,fullfile(pngfld, 'MRI_nVoxSign_ROI.png')); saveas(ff,fullfile(svgfld, 'MRI_nVoxSign_ROI.svg')); end if CloseFigs; close(ff); end %% SFig 1: prop significant voxels per ROI (split by monkey and model) ==== ff2 = figure; set(ff2,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap)); set(ff2,'Position',[10 10 1600 1000]); paneln=1; for s = 1: length(SUBS) fprintf(['SUB: ' SUBS{s} '\n']) for rowmod=1:4 subplot(4,2,paneln); hold on; nvox=[]; for r=1:length(roi) nvox = [nvox; sum(... T(modidx.(MRI_MODELS{rowmod,1})).mod.R2 > RTHRES & ... strcmp(T(modidx.(MRI_MODELS{rowmod,1})).mod.Monkey, cSUBS(s)) & ... ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,... ck_GetROIidx(roilabels(r),rois) ) ) ... sum(strcmp(T(modidx.(MRI_MODELS{rowmod,1})).mod.Monkey, cSUBS(s)) & ... ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,... ck_GetROIidx(roilabels(r),rois) ) ) ]; end for xval=1:size(nvox,1) bar(xval,nvox(xval,1)./nvox(xval,2)); end title(['SUB ' SUBS{s} ', ' MMS{rowmod,1} ' proportion Vox with R2 > ' ... num2str(RTHRES)],'interpreter','none'); xlabel('ROI'); ylabel('Proportion Vox','interpreter','none'); set(gca, 'Box','off','ylim',[0 .6]); set(gca,'xtick',1:length(nvox),'xticklabels',roilabels,'TickDir','out'); xtickangle(45) paneln=paneln+1; end end if SaveFigs saveas(ff2,fullfile(pngfld, 'MRI_propVoxSign_ROI.png')); saveas(ff2,fullfile(svgfld, 'MRI_propVoxSign_ROI.svg')); end if CloseFigs; close(ff2); end %% MRI scatter plots & differences R2 per ROI ============================= RTHRES = 0; % include all voxels for sidx = 0:2 % both monkeys and individuals % scatter plots --- % Different ROIS in different colors f=figure; set(f,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap)); set(f,'Position',[10 10 1300 600]); % Tell us what's happening if sidx fprintf(['Monkey ' SUBS{sidx} '\n']); marker = [' ' SUBS{sidx}]; else fprintf(['Both monkeys\n']); marker = ' both'; end % plot scatter per area figure(f); paneln=1; for rowmod=1:4 for colmod=rowmod+1:4 subplot(2,3,paneln); hold on; PerROI(sidx+1,paneln).Models = {MRI_MODELS{rowmod,1},MRI_MODELS{colmod,1}}; PerROI(sidx+1,paneln).mR2 = []; PerROI(sidx+1,paneln).seR2 = []; for r=1:length(roi) s_R2 = T(modidx.(MRI_MODELS{rowmod,1})).mod.R2 >= RTHRES | ... T(modidx.(MRI_MODELS{colmod,1})).mod.R2 >= RTHRES; if sidx s_Monkey = strcmp(T(modidx.(MRI_MODELS{rowmod,1})).mod.Monkey,cSUBS(sidx)); SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,... ck_GetROIidx(roilabels(r),rois)) & s_Monkey; else SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,... ck_GetROIidx(roilabels(r),rois)); end scatter(... T(modidx.(MRI_MODELS{rowmod,1})).mod.R2(SSS),... T(modidx.(MRI_MODELS{colmod,1})).mod.R2(SSS),... 100,'Marker','.'); PerROI(sidx+1,paneln).roi(r).sub = marker(2:end); PerROI(sidx+1,paneln).roi(r).name = roilabels(r); PerROI(sidx+1,paneln).roi(r).R2 = [... T(modidx.(MRI_MODELS{rowmod,1})).mod.R2(SSS),... T(modidx.(MRI_MODELS{colmod,1})).mod.R2(SSS)]; PerROI(sidx+1,paneln).mR2 = [PerROI(sidx+1,paneln).mR2; ... nanmean(T(modidx.(MRI_MODELS{rowmod,1})).mod.R2(SSS)),... nanmean(T(modidx.(MRI_MODELS{colmod,1})).mod.R2(SSS))]; PerROI(sidx+1,paneln).seR2 = [PerROI(sidx+1,paneln).seR2; ... nanstd(T(modidx.(MRI_MODELS{rowmod,1})).mod.R2(SSS))./... sqrt(length(T(modidx.(MRI_MODELS{rowmod,1})).mod.R2(SSS))),... nanstd(T(modidx.(MRI_MODELS{colmod,1})).mod.R2(SSS))./... sqrt(length(T(modidx.(MRI_MODELS{colmod,1})).mod.R2(SSS)))]; % stats per ROI [PerROI(sidx+1,paneln).roi(r).p,... PerROI(sidx+1,paneln).roi(r).h,... PerROI(sidx+1,paneln).roi(r).stats] = signrank(... PerROI(sidx+1,paneln).roi(r).R2(:,1),... PerROI(sidx+1,paneln).roi(r).R2(:,2)); PerROI(sidx+1,paneln).roi(r).dir = ... diff(PerROI(sidx+1,paneln).mR2(r,:))./... abs(diff(PerROI(sidx+1,paneln).mR2(r,:))); end plot([-2 100],[-2 100],'k','Linewidth',1); title([MMS{rowmod,1} ' vs ' MMS{colmod,1}],'interpreter','none'); xlabel(MMS{rowmod,1},'interpreter','none'); ylabel(MMS{colmod,1},'interpreter','none'); set(gca, 'Box','off', 'xlim', [-2 100], 'ylim',[-2 100]); paneln=paneln+1; end end sgtitle(['SUBJECT' marker]) if SaveFigs; saveas(f,fullfile(pngfld, ... ['MRI_ModelComparison_ROI_R2' marker '.png'])); end if SaveFigs; saveas(f,fullfile(svgfld, ... ['MRI_ModelComparison_ROI_R2' marker '.svg'])); end if CloseFigs; close(f); end end % Report stats per ROI ---- fprintf('Compare CSS against P-LIN per ROI -----\n'); nBetter = sum([PerROI(1,2).roi(:).p] < 0.05 & [PerROI(1,2).roi(:).dir] > 0); nWorse = sum([PerROI(1,2).roi(:).p] < 0.05 & [PerROI(1,2).roi(:).dir] < 0); nSimilar = sum([PerROI(1,2).roi(:).p] > 0.05); nROI = length(PerROI(1,2).roi); fprintf('- Both together -\n'); fprintf([num2str(nBetter) '/' num2str(nROI) ' CSS better\n']); fprintf([num2str(nWorse) '/' num2str(nROI) ' CSS worse\n']); fprintf([num2str(nSimilar) '/' num2str(nROI) ' no difference\n']); nBetter = sum([PerROI(2,2).roi(:).p] < 0.05 & [PerROI(2,2).roi(:).dir] > 0); nWorse = sum([PerROI(2,2).roi(:).p] < 0.05 & [PerROI(2,2).roi(:).dir] < 0); nSimilar = sum([PerROI(2,2).roi(:).p] > 0.05); nROI = length(PerROI(2,2).roi); fprintf(['- ' SUBS{1} ' -\n']); fprintf([num2str(nBetter) '/' num2str(nROI) ' CSS better\n']); fprintf([num2str(nWorse) '/' num2str(nROI) ' CSS worse\n']); fprintf([num2str(nSimilar) '/' num2str(nROI) ' no difference\n']); nBetter = sum([PerROI(3,2).roi(:).p] < 0.05 & [PerROI(3,2).roi(:).dir] > 0); nWorse = sum([PerROI(3,2).roi(:).p] < 0.05 & [PerROI(3,2).roi(:).dir] < 0); nSimilar = sum([PerROI(3,2).roi(:).p] > 0.05); nROI = length(PerROI(3,2).roi); fprintf(['- ' SUBS{2} ' -\n']); fprintf([num2str(nBetter) '/' num2str(nROI) ' CSS better\n']); fprintf([num2str(nWorse) '/' num2str(nROI) ' CSS worse\n']); fprintf([num2str(nSimilar) '/' num2str(nROI) ' no difference\n']); %% Fig 4A / SFig 3C: MRI scatter plots & differences R2 all voxels ======== RTHRES = 0; % include all voxels for sidx = 0:2 % both monkeys and individuals % All ROIS together in black fsc=figure; set(fsc,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap)); set(fsc,'Position',[10 10 1300 600]); % Tell us what's happening if sidx fprintf(['Monkey ' SUBS{sidx} '\n']); marker = [' ' SUBS{sidx}]; else fprintf(['Both monkeys\n']); marker = ' both'; end % plot scatter over all voxels figure(fsc); paneln=1; for rowmod=1:4 for colmod=rowmod+1:4 subplot(2,3,paneln); hold on; XY=[]; for r=1:length(roi) s_R2 = T(modidx.(MRI_MODELS{rowmod,1})).mod.R2 >= RTHRES | ... T(modidx.(MRI_MODELS{colmod,1})).mod.R2 >= RTHRES; if sidx s_Monkey = strcmp(T(modidx.(MRI_MODELS{rowmod,1})).mod.Monkey,cSUBS(sidx)); SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,... ck_GetROIidx(roilabels(r),rois)) & s_Monkey; else SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,... ck_GetROIidx(roilabels(r),rois)); end XY=[XY; [T(modidx.(MRI_MODELS{rowmod,1})).mod.R2(SSS) ... T(modidx.(MRI_MODELS{colmod,1})).mod.R2(SSS)] ]; end binscatter(XY(:,1),XY(:,2),100); colorbar; set(gca,'ColorScale','log'); colormap(inferno) plot([-2 100],[-2 100],'k','Linewidth',1); title([MMS{rowmod,1} ' vs ' MMS{colmod,1}],'interpreter','none'); xlabel(MMS{rowmod,1},'interpreter','none'); ylabel(MMS{colmod,1},'interpreter','none'); set(gca, 'Box','off', 'xlim', [-2 100], 'ylim',[-2 100]); caxis([1 1e4]) paneln=paneln+1; end end sgtitle(['SUBJECT' marker]) if SaveFigs saveas(fsc,fullfile(pngfld,['MRI_ModelComparison_R2' marker '.png'])); saveas(fsc,fullfile(svgfld,['MRI_ModelComparison_R2' marker '.svg'])); end if CloseFigs; close(fsc); end end % stats over both subs MR2 = [T(2).mod.R2 T(4).mod.R2 T(7).mod.R2 T(8).mod.R2]; sel=logical(sum(MR2>0,2)); [p,tbl,stats] = kruskalwallis(MR2(sel,:),{'css','dog','p-lin','u-lin'}); fprintf(['Kruskal Wallis ----- p = ' num2str(p) '\n']); [c,m,h,gnames] = multcompare(stats); fprintf('Tukey HSD -----\n') for i=1:size(c,1) fprintf([gnames{c(i,1)} ' vs ' gnames{c(i,2)} ... ', p = ' num2str(c(i,6)) '\n']) end %% SFig 3A-B Difference in R2 across models per ROI ======================= for sidx = 0:2 % both monkeys and individuals % Tell us what's happening if sidx fprintf(['Monkey ' SUBS{sidx} '\n']); marker = [' ' SUBS{sidx}]; else fprintf(['Both monkeys\n']); marker = ' both'; end idx=1; for rowmod=1:4 for colmod=rowmod+1:4 diffmat{sidx+1,idx}=[]; for r=1:length(roi) if sidx s_Monkey = strcmp(T(modidx.(MRI_MODELS{rowmod,1})).mod.Monkey,cSUBS(sidx)); SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,... ck_GetROIidx(roilabels(r),rois)) & s_Monkey; else SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,... ck_GetROIidx(roilabels(r),rois)); end [n,x] = hist(... T(modidx.(MRI_MODELS{colmod,1})).mod.R2(SSS)-... T(modidx.(MRI_MODELS{rowmod,1})).mod.R2(SSS), 100); f = n./sum(SSS); m = mean(T(modidx.(MRI_MODELS{colmod,1})).mod.R2(SSS)-... T(modidx.(MRI_MODELS{rowmod,1})).mod.R2(SSS)); sd = std(T(modidx.(MRI_MODELS{colmod,1})).mod.R2(SSS)-... T(modidx.(MRI_MODELS{rowmod,1})).mod.R2(SSS)); se = sd ./ sqrt(sum(SSS)); diffmat{sidx+1,idx} = [diffmat{sidx+1,idx}; m sd se]; end idx=idx+1; end end end f2=figure; set(f2,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap)); set(f2,'Position',[100 100 1800 1000]); cc=1; for rowmod=1:4 for colmod=rowmod+1:4 subplot(3,2,cc); hold on; for xval=1:length(diffmat{1,cc}) bar(xval,diffmat{1,cc}(xval,1)); end for xval=1:length(diffmat{1,cc}) errorbar(xval,diffmat{1,cc}(xval,1),diffmat{cc}(xval,3),... 'k-','Linestyle','none'); end set(gca,'xtick',1:length(diffmat{cc}),... 'xticklabels',roilabels,'ylim',[-2 3],'TickDir','out'); xlabel('ROI'); ylabel('Diff R2'); title([MMS{colmod,1} ' - ' MMS{rowmod,1}],... 'interpreter','none'); xtickangle(45) cc=cc+1; end end if SaveFigs saveas(f2,fullfile(pngfld, 'MRI_ModelComparison_ROI_R2.png')); saveas(f2,fullfile(svgfld, 'MRI_ModelComparison_ROI_R2.svg')); end if CloseFigs; close(f2); end %% only CSS and DoG vs P-LIN per monkey =================================== f3=figure; set(f3,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap)); set(f3,'Position',[100 100 1200 1000]); subplot(2,2,1); hold on; for xval=1:length(diffmat{2,2}) bar(xval,diffmat{2,2}(xval,1)); end for xval=1:length(diffmat{2,2}) errorbar(xval,diffmat{2,2}(xval,1),diffmat{2,2}(xval,3),... 'k-','Linestyle','none'); end set(gca,'xtick',1:length(diffmat{2,2}),... 'xticklabels',roilabels,'ylim',[-0.1 2.5],'TickDir','out'); xlabel('ROI'); ylabel('Diff R2'); title('M1: CSS - P-LIN','interpreter','none'); xtickangle(45) subplot(2,2,2); hold on; for xval=1:length(diffmat{2,3}) bar(xval,diffmat{2,3}(xval,1)); end for xval=1:length(diffmat{2,3}) errorbar(xval,diffmat{2,3}(xval,1),diffmat{2,3}(xval,3),... 'k-','Linestyle','none'); end set(gca,'xtick',1:length(diffmat{2,3}),... 'xticklabels',roilabels,'ylim',[-0.1 2.5],'TickDir','out'); xlabel('ROI'); ylabel('Diff R2'); title('M1: DoG - P-LIN','interpreter','none'); xtickangle(45) subplot(2,2,3); hold on; for xval=1:length(diffmat{3,2}) bar(xval,diffmat{3,2}(xval,1)); end for xval=1:length(diffmat{3,2}) errorbar(xval,diffmat{3,2}(xval,1),diffmat{3,2}(xval,3),... 'k-','Linestyle','none'); end set(gca,'xtick',1:length(diffmat{3,2}),... 'xticklabels',roilabels,'ylim',[-0.1 2.5],'TickDir','out'); xlabel('ROI'); ylabel('Diff R2'); title('M2: CSS - P-LIN','interpreter','none'); xtickangle(45) subplot(2,2,4); hold on; for xval=1:length(diffmat{3,3}) bar(xval,diffmat{3,3}(xval,1)); end for xval=1:length(diffmat{3,3}) errorbar(xval,diffmat{3,3}(xval,1),diffmat{3,3}(xval,3),... 'k-','Linestyle','none'); end set(gca,'xtick',1:length(diffmat{3,3}),... 'xticklabels',roilabels,'ylim',[-0.1 2.5],'TickDir','out'); xlabel('ROI'); ylabel('Diff R2'); title('M2: DoG - P-LIN','interpreter','none'); xtickangle(45) if SaveFigs saveas(f3,fullfile(pngfld, 'MRI_SelModelComparison_ROI_R2.png')); saveas(f3,fullfile(svgfld, 'MRI_SelModelComparison_ROI_R2.svg')); end if CloseFigs; close(f3); end %% Fig 4B NegGain Fits: Characterize ====================================== R2th = 5; % minimum R2 R2enh = 5; % R2 improvement DoG = tMRI(... strcmp(tMRI.Model,'doghrf_cv1_mhrf'),:); lin_n = tMRI(... strcmp(tMRI.Model,'linhrf_cv1_mhrf_neggain'),:); lin = tMRI(... strcmp(tMRI.Model,'linhrf_cv1_mhrf'),:); f_neg2 = figure; set(f_neg2,'Position',[10 10 1200 1000]); vox_sel = lin_n.R2>R2th & lin_n.R2>lin.R2+R2enh & ... (lin_n.ecc<25 & lin.ecc<25 & DoG.ecc<25) & .... (lin_n.rfs<20 & lin.rfs<20 & DoG.rfs<20); fprintf('-------------\n'); subplot(2,3,1); hold on; histogram(lin_n.gain(vox_sel),-100:0.1:100,... 'Normalization','probability','FaceColor','k','FaceAlpha',0.75); xlabel('gain LIN-POSNEG');ylabel('nvoxels'); set(gca,'xlim',[-2.6 2.6]); MM=median(lin_n.gain(vox_sel)); yy=get(gca,'ylim'); plot([MM MM], [0 1],'k','Linewidth',2) set(gca,'ylim',[0 0.3],'TickDir','out'); title('Gain SELECTED voxels') fprintf(['MEDIAN GAIN SEL: ' num2str(MM) '\n']) % Wilcoxon 1-tailed < 1 [p,h,stats] = signrank(lin_n.gain(vox_sel),0,'tail','left'); fprintf(['Gain SEL < 0: Wilcoxon z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); % all vox subplot(2,3,4); hold on; histogram(lin_n.gain,-100:0.1:100,... 'Normalization','probability','FaceColor','k','FaceAlpha',0.75); xlabel('gain LIN-POSNEG');ylabel('nvoxels'); set(gca,'xlim',[-2.6 2.6]); MM=median(lin_n.gain); yy=get(gca,'ylim'); plot([MM MM], [0 1],'k','Linewidth',2) set(gca,'ylim',[0 0.3],'TickDir','out'); title('Gain ALL voxels') fprintf(['MEDIAN GAIN ALL: ' num2str(MM) '\n']) % Wilcoxon 1-tailed < 1 [p,h,stats] = signrank(lin_n.gain,0,'tail','right'); fprintf(['Gain ALL < 0: Wilcoxon z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); fprintf('-------------\n'); subplot(2,3,2); hold on; bb = [lin_n.ecc(vox_sel) lin.ecc(vox_sel)]; plot([1 2],mean(bb),'o','Linewidth',2) errorbar([1 2],mean(bb),std(bb),'k','Linewidth',2) plot([1 2],mean(bb),'o','MarkerSize',10,'MarkerFaceColor','k','MarkerEdgeColor','k') set(gca,'xtick',1:2,'xticklabels',{'LIN-N','LIN'},... 'ylim',[0 15],'xlim',[0.8 2.2],'TickDir','out') ylabel('Eccentricity'); title('Ecc') subplot(2,3,3); hold on; histogram(lin.ecc(vox_sel)-lin_n.ecc(vox_sel),-10:0.5:10,... 'FaceColor','k','FaceAlpha',0.5); xlabel('Ecc. Diff (POS-POSNEG)');ylabel('nvoxels'); MM=median(bb(:,2)-bb(:,1)); yy=get(gca,'ylim'); plot([MM MM], [0 yy(2)+40],'k','Linewidth',5) set(gca,'ylim',[0 yy(2)+100],'TickDir','out'); title('Ecc Diff') fprintf(['MEDIAN ECC DIFF: ' num2str(MM) '\n']) bbecc=bb; subplot(2,3,5); hold on; bb = [lin_n.rfs(vox_sel) lin.rfs(vox_sel)]; plot([1 2],mean(bb),'o','Linewidth',2) errorbar([1 2],mean(bb),std(bb),'k','Linewidth',2) plot([1 2],mean(bb),'o','MarkerSize',10,'MarkerFaceColor','k','MarkerEdgeColor','k') set(gca,'xtick',1:2,'xticklabels',{'LIN-N','LIN'},... 'ylim',[0 5.1],'xlim',[0.8 2.2],'TickDir','out') ylabel('Size'); title('Size') subplot(2,3,6); hold on; histogram(lin.rfs(vox_sel)-lin_n.rfs(vox_sel),-10:0.25:10,... 'FaceColor','k','FaceAlpha',0.5); xlabel('Size Diff (POS-POSNEG)');ylabel('nvoxels'); MM=median(bb(:,2)-bb(:,1)); yy=get(gca,'ylim'); plot([MM MM], [0 yy(2)+40],'k','Linewidth',5) set(gca,'ylim',[0 yy(2)+110]); set(gca,'xlim',[-5.1 5.1],'TickDir','out'); title('Size Diff') fprintf(['MEDIAN SZ DIFF: ' num2str(MM) '\n']) bbsz=bb; sgtitle('pRFs POS LINEAR vs POSNEG LINEAR model') if SaveFigs saveas(f_neg2,fullfile(pngfld, 'MRI_NEG-PRF2.png')); saveas(f_neg2,fullfile(svgfld, 'MRI_NEG-PRF2.svg')); end if CloseFigs; close(f_neg2); end % Wilcoxon 1-tailed < 0 [p,h,stats] = signrank(bbecc(:,1),bbecc(:,2)); fprintf(['Wilcoxon z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); %% DoG Fits: Characterize ================================================= f_neg3 = figure; set(f_neg3,'Position',[10 10 1300 1100]); vox_sel = DoG.R2>R2th & DoG.R2>lin.R2+R2enh; subplot(1,3,1); hold on; histogram(DoG.normamp(vox_sel),-10:0.1:10,'FaceColor','k','FaceAlpha',0.5); xlabel('INH nAMP');ylabel('nvoxels'); set(gca,'xlim',[-0.2 3.1]); MM=median(DoG.normamp(vox_sel)); yy=get(gca,'ylim'); plot([MM MM], [0 yy(2)+40],'k','Linewidth',5) set(gca,'ylim',[0 yy(2)+30],'TickDir','out'); title('NORMAMP') fprintf(['MEDIAN NAMP: ' num2str(MM) '\n']) % Wilcoxon 1-tailed > 1 [p,h,stats] = signrank(DoG.normamp(vox_sel),1,'tail','right'); fprintf(['nAmp > 1: Wilcoxon z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); fprintf(['Median nAMP: ' num2str(median(DoG.normamp(vox_sel))) ', IQR: '... num2str(iqr(DoG.normamp(vox_sel))) '\n']) subplot(1,3,2); hold on; bb = [DoG.ecc(vox_sel) lin.ecc(vox_sel)]; bbecc=bb; plot([1 2],mean(bb),'o','Linewidth',2) errorbar([1 2],mean(bb),std(bb),'k','Linewidth',2) plot([1 2],mean(bb),'o','MarkerSize',10,'MarkerFaceColor','k','MarkerEdgeColor','k') set(gca,'xtick',1:2,'xticklabels',{'DoG','LIN'},... 'ylim',[0 20],'xlim',[0.8 2.2],'TickDir','out') ylabel('Eccentricity'); title('Ecc Diff') subplot(1,3,3); hold on; histogram(lin.ecc(vox_sel)-DoG.ecc(vox_sel),-10:0.5:10,... 'FaceColor','k','FaceAlpha',0.5); xlabel('Ecc. Diff (POS-DoG)');ylabel('nvoxels'); MM=median(bb(:,2)-bb(:,1)); yy=get(gca,'ylim'); plot([MM MM], [0 yy(2)+40],'k','Linewidth',5) set(gca,'ylim',[0 yy(2)+110],'TickDir','out'); title('Ecc Diff') fprintf(['MEDIAN ECC DIFF: ' num2str(MM) '\n']) sgtitle('pRFs POS LINEAR vs DoG model') if SaveFigs saveas(f_neg3,fullfile(pngfld, 'MRI_NEG-PRF3.png')); saveas(f_neg3,fullfile(svgfld, 'MRI_NEG-PRF3.svg')); end if CloseFigs; close(f_neg3); end % Wilcoxon 1-tailed < 1 [p,h,stats] = signrank(bbecc(:,1),bbecc(:,2)); fprintf(['Wilcoxon z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); %% Value of exponential parameter for CSS across ROIs ===================== RTHRES = 5; fexp = figure; set(fexp,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap)); set(fexp,'Position',[10 10 1600 300]); hold on; for s = 1:length(SUBS) fprintf(['SUB: ' SUBS{s} '\n']) exptv(s).roi={}; for r=1:length(roi) exptvals = T(modidx.csshrf_cv1_mhrf).mod.expt(... T(modidx.csshrf_cv1_mhrf).mod.R2 > RTHRES & ... strcmp(T(modidx.csshrf_cv1_mhrf).mod.Monkey,cSUBS(s)) & ... ismember(T(modidx.csshrf_cv1_mhrf).mod.ROI,... ck_GetROIidx(roilabels(r),rois) ) ); exptv(s).roi{r}=exptvals; end end all_expt=[]; for xval=1:length(roi) if size([exptv(1).roi{xval};exptv(2).roi{xval}],1) > 4 bar(xval,mean([exptv(1).roi{xval};exptv(2).roi{xval}])); all_expt = [all_expt; exptv(1).roi{xval}; exptv(2).roi{xval}]; else bar(xval,0); all_expt = [all_expt; exptv(1).roi{xval}; exptv(2).roi{xval}]; end end for xval=1:length(roi) if size([exptv(1).roi{xval};exptv(2).roi{xval}],1) > 4 errorbar(xval,... mean([exptv(1).roi{xval};exptv(2).roi{xval}]),... std([exptv(1).roi{xval};exptv(2).roi{xval}]),... 'k','Linestyle','none'); % errorbar(xval,... % mean([exptv(1).roi{xval};exptv(2).roi{xval}]),... % std([exptv(1).roi{xval};exptv(2).roi{xval}])./... % sqrt(length([exptv(1).roi{xval};exptv(2).roi{xval}])),... % 'k','Linestyle','none'); end end title(['Avg EXPT parameter, R2 > ' ... num2str(RTHRES)],'interpreter','none'); ylabel('EXPT PM','interpreter','none'); set(gca,'xtick',1:length(roi),'xticklabels',roilabels,'TickDir','out'); xtickangle(45) if SaveFigs saveas(fexp,fullfile(pngfld, 'MRI_CSS_EXPT-PM_ROI.png')); saveas(fexp,fullfile(svgfld, 'MRI_CSS_EXPT-PM_ROI.svg')); end if CloseFigs; close(fexp); end % Wilcoxon 1-tailed < 1 [p,h,stats] = signrank(all_expt,1,'tail','left'); fprintf(['ALL VOXELS - Wilcoxon EXPT < 1 : z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); fprintf('\n------\nSTATS\n------\n') for xval=1:length(roi) if size([exptv(1).roi{xval};exptv(2).roi{xval}],1) > 4 [p,h,stats] = signrank([exptv(1).roi{xval};exptv(2).roi{xval}],1,'tail','left'); if isfield(stats,'zval') fprintf([roilabels{xval} ' VOXELS - Wilcoxon EXPT < 1 : z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); end else fprintf([roilabels{xval} ': Not enough values for statistics\n']) end end %% SFig 13: MRI scatter plot HRF & differences ============================ RTHRES = 0; f3a=figure; set(f3a,'Position',[100 100 1500 400]); set(f3a,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap)); s_R2 = T(modidx.linhrf_cv1_mhrf).mod.R2>RTHRES; f3b=figure; set(f3b,'Position',[100 100 1300 1200]); set(f3b,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap)); s_R2 = T(modidx.linhrf_cv1_mhrf).mod.R2>RTHRES; for mm = 1:size(MRI_MODELS,1) figure(f3a) subplot(1,4,mm); hold on; XY=[]; for r=1:length(roi) SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,... ck_GetROIidx(roilabels(r),rois) ); % scatter(T(modidx.(MRI_MODELS{mm,1})).mod.R2(SSS),... % T(modidx.(MRI_MODELS{mm,2})).mod.R2(SSS),100,'Marker','.'); XY=[XY; T(modidx.(MRI_MODELS{mm,1})).mod.R2(SSS) ... T(modidx.(MRI_MODELS{mm,2})).mod.R2(SSS)]; end binscatter(XY(:,1),XY(:,2),100); colorbar; set(gca,'ColorScale','log'); colormap(inferno) caxis([1 1e4]) plot([-2 100],[-2 100],'k','Linewidth',1); title(['mHRF vs cHRF (' MRI_MODELS{mm,1} ')'],'interpreter','none'); xlabel('Monkey HRF R2'); ylabel('Canonical HRF R2'); set(gca, 'Box','off', 'xlim', [-2 100], 'ylim',[-2 100]); diffmat2{1}=[]; dH=[]; for r=1:length(roi) SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,... ck_GetROIidx(roilabels(r),rois) ); [n,x] = hist(... T(modidx.(MRI_MODELS{mm,1})).mod.R2(SSS)-... T(modidx.(MRI_MODELS{mm,2})).mod.R2(SSS),100); f = n./sum(SSS); m = mean(T(modidx.(MRI_MODELS{mm,1})).mod.R2(SSS)-... T(modidx.(MRI_MODELS{mm,2})).mod.R2(SSS)); sd = std(T(modidx.(MRI_MODELS{mm,1})).mod.R2(SSS)-... T(modidx.(MRI_MODELS{mm,2})).mod.R2(SSS)); nvox = sum(SSS); se = sd ./ sqrt(nvox); diffmat2{1} = [diffmat2{1}; m sd se nvox]; dH=[dH;... T(modidx.(MRI_MODELS{mm,1})).mod.R2(SSS) - ... T(modidx.(MRI_MODELS{mm,2})).mod.R2(SSS) ]; end figure(f3b) subplot(4,1,mm); hold on for xval=1:length(diffmat2{1}) bar(xval,diffmat2{1}(xval,1)); end for xval=1:length(diffmat2{1}) errorbar(xval,diffmat2{1}(xval,1),diffmat2{1}(xval,3),... 'k-','Linestyle','none') end % mean (weighted mean taking nvox into account) mAll = sum((diffmat2{1}(:,1).*diffmat2{1}(:,4)))./... sum(diffmat2{1}(:,4)); mAll = mean(dH); stdAll = std(dH); % Wilcoxon 1-tailed < 1 [p,h,stats] = signrank(dH,0); fprintf(['ALL VOXELS - Wilcoxon dHRF (' MRI_MODELS{mm,1} ') < 1 : z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); text(0.5, -1.2,[num2str(mAll) ' +/- ' num2str(stdAll)]) set(gca,'xticklabels',[],'ylim',[-1.6 1.6],'TickDir','out'); %xlabel('ROI'); ylabel('Diff R2'); title(['mHRF - cHRF (' MRI_MODELS{mm,1} ')'],'interpreter','none'); %legend(roilabels,'interpreter','none','Location','NorthEast'); set(gca,'xtick',1:length(diffmat2{1}),'xticklabels',roilabels); xtickangle(45) end if SaveFigs saveas(f3a,fullfile(pngfld, 'HRF_Comparison.png')); saveas(f3a,fullfile(svgfld, 'HRF_Comparison.svg')); saveas(f3b,fullfile(pngfld, 'HRF_Comparison_binned.png')); saveas(f3b,fullfile(svgfld, 'HRF_Comparison_binned.svg')); end if CloseFigs; close(f3a); close(f3b); close all; end %% MRI rf size depending on HRF =========================================== R2th=5; s_R2 = T(modidx.csshrf_cv1_mhrf).mod.R2 > R2th; for mm = 1:size(MRI_MODELS,1) xy_sz{mm}=[]; xy_ecc{mm}=[]; for r=1:length(roi) SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,... ck_GetROIidx(roilabels(r),rois) ); xy_sz{mm} = [xy_sz{mm};... T(modidx.(MRI_MODELS{mm,1})).mod.rfs(SSS) ... T(modidx.(MRI_MODELS{mm,2})).mod.rfs(SSS)]; xy_ecc{mm} = [xy_ecc{mm};... T(modidx.(MRI_MODELS{mm,1})).mod.ecc(SSS) ... T(modidx.(MRI_MODELS{mm,2})).mod.ecc(SSS)]; end xy_sz{mm}( isinf(xy_sz{mm}) ) = nan; xy_ecc{mm}( isinf(xy_ecc{mm}) ) = nan; % Wilcoxon fprintf(['--' MMS{mm,1} '--\n']) [p,h,stats] = signrank(xy_sz{mm}(:,1),xy_sz{mm}(:,2),'method','approximate'); fprintf(['ALL VOXELS R2>TH - Wilcoxon HRF sz: z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); fprintf(['Mean:' num2str(nanmean(diff(xy_sz{mm},1,2))) ', std ' ... num2str(nanstd(diff(xy_sz{mm},1,2))) '\n']) [p,h,stats] = signrank(xy_ecc{mm}(:,1),xy_ecc{mm}(:,2),'method','approximate'); fprintf(['ALL VOXELS R2>TH - Wilcoxon HRF ecc: z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); fprintf(['Mean:' num2str(nanmean(diff(xy_ecc{mm},1,2))) ', std ' ... num2str(nanstd(diff(xy_ecc{mm},1,2))) '\n']) diffmat2{1}=[]; for r=1:length(roi) SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{rowmod,1})).mod.ROI,... ck_GetROIidx(roilabels(r),rois) ); [n,x] = hist(... T(modidx.(MRI_MODELS{mm,1})).mod.rfs(SSS) - ... T(modidx.(MRI_MODELS{mm,2})).mod.rfs(SSS), 100); f = n./sum(SSS); dsz = T(modidx.(MRI_MODELS{mm,1})).mod.rfs - ... T(modidx.(MRI_MODELS{mm,2})).mod.rfs; dsz = dsz(SSS); dsz = dsz(isfinite(dsz)); m = mean(dsz); sd = std(dsz); se = sd ./ sqrt(length(dsz)); diffmat2{1} = [diffmat2{1}; m sd se size(dsz,1)]; end end %% MRI ECC vs PRF Size ==================================================== Rth=5; f5=figure; set(f5,'Position',[100 100 800 1200]); set(f5,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap)); for m=1:length(MRI_MODELS)-1 s_R2 = T(modidx.(MRI_MODELS{m,1})).mod.R2 > Rth; EccBin = 0.5:1:16.5; subplot(3,2,(m-1)*2 +1);hold on; for r=1:length(roi) SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{m,1})).mod.ROI,... ck_GetROIidx(roilabels(r),rois) ); ES{r}=[]; scatter(T(modidx.(MRI_MODELS{m,1})).mod.ecc(SSS),... T(modidx.(MRI_MODELS{m,1})).mod.rfs(SSS),100,'Marker','.'); for b=1:length(EccBin) bb=[EccBin(b)-0.5 EccBin(b)+0.5]; PSZ=T(modidx.(MRI_MODELS{m,1})).mod.rfs(SSS); ECC=T(modidx.(MRI_MODELS{m,1})).mod.ecc(SSS); ES{r}=[ES{r}; EccBin(b) median(PSZ(ECC>=bb(1) & ECC<=bb(2)))]; end end title(['Ecc vs pRF size [' MMS{m} ', R>' num2str(Rth) ']'],... 'interpreter','none'); xlabel('Eccentricity');ylabel('pRF size'); set(gca, 'Box','off', 'xlim', [0 10], 'ylim',[0 10]); subplot(3,2,(m-1)*2 +2);hold on; for r=1:length(roi) h=plot(ES{r}(:,1),ES{r}(:,2),'o'); set(h,'MarkerSize',6,'markerfacecolor', get(h, 'color')); end title(['Ecc vs pRF size [' MMS{m} ', R>' num2str(Rth) ']'],... 'interpreter','none'); xlabel('Eccentricity');ylabel('pRF size'); set(gca, 'Box','off', 'xlim', [0 10], 'ylim',[0 10]); end if SaveFigs saveas(f5,fullfile(pngfld, 'MRI_Ecc_vs_Size.png')); saveas(f5,fullfile(svgfld, 'MRI_Ecc_vs_Size.svg')); end if CloseFigs; close(f5); end %% Fig 5 & SFig 5: Ecc-Sz for CSS per area including CI =================== SaveFigs=false; for Rth=[5 3] clear tbl lm; max_ecc = 100; bins=[0:18; 2:20]'; f5css=figure; set(f5css,'Position',[100 100 2000 1500]); set(f5css,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap)); sgtitle(['Voxels & fit, R>' num2str(Rth)]) m=3; % css s_R2 = T(modidx.(MRI_MODELS{m,1})).mod.R2 > Rth & ... T(modidx.(MRI_MODELS{m,1})).mod.ecc < max_ecc; nsp=length(roi); nrc=ceil(sqrt(nsp)); nvox_roi = []; binned_data=[]; for r=1:length(roi) ES{r} = []; ES_m{r} = []; ESb{r} = []; ESf{r} = []; SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{m,1})).mod.ROI,... ck_GetROIidx(roilabels(r),rois) ); s_Monkey = strcmp(T(modidx.(MRI_MODELS{m,1})).mod.Monkey,'M01'); SSS_m1 = s_R2 & ismember( T(modidx.(MRI_MODELS{m,1})).mod.ROI,... ck_GetROIidx(roilabels(r),rois)) & s_Monkey; s_Monkey = strcmp(T(modidx.(MRI_MODELS{m,1})).mod.Monkey,'M02'); SSS_m2 = s_R2 & ismember( T(modidx.(MRI_MODELS{m,1})).mod.ROI,... ck_GetROIidx(roilabels(r),rois)) & s_Monkey; nvox_roi = [nvox_roi; sum(SSS_m1) sum(SSS_m2)]; roi_n{r,1} = rois{r}; roi_n{r,2} = sum(SSS_m1); roi_n{r,3} = sum(SSS_m2); ES{r}=[T(modidx.(MRI_MODELS{m,1})).mod.ecc(SSS) ... T(modidx.(MRI_MODELS{m,1})).mod.rfs(SSS) ]; ES{r}( ES{r}(:,2) > 50,: ) = []; % remove insanely large pRF's ES_m{1,r}=[T(modidx.(MRI_MODELS{m,1})).mod.ecc(SSS_m1) ... T(modidx.(MRI_MODELS{m,1})).mod.rfs(SSS_m1) ]; ES_m{2,r}=[T(modidx.(MRI_MODELS{m,1})).mod.ecc(SSS_m2) ... T(modidx.(MRI_MODELS{m,1})).mod.rfs(SSS_m2) ]; %subplot(nrc,nrc,r); hold on; subplot(5,8,r); hold on; scatter(ES_m{1,r}(:,1),ES_m{1,r}(:,2),10,'Marker','o',... 'MarkerFaceColor','b','MarkerFaceAlpha',.2,... 'MarkerEdgeColor','none'); scatter(ES_m{2,r}(:,1),ES_m{2,r}(:,2),10,'Marker','o',... 'MarkerFaceColor','k','MarkerFaceAlpha',.2,... 'MarkerEdgeColor','none'); % make table tbl = table(ES{r}(:,1),ES{r}(:,2),'VariableNames',{'ecc','sz'}); if ~isempty(tbl) % remove inf's tbl(isinf(tbl.sz),:)=[]; % fit lm{r} = fitlm(tbl,'sz ~ 1 + ecc'); % CI lmCI{r}=coefCI(lm{r},0.05); if ~isnan(lm{r}.Coefficients.pValue(1)) % Prediction x=[0;20]; [ypred,yci] = predict(lm{r},x); yci2=yci(:,2)-ypred; if lm{r}.Coefficients.pValue(2) < 0.05 boundedline(x,ypred',yci2,'r','alpha') else boundedline(x,ypred',yci2,'k','alpha') end %fprintf([roilabels{r} ' n = ' num2str(size(tbl,1)) '\n']) %fprintf([roilabels{r} ' p = ' num2str(lm{r}.Coefficients.pValue(2)) '\n']) fprintf([roilabels{r} ' slope = ' num2str(lm{r}.Coefficients.Estimate(2)) '\n']) end % binned for b=1:size(bins,1) seldata = (tbl.ecc>bins(b,1) & tbl.ecc<=bins(b,2)) ; if sum(seldata)>1 binned_data=[binned_data;... r b mean(tbl.sz(seldata)) std(tbl.sz(seldata))./sqrt(sum(seldata))]; end end end title(['Ecc vs pRF size [' roilabels{r} ', R>' num2str(Rth) ']'],... 'interpreter','none'); title(roilabels{r}); xlabel('Eccentricity');ylabel('pRF size'); set(gca, 'Box','off', 'xlim', [0 20], 'ylim',[0 25]); end % binned f5css_binned=figure; set(f5css_binned,'Position',[100 100 2000 1500]); set(f5css_binned,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap)); sgtitle(['Binned & voxel-based fit, R>' num2str(Rth)]) for r=1:length(roi) subplot(5,8,r); hold on; errorbar(... binned_data(binned_data(:,1)==r,2),... binned_data(binned_data(:,1)==r,3),... binned_data(binned_data(:,1)==r,4),... 'LineStyle','none','Marker','o','MarkerSize',6,... 'Color','k','MarkerFaceColor',[.3 .3 .3]); plot([8 8],[0 25],'k--') if ~isempty(lm{r}) && ~isnan(lm{r}.Coefficients.pValue(1)) x=[0;20]; [ypred,yci] = predict(lm{r},x); yci2=yci(:,2)-ypred; if lm{r}.Coefficients.pValue(2) < 0.05 boundedline(x,ypred',yci2,'r','alpha') else boundedline(x,ypred',yci2,'k','alpha') end end title([roilabels{r} ', R>' num2str(Rth) ', nVox [' ... num2str(roi_n{r,2}) ' ' num2str(roi_n{r,3}) ']'],... 'interpreter','none'); %title(roilabels{r}); xlabel('Ecc (dva)');ylabel('Sz (dva)'); set(gca, 'Box','off', 'xlim', [0 20], 'ylim',[0 25]); end if SaveFigs saveas(f5css,fullfile(pngfld, ['MRI_Ecc_vs_Size_R' num2str(Rth) '.png'])); saveas(f5css,fullfile(svgfld, ['MRI_Ecc_vs_Size_R' num2str(Rth) '.svg'])); saveas(f5css_binned,fullfile(pngfld, ['MRI_Ecc_vs_Size_binned_R' num2str(Rth) '.png'])); saveas(f5css_binned,fullfile(svgfld, ['MRI_Ecc_vs_Size_binned_R' num2str(Rth) '.svg'])); end if CloseFigs; close(f5css_R5); close(f5css_binned_R5); end % Select areas with data in both monkeys ================================= if false minVox_both=25; % which ROIs have a minimum number of voxels in both animals? s_roi = mean(nvox_roi>=minVox_both,2)==1; fprintf(['ROIs with more than ' num2str(minVox_both) ' sign. voxels in BOTH:\n']) rois(s_roi,1) minVox_s1=minVox_both; s1_roi = nvox_roi(:,1)>minVox_s1; fprintf(['ROIs with more than ' num2str(minVox_s1) ' sign. voxels in S1:\n']) rois(s1_roi,1) minVox_s2=minVox_both; s2_roi = nvox_roi(:,2)>minVox_s2; fprintf(['ROIs with more than ' num2str(minVox_s2) ' sign. voxels in S2:\n']) rois(s2_roi,1) figure; hold on; lidx=1; clear L for idx=find(s_roi==1)' x=[0;20]; [ypred,yci] = predict(lm{idx},x); yci2=yci(:,2)-ypred; plot(x,ypred','LineWidth',3) L{lidx}=roilabels{idx}; lidx=lidx+1; end legend(L) end % NEW Fig 5: Ecc-Sz for CSS, selected ROIs =============================== A_idx = [1:5 8:9]; B_idx = [10:13 23 25:26]; C_idx = [19 18 20]; D_idx = [29:31 33 37:38]; % precreate color matrices for plotting A_cols = brewermap(length(A_idx),def_cmap); B_cols = brewermap(length(B_idx),def_cmap); C_cols = brewermap(length(C_idx),def_cmap); D_cols = brewermap(length(D_idx),def_cmap); f_eccsz = figure; set(f_eccsz,'Position',[100 100 2400 2000]); sgtitle(['Ecc vs Sz, R>' num2str(Rth)]) subplot(2,2,1); hold on; lidx=1; clear L for r=A_idx errorbar(... binned_data(binned_data(:,1)==r,2),... binned_data(binned_data(:,1)==r,3),... binned_data(binned_data(:,1)==r,4),... 'LineStyle','none','Marker','o','MarkerSize',6,... 'Color',A_cols(lidx,:),'MarkerFaceColor',A_cols(lidx,:)); lidx=lidx+1; end lidx=1; for r=A_idx x=[0;20]; [ypred,yci] = predict(lm{r},x); yci2=yci(:,2)-ypred; %boundedline(x,ypred',yci2,'cmap',A_cols(lidx,:),'alpha') plot(x,ypred','-','LineWidth',2,'Color',A_cols(lidx,:)) L{lidx}=roilabels{r}; lidx=lidx+1; end legend(L,'Location','NorthWest') plot([8 8],[0 25],'k--') set(gca,'ylim',[0 15], 'TickDir','out'); subplot(2,2,2); hold on; lidx=1; clear L for r=B_idx errorbar(... binned_data(binned_data(:,1)==r,2),... binned_data(binned_data(:,1)==r,3),... binned_data(binned_data(:,1)==r,4),... 'LineStyle','none','Marker','o','MarkerSize',6,... 'Color',B_cols(lidx,:),'MarkerFaceColor',B_cols(lidx,:)); lidx=lidx+1; end lidx=1; for r=B_idx x=[0;20]; [ypred,yci] = predict(lm{r},x); yci2=yci(:,2)-ypred; %boundedline(x,ypred',yci2,'cmap',B_cols(lidx,:),'alpha') plot(x,ypred','-','LineWidth',2,'Color',B_cols(lidx,:)) L{lidx}=roilabels{r}; lidx=lidx+1; end legend(L,'Location','NorthWest') plot([8 8],[0 25],'k--') set(gca,'ylim',[0 20], 'TickDir','out'); C_idx = [19 18 20]; D_idx = [29:31 33 37 38]; subplot(2,2,3); hold on; lidx=1; clear L for r=C_idx errorbar(... binned_data(binned_data(:,1)==r,2),... binned_data(binned_data(:,1)==r,3),... binned_data(binned_data(:,1)==r,4),... 'LineStyle','none','Marker','o','MarkerSize',6,... 'Color',C_cols(lidx,:),'MarkerFaceColor',C_cols(lidx,:)); lidx=lidx+1; end lidx=1; for r=C_idx x=[0;20]; [ypred,yci] = predict(lm{r},x); yci2=yci(:,2)-ypred; %boundedline(x,ypred',yci2,'cmap',C_cols(lidx,:),'alpha') plot(x,ypred','-','LineWidth',2,'Color',C_cols(lidx,:)) L{lidx}=roilabels{r}; lidx=lidx+1; end legend(L,'Location','NorthWest') plot([8 8],[0 25],'k--') set(gca,'ylim',[0 20], 'TickDir','out'); subplot(2,2,4); hold on; lidx=1; clear L for r=D_idx errorbar(... binned_data(binned_data(:,1)==r,2),... binned_data(binned_data(:,1)==r,3),... binned_data(binned_data(:,1)==r,4),... 'LineStyle','none','Marker','o','MarkerSize',6,... 'Color',D_cols(lidx,:),'MarkerFaceColor',D_cols(lidx,:)); lidx=lidx+1; end lidx=1; for r=D_idx x=[0;20]; [ypred,yci] = predict(lm{r},x); yci2=yci(:,2)-ypred; %boundedline(x,ypred',yci2,'cmap',D_cols(lidx,:),'alpha') plot(x,ypred','-','LineWidth',2,'Color',D_cols(lidx,:)) L{lidx}=roilabels{r}; lidx=lidx+1; end legend(L,'Location','NorthWest') plot([8 8],[0 25],'k--') set(gca,'ylim',[0 20], 'TickDir','out'); if SaveFigs saveas(f_eccsz,fullfile(pngfld, ['MRI_Ecc_vs_Size_ROIs_R' num2str(Rth) '.png'])); saveas(f_eccsz,fullfile(svgfld, ['MRI_Ecc_vs_Size_ROIs_R' num2str(Rth) '.svg'])); end if CloseFigs; close(f_eccsz); end end %% SFig 5A/B: Ecc-Sz per area including CI and binned data with SEM ======= reroi = [19 18 20 1:5 7:14 21:38]; for Rth=[5 3] clear tbl lm; bins=[0:18; 2:20]'; sf5css=figure; set(sf5css,'Position',[100 100 1600 1200]); set(sf5css,'DefaultAxesColorOrder',brewermap(length(roi),def_cmap)); sgtitle(['R>' num2str(Rth)]) m=3; % css s_R2 = T(modidx.(MRI_MODELS{m,1})).mod.R2 > Rth; nsp=length(reroi); nvox_roi = []; binned_data=[]; for r=1:length(reroi) rr=reroi(r); ES2{r} = []; ES2_m{r} = []; ES2b{r} = []; ES2f{r} = []; SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{m,1})).mod.ROI,... ck_GetROIidx(roilabels(rr),rois) ); s_Monkey = strcmp(T(modidx.(MRI_MODELS{m,1})).mod.Monkey,'M01'); SSS_m1 = s_R2 & ismember( T(modidx.(MRI_MODELS{m,1})).mod.ROI,... ck_GetROIidx(roilabels(rr),rois)) & s_Monkey; s_Monkey = strcmp(T(modidx.(MRI_MODELS{m,1})).mod.Monkey,'M02'); SSS_m2 = s_R2 & ismember( T(modidx.(MRI_MODELS{m,1})).mod.ROI,... ck_GetROIidx(roilabels(rr),rois)) & s_Monkey; nvox_roi = [nvox_roi; sum(SSS_m1) sum(SSS_m2)]; roi_n{r,1} = rois{rr}; roi_n{r,2} = sum(SSS_m1); roi_n{r,3} = sum(SSS_m2); ES2{r}=[T(modidx.(MRI_MODELS{m,1})).mod.ecc(SSS) ... T(modidx.(MRI_MODELS{m,1})).mod.rfs(SSS) ]; ES2{r}( ES2{r}(:,2) > 50,: ) = []; % remove insanely large pRF's ES2_m{1,r}=[T(modidx.(MRI_MODELS{m,1})).mod.ecc(SSS_m1) ... T(modidx.(MRI_MODELS{m,1})).mod.rfs(SSS_m1) ]; ES2_m{2,r}=[T(modidx.(MRI_MODELS{m,1})).mod.ecc(SSS_m2) ... T(modidx.(MRI_MODELS{m,1})).mod.rfs(SSS_m2) ]; subplot(5,7,r); hold on; scatter(ES2_m{1,r}(:,1),ES2_m{1,r}(:,2),20,'Marker','o',... 'MarkerFaceColor','b','MarkerFaceAlpha',.2,... 'MarkerEdgeColor','none'); scatter(ES2_m{2,r}(:,1),ES2_m{2,r}(:,2),20,'Marker','o',... 'MarkerFaceColor','k','MarkerFaceAlpha',.2,... 'MarkerEdgeColor','none'); % make table tbl2 = table(ES2{r}(:,1),ES2{r}(:,2),'VariableNames',{'ecc','sz'}); if ~isempty(tbl2) % remove inf's tbl2(isinf(tbl2.sz),:)=[]; % fit lm2{r} = fitlm(tbl2,'sz ~ 1 + ecc'); % CI lmCI2{r}=coefCI(lm2{r},0.05); if ~isnan(lm2{r}.Coefficients.pValue(1)) % Prediction x=[0;20]; [ypred,yci] = predict(lm2{r},x); yci2=yci(:,2)-ypred; if lm2{r}.Coefficients.pValue(2) < 0.05 boundedline(x,ypred',yci2,'r','alpha') else boundedline(x,ypred',yci2,'k','alpha') end fprintf([roilabels{rr} ' n = ' num2str(size(tbl2,1)) '\n']) fprintf([roilabels{rr} ' p = ' num2str(lm2{r}.Coefficients.pValue(2)) '\n']) fprintf([roilabels{rr} ' slope = ' num2str(lm2{r}.Coefficients.Estimate(2)) '\n']) end % binned for b=1:size(bins,1) seldata = (tbl2.ecc>bins(b,1) & tbl2.ecc<=bins(b,2)) ; if sum(seldata)>1 binned_data=[binned_data;... rr b mean(tbl2.sz(seldata)) std(tbl2.sz(seldata))./sqrt(sum(seldata))]; end end end errorbar(... binned_data(binned_data(:,1)==rr,2),... binned_data(binned_data(:,1)==rr,3),... binned_data(binned_data(:,1)==rr,4),... 'LineStyle','none','Marker','o','MarkerSize',6,... 'Color','k','MarkerFaceColor',[.8 .8 .8]); plot([8 8],[0 25],'k--') %xlabel('Eccentricity');ylabel('pRF size'); set(gca, 'Box','off', 'xlim', [0 20], 'ylim',[0 25],... 'FontName','Myriad Pro','FontSize',11,'TitleHorizontalAlignment','left'); title(['Ecc vs pRF size [' roilabels{rr} ', R>' num2str(Rth) ']'],... 'interpreter','none'); title(roilabels{rr}); text(2,20,sprintf(['n = [' num2str(roi_n{r,2}) ' ' num2str(roi_n{r,3}) ']\n'... 'p = ' num2str(lm2{r}.Coefficients.pValue(2)) '\n'... 'slope = ' num2str(lm2{r}.Coefficients.Estimate(2))])); end if SaveFigs saveas(sf5css,fullfile(pngfld, ['SF5_MRI_Ecc_vs_Size_R' num2str(Rth) '.png'])); saveas(sf5css,fullfile(svgfld, ['SF5_MRI_Ecc_vs_Size_R' num2str(Rth) '.svg'])); end if CloseFigs; close(sf5css); end end %% ======================================================================== % EPHYS PRF ANALYSIS ----------------------------------------------------- % ======================================================================== %% Fig 9: EPHYS ECC vs PRF Size =========================================== Rth=25; SNRth=3; m=unique(tMUA.Model); R2=[]; for i=2%1:length(m) R2 = [R2 tMUA.R2(strcmp(tMUA.Model,m{i}))]; end subs = tMUA.Monkey; for i=1:length(subs) if strcmp(subs{i},'M03') subs{i} = 1; elseif strcmp(subs{i},'M04') subs{i} = 2; end end subs=cell2mat(subs); names={'M03','M04'}; mm=unique(tMUA.Model); for m=2%1:length(mm) % collect if ~strcmp(mm{m},'classicRF') eccsz = [tMUA.R2(strcmp(tMUA.Model,mm{m})) ... tMUA.ecc(strcmp(tMUA.Model,mm{m})) ... tMUA.rfs(strcmp(tMUA.Model,mm{m})) ... tMUA.Area(strcmp(tMUA.Model,mm{m})) ... tMUA.Array(strcmp(tMUA.Model,mm{m})) ... subs(strcmp(tMUA.Model,mm{m})) ... tMUA.gain(strcmp(tMUA.Model,mm{m})) ... ]; else eccsz = [tMUA.SNR(strcmp(tMUA.Model,mm{m})) ... tMUA.ecc(strcmp(tMUA.Model,mm{m})) ... tMUA.rfs(strcmp(tMUA.Model,mm{m})) ... tMUA.Area(strcmp(tMUA.Model,mm{m})) ... tMUA.Array(strcmp(tMUA.Model,mm{m})) ... subs(strcmp(tMUA.Model,mm{m})) ... ]; end f_eccsz_arr = figure; set(f_eccsz_arr,'Position',[100 100 1800 600]); for ss=1:2 a=[1 4]; for aa=1:length(a) leg={}; if strcmp(mm{m},'classicRF') sel = eccsz(:,1)>SNRth & eccsz(:,4)==a(aa) & eccsz(:,6)==ss; else sel = eccsz(:,1)>Rth & eccsz(:,4)==a(aa) & eccsz(:,6)==ss; end es1=eccsz(sel,:); subplot(2,2,(ss-1)*2 + aa); hold on; arr=unique(es1(:,5)); for aaa = 1:length(arr) if arr(aaa)~=99 sel2=es1(:,5)==arr(aaa); scatter(es1(sel2,2),es1(sel2,3)); leg=[leg {num2str(arr(aaa))}]; end end if a(aa)==1 set(gca,'xlim',[0 8],'ylim',[0 3]); else set(gca,'xlim',[0 8],'ylim',[0 5]); end legend(leg,'Location','NorthEastOutside'); title([names{ss} ' V' num2str(a(aa))]); xlabel('Ecc (dva)'); ylabel('RF-size (sd, dva)') end end sgtitle(['MODEL: ' mm{m} ', R2th: ' num2str(Rth)],... 'interpreter','none'); if SaveFigs saveas(f_eccsz_arr,fullfile(svgfld, ... ['EPHYS_MUA_Ecc_vs_Size_' mm{m} '.svg'])); saveas(f_eccsz_arr,fullfile(pngfld, ... ['EPHYS_MUA_Ecc_vs_Size_' mm{m} '.png'])); end if CloseFigs; close(f_eccsz_arr); end f_eccsz_v1 = figure; set(f_eccsz_v1,'Position',[100 100 1600 400]); for ss=1:2 leg={}; if strcmp(mm{m},'classicRF') sel = eccsz(:,1)>SNRth & eccsz(:,4)==1 & eccsz(:,6)==ss; else sel = eccsz(:,1)>Rth & eccsz(:,4)==1 & eccsz(:,6)==ss; end es1=eccsz(sel,:); subplot(1,2,ss); hold on; arr=unique(es1(:,5)); for aaa = 1:length(arr) sel2=es1(:,5)==arr(aaa); scat=scatter(es1(sel2,2),es1(sel2,3)); %if (ss==1 && (arr(aaa)==11 || arr(aaa)==13)) || ... % (ss==2 && (arr(aaa)==10 || arr(aaa)==12)) if (ss==1 && arr(aaa)==11) || ... (ss==2 && arr(aaa)==10) set(scat,'MarkerFaceColor','r','MarkerEdgeColor','r'); elseif (ss==1 && arr(aaa)==13) || ... (ss==2 && arr(aaa)==12) set(scat,'MarkerFaceColor','b','MarkerEdgeColor','b'); else set(scat,'MarkerFaceColor','k','MarkerEdgeColor','k'); end leg=[leg {num2str(arr(aaa))}]; end title([names{ss} ' V1']); xlabel('Ecc (dva)'); ylabel('RF-size (sd, dva)') set(gca,'xlim',[0 8],'ylim',[0 3]); %legend(leg,'Location','NorthEastOutside'); end sgtitle(['MODEL: ' mm{m} ', R2th: ' num2str(Rth)],... 'interpreter','none'); if SaveFigs saveas(f_eccsz_v1,fullfile(pngfld, ... ['EPHYS_MUA_Ecc_vs_Size_' mm{m} '.png'])); saveas(f_eccsz_v1,fullfile(svgfld, ... ['EPHYS_MUA_Ecc_vs_Size_' mm{m} '.svg'])); end if CloseFigs; close(f_eccsz_v1); end end %close all; %% Fig 6 / SFig 6 EPHYS VFC MUA (scatter and heatmap) ===================== % scatter MUA ============================================================= R2TH = 50; s=tMUA.R2>R2TH; model='css_ephys_cv1'; fm=figure; set(fm,'Position',[100 100 1400 1000]); [cmap,~,~] = brewermap(length(unique(tMUA.Array)),'spectral'); set(fm,'DefaultAxesColorOrder',cmap); msz = 5; % M03 V1 ---- m=strcmp(tMUA.Monkey,'M03') & strcmp(tMUA.Model,model); v=tMUA.Area==1; subplot(2,2,1);hold on; plot([-10 20],[0 0],'k'); plot([0 0],[-30 10],'k'); lt={'meridian','meridian'}; nelec=0; for r=unique(tMUA.Array)' a=tMUA.Array==r; currcol = cmap(r,:); if sum(s & m & a & v) > 0 nelec=nelec+sum(s & m & a & v); p=plot(tMUA.X(s & m & a & v),... tMUA.Y(s & m & a & v),'o','Color',currcol,... 'LineStyle','none','LineWidth',1,... 'MarkerSize',msz,'MarkerFaceColor',currcol); lt=[lt num2str(r)]; end end set(gca, 'Box','off', 'xlim', [-1 8], 'ylim',[-6 1]); title(['M03 V1 MUA (n = ' num2str(nelec) ')']) l=legend(lt); set(l,'Location','NorthEastOutside'); clear lt; % M03 V4 ---- v=tMUA.Area==4; subplot(2,2,3);hold on; plot([-10 20],[0 0],'k'); plot([0 0],[-30 10],'k'); lt={'meridian','meridian'}; nelec=0; for r=unique(tMUA.Array)' a=tMUA.Array==r; currcol = cmap(r,:); if sum(s & m & a & v) > 0 nelec=nelec+sum(s & m & a & v); p=plot(tMUA.X(s & m & a & v),... tMUA.Y(s & m & a & v),'o','Color',currcol,... 'LineStyle','none','LineWidth',1,... 'MarkerSize',msz,'MarkerFaceColor',currcol); lt=[lt num2str(r)]; end end set(gca, 'Box','off', 'xlim', [-1 8], 'ylim',[-6 1]); % set(gca, 'Box','off', 'xlim', [-2 25], 'ylim',[-25 2]); title(['M03 V4 MUA (n = ' num2str(nelec) ')']) l=legend(lt); set(l,'Location','NorthEastOutside'); clear lt; % M04 V1 ---- m=strcmp(tMUA.Monkey,'M04') & strcmp(tMUA.Model,model); v=tMUA.Area==1; subplot(2,2,2);hold on; plot([-10 20],[0 0],'k'); plot([0 0],[-30 10],'k'); lt={'meridian','meridian'}; nelec=0; for r=unique(tMUA.Array)' a=tMUA.Array==r; currcol = cmap(r,:); if sum(s & m & a & v) > 0 nelec=nelec+sum(s & m & a & v); p=plot(tMUA.X(s & m & a & v),... tMUA.Y(s & m & a & v),'o','Color',currcol,... 'LineStyle','none','LineWidth',1,... 'MarkerSize',msz,'MarkerFaceColor',currcol); lt=[lt num2str(r)]; end end set(gca, 'Box','off', 'xlim', [-1 8], 'ylim',[-6 1]); title(['M04 V1 MUA (n = ' num2str(nelec) ')']) l=legend(lt); set(l,'Location','NorthEastOutside'); clear lt; % M04 V4 ---- v=tMUA.Area==4; subplot(2,2,4);hold on; plot([-10 30],[0 0],'k'); plot([0 0],[-30 10],'k'); lt={'meridian','meridian'}; nelec=0; for r=unique(tMUA.Array)' a=tMUA.Array==r; currcol = cmap(r,:); if sum(s & m & a & v) > 0 nelec=nelec+sum(s & m & a & v); p=plot(tMUA.X(s & m & a & v),... tMUA.Y(s & m & a & v),'o','Color',currcol,... 'LineStyle','none','LineWidth',1,... 'MarkerSize',msz,'MarkerFaceColor',currcol); lt=[lt num2str(r)]; end end set(gca, 'Box','off', 'xlim', [-2 25], 'ylim',[-25 2]); title(['M04 V4 MUA (n = ' num2str(nelec) ')']) l=legend(lt); set(l,'Location','NorthEastOutside'); clear lt; if SaveFigs saveas(fm,fullfile(pngfld, 'MUA_VFC.png')); saveas(fm,fullfile(svgfld, 'MUA_VFC.svg')); end if CloseFigs; close(fm); end % Heatmap MUA ============================================================= fhm=figure; set(fhm,'Position',[100 100 1800 600]); settings.PixPerDeg = 29.5032; settings.meshsize = 2000; colormap(inferno) % M03 V1 ---- m=strcmp(tMUA.Monkey,'M03') & strcmp(tMUA.Model,model); v=tMUA.Area==1; subplot(2,4,1);hold on; allprf=[]; for r=unique(tMUA.Array)' a=tMUA.Array==r; if sum(s & m & a & v) > 0 allprf = [allprf; ... tMUA.X(s & m & a & v) ... tMUA.Y(s & m & a & v) ... tMUA.rfs(s & m & a & v)]; end end res = ck_2dPRF_ephys(allprf(:,1),allprf(:,2),allprf(:,3)); img=flipud(res.img); res.ymesh2 = fliplr(res.ymesh); % zoom in on plot xrange = [-3 8]; yrange = [-6 3]; xrange_idx = [... find(res.xmesh >= xrange(1),1,'first') find(res.xmesh <= xrange(2),1,'last')]; yrange_idx = [... find(res.ymesh2 >= yrange(1),1,'first') find(res.ymesh2 <= yrange(2),1,'last')]; xrr=res.xmesh(xrange_idx); yrr=res.ymesh2(yrange_idx); % plot sumimg=sum(img,3); imagesc(sumimg); set(gca,'xlim',xrange_idx,'ylim',... yrange_idx,'Color','k','xtick',[],'ytick',[]) colorbar; subplot(2,4,5); hold on; plot(1.2*res.xr,[0 0],'w'); plot([0 0],1.2*res.yr,'w'); set(gca,'xlim',xrr,'ylim',yrr,'Color','k') colorbar; clear('res','img','sumimg'); title('M03 V1 MUA') % M03 V4 ---- v=tMUA.Area==4; subplot(2,4,2);hold on; allprf=[]; for r=unique(tMUA.Array)' a=tMUA.Array==r; if sum(s & m & a & v) > 0 allprf = [allprf; ... tMUA.X(s & m & a & v) ... tMUA.Y(s & m & a & v) ... tMUA.rfs(s & m & a & v)]; end end res = ck_2dPRF_ephys(allprf(:,1),allprf(:,2),allprf(:,3)); img=flipud(res.img); res.ymesh2 = fliplr(res.ymesh); % zoom in on plot xrange = [-3 8]; yrange = [-6 3]; xrange_idx = [... find(res.xmesh >= xrange(1),1,'first') find(res.xmesh <= xrange(2),1,'last')]; yrange_idx = [... find(res.ymesh2 >= yrange(1),1,'first') find(res.ymesh2 <= yrange(2),1,'last')]; xrr=res.xmesh(xrange_idx); yrr=res.ymesh2(yrange_idx); % plot sumimg=sum(img,3); imagesc(sumimg); set(gca,'xlim',xrange_idx,'ylim',... yrange_idx,'Color','k','xtick',[],'ytick',[]) colorbar; subplot(2,4,6); hold on; plot(1.2*res.xr,[0 0],'w'); plot([0 0],1.2*res.yr,'w'); set(gca,'xlim',xrr,'ylim',yrr,'Color','k') colorbar; clear('res','img','sumimg'); title('M03 V4 MUA') % M04 V1 ---- m=strcmp(tMUA.Monkey,'M04') & strcmp(tMUA.Model,model); v=tMUA.Area==1; subplot(2,4,3);hold on; allprf=[]; for r=unique(tMUA.Array)' a=tMUA.Array==r; if sum(s & m & a & v) > 0 allprf = [allprf; ... tMUA.X(s & m & a & v) ... tMUA.Y(s & m & a & v) ... tMUA.rfs(s & m & a & v)]; end end res = ck_2dPRF_ephys(allprf(:,1),allprf(:,2),allprf(:,3)); img=flipud(res.img); res.ymesh2 = fliplr(res.ymesh); % zoom in on plot xrange = [-3 8]; yrange = [-6 3]; xrange_idx = [... find(res.xmesh >= xrange(1),1,'first') find(res.xmesh <= xrange(2),1,'last')]; yrange_idx = [... find(res.ymesh2 >= yrange(1),1,'first') find(res.ymesh2 <= yrange(2),1,'last')]; xrr=res.xmesh(xrange_idx); yrr=res.ymesh2(yrange_idx); % plot sumimg=sum(img,3); imagesc(sumimg); set(gca,'xlim',xrange_idx,'ylim',... yrange_idx,'Color','k','xtick',[],'ytick',[]) colorbar; subplot(2,4,7); hold on; plot(1.2*res.xr,[0 0],'w'); plot([0 0],1.2*res.yr,'w'); set(gca,'xlim',xrr,'ylim',yrr,'Color','k') colorbar; clear('res','img','sumimg'); title('M04 V1 MUA') % M04 V4 ---- v=tMUA.Area==4; subplot(2,4,4);hold on; allprf=[]; for r=unique(tMUA.Array)' a=tMUA.Array==r; if sum(s & m & a & v) > 0 allprf = [allprf; ... tMUA.X(s & m & a & v) ... tMUA.Y(s & m & a & v) ... tMUA.rfs(s & m & a & v)]; end end res = ck_2dPRF_ephys(allprf(:,1),allprf(:,2),allprf(:,3)); img=flipud(res.img); res.ymesh2 = fliplr(res.ymesh); % zoom in on plot xrange = [-5 25]; yrange = [-25 5]; xrange_idx = [... find(res.xmesh >= xrange(1),1,'first') find(res.xmesh <= xrange(2),1,'last')]; yrange_idx = [... find(res.ymesh2 >= yrange(1),1,'first') find(res.ymesh2 <= yrange(2),1,'last')]; xrr=res.xmesh(xrange_idx); yrr=res.ymesh2(yrange_idx); % plot sumimg=sum(img,3); imagesc(sumimg); set(gca,'xlim',xrange_idx,'ylim',... yrange_idx,'Color','k','xtick',[],'ytick',[]) colorbar; subplot(2,4,8); hold on; plot(1.2*res.xr,[0 0],'w'); plot([0 0],1.2*res.yr,'w'); set(gca,'xlim',xrr,'ylim',yrr,'Color','k') colorbar; clear('res','img','sumimg'); title('M04 V4 MUA') if SaveFigs saveas(fhm,fullfile(pngfld, 'MUA_VFC_hm.png')); saveas(fhm,fullfile(svgfld, 'MUA_VFC_hm.svg')); end if CloseFigs; close(fhm); end %% Ephys location difference & size difference ============================ rth=25; snrth=3; c_rth=25; ephys_MMS = MMS(:,1); clear C R2m SZ for m=1:length(ephys_MOD) C{m}=[];R2m{m}=[];SZ{m}=[]; model=ephys_MOD{m}; s = strcmp(tMUA.Model,model); C{m}=[C{m} tMUA.R2(s) tMUA.X(s) tMUA.Y(s) tMUA.rfs(s)]; R2m{m}=[R2m{m} tMUA.R2(s)]; SZ{m}=[SZ{m} tMUA.R2(s) tMUA.rfs(s)]; PRF_EST(m,1).sig = 'MUA'; PRF_EST(m,1).R2 = tMUA.R2(s); PRF_EST(m,1).X = tMUA.X(s); PRF_EST(m,1).Y = tMUA.Y(s); PRF_EST(m,1).S = tMUA.rfs(s); PRF_EST(m,1).A = tMUA.Area(s); PRF_EST(m,1).G = tMUA.gain(s); PRF_EST(m,1).M = tMUA.Monkey(s); PRF_EST(m,1).ARR = tMUA.Array(s); s = strcmp(tMUA.Model,'classicRF'); %C{m}=[C{m} tMUA.X(s)./668.745 tMUA.Y(s)./668.745 tMUA.rfs(s)./2]; PRF_EST(m,2).sig = 'MUACLASSIC'; PRF_EST(m,2).R2 = tMUA.R2(s); PRF_EST(m,2).X = tMUA.X(s); PRF_EST(m,2).Y = tMUA.Y(s); PRF_EST(m,2).S = tMUA.rfs(s); PRF_EST(m,2).A = tMUA.Area(s); PRF_EST(m,2).G = tMUA.gain(s); PRF_EST(m,2).M = tMUA.Monkey(s); PRF_EST(m,2).SNR = tMUA.SNR(s); PRF_EST(m,2).ARR = tMUA.Array(s); s = strcmp(tLFP.Model,model); sig=unique(tLFP.SigType); lfp_order = [3 1 2 5 4]; cnt=1; for i=lfp_order b=strcmp(tLFP.SigType,sig{i}); C{m}=[C{m} tLFP.R2(s & b) tLFP.X(s & b) tLFP.Y(s & b) tLFP.rfs(s & b)]; R2m{m}=[R2m{m} tLFP.R2(s & b)]; SZ{m}=[SZ{m} tLFP.R2(s & b) tLFP.rfs(s & b)]; PRF_EST(m,2+cnt).sig = sig{i}; PRF_EST(m,2+cnt).R2 = tLFP.R2(s & b); PRF_EST(m,2+cnt).X = tLFP.X(s & b); PRF_EST(m,2+cnt).Y = tLFP.Y(s & b); PRF_EST(m,2+cnt).S = tLFP.rfs(s & b); PRF_EST(m,2+cnt).A = tLFP.Area(s & b); PRF_EST(m,2+cnt).G = tLFP.gain(s & b); PRF_EST(m,2+cnt).M = tLFP.Monkey(s & b); cnt=cnt+1; end s= sum(R2m{m}>rth,2)==size(R2m{m},2); distRF{m} = [... sqrt(((C{m}(s,2)-C{m}(s,5)).^2) + ((C{m}(s,3)-C{m}(s,6)).^2)) ... sqrt(((C{m}(s,2)-C{m}(s,9)).^2) + ((C{m}(s,3)-C{m}(s,10)).^2)) ... sqrt(((C{m}(s,2)-C{m}(s,13)).^2) + ((C{m}(s,3)-C{m}(s,14)).^2)) ... sqrt(((C{m}(s,2)-C{m}(s,17)).^2) + ((C{m}(s,3)-C{m}(s,18)).^2)) ... sqrt(((C{m}(s,2)-C{m}(s,21)).^2) + ((C{m}(s,3)-C{m}(s,22)).^2)) ]; distSZ{m} = [... C{m}(s,4)-C{m}(s,7) ... C{m}(s,4)-C{m}(s,11) ... C{m}(s,4)-C{m}(s,15) ... C{m}(s,4)-C{m}(s,19) ... C{m}(s,4)-C{m}(s,23) ]; % normalize by MUA-pRF normSz{m} = [... C{m}(s,7)./C{m}(s,4) ... C{m}(s,11)./C{m}(s,4) ... C{m}(s,15)./C{m}(s,4) ... C{m}(s,19)./C{m}(s,4) ... C{m}(s,23)./C{m}(s,4) ]; % % normalize by cRF % normSz{m} = [... % C{m}(s,7)./PRF_EST(m,2).S(s,:) ... % C{m}(s,11)./PRF_EST(m,2).S(s,:) ... % C{m}(s,15)./PRF_EST(m,2).S(s,:) ... % C{m}(s,19)./PRF_EST(m,2).S(s,:) ... % C{m}(s,23)./PRF_EST(m,2).S(s,:) ]; end %% Fig 7: compare MUA with classic ======================================== for m=[1 3] mlabel=ephys_MMS{m}; cRF_coll = []; pRF_coll = []; f_cRFpRF=figure; set(f_cRFpRF,'Position',[100 100 1500 600]); pn=0; for area = [1 4] pn=pn+1; fprintf(['=== AREA V' num2str(area) ' ===\n']) % location sel = PRF_EST(m,1).A == area & ... PRF_EST(m,1).R2 > rth & ... PRF_EST(m,2).R2 >= c_rth & ... PRF_EST(m,2).SNR >= snrth & ... ( strcmp(PRF_EST(m,1).M,'M03') & (PRF_EST(m,1).ARR ~= 11 & PRF_EST(m,1).ARR ~= 13) | ... strcmp(PRF_EST(m,1).M,'M04') & (PRF_EST(m,1).ARR ~= 10 & PRF_EST(m,1).ARR ~= 12)); dLOCATION = sqrt(... (PRF_EST(m,1).X(sel)-PRF_EST(m,2).X(sel)).^2 + ... (PRF_EST(m,1).Y(sel)-PRF_EST(m,2).Y(sel)).^2); mselect = PRF_EST(m,1).M(sel); m3idx = strcmp(mselect,'M03'); m4idx = strcmp(mselect,'M04'); fprintf(['MODEL ' ephys_MOD{m} ', MUA vs CLASSIC distance ----\n']) fprintf(['Mean ' num2str(nanmean(dLOCATION)) ', STD ' num2str(nanstd(dLOCATION)) '\n']) fprintf(['Median ' num2str(nanmedian(dLOCATION)) ', IQR ' num2str(iqr(dLOCATION)) '\n']) % size SIZE_MUA = [PRF_EST(m,1).S(sel) PRF_EST(m,2).S(sel)]; rmINF = isinf(SIZE_MUA(:,1)); SIZE_MUA(rmINF,:)=[]; m3idx(rmINF,:)=[]; m4idx(rmINF,:)=[]; [p,h,stats] = signrank(SIZE_MUA(:,1),SIZE_MUA(:,2)); fprintf(['MODEL ' ephys_MOD{m} ', MUA vs CLASSIC size ----\n']) fprintf(['dSZ: Wilcoxon z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); fprintf(['Mean (mua-classic) ' num2str(nanmean(SIZE_MUA(:,1)-SIZE_MUA(:,2))) ... ', STD ' num2str(nanstd(SIZE_MUA(:,1)-SIZE_MUA(:,2))) '\n']) fprintf(['Median ' num2str(nanmedian(SIZE_MUA(:,1)-SIZE_MUA(:,2))) ... ', IQR ' num2str(iqr(SIZE_MUA(:,1)-SIZE_MUA(:,2))) '\n']) subplot(1,2,pn);hold on; splim=[3 8]; plot([0 10],[0 10],'-r'); % split by monkey pRF3=SIZE_MUA(m3idx,1); cRF3=SIZE_MUA(m3idx,2); scatter(pRF3,cRF3,50,'Marker','o',... 'MarkerEdgeColor','none',... 'MarkerFaceColor',[.3 .3 .3],'MarkerFaceAlpha',0.3); idx = ~isnan(pRF3) & ~isnan(cRF3); polyfit(pRF3(idx),cRF3(idx),1) pRF4=SIZE_MUA(m4idx,1); cRF4=SIZE_MUA(m4idx,2); scatter(pRF4,cRF4,50,'Marker','o',... 'MarkerEdgeColor','none',... 'MarkerFaceColor',[.0 .0 .5],'MarkerFaceAlpha',0.3); idx = ~isnan(pRF4) & ~isnan(cRF4); polyfit(pRF4(idx),cRF4(idx),1) pRF=[pRF3;pRF4];cRF=[cRF3;cRF4]; n_electrodes = [length(pRF3) length(pRF4)] plot([0 10],[0 10],'-r'); set(gca, 'xlim',[0 splim(pn)],'ylim',[0 splim(pn)]); legend({'x=y','M3','M4'}) title(['Area V' num2str(area)]) xlabel('MUA pRF size'); ylabel('Classic RF size') % Stats per area pRF_both = [pRF3;pRF4]; cRF_both = [cRF3;cRF4]; [R,p] = corr(cRF_both,pRF_both,'rows','complete'); fprintf(['Pearson corr for V' num2str(area) ': R=' num2str(R) ', p=' num2str(p) '\n']) cRF_coll = [cRF_coll; ... cRF3 3*ones(size(cRF3)) area*ones(size(cRF3));... cRF4 4*ones(size(cRF4)) area*ones(size(cRF4))]; pRF_coll = [pRF_coll; ... pRF3 3*ones(size(cRF3)) area*ones(size(pRF3));... pRF4 4*ones(size(cRF4)) area*ones(size(pRF4))]; end [R,p] = corr(cRF_coll(:,1),pRF_coll(:,1),'rows','complete'); fprintf(['Pearson corr: R=' num2str(R) ', p=' num2str(p) '\n']) %% f_cRFpRF2 = figure; set(f_cRFpRF2,'Position',[0 0 2000 400]); subplot(1,3,1); hold on plot([0 10],[0 10],'-r'); % M3 V1 s = pRF_coll(:,2) == 3 & pRF_coll(:,3) == 1; scatter(pRF_coll(s,1),cRF_coll(s,1),55,'Marker','o',... 'MarkerEdgeColor','none',... 'MarkerFaceColor','k','MarkerFaceAlpha',0.3); % M4 V1 s = pRF_coll(:,2) == 4 & pRF_coll(:,3) == 1; scatter(pRF_coll(s,1),cRF_coll(s,1),55,'Marker','o',... 'MarkerEdgeColor','none',... 'MarkerFaceColor',[.0 .0 .5],'MarkerFaceAlpha',0.3); % M3 V4 s = pRF_coll(:,2) == 3 & pRF_coll(:,3) == 4; scatter(pRF_coll(s,1),cRF_coll(s,1),50,'Marker','o',... 'MarkerFaceColor',[1 1 1],'LineWidth',1,... 'MarkerEdgeColor','k','MarkerFaceAlpha',0.3); % M4 V4 s = pRF_coll(:,2) == 4 & pRF_coll(:,3) == 4; scatter(pRF_coll(s,1),cRF_coll(s,1),50,'Marker','o',... 'MarkerFaceColor',[1 1 1],'LineWidth',1,... 'MarkerEdgeColor',[.0 .0 .5],'MarkerFaceAlpha',0.3); set(gca, 'xlim',[0 splim(pn)],'ylim',[0 splim(pn)]); legend({'x=y','M3 V1','M4 V1','M3 V4','M4 V4'}) xlabel('MUA pRF size'); ylabel('Classic RF size') ppp=pRF_coll(:,1);ccc=cRF_coll(:,1); title(['pRF model: ' mlabel]); subplot(1,3,2); hold on; histogram(ppp-ccc,-5.05:0.1:5.05,'FaceColor',[.5 .5 .5],'LineStyle','none') xlabel('pRF size - cRF size'); ylabel('n channels'); med_diff = nanmedian(ppp-ccc); iqr_diff = [med_diff-(iqr(ppp-ccc)/2) med_diff+(iqr(ppp-ccc)/2)]; m_diff = nanmean(ppp-ccc); sem_diff = nanstd(ppp-ccc)./sqrt(sum(~isnan(ppp.*ccc))); [p,h,stats] = signrank(ppp-ccc,0); fprintf(['Difference pRF-cRF: Wilcoxon z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); plot([med_diff med_diff],[0 100],'r'); plot([m_diff m_diff],[0 100],'r--'); plot([0 0],[0 100],'k'); legend({'histogram', ['median ' num2str(med_diff)], ['mean ' num2str(m_diff)]}) title(['IQR ' num2str(iqr_diff(1)) ' ' num2str(iqr_diff(2)) ... '; SEM ' num2str(sem_diff)]) subplot(1,3,3); hold on; histogram(ppp./ccc,-0.05:0.1:10.05,'FaceColor',[.5 .5 .5],'LineStyle','none') xlabel('pRF size / cRF size'); ylabel('n channels'); med_rat = nanmedian(ppp./ccc); iqr_rat = [med_rat-(iqr(ppp./ccc)/2) med_rat+(iqr(ppp./ccc)/2)]; m_rat = nanmean(ppp./ccc); sem_rat = nanstd(ppp./ccc)./sqrt(sum(~isnan(ppp.*ccc))); [p,h,stats] = signrank(ppp./ccc,1); fprintf(['Ratio pRF/cRF: Wilcoxon z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); plot([med_rat med_rat],[0 100],'r'); plot([m_rat m_rat],[0 100],'r--'); plot([1 1],[0 100],'k'); legend({'histogram', ['median ' num2str(med_rat)], ['mean ' num2str(m_rat)]}) title(['IQR ' num2str(iqr_rat(1)) ' ' num2str(iqr_rat(2)) ... '; SEM ' num2str(sem_rat)]) if SaveFigs saveas(f_cRFpRF2,fullfile(pngfld, ['MUA_cRF-pRF_' mlabel '.png'])); saveas(f_cRFpRF2,fullfile(svgfld, ['MUA_cRF-pRF_' mlabel '.svg'])); saveas(f_cRFpRF,fullfile(pngfld, ['MUA_cRF-pRF' mlabel '.png'])); saveas(f_cRFpRF,fullfile(svgfld, ['MUA_cRF-pRF' mlabel '.svg'])); end if CloseFigs; close(f_cRFpRF); close(f_cRFpRF2); end end %% Fig 8 MUA model comparison ============================================= m=unique(tMUA.Model); R2=[]; RTH = 0; for i=1:length(m) R2 = [R2 tMUA.R2(strcmp(tMUA.Model,m{i}))]; end v1=tMUA.Area(strcmp(tMUA.Model,m{1}))==1; v4=tMUA.Area(strcmp(tMUA.Model,m{1}))==4; sub = 'both'; if strcmp(sub,'M03') || strcmp(sub,'M04') v1=v1 & strcmp(tMUA.Monkey(strcmp(tMUA.Model,m{1})),sub); v4=v4 & strcmp(tMUA.Monkey(strcmp(tMUA.Model,m{1})),sub); end f6=figure; msz=15; set(f6,'Position',[100 100 1200 1200]); for row=1:4 for column=1:4 subplot(4,4,((row-1)*4)+column); hold on; plot([0 100],[0 100],'k'); XY = [R2(v1,row+1), R2(v1,column+1)]; XYs = XY(XY(:,1)>RTH & XY(:,2)>RTH,:); scatter(XYs(:,1), XYs(:,2),msz,'Marker','o',... 'MarkerEdgeColor',[.3 .3 .3],'MarkerFaceColor',[.3 .3 .3],... 'MarkerFaceAlpha',0.5); XY = [R2(v4,row+1), R2(v4,column+1)]; XYs = XY(XY(:,1)>RTH & XY(:,2)>RTH,:); scatter(XYs(:,1), XYs(:,2),msz,'Marker','o',... 'MarkerEdgeColor',[.3 .8 .3],'MarkerFaceColor',[.3 .8 .3],... 'MarkerFaceAlpha',0.5); set(gca, 'Box','off', 'xlim', [-2 100], 'ylim',[-2 100]); xlabel(m{row+1},'interpreter','none'); ylabel(m{column+1},'interpreter','none'); title('MUA'); if row==1 && column==1 legend({'','V1','V4'},'location','NorthWest'); end end end if SaveFigs saveas(f6,fullfile(pngfld, 'EPHYS_MUA_ModelComparison.png')); saveas(f6,fullfile(svgfld, 'EPHYS_MUA_ModelComparison.svg')); end if CloseFigs; close(f6); end %% Stats model comparison R2 MUA ========================================== for RTH = [0 25] % stats V1 MR2=R2(v1,2:5); sel=logical(sum(MR2>RTH,2)); [p,tbl,stats] = kruskalwallis(MR2(sel,:),{'css','dog','p-lin','u-lin'}); [c,m,h,gnames] = multcompare(stats); for i=1:size(c,1) fprintf(['V1 RTH-' num2str(RTH) ': ' gnames{c(i,1)} ' vs ' gnames{c(i,2)} ... ', p = ' num2str(c(i,6)) '\n']) end % stats V4 MR2=R2(v4,2:5); sel=logical(sum(MR2>RTH,2)); [p,tbl,stats] = kruskalwallis(MR2(sel,:),{'css','dog','p-lin','u-lin'}); [c,m,h,gnames] = multcompare(stats); for i=1:size(c,1) fprintf(['V4 RTH-' num2str(RTH) ': ' gnames{c(i,1)} ' vs ' gnames{c(i,2)} ... ', p = ' num2str(c(i,6)) '\n']) end end %% SFig 7 & SFig 8: LFP model comparison ================================== f7=figure; set(f7,'Position',[100 100 1600 1200]); msz=15; m=unique(tMUA.Model); v1=tMUA.Area(strcmp(tMUA.Model,m{1}))==1; v4=tMUA.Area(strcmp(tMUA.Model,m{1}))==4; sub = 'both'; if strcmp(sub,'M03') || strcmp(sub,'M04') v1=v1 & strcmp(tMUA.Monkey(strcmp(tMUA.Model,m{1})),sub); v4=v4 & strcmp(tMUA.Monkey(strcmp(tMUA.Model,m{1})),sub); end % scatter dots sig=unique(tLFP.SigType); lfp_order = [3 1 2 5 4]; spn=1; fbn=1; for fb=lfp_order lfpmods{fbn,1}=[]; lfpmods{fbn,2}=[]; m=unique(tLFP.Model); % reorder --------- m=m([3 4 1 2]); modname={'P-LIN','U-LIN','CSS','DoG'}; % ----------------- R2=[]; for i=1:length(m) R2 = [R2 tLFP.R2(... strcmp(tLFP.Model,m{i}) & ... strcmp(tLFP.SigType,sig{fb}))]; end for m1=1:4 lfpmods{fbn,1}=[lfpmods{fbn,1} R2(v1,m1)]; lfpmods{fbn,2}=[lfpmods{fbn,2} R2(v4,m1)]; for m2=m1+1:4 subplot(length(sig),6,spn); hold on; plot([-5 100],[-5 100],'k'); % M1R2=R2(v1,m1); % M2R2=R2(v1,m2); % M1R2=M1R2(M1R2>0 & M2R2>0); % M2R2=M2R2(M1R2>0 & M2R2>0); scatter(R2(v1,m1), R2(v1,m2),msz,'Marker','o',... 'MarkerEdgeColor',[.3 .3 .3],'MarkerEdgeAlpha',0,... 'MarkerFaceColor',[.3 .3 .3],'MarkerFaceAlpha',0.25); scatter(R2(v4,m1), R2(v4,m2),msz,'Marker','o',... 'MarkerEdgeColor',[.3 .8 .3],'MarkerEdgeAlpha',0,... 'MarkerFaceColor',[.3 .8 .3],'MarkerFaceAlpha',0.25); % M1R2=R2(v4,m1); M1R2=M1R2(M1R2>0 & M2R2>0); % M2R2=R2(v4,m1); M1R2=M1R2(M1R2>0 & M2R2>0); % scatter(R2(v4,m1), R2(v4,m2),msz,'Marker','o',... % 'MarkerEdgeColor',[.3 .3 .3],'MarkerEdgeAlpha',0,... % 'MarkerFaceColor',[.3 .3 .3],'MarkerFaceAlpha',0.25); set(gca, 'Box','off', 'xlim', [-5 100], 'ylim',[-5 100],... 'xticklabels',{},'yticklabels',{},'TickDir','out'); xlabel(modname{m1},'interpreter','none');ylabel(modname{m2},'interpreter','none') title(sig{fb}) if spn==1 %legend({'','V1','V4'},'location','SouthEast'); end spn=spn+1; end end fbn=fbn+1; end if SaveFigs saveas(f7,fullfile(pngfld, 'EPHYS_LFP_ModelComparison.png')); saveas(f7,fullfile(svgfld, 'EPHYS_LFP_ModelComparison.svg')); end if CloseFigs; close(f7); end %% stats ================================================================== FreqNames = {'Theta','Alpha','Beta','Gamma-low','Gamma-high'}; for fb=1:5 % stats V1 fprintf(['\n= V1 LFP-' FreqNames{fb} ' ================\n']) [p,tbl,stats] = kruskalwallis(lfpmods{fb,1},modname); fprintf(['H =' num2str(tbl{2,5}(1)) ', df = ' num2str(tbl{2,3}(1)) ', p = ' num2str(tbl{2,6}(1)) '\n']) [c,m,h,gnames] = multcompare(stats); for i=1:size(c,1) fprintf([gnames{c(i,1)} ' vs ' gnames{c(i,2)} ... ', p = ' num2str(c(i,6)) '\n']) end % stats V4 fprintf(['\n= V4 LFP-' FreqNames{fb} ' ================\n']) [p,tbl,stats] = kruskalwallis(lfpmods{fb,2},modname); fprintf(['H =' num2str(tbl{2,5}(1)) ', df = ' num2str(tbl{2,3}(1)) ', p = ' num2str(tbl{2,6}(1)) '\n']) [c,m,h,gnames] = multcompare(stats); for i=1:size(c,1) fprintf([gnames{c(i,1)} ' vs ' gnames{c(i,2)} ... ', p = ' num2str(c(i,6)) '\n']) end end %% SFig 9: R2 for different ephys signals ================================= r2th=0; ephys_MOD={'linear_ephys_cv1','linear_ephys_cv1_neggain',... 'css_ephys_cv1','dog_ephys_cv1'}; ephys_MMS = MMS(:,1); for m=1:length(ephys_MOD) RR=[]; fprintf(['\n============ ' ephys_MOD{m} ' ===========\n']) model=ephys_MOD{m}; s = strcmp(tMUA.Model,model); RR=[RR tMUA.R2(s)]; s = strcmp(tLFP.Model,model); sig=unique(tLFP.SigType); lfp_order = [3 1 2 5 4]; for i=lfp_order b=strcmp(tLFP.SigType,sig{i}); RR=[RR tLFP.R2(s & b)]; end LAB=['MUA';sig(lfp_order)]; f8=figure; set(f8,'Position',[100 100 1900 1350]); sgtitle(['R2 per Model: ' model],'interpreter','none'); c=0;d=0; for ref=1:6 c=c+1; for fb=1:6 d=d+1; s=(RR(:,ref)>r2th & RR(:,fb)>r2th); subplot(6,6,d); hold on; %scatter(RR(s,ref),RR(s,fb),120,[0.3 0.3 0.3],'Marker','.'); binscatter(RR(s,ref),RR(s,fb),25,... 'XLimits', [0 100],... 'YLimits', [0 100],... 'ShowEmptyBins', 'off') colorbar; set(gca,'ColorScale','log'); colormap(inferno) caxis([1 256]) plot([0 100],[0 100],'Color',[.7 .7 .7],'LineWidth',2); xlabel(LAB{ref});ylabel(LAB{fb}); %title(['R2 ' model],'interpreter','none'); set(gca,'xlim',[0 100],'ylim',[0 100]); set(gca,'TickDir','out','xtick',[0 50 100],'ytick',[0 50 100]); end end if SaveFigs saveas(f8,fullfile(pngfld, ['EPHYS_MUA_R2_' ephys_MOD{m} '.png'])); saveas(f8,fullfile(svgfld, ['EPHYS_MUA_R2_' ephys_MOD{m} '.svg'])); end if CloseFigs; close(f8); end % Stats [p,tbl,stats] = kruskalwallis(RR,LAB'); fprintf(['H =' num2str(tbl{2,5}(1)) ', df = ' num2str(tbl{2,3}(1)) ', p = ' num2str(tbl{2,6}(1)) '\n']) [c,m,h,gnames] = multcompare(stats); for i=1:size(c,1) fprintf([gnames{c(i,1)} ' (' num2str(m(c(i,1),1)) ' +/- ' num2str(m(c(i,1),2)) ') vs ' ... gnames{c(i,2)} ' (' num2str(m(c(i,2),1)) ' +/- ' num2str(m(c(i,2),2)) '), p = ' num2str(c(i,6)) '\n']) end close all end %% Fig 11A,B: location and size across ephys channels ===================== RTH=25; SNRTH = 3; % CSS ----- midx = 3; % 2 = U-LIN, 3 = CSS fprintf(['MODEL ' ephys_MOD{midx} '\n']); % 1 = MUA, 2 = ClasRF, 3 = Theta, 4 = Alpha, 5 = Beta, 6 = Gamma-low, 7 = % Gamma-high SigCompNames = {}; SigComp_nElec =[]; SigCompDist = []; SigCompDist_columns = {'mean','std','median','iqr'}; SigCompSz = []; SigCompSz_columns = {'mean_nS1','std_nS1','median_nS1','iqr_nS1',... 'mean_nS2','std_nS2','median_nS2','iqr_nS2',... 'mean_nS2/nS1','std_nS2/nS1','median_nS2/nS1','iqr_nS2/nS1'}; sub = 'both'; rown=1; for sigidx1 = 1:7 for sigidx2 = 1:7 SigCompNames{rown,1} = PRF_EST(midx,sigidx1).sig; SigCompNames{rown,2} = PRF_EST(midx,sigidx2).sig; % threshold if strcmp(sub,'M03') || strcmp(sub,'M04') elec = PRF_EST(midx,sigidx1).R2 > RTH & ... PRF_EST(midx,sigidx2).R2 > RTH & ... strcmp(PRF_EST(midx,sigidx2).M,sub); else elec = PRF_EST(midx,sigidx1).R2 > RTH & PRF_EST(midx,sigidx2).R2 > RTH; end % calculate distance DIST = sqrt((PRF_EST(midx,sigidx1).X(elec) - PRF_EST(midx,sigidx2).X(elec)).^2 + ... (PRF_EST(midx,sigidx1).Y(elec) - PRF_EST(midx,sigidx2).Y(elec)).^2); SigCompDist = [SigCompDist; nanmean(DIST) nanstd(DIST) nanmedian(DIST) iqr(DIST)]; % calculate normalized size % normalize by MUA pRF % nS1 = PRF_EST(midx,sigidx1).S(elec)./PRF_EST(midx,1).S(elec); % nS2 = PRF_EST(midx,sigidx2).S(elec)./PRF_EST(midx,1).S(elec); % normalize by cRF nS1 = PRF_EST(midx,sigidx1).S(elec)./PRF_EST(midx,2).S(elec); nS2 = PRF_EST(midx,sigidx2).S(elec)./PRF_EST(midx,2).S(elec); SigCompSz = [SigCompSz;... nanmean(nS1(~isinf(nS1))) nanstd(nS1(~isinf(nS1))) ... nanmedian(nS1(~isinf(nS1))) iqr(nS1(~isinf(nS1))) ... nanmean(nS2(~isinf(nS2))) nanstd(nS2(~isinf(nS2))) ... nanmedian(nS2(~isinf(nS2))) iqr(nS2(~isinf(nS2))) ... nanmean(nS2(~isinf(nS1) & ~isinf(nS2))./nS1(~isinf(nS1) & ~isinf(nS2))) ... nanstd(nS2(~isinf(nS1) & ~isinf(nS2))./nS1(~isinf(nS1) & ~isinf(nS2))) ... nanmedian(nS2(~isinf(nS1) & ~isinf(nS2))./nS1(~isinf(nS1) & ~isinf(nS2))) ... iqr(nS2(~isinf(nS1) & ~isinf(nS2))./nS1(~isinf(nS1) & ~isinf(nS2)))]; SigComp_nElec = [SigComp_nElec;sum(elec)]; rown = rown+1; end end [sorted_names,sidx] = sortrows(SigCompNames); sorted_n = SigComp_nElec(sidx,:); sorted_dist = SigCompDist(sidx,:); sorted_sz = SigCompSz(sidx,:); uSig = unique(SigCompNames); nSig = size(unique(SigCompNames),1); DistMat_median = zeros(nSig); DistMat_iqr = zeros(nSig); DistMat_n = zeros(nSig); SzMat_median = zeros(nSig); SzMat_iqr = zeros(nSig); for i=1:nSig ii=((i-1)*nSig)+1; DistMat_median(:,i) = sorted_dist(ii:ii+nSig-1,3); DistMat_iqr(:,i) = sorted_dist(ii:ii+nSig-1,4)./2; DistMat_n(:,i) = sorted_n(ii:ii+nSig-1); SzMat_median(:,i) = sorted_sz(ii:ii+nSig-1,11); SzMat_iqr(:,i) = sorted_sz(ii:ii+nSig-1,12)./2; end % re-order for plot order=[4 3 5 1 2 7 6]; DistMat_median = DistMat_median(order,:); DistMat_median = DistMat_median(:,order); DistMat_iqr = DistMat_iqr(order,:); DistMat_iqr = DistMat_iqr(:,order); SzMat_median = SzMat_median(order,:); SzMat_median = SzMat_median(:,order); SzMat_iqr = SzMat_iqr(order,:); SzMat_iqr = SzMat_iqr(:,order); DistMat_n = DistMat_n(order,:); DistMat_n = DistMat_n(:,order); uSig = uSig(order); % fcss=figure; set(fcss,'Position',[10 10 2200 1200],'Renderer','painters'); colormap(viridis) %colormap(brewermap([],'RdBu')); subplot(2,3,1); imagesc(DistMat_median) set(gca,'TickDir','out','xtick',1:11,'xticklabels',uSig, 'yticklabels',uSig,... 'XTickLabelRotation',45) caxis([0 1.5]);colorbar; title('pRF distance median'); subplot(2,3,2); imagesc(DistMat_iqr) set(gca,'TickDir','out','xtick',1:11,'xticklabels',uSig, 'yticklabels',uSig,... 'XTickLabelRotation',45) caxis([0 0.8]);colorbar; title('pRF distance iqr'); subplot(2,3,3); imagesc(DistMat_n) set(gca,'TickDir','out','xtick',1:11,'xticklabels',uSig, 'yticklabels',uSig,... 'XTickLabelRotation',45,'ColorScale','log'); colorbar; caxis([1 1700]) title('n'); subplot(2,3,4); imagesc(SzMat_median) set(gca,'TickDir','out','xtick',1:11,'xticklabels',uSig, 'yticklabels',uSig,... 'XTickLabelRotation',45) caxis([0 2]);colorbar; title('pRF relative sz median'); subplot(2,3,5); imagesc(SzMat_iqr) set(gca,'TickDir','out','xtick',1:11,'xticklabels',uSig, 'yticklabels',uSig,... 'XTickLabelRotation',45) caxis([0 2]);colorbar; title('pRF relative sz iqr'); subplot(2,3,6); nMask = DistMat_n >= 10; imagesc(nMask) set(gca,'TickDir','out','xtick',1:11,'xticklabels',uSig, 'yticklabels',uSig,... 'XTickLabelRotation',45); colorbar; caxis([0 1]) title('mask n > 10'); sgtitle(['MODEL ' ephys_MOD{midx}], 'interpreter','none') if SaveFigs saveas(fcss,fullfile(pngfld, 'EPHYS_xsig_dist.png')); saveas(fcss,fullfile(svgfld, 'EPHYS_xsig_dist.svg')); end if CloseFigs; close(fcss); end fcss_sz=figure;hold on; %errorbar(1:7,SzMat_median(:,2),SzMat_iqr(:,2),'ko','LineStyle','none'); % norm by MUA pRF errorbar(1:7,SzMat_median(:,1),SzMat_iqr(:,1),'ko','LineStyle','none'); % norm by MUA cRF ylabel('Normalized pRF size') set(gca,'xlim',[0.5 7.5], 'xtick',1:7, 'xticklabels',... {'cRF-MUA','pRF-MUA','Theta','Alpha', 'Beta', 'Gam-low', 'Gam-high'},... 'XTickLabelRotation',45); if SaveFigs saveas(fcss_sz,fullfile(pngfld, 'EPHYS_xsig_sz.png')); saveas(fcss_sz,fullfile(svgfld, 'EPHYS_xsig_sz.svg')); end if CloseFigs; close(fcss_sz); end %% Fig 10: Split Alpha and Beta LFP by positive negative gain ============= R2th = 20; % minimum R2 R2enh = 5; % R2 improvement fb = {'Alpha','Beta'}; %% ALPHA Fig 10, SFig 11 ================================================== fidx = 1; DoG = tLFP(... strcmp(tLFP.Model,'dog_ephys_cv1') & strcmp(tLFP.SigType,fb{fidx}),:); lin_n = tLFP(... strcmp(tLFP.Model,'linear_ephys_cv1_neggain') & strcmp(tLFP.SigType,fb{fidx}),:); lin = tLFP(... strcmp(tLFP.Model,'linear_ephys_cv1') & strcmp(tLFP.SigType,fb{fidx}),:); css = tLFP(... strcmp(tLFP.Model,'css_ephys_cv1') & strcmp(tLFP.SigType,fb{fidx}),:); lin_nGAM1 = tLFP(... strcmp(tLFP.Model,'linear_ephys_cv1_neggain') & strcmp(tLFP.SigType,'lGamma'),:); lin_nGAM2 = tLFP(... strcmp(tLFP.Model,'linear_ephys_cv1_neggain') & strcmp(tLFP.SigType,'hGamma'),:); % U-LIN ---- f_neg2 = figure; roi=1; % only do this for V1 channels set(f_neg2,'Position',[10 10 900 800]); chan_sel = lin_n.Area==roi & lin_n.R2>R2th & lin_n.R2>lin.R2+R2enh; chan_sel2 = lin_n.Area==roi & lin_n.R2>R2th; chan_pgain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain>0; chan_ngain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain<0; chan_gam1 = lin_n.Area==roi & lin_nGAM1.R2>R2th; chan_gam2 = lin_n.Area==roi & lin_nGAM2.R2>R2th; subplot(2,2,1); hold on; % gain alpha U-LIN histogram(lin_n.gain(chan_sel2),-2000:50:2000,'FaceColor','k','FaceAlpha',0.5); histogram(lin_n.gain(chan_ngain),-2000:50:2000,'FaceColor','r','FaceAlpha',0.5); histogram(lin_n.gain(chan_pgain),-2000:50:2000,'FaceColor','b','FaceAlpha',0.5); xlabel('gain U-LIN - ALL ELEC');ylabel('nChannels'); set(gca,'xlim',[-800 1800],'TickDir','out'); MM=median(lin_n.gain(chan_sel2)); yy=get(gca,'ylim'); plot([MM MM], [0 yy(2)+40],'k','Linewidth',5) set(gca,'ylim',[0 yy(2)+10]); title('Gain') fprintf(['UNSELECTED - ALPHA MEDIAN GAIN: ' num2str(MM) ', IQR ' num2str(iqr(lin_n.gain(chan_sel))) '\n']) subplot(2,2,2); hold on; % gain alpha U-LIN histogram(lin_n.gain(chan_sel),-1000:50:1000,'FaceColor','k','FaceAlpha',0.5); xlabel('gain LIN-POSNEG');ylabel('nChannels'); set(gca,'xlim',[-800 1700],'TickDir','out'); MM=median(lin_n.gain(chan_sel)); yy=get(gca,'ylim'); plot([MM MM], [0 yy(2)+40],'k','Linewidth',5) set(gca,'ylim',[0 yy(2)+10]); title('Gain selected channels (based on R2 ULIN>PLIN') fprintf(['ALPHA MEDIAN GAIN: ' num2str(MM) ', IQR ' num2str(iqr(lin_n.gain(chan_sel))) '\n']) % Wilcoxon 1-tailed < 1 [p,h,stats] = signrank(lin_n.gain(chan_sel),0,'tail','left'); fprintf(['Gain < 0: Wilcoxon z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); subplot(2,2,3); hold on; bb = [lin_n.ecc(chan_sel) lin.ecc(chan_sel)]; mEcc = [mean(lin_n.ecc(chan_ngain)) mean(lin_n.ecc(chan_pgain))]; sdEcc = [std(lin_n.ecc(chan_ngain)) std(lin_n.ecc(chan_pgain))]; mdEcc = [median(lin_n.ecc(chan_ngain)) median(lin_n.ecc(chan_pgain)) ]; iqrEcc = [iqr(lin_n.ecc(chan_ngain)) iqr(lin_n.ecc(chan_pgain)) ]./2; errorbar(1:2,mdEcc,iqrEcc,... 'ko','MarkerSize',4,'MarkerFaceColor','k','Linewidth',2) set(gca,'xtick',1:2,'xticklabels',{'ULIN-N','ULIN-P'},... 'xlim',[0.8 2.2],'TickDir','out') ylabel('Eccentricity'); title('Ecc') % Ecc difference fprintf('-- STATS Ecc --\n'); fprintf(['mEcc nGain U-LIN: ' num2str(mEcc(1)) ', STD ' num2str(sdEcc(1)) '\n']) fprintf(['mEcc pGain U-LIN: ' num2str(mEcc(2)) ', STD ' num2str(sdEcc(2)) '\n']) fprintf(['mdEcc nGain U-LIN: ' num2str(mdEcc(1)) ', IQR ' num2str(iqrEcc(1)) '\n']) fprintf(['mdEcc pGain U-LIN: ' num2str(mdEcc(2)) ', IQR ' num2str(iqrEcc(2)) '\n']) [p,h,stats] = ranksum(lin_n.ecc(chan_ngain),lin_n.ecc(chan_pgain)); fprintf(['Ecc difference pos gain U-LIN vs neg gain U-LIN: Wilcoxon z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); subplot(2,2,4); hold on; bb = [lin_n.rfs(chan_sel) lin.rfs(chan_sel)]; mSz = [mean(lin_n.rfs(chan_ngain)) mean(lin_n.rfs(chan_pgain)) ]; sdSz = [std(lin_n.rfs(chan_ngain)) std(lin_n.rfs(chan_pgain)) ]; mdSz = [median(lin_n.rfs(chan_ngain)) median(lin_n.rfs(chan_pgain))]; iqrSz = [iqr(lin_n.rfs(chan_ngain)) iqr(lin_n.rfs(chan_pgain)) ]./2; errorbar(1:2,mdSz,iqrSz,'ko','MarkerSize',4,'MarkerFaceColor','k','Linewidth',2) set(gca,'xtick',1:2,'xticklabels',{'ULIN-N','ULIN-P'},... 'xlim',[0.8 2.2],'TickDir','out') ylabel('Size'); title('Size') % Sz difference fprintf('-- STATS Sz --\n'); fprintf(['mSz nGain U-LIN: ' num2str(mSz(1)) ', STD ' num2str(sdSz(1)) '\n']) fprintf(['mSz pGain U-LIN: ' num2str(mSz(2)) ', STD ' num2str(sdSz(2)) '\n']) fprintf(['mdSz nGain U-LIN: ' num2str(mdSz(1)) ', IQR ' num2str(iqrSz(1)) '\n']) fprintf(['mdSz pGain U-LIN: ' num2str(mdSz(2)) ', IQR ' num2str(iqrSz(2)) '\n']) [p,h,stats] = ranksum(lin_n.rfs(chan_ngain),lin_n.rfs(chan_pgain)); fprintf(['Sz difference pos gain U-LIN vs neg gain U-LIN: Wilcoxon z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); if SaveFigs saveas(f_neg2,fullfile(pngfld, 'ALPHA_posneg_ecc_sz.png')); saveas(f_neg2,fullfile(svgfld, 'ALPHA_posneg_ecc_sz.svg')); end if CloseFigs; close(f_neg2); end %% Compare eccentricities between alpha and gamma prfs ==================== fextra=figure; set(fextra,'Position',[100 100 600 1500]); dE=[]; subplot(4,2,1);hold on; plot([0 15],[0 15]); scatter(lin_n.ecc(chan_ngain),lin_nGAM1.ecc(chan_ngain)); set(gca,'xlim',[0 15],'ylim',[0 15]); xlabel('A- CHAN: Ecc LOW GAM');ylabel('A- CHAN: Ecc ALPHA'); fprintf('Are Negative ALPHA pRF shifted relative to fixation from Low GAMMA?\n'); [p,h,stats] = signrank(lin_n.ecc(chan_ngain),lin_nGAM1.ecc(chan_ngain)); fprintf(['Median dEcc is ' ... num2str(median(lin_nGAM1.ecc(chan_ngain)-lin_n.ecc(chan_ngain))) ... ' IQR ' ... num2str(iqr(lin_nGAM1.ecc(chan_ngain)-lin_n.ecc(chan_ngain))./2) ... ' (POS --> towards fixation) \n']); fprintf(['p = ' num2str(p) '\n']); dE=[dE; ... median(lin_nGAM1.ecc(chan_ngain)-lin_n.ecc(chan_ngain)) ... iqr(lin_nGAM1.ecc(chan_ngain)-lin_n.ecc(chan_ngain))./2]; subplot(4,2,2);hold on; plot([0 15],[0 15]); scatter(lin_n.ecc(chan_pgain),lin_nGAM1.ecc(chan_pgain)); set(gca,'xlim',[0 15],'ylim',[0 15]); xlabel('A+ CHAN: Ecc LOW GAM');ylabel('A+ CHAN: Ecc ALPHA'); fprintf('Are Positive ALPHA pRF shifted relative to fixation from Low GAMMA?\n'); [p,h,stats] = signrank(lin_n.ecc(chan_pgain),lin_nGAM1.ecc(chan_pgain)); fprintf(['Median dEcc is ' ... num2str(median(lin_nGAM1.ecc(chan_pgain)-lin_n.ecc(chan_pgain))) ... ' IQR ' ... num2str(iqr(lin_nGAM1.ecc(chan_pgain)-lin_n.ecc(chan_pgain))./2) ... ' (POS --> towards fixation) \n']); fprintf(['p = ' num2str(p) '\n']); dE=[dE; ... median(lin_nGAM1.ecc(chan_pgain)-lin_n.ecc(chan_pgain)) ... iqr(lin_nGAM1.ecc(chan_pgain)-lin_n.ecc(chan_pgain))./2]; subplot(4,2,3);hold on; plot([0 15],[0 15]); scatter(lin_n.ecc(chan_ngain),lin_nGAM2.ecc(chan_ngain)); set(gca,'xlim',[0 15],'ylim',[0 15]); xlabel('A- CHAN: Ecc HIGH GAM');ylabel('A- CHAN: Ecc ALPHA'); fprintf('Are Negative ALPHA pRF shifted relative to fixation from High GAMMA?\n'); [p,h,stats] = signrank(lin_n.ecc(chan_ngain),lin_nGAM2.ecc(chan_ngain)); fprintf(['Median dEcc is ' ... num2str(median(lin_nGAM2.ecc(chan_ngain)-lin_n.ecc(chan_ngain))) ... ' IQR ' ... num2str(iqr(lin_nGAM2.ecc(chan_ngain)-lin_n.ecc(chan_ngain))./2) ... ' (POS --> towards fixation) \n']); fprintf(['p = ' num2str(p) '\n']); dE=[dE; ... median(lin_nGAM2.ecc(chan_ngain)-lin_n.ecc(chan_ngain)) ... iqr(lin_nGAM2.ecc(chan_ngain)-lin_n.ecc(chan_ngain))./2]; subplot(4,2,4);hold on; plot([0 15],[0 15]); scatter(lin_n.ecc(chan_pgain),lin_nGAM2.ecc(chan_pgain)); set(gca,'xlim',[0 15],'ylim',[0 15]); xlabel('A+ CHAN: Ecc HIGH GAM');ylabel('A+ CHAN: Ecc ALPHA'); fprintf('Are Positive ALPHA pRF shifted relative to fixation from High GAMMA?\n'); [p,h,stats] = signrank(lin_n.ecc(chan_pgain),lin_nGAM2.ecc(chan_pgain)); fprintf(['Median dEcc is ' ... num2str(median(lin_nGAM2.ecc(chan_pgain)-lin_n.ecc(chan_pgain))) ... ' IQR ' ... num2str(iqr(lin_nGAM2.ecc(chan_pgain)-lin_n.ecc(chan_pgain))./2) ... ' (POS --> towards fixation) \n']); fprintf(['p = ' num2str(p) '\n']); dE=[dE; ... median(lin_nGAM2.ecc(chan_pgain)-lin_n.ecc(chan_pgain)) ... iqr(lin_nGAM2.ecc(chan_pgain)-lin_n.ecc(chan_pgain))./2]; subplot(4,2,5:8);hold on; hold on; bar(1:4, dE(:,1)); errorbar(1:4, dE(:,1), dE(:,2), 'k','Linestyle', 'none'); ylabel('dEcc Alpha vs Gam') set(gca,'xtick',1:4,'xticklabels',{'A- v lG','A+ v lG','A- v hG','A+ v hG'}) %% plot the diff in within channel pRF location for low/high LFP freq ===== fshiftprf = figure; set(fshiftprf,'Position',[100 100 1600 1600]); Xa1 = lin_n.X(chan_ngain); Ya1 = lin_n.Y(chan_ngain); Xlg1 = lin_nGAM1.X(chan_ngain); Ylg1 = lin_nGAM1.Y(chan_ngain); Xa2 = lin_n.X(chan_pgain); Ya2 = lin_n.Y(chan_pgain); Xlg2 = lin_nGAM1.X(chan_pgain); Ylg2 = lin_nGAM1.Y(chan_pgain); % only prfs within 8 deg for display csel1 = Xa1<8 & Ya1<2 & Xlg1<8 & Ylg1<2 & ... Xa1>-2 & Ya1>-8 & Xlg1>-2 & Ylg1>-8; csel2 = Xa2<8 & Ya2<2 & Xlg2<8 & Ylg2<2 & ... Xa2>-2 & Ya2>-8 & Xlg2>-2 & Ylg2>-8; subplot(2,2,1);hold on; plot([0 0],[-10 10],'--','LineWidth',2,'Color',[0.5 0.5 0.5]); plot([-10 10],[0 0],'--','LineWidth',2,'Color',[0.5 0.5 0.5]); quiver(Xlg1(csel1),Ylg1(csel1),Xa1(csel1)-Xlg1(csel1),Ya1(csel1)-Ylg1(csel1),... 'Color',[0.8 0.8 0.8],... 'AutoScale',0,'ShowArrowHead',0); % lines from gam to alpha pRF scatter(Xa1(csel1),Ya1(csel1),20,'ko','MarkerFaceColor','r'); % alpha pRF center scatter(Xlg1(csel1),Ylg1(csel1),20,'ko','MarkerFaceColor','b'); % lgam prf center % quiver(Xlg1(csel1),Ylg1(csel1),Xa1(csel1)-Xlg1(csel1),Ya1(csel1)-Ylg1(csel1),... % 'k','LineWidth',1,'AutoScale',1,'ShowArrowHead',1); % arrows from gam to alpha pRF plot(0,0,'go','MarkerSize',8,'MarkerFaceColor','g'); % fovea set(gca,'xlim',[-1.5 8.2],'ylim',[-6.5 2.2],'FontSize',14,... 'FontName','Myriad Pro','TickDir','out'); xlabel('Horizontal dva');ylabel('Vertical dva'); shiftsz = sqrt((Xa1(csel1)-Xlg1(csel1)).^2 + (Ya1(csel1)-Ylg1(csel1)).^2); m_ss = mean(shiftsz); se_ss = std(shiftsz)./sqrt(length(shiftsz)); fprintf(['Average negALPHA shift size: ' num2str(m_ss) ' +/- ' num2str(se_ss) ' SEM\n']); subplot(2,2,2);hold on; quiver(Xlg1(csel1),Ylg1(csel1),Xa1(csel1)-Xlg1(csel1),Ya1(csel1)-Ylg1(csel1),... 'k','AutoScale',1,'ShowArrowHead',1) plot(0,0,'go','MarkerSize',8,'MarkerFaceColor','g') set(gca,'xlim',[-1.5 8.2],'ylim',[-6.5 2.2],'FontSize',14,... 'FontName','Myriad Pro','TickDir','out'); xlabel('Horizontal dva');ylabel('Vertical dva'); subplot(2,2,3);hold on; plot([0 0],[-10 10],'--','LineWidth',2,'Color',[0.5 0.5 0.5]); plot([-10 10],[0 0],'--','LineWidth',2,'Color',[0.5 0.5 0.5]); quiver(Xlg2(csel2),Ylg2(csel2),Xa2(csel2)-Xlg2(csel2),Ya2(csel2)-Ylg2(csel2),... 'Color',[0.8 0.8 0.8],... 'AutoScale',0,'ShowArrowHead',0); % lines from gam to alpha pRF scatter(Xa2(csel2),Ya2(csel2),20,'ko','MarkerFaceColor','r'); % alpha pRF center scatter(Xlg2(csel2),Ylg2(csel2),20,'ko','MarkerFaceColor','b'); % lgam prf center % quiver(Xlg2(csel2),Ylg2(csel2),Xa2(csel2)-Xlg2(csel2),Ya2(csel2)-Ylg2(csel2),... % 'k','LineWidth',1,'AutoScale',1,'ShowArrowHead',1); % arrows from gam to alpha pRF plot(0,0,'go','MarkerSize',8,'MarkerFaceColor','g'); % fovea set(gca,'xlim',[-1.5 8.2],'ylim',[-6.5 2.2],'FontSize',14,... 'FontName','Myriad Pro','TickDir','out'); xlabel('Horizontal dva');ylabel('Vertical dva'); shiftsz = sqrt((Xa2(csel2)-Xlg2(csel2)).^2 + (Ya2(csel2)-Ylg2(csel2)).^2); m_ss = mean(shiftsz); se_ss = std(shiftsz)./sqrt(length(shiftsz)); fprintf(['Average posALPHA shift size: ' num2str(m_ss) ' +/- ' num2str(se_ss) ' SEM\n']); subplot(2,2,4);hold on; quiver(Xlg2(csel2),Ylg2(csel2),Xa2(csel2)-Xlg2(csel2),Ya2(csel2)-Ylg2(csel2),... 'k','AutoScale',1,'ShowArrowHead',1) plot(0,0,'go','MarkerSize',8,'MarkerFaceColor','g') set(gca,'xlim',[-1.5 8.2],'ylim',[-6.5 2.2],'FontSize',14,... 'FontName','Myriad Pro','TickDir','out'); xlabel('Horizontal dva');ylabel('Vertical dva'); if SaveFigs saveas(fextra,fullfile(pngfld, 'ALPHAvGAM_posneg_ecc_sz.png')); saveas(fextra,fullfile(svgfld, 'ALPHAvGAM_posneg_ecc_sz.svg')); end if CloseFigs; close(fextra); end %% ALPHA vs GAMMA pRF separation index ==================================== lin_n = tLFP(... strcmp(tLFP.Model,'linear_ephys_cv1_neggain') & strcmp(tLFP.SigType,'Alpha'),:); DoG = tMUA(... strcmp(tMUA.Model,'dog_ephys_cv1'),:); roi=1; % only do this for V1 channels chan_sel = lin_n.Area==roi & lin_n.R2>R2th & lin_n.R2>lin.R2+R2enh; chan_sel2 = lin_n.Area==roi & lin_n.R2>R2th; chan_pgain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain>0; chan_ngain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain<0; chan_gam1 = lin_n.Area==roi & lin_nGAM1.R2>R2th; chan_gam2 = lin_n.Area==roi & lin_nGAM2.R2>R2th; [mean(sqrt((lin_n.X(chan_ngain)-lin_nGAM1.X(chan_ngain)).^2 + ... (lin_n.Y(chan_ngain)-lin_nGAM1.Y(chan_ngain)).^2)) ... std(sqrt((lin_n.X(chan_ngain)-lin_nGAM1.X(chan_ngain)).^2 + ... (lin_n.Y(chan_ngain)-lin_nGAM1.Y(chan_ngain)).^2))./sqrt(sum(chan_ngain))] [mean(sqrt((lin_n.X(chan_pgain)-lin_nGAM1.X(chan_pgain)).^2 + ... (lin_n.Y(chan_pgain)-lin_nGAM1.Y(chan_pgain)).^2)) ... std(sqrt((lin_n.X(chan_pgain)-lin_nGAM1.X(chan_pgain)).^2 + ... (lin_n.Y(chan_pgain)-lin_nGAM1.Y(chan_pgain)).^2))./sqrt(sum(chan_pgain))] [mean(sqrt((lin_n.X(chan_ngain)-lin_nGAM2.X(chan_ngain)).^2 + ... (lin_n.Y(chan_ngain)-lin_nGAM1.Y(chan_ngain)).^2)) ... std(sqrt((lin_n.X(chan_ngain)-lin_nGAM2.X(chan_ngain)).^2 + ... (lin_n.Y(chan_ngain)-lin_nGAM2.Y(chan_ngain)).^2))./sqrt(sum(chan_ngain))] [mean(sqrt((lin_n.X(chan_pgain)-lin_nGAM2.X(chan_pgain)).^2 + ... (lin_n.Y(chan_pgain)-lin_nGAM2.Y(chan_pgain)).^2)) ... std(sqrt((lin_n.X(chan_pgain)-lin_nGAM2.X(chan_pgain)).^2 + ... (lin_n.Y(chan_pgain)-lin_nGAM2.Y(chan_pgain)).^2))./sqrt(sum(chan_pgain))] fprintf('=== ALPHA ===\n') % size fprintf('Size diff (a-lG)\n') fprintf('NGAIN\n') [mean(lin_n.rfs(chan_ngain)-lin_nGAM1.rfs(chan_ngain)) ... std(lin_n.rfs(chan_ngain)-lin_nGAM1.rfs(chan_ngain))./sqrt(sum(chan_ngain))] fprintf('PGAIN\n') [mean(lin_n.rfs(chan_pgain)-lin_nGAM1.rfs(chan_pgain)) ... std(lin_n.rfs(chan_pgain)-lin_nGAM1.rfs(chan_pgain))./sqrt(sum(chan_pgain))] % normalized size fprintf('Norm Size diff (a/lG)\n') fprintf('NGAIN\n') [mean(lin_n.rfs(chan_ngain)./lin_nGAM1.rfs(chan_ngain)) ... std(lin_n.rfs(chan_ngain)./lin_nGAM1.rfs(chan_ngain))./sqrt(sum(chan_ngain))] fprintf('PGAIN\n') [mean(lin_n.rfs(chan_pgain)./lin_nGAM1.rfs(chan_pgain)) ... std(lin_n.rfs(chan_pgain)./lin_nGAM1.rfs(chan_pgain))./sqrt(sum(chan_pgain))] fprintf('Norm Size diff (a/hG)\n') fprintf('NGAIN\n') [mean(lin_n.rfs(chan_ngain)./lin_nGAM2.rfs(chan_ngain)) ... std(lin_n.rfs(chan_ngain)./lin_nGAM2.rfs(chan_ngain))./sqrt(sum(chan_ngain))] fprintf('PGAIN\n') [mean(lin_n.rfs(chan_pgain)./lin_nGAM2.rfs(chan_pgain)) ... std(lin_n.rfs(chan_pgain)./lin_nGAM2.rfs(chan_pgain))./sqrt(sum(chan_pgain))] % distance in relation to size distn = sqrt((lin_n.X(chan_ngain)-lin_nGAM1.X(chan_ngain)).^2 + ... (lin_n.Y(chan_ngain)-lin_nGAM1.Y(chan_ngain)).^2); distp = sqrt((lin_n.X(chan_pgain)-lin_nGAM1.X(chan_pgain)).^2 + ... (lin_n.Y(chan_pgain)-lin_nGAM1.Y(chan_pgain)).^2); sz2n=lin_n.rfs(chan_ngain)+lin_nGAM1.rfs(chan_ngain); sz2p=lin_n.rfs(chan_pgain)+lin_nGAM1.rfs(chan_pgain); ff=figure; set(ff,'Position',[10 10 600 1500]); fprintf('DIST REL/Sz (dist/(s1+s2))\n') subplot(4,2,1);hold on; scatter(distn,sz2n); scatter(distp,sz2p); plot([0 20],[0 20]); xlabel('distance A-G pRF');ylabel('sumsize A-G pRF') set(gca,'xlim',[0 20],'ylim',[0 20]); legend({'negative', 'positive'}); subplot(4,2,2);hold on; fprintf('lG NEG\n') fprintf('mean std median iqr\n') [mean(distn./sz2n) std(distn./sz2n)./sqrt(sum(chan_ngain))... median(distn./sz2n) iqr(distn./sz2n)./2] histogram(distn./sz2n,100,'Normalization','probability') fprintf('lG POS\n') [mean(distp./sz2p) std(distp./sz2p)./sqrt(sum(chan_pgain))... median(distp./sz2p) iqr(distp./sz2p)./2] histogram(distp./sz2p,100,'Normalization','probability') title('Low Gamma'); xlabel('Separation Index'); BC=[median(distn./sz2n) median(distp./sz2p)]; BCE=[iqr(distn./sz2n)/2 iqr(distp./sz2p)/2]; [p,h,stats] = ranksum(distn./sz2n, distp./sz2p); fprintf('SEPARATION INDEX [Alpha-lGam] - Neg vs Pos channels\n'); fprintf(['p = ' num2str(p) '\n']); distn = sqrt((lin_n.X(chan_ngain)-lin_nGAM2.X(chan_ngain)).^2 + ... (lin_n.Y(chan_ngain)-lin_nGAM2.Y(chan_ngain)).^2); distp = sqrt((lin_n.X(chan_pgain)-lin_nGAM2.X(chan_pgain)).^2 + ... (lin_n.Y(chan_pgain)-lin_nGAM2.Y(chan_pgain)).^2); sz2n=lin_n.rfs(chan_ngain)+lin_nGAM2.rfs(chan_ngain); sz2p=lin_n.rfs(chan_pgain)+lin_nGAM2.rfs(chan_pgain); subplot(4,2,3);hold on; scatter(distn,sz2n); scatter(distp,sz2p); plot([0 20],[0 20]); xlabel('distance A-G pRF');ylabel('sumsize A-G pRF') set(gca,'xlim',[0 20],'ylim',[0 20]); legend({'negative', 'positive'}); subplot(4,2,4);hold on; fprintf('hG NEG\n') fprintf('mean std median iqr\n') [mean(distn./sz2n) std(distn./sz2n)./sqrt(sum(chan_ngain))... median(distn./sz2n) iqr(distn./sz2n)./2] histogram(distn./sz2n,100,'Normalization','probability') fprintf('hG POS\n') [mean(distp./sz2p) std(distp./sz2p)./sqrt(sum(chan_pgain))... median(distp./sz2p) iqr(distp./sz2p)./2] histogram(distp./sz2p,100,'Normalization','probability') title('High Gamma'); xlabel('Separation Index'); BC=[BC median(distn./sz2n) median(distp./sz2p)]; BCE=[BCE iqr(distn./sz2n)/2 iqr(distp./sz2p)/2]; [p,h,stats] = ranksum(distn./sz2n, distp./sz2p); fprintf('SEPARATION INDEX [Alpha-hGam] - Neg vs Pos channels\n'); fprintf(['p = ' num2str(p) '\n']); subplot(4,2,5:8); hold on; bar(1:4, BC); errorbar(1:4, BC, BCE, 'k', 'Linestyle','none') ylabel('Separation Index') set(gca, 'xtick',1:4,'xticklabels',{'lG-','lG+','hG-','hG+'}) if SaveFigs saveas(ff,fullfile(pngfld, 'ALPHAvGAM_sepidx.png')); saveas(ff,fullfile(svgfld, 'ALPHAvGAM_sepidx.svg')); end if CloseFigs; close(ff); end %% ALPHA DoG ============================================================== chan_sel = DoG.R2>R2th & DoG.R2>lin.R2+R2enh & ... DoG.normamp~=0 & DoG.ecc<16; MM=median(DoG.normamp(chan_sel)); fprintf(['ALPHA MEDIAN NAMP: ' num2str(MM) ', IQR ' num2str(iqr(DoG.normamp(chan_sel))) '\n']) % Wilcoxon 1-tailed < 1 [p,h,stats] = signrank(lin_n.gain(chan_sel),0,'tail','left'); fprintf(['Gain < 0: Wilcoxon z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); bb = [DoG.ecc(chan_sel) lin.ecc(chan_sel)]; % Wilcoxon [p,h,stats] = signrank(bb(:,1),bb(:,2)); fprintf(['ECC diff Wilcoxon z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); MM=median(bb(:,2)-bb(:,1)); fprintf(['ALPHA MEDIAN ECC DIFF: ' num2str(median(lin.ecc(chan_sel)-DoG.ecc(chan_sel))) ... ', IQR ' num2str(iqr(lin.ecc(chan_sel)-DoG.ecc(chan_sel))) '\n']) %% BETA SFig 10, 11 ======================================================= fidx = 2; DoG = tLFP(... strcmp(tLFP.Model,'dog_ephys_cv1') & strcmp(tLFP.SigType,fb{fidx}),:); lin_n = tLFP(... strcmp(tLFP.Model,'linear_ephys_cv1_neggain') & strcmp(tLFP.SigType,fb{fidx}),:); lin = tLFP(... strcmp(tLFP.Model,'linear_ephys_cv1') & strcmp(tLFP.SigType,fb{fidx}),:); css = tLFP(... strcmp(tLFP.Model,'css_ephys_cv1') & strcmp(tLFP.SigType,fb{fidx}),:); % U-LIN ---- f_neg2 = figure; roi=1; % only do this for V1 channels set(f_neg2,'Position',[10 10 900 1200]); chan_sel = lin_n.Area==roi & lin_n.R2>R2th & lin_n.R2>lin.R2+R2enh; chan_sel2 = lin_n.Area==roi & lin_n.R2>R2th; chan_pgain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain>0; chan_ngain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain<0; subplot(2,2,1); hold on; % gain alpha U-LIN histogram(lin_n.gain(chan_sel2),-2000:50:2000,'FaceColor','k','FaceAlpha',0.5); histogram(lin_n.gain(chan_ngain),-2000:50:2000,'FaceColor','r','FaceAlpha',0.5); histogram(lin_n.gain(chan_pgain),-2000:50:2000,'FaceColor','b','FaceAlpha',0.5); xlabel('gain U-LIN - ALL ELEC');ylabel('nChannels'); set(gca,'xlim',[-800 1800],'TickDir','out'); MM=median(lin_n.gain(chan_sel2)); yy=get(gca,'ylim'); plot([MM MM], [0 yy(2)+40],'k','Linewidth',5) set(gca,'ylim',[0 yy(2)+10]); title('Gain') fprintf(['UNSELECTED - BETA MEDIAN GAIN: ' num2str(MM) ', IQR ' num2str(iqr(lin_n.gain(chan_sel))) '\n']) subplot(2,2,2); hold on; % gain alpha U-LIN histogram(lin_n.gain(chan_sel),-1000:50:1000,'FaceColor','k','FaceAlpha',0.5); xlabel('gain LIN-POSNEG');ylabel('nChannels'); set(gca,'xlim',[-800 1700],'TickDir','out'); MM=median(lin_n.gain(chan_sel)); yy=get(gca,'ylim'); plot([MM MM], [0 yy(2)+40],'k','Linewidth',5) set(gca,'ylim',[0 yy(2)+10]); title('Gain selected channels (based on R2 ULIN>PLIN') fprintf(['BETA MEDIAN GAIN: ' num2str(MM) ', IQR ' num2str(iqr(lin_n.gain(chan_sel))) '\n']) % Wilcoxon 1-tailed < 1 [p,h,stats] = signrank(lin_n.gain(chan_sel),0,'tail','left'); fprintf(['Gain < 0: Wilcoxon z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); subplot(2,2,3); hold on; bb = [lin_n.ecc(chan_sel) lin.ecc(chan_sel)]; mEcc = [mean(lin_n.ecc(chan_ngain)) mean(lin_n.ecc(chan_pgain))]; sdEcc = [std(lin_n.ecc(chan_ngain)) std(lin_n.ecc(chan_pgain))]; mdEcc = [median(lin_n.ecc(chan_ngain)) median(lin_n.ecc(chan_pgain))]; iqrEcc = [iqr(lin_n.ecc(chan_ngain)) iqr(lin_n.ecc(chan_pgain))]./2; errorbar(1:2,mdEcc,iqrEcc,... 'ko','MarkerSize',4,'MarkerFaceColor','k','Linewidth',2) set(gca,'xtick',1:2,'xticklabels',{'ULIN-N','ULIN-P'},... 'xlim',[0.8 2.2],'TickDir','out') ylabel('Eccentricity'); title('Ecc') % Ecc difference fprintf('-- STATS Ecc --\n'); fprintf(['mEcc nGain U-LIN: ' num2str(mEcc(1)) ', STD ' num2str(sdEcc(1)) '\n']) fprintf(['mEcc pGain U-LIN: ' num2str(mEcc(2)) ', STD ' num2str(sdEcc(2)) '\n']) fprintf(['mdEcc nGain U-LIN: ' num2str(mdEcc(1)) ', IQR ' num2str(iqrEcc(1)) '\n']) fprintf(['mdEcc pGain U-LIN: ' num2str(mdEcc(2)) ', IQR ' num2str(iqrEcc(2)) '\n']) [p,h,stats] = ranksum(lin_n.ecc(chan_ngain),lin_n.ecc(chan_pgain)); fprintf(['Ecc difference pos gain U-LIN vs neg gain U-LIN: Wilcoxon z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); subplot(2,2,4); hold on; bb = [lin_n.rfs(chan_sel) lin.rfs(chan_sel)]; mSz = [mean(lin_n.rfs(chan_ngain)) mean(lin_n.rfs(chan_pgain))]; sdSz = [std(lin_n.rfs(chan_ngain)) std(lin_n.rfs(chan_pgain))]; mdSz = [median(lin_n.rfs(chan_ngain)) median(lin_n.rfs(chan_pgain))]; iqrSz = [iqr(lin_n.rfs(chan_ngain)) iqr(lin_n.rfs(chan_pgain))]./2; errorbar(1:2,mdSz,iqrSz,'ko','MarkerSize',4,'MarkerFaceColor','k','Linewidth',2) set(gca,'xtick',1:4,'xticklabels',{'ULIN-N','ULIN-P'},... 'xlim',[0.8 2.2],'TickDir','out') ylabel('Size'); title('Size') % Sz difference fprintf('-- STATS Sz --\n'); fprintf(['mSz nGain U-LIN: ' num2str(mSz(1)) ', STD ' num2str(sdSz(1)) '\n']) fprintf(['mSz pGain U-LIN: ' num2str(mSz(2)) ', STD ' num2str(sdSz(2)) '\n']) fprintf(['mdSz nGain U-LIN: ' num2str(mdSz(1)) ', IQR ' num2str(iqrSz(1)) '\n']) fprintf(['mdSz pGain U-LIN: ' num2str(mdSz(2)) ', IQR ' num2str(iqrSz(2)) '\n']) [p,h,stats] = ranksum(lin_n.rfs(chan_ngain),lin_n.rfs(chan_pgain)); fprintf(['Sz difference pos gain U-LIN vs neg gain U-LIN: Wilcoxon z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); if SaveFigs saveas(f_neg2,fullfile(pngfld, 'BETA_posneg_ecc_sz.png')); saveas(f_neg2,fullfile(svgfld, 'BETA_posneg_ecc_sz.svg')); end if CloseFigs; close(f_neg2); end %% Compare eccentricities between beta and gamma prfs ===================== lin_n = tLFP(... strcmp(tLFP.Model,'linear_ephys_cv1_neggain') & strcmp(tLFP.SigType,'Beta'),:); roi=1; % only do this for V1 channels chan_sel = lin_n.Area==roi & lin_n.R2>R2th & lin_n.R2>lin.R2+R2enh; chan_sel2 = lin_n.Area==roi & lin_n.R2>R2th; chan_pgain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain>0; chan_ngain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain<0; chan_gam1 = lin_n.Area==roi & lin_nGAM1.R2>R2th; chan_gam2 = lin_n.Area==roi & lin_nGAM2.R2>R2th; chan_gam1 = lin_n.Area==roi & lin_nGAM1.R2>R2th; chan_gam2 = lin_n.Area==roi & lin_nGAM2.R2>R2th; fextra=figure; set(fextra,'Position',[10 10 600 1500]); dE=[]; subplot(4,2,1);hold on; plot([0 15],[0 15]); scatter(lin_n.ecc(chan_ngain),lin_nGAM1.ecc(chan_ngain)); set(gca,'xlim',[0 15],'ylim',[0 15]); xlabel('B- CHAN: Ecc LOW GAM');ylabel('B- CHAN: Ecc BETA'); fprintf('Are Negative BETA pRF shifted relative to fixation from Low GAMMA?\n'); [p,h,stats] = signrank(lin_n.ecc(chan_ngain),lin_nGAM1.ecc(chan_ngain)); fprintf(['Median dEcc is ' ... num2str(median(lin_nGAM1.ecc(chan_ngain)-lin_n.ecc(chan_ngain))) ... ' IQR ' ... num2str(iqr(lin_nGAM1.ecc(chan_ngain)-lin_n.ecc(chan_ngain))./2) ... ' (POS --> towards fixation) \n']); fprintf(['p = ' num2str(p) '\n']); dE=[dE; ... median(lin_nGAM1.ecc(chan_ngain)-lin_n.ecc(chan_ngain)) ... iqr(lin_nGAM1.ecc(chan_ngain)-lin_n.ecc(chan_ngain))./2]; subplot(4,2,2);hold on; plot([0 15],[0 15]); scatter(lin_n.ecc(chan_pgain),lin_nGAM1.ecc(chan_pgain)); set(gca,'xlim',[0 15],'ylim',[0 15]); xlabel('B+ CHAN: Ecc LOW GAM');ylabel('B+ CHAN: Ecc BETA'); fprintf('Are Positive BETA pRF shifted relative to fixation from Low GAMMA?\n'); [p,h,stats] = signrank(lin_n.ecc(chan_pgain),lin_nGAM1.ecc(chan_pgain)); fprintf(['Median dEcc is ' ... num2str(median(lin_nGAM1.ecc(chan_pgain)-lin_n.ecc(chan_pgain))) ... ' IQR ' ... num2str(iqr(lin_nGAM1.ecc(chan_pgain)-lin_n.ecc(chan_pgain))./2) ... ' (POS --> towards fixation) \n']); fprintf(['p = ' num2str(p) '\n']); dE=[dE; ... median(lin_nGAM1.ecc(chan_pgain)-lin_n.ecc(chan_pgain)) ... iqr(lin_nGAM1.ecc(chan_pgain)-lin_n.ecc(chan_pgain))./2]; subplot(4,2,3);hold on; plot([0 15],[0 15]); scatter(lin_n.ecc(chan_ngain),lin_nGAM2.ecc(chan_ngain)); set(gca,'xlim',[0 15],'ylim',[0 15]); xlabel('B- CHAN: Ecc HIGH GAM');ylabel('B- CHAN: Ecc BETA'); fprintf('Are Negative BETA pRF shifted relative to fixation from High GAMMA?\n'); [p,h,stats] = signrank(lin_n.ecc(chan_ngain),lin_nGAM2.ecc(chan_ngain)); fprintf(['Median dEcc is ' ... num2str(median(lin_nGAM2.ecc(chan_ngain)-lin_n.ecc(chan_ngain))) ... ' IQR ' ... num2str(iqr(lin_nGAM2.ecc(chan_ngain)-lin_n.ecc(chan_ngain))./2) ... ' (POS --> towards fixation) \n']); fprintf(['p = ' num2str(p) '\n']); dE=[dE; ... median(lin_nGAM2.ecc(chan_ngain)-lin_n.ecc(chan_ngain)) ... iqr(lin_nGAM2.ecc(chan_ngain)-lin_n.ecc(chan_ngain))./2]; subplot(4,2,4);hold on; plot([0 15],[0 15]); scatter(lin_n.ecc(chan_pgain),lin_nGAM2.ecc(chan_pgain)); set(gca,'xlim',[0 15],'ylim',[0 15]); xlabel('B+ CHAN: Ecc HIGH GAM');ylabel('B+ CHAN: Ecc BETA'); fprintf('Are Positive BETA pRF shifted relative to fixation from High GAMMA?\n'); [p,h,stats] = signrank(lin_n.ecc(chan_pgain),lin_nGAM2.ecc(chan_pgain)); fprintf(['Median dEcc is ' ... num2str(median(lin_nGAM2.ecc(chan_pgain)-lin_n.ecc(chan_pgain))) ... ' IQR ' ... num2str(iqr(lin_nGAM2.ecc(chan_pgain)-lin_n.ecc(chan_pgain))./2) ... ' (POS --> towards fixation) \n']); fprintf(['p = ' num2str(p) '\n']); dE=[dE; ... median(lin_nGAM2.ecc(chan_pgain)-lin_n.ecc(chan_pgain)) ... iqr(lin_nGAM2.ecc(chan_pgain)-lin_n.ecc(chan_pgain))./2]; subplot(4,2,5:8);hold on; hold on; bar(1:4, dE(:,1)); errorbar(1:4, dE(:,1), dE(:,2), 'k','Linestyle', 'none'); ylabel('dEcc Beta vs Gam') set(gca,'xtick',1:4,'xticklabels',{'B- v lG','B+ v lG','B- v hG','B+ v hG'}) %% plot the diff in within channel pRF location for low/high LFP freq ===== fshiftprf = figure; set(fshiftprf,'Position',[100 100 1600 1600]); Xa1 = lin_n.X(chan_ngain); Ya1 = lin_n.Y(chan_ngain); Xlg1 = lin_nGAM1.X(chan_ngain); Ylg1 = lin_nGAM1.Y(chan_ngain); Xa2 = lin_n.X(chan_pgain); Ya2 = lin_n.Y(chan_pgain); Xlg2 = lin_nGAM1.X(chan_pgain); Ylg2 = lin_nGAM1.Y(chan_pgain); % only prfs within 8 deg for display csel1 = Xa1<8 & Ya1<2 & Xlg1<8 & Ylg1<2 & ... Xa1>-2 & Ya1>-8 & Xlg1>-2 & Ylg1>-8; csel2 = Xa2<8 & Ya2<2 & Xlg2<8 & Ylg2<2 & ... Xa2>-2 & Ya2>-8 & Xlg2>-2 & Ylg2>-8; subplot(2,2,1);hold on; plot([0 0],[-10 10],'--','LineWidth',2,'Color',[0.5 0.5 0.5]); plot([-10 10],[0 0],'--','LineWidth',2,'Color',[0.5 0.5 0.5]); quiver(Xlg1(csel1),Ylg1(csel1),Xa1(csel1)-Xlg1(csel1),Ya1(csel1)-Ylg1(csel1),... 'Color',[0.8 0.8 0.8],... 'AutoScale',0,'ShowArrowHead',0); % lines from gam to alpha pRF scatter(Xa1(csel1),Ya1(csel1),20,'ko','MarkerFaceColor','r'); % alpha pRF center scatter(Xlg1(csel1),Ylg1(csel1),20,'ko','MarkerFaceColor','b'); % lgam prf center % quiver(Xlg1(csel1),Ylg1(csel1),Xa1(csel1)-Xlg1(csel1),Ya1(csel1)-Ylg1(csel1),... % 'k','LineWidth',1,'AutoScale',1,'ShowArrowHead',1); % arrows from gam to alpha pRF plot(0,0,'go','MarkerSize',8,'MarkerFaceColor','g'); % fovea set(gca,'xlim',[-1.5 8.2],'ylim',[-6.5 2.2],'FontSize',14,... 'FontName','Myriad Pro','TickDir','out'); xlabel('Horizontal dva');ylabel('Vertical dva'); shiftsz = sqrt((Xa1(csel1)-Xlg1(csel1)).^2 + (Ya1(csel1)-Ylg1(csel1)).^2); m_ss = mean(shiftsz); se_ss = std(shiftsz)./sqrt(length(shiftsz)); fprintf(['Average negALPHA shift size: ' num2str(m_ss) ' +/- ' num2str(se_ss) ' SEM\n']); subplot(2,2,2);hold on; quiver(Xlg1(csel1),Ylg1(csel1),Xa1(csel1)-Xlg1(csel1),Ya1(csel1)-Ylg1(csel1),... 'k','AutoScale',1,'ShowArrowHead',1) plot(0,0,'go','MarkerSize',8,'MarkerFaceColor','g') set(gca,'xlim',[-1.5 8.2],'ylim',[-6.5 2.2],'FontSize',14,... 'FontName','Myriad Pro','TickDir','out'); xlabel('Horizontal dva');ylabel('Vertical dva'); subplot(2,2,3);hold on; plot([0 0],[-10 10],'--','LineWidth',2,'Color',[0.5 0.5 0.5]); plot([-10 10],[0 0],'--','LineWidth',2,'Color',[0.5 0.5 0.5]); quiver(Xlg2(csel2),Ylg2(csel2),Xa2(csel2)-Xlg2(csel2),Ya2(csel2)-Ylg2(csel2),... 'Color',[0.8 0.8 0.8],... 'AutoScale',0,'ShowArrowHead',0); % lines from gam to alpha pRF scatter(Xa2(csel2),Ya2(csel2),20,'ko','MarkerFaceColor','r'); % alpha pRF center scatter(Xlg2(csel2),Ylg2(csel2),20,'ko','MarkerFaceColor','b'); % lgam prf center % quiver(Xlg2(csel2),Ylg2(csel2),Xa2(csel2)-Xlg2(csel2),Ya2(csel2)-Ylg2(csel2),... % 'k','LineWidth',1,'AutoScale',1,'ShowArrowHead',1); % arrows from gam to alpha pRF plot(0,0,'go','MarkerSize',8,'MarkerFaceColor','g'); % fovea set(gca,'xlim',[-1.5 8.2],'ylim',[-6.5 2.2],'FontSize',14,... 'FontName','Myriad Pro','TickDir','out'); xlabel('Horizontal dva');ylabel('Vertical dva'); shiftsz = sqrt((Xa2(csel2)-Xlg2(csel2)).^2 + (Ya2(csel2)-Ylg2(csel2)).^2); m_ss = mean(shiftsz); se_ss = std(shiftsz)./sqrt(length(shiftsz)); fprintf(['Average posALPHA shift size: ' num2str(m_ss) ' +/- ' num2str(se_ss) ' SEM\n']); subplot(2,2,4);hold on; quiver(Xlg2(csel2),Ylg2(csel2),Xa2(csel2)-Xlg2(csel2),Ya2(csel2)-Ylg2(csel2),... 'k','AutoScale',1,'ShowArrowHead',1) plot(0,0,'go','MarkerSize',8,'MarkerFaceColor','g') set(gca,'xlim',[-1.5 8.2],'ylim',[-6.5 2.2],'FontSize',14,... 'FontName','Myriad Pro','TickDir','out'); xlabel('Horizontal dva');ylabel('Vertical dva'); if SaveFigs saveas(fextra,fullfile(pngfld, 'BETAvGAM_posneg_ecc_sz.png')); saveas(fextra,fullfile(svgfld, 'BETAvGAM_posneg_ecc_sz.svg')); end if CloseFigs; close(fextra); end %% Beta vs Gamma pRF separation index ===================================== lin_n = tLFP(... strcmp(tLFP.Model,'linear_ephys_cv1_neggain') & strcmp(tLFP.SigType,'Beta'),:); roi=1; % only do this for V1 channels chan_sel = lin_n.Area==roi & lin_n.R2>R2th & lin_n.R2>lin.R2+R2enh; chan_sel2 = lin_n.Area==roi & lin_n.R2>R2th; chan_pgain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain>0; chan_ngain = lin_n.Area==roi & lin_n.R2>R2th & lin_n.gain<0; chan_gam1 = lin_n.Area==roi & lin_nGAM1.R2>R2th; chan_gam2 = lin_n.Area==roi & lin_nGAM2.R2>R2th; chan_gam1 = lin_n.Area==roi & lin_nGAM1.R2>R2th; chan_gam2 = lin_n.Area==roi & lin_nGAM2.R2>R2th; [mean(sqrt((lin_n.X(chan_ngain)-lin_nGAM1.X(chan_ngain)).^2 + ... (lin_n.Y(chan_ngain)-lin_nGAM1.Y(chan_ngain)).^2)) ... std(sqrt((lin_n.X(chan_ngain)-lin_nGAM1.X(chan_ngain)).^2 + ... (lin_n.Y(chan_ngain)-lin_nGAM1.Y(chan_ngain)).^2))./sqrt(sum(chan_ngain))] [mean(sqrt((lin_n.X(chan_pgain)-lin_nGAM1.X(chan_pgain)).^2 + ... (lin_n.Y(chan_pgain)-lin_nGAM1.Y(chan_pgain)).^2)) ... std(sqrt((lin_n.X(chan_pgain)-lin_nGAM1.X(chan_pgain)).^2 + ... (lin_n.Y(chan_pgain)-lin_nGAM1.Y(chan_pgain)).^2))./sqrt(sum(chan_pgain))] [mean(sqrt((lin_n.X(chan_ngain)-lin_nGAM2.X(chan_ngain)).^2 + ... (lin_n.Y(chan_ngain)-lin_nGAM1.Y(chan_ngain)).^2)) ... std(sqrt((lin_n.X(chan_ngain)-lin_nGAM2.X(chan_ngain)).^2 + ... (lin_n.Y(chan_ngain)-lin_nGAM2.Y(chan_ngain)).^2))./sqrt(sum(chan_ngain))] [mean(sqrt((lin_n.X(chan_pgain)-lin_nGAM2.X(chan_pgain)).^2 + ... (lin_n.Y(chan_pgain)-lin_nGAM2.Y(chan_pgain)).^2)) ... std(sqrt((lin_n.X(chan_pgain)-lin_nGAM2.X(chan_pgain)).^2 + ... (lin_n.Y(chan_pgain)-lin_nGAM2.Y(chan_pgain)).^2))./sqrt(sum(chan_pgain))] fprintf('=== BETA ===\n') % size fprintf('Size diff (a-lG)\n') fprintf('NGAIN\n') [mean(lin_n.rfs(chan_ngain)-lin_nGAM1.rfs(chan_ngain)) ... std(lin_n.rfs(chan_ngain)-lin_nGAM1.rfs(chan_ngain))./sqrt(sum(chan_ngain))] fprintf('PGAIN\n') [mean(lin_n.rfs(chan_pgain)-lin_nGAM1.rfs(chan_pgain)) ... std(lin_n.rfs(chan_pgain)-lin_nGAM1.rfs(chan_pgain))./sqrt(sum(chan_pgain))] % normalized size fprintf('Norm Size diff (a/lG)\n') fprintf('NGAIN\n') [mean(lin_n.rfs(chan_ngain)./lin_nGAM1.rfs(chan_ngain)) ... std(lin_n.rfs(chan_ngain)./lin_nGAM1.rfs(chan_ngain))./sqrt(sum(chan_ngain))] fprintf('PGAIN\n') [mean(lin_n.rfs(chan_pgain)./lin_nGAM1.rfs(chan_pgain)) ... std(lin_n.rfs(chan_pgain)./lin_nGAM1.rfs(chan_pgain))./sqrt(sum(chan_pgain))] fprintf('Norm Size diff (a/hG)\n') fprintf('NGAIN\n') [mean(lin_n.rfs(chan_ngain)./lin_nGAM2.rfs(chan_ngain)) ... std(lin_n.rfs(chan_ngain)./lin_nGAM2.rfs(chan_ngain))./sqrt(sum(chan_ngain))] fprintf('PGAIN\n') [mean(lin_n.rfs(chan_pgain)./lin_nGAM2.rfs(chan_pgain)) ... std(lin_n.rfs(chan_pgain)./lin_nGAM2.rfs(chan_pgain))./sqrt(sum(chan_pgain))] % distance in relation to size distn = sqrt((lin_n.X(chan_ngain)-lin_nGAM1.X(chan_ngain)).^2 + ... (lin_n.Y(chan_ngain)-lin_nGAM1.Y(chan_ngain)).^2); distp = sqrt((lin_n.X(chan_pgain)-lin_nGAM1.X(chan_pgain)).^2 + ... (lin_n.Y(chan_pgain)-lin_nGAM1.Y(chan_pgain)).^2); sz2n=lin_n.rfs(chan_ngain)+lin_nGAM1.rfs(chan_ngain); sz2p=lin_n.rfs(chan_pgain)+lin_nGAM1.rfs(chan_pgain); fextra=figure; set(fextra,'Position',[10 10 600 1500]); subplot(4,2,1);hold on; scatter(distn,sz2n); scatter(distp,sz2p); plot([0 20],[0 20]); set(gca,'xlim',[0 20],'ylim',[0 20]); xlabel('distance B-G pRF');ylabel('sumsize B-G pRF') legend({'negative', 'positive'}); fprintf('DIST REL/Sz (dist/(s1+s2))\n') subplot(4,2,2);hold on; fprintf('lG NEG\n') fprintf('mean std median iqr\n') [mean(distn./sz2n) std(distn./sz2n)./sqrt(sum(chan_ngain))... median(distn./sz2n) iqr(distn./sz2n)./2] histogram(distn./sz2n,100,'Normalization','probability') fprintf('lG POS\n') [mean(distp./sz2p) std(distp./sz2p)./sqrt(sum(chan_pgain))... median(distp./sz2p) iqr(distp./sz2p)./2] histogram(distp./sz2p,100,'Normalization','probability') title('Low Gamma'); xlabel('Separation Index'); BC=[median(distn./sz2n) median(distp./sz2p)]; BCE=[iqr(distn./sz2n)/2 iqr(distp./sz2p)/2]; distn = sqrt((lin_n.X(chan_ngain)-lin_nGAM2.X(chan_ngain)).^2 + ... (lin_n.Y(chan_ngain)-lin_nGAM2.Y(chan_ngain)).^2); distp = sqrt((lin_n.X(chan_pgain)-lin_nGAM2.X(chan_pgain)).^2 + ... (lin_n.Y(chan_pgain)-lin_nGAM2.Y(chan_pgain)).^2); sz2n=lin_n.rfs(chan_ngain)+lin_nGAM2.rfs(chan_ngain); sz2p=lin_n.rfs(chan_pgain)+lin_nGAM2.rfs(chan_pgain); subplot(4,2,3);hold on; scatter(distn,sz2n); scatter(distp,sz2p); plot([0 20],[0 20]); xlabel('distance B-G pRF');ylabel('sumsize B-G pRF') set(gca,'xlim',[0 20],'ylim',[0 20]); legend({'negative', 'positive'}); subplot(4,2,4);hold on; fprintf('hG NEG\n') fprintf('mean std median iqr\n') [mean(distn./sz2n) std(distn./sz2n)./sqrt(sum(chan_ngain))... median(distn./sz2n) iqr(distn./sz2n)./2] histogram(distn./sz2n,100,'Normalization','probability') fprintf('hG POS\n') [mean(distp./sz2p) std(distp./sz2p)./sqrt(sum(chan_pgain))... median(distp./sz2p) iqr(distp./sz2p)./2] histogram(distp./sz2p,100,'Normalization','probability') title('High Gamma'); xlabel('Separation Index'); BC=[BC median(distn./sz2n) median(distp./sz2p)]; BCE=[BCE iqr(distn./sz2n)/2 iqr(distp./sz2p)/2]; subplot(4,2,5:8); hold on; bar(1:4, BC); ylabel('Separation Index') errorbar(1:4, BC, BCE, 'k', 'Linestyle','none') set(gca, 'xtick',1:4,'xticklabels',{'lG-','lG+','hG-','hG+'}) if SaveFigs saveas(fextra,fullfile(pngfld, 'BETAvGAM_sepidx.png')); saveas(fextra,fullfile(svgfld, 'BETAvGAM_sepidx.svg')); end if CloseFigs; close(fextra); end %% Beta DoG =============================================================== chan_sel = DoG.R2>R2th & DoG.R2>lin.R2+R2enh & ... DoG.normamp~=0 & DoG.ecc<16; MM=median(DoG.normamp(chan_sel)); fprintf(['BETA MEDIAN NAMP: ' num2str(MM) ', IQR ' num2str(iqr(DoG.normamp(chan_sel))) '\n']) % Wilcoxon 1-tailed < 1 [p,h,stats] = signrank(lin_n.gain(chan_sel),0,'tail','left'); fprintf(['Gain < 0: Wilcoxon z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); bb = [DoG.ecc(chan_sel) lin.ecc(chan_sel)]; % Wilcoxon [p,h,stats] = signrank(bb(:,1),bb(:,2)); fprintf(['ECC diff Wilcoxon z = ' ... num2str(stats.zval) ', p = ' num2str(p) '\n']); MM=median(bb(:,2)-bb(:,1)); fprintf(['BETA MEDIAN ECC DIFF: ' num2str(median(lin.ecc(chan_sel)-DoG.ecc(chan_sel))) ... ', IQR ' num2str(iqr(lin.ecc(chan_sel)-DoG.ecc(chan_sel))) '\n']) %% ======================================================================== % FMRI-EPHYS PRF ANALYSIS ------------------------------------------------ % ======================================================================== %% Fig 12 & SFig 12: Correlate MRI-ephys REGRESSION BASED ================= % based on the ecc vs size correlation % PER SIGNAL SOURCE % - calculate linear regression on random (filtered) subset of data % - repeat % - filter results % ACROSS SIGNAL SOURCES % - correlate the bootstrapped fit-result parameters % - optional >> do this in a bootstrapped way rng(1); % seed the random number generator % adjust these to check robustness ---------------------------------------- Rth_mri = 5; % R2 threshold MRI Rth_ephys = 50; % R2 threshold ephys selVFC = 1; % if true, only use lower right visual quadrant selANAT = 0; % if true, only use approximate area of electrodes ephysSUB = 'Both'; % ------------------------------------------------------------------------- MaxECC = [10 30]; % max ecc to use for fitting [V1 V4] MODS = {... 'linhrf_cv1_mhrf','linear_ephys_cv1';... 'linhrf_cv1_mhrf_neggain','linear_ephys_cv1_neggain';... 'csshrf_cv1_mhrf','css_ephys_cv1';... 'doghrf_cv1_mhrf','dog_ephys_cv1';... }; MRI_MODEL = MODS(:,1); EPHYS_MODEL = MODS(:,2); MMS={'linear','linear_ng','css','dog'}; warning off; cmROI = {'V1','V4'}; fprintf('=======================\n'); for m = 3 % select which model fprintf(['\nCrossmodal Correlation for Model: ' MODS{m} '\n']); s_R2 = T(modidx.(MRI_MODEL{m})).mod.R2 > Rth_mri & ... T(modidx.(MRI_MODEL{m})).mod.rfs < 1000; % collect mri prfs for r = 1:size(cmROI,2) if ~selANAT SSS = s_R2 & ismember( T(modidx.(MRI_MODELS{m})).mod.ROI,... ck_GetROIidx(cmROI(r),rois) ); else if strcmp(cmROI{r},'V1') SSS = s_R2 & T(modidx.(MRI_MODEL{m})).mod.ELEC_V1 > 0 & ... ismember( T(modidx.(MRI_MODELS{m})).mod.ROI,... ck_GetROIidx(cmROI(r),rois) ); elseif strcmp(cmROI{r},'V4') SSS = s_R2 & T(modidx.(MRI_MODEL{m})).mod.ELEC_V4 > 0 & ... ismember( T(modidx.(MRI_MODELS{m})).mod.ROI,... ck_GetROIidx(cmROI(r),rois) ); end end if selVFC VFSEL = T(modidx.(MRI_MODEL{m})).mod.X >= 0 & ... T(modidx.(MRI_MODEL{m})).mod.Y <=0; SSS = SSS & VFSEL; end % if selANAT && strcmp(cmROI{r},'V1') % ANATSEL = T(modidx.(MRI_MODEL{m})).mod.ELEC_V1 > 0; % SSS = SSS & ANATSEL; % elseif selANAT && strcmp(cmROI{r},'V4') % ANATSEL = T(modidx.(MRI_MODEL{m})).mod.ELEC_V4 > 0; % SSS = SSS & ANATSEL; % end if strcmp(cmROI{r},'V1') % V1 mri1(m).ECC = T(modidx.(MRI_MODEL{m})).mod.ecc(SSS); mri1(m).S = T(modidx.(MRI_MODEL{m})).mod.rfs(SSS); mri1(m).G = T(modidx.(MRI_MODEL{m})).mod.gain(SSS); elseif strcmp(cmROI{r},'V4') % V4 mri4(m).ECC = T(modidx.(MRI_MODEL{m})).mod.ecc(SSS); mri4(m).S = T(modidx.(MRI_MODEL{m})).mod.rfs(SSS); mri4(m).G = T(modidx.(MRI_MODEL{m})).mod.gain(SSS); end end % collect ephys prfs % MUA V1 if strcmp(ephysSUB,'M03') || strcmp(ephysSUB,'M04') s = strcmp(tMUA.Model,EPHYS_MODEL{m}) & ... strcmp(tMUA.Monkey,ephysSUB) & ... tMUA.Area == 1 & tMUA.R2 > Rth_ephys & tMUA.rfs < 1000; else s = strcmp(tMUA.Model,EPHYS_MODEL{m}) & ... tMUA.Area == 1 & tMUA.R2 > Rth_ephys & tMUA.rfs < 1000; end % Exclude 2 outlier V1 arrays in both monkeys s2 = ( strcmp(tMUA.Monkey,'M03') & (tMUA.Array == 11 | tMUA.Array == 13)) | ... ( strcmp(tMUA.Monkey,'M04') & (tMUA.Array == 10 | tMUA.Array == 12)); s(s2) = false; mua1(m).ECC = tMUA.ecc(s); mua1(m).S = tMUA.rfs(s); mua1(m).G = tMUA.gain(s); % MUA V4 s = strcmp(tMUA.Model,EPHYS_MODEL{m}) & ... tMUA.Area == 4 & tMUA.R2 > Rth_ephys & tMUA.rfs < 1000; mua4(m).ECC = tMUA.ecc(s); mua4(m).S = tMUA.rfs(s); mua4(m).G = tMUA.gain(s); % LFP freqband=unique(tLFP.SigType); for fb = 1: length(freqband) % V1 if strcmp(ephysSUB,'M03') || strcmp(ephysSUB,'M04') s = strcmp(tLFP.Model,EPHYS_MODEL{m}) & ... strcmp(tLFP.Monkey,ephysSUB) & ... tLFP.Area == 1 & tLFP.R2 > Rth_ephys & tLFP.rfs < 1000 & ... strcmp(tLFP.SigType, freqband{fb}); else s = strcmp(tLFP.Model,EPHYS_MODEL{m}) & ... tLFP.Area == 1 & tLFP.R2 > Rth_ephys & tLFP.rfs < 1000 & ... strcmp(tLFP.SigType, freqband{fb}); end s2 = ( strcmp(tMUA.Monkey,'M03') & (tMUA.Array == 11 | tMUA.Array == 13)) | ... ( strcmp(tMUA.Monkey,'M04') & (tMUA.Array == 10 | tMUA.Array == 12)); s(s2) = false; lfp1(fb,m).freqband = freqband{fb}; lfp1(fb,m).ECC = tLFP.ecc(s); lfp1(fb,m).S = tLFP.rfs(s); lfp1(fb,m).G = tLFP.gain(s); % V4 s = strcmp(tLFP.Model,EPHYS_MODEL{m}) & ... tLFP.Area == 4 & tLFP.R2 > Rth_ephys & tLFP.rfs < 1000 & ... strcmp(tLFP.SigType, freqband{fb}); lfp4(fb,m).freqband = freqband{fb}; lfp4(fb,m).ECC = tLFP.ecc(s); lfp4(fb,m).S = tLFP.rfs(s); lfp4(fb,m).G = tLFP.gain(s); end % Calculate & bootstrap linear regressions % ============================================= fprintf('Performing regressions on ALL data\n') MRI_stats1 = regstats(mri1(m).S,mri1(m).ECC); MRI_stats4 = regstats(mri4(m).S,mri4(m).ECC); MUA_stats1 = regstats(mua1(m).S,mua1(m).ECC); MUA_stats4 = regstats(mua4(m).S,mua4(m).ECC); for fb=1:length(freqband) try LFP_stats1{fb} = regstats(lfp1(fb,m).S,lfp1(fb,m).ECC); catch LFP_stats1{fb} = []; end try LFP_stats4{fb} = regstats(lfp4(fb,m).S,lfp4(fb,m).ECC); catch LFP_stats4{fb} = []; end end % Do statistics on model comparisons with fitlm ck_xmod_stats; XMOD(m).model = MMS{m}; XMOD(m).stats = stats; if m==2 ck_xmod_stats_lowfreqlfp; end end warning on; %% Fig 11C: Cross-signal comparison exponential parameter ================= RTHRES = 25; clear exptvals_MUA exptvals_LFP exptvals_MRI kw rois2 = [1 4]; f=figure;set(f,'Position',[10 10 1200 500]); ephysSub='M04'; for r=1:length(rois2) kw{r}=[]; exptvals_MRI{r} = [exptv(1).roi{rois2(r)};exptv(2).roi{rois2(r)}]; kw{r}=[kw{r}; exptvals_MRI{r} ones(size(exptvals_MRI{r}))]; if strcmp(ephysSub,'M03') || strcmp(ephysSub,'M04') exptvals_MUA{r} = tMUA.expt(strcmp(tMUA.Model,'css_ephys_cv1') & ... strcmp(tMUA.Monkey,ephysSub) & strcmp(tMUA.SigType,'MUA') & ... tMUA.Area==rois2(r) & tMUA.R2 > RTHRES ); fprintf(['Only ' ephysSub '\n']) else exptvals_MUA{r} = tMUA.expt(strcmp(tMUA.Model,'css_ephys_cv1') & ... strcmp(tMUA.SigType,'MUA') & tMUA.Area==rois2(r) & tMUA.R2 > RTHRES ); end kw{r}=[kw{r}; exptvals_MUA{r} 2*ones(size(exptvals_MUA{r}))]; sig=unique(tLFP.SigType); lfp_order = [3 1 2 5 4]; spn=1; cl=2; for fb=lfp_order cl=cl+1; if strcmp(ephysSub,'M03') || strcmp(ephysSub,'M04') exptvals_LFP{r,spn} = tLFP.expt(strcmp(tLFP.Model,'css_ephys_cv1') & ... strcmp(tLFP.Monkey,ephysSub) & strcmp(tLFP.SigType,sig{fb}) & ... tLFP.Area==rois2(r) & tLFP.R2 > RTHRES ); fprintf(['Only ' ephysSub '\n']) else exptvals_LFP{r,spn} = tLFP.expt(strcmp(tLFP.Model,'css_ephys_cv1') & ... strcmp(tLFP.SigType,sig{fb}) & tLFP.Area==rois2(r) & tLFP.R2 > RTHRES ); end if size(exptvals_LFP{r,spn},1)>3 kw{r}=[kw{r}; exptvals_LFP{r,spn} cl*ones(size(exptvals_LFP{r,spn}))]; end spn=spn+1; end figure(f); subplot(1,2,r); hold on; bar(1:7,[mean(exptvals_MRI{r}) ... mean(exptvals_MUA{r}) ... mean(exptvals_LFP{r,1}) ... mean(exptvals_LFP{r,2}) ... mean(exptvals_LFP{r,3}) ... mean(exptvals_LFP{r,4}) ... mean(exptvals_LFP{r,5})]); errorbar(1:7,[mean(exptvals_MRI{r}) ... mean(exptvals_MUA{r}) ... mean(exptvals_LFP{r,1}) ... mean(exptvals_LFP{r,2}) ... mean(exptvals_LFP{r,3}) ... mean(exptvals_LFP{r,4}) ... mean(exptvals_LFP{r,5})],... [std(exptvals_MRI{r}) ... std(exptvals_MUA{r}) ... std(exptvals_LFP{r,1}) ... std(exptvals_LFP{r,2}) ... std(exptvals_LFP{r,3}) ... std(exptvals_LFP{r,4}) ... std(exptvals_LFP{r,5})],'ko','Linestyle','none') set(gca,'xlim',[.5 7.5],'ylim',[0 1.1],'xtick',1:7,'xticklabels',... {'MRI','MUA', 'THETA','ALPHA','BETA','lGAM','hGAM'}) title(['V' num2str(rois2(r))]) ylabel('pRF exponent') fprintf('-- Compare < 1 --\n'); fprintf(['AREA V' num2str(rois2(r)) '\n']); [p,h,stats] = signrank(exptvals_MRI{r},1,'tail','left'); fprintf(['MRI, Wilcoxon EXPT < 1 : z = ' num2str(stats.zval) ', p = ' num2str(p) '\n']); [p,h,stats] = signrank(exptvals_MUA{r},1,'tail','left'); fprintf(['MUA, Wilcoxon EXPT < 1 : z = ' num2str(stats.zval) ', p = ' num2str(p) '\n']); for i=1:5 if size(exptvals_LFP{r,i},1) > 5 [p,h,stats] = signrank(exptvals_LFP{r,i},1,'tail','left'); if isfield(stats,'zval') fprintf(['LFP-' num2str(i) ', Wilcoxon EXPT < 1 : z = ' num2str(stats.zval) ', p = ' num2str(p) '\n']); end end end fprintf('-- Compare across signals --\n'); fprintf(['AREA V' num2str(rois2(r)) '\n']); [expcss(r).p,expcss(r).tbl,expcss(r).stats] = kruskalwallis(kw{r}(:,1), kw{r}(:,2)); [expcss(r).c,expcss(r).m,expcss(r).h,expcss(r).gnames] = multcompare(expcss(r).stats); end if SaveFigs saveas(f,fullfile(pngfld, 'xsig_prfexp.png')); saveas(f,fullfile(svgfld, 'xsig_prfexp.svg')); end if CloseFigs; close(f); end %% ======================================================================== % CLEAN UP --------------------------------------------------------------- % ======================================================================== %% clear and close close all; clc;