Microstim_FF-FB/SCRIPTS/FUNCTIONS/FF_FIGURES.m

%% Generate the FEEDFORWARD MS figure panels ============================
fprintf('== GENERATING FEEDFORWARD MS EFFECT FIGURES ==================\n')

% Effect Windows
EffWin.V1toV4.Exc = [0 0.050];
EffWin.V1toV4.Sup = [0.050 0.120];
EffWin.V4toV1.Exc = [0 0.030];
EffWin.V4toV1.Sup = [0.030 0.150];

% FEEDFORWARD V1 to V4
EffectWindow.exc = EffWin.V1toV4.Exc;
EffectWindow.sup = EffWin.V1toV4.Sup;

fprintf(['Calculating V4 effect sizes in the following post-stimulation windows:\n'...
    'Excitation: ' num2str(EffectWindow.exc(1)*1e3) ' to ' num2str(EffectWindow.exc(2)*1e3) ' ms\n'...
    'Suppression: ' num2str(EffectWindow.sup(1)*1e3) ' to ' num2str(EffectWindow.sup(2)*1e3) ' ms\n']);

V1MS = preProcessV1MS(V1MS,EffectWindow);

%% FIGURE 3D - Example unit V1 ==========================================
fprintf('\n -- FIGURE 3D ----------\n');

% Set unit to plot as example. 
INDS = find(V1MS.RFDATA(:,5)>-1);
indCH = 31; soa = 1; S = 8;

% Normalized based on peak of blank (no MS) average. 
BL_EX = smooth(V1MS.CLEAN.UNITS(:,INDS(indCH),1),S );
MAX = nanmax(BL_EX); BL_EX = BL_EX/MAX;
MS_EX = smooth(V1MS.CLEAN.UNITS(:,INDS(indCH),soa+1),S )/MAX;

FSZ = 12;
figure (44);clf; hold all;
set(gca,'FontSize',FSZ)
xlim([-0.1 0.35]); ylim([-0.15 1.6]); YL = ylim;
p = patch([-0.05 -0.05 -0.05+0.025 -0.05+0.025],[YL(1) ...
    YL(2) YL(2) YL(1)],0.8*[1 1 1],'EdgeColor','none','FaceColor',[230 220 110]/255);
hms = plot(V1MS.TIME,smooth(MS_EX,1),'Color',V1MS.META.COLS(soa+1,:),'LineWidth',2);
hbl = plot(V1MS.TIME,smooth(BL_EX,1),'Color',V1MS.META.COLS(1,:)*0,'LineWidth',2);
legend([hbl hms],{'MS' 'BL'})
xlabel('Time [s]'); ylabel('MUA');
print_figure('figure3D',[0 0 7 7],fullfile(fld.figs,'FIG_3'),save_figs)
cd(fld.scripts);

%% FIGURE 3B - Make RF figure FF ========================================
fprintf('\n -- FIGURE 3B ----------\n');
load(fullfile(fld.data,'MAPPINGS','RF_B.mat'));

% Compute and plot V4 RF.
figure(18);
set(gca,'FontSize',FSZ)
V4Lims=[0 9; -9 0];
NV1RFs=[3]; % Say how many V1 RF there are.
v1ch = 3; % Say which one do you want to plot 
% (if you want only RF 1 and 3 for example you can use v1ch=[1 3] )
NV4CH=27; % Set the V4 channel that you want to plot.
YTicks=[-9:2:0];
XTicks=[0:2:9];

[X,Y]=meshgrid(V4Lims(1,1):V4Lims(1,2),V4Lims(2,1):V4Lims(2,2));
ST=0.01;
[XI,YI]=meshgrid([V4Lims(1,1):ST:V4Lims(1,2)],...
    [V4Lims(2,1):ST:V4Lims(2,2)]);
SM=3; count=0;
wh=[0.5 0.5]; POS=[0.5 1; 1 1; 0.5 0.5; 1 0.5];
count=count+1;

TEMP=squeeze(RESPMATV4(NV4CH,:,:));
if SM~=0
    TEMP3=TEMP;
    for j=1:size(TEMP,1)
        TEMP3(j,:)=smooth(TEMP(j,:),SM);
    end
    for i=1:size(TEMP3,2)
        TEMP3(:,i)=smooth(TEMP3(:,i),SM);
    end
