Microstim_FF-FB/SCRIPTS/FUNCTIONS/V1V1_spikes_ms1p.m

%% Population analysis of the straightforward MS data SPIKES ============
% no MS vs MS @ 150ms, visual stimulus: homogenous texture
% ------------------------------------------------------------
% Loads an overview table of data to potentially include, imports the 
% bandpass filtered data and saves it together in a population data file. 
% Once this has been done, it can also directly load this file and work 
% with that.

% single pulses

presaved = true; % 
set_spike_window = false; % allows to set the spike thresholding windows
MeanAlreadySubtracted = true;

%% load info ============================================================
load(fullfile(fld.data,'V1V1.mat'));
if ~set_spike_window 
    load(fullfile(fld.data,'V1V1_spikes.mat'));
end

%% load data ============================================================
fprintf('Loading data..');
if ~presaved
    curr_sel = 0;
    for ff=1:length(V1V1.Day)
        indstr = V1V1.indeces(ff);
        sep = find(indstr{1}=='_'==1);
        sel_files = ...
            [str2double(indstr{1}(1:sep-1)) str2double(indstr{1}(sep+1:end))];
        if curr_sel ~= sel_files(1) % only load files once
            curr_sel=sel_files(1);
            % Wav2 = CH1-16
            % Wave = CH17-32
            if MeanAlreadySubtracted
                % << artefact removal in raw
                cd(fullfile(fld.data,'BANDPASS','BP3','500to5000'));
            else
                % << first BP version (no artefact removal in raw)
                cd(fullfile(fld.data,'BANDPASS','BP','500to5000'));
            end
            
            CHMap = [1:16;17:32];
            CHi = [ones(1,16) ones(1,16).*2 ; 1:16 1:16];
            
            for df=1:2
                fprintf(['Loading ' filenames{sel_files(df)} '\n']);
                load(filenames{sel_files(df)},'Data')
                D(df)=Data; clear Data; %#ok<*SAGROW>
                t_vis = D(df).Time+0.150;
            end
        end
        % Resave selected channels
        % use 'CHi'
        sC = V1V1.rec_channel(ff);
        V1_V1(ff).info = V1V1(ff,:);
        V1_V1(ff).ms_trials = ...
            D(CHi(1,sC)).Trials(:,5)==1 & D(CHi(1,sC)).Trials(:,7)==1;
        V1_V1(ff).ms_i = (find(V1_V1(ff).ms_trials==1))';
        V1_V1(ff).noms_trials = ...
            D(CHi(1,sC)).Trials(:,5)==1 & D(CHi(1,sC)).Trials(:,7)==2;
        V1_V1(ff).noms_i = (find(V1_V1(ff).noms_trials==1))';
        V1_V1(ff).t_vis = t_vis;
        V1_V1(ff).d = D(CHi(1,sC)).Signal(:,CHi(2,sC),:);
    end
    cd(fld.scripts);
    V1_V1_info = V1V1;
    if MeanAlreadySubtracted
        save(fullfile(fld.data,'V1V1_Selected2_1p'),...
            'V1_V1','V1_V1_info','-v7.3');
    else
        save(fullfile(fld.data,'V1V1_Selected_1p'),...
            'V1_V1','V1_V1_info','-v7.3');
    end
else
    if MeanAlreadySubtracted
        load(fullfile(fld.data,'V1V1_Selected2_1p.mat'));
    else
        load(fullfile(fld.data,'V1V1_Selected_1p.mat'));
    end
end
fprintf('DONE!\n');

