%% Clear Workspace and go to folder clc; clear; close all figures = fullfile(cd,'figures'); set(0, 'DefaultFigureRenderer', 'painters'); % Save figures in vector format load('PupilData2.mat'); % Load the data % exclude subject 9 because of lack of data (only single session) Pupil(9,:) = []; for j = 1: size(Pupil,1); Pupil.AgeYears(j) = years(Pupil.T1_date(j)-Pupil.DOB(j)); end % Add age %% Analyze sessions close all CSAlone_Sessionavg = cell(size(Pupil,1),12); Paired_Sessionavg = cell(size(Pupil,1),12); for j = 1: size(Pupil,1) Subject = j; for i = 1: Pupil.Sessions(Subject) % 12= Max number of sessions Session = ['T',num2str(i)]; Time = Pupil.([Session '_Timestamps']){Subject}; % Timestamps Events = Pupil.([Session '_triggers']){Subject}; % Trigger annotations Events_Time = Pupil.([Session '_triggers_Timestamps']){Subject};% Timestamps for triggers dat.RD = Pupil.([Session '_Pupil_Right']){Subject}; % Size of Right Pupil dat.LD = Pupil.([Session '_Pupil_Left']){Subject}; % Size of Left Pupil % check and remove unrealistically small or large pupil sizes vars = {'RD','LD'}; for v=1:length(vars) qRem = dat.(vars{v})<1 | dat.(vars{v})>7; dat.(vars{v})(qRem) = nan; if any(qRem(:)) fprintf('subject %d, session %s, RD: %d/%d samples removed\n',Subject,Session,sum(qRem(:)),numel(dat.(vars{v}))); end end % find and remove blinks for v=1:length(vars) [on,off] = bool2bounds(isnan(dat.(vars{v}))); if 0 % debug plot clf plot(dat.(vars{v})) hold on plot(on -1,dat.(vars{v})( on-1),'go') plot(off+1,dat.(vars{v})(min(off+1,end)),'r*') end % remove 5 samples before and after missing to deal with blink % artifacts on = max( on-5,1); off = min(off+5,length(dat.(vars{v}))); qRem= bounds2bool(on,off,length(dat.(vars{v}))); dat.(vars{v})(qRem) = nan; if 0 % debug plot plot(dat.(vars{v}),'r') end end LD = dat.LD; RD = dat.RD; % Find CS only and paired LD_Paired = nan(75,300); RD_Paired = nan(75,300); LD_CSAlone = nan(25,300); RD_CSAlone = nan(25,300); x=1; y=1; qSkipNext = false; for k=1:numel(Events)-1 if ismember(char(Events(k)),{'TRIAL 1','TRIAL 2'}) qSkipNext = true; % skip the CS alone trials at start end if strcmp(char(Events(k)),'SOUND ON (duration 1.000 s)') && ~strcmp(char(Events(k+1)),'BRIGHT ON') % CS alone criteria if ~qSkipNext qTime= Time > Events_Time(k)-1000000 & Time < Events_Time(k) + 4000000; temp = LD(qTime); LD_CSAlone(x,1:min(numel(temp),300)) = temp(1:min(end,300)); temp = RD(qTime); RD_CSAlone(x,1:min(numel(temp),300)) = temp(1:min(end,300)); x=x+1; end qSkipNext = false; elseif strcmp(char(Events(k)),'SOUND ON (duration 1.000 s)') && strcmp(char(Events(k+1)),'BRIGHT ON') % Paired criteria qTime= Time > Events_Time(k)-1000000 & Time < Events_Time(k) + 4000000; temp = LD(qTime); LD_Paired(y,1:min(numel(temp),300)) = temp(1:min(end,300)); temp = RD(qTime); RD_Paired(y,1:min(numel(temp),300)) = temp(1:min(end,300)); y=y+1; end end assert(x==26) % 25 trials expected assert(y==76) % 75 trials expected % baseline correction (subtract median of first 60 samples (1s % before CS onset till CS onset) RD_CSAlone = RD_CSAlone-median(RD_CSAlone(:,1:60),2,'omitnan'); LD_CSAlone = LD_CSAlone-median(LD_CSAlone(:,1:60),2,'omitnan'); RD_Paired = RD_Paired -median(RD_Paired (:,1:60),2,'omitnan'); LD_Paired = LD_Paired -median(LD_Paired (:,1:60),2,'omitnan'); CSAlone_Sessionavg{j,i} = (RD_CSAlone + LD_CSAlone)/2; Paired_Sessionavg {j,i} = (RD_Paired + LD_Paired )/2; end end % find minimum pupil value for each trial during interval (0--2 s w.r.t. % CS) [CSAlone_all,Paired_all] = deal(cell(size(CSAlone_Sessionavg))); [CSAlone,Paired] = deal(nan(size(CSAlone_Sessionavg))); qProc = ~cellfun(@isempty,CSAlone_Sessionavg); CSAlone_all(qProc) = cellfun(@(x) min(x(:,61:180),[],2),CSAlone_Sessionavg(qProc),'uni',false); Paired_all(qProc) = cellfun(@(x) min(x(:,61:180),[],2),Paired_Sessionavg(qProc),'uni',false); CSAlone(qProc) = cellfun(@(x) mean(min(x(:,61:180),[],2),'omitnan'),CSAlone_Sessionavg(qProc)); Paired(qProc) = cellfun(@(x) mean(min(x(:,61:180),[],2),'omitnan'),Paired_Sessionavg(qProc)); % average traces CSAlone_Sessionavg_mean = cellfun(@(x) mean(x,'omitnan'),CSAlone_Sessionavg,'uni',false); Paired_Sessionavg_mean = cellfun(@(x) mean(x,'omitnan'), Paired_Sessionavg,'uni',false); % RMS % per session RMS_Session = cellfun(@(x) mean(sqrt(mean(diff(x,1,2).^2,'omitnan')),'omitnan'),CSAlone_Sessionavg); RMS_Subject = mean(RMS_Session,2,'omitnan'); fprintf('min RMS: %.4fmm, max RMS: %.4fmm\n',min(RMS_Subject),max(RMS_Subject)); %% Stats First = CSAlone(:,1); % Get Data from first session for j = 1: size(Pupil,1) Last(j,1) = CSAlone(j,find(~isnan(CSAlone(j,:)),1,'last')); % Get data from last session end jbtest(Last-First); % Test normality [h,p,CI, stats] = ttest(First-Last); fprintf('Mean First= %.5f +/- %.5f\n',mean(First,'omitnan'), std(First,'omitnan')); % Print out the mean for group A fprintf('Mean Last= %.5f +/- %.5f\n',mean( Last,'omitnan'), std( Last,'omitnan')); % Print out the mean for group B fprintf('Increase= %.5f\n',mean(Last,'omitnan')-mean(First,'omitnan')); fprintf('t = %.5f df = %d, p = %.5f\n',stats.tstat,stats.df,p); % Print out t-value, df, and p-value A = First; AX = ones(numel(A),1); B = Last; BX = ones(numel(B),1)*2; Box = [A; B]; X = [AX; BX]; lab = {'First', 'Last'}; % Labels for boxplot boxplot(Box',X',"Colors",'br',"Labels",lab); hold on % Make the boxplot plot([AX BX].',[A B].','-o','linewidth',0.1,'Color', [85 85 85]/255) ylabel('Amplitude (difference from baseline, mm)') % Set y-label %ylim([2 6]) set(gca, 'TickDir', 'out'); box off; set(gca,'color','none'); % Set tics out and remove white box str = ['Pupilsize, Before After']; title(str); % Add title %% Plots %% Plot average traces Trace_First = cat(1,CSAlone_Sessionavg_mean{:,1}); Paired_Trace_First = cat(1,Paired_Sessionavg_mean{:,1}); for j = 1:size(Pupil,1) Trace_Last(j,1:300) = CSAlone_Sessionavg_mean{j,Pupil.Sessions(j)}; Paired_Trace_Last(j,1:300) = Paired_Sessionavg_mean{j,Pupil.Sessions(j)}; end close subplot(2,2,1) errorbar(linspace(-1,4,300)',mean(Trace_First),std(Trace_First) / sqrt(10)); hold on errorbar(linspace(-1,4,300)',mean(Trace_Last),std(Trace_Last) / sqrt(10)); hold on errorbar(linspace(-1,4,300)',mean(Paired_Trace_First),std(Paired_Trace_First) / sqrt(10)); hold on errorbar(linspace(-1,4,300)',mean(Paired_Trace_Last),std(Paired_Trace_Last) / sqrt(10)); hold on %ylim([2 6]) xlim([-1 4]) set(gca,'YDir','normal') set(gca, 'TickDir', 'out') title('First and last session, average') set(gca, 'TickDir', 'out'); box off; set(gca,'color','none'); subplot(2,2,3) CSperiod = find(linspace(-1,4,300)'>0 & linspace(-1,4,300)'<1); X = linspace(-1,4,300)'; X=X(CSperiod); errorbar(X,mean(Trace_First(:,CSperiod)),std(Trace_First(:,CSperiod)) / sqrt(10)); hold on errorbar(X,mean(Trace_Last(:,CSperiod)),std(Trace_Last(:,CSperiod)) / sqrt(10)); hold on errorbar(X,mean(Paired_Trace_First(:,CSperiod)),std(Paired_Trace_First(:,CSperiod)) / sqrt(10)); hold on errorbar(X,mean(Paired_Trace_Last(:,CSperiod)),std(Paired_Trace_Last(:,CSperiod)) / sqrt(10)); hold on %ylim([4.2 5.1]) set(gca,'YDir','normal') set(gca, 'TickDir', 'out') title('First and last session, CS average') set(gca, 'TickDir', 'out'); box off; set(gca,'color','none'); str = ['Average traces']; sgtitle(str); % Add grand title cd(figures); saveas(gcf, str, 'pdf'); cd('..'); % Save the figure %% Heat figure() subplot(2,2,1) imagesc(Paired_Trace_First) colormap("summer") % Change colormap colorbar() % Add colorbar caxis([-2.5 0]) set(gca,'YDir','normal') set(gca, 'TickDir', 'out') subplot(2,2,2) imagesc(Trace_First) colormap("summer") % Change colormap colorbar() % Add colorbar caxis([-2.5 0]) set(gca,'YDir','normal') set(gca, 'TickDir', 'out') subplot(2,2,3) imagesc(Paired_Trace_Last) colormap("summer") % Change colormap colorbar() % Add colorbar caxis([-2.5 0]) set(gca,'YDir','normal') set(gca, 'TickDir', 'out') subplot(2,2,4) imagesc(Trace_Last) colormap("summer") % Change colormap colorbar() % Add colorbar caxis([-2.5 0]) set(gca,'YDir','normal') set(gca, 'TickDir', 'out') str = ['Heat First and last']; sgtitle(str); % Add grand title cd(figures); saveas(gcf, str, 'pdf'); cd('..'); % Save the figure %% Correlation clf errorbar(nanmean(CSAlone),nanstd(CSAlone) ./ sqrt(sum(~isnan(CSAlone))) ,'o') hold on errorbar(nanmean(Paired),nanstd(Paired) ./ sqrt(sum(~isnan(Paired))) ,'o') [RHO,PVAL] = corr(nanmean(CSAlone)',[1:12]'); % Calculate correlation [RHO,PVAL] = corr(nanmean(Paired)',[1:12]'); %ylim([2 6]) set(gca,'YDir','normal') set(gca, 'TickDir', 'out'); box off; set(gca,'color','none'); str = ['Minimal pupil size on probe trials']; title(str); % Add grand title cd(figures); saveas(gcf, str, 'pdf'); cd('..'); % Save the figure