else
    TEMP3=TEMP;
end

Z=interp2(X,Y,TEMP3,XI,YI);
colormap('default')
imagesc(V4Lims(1,:),V4Lims(2,:),Z/max(max(Z)));
set(gca,'YDir','normal')
hold on; axis('square')

% Plot v1 channel RF.
COLORS=[159/255, 29/255, 53/255; 0 128/255 0; ...
    125/255, 249/255, 1; 244/255, 196/255, 48/255];
for v1=v1ch
    TEMP=DATAV1(v1);
    
    DAT=TEMP.DAT;
    sampf=TEMP.Sampf;
    MP=TEMP.MP;
    DegperS=TEMP.DegperS;
    OFF.rho=TEMP.RFPecc;
    OFF.angle=TEMP.RFPangle;
    
    RF_INFO(v1)= RFpos(DAT, sampf, MP, DegperS, OFF);
    
    Mid=RF_INFO(v1).Mid;
    WH=RF_INFO(v1).WH;
    clr=RF_INFO(v1).clr;
    Ic=RF_INFO(v1).Ic;
    if v1==1;channs=[14];
    elseif v1==2;channs=[1 2 3 4 6 11 12 14 15];
    elseif v1==3;channs=[2];
    end
    for chan = channs
        M = [Mid.x(chan) Mid.y(chan)];
        P = patch([Mid.x(chan)-WH(chan).w2 Mid.x(chan)+WH(chan).w2 ...
            Mid.x(chan)+WH(chan).w2 Mid.x(chan)-WH(chan).w2], ...
            [Mid.y(chan)-WH(chan).h2 Mid.y(chan)-WH(chan).h2 ...
            Mid.y(chan)+WH(chan).h2 Mid.y(chan)+WH(chan).h2],COLORS(v1,:), ...
            'EdgeColor', 'w','LineWidth',4); 
        set(P,'FaceColor','none')
    end
end

set(gca,'XAxisLocation','top'); set(gca,'Color',[0 0 1]);
xlabel('Horizontal Position [deg]'); ylabel('Vertical Position [deg]');

print_figure('figure3B',[0 0 7 7],fullfile(fld.figs,'FIG_3'),save_figs)
cd(fld.scripts);

%% FIGURE 4A - Test Overlap RF ==========================================
fprintf('\n -- FIGURE 4A ----------\n');
OVER = V1MS.RFDATA(:,5); % Overlapp between V1 anv V4 RFs. 
EFFE  = 50*V1MS.EFFECT(1,:); % Effect size for ALL channels. 
AMP = V1MS.RFDATA(:,8); % MS current amplitude.

% Same as before but only for unit with a significant EFFECT and there is
% RF data for both V1 and V4 electrode.
OVER_SIG = OVER(logical(V1MS.CLEAN.UNITSIND) & ~isnan(OVER));
EFFE_SIG = EFFE(logical(V1MS.CLEAN.UNITSIND) & ~isnan(OVER));
AMP_SIG  = AMP(logical(V1MS.CLEAN.UNITSIND) & ~isnan(OVER));

% Compute correlation for current normalized effect size. 
eff = EFFE_SIG./AMP_SIG';
over = OVER_SIG;
iNAN = isnan(eff);
eff(iNAN) = [];
over(iNAN) = [];
X = [ones(length(eff),1) over];
y = eff';
[b,bint,r,rint,stats] = regress(y,X);

% Same as before with both significant and non significant units. 
OVER = OVER(~isnan(OVER));
EFFE = EFFE(~isnan(OVER));
AMP = AMP(~isnan(OVER));
eff = EFFE./AMP';
over = OVER;
iNAN = isnan(eff);
eff(iNAN) = [];
over(iNAN) = [];
X = [ones(length(eff),1) over];
y = eff';
[b,bint,r,rint,stats] = regress(y,X);

%Significante of the correlations. 
fprintf('\n-- Regression RF vs Exc Effect --\n');
[r,n,p,t] = regressionTTest(y,stats);
    disp(['R = ' num2str(r)])
    disp(['N = ' num2str(n)])
    disp(['P = ' num2str(p)])
    disp(['T = ' num2str(t)])