%% Set the spike thresholding window ====================================
if set_spike_window
    % select a channel
    sC = input(['Which data? (max ' num2str(length(V1_V1)) '): ']);
    msw = [0.148 0.152];
    MSmask = V1_V1(sC).t_vis>msw(1) & V1_V1(sC).t_vis<msw(2);
    
    figure
    subplot(2,1,1);hold on;
    for i=V1_V1(sC).noms_i
        plot(V1_V1(sC).t_vis,V1_V1(sC).d(:,:,i));
    end
    set(gca,'xlim',[-.05 0.3],'ylim',[-500 400])
    title(['All No MS Trials DataNum ' num2str(sC) ])
    subplot(2,1,2);hold on;
    for i=V1_V1(sC).ms_i
        V1_V1(sC).md(:,:,i)=V1_V1(sC).d(:,:,i);
        V1_V1(sC).md(MSmask,:,i)=0;
        plot(V1_V1(sC).t_vis,V1_V1(sC).md(:,:,i));
    end
    set(gca,'xlim',[-.05 0.3],'ylim',[-500 400])
    title(['All MS Trials DataNum ' num2str(sC) ])
    
    % set spike threshold window
    if isfield(V1_V1(sC),'TH') && ~isempty(V1_V1(sC).TH)
        TH = V1_V1(sC).TH;
    else
        V1_V1(sC).TH=[];
        while isempty(V1_V1(sC).TH)
            TH = input('Set spike threshold window: ');
            % plot the histograms
            PosNeg = -1; % neg peak detection
            spikes.ms=[];spikes.noms=[];
            n=0;
            for i=V1_V1(sC).noms_i(1:end)
                n=n+1;
                y = V1_V1(sC).d(:,:,i);
                [yp,xp] = findpeaks(double(PosNeg.*y),V1_V1(sC).t_vis,...
                    'MinPeakHeight',TH(1));
                xp(yp>TH(2))=[];
                yp(yp>TH(2))=[];
                spikes.noms=[spikes.noms; xp' yp ones(size(yp))*n];
            end
            n=0;
            for i=V1_V1(sC).ms_i(1:end)
                n=n+1;
                y = V1_V1(sC).md(:,:,i);
                [yp,xp] = findpeaks(double(PosNeg.*y),V1_V1(sC).t_vis,...
                    'MinPeakHeight',TH(1));
                xp(yp>TH(2))=[];
                yp(yp>TH(2))=[];
                spikes.ms=[spikes.ms; xp' yp ones(size(yp))*n];
            end
            
            BinSpecs=[-200 5 400]./1000;
            [spikes.noms_n, spikes.noms_bin] = ...
                hist(spikes.noms(:,1),BinSpecs(1):BinSpecs(2):BinSpecs(3));
            spikes.noms_n=spikes.noms_n./length(V1_V1(sC).noms_i);
            [spikes.ms_n, spikes.ms_bin] = hist(spikes.ms(:,1),...
                BinSpecs(1):BinSpecs(2):BinSpecs(3));
            spikes.ms_n=spikes.ms_n./length(V1_V1(sC).ms_i);
            
            figure;
            % Histograms -------
            subplot(2,3,4);hold on;
            bar(spikes.noms_bin(2:end-1),spikes.noms_n(2:end-1),'FaceColor','k','Barwidth',1)
            peak = max(spikes.noms_n(spikes.noms_bin>0 & spikes.noms_bin<0.150));
            plot([msw(1) msw(1)],[0 1000],'b--');
            plot([msw(2) msw(2)],[0 1000],'b--')
            set(gca,'xlim',[-0.1 0.350],'ylim',[0 peak*1.1])
            title('No MS');
            
            subplot(2,3,5);hold on;
            bar(spikes.ms_bin(2:end-1),spikes.ms_n(2:end-1),'FaceColor','r','EdgeColor','r','Barwidth',1)
            peak = max(spikes.ms_n(spikes.ms_bin>0 & spikes.ms_bin<0.150));
            plot([msw(1) msw(1)],[0 1000],'b--');
            plot([msw(2) msw(2)],[0 1000],'b--')
            set(gca,'xlim',[-0.1 0.350],'ylim',[0 peak*1.1])
            title('MS');
            
            % Raster plots ------
            subplot(2,3,1);hold on;
            for i=1:length(V1_V1(sC).noms_i)
                s=spikes.noms(:,3)==i;
                plot(spikes.noms(s,1),spikes.noms(s,3),'ko',...
                    'MarkerSize',3,'MarkerFaceColor','k')
            end
            plot([msw(1) msw(1)],[0 1000],'b--');
            plot([msw(2) msw(2)],[0 1000],'b--')
            set(gca,'xlim',[-0.1 0.35],'ylim',[0 i]);
            title('No MS');
            
            subplot(2,3,2);hold on
            for i=1:length(V1_V1(sC).ms_i)
                s=spikes.ms(:,3)==i;
                plot(spikes.ms(s,1),spikes.ms(s,3),'ro',...
                    'MarkerSize',3,'MarkerFaceColor','r')
            end
            plot([msw(1) msw(1)],[0 1000],'b--');
            plot([msw(2) msw(2)],[0 1000],'b--')
            set(gca,'xlim',[-0.1 0.35],'ylim',[0 i]);
            title('MS');
            
            % 2 histograms overlaid -----
            subplot(2,3,3);hold on;
            bar(spikes.noms_bin(2:end-1),spikes.noms_n(2:end-1),...
                'EdgeColor','k','LineWidth',2,'FaceColor',...
                'none','Barwidth',1)
            bar(spikes.ms_bin(2:end-1),spikes.ms_n(2:end-1),...
                'EdgeColor','r','LineWidth',2,'FaceColor',...
                'none','Barwidth',1)
            plot([msw(1) msw(1)],[0 1000],'b--');
            plot([msw(2) msw(2)],[0 1000],'b--')
            legend({'No MS','MS'})
            set(gca,'xlim',[-0.05 0.350],'ylim',[0 peak*1.1])
            title('No Ms vs. MS');
            
            % Difference zoomed in ------
            subplot(2,3,6);hold on;
            bar(spikes.noms_bin(2:end-1),...
                spikes.ms_n(2:end-1)-spikes.noms_n(2:end-1),...
                'EdgeColor','k','FaceColor','k','Barwidth',1)
            plot(spikes.noms_bin(2:end-1),...
                smooth(spikes.ms_n(2:end-1)-spikes.noms_n(2:end-1),5),...
                'r','LineWidth',1)
            plot([msw(1) msw(1)],[0 1000],'b--');
            plot([msw(2) msw(2)],[0 1000],'b--')
            set(gca,'xlim',[0.100 0.350],'ylim',[-0.08 0.3])
            title('Difference MS-NoMS');
            
            % ask to approve
            keepTH = input('Keep this spike window? (y/n): ','s');
            if strcmp(keepTH,'y') || strcmp(keepTH,'Y')
                V1_V1(sC).TH = TH;
                close all;
            end
        end
    end
else
    V1V1.TH1 = [V1V1_info{2:end,13}]';
    V1V1.TH2 = [V1V1_info{2:end,14}]';
    for sC = 1:size(V1_V1,2)
        V1_V1(sC).TH = [V1V1.TH1(sC) V1V1.TH2(sC)];
    end
end

%% Average ==============================================================
fprintf('Processing spikes\n');
if ~presaved
    TrialSel = [...
        3 10 40 10 40;...
        6 1 70 1 55;...
        10 70 215 70 215 ];
    warning off %#ok<*WNOFF>
    msw = [0.148 0.152];
    MaskStimArtefact = false;
    MaskStimBinOnly =true; % extrapolates previous timepoint
    
    BinSpecs=[-200 5 400]./1000;
    PosNeg = -1; % neg peak detection
    psth_nms=[];psth_ms=[];
    fprintf('Making PSTHs\n')
    for sC = 1:length(V1_V1)
        fprintf([num2str(sC) '\n']);
        spk(sC).ms=[];spk(sC).noms=[];
        
        if ismember(sC,TrialSel(:,1))
            r=find(TrialSel(:,1)==sC);
            nm_ti = TrialSel(r,2:3);
            m_ti = TrialSel(r,4:5);
        else
            nm_ti = [1 length(V1_V1(sC).noms_i)];
            m_ti = [1 length(V1_V1(sC).ms_i)];
        end
        n=0;
        
        % subtract average artefact shape ------
        MSTi = V1_V1(sC).ms_i(m_ti(1):m_ti(2));
        av_msTrace = nanmean(V1_V1(sC).d(:,:,MSTi),3);
        SubtractAverage=false;
        % --------------------------------------
        
        for i=V1_V1(sC).noms_i(nm_ti(1):nm_ti(2))
            n=n+1;
            y = V1_V1(sC).d(:,:,i);
            [yp,xp] = findpeaks(double(PosNeg.*y),V1_V1(sC).t_vis,...
                'MinPeakHeight',V1_V1(sC).TH(1));
            xp(yp>V1_V1(sC).TH(2))=[];
            yp(yp>V1_V1(sC).TH(2))=[];
            spk(sC).noms=[spk(sC).noms; xp' yp ones(size(yp))*n];
        end
        n=0;
        for i=V1_V1(sC).ms_i(m_ti(1):m_ti(2))
            n=n+1;
            if SubtractAverage
                y = V1_V1(sC).d(:,:,i)-av_msTrace;
                if MaskStimArtefact
                    MSmask = V1_V1(sC).t_vis>msw(1) & V1_V1(sC).t_vis<msw(2);
                    y(MSmask,:,:) = 0;
                end
            else
                y = V1_V1(sC).d(:,:,i); %use masked trace
            end
            [yp,xp] = findpeaks(double(PosNeg.*y),V1_V1(sC).t_vis,...
                'MinPeakHeight',V1_V1(sC).TH(1));
            xp(yp>V1_V1(sC).TH(2))=[];
            yp(yp>V1_V1(sC).TH(2))=[];
            % remove spikes during pulse window
            yp(xp>msw(1) & xp<msw(2))=[];
            xp(xp>msw(1) & xp<msw(2))=[];
            spk(sC).ms=[spk(sC).ms; xp' yp ones(size(yp))*n];
        end
                   
        % sliding window approach ----
        SlidWin.width = 0.010;
        SlidWin.step = 0.001;
        SlidWin.Min = -0.200;
        SlidWin.Max = 0.350;
        
        spk(sC).noms_n = []; spk(sC).noms_bin = [];
        for tstep = (SlidWin.Min-(SlidWin.width/2)):...
                SlidWin.step:(SlidWin.Max-(SlidWin.width/2))
            t1 = tstep; t2 = tstep+SlidWin.width;
            y=sum(spk(sC).noms(:,1) >= t1 & spk(sC).noms(:,1) <= t2)/...
                ( (t2-t1)*max(spk(sC).noms(:,3)));
            spk(sC).noms_bin=[spk(sC).noms_bin tstep];
            spk(sC).noms_n=[spk(sC).noms_n; y];
        end
        spk(sC).ms_n = []; spk(sC).ms_bin = [];
        for tstep = (SlidWin.Min-(SlidWin.width/2)):...
                SlidWin.step:(SlidWin.Max-(SlidWin.width/2))
            t1 = tstep; t2 = tstep+SlidWin.width;
            if t1>msw(1) && t1<=msw(2); t1 = msw(2); end
            if t2>=msw(1) && t2<msw(2); t2 = msw(1); end
            if t1<msw(1) && t2>=msw(2); t2 = msw(1); end
            if abs(t1-msw(1))<SlidWin.step 
                % don't change y
            else
                y=sum(spk(sC).ms(:,1) >= t1 & spk(sC).ms(:,1) <= t2)/...
                    ((t2-t1)*max(spk(sC).ms(:,3)));
            end
            spk(sC).ms_bin=[spk(sC).ms_bin tstep];
            spk(sC).ms_n=[spk(sC).ms_n; y];
        end 
        psth_nms=[psth_nms;spk(sC).noms_n']; %#ok<*AGROW>
        psth_ms=[psth_ms;spk(sC).ms_n'];
    end
    if MeanAlreadySubtracted
        save(fullfile(fld.data,'V1V1_spikes'),...
            'spk','psth_nms','psth_ms', 'V1V1');
    else
        save(fullfile(fld.data,'V1V1_spikes2'),...
            'spk','psth_nms','psth_ms', 'V1V1');
    end
    fprintf('DONE!\n')
else
    msw = [0.148 0.152];
    if MeanAlreadySubtracted
        load(fullfile(fld.data,'V1V1_spikes.mat'));
    else
        load(fullfile(fld.data,'V1V1_spikes2.mat'));
    end
    fprintf('DONE!\n')
end

%% subtract baseline, peak normalize & plot =============================
t=spk(1).ms_bin;
for s=1:size(spk,2)
    y1=spk(s).noms_n; y2=spk(s).ms_n; 
    y1_bl= y1-nanmean([y1(t>-0.100 & t<0); y2(t>-0.100 & t<0)]);
    y2_bl= y2-nanmean([y1(t>-0.100 & t<0); y2(t>-0.100 & t<0)]);
    yp = max(y1_bl(t>0 & t<0.100));
    y1n=y1_bl./yp; y2n=y2_bl./yp;
    spk(s).noms_n2=y1n; 
    spk(s).ms_n2=y2n; 
end

for s=1:size(spk,2)
    psth_nms(s,:) = spk(s).noms_n2; % Normalized
    %psth_nms(s,:) = spk(s).noms_n;
    psth_ms(s,:) = spk(s).ms_n2; % Normalized
    %psth_ms(s,:) = spk(s).ms_n;
end

%% plot the traces ======================================================
% select a subset of combinations
sel = V1V1.ManSel==1;
msw = [0.148 0.152];
% sel = V1V1.current>0; % all

N=sum(sel);
fprintf(['Number of stim-rec combinations in average: ' num2str(N) '\n']);

figure;
% Histograms -------
subplot(2,2,1);hold on;
BL_nms = mean2(psth_nms(sel,spk(1).noms_bin(2:end-1)<0 & ...
    spk(1).noms_bin(2:end-1)>-0.100));
mPSTH_nms = mean(psth_nms(sel,2:end-1),1);%-BL_nms;
sePSTH_nms = std(psth_nms(sel,2:end-1),1)./sqrt(length(sel));%-BL_nms;
bar(spk(1).noms_bin(2:end-1),mPSTH_nms,'FaceColor','k','Barwidth',1)
plot([msw(1) msw(1)],[-100 100],'b--');
plot([msw(2) msw(2)],[-100 100],'b--')
set(gca,'xlim',[-0.05 0.350],'ylim',[-0.1 1.1])
title('No MS');
text(0.25,0.8,['N = ' num2str(N)],'fontsize',12);

subplot(2,2,3);hold on;
BL_ms = mean2(psth_ms(sel,spk(1).ms_bin(2:end-1)<0 & ...
    spk(1).ms_bin(2:end-1)>-0.100));
mPSTH_ms = mean(psth_ms(sel,2:end-1),1);%-BL_ms;
sePSTH_ms = std(psth_ms(sel,2:end-1),1)./sqrt(length(sel));%-BL_nms;
bar(spk(1).ms_bin(2:end-1),mPSTH_ms,'FaceColor','r',...
    'EdgeColor','r','Barwidth',1)
plot([msw(1) msw(1)],[-100 100],'b--');
plot([msw(2) msw(2)],[-100 100],'b--')
set(gca,'xlim',[-0.05 0.350],'ylim',[-0.1 1.1])
title('MS');
text(0.25,0.8,['N = ' num2str(N)],'fontsize',12);

% 2 histograms overlaid -----
subplot(2,2,2);hold on;
plot(spk(1).noms_bin(2:end-1),mPSTH_nms,'Color','k','LineWidth',2)
plot(spk(1).ms_bin(2:end-1),mPSTH_ms,'Color','r','LineWidth',2)
plot([msw(1) msw(1)],[-100 100],'b--');
plot([msw(2) msw(2)],[-100 100],'b--')
legend({'No MS','MS'})
set(gca,'xlim',[-0.05 0.350],'ylim',[-0.1 1.1])
title('No Ms vs. MS');
text(0.25,0.8,['N = ' num2str(N)],'fontsize',12);

% Difference zoomed in ------
subplot(2,2,4);hold on;
bar(spk(1).noms_bin(2:end-1),mPSTH_ms-mPSTH_nms,...
    'EdgeColor','k','FaceColor','k','Barwidth',1)
plot([msw(1) msw(1)],[-100 100],'b--');
plot([msw(2) msw(2)],[-100 100],'b--')
set(gca,'xlim',[0.100 0.30],'ylim',[-0.20 0.3])
title('Difference MS-NoMS');
text(0.26,0.20,['N = ' num2str(N)],'fontsize',12);

print_figure('Fig2D' ,[0 0 32 20],fullfile(fld.figs,'FIG_2'),save_figs);
cd(fld.scripts);

%% sem-plot =============================================================
figure;
smw = 7;
subplot(1,3,1);hold on;
plot(spk(1).noms_bin(2:end-1),smooth(mPSTH_nms,smw),'k','LineWidth',2)
plot(spk(1).ms_bin(2:end-1),smooth(mPSTH_ms,smw),'r','LineWidth',2)
plot([-100 100],[0 0],'k');
plot([msw(1) msw(1)],[-100 100],'b--');
plot([msw(2) msw(2)],[-100 100],'b--')
legend({'No MS','MS'})
set(gca,'xlim',[-0.05 0.350],'ylim',[-.1 1.1])
title('No Ms vs. MS');
text(0.25,0.8,['N = ' num2str(N)],'fontsize',12);

subplot(1,3,2);hold on;
x=spk(1).noms_bin(2:end-1);
y=[mPSTH_nms-sePSTH_nms; mPSTH_nms; mPSTH_nms+sePSTH_nms];
px=[x,fliplr(x)];
py=[y(1,:), fliplr(y(3,:))];
patch(px,smooth(py,smw),1,'FaceColor','k','FaceAlpha',0.3,'EdgeColor','none');
plot(x,smooth(y(2,:),smw),'k');

x=spk(1).ms_bin(2:end-1);
y=[mPSTH_ms-sePSTH_ms; mPSTH_ms; mPSTH_ms+sePSTH_ms];
px=[x,fliplr(x)];
py=[y(1,:), fliplr(y(3,:))];
patch(px,smooth(py,smw),1,'FaceColor','r','FaceAlpha',0.3,'EdgeColor','none');
plot(x,smooth(y(2,:),smw),'r');
plot([-100 100],[0 0],'k')
plot([msw(1) msw(1)],[-100 100],'b--');
plot([msw(2) msw(2)],[-100 100],'b--')
set(gca,'xlim',[-0.05 0.350],'ylim',[-.1 1.1])
title('No Ms vs. MS (SEM)');
text(0.25,0.8,['N = ' num2str(N)],'fontsize',12);

subplot(1,3,3);hold on;
area(spk(1).noms_bin(2:end-1),smooth(mPSTH_ms-mPSTH_nms,smw), 'FaceColor','k')
plot([msw(1) msw(1)],[-100 100],'b--');
plot([msw(2) msw(2)],[-100 100],'b--')
set(gca,'xlim',[-0.1 0.350],'ylim',[-0.2 0.3])
title('Ms-NoMS');
text(0.25,0.25,['N = ' num2str(N)],'fontsize',12);

%% Statistics ===========================================================
ModWin = [0.150 0.160; 0.160 0.250];
t=spk(1).noms_bin(2:end-1);
MS = psth_ms(sel,2:end-1);
NMS = psth_nms(sel,2:end-1);
ExcWin = t>=ModWin(1,1) & t<=ModWin(1,2);
SupWin = t>=ModWin(2,1) & t<=ModWin(2,2);
ExcEffect = [mean(NMS(:,ExcWin),2) mean(MS(:,ExcWin),2)];
SupEffect = [mean(NMS(:,SupWin),2) mean(MS(:,SupWin),2)];

% 2-tailed
[H_exc_2t,p_exc_2t, ci_exc_2t, stats_exc_2t] = ...
    ttest(ExcEffect(:,1),ExcEffect(:,2),'tail','both');
[H_sup_2t,p_sup_2t, ci_sup_2t, stats_sup_2t] = ...
    ttest(SupEffect(:,1),SupEffect(:,2),'tail','both');

% 1-tailed >> test specifically excitation and suppression
[H_exc_1t,p_exc_1t, ci_exc_1t, stats_exc_1t] = ...
    ttest(ExcEffect(:,1),ExcEffect(:,2),'tail','left');
[H_sup_1t,p_sup_1t, ci_sup_1t, stats_sup_1t] = ...
    ttest(SupEffect(:,1),SupEffect(:,2),'tail','right');

figure;
subplot(1,2,1);hold on;
for s=1:size(ExcEffect,1)
    plot(1:2,ExcEffect(s,:),'o-','Markersize',6);
end
errorbar(1:2,mean(ExcEffect),std(ExcEffect)./sqrt(size(ExcEffect,1)),...
    'ko-','MarkerSize',10,'MarkerFaceColor','k','LineWidth',2)
set(gca,'xtick',1:2,'xticklabel',{'No MS','MS'},'xlim',[.5 2.5]);
title(['Excitatory Effect of Microstim (' num2str( 1000*(ModWin(1,1)-0.150) ) '-' ...
    num2str( 1000*(ModWin(1,2)-0.150)) ' ms)']);

subplot(1,2,2);hold on;
for s=1:size(SupEffect,1)
    plot(1:2,SupEffect(s,:),'o-','Markersize',6);
end
errorbar(1:2,mean(SupEffect),std(SupEffect)./sqrt(size(SupEffect,1)),...
    'ko-','MarkerSize',10,'MarkerFaceColor','k','LineWidth',2)
set(gca,'xtick',1:2,'xticklabel',{'No MS','MS'},'xlim',[.5 2.5]);
title(['Suppressive Effect of Microstim (' num2str( 1000*(ModWin(2,1)-0.150) ) '-' ...
    num2str( 1000*(ModWin(2,2)-0.150)) ' ms)']);

fprintf(['\nMicrostimulation induced excitation:\n' ...
    't(' num2str(stats_exc_1t.df) ') = ' num2str(stats_exc_1t.tstat) ...
    ', p = ' num2str(p_exc_1t) '\n']);
fprintf(['Microstimulation induced suppression:\n'...
    't(' num2str(stats_sup_1t.df) ') = ' num2str(stats_sup_1t.tstat) ...
    ', p = ' num2str(p_sup_1t) '\n']);

%% Example plot for FIGURE 2A-C =========================================
Check_IndV1;