figure(132);clf;hold all
set(gca,'FontSize',FSZ)
scatter(OVER,EFFE./AMP',60,[0.75 0.75 0.9],'filled','LineWidth',1)
scatter(OVER_SIG,EFFE_SIG./AMP_SIG',60,V1MS.META.COLS(4,:),...
    'filled','LineWidth',1)

xlim([0 1]); ylim([-0.01 0.055]*50);
plot([0.1 0.97],b(2)*[0.15 0.99]+b(1),'Color',0*V1MS.META.COLS(2,:))
ylabel('Effect Size'); plot(xlim,[0 0],'k:','LineWidth',2);
xlabel('RFs Overlap');
print_figure('figure4A',[0 0 7 7],fullfile(fld.figs,'FIG_4'),save_figs)
cd(fld.scripts);

%% FIGURE S3A - SOA per Monkey ==========================================
fprintf('\n -- SUPP FIGURE 3A ----------\n');

S = 10; SS = [-0.05 0.05 .15];
monkeys = {'B' 'C'};
for M = 1:2
    figure (800+M);clf;
    iM = strcmp(V1MS.MONKEY,monkeys{M});
    temp = nanmean(nanmean(V1MS.CLEAN.UNITS(:,iM,:),3),1);
    ncount=sum(~isnan(temp));
    fprintf(['Monkey ' monkeys{M} ', n = ' num2str(ncount) '\n']);
    AVE =  squeeze(nanmean(V1MS.CLEAN.UNITS(:,iM,:),2));
    
    for soa = 1:3
        subplot(1,3,soa);hold all;
        set(gca,'FontSize',FSZ)
        xlim([-0.1 0.35])
        if strcmp(monkeys{M},'C')
            ylim([-1 3.3])
        else
            ylim([-0.2 1.7])
        end
        YL = ylim;
        p = patch([SS(soa) SS(soa) SS(soa)+0.025 ...
            SS(soa)+0.025],[YL(1) YL(2) YL(2) YL(1)],0.8*[1 1 1],...
            'EdgeColor','none','FaceColor',[230 220 110]/255);
        hms = plot(V1MS.TIME,smooth(AVE(:,soa+1)/nanmax(smooth(AVE(:,1),S)),S),...
            'Color',V1MS.META.COLS(soa+1,:), 'LineWidth',2);
        hbl = plot(V1MS.TIME,smooth(AVE(:,1),S)/nanmax(smooth(AVE(:,1),S)),...
            'Color',V1MS.META.COLS(1,:)*0,'LineWidth',2);
        legend([hbl hms],{'no MS' 'MS'}); xlabel('Time [s]');
        title(['SOA ' num2str(SS(soa)*1000) ' ms']);
    end
    print_figure(['suppfigure3A_Monkey' monkeys{M}],[0 0 7 7],...
        fullfile(fld.figs,'SUPPFIG_3'),save_figs);
    cd(fld.scripts);
end

%% FIGURE 3E - SOA both monkeys =========================================
fprintf('\n -- FIGURE 3E ----------\n');

S = 10; SS = [-0.05 0.05 .15];
figure(55);clf;
for soa = 1:3
    subplot(1,3,soa);hold all;
    set(gca,'FontSize',FSZ)
    xlim([-0.1 0.35])
    ylim([-0.7 3.4])
    YL = ylim; 
    p = patch([SS(soa) SS(soa) SS(soa)+0.025 SS(soa)+0.025],...
        [YL(1) YL(2) YL(2) YL(1)],0.8*[1 1 1],'EdgeColor','none','FaceColor',...
        [230 220 110]/255);
    hms = plot(V1MS.TIME,smooth(V1MS.CLEAN.AVE(:,soa+1)/...
        nanmax(smooth(V1MS.CLEAN.AVE(:,1),S)),S),...
        'Color',V1MS.META.COLS(soa+1,:), 'LineWidth',2);
    hbl = plot(V1MS.TIME,smooth(V1MS.CLEAN.AVE(:,1),S)/...
        nanmax(smooth(V1MS.CLEAN.AVE(:,1),S)),...
        'Color',V1MS.META.COLS(1,:)*0,'LineWidth',2);
    legend([hbl hms],{'no MS' 'MS'})
    title(['SOA ' num2str(SS(soa)*1000) ' ms']);
    xlabel('Time [s]')  
end
print_figure('figure3E' ,[0 0 7 7],fullfile(fld.figs,'FIG_3'),save_figs);
cd(fld.scripts);

%% FIGURE 3F - Difference for every SOA locked to MS ====================
fprintf('\n -- FIGURE 3F ----------\n');

figure(8);clf;hold all;
set(gca,'FontSize',FSZ)
ylim([-0.8 3.45])
YL = ylim;
p = patch([0 0 0.025 0.025],[YL(1) YL(2) YL(2) YL(1)],0.8*[1 1 1],...
    'EdgeColor','none','FaceColor',[230 220 110]/255);
S = 10; SAOS = [-0.05 0.05 0.15];
hds = [];
for soa = 1:3
    hds(end+1) = plot(V1MS.TIME-SAOS(soa),...
        smooth(-V1MS.CLEAN.AVE(:,1)+V1MS.CLEAN.AVE(:,soa+1),S)/...
        nanmax(smooth(V1MS.CLEAN.AVE(:,1),S)),...
        'Color',V1MS.META.COLS(soa+1,:),'LineWidth',2);
    xlim([-0.1 0.2])
end

plot(xlim,[0 0],'k:')
legend(hds,{'-50ms' '50ms' '150ms'},'Location','NorthEast')
xlabel('Time [s]')
print_figure('figure3F' ,[0 0 7 7],fullfile(fld.figs,'FIG_3'),save_figs);
    cd(fld.scripts);

%% FIGURE S4B - 25% of electrodes with smalles effect size at SOA 2 =====
fprintf('\n -- SUPP FIGURE 4B ----------\n');

eff = V1MS.CLEAN.EFFECT(2,:)';
amp = V1MS.RFDATA(:,8);
units = V1MS.CLEAN.UNITS;

iNotNan = ~isnan(eff) & ~isnan(amp);

eff = eff(iNotNan);amp = amp(iNotNan);units = units(:,iNotNan,:);
[v,indSort] = sort(abs(eff));

i25percent = 1:round(0.25*length(amp));

avUnits   = squeeze(nanmean(units,2));
av25units = squeeze(nanmean(units(:,indSort(i25percent),:),2));

subUnits = avUnits(:,2:4) - repmat(avUnits(:,1),1,3);
sub25units = av25units(:,2:4) - repmat(av25units(:,1),1,3);

figure(856);clf;hold all;
set(gca,'FontSize',FSZ)
ylim([-0.6 1.8])
YL = ylim;p = patch([0 0 0.025 0.025],[YL(1) YL(2) YL(2) YL(1)],...
    0.8*[1 1 1],'EdgeColor','none','FaceColor',[230 220 110]/255);
S = 3;
SAOS = [-0.05 0.05 0.15];
hds = [];
for soa = 1:3
    hds(end+1) = plot(V1MS.TIME-SAOS(soa),...
        smooth(-av25units(:,1)+av25units(:,soa+1),S)/...
        nanmax(smooth(av25units(:,1),S)),'Color',...
        V1MS.META.COLS(soa+1,:),'LineWidth',2);
    xlim([-0.1 0.3])
end

plot(xlim,[0 0],'k:')
legend(hds,{'-50ms' '50ms' '150ms'},'Location','NorthEast')
ylabel('MS - BL'); xlabel('Time [s]')

print_figure('suppfigure4B' ,[0 0 7 7],...
    fullfile(fld.figs,'SUPPFIG_4'),save_figs);
    cd(fld.scripts);

%% FIGURE S4A - MUA for 25% electrodes with smallest effect at SOA 2 ====
fprintf('\n -- SUPP FIGURE 4A ----------\n');

S = 10; SS = [-0.05 0.05 .15];
figure (soa+3737);clf;
for soa = 1:3
    subplot(1,3,soa);hold all;
    set(gca,'FontSize',FSZ)
    xlim([-0.1 0.35])
    ylim([-0.3 1.7])
    YL = ylim;p = patch([SS(soa) SS(soa) SS(soa)+0.025 SS(soa)+0.025],...
        [YL(1) YL(2) YL(2) YL(1)],0.8*[1 1 1],...
        'EdgeColor','none','FaceColor',[230 220 110]/255);
    hms = plot(V1MS.TIME,smooth(av25units(:,soa+1)/...
        nanmax(smooth(av25units(:,1),S)),S),'Color',...
        V1MS.META.COLS(soa+1,:), 'LineWidth',2);
    hbl = plot(V1MS.TIME,smooth(av25units(:,1),S)/...
        nanmax(smooth(av25units(:,1),S)),'Color',...
        V1MS.META.COLS(1,:)*0,'LineWidth',2);
    legend([hbl hms],{'no MS' 'MS'})
    xlabel('Time [s]')
    title(['SOA ' num2str(SS(soa)*1000) ' ms']);
end
print_figure('suppfigure4A' ,[0 0 7 7],...
    fullfile(fld.figs,'SUPPFIG_4'),save_figs);
cd(fld.scripts);

%% FIGURE S4C - Scatter plots for these conditions ======================
fprintf('\n -- SUPP FIGURE 4C ----------\n');

EFF = V1MS.CLEAN.EFFECT(:,iNotNan);
SS = [-0.05 0.05 .15];
COMBS = [1 2; 1 3; 2 3];
LIMS = [-0.3 1.5];
figure(178);clf;
plotCount = 0;
for comb = [1 3]
    plotCount = plotCount + 1;
    subplot(1,2,plotCount)
    hold all;
    scatter(V1MS.CLEAN.EFFECT(COMBS(comb,1),indSort((i25percent+1):end)),...
        V1MS.CLEAN.EFFECT(COMBS(comb,2),...
        indSort((i25percent+1):end)),10,[0 0 0],'filled');
    scatter(EFF(COMBS(comb,1),indSort(i25percent)),EFF(COMBS(comb,2),...
        indSort(i25percent)),10,[0.8 0.2 0.5],'filled');
    plot(LIMS,LIMS,':k');
    xlim(LIMS);ylim(LIMS);
    plot([0 0],ylim,':','Color',0.5*[1 1 1])
    plot(xlim,[0 0],':','Color',0.5*[1 1 1])
    xlabel(['SOA ' num2str(1000*SS(COMBS(comb,1))) ' ms']);
    ylabel(['SOA ' num2str(1000*SS(COMBS(comb,2))) ' ms']);
    axis square
end
print_figure('suppfigure4C' ,[0 0 7 7],...
    fullfile(fld.figs,'SUPPFIG_4'),save_figs);
cd(fld.scripts);

%% FIGURE 3G - Average Effect Size per SOA ==============================
fprintf('\n -- SUPP FIGURE 3G ----------\n');

ANOVA_DATA = [];ANOVA_DATA2 = [];

figure(103);clf;hold all
subplot(1,2,1);hold on;
set(gca,'FontSize',FSZ)
for soa = 1:3
    EFFE  = V1MS.EFFECT(soa,:);
    EFFE_SIG = EFFE(logical(V1MS.CLEAN.UNITSIND))/...
        nanmax(smooth(V1MS.CLEAN.AVE(:,1),S));
    mEFF    = nanmean(EFFE_SIG);
    sEFF    = nanstd(EFFE_SIG)/sqrt(sum(~isnan(EFFE_SIG)));
    errorbar(soa,mEFF,sEFF,'sk','LineWidth',2);
    bar(soa,mEFF,'FaceColor',V1MS.META.COLS(soa+1,:))
    ANOVA_DATA = [ANOVA_DATA EFFE_SIG'];
end
ylim([0 1.65]); xlim([0.5 3.5]);
set(gca,'XTick',[1 2 3],'XTickLabel',{'-20' '50' '150'})
ylabel('Mean Difference'); xlabel('SOA [ms]');
title('Excitation');

subplot(1,2,2);hold on;
set(gca,'FontSize',FSZ)
for soa = 1:3
    EFFE2  = V1MS.SUPPR_EFFECT(soa,:);
    EFFE_SIG2 = EFFE2(logical(V1MS.CLEAN.UNITSIND))/...
        nanmax(smooth(V1MS.CLEAN.AVE(:,1),S));
    mEFF2    = nanmean(EFFE_SIG2);
    sEFF2    = nanstd(EFFE_SIG2)/sqrt(sum(~isnan(EFFE_SIG2)));
    errorbar(soa,mEFF2,sEFF2,'sk','LineWidth',2);
    bar(soa,mEFF2,'FaceColor',V1MS.META.COLS(soa+1,:))
    ANOVA_DATA2 = [ANOVA_DATA2 EFFE_SIG2'];
end
ylim([-0.55 0.05])
xlim([0.5 3.5])
set(gca,'XTick',[1 2 3],'XTickLabel',{'-20' '50' '150'})
ylabel('Mean Difference');xlabel('SOA [ms]');
title('Suppression');

print_figure('figure3G' ,[0 0 7 7],fullfile(fld.figs,'FIG_3'),save_figs);
cd(fld.scripts);

%% STATS for FIG 3G =====================================================
fprintf('Statistics for FIG 3G\n');

fprintf('- Excitation -\n');
anova_data = ANOVA_DATA;
anova_data(isnan(ANOVA_DATA(:,1)),:) = [];
[P,ANOVATAB,STATS] = anova1(anova_data);
[c,m,h,names] = multcompare(STATS,'Alpha',10^-6,'CType','tukey-kramer');
for i=1:3
    fprintf(['SOA ' num2str(c(i,1)) ' vs ' num2str(c(i,2)) ': p = ' ...
        num2str(c(i,6)) '\n']);
end

fprintf('- Suppression -\n');
anova_data2 = ANOVA_DATA2;
anova_data2(isnan(ANOVA_DATA2(:,1)),:) = [];
[P2,ANOVATAB2,STATS2] = anova1(anova_data2);
[c2,m2,h2,names2] = multcompare(STATS2,'Alpha',10^-6,'CType','tukey-kramer');
for i=1:3
    fprintf(['SOA ' num2str(c2(i,1)) ' vs ' num2str(c2(i,2)) ': p = ' ...
        num2str(c2(i,6)) '\n']);
end
%% FIGURE 4B - Exc effect size vs MUA of Blank trials during stim. ======
fprintf('\n -- FIGURE 4B ----------\n');
S = 10; SS = [-0.05 0.05 .15];

% Compute activity during blank trials.
BASES = nan(size(V1MS.CLEAN.UNITS,2),3);
for soa = 1:3
    indT = V1MS.TIME > (SS(soa)+0) & V1MS.TIME < (SS(soa)+0.05);
    BASES(:,soa) = squeeze(nanmean(V1MS.UNITS(indT,:,1),1));
end

% Get effect size per soa per unit
EFFS = nan(size(V1MS.CLEAN.UNITS,2),3);
for soa = 1:3
    EFFE  = V1MS.EFFECT(soa,:);
    EFFS(:,soa) = EFFE;
end

figure(104);clf;hold on
set(gca,'FontSize',FSZ)

for soa = 3 % only the 150 ms SOA
    scatter(BASES(:,soa),EFFS(:,soa),60,[0.75 0.75 0.9],...
        'filled','LineWidth',1),
    scatter(BASES(logical(V1MS.CLEAN.UNITSIND),soa),...
        EFFS(logical(V1MS.CLEAN.UNITSIND),soa),60,...
        V1MS.META.COLS(soa+1,:),'filled','LineWidth',1);
    % compute regressions
    eff = EFFS(:,soa);
    over = BASES(:,soa);
    iNAN = isnan(eff);eff(iNAN) = [];over(iNAN) = [];
    iNAN = isnan(over);eff(iNAN) = [];over(iNAN) = [];
    X = [ones(length(eff),1) over];
    y = eff;
    [b,bint,r,rint,stats] = regress(y,X);
    disp(['Regression stats for SOA #' num2str(soa)])
    [r,n,p,t] = regressionTTest(y,stats);
    disp(['R = ' num2str(r)])
    disp(['N = ' num2str(n)])
    disp(['P = ' num2str(p)])
    disp(['T = ' num2str(t)])   
end

h = lsline;

ylim([-0.3 1.4]); xlim([-0.28 0.62]);
set(h(1),'Color','k','LineStyle','--')
set(h(2),'LineStyle','none')
xlabel('Average MUA'); ylabel('Effect Size')

print_figure('figure4B' ,[0 0 7 7],fullfile(fld.figs,'FIG_4'),save_figs);
cd(fld.scripts);