Transférer les fichiers vers 'scripts'

scripts/E_generateFigure2.m

@@ -0,0 +1,356 @@
+%% generate Figure 2.
+% Plotting Heatmaps all neurons 1st/last session early training and  extended training
+% need R_light datafiles, Celltype datafiles to only include MSN, and TRN
+% datafiles to restrict analysis on TRN.
+
+
+clc;clear;close all
+Task=1;
+load('Rextendedtraining_light.mat');
+load('Celltype_extendedTraining.mat');
+load('TRN_DT5_extendedtrain.mat')
+
+%% --- --- Plotting color-coded and average PSTHs --- ---
+BrainRegion=[10 20]; %10 for DLS, 20 for DMS
+ExcInh=1;        % 1 for excitations, 2 for inhibitions -- Used to find Exc or Inh onsets in R.Ev(.).Onsets
+
+Pstat=0.05;
+stepsize=0.3;      % Used for the scatter plot projections
+
+% For COLOR-CODED PSTH only (not Average PSTH)
+Norm=0;          % 0 = no norm - max accross all the selected neurons
+                 % 1 each neuron is normalized wh its own max
+                 % 2 for baseline substraction
+
+Sorting=4;       % 1 to sort neurons by onset
+                 % 2 to sort neurons by duration
+                 % 3 to sort neurons by PeakTime
+                 % 4 to sort by class of neuron
+%for RExt ONLY!  % 5 to sort by the magnitude of the cue response
+%for Cues ONLY!  % 6 to sort by the magnitude of the RExt response
+
+Xaxis=[0 357];Yaxis=[-1 3];
+Xaxis1=[-0.25 0.25];
+Ishow=find(Tm>=Xaxis1(1) & Tm<=Xaxis1(2));
+
+% ----------- COLORS ---------
+myblue=[0 113/255 188/255];
+mypurple=[200/255 0 255/255];
+mydarkblue=[0 0/255 255/255];
+TickSize=[0.015 0.02];
+
+Eventnames={'';'LI';'';'LP1';'';'LP2';'';'LP3';'';'LP4';'';'LP5';'';'PE'};
+
+time=-0.25:0.25:3.25;time2=1:1:357;
+
+%% Indiv Neuron
+%analysis restricted on putative MSN. nonTRN included
+
+%concat with normalized data
+DSconcat_OT = Ev(7).PSTHz(:,Ishow);
+for i=1:6
+    DSconcat_OT = cat(2,DSconcat_OT,Ev(i).PSTHz(:,Ishow));
+end
+
+DLSselection=Coord(:,4)==10 & Celltype(:,1)==1 & TRN(:,1)~=0;
+DMSselection=Coord(:,4)==20 & Celltype(:,1)==1 & TRN(:,1)~=0;
+
+LPwin=[52:280];
+binLP=mean(DSconcat_OT(:,LPwin),2);
+SEQwin=[26:331];
+binSEQ=mean(DSconcat_OT(TRN(:,1)~=0,SEQwin),2);
+TRN(TRN(:,1)~=0,2)=sign(binSEQ);
+
+DLSconcat_OT=DSconcat_OT(DLSselection,:);
+DMSconcat_OT=DSconcat_OT(DMSselection,:);
+
+%TMPdls=R.Class(DLSselection,1); TMPdms=R.Class(DMSselection,1);
+TMPdls=binLP(DLSselection,1); 
+TMPdms=binLP(DMSselection,1);
+TMPdls(isnan(TMPdls))=0;
+TMPdms(isnan(TMPdms))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
+[~,SORTimgDLS]=sort(TMPdls);[~,SORTimgDMS]=sort(TMPdms);
+
+DLSconcat_OT=DLSconcat_OT(SORTimgDLS,:);
+DMSconcat_OT=DMSconcat_OT(SORTimgDMS,:);
+
+
+subplot(9,8,[29 30 37 38 45 46 53 54 61 62 69 70])
+Clim=[-5 5];
+imagesc(time,[1 size(DLSconcat_OT,1)],DLSconcat_OT,Clim); %colorbar;axis tight;
+colormap(jet);
+hold on
+plot([0.25 0.25],[size(DLSconcat_OT,1)+0.5 0.5],'k')
+plot([0.75 0.75],[size(DLSconcat_OT,1)+0.5 0.5],'k')
+plot([1.25 1.25],[size(DLSconcat_OT,1)+0.5 0.5],'k')
+plot([1.75 1.75],[size(DLSconcat_OT,1)+0.5 0.5],'k')
+plot([2.25 2.25],[size(DLSconcat_OT,1)+0.5 0.5],'k')
+plot([2.75 2.75],[size(DLSconcat_OT,1)+0.5 0.5],'k')
+set(gca,'XTick',[0:0.25:3.25]);
+set(gca,'xticklabel',Eventnames)
+title('Extended training DLS');
+ylabel('neurons');
+hold off
+
+subplot(9,8,[31 32 39 40 47 48 55 56 63 64 71 72])
+imagesc(time,[1 size(DMSconcat_OT,1)],DMSconcat_OT,Clim); %colorbar;axis tight;
+colormap(jet);
+hold on
+plot([0.25 0.25],[size(DMSconcat_OT,1)+0.5 0.5],'k')
+plot([0.75 0.75],[size(DMSconcat_OT,1)+0.5 0.5],'k')
+plot([1.25 1.25],[size(DMSconcat_OT,1)+0.5 0.5],'k')
+plot([1.75 1.75],[size(DMSconcat_OT,1)+0.5 0.5],'k')
+plot([2.25 2.25],[size(DMSconcat_OT,1)+0.5 0.5],'k')
+plot([2.75 2.75],[size(DMSconcat_OT,1)+0.5 0.5],'k')
+set(gca,'XTick',[0:0.25:3.25]);
+set(gca,'xticklabel',Eventnames)
+title('Extended training DMS');
+ylabel('neurons');
+hold off
+
+
+%% Average OT PSTH
+
+subplot(9,8,[23 24])
+
+DLSpsth=nanmean(DLSconcat_OT,1);
+DLSsem=nanste(DLSconcat_OT,1);
+DLSup=DLSpsth+DLSsem;
+DLSdown=DLSpsth-DLSsem;
+hold on;
+patch([time2,time2(end:-1:1)],[DLSup,DLSdown(end:-1:1)],myblue,'EdgeColor','none');alpha(0.5);
+plot(time2,DLSpsth,'Color',myblue,'linewidth',1.5);
+
+DMSpsth=nanmean(DMSconcat_OT,1);
+DMSsem=nanste(DMSconcat_OT,1);
+DMSup=DMSpsth+DMSsem;
+DMSdown=DMSpsth-DMSsem;
+patch([time2,time2(end:-1:1)],[DMSup,DMSdown(end:-1:1)],'r','EdgeColor','none');alpha(0.5);
+plot(time2,DMSpsth,'r','linewidth',1.5);
+
+%plot([0 0], [-2 2.5],'LineStyle',':','Color','k');
+plot([0 357], [0 0],'LineStyle',':','Color','k');
+plot([51 51],Yaxis,'LineStyle',':','Color','k')
+plot([102 102],Yaxis,'LineStyle',':','Color','k')
+plot([153 153],Yaxis,'LineStyle',':','Color','k')
+plot([204 204],Yaxis,'LineStyle',':','Color','k')
+plot([255 255],Yaxis,'LineStyle',':','Color','k')
+plot([306 306],Yaxis,'LineStyle',':','Color','k')
+axis([Xaxis,Yaxis]);
+set(gca,'XTick',[0:25.5:357]);
+set(gca,'xticklabel',Eventnames)
+title('Extended training');
+ylabel('z-score');
+
+%% Pie Plot
+%proportion nonTRN excited inhibited
+%analysis restricted on putative MSN
+
+DLSselectionNonTRN=Coord(:,4)==10 & TRN(:,1)==0 & Celltype(:,1)==1;
+DLSselectionTRN_excited=Coord(:,4)==10 & TRN(:,1)~=0 & TRN(:,2)>0 & Celltype(:,1)==1;
+DLSselectionTRN_Inhibited=Coord(:,4)==10 & TRN(:,1)~=0 & TRN(:,2)<0 & Celltype(:,1)==1;
+
+DMSselectionNonTRN=Coord(:,4)==20 & TRN(:,1)==0 & Celltype(:,1)==1;
+DMSselectionTRN_excited=Coord(:,4)==20 & TRN(:,1)~=0 & TRN(:,2)>0 & Celltype(:,1)==1;
+DMSselectionTRN_Inhibited=Coord(:,4)==20 & TRN(:,1)~=0 & TRN(:,2)<0 & Celltype(:,1)==1;
+
+DLS_X(1)=sum(DLSselectionTRN_excited);
+DLS_X(2)=sum(DLSselectionTRN_Inhibited);
+DLS_X(3)=sum(DLSselectionNonTRN);
+
+DMS_X(1)=sum(DMSselectionTRN_excited);
+DMS_X(2)=sum(DMSselectionTRN_Inhibited);
+DMS_X(3)=sum(DMSselectionNonTRN);
+
+DLS_Percent=DLS_X*100./(sum(Coord(:,4)==10 & Celltype(:,1)==1));
+DMS_Percent=DMS_X*100./(sum(Coord(:,4)==20 & Celltype(:,1)==1));
+
+labels={'EXC','INH','non-TRN'};
+subplot(9,8,[11 12])
+bar(sum(DLS_Percent),'FaceColor',[0.9290 0.6940 0.1250]);
+hold on
+bar(sum(DLS_Percent(1:2)),'FaceColor',[0.8500 0.3250 0.0980]);
+bar(DLS_Percent(1),'FaceColor',[0 0.4470 0.7410]);
+axis([0.5 10.5 0 100 ])
+
+subplot(9,8,[19 20])
+bar(sum(DMS_Percent),'FaceColor',[0.9290 0.6940 0.1250]);
+hold on
+bar(sum(DMS_Percent(1:2)),'FaceColor',[0.8500 0.3250 0.0980]);
+bar(DMS_Percent(1),'FaceColor',[0 0.4470 0.7410]);
+axis([0.5 10.5 0 100 ])
+
+
+ save('Rextendedtraining_light.mat','Coord','Erefnames','Ev','Ninfo','Tm','TRN')
+
+%% ACQUISITION DATASET
+
+clear Ev 
+load('Rearlytraining_light.mat');
+load('Celltype_earlyTraining.mat');
+load('TRN_DT5_earlytrain.mat')
+
+for i=1:max(Celltype(:,1))
+    selectionIN=Ses(i).Coord(:,4)~=0;
+    selectsession=Celltype(:,1)==i;
+
+    Ses(i).Celltype(1:length(selectionIN),1)=0;
+    Ses(i).Celltype(selectionIN,1)=Celltype(selectsession,2);  
+end
+
+
+%analysis restricted on putative MSN. nonTRN included
+%concat with normalized data
+titlegraph={'Early training 1st session','Early training last session'};
+
+for j=1:2
+    vars={'TMPdls','TMPdms','binSEQ','binLP','DLSselection','DMSselection','DSconcat_ACQ','DLSconcat_ACQ','DMSconcat_ACQ'};
+    clear(vars{:});
+    
+    if j==1
+        k=4;
+    else
+        k=13;
+    end
+    
+    DSconcat_ACQ = Ses(k).Ev(7).PSTHz(:,Ishow);
+    for i=1:6
+        DSconcat_ACQ = cat(2,DSconcat_ACQ,Ses(k).Ev(i).PSTHz(:,Ishow));
+    end
+    
+    DLSselection=Ses(k).Coord(:,4)==10 & Ses(k).Celltype(:,1)==1 & sesTRN(k).TRN(:,1)~=0;
+    DMSselection=Ses(k).Coord(:,4)==20 & Ses(k).Celltype(:,1)==1 & sesTRN(k).TRN(:,1)~=0;
+    
+    LPwin=[52:280];
+    binLP=mean(DSconcat_ACQ(:,LPwin),2);
+    SEQwin=[26:331];
+    binSEQ=mean(DSconcat_ACQ(:,SEQwin),2);
+    sesTRN(k).TRN(:,2)=sign(binSEQ);
+    
+    DLSconcat_ACQ=DSconcat_ACQ(DLSselection,:);
+    DMSconcat_ACQ=DSconcat_ACQ(DMSselection,:);
+    
+    %TMPdls=R.Class(DLSselection,1); TMPdms=R.Class(DMSselection,1);
+    TMPdls=binLP(DLSselection,1);
+    TMPdms=binLP(DMSselection,1);
+    TMPdls(isnan(TMPdls))=0;
+    TMPdms(isnan(TMPdms))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
+    [~,SORTimgDLS]=sort(TMPdls);[~,SORTimgDMS]=sort(TMPdms);
+    
+    DLSconcat_ACQ=DLSconcat_ACQ(SORTimgDLS,:);
+    DMSconcat_ACQ=DMSconcat_ACQ(SORTimgDMS,:);
+    
+    
+    subplot(9,8,[25+(j-1)*24 26+(j-1)*24 33+(j-1)*24 34+(j-1)*24 41+(j-1)*24 42+(j-1)*24])
+    Clim=[-5 5];
+    imagesc(time,[1 size(DLSconcat_ACQ,1)],DLSconcat_ACQ,Clim); %colorbar;axis tight;
+    colormap(jet);
+    hold on
+    plot([0.25 0.25],[size(DLSconcat_ACQ,1)+0.5 0.5],'k')
+    plot([0.75 0.75],[size(DLSconcat_ACQ,1)+0.5 0.5],'k')
+    plot([1.25 1.25],[size(DLSconcat_ACQ,1)+0.5 0.5],'k')
+    plot([1.75 1.75],[size(DLSconcat_ACQ,1)+0.5 0.5],'k')
+    plot([2.25 2.25],[size(DLSconcat_ACQ,1)+0.5 0.5],'k')
+    plot([2.75 2.75],[size(DLSconcat_ACQ,1)+0.5 0.5],'k')
+    set(gca,'XTick',[0:0.25:3.25]);
+    set(gca,'xticklabel',Eventnames)
+    title('DLS');
+    ylabel('neurons');
+    hold off
+
+    subplot(9,8,[27+(j-1)*24 28+(j-1)*24 35+(j-1)*24 36+(j-1)*24 43+(j-1)*24 44+(j-1)*24])
+    imagesc(time,[1 size(DMSconcat_ACQ,1)],DMSconcat_ACQ,Clim); %colorbar;axis tight;
+    colormap(jet);
+    hold on
+    plot([0.25 0.25],[size(DMSconcat_ACQ,1)+0.5 0.5],'k')
+    plot([0.75 0.75],[size(DMSconcat_ACQ,1)+0.5 0.5],'k')
+    plot([1.25 1.25],[size(DMSconcat_ACQ,1)+0.5 0.5],'k')
+    plot([1.75 1.75],[size(DMSconcat_ACQ,1)+0.5 0.5],'k')
+    plot([2.25 2.25],[size(DMSconcat_ACQ,1)+0.5 0.5],'k')
+    plot([2.75 2.75],[size(DMSconcat_ACQ,1)+0.5 0.5],'k')
+    set(gca,'XTick',[0:0.25:3.25]);
+    set(gca,'xticklabel',Eventnames)
+    title('DMS');
+    ylabel('neurons');
+    hold off
+    
+    
+    %% Average OT PSTH
+    
+    subplot(9,8,[7+(j-1)*8 8+(j-1)*8])
+    
+    DLSpsth=nanmean(DLSconcat_ACQ,1);
+    DLSsem=nanste(DLSconcat_ACQ,1);
+    DLSup=DLSpsth+DLSsem;
+    DLSdown=DLSpsth-DLSsem;
+    hold on;
+    patch([time2,time2(end:-1:1)],[DLSup,DLSdown(end:-1:1)],myblue,'EdgeColor','none');alpha(0.5);
+    plot(time2,DLSpsth,'Color',myblue,'linewidth',1.5);
+    
+    DMSpsth=nanmean(DMSconcat_ACQ,1);
+    DMSsem=nanste(DMSconcat_ACQ,1);
+    DMSup=DMSpsth+DMSsem;
+    DMSdown=DMSpsth-DMSsem;
+    patch([time2,time2(end:-1:1)],[DMSup,DMSdown(end:-1:1)],'r','EdgeColor','none');alpha(0.5);
+    plot(time2,DMSpsth,'r','linewidth',1.5);
+    
+    %plot([0 0], [-2 2.5],'LineStyle',':','Color','k');
+    plot([0 357], [0 0],'LineStyle',':','Color','k');
+    plot([51 51],Yaxis,'LineStyle',':','Color','k')
+    plot([102 102],Yaxis,'LineStyle',':','Color','k')
+    plot([153 153],Yaxis,'LineStyle',':','Color','k')
+    plot([204 204],Yaxis,'LineStyle',':','Color','k')
+    plot([255 255],Yaxis,'LineStyle',':','Color','k')
+    plot([306 306],Yaxis,'LineStyle',':','Color','k')
+    axis([Xaxis,Yaxis]);
+    set(gca,'XTick',[0:25.5:357]);
+    set(gca,'xticklabel',Eventnames)
+    title(titlegraph(j));
+    ylabel('z-score');
+    
+end
+
+%% Pie Plot
+%proportion nonTRN excited inhibited
+%analysis restricted on putative MSN
+
+for i=4:13
+    
+    DSconcat_ACQ = Ses(i).Ev(7).PSTHz(:,Ishow);
+    for j=1:6
+        DSconcat_ACQ = cat(2,DSconcat_ACQ,Ses(i).Ev(j).PSTHz(:,Ishow));
+    end
+    SEQwin=[26:331];
+    binSEQ=nanmean(DSconcat_ACQ(:,SEQwin),2);
+    sesTRN(i).TRN(:,2)=sign(binSEQ);
+    
+    DLSselectionNonTRN_ACQ=Ses(i).Coord(:,4)==10 & sesTRN(i).TRN(:,1)==0 & Ses(i).Celltype(:,1)==1;
+    DLSselectionTRN_excited_ACQ=Ses(i).Coord(:,4)==10 & sesTRN(i).TRN(:,1)~=0 & sesTRN(i).TRN(:,2)>0 & Ses(i).Celltype(:,1)==1;
+    DLSselectionTRN_Inhibited_ACQ=Ses(i).Coord(:,4)==10 & sesTRN(i).TRN(:,1)~=0 & sesTRN(i).TRN(:,2)<0 & Ses(i).Celltype(:,1)==1;
+    
+    DMSselectionNonTRN_ACQ=Ses(i).Coord(:,4)==20 & sesTRN(i).TRN(:,1)==0 & Ses(i).Celltype(:,1)==1;
+    DMSselectionTRN_excited_ACQ=Ses(i).Coord(:,4)==20 & sesTRN(i).TRN(:,1)~=0 & sesTRN(i).TRN(:,2)>0 & Ses(i).Celltype(:,1)==1;
+    DMSselectionTRN_Inhibited_ACQ=Ses(i).Coord(:,4)==20 & sesTRN(i).TRN(:,1)~=0 & sesTRN(i).TRN(:,2)<0 & Ses(i).Celltype(:,1)==1;
+    
+    DLS_X_ACQ(i-3,1)=sum(DLSselectionTRN_excited_ACQ);
+    DLS_X_ACQ(i-3,2)=sum(DLSselectionTRN_Inhibited_ACQ);
+    DLS_X_ACQ(i-3,3)=sum(DLSselectionNonTRN_ACQ);
+    
+    DMS_X_ACQ(i-3,1)=sum(DMSselectionTRN_excited_ACQ);
+    DMS_X_ACQ(i-3,2)=sum(DMSselectionTRN_Inhibited_ACQ);
+    DMS_X_ACQ(i-3,3)=sum(DMSselectionNonTRN_ACQ);
+    
+    DLS_Percent_ACQ(i-3,:)=DLS_X_ACQ(i-3,:)*100./(sum(Ses(i).Coord(:,4)==10 & Ses(i).Celltype(:,1)==1));
+    DMS_Percent_ACQ(i-3,:)=DMS_X_ACQ(i-3,:)*100./(sum(Ses(i).Coord(:,4)==20 & Ses(i).Celltype(:,1)==1));
+end
+
+labels={'EXC','INH','non-TRN'};
+subplot(9,8,[9 10])
+bar(DLS_Percent_ACQ,1,'stacked');
+axis([0.5 10.5 0 100 ])
+subplot(9,8,[17 18])
+bar(DMS_Percent_ACQ,1,'stacked');
+axis([0.5 10.5 0 100 ])
+
+save('Rearlytraining_light.mat','Erefnames','Ses','Tm') 
+save('TRN_DT5_earlytrain.mat','sesTRN')

scripts/F_generateFigure3and5.m

@@ -0,0 +1,706 @@
+%% generate figure 3, 5 and Figure 3 - supplement figure 1
+
+% fourrier analysis - hierarchical clustering - Classification approach on
+% extended training dataset
+% statistics for the extended training dataset
+
+clc;clear;close all;
+load('Rextendedtraining_light.mat');
+load('RatID_extendedtraining.mat'); %% column 1 = rat ID, column 2 = habit ID; 0=goal-directed / 1=habit
+load('Celltype_extendedTraining.mat');
+load('TRN_DT5_extendedtrain.mat')
+
+timeconcat=[1:1:357];
+BrainRegion=[10 20];
+Erefnames=Erefnames(1:7);
+Xaxis1=[-0.25 0.25]; Yaxis=[-4 8];
+Xaxis2=[1 357]; Yaxis2=[-2 6];
+Ishow=find(Tm>=Xaxis1(1) & Tm<=Xaxis1(2));
+Eventnames={'';'LI';'';'LP1';'';'LP2';'';'LP3';'';'LP4';'';'LP5';'';'PE'};
+
+% ----------- COLORS ---------
+myblue=[0 113/255 188/255];
+mypurple=[200/255 0 255/255];
+mydarkblue=[0 0/255 255/255];
+TickSize=[0.015 0.02];
+Clim=[-5 5];
+
+
+%% preparation of data
+%normalization based on max
+for k=1:length(Erefnames)
+    Ev(k).PSTH_norm1=normalize(Ev(k).PSTHraw,Ev(k).PSTHrawBL,1);
+end
+
+%concat with normalized data
+DSconcat_OT = Ev(7).PSTH_norm1(:,Ishow);
+for i=1:6
+    DSconcat_OT = cat(2,DSconcat_OT,Ev(i).PSTH_norm1(:,Ishow));
+end
+
+% concat with Zscore
+DSconcat_OTz = Ev(7).PSTHz(:,Ishow);
+for i=1:6
+    DSconcat_OTz = cat(2,DSconcat_OTz,Ev(i).PSTHz(:,Ishow));
+end
+
+selection=Coord(:,4)~=0 & Celltype(:,1)==1;% & TRN(:,1)~=0; %exclude NaN values
+selectionID=find(Coord(:,4)~=0 & Celltype(:,1)==1);% & TRN(:,1)~=0);
+outID=find(Coord(:,4)==0 | Celltype(:,1)~=1);% & TRN(:,1)==0);
+DSconcat_OTshort=DSconcat_OTz(selection,:);
+
+
+%% Clustering on phasic and non phasic neurons 
+
+close all;
+f_low=1;f_high=4;
+
+%fourier analysis
+[X_all,f_all,P_all,Nt_all]=myfft(DSconcat_OTshort',0.01,1);
+total_power=sum(P_all);
+P_all_norm = bsxfun(@rdivide,P_all, total_power);
+
+
+a=find(f_all>f_low);b=find(f_all<f_high);
+lower_window_sum_all=sum(P_all_norm(1:a(1),:));
+upper_window_sum_all=sum(P_all_norm(a(1):b(end),:));
+
+ratio_all=lower_window_sum_all./upper_window_sum_all;
+
+%hierarchical clustering
+X_ratio=[[lower_window_sum_all]' [upper_window_sum_all]']; %#ok<NBRAK>
+eva_Ch_ratio = evalclusters(X_ratio,'linkage','CalinskiHarabasz','KList',[1:10]);
+c_ratio=eva_Ch_ratio.OptimalY; % if wat to use optimal number of cluster
+
+ Z = linkage(X_ratio,'ward'); % X is the data set, neuron in rows, features in columns
+ c = cluster(Z,'Maxclust',3); % determine the number cluster here, if you don't use optimalK
+ D = pdist(X_ratio);
+leafOrder = optimalleaforder(Z,D);
+cutoff = median([Z(end-2,3) Z(end-1,3)]);
+
+ phasicID(selectionID,1)=c_ratio;
+ phasicID(outID,1)=NaN;
+
+
+%% heatmap results clustering
+figure(1)
+
+for i=1:2
+    ind=find(phasicID==i);
+    t=num2str(i);
+    Clim=[-5 5];
+    DSconcat_OT_class=[];SORTimg=[];MaxVal=[];MeanVal=[]; MaxTime=[];TMP=[];NNid=[];
+    DSconcat_OT_class=DSconcat_OT(ind,:);
+    
+    for NN=1:size(DSconcat_OT_class,1)
+        MeanVal(NN,1)=mean(DSconcat_OT_class(NN,52:306),2);
+        [MaxVal(NN,1),MaxInd]=max(DSconcat_OT_class(NN,:));
+        MaxTime(NN,1)=timeconcat(timeconcat(1)+MaxInd-1);
+        if i==2
+            MaxTimePhasic=MaxTime;
+        end
+    end
+    if i==1
+        TMP=MeanVal;
+        subplot(5,3,[4 7 10 13])
+    elseif i==2 
+        TMP=MeanVal;
+        subplot(5,3,[5 8 11 14])
+    end
+    TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
+    [~,SORTimg]=sort(TMP);
+    %NNid=cell2mat(Ninfo(phasic_neuron_selection,3));NNid2=NNid(SORTimg);
+    DSconcat_OT_class=DSconcat_OT_class(SORTimg,:);
+    
+    
+    imagesc(timeconcat,[1 size(DSconcat_OT_class,1)],DSconcat_OT_class); %colorbar;axis tight;
+    colormap(jet);
+    hold on
+    plot([51 51],[size(DSconcat_OT_class,1)+0.5 0.5],'k')
+    plot([102 102],[size(DSconcat_OT_class,1)+0.5 0.5],'k')
+    plot([153 153],[size(DSconcat_OT_class,1)+0.5 0.5],'k')
+    plot([204 204],[size(DSconcat_OT_class,1)+0.5 0.5],'k')
+    plot([255 255],[size(DSconcat_OT_class,1)+0.5 0.5],'k')
+    plot([306 306],[size(DSconcat_OT_class,1)+0.5 0.5],'k')
+    set(gca,'XTick',[0:25.5:357]);
+    set(gca,'xticklabel',Eventnames)
+    hold off
+end
+
+%% heatmap phasic versus non phasic neurons
+
+phasicID(phasicID(:,1)==3,1)=1;
+subplot(5,3,2)
+dendrogram(Z,'ColorThreshold',2, 'Reorder',leafOrder);
+
+
+for j=1:2 % loop thru structure
+    for k=1:max(phasicID(:,1)) % loop thru class
+        DSselection(:,1)=Coord(:,4)==BrainRegion(j) & phasicID(:,1)==k;
+        DSselectionNumber(j,k)=sum(DSselection);
+        DSselectionProp(j,k)=100*sum(DSselection)/sum(Coord(:,4)==BrainRegion(j) & Celltype(:,1)==1);
+    end
+end
+subplot(5,3,3)
+X = categorical({'DLS','DMS'});
+b2=bar(X,DSselectionProp,'stacked');
+ylabel('% neurons');
+legend('phasic', 'nonphasic')
+
+
+%% Figure scatter plot Hierarchical clustering based on fourier
+
+for i=1:length(c_ratio)
+    if c_ratio(i,1)==2
+        color(i,:)=[255/255 0/255 0/255]; %phasic red
+    elseif c_ratio(i,1)==1
+        color(i,:)=[0 255/255 0/255]; %sustained green
+    end
+end
+subplot(5,3,1)
+scatter(lower_window_sum_all,upper_window_sum_all,[],color)
+hold on
+title(strcat('hierarchical clustering'))
+xlabel('Power low freq') % x-axis label
+ylabel('Power int freq')% y-axis label
+legend('phasic', 'nonphasic')
+hold off
+
+%% Plot heatmap phasic neurons, separation DMS/DLS, 
+PhasicNeurons_ID(:,1)=find(phasicID==2);
+PhasicNeurons_ID(:,2)=Coord(PhasicNeurons_ID(:,1),4);
+PhasicNeurons_ID(:,3)=RatID(PhasicNeurons_ID(:,1),1);
+
+DSconcat_OT_phasic=DSconcat_OTz(PhasicNeurons_ID(:,1),:);
+DLSselection=PhasicNeurons_ID(:,2)~=0 & PhasicNeurons_ID(:,2)==10;
+DMSselection=PhasicNeurons_ID(:,2)~=0 & PhasicNeurons_ID(:,2)==20;
+
+% sorting DLS
+DSconcat_OT_classDLS=DSconcat_OT_phasic(DLSselection,:);
+for NN=1:size(DSconcat_OT_classDLS,1)
+    [MaxValDLS(NN,1),MaxInd]=max(DSconcat_OT_classDLS(NN,:));
+    MaxTimeDLS(NN,1)=timeconcat(timeconcat(1)+MaxInd-1);
+end
+TMPdls=MaxTimeDLS;
+TMPdls(isnan(TMPdls))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
+[~,SORTimgDLS]=sort(TMPdls);
+DSconcat_OT_classDLS=DSconcat_OT_classDLS(SORTimgDLS,:);
+
+% sorting DMS
+DSconcat_OT_classDMS=DSconcat_OT_phasic(DMSselection,:);
+for NN=1:size(DSconcat_OT_classDMS,1)
+    [MaxValDMS(NN,1),MaxInd]=max(DSconcat_OT_classDMS(NN,:));
+    MaxTimeDMS(NN,1)=timeconcat(timeconcat(1)+MaxInd-1);
+end
+TMPdms=MaxTimeDMS;
+TMPdms(isnan(TMPdms))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
+[~,SORTimgDMS]=sort(TMPdms);
+DSconcat_OT_classDMS=DSconcat_OT_classDMS(SORTimgDMS,:);
+
+% plot heatmap
+figure(2)
+subplot(9,4,[1 2 5 6 9 10])
+imagesc(timeconcat,[1 size(DSconcat_OT_classDLS,1)],DSconcat_OT_classDLS,Clim); %colorbar;axis tight;
+colormap(jet);
+hold on
+plot([51 51],[size(DSconcat_OT_classDLS,1)+0.5 0.5],'k')
+plot([102 102],[size(DSconcat_OT_classDLS,1)+0.5 0.5],'k')
+plot([153 153],[size(DSconcat_OT_classDLS,1)+0.5 0.5],'k')
+plot([204 204],[size(DSconcat_OT_classDLS,1)+0.5 0.5],'k')
+plot([255 255],[size(DSconcat_OT_classDLS,1)+0.5 0.5],'k')
+plot([306 306],[size(DSconcat_OT_classDLS,1)+0.5 0.5],'k')
+set(gca,'XTick',[0:25.5:357]);
+set(gca,'xticklabel',Eventnames)
+ylabel('DLS neurons')
+hold off
+
+subplot(9,4,[17 18 21 22 25 26])
+imagesc(timeconcat,[1 size(DSconcat_OT_classDMS,1)],DSconcat_OT_classDMS,Clim); %colorbar;axis tight;
+colormap(jet);
+hold on
+plot([51 51],[size(DSconcat_OT_classDMS,1)+0.5 0.5],'k')
+plot([102 102],[size(DSconcat_OT_classDMS,1)+0.5 0.5],'k')
+plot([153 153],[size(DSconcat_OT_classDMS,1)+0.5 0.5],'k')
+plot([204 204],[size(DSconcat_OT_classDMS,1)+0.5 0.5],'k')
+plot([255 255],[size(DSconcat_OT_classDMS,1)+0.5 0.5],'k')
+plot([306 306],[size(DSconcat_OT_classDMS,1)+0.5 0.5],'k')
+set(gca,'XTick',[0:25.5:357]);
+set(gca,'xticklabel',Eventnames)
+ylabel('DMS neurons')
+hold off
+
+
+%% Hierarchical clustering of phasic neurons based on Sequence-related activity
+
+LIwin=[26:51];%post lever insertion window
+LPwin=[52:280];%pre-1stLP to pre-lastLP window
+PEwin=[307:332];% pre-PE window
+
+binLI=mean(DSconcat_OT(:,LIwin),2);
+binLP=mean(DSconcat_OT(:,LPwin),2);
+binPE=mean(DSconcat_OT(:,PEwin),2);
+binnedDS=cat(2,binLI,binLP,binPE);   
+
+sel=Coord(:,4)~=0 & Celltype==1 & phasicID==2;
+X=binnedDS(sel,:);
+
+Z = linkage(X,'ward');
+c = cluster(Z,'Maxclust',4);
+eva = evalclusters(X,'linkage','CalinskiHarabasz','KList',[1:10]); %% You can change the criterion here. There are 3 other options. Gap, 'Silhouette' , 'DaviesBouldin'. You can do help evalclusters to use this code.
+D = pdist(X);
+leafOrder = optimalleaforder(Z,D);
+cutoff = median([Z(end-2,3) Z(end-1,3)]);
+subplot(9,4,[29 33])
+dendrogram(Z,'ColorThreshold',cutoff, 'Reorder',leafOrder);
+
+ClassPHAS(:,1)=eva.OptimalY;
+ClassPHAS(:,2:4)=X;
+
+
+%% eva has the results of fitting the clusters. eva.OptimalK tells us how many clusters fit this data best based on the criterion used.
+%% eva.OptimalY tells us the cluster identities of each point.
+%% CH criterion clustering labels
+ind=[];
+DSconcat_OT_Phasic=DSconcat_OTz(sel,:);
+
+for i=1:max(eva.OptimalY)
+    DSconcat_OT_classDLS=[]; SORTimgDLS=[];MaxValDLS=[];TMPdls=[];
+    DSconcat_OT_classDMS=[]; SORTimgDMS=[];MaxValDMS=[];TMPdms=[];
+    ind=find(eva.OptimalY==i);
+    
+    PhasicNeurons_ID(:,1)=find(phasicID~=1 & ~isnan(phasicID));
+    PhasicNeurons_ID(:,2)=Coord(PhasicNeurons_ID(:,1),4);
+    PhasicNeurons_ID(:,3)=eva.OptimalY;
+    
+    DLSselectionP=PhasicNeurons_ID(:,2)==10 &  PhasicNeurons_ID(:,3)==i;
+    DMSselectionP=PhasicNeurons_ID(:,2)==20 &  PhasicNeurons_ID(:,3)==i;
+
+    
+    % sorting DLS
+    DSconcat_OT_classDLS=DSconcat_OT_Phasic(DLSselectionP,:);
+    for NN=1:size(DSconcat_OT_classDLS,1)
+        MaxValDLS(NN,1)=mean(DSconcat_OT_classDLS(NN,52:306),2);
+    end
+    TMPdls=MaxValDLS;
+    TMPdls(isnan(TMPdls))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
+    [~,SORTimgDLS]=sort(TMPdls);
+    DSconcat_OT_classDLS=DSconcat_OT_classDLS(SORTimgDLS,:);
+    
+    % sorting DMS
+    DSconcat_OT_classDMS=DSconcat_OT_Phasic(DMSselectionP,:);
+    for NN=1:size(DSconcat_OT_classDMS,1)
+        MaxValDMS(NN,1)=mean(DSconcat_OT_classDMS(NN,52:306),2);
+    end
+    TMPdms=MaxValDMS;
+    TMPdms(isnan(TMPdms))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
+    [~,SORTimgDMS]=sort(TMPdms);
+    DSconcat_OT_classDMS=DSconcat_OT_classDMS(SORTimgDMS,:);
+    
+    for j=1:2
+        if j==1
+            DSconcat_OT_class_DLSDMS=DSconcat_OT_classDLS;
+        else
+            DSconcat_OT_class_DLSDMS=DSconcat_OT_classDMS;
+        end
+    subplot(9,4,[3+(i-1)*12+(j-1)*4 4+(i-1)*12+(j-1)*4])
+    imagesc(timeconcat,[1 size(DSconcat_OT_class_DLSDMS,1)],DSconcat_OT_class_DLSDMS,Clim);
+    colormap(jet);
+    hold on
+    plot([51 51],[size(DSconcat_OT_class_DLSDMS,1)+0.5 0.5],'k')
+    plot([102 102],[size(DSconcat_OT_class_DLSDMS,1)+0.5 0.5],'k')
+    plot([153 153],[size(DSconcat_OT_class_DLSDMS,1)+0.5 0.5],'k')
+    plot([204 204],[size(DSconcat_OT_class_DLSDMS,1)+0.5 0.5],'k')
+    plot([255 255],[size(DSconcat_OT_class_DLSDMS,1)+0.5 0.5],'k')
+    plot([306 306],[size(DSconcat_OT_class_DLSDMS,1)+0.5 0.5],'k')
+    set(gca,'XTick',[0:25.5:357]);
+    set(gca,'xticklabel',Eventnames)
+    hold off
+    
+    % plot average PSTH
+    subplot(9,4,[11+(i-1)*12 12+(i-1)*12])
+    DLSpsth=nanmean(DSconcat_OT_classDLS(:,timeconcat),1);
+    DLSsem=nanste(DSconcat_OT_classDLS(:,timeconcat),1);
+    DLSup=DLSpsth+DLSsem;
+    DLSdown=DLSpsth-DLSsem;
+    
+    DMSpsth=nanmean(DSconcat_OT_classDMS(:,timeconcat),1);
+    DMSsem=nanste(DSconcat_OT_classDMS(:,timeconcat),1);
+    DMSup=DMSpsth+DMSsem;
+    DMSdown=DMSpsth-DMSsem;
+    hold on
+    patch([timeconcat,timeconcat(end:-1:1)],[DLSup,DLSdown(end:-1:1)],myblue,'EdgeColor','none');alpha(0.5);
+    g1=plot(timeconcat,DLSpsth,'Color',myblue,'linewidth',1.5);
+    patch([timeconcat,timeconcat(end:-1:1)],[DMSup,DMSdown(end:-1:1)],'r','EdgeColor','none');alpha(0.5);
+    g2=plot(timeconcat,DMSpsth,'r','linewidth',1.5);
+    
+    plot([0 357], [0 0],'LineStyle',':','Color','k');
+    plot([51 51],Yaxis,'LineStyle',':','Color','k')
+    plot([102 102],Yaxis,'LineStyle',':','Color','k')
+    plot([153 153],Yaxis,'LineStyle',':','Color','k')
+    plot([204 204],Yaxis,'LineStyle',':','Color','k')
+    plot([255 255],Yaxis,'LineStyle',':','Color','k')
+    plot([306 306],Yaxis,'LineStyle',':','Color','k')
+    axis([Xaxis2,Yaxis]);
+    set(gca,'XTick',[0:25.5:357]);
+    set(gca,'xticklabel',Eventnames)
+    hold off
+    end
+end
+
+
+%% proportion of neurons in each classes / plot stack bars / 2*3 classes with separation phasic / non phasic neurons
+
+DSselection=[];
+
+for j=1:2 % loop thru structure
+    for k=1:3 % loop thru class
+        DSselection(:,1)=PhasicNeurons_ID(:,2)==BrainRegion(j) & PhasicNeurons_ID(:,3)==k;
+        DSselectionNumber(j,k)=sum(DSselection);
+        DSselectionProp(j,k)=100*sum(DSselection)/sum(PhasicNeurons_ID(:,2)==BrainRegion(j));
+    end
+end
+
+subplot(9,4,[30 34])
+X=categorical({'DLS','DMS'});
+b1=bar(X,DSselectionProp,'stacked');
+legend('Start','Stop','Middle')
+
+
+%% Hierarchical clustering of Nonphasic neurons based on Seqeunce-related activity
+figure
+sel=Coord(:,4)~=0 & Celltype==1 & phasicID==1;
+X=binnedDS(sel,:);
+
+Z = linkage(X,'ward');
+c = cluster(Z,'Maxclust',4);
+eva = evalclusters(X,'linkage','CalinskiHarabasz','KList',[1:10]); %% You can change the criterion here. There are 3 other options. Gap, 'Silhouette' , 'DaviesBouldin'. You can do help evalclusters to use this code.
+D = pdist(X);
+leafOrder = optimalleaforder(Z,D);
+cutoff = median([Z(end-2,3) Z(end-1,3)]);
+subplot(6,4,[17 21])
+dendrogram(Z,'ColorThreshold',5, 'Reorder',leafOrder);
+
+ClassNONPHASIC(:,1)=eva.OptimalY;
+ClassNONPHASIC(:,2:4)=X;
+
+save('ClassHC.mat','ClassPHAS','ClassNONPHASIC')
+
+
+%% eva has the results of fitting the clusters. eva.OptimalK tells us how many clusters fit this data best based on the criterion used.
+%% eva.OptimalY tells us the cluster identities of each point.
+%% CH criterion clustering labels
+ind=[];
+DSconcat_OT_NonPhasic=DSconcat_OTz(sel,:);
+DLSsel=Coord(:,4)==10 & Celltype==1 & phasicID==1;
+DMSsel=Coord(:,4)==20 & Celltype==1 & phasicID==1;
+
+% sorting all neurons
+DSconcat_OT_NPclassDLS=DSconcat_OTz(DLSsel,:);
+
+
+for NN=1:size(DSconcat_OT_NPclassDLS,1)
+    MeanValdls(NN,1)=mean(DSconcat_OT_NPclassDLS(NN,52:306),2);
+end
+TMPdls=MeanValdls;
+TMPdls(isnan(TMPdls))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
+[~,SORTimgDLS]=sort(TMPdls);
+DSconcat_OT_NPclassDLS=DSconcat_OT_NPclassDLS(SORTimgDLS,:);
+
+DSconcat_OT_NPclassDMS=DSconcat_OTz(DMSsel,:);
+for NN=1:size(DSconcat_OT_NPclassDMS,1)
+    MeanValdms(NN,1)=mean(DSconcat_OT_NPclassDMS(NN,52:306),2);
+end
+TMPdms=MeanValdms;
+TMPdms(isnan(TMPdms))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
+[~,SORTimgDMS]=sort(TMPdms);
+DSconcat_OT_NPclassDMS=DSconcat_OT_NPclassDMS(SORTimgDMS,:);
+
+% plot heatmap
+subplot(6,4,[1 2 5 6])
+imagesc(timeconcat,[1 size(DSconcat_OT_NPclassDLS,1)],DSconcat_OT_NPclassDLS,Clim); %colorbar;axis tight;
+colormap(jet);
+hold on
+plot([51 51],[size(DSconcat_OT_NPclassDLS,1)+0.5 0.5],'k')
+plot([102 102],[size(DSconcat_OT_NPclassDLS,1)+0.5 0.5],'k')
+plot([153 153],[size(DSconcat_OT_NPclassDLS,1)+0.5 0.5],'k')
+plot([204 204],[size(DSconcat_OT_NPclassDLS,1)+0.5 0.5],'k')
+plot([255 255],[size(DSconcat_OT_NPclassDLS,1)+0.5 0.5],'k')
+plot([306 306],[size(DSconcat_OT_NPclassDLS,1)+0.5 0.5],'k')
+set(gca,'XTick',[0:25.5:357]);
+set(gca,'xticklabel',Eventnames)
+hold off
+
+subplot(6,4,[9 10 13 14])
+imagesc(timeconcat,[1 size(DSconcat_OT_NPclassDMS,1)],DSconcat_OT_NPclassDMS,Clim); %colorbar;axis tight;
+colormap(jet);
+hold on
+plot([51 51],[size(DSconcat_OT_NPclassDMS,1)+0.5 0.5],'k')
+plot([102 102],[size(DSconcat_OT_NPclassDMS,1)+0.5 0.5],'k')
+plot([153 153],[size(DSconcat_OT_NPclassDMS,1)+0.5 0.5],'k')
+plot([204 204],[size(DSconcat_OT_NPclassDMS,1)+0.5 0.5],'k')
+plot([255 255],[size(DSconcat_OT_NPclassDMS,1)+0.5 0.5],'k')
+plot([306 306],[size(DSconcat_OT_NPclassDMS,1)+0.5 0.5],'k')
+set(gca,'XTick',[0:25.5:357]);
+set(gca,'xticklabel',Eventnames)
+hold off
+
+for i=1:max(eva.OptimalY)
+    DSconcat_OT_classDLS=[]; SORTimgDLS=[];MaxValDLS=[];TMPdls=[];
+    DSconcat_OT_classDMS=[]; SORTimgDMS=[];MaxValDMS=[];TMPdms=[];
+    ind=find(eva.OptimalY==i);
+    
+    NPhasicNeurons_ID(:,1)=find(phasicID==1 & ~isnan(phasicID));
+    NPhasicNeurons_ID(:,2)=Coord(NPhasicNeurons_ID(:,1),4);
+    NPhasicNeurons_ID(:,3)=eva.OptimalY;
+    
+    DLSselectionNP=NPhasicNeurons_ID(:,2)==10 &  NPhasicNeurons_ID(:,3)==i;
+    DMSselectionNP=NPhasicNeurons_ID(:,2)==20 &  NPhasicNeurons_ID(:,3)==i;
+    
+    
+    % sorting DLS
+    DSconcat_OT_classDLS=DSconcat_OT_NonPhasic(DLSselectionNP,:);
+    for NN=1:size(DSconcat_OT_classDLS,1)
+        MaxValDLS(NN,1)=mean(DSconcat_OT_classDLS(NN,52:306),2);
+    end
+    TMPdls=MaxValDLS;
+    TMPdls(isnan(TMPdls))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
+    [~,SORTimgDLS]=sort(TMPdls);
+    DSconcat_OT_classDLS=DSconcat_OT_classDLS(SORTimgDLS,:);
+    
+    % sorting DMS
+    DSconcat_OT_classDMS=DSconcat_OT_NonPhasic(DMSselectionNP,:);
+    for NN=1:size(DSconcat_OT_classDMS,1)
+        MaxValDMS(NN,1)=mean(DSconcat_OT_classDMS(NN,52:306),2);
+    end
+    TMPdms=MaxValDMS;
+    TMPdms(isnan(TMPdms))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
+    [~,SORTimgDMS]=sort(TMPdms);
+    DSconcat_OT_classDMS=DSconcat_OT_classDMS(SORTimgDMS,:);
+    
+    for j=1:2
+        if j==1
+            DSconcat_OT_class_DLSDMS= DSconcat_OT_classDLS;
+        else
+            DSconcat_OT_class_DLSDMS= DSconcat_OT_classDMS;
+        end
+        
+        subplot(6,4,[3+(i-1)*12+(j-1)*4 4+(i-1)*12+(j-1)*4])
+        imagesc(timeconcat,[1 size(DSconcat_OT_class_DLSDMS,1)],DSconcat_OT_class_DLSDMS,Clim);
+        colormap(jet);
+        hold on
+        plot([51 51],[size(DSconcat_OT_class_DLSDMS,1)+0.5 0.5],'k')
+        plot([102 102],[size(DSconcat_OT_class_DLSDMS,1)+0.5 0.5],'k')
+        plot([153 153],[size(DSconcat_OT_class_DLSDMS,1)+0.5 0.5],'k')
+        plot([204 204],[size(DSconcat_OT_class_DLSDMS,1)+0.5 0.5],'k')
+        plot([255 255],[size(DSconcat_OT_class_DLSDMS,1)+0.5 0.5],'k')
+        plot([306 306],[size(DSconcat_OT_class_DLSDMS,1)+0.5 0.5],'k')
+        set(gca,'XTick',[0:25.5:357]);
+        set(gca,'xticklabel',Eventnames)
+        hold off
+        
+        if i==1
+            Yaxis=[-3 2];
+        else
+            Yaxis=[-1 15];
+        end
+        % plot average PSTH
+        subplot(6,4,[11+(i-1)*12 12+(i-1)*12])
+        DLSpsth=nanmean(DSconcat_OT_classDLS(:,timeconcat),1);
+        DLSsem=nanste(DSconcat_OT_classDLS(:,timeconcat),1);
+        DLSup=DLSpsth+DLSsem;
+        DLSdown=DLSpsth-DLSsem;
+        
+        DMSpsth=nanmean(DSconcat_OT_classDMS(:,timeconcat),1);
+        DMSsem=nanste(DSconcat_OT_classDMS(:,timeconcat),1);
+        DMSup=DMSpsth+DMSsem;
+        DMSdown=DMSpsth-DMSsem;
+        hold on
+        patch([timeconcat,timeconcat(end:-1:1)],[DLSup,DLSdown(end:-1:1)],myblue,'EdgeColor','none');alpha(0.5);
+        g1=plot(timeconcat,DLSpsth,'Color',myblue,'linewidth',1.5);
+        patch([timeconcat,timeconcat(end:-1:1)],[DMSup,DMSdown(end:-1:1)],'r','EdgeColor','none');alpha(0.5);
+        g2=plot(timeconcat,DMSpsth,'r','linewidth',1.5);
+        
+        plot([0 357], [0 0],'LineStyle',':','Color','k');
+        plot([51 51],Yaxis,'LineStyle',':','Color','k')
+        plot([102 102],Yaxis,'LineStyle',':','Color','k')
+        plot([153 153],Yaxis,'LineStyle',':','Color','k')
+        plot([204 204],Yaxis,'LineStyle',':','Color','k')
+        plot([255 255],Yaxis,'LineStyle',':','Color','k')
+        plot([306 306],Yaxis,'LineStyle',':','Color','k')
+        axis([Xaxis2,Yaxis]);
+        set(gca,'XTick',[0:25.5:357]);
+        set(gca,'xticklabel',Eventnames)
+        hold off
+    end
+end
+
+
+%% proportion of neurons in each classes / plot stack bars / 2*3 classes with separation phasic / non phasic neurons
+
+DSselection=[];
+
+for j=1:2 % loop thru structure
+    for k=1:6 % loop thru class
+        DSselection(:,1)=NPhasicNeurons_ID(:,2)==BrainRegion(j) & NPhasicNeurons_ID(:,3)==k;
+        DSselectionNumberNP(j,k)=sum(DSselection);
+        DSselectionPropNP(j,k)=100*sum(DSselection)/sum(NPhasicNeurons_ID(:,2)==BrainRegion(j));
+    end
+end
+
+subplot(6,4,[18 22])
+X=categorical({'DLS','DMS'});
+b1=bar(X,DSselectionPropNP,'stacked');
+legend('INH','EXC');
+
+
+
+%% make table for statistics
+load('RatID_extendedtraining.mat');
+ Class_OT(1:length(Coord(:,4)),1)=phasicID;
+ Class_OT(PhasicNeurons_ID(:,1),2)=PhasicNeurons_ID(:,3);
+ Class_OT(NPhasicNeurons_ID(:,1),2)=NPhasicNeurons_ID(:,3)+3;
+ DSmeanz_OT=[];
+
+ DSmeanz_OT = cat(2,DSmeanz_OT,Ev(7).Meanz(:,1)); % use meanZ instead of average of bin
+  for k=1:5
+      DSmeanz_OT = cat(2,DSmeanz_OT,Ev(k).MeanzPRE(:,1));
+      DSmeanz_OT = cat(2,DSmeanz_OT,Ev(k).Meanz(:,1));
+  end
+ DSmeanz_OT = cat(2,DSmeanz_OT,Ev(6).MeanzPRE(:,1));
+
+
+ table_OT=cat(2,RatID,Coord(:,4),Celltype(:,1),Class_OT,DSmeanz_OT, TRN);
+ save('table_OT.mat','table_OT')
+
+  nbevent=[1:12];
+ 
+ selectionSTAT=table_OT(:,3)~=0 & table_OT(:,4)==1 & TRN(:,1)~=0;
+ tableSTAT=array2table(table_OT(selectionSTAT,1:18),'VariableNames',{'rat','habit','region','celltype','phasicID','class','win1','win2','win3','win4','win5','win6','win7','win8','win9','win10','win11','win12'});
+
+ rm = fitrm(tableSTAT,'win1-win12~region','WithinDesign',nbevent');
+ ranovatbl_event = ranova(rm);   
+ Btw_ranovatbl_event = anova(rm);
+ 
+ %control for effect of habit
+ rm = fitrm(tableSTAT,'win1-win12~region*habit','WithinDesign',nbevent');
+ ranovatbl_event_habit = ranova(rm);   
+ Btw_ranovatbl_event_habit = anova(rm);
+ 
+ for i=1:max(table_OT(:,6))
+    selectionSTAT=table_OT(:,3)~=0 & table_OT(:,4)==1 & table_OT(:,6)==i; 
+    tableSTAT=array2table(table_OT(selectionSTAT,1:18),'VariableNames',{'rat','habit','region','celltype','phasicID','class','win1','win2','win3','win4','win5','win6','win7','win8','win9','win10','win11','win12'});
+
+    rm = fitrm(tableSTAT,'win1-win12~region','WithinDesign',nbevent');
+    stat(i).ranovatbl_event = ranova(rm);
+    stat(i).Btw_ranovatbl_event = anova(rm);  
+    
+    rm = fitrm(tableSTAT,'win1-win12~region*habit','WithinDesign',nbevent');
+    stat(i).ranovatbl_event_Habit = ranova(rm);
+    stat(i).Btw_ranovatbl_event_Habit = anova(rm);   
+ end
+
+
+%% permutation test
+allsamples=X_ratio;
+
+for p=1:10000
+    shuffsamples=c_ratio(randperm(length(c_ratio)));
+    MeanCluster1(p,:)=mean(allsamples(shuffsamples==1,:),1);
+    MeanCluster2(p,:)=mean(allsamples(shuffsamples==2,:),1);
+    DistCluster(p,1) = pdist([MeanCluster1(p,:);MeanCluster2(p,:)],'euclidean');
+end
+MeanP=mean(allsamples(c_ratio==2,:),1);
+MeanNP=mean(allsamples(c_ratio==1,:),1);
+RealDist1=pdist([MeanP;MeanNP],'euclidean');
+Classpval=(sum(RealDist1>DistCluster)+1)/(length(DistCluster)+1);
+
+figure;
+hist(DistCluster);
+hold on;
+plot([RealDist1 RealDist1],[0 4000]);
+xlabel('Between-cluster distance')
+ylabel('# iterations')
+text(0.225,100,'Phasic vs Non-phasic classification','color','b','rotation',90);
+
+%% Clustering on phasic and non phasic neurons 
+
+rangeHighFreq=[2:0.1:8];
+rangeLowFreq=[0.3:0.1:1.1];
+nbClassif2=0;
+for i=1:length(rangeHighFreq)
+    for j=1:length(rangeLowFreq)
+        
+        f_low=rangeLowFreq(j);
+        f_high=rangeHighFreq(i);
+        
+        %fourier analysis
+        [X_all,f_all,P_all,Nt_all]=myfft(DSconcat_OTshort',0.01,0);
+        total_power=sum(P_all);
+        P_all_norm = bsxfun(@rdivide,P_all, total_power);
+        
+        a=find(f_all>f_low);b=find(f_all<f_high);
+        lower_window_sum_all=sum(P_all_norm(1:a(1),:));
+        upper_window_sum_all=sum(P_all_norm(a(1):b(end),:));
+        
+        ratio_all=lower_window_sum_all./upper_window_sum_all;
+        
+        %hierarchical clustering
+        X_ratio=[[lower_window_sum_all]' [upper_window_sum_all]']; %#ok<NBRAK>
+        eva_Ch_ratio = evalclusters(X_ratio,'linkage','CalinskiHarabasz','KList',[1:10]);
+        OptimalK(i,j)=eva_Ch_ratio.OptimalK;
+        c_ratio3=eva_Ch_ratio.OptimalY;
+        if eva_Ch_ratio.OptimalK==2
+            MeanP=mean(X_ratio(c_ratio3==2,:),1);
+            MeanNP=mean(X_ratio(c_ratio3==1,:),1);
+            RealDist(i,j)=pdist([MeanP;MeanNP],'euclidean');
+            nbClassif2=nbClassif2+1;
+            OptimalY(:,nbClassif2)=eva_Ch_ratio.OptimalY;
+        end
+    end
+end
+
+figure
+[X,Y] = meshgrid(rangeLowFreq,rangeHighFreq);
+pcolor(X,Y,OptimalK)
+colorbar
+xlabel('Lower frequency limit (Hz)')
+ylabel('Higher frequency limit (Hz)')
+
+corrOptimalK=OptimalK(:,1:8);
+percentK2=100*sum(sum(corrOptimalK==2))/(size(corrOptimalK,1)*size(corrOptimalK,2));
+percentK3=100*sum(sum(corrOptimalK==3))/(size(corrOptimalK,1)*size(corrOptimalK,2));
+
+distancesDist=reshape(RealDist,size(RealDist,1)*size(RealDist,2),1);
+figure;
+hist(distancesDist(distancesDist(:,1)~=0,1));
+hold on
+hist(DistCluster);
+plot([RealDist1 RealDist1],[0 4000]);
+xlabel('Between-cluster distance')
+ylabel('# iterations')
+text(0.225,100,'Phasic vs Non-phasic classification','color','b','rotation',90);
+
+A(1,1)=length(find(OptimalK==2));
+A(1,2)=length(find(OptimalK==3));
+A(1,3)=size(OptimalK,1)*size(OptimalK,2);
+A(2,1)=100*A(1,1)./A(1,3);
+A(2,2)=100*A(1,2)./A(1,3);
+A(2,3)=100*(A(1,1)+A(1,2))./A(1,3);
+
+%% look at 2 classes when OptimalK=2
+
+for j=1:size(OptimalY,2)
+    for i=1:2
+        NbNeuronClass(j,i)=sum(OptimalY(:,j)==i);
+    end
+    NbNeuronClass(j,:)=sort(NbNeuronClass(j,:),2);
+end
+
+figure
+histogram(NbNeuronClass(:,1))
+hold on
+histogram(NbNeuronClass(:,2))
+hold off
+xlabel('number of neurons in each class')
+ylabel('# iterations')
+legend('Non-Phasic','Phasic')

scripts/G_generateFigure4.m

@@ -0,0 +1,210 @@
+%% generate figure 4
+% Analysis of phasic Start and Stop neurons in the extended training dataset
+
+clc;clear;close all
+Task=1;
+load('table_OT.mat')
+load('Rextendedtraining.mat');
+load('RatID_extendedtraining.mat');
+load('Celltype_extendedTraining.mat');
+load('TRN_DT5_extendedtrain.mat')
+%% --- --- Plotting color-coded and average PSTHs --- ---
+
+Yaxis=[-2 6];
+Xaxis1=[-0.25 1];
+Xaxis2=[-2 0.25];
+
+Ishow1=find(R.Param.Tm>=Xaxis1(1) & R.Param.Tm<=Xaxis1(2));
+Ishow2=find(R.Param.Tm>=Xaxis2(1) & R.Param.Tm<=Xaxis2(2));
+
+Eventnames={'';'LI';'';'LP1';'';'LP2';'';'LP3';'';'LP4';'';'LP5';'';'PE'};
+namebars=categorical({'All START','EXC pre-LP1'});
+labelxaxis={'time from LI','time from LP1'};
+
+% ----------- COLORS ---------
+myblue=[0 113/255 188/255];
+mypurple=[200/255 0 255/255];
+mydarkblue=[0 0/255 255/255];
+TickSize=[0.015 0.02];
+
+time1=-0.25:0.01:1;
+time2=-2:0.01:0.25;
+
+%% Pie Plot
+%proportion of start neurons responding before the 1st lever press
+
+DLSselectionPRELP_start=R.Coord(:,4)==10 & Celltype(:,1)==1 & TRNmatrice(:,1)==1 & table_OT(:,5)==2 & table_OT(:,6)==1;
+DLSselectionSTART=R.Coord(:,4)==10 & Celltype(:,1)==1 & table_OT(:,5)==2 & table_OT(:,6)==1;
+
+DMSselectionPRELP_start=R.Coord(:,4)==20 & Celltype(:,1)==1 & TRNmatrice(:,1)==1 & table_OT(:,5)==2 & table_OT(:,6)==1;
+DMSselectionSTART=R.Coord(:,4)==20 & Celltype(:,1)==1 & table_OT(:,5)==2 & table_OT(:,6)==1;
+
+DSselectionPRELP_start=R.Coord(:,4)~=0 & Celltype(:,1)==1 & TRNmatrice(:,1)==1 & table_OT(:,5)==2 & table_OT(:,6)==1;
+
+DLS_X(1)=sum(DLSselectionSTART);
+DMS_X(1)=sum(DMSselectionSTART);
+DLS_X(2)=sum(DLSselectionPRELP_start);
+DMS_X(2)=sum(DMSselectionPRELP_start);
+DS_X=sum(DSselectionPRELP_start);
+
+DLS_tot(1)=sum(R.Coord(:,4)==10 & Celltype(:,1)==1 & table_OT(:,6)==1);
+DMS_tot(1)=sum(R.Coord(:,4)==20 & Celltype(:,1)==1 & table_OT(:,6)==1);
+DS_tot=sum(R.Coord(:,4)~=0 & Celltype(:,1)==1 & table_OT(:,6)==1);
+
+DS_Percent(1,:)=DS_X*100./DS_tot;
+DS_Percent(2,:)=100-DS_Percent(1,:);
+
+figure
+subplot(2,4,1)
+pie(DS_Percent)
+title('Phasic Start');
+legend({'EXC pre-LP1','All Start'})
+
+
+%% Average PSTH start neurons
+
+nplot=0;
+for i=1:2 %loop thru event
+    nplot=nplot+1;
+    
+    selectionDLS=R.Coord(:,4)==10 & Celltype(:,1)==1 & table_OT(:,6)==1;
+    selectionDMS=R.Coord(:,4)==20 & Celltype(:,1)==1 & table_OT(:,6)==1;
+    
+    if i==1
+        DSconcat_OTz = R.Ev(7).PSTHz(:,Ishow1);
+        Xaxis=Xaxis1;time=time1;
+        tickdist=0.25;
+    else
+        DSconcat_OTz = R.Ev(1).PSTHz(:,Ishow2);
+        Xaxis=Xaxis2;time=time2;
+        tickdist=0.5;
+    end
+    
+    subplot(2,4,[3+(nplot-1)])
+    
+    DLSconcat_OT=DSconcat_OTz(selectionDLS,:);
+    DMSconcat_OT=DSconcat_OTz(selectionDMS,:);
+    
+    DLSpsth=nanmean(DLSconcat_OT,1);
+    DLSsem=nanste(DLSconcat_OT,1);
+    DLSup=DLSpsth+DLSsem;
+    DLSdown=DLSpsth-DLSsem;
+    hold on;
+    patch([time,time(end:-1:1)],[DLSup,DLSdown(end:-1:1)],myblue,'EdgeColor','none');alpha(0.5);
+    plot(time,DLSpsth,'Color',myblue,'linewidth',1.5);
+    
+    DMSpsth=nanmean(DMSconcat_OT,1);
+    DMSsem=nanste(DMSconcat_OT,1);
+    DMSup=DMSpsth+DMSsem;
+    DMSdown=DMSpsth-DMSsem;
+    patch([time,time(end:-1:1)],[DMSup,DMSdown(end:-1:1)],'r','EdgeColor','none');alpha(0.5);
+    plot(time,DMSpsth,'r','linewidth',1.5);
+    
+    plot([Xaxis(1) Xaxis(end)], [0 0],'LineStyle',':','Color','k');
+    plot([0 0], [Yaxis(1) Yaxis(end)],'LineStyle',':','Color','k');
+    axis([Xaxis,Yaxis]);
+    set(gca,'XTick',[Xaxis(1):tickdist:Xaxis(end)]);
+    title('START neurons');
+    ylabel('z-score');
+    xlabel(labelxaxis(i));
+end
+ 
+%% Average PSTH stop neurons
+Xaxis2=[-1 1];
+Ishow2=find(R.Param.Tm>=Xaxis2(1) & R.Param.Tm<=Xaxis2(2));
+time1=-0.25:0.01:1;
+time2=-1:0.01:1;
+labelxaxis={'time from LLP','time from PE'};
+
+nplot=0;
+for i=1:2 %loop thru event
+    nplot=nplot+1;
+    
+    selectionDLS=R.Coord(:,4)==10 & Celltype(:,1)==1 & table_OT(:,6)==2;
+    selectionDMS=R.Coord(:,4)==20 & Celltype(:,1)==1 & table_OT(:,6)==2;
+    
+    if i==1
+        DSconcat_OTz = R.Ev(5).PSTHz(:,Ishow1);
+        Xaxis=Xaxis1;time=time1;
+        tickdist=0.25;
+    else
+        DSconcat_OTz = R.Ev(6).PSTHz(:,Ishow2);
+        Xaxis=Xaxis2;time=time2;
+        tickdist=0.5;
+    end
+    
+    subplot(2,4,[5+(nplot-1)])
+    
+    DLSconcat_OT=DSconcat_OTz(selectionDLS,:);
+    DMSconcat_OT=DSconcat_OTz(selectionDMS,:);
+    
+    DLSpsth=nanmean(DLSconcat_OT,1);
+    DLSsem=nanste(DLSconcat_OT,1);
+    DLSup=DLSpsth+DLSsem;
+    DLSdown=DLSpsth-DLSsem;
+    hold on;
+    patch([time,time(end:-1:1)],[DLSup,DLSdown(end:-1:1)],myblue,'EdgeColor','none');alpha(0.5);
+    plot(time,DLSpsth,'Color',myblue,'linewidth',1.5);
+    
+    DMSpsth=nanmean(DMSconcat_OT,1);
+    DMSsem=nanste(DMSconcat_OT,1);
+    DMSup=DMSpsth+DMSsem;
+    DMSdown=DMSpsth-DMSsem;
+    patch([time,time(end:-1:1)],[DMSup,DMSdown(end:-1:1)],'r','EdgeColor','none');alpha(0.5);
+    plot(time,DMSpsth,'r','linewidth',1.5);
+    
+    %plot([0 0], [-2 2.5],'LineStyle',':','Color','k');
+    plot([Xaxis(1) Xaxis(end)], [0 0],'LineStyle',':','Color','k');
+    plot([0 0], [Yaxis(1) Yaxis(end)],'LineStyle',':','Color','k');
+    axis([Xaxis,Yaxis]);
+    set(gca,'XTick',[Xaxis(1):tickdist:Xaxis(end)]);
+    title('STOP neurons');
+    ylabel('z-score');
+    xlabel(labelxaxis(i));
+end
+
+ %% Stop neurons, activity around PE at the end of trials VS during inter-trial intervals
+
+titlelabel={'STOP neurons DLS' 'STOP neurons DMS'};
+selectionDLS=R.Coord(:,4)==10 & Celltype(:,1)==1 & table_OT(:,6)==2;
+selectionDMS=R.Coord(:,4)==20 & Celltype(:,1)==1 & table_OT(:,6)==2;
+
+Xaxis1=[-1 1]; Yaxis=[-2 6];
+Ishow=find(R.Param.Tm>=Xaxis1(1) & R.Param.Tm<=Xaxis1(2));
+
+time=[-1:0.01:1];
+colors{1,1}=[0.00  0.75  0.75];
+colors{1,2}=[0.3  0.1  0.8];
+
+for i=1:2
+    if i==1
+        selection=selectionDLS;
+    else
+        selection=selectionDMS;
+    end
+    meanPostRew_PE=nanmean(R.Ev(6).PSTHz(selection,Ishow),1);
+    stePostRew_PE=nanste(R.Ev(6).PSTHz(selection,Ishow),1);
+    meanITI_PE=nanmean(R.Ev(12).PSTHz(selection,Ishow),1);
+    steITI_PE=nanste(R.Ev(12).PSTHz(selection,Ishow),1);
+    
+    PRPE_up=meanPostRew_PE+stePostRew_PE;
+    PRPE_down=meanPostRew_PE-stePostRew_PE;
+    ITIPE_up=meanITI_PE+steITI_PE;
+    ITIPE_down=meanITI_PE-steITI_PE;
+    
+    subplot(2,4,7+i-1)
+    patch([time,time(end:-1:1)],[PRPE_up,PRPE_down(end:-1:1)],colors{1,1},'EdgeColor','none');alpha(0.5);
+    hold on
+    g1=plot(time,meanPostRew_PE,'Color',colors{1,1},'linewidth',1.5);
+    patch([time,time(end:-1:1)],[ITIPE_up,ITIPE_down(end:-1:1)],colors{1,2},'EdgeColor','none');alpha(0.5);
+    g2=plot(time,meanITI_PE,'Color',colors{1,2},'linewidth',1.5);
+    
+    plot([-1 1], [0 0],'LineStyle',':','Color','k');
+    plot([0 0],Yaxis,'LineStyle',':','Color','k')
+    axis([Xaxis1,Yaxis]);
+    title(titlelabel(i));
+    set(gca,'XTick',[-1:0.5:1]);
+    ylabel('z-score');
+    xlabel('time from PE');
+    hold off
+end

scripts/H_generateFigure6.m

@@ -0,0 +1,527 @@
+%% generate figure 6 and figure 6 - supplement figure 1
+% fourrier analysis - hierarchical clustering - Random forest approach on
+% early training dataset
+% statistics for the early training dataset
+
+
+clc;clear;close all;
+load ('Rearlytraining_light.mat');
+load ('Celltype_earlyTraining.mat');
+load('TRN_DT5_earlytrain.mat')
+
+Xevent=[-0.25 0.25];
+Yaxis=[-5 10];Xaxis=[1 357];
+Ishow=find(Tm>=Xevent(1) & Tm<=Xevent(2));
+time=[1:1:357];
+BrainRegion=[10 20];
+% ----------- COLORS ---------
+myblue=[0 113/255 188/255];
+mypurple=[200/255 0 255/255];
+mydarkblue=[0 0/255 255/255];
+TickSize=[0.015 0.02];
+Clim=[-5 5];
+Eventnames={'';'LI';'';'LP1';'';'LP2';'';'LP3';'';'LP4';'';'LP5';'';'PE'};
+
+DSconcat_ACQ_allsession=[];Coord_allsession=[];Ninfo_allsession=[];DSmeanz_ACQ_allsession=[];
+TRN_allsession=[];EXCINH_allsession=[];sessionID_allsession=[];
+for i=1:length(Ses)
+    
+    %normalization based on max
+    for k=1:length(Erefnames)
+        Ev(k).PSTH_norm1=normalize(Ses(i).Ev(k).PSTHraw,Ses(i).Ev(k).PSTHrawBL,1);
+    end
+    
+    DSconcat_ACQ=[];DSmeanz_ACQ=[];
+    %concat with normalized data
+    DSconcat_ACQ = Ev(7).PSTH_norm1(:,Ishow);
+    DSmeanz_ACQ = Ses(i).Ev(7).Meanz(:,1);
+    for k=1:5
+        DSconcat_ACQ = cat(2,DSconcat_ACQ,Ev(k).PSTH_norm1(:,Ishow));
+        DSmeanz_ACQ = cat(2,DSmeanz_ACQ,Ses(i).Ev(k).MeanzPRE(:,1));
+        DSmeanz_ACQ = cat(2,DSmeanz_ACQ,Ses(i).Ev(k).Meanz(:,1));
+    end
+    DSconcat_ACQ = cat(2,DSconcat_ACQ,Ev(6).PSTH_norm1(:,Ishow));
+    DSmeanz_ACQ = cat(2,DSmeanz_ACQ,Ses(i).Ev(6).MeanzPRE(:,1));
+    
+    DSconcat_ACQ_allsession=cat(1,DSconcat_ACQ_allsession,DSconcat_ACQ);
+    Coord_allsession=cat(1,Coord_allsession,Ses(i).Coord(:,4));
+    Ninfo_allsession=cat(1,Ninfo_allsession,Ses(i).Ninfo);
+    DSmeanz_ACQ_allsession=cat(1,DSmeanz_ACQ_allsession,DSmeanz_ACQ);
+    TRN_allsession=cat(1,TRN_allsession,sesTRN(i).TRN(:,1));
+    sessionID_allsession=cat(1,sessionID_allsession,repmat(i,length(Ses(i).TRN(:,1)),1));
+    if i>3
+        EXCINH_allsession=cat(1,EXCINH_allsession,sesTRN(i).TRN(:,2));
+    end
+end
+
+s2={'YVG04','YVG05','YVG08','YVJ11','YVJ15'};
+for j=1:5
+    for i=1:length(Ninfo_allsession)
+        RatID(i,j)=strncmp(Ninfo_allsession{i,1},s2(j),5);
+    end
+end
+RatID_ACQ=sum(RatID,2);
+% concat with Zscore
+DSconcatZ_ACQ_allsession=[];DSevent_ACQ_allsession=[];
+for i=1:length(Ses)
+    DSconcat_ACQ=[];DSevent_ACQ=[];
+    
+    DSconcat_ACQ=Ses(i).Ev(7).PSTHz(:,Ishow);
+    DSevent_ACQ=Ses(i).Ev(7).Meanz(:,1);
+    for k=1:6 %loop through event
+        DSconcat_ACQ=cat(2,DSconcat_ACQ,Ses(i).Ev(k).PSTHz(:,Ishow));
+        DSevent_ACQ = cat(2,DSevent_ACQ,Ses(i).Ev(k).Meanz(:,1));
+    end
+    DSconcatZ_ACQ_allsession=cat(1,DSconcatZ_ACQ_allsession,DSconcat_ACQ);
+    DSevent_ACQ_allsession=cat(1,DSevent_ACQ_allsession,DSevent_ACQ);
+end
+
+selection=Coord_allsession~=0 & Celltype(:,2)==1 & ~isnan(mean(DSconcat_ACQ_allsession,2));
+DSconcat_ACQshort=DSconcat_ACQ_allsession(selection,:);
+
+DSevent_ACQshort=DSevent_ACQ_allsession(selection,:);
+Coord_short=Coord_allsession(selection,:);
+Ninfo_short=Ninfo_allsession(selection,:);
+
+%% Clustering on phasic and non phasic neurons with Fourier analysis 
+
+close all;
+f_low=1;f_high=4;
+
+[X_all,f_all,P_all,Nt_all]=myfft(DSconcat_ACQshort',0.01,1);
+total_power=sum(P_all);
+P_all_norm = bsxfun(@rdivide,P_all, total_power);
+%figure;plot(f_all,P_all_norm);set(gca,'ylim',[0 0.6]);xlabel('Hz');ylabel('Normalized power');
+%title('Normalized FFT of all neurons');
+
+a=find(f_all>f_low);b=find(f_all<f_high);
+interval_to_sum=[a(1) b(end)];
+lower_window_sum_all=sum(P_all_norm(1:a(1),:));
+upper_window_sum_all=sum(P_all_norm(a(1):b(end),:));
+
+ratio_all=lower_window_sum_all./upper_window_sum_all;
+close all
+
+X_ratio=[[lower_window_sum_all]' [upper_window_sum_all]']; %#ok<NBRAK>
+eva_Ch_ratio = evalclusters(X_ratio,'linkage','CalinskiHarabasz','KList',[1:10]);
+phasicID=eva_Ch_ratio.OptimalY;
+
+Z = linkage(X_ratio,'ward'); % X is the data set, neuron in rows, features in columns
+c = cluster(Z,'Maxclust',3); % determine the number cluster here)
+D = pdist(X_ratio);
+leafOrder = optimalleaforder(Z,D);
+cutoff = median([Z(end-2,3) Z(end-1,3)]);
+% 
+% figure
+% subplot(9,4,[29 33])
+% dendrogram(Z,'ColorThreshold',cutoff, 'Reorder',leafOrder);
+%% Figure Fourier transform
+
+c_ratio=eva_Ch_ratio.OptimalY;
+
+for i=1:length(c_ratio)
+    if c_ratio(i,1)==1
+        color(i,:)=[255/255 0/255 0/255]; %phasic red
+    elseif c_ratio(i,1)==2
+        color(i,:)=[0 255/255 0/255]; %sustained green
+    end
+end
+
+figure
+subplot(5,3,1)
+scatter(lower_window_sum_all,upper_window_sum_all,[],color)
+hold on
+title(strcat('hierarchical clustering'))
+xlabel('power in low freq') % x-axis label
+ylabel('power in int freq')% y-axis label
+legend('phasic','non-phasic')
+hold off
+
+%% visualization classes
+
+for i=1:2 
+    ind=find(phasicID==i);
+    t=num2str(i);
+    Clim=[-5 5];
+    DSconcat_ACQ_class=[];SORTimg=[];MaxValDLS=[];MeanValDLS=[]; MaxTimeDLS=[];TMP=[];NNid=[];
+    DSconcat_ACQ_class=DSconcat_ACQshort(ind,:);
+    
+    for NN=1:size(DSconcat_ACQ_class,1)
+        MeanValDLS(NN,1)=mean(DSconcat_ACQ_class(NN,52:306),2);
+        [MaxValDLS(NN,1),MaxInd]=max(DSconcat_ACQ_class(NN,:));
+        MaxTimeDLS(NN,1)=time(time(1)+MaxInd-1);
+        if i==2
+            MaxTimePhasic=MaxTimeDLS;
+        end
+    end
+    if i==1
+        TMP=MeanValDLS;
+        subplot(5,3,[4 7 10 13])
+    elseif i==2 
+        TMP=MeanValDLS;
+        subplot(5,3,[5 8 11 14])
+    end
+    TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
+    [~,SORTimg]=sort(TMP);
+    %NNid=cell2mat(Ninfo(phasic_neuron_selection,3));NNid2=NNid(SORTimg);
+    DSconcat_ACQ_class=DSconcat_ACQ_class(SORTimg,:);
+    
+    
+    imagesc(time,[1 size(DSconcat_ACQ_class,1)],DSconcat_ACQ_class);%colorbar;axis tight;
+    colormap(jet);
+    hold on
+    plot([51 51],[size(DSconcat_ACQ_class,1)+0.5 0.5],'k')
+    plot([102 102],[size(DSconcat_ACQ_class,1)+0.5 0.5],'k')
+    plot([153 153],[size(DSconcat_ACQ_class,1)+0.5 0.5],'k')
+    plot([204 204],[size(DSconcat_ACQ_class,1)+0.5 0.5],'k')
+    plot([255 255],[size(DSconcat_ACQ_class,1)+0.5 0.5],'k')
+    plot([306 306],[size(DSconcat_ACQ_class,1)+0.5 0.5],'k')
+    set(gca,'XTick',[0:25.5:357]);
+    set(gca,'xticklabel',Eventnames)
+    hold off
+    if i==1
+        ylabel('Non-Phasic neurons')
+    else
+        ylabel('Phasic neurons')
+    end
+end
+
+ClassACQ(1:length(Celltype),1)=0;
+ClassACQ(selection,1)=phasicID;
+
+
+%% Run random forest for phasic neurons
+
+load('ClassHC.mat' );
+XTrainingDATA=ClassPHAS(:,2:4);
+YTrainingDATA=ClassPHAS(:,1);
+
+LIwin=[26:51];%post lever insertion window
+LPwin=[52:280];%pre-1stLP to pre-lastLP window
+PEwin=[307:332];% pre-PE window
+
+binLI=mean(DSconcat_ACQ_allsession(:,LIwin),2);
+binLP=mean(DSconcat_ACQ_allsession(:,LPwin),2);
+binPE=mean(DSconcat_ACQ_allsession(:,PEwin),2);
+binnedDS=cat(2,binLI,binLP,binPE);   
+XTestingDATA=binnedDS(ClassACQ(:,1)==2,:);
+
+model = classRF_train(XTrainingDATA, YTrainingDATA); % train Random Forest model
+
+% prepare the data for FOR loop
+X = {XTrainingDATA, XTestingDATA};
+Y = {YTrainingDATA};
+classifier = {'Raw Data','Random Forest'};
+mc = ['r','g']; msym = ['o','o']; msiz = [3,15];
+
+
+for i = 1 : 2 % iterate through training and test data
+    x = X{i}; y = Y;
+    for j = 1 : 2 % iterate through different classifiers
+        
+        switch(j) % y_hat = running models ? model predictions : y
+            case 1 % raw data
+                y_hat = y;
+            case 2 % Random Forest
+                y_hat = classRF_predict(x ,model); % Random Forest predictions
+                if i==2
+                    Class_Acq(:,1)=y_hat;
+                end
+        end
+    end
+end
+
+ClassACQ(1:length(Celltype),2)=0;
+ClassACQ(ClassACQ(:,1)==2,2)=Class_Acq;
+
+%% Random forest for non phasic neurons
+
+XTrainingDATA_NP=ClassNONPHASIC(:,2:4);
+YTrainingDATA_NP=ClassNONPHASIC(:,1);
+XTestingDATA_NP=binnedDS(ClassACQ(:,1)==1,:);
+
+model = classRF_train(XTrainingDATA_NP, YTrainingDATA_NP); % train Random Forest model
+
+% prepare the data for FOR loop
+X = {XTrainingDATA_NP, XTestingDATA_NP};
+Y = {YTrainingDATA_NP};
+classifier = {'Raw Data','Random Forest'};
+mc = ['r','g']; msym = ['o','o']; msiz = [3,15];
+
+
+for i = 1 : 2 % iterate through training and test data
+    x = X{i}; y = Y;
+    for j = 1 : 2 % iterate through different classifiers
+        
+        switch(j) % y_hat = running models ? model predictions : y
+            case 1 % raw data
+                y_hat = y;
+            case 2 % Random Forest
+                y_hat = classRF_predict(x ,model); % Random Forest predictions
+                if i==2
+                    Class_Acq_NP(:,1)=y_hat;
+                end
+        end
+    end
+end
+
+ClassACQ(ClassACQ(:,1)==1,2)=Class_Acq_NP+3;
+
+%% selection of neurons per session per class
+% 1-3 is for phasic and 4-5 for non phasic
+
+for i=1:max(Celltype(:,1))
+    InDSsel=Ses(i).Coord(:,4)~=0;
+    Ses(i).Class(InDSsel,1)=ClassACQ(Celltype(:,1)==i,2);
+    Ses(i).Celltype(InDSsel,1)=Celltype(Celltype(:,1)==i,2);
+end
+
+for k=1:max(ClassACQ(:,2))
+    DLSconcat_ACQ_3ds(k).DLS=[];
+    DMSconcat_ACQ_3ds(k).DMS=[];
+    for i=4:13
+        DLSselection=Ses(i).Coord(:,4)==10 & Ses(i).Class(:,1)==k & Ses(i).Celltype(:,1)==1;
+        DMSselection=Ses(i).Coord(:,4)==20 & Ses(i).Class(:,1)==k & Ses(i).Celltype(:,1)==1;
+        
+        DLSconcat_ACQ=[];DMSconcat_ACQ=[];
+        DLSconcat_ACQ = cat(2,DLSconcat_ACQ,Ses(i).Ev(7).PSTHz(DLSselection,Ishow));
+        DMSconcat_ACQ = cat(2,DMSconcat_ACQ,Ses(i).Ev(7).PSTHz(DMSselection,Ishow));
+        for j=1:6
+            DLSconcat_ACQ = cat(2,DLSconcat_ACQ,Ses(i).Ev(j).PSTHz(DLSselection,Ishow));
+            DMSconcat_ACQ = cat(2,DMSconcat_ACQ,Ses(i).Ev(j).PSTHz(DMSselection,Ishow));
+        end
+        MeanDLS(k).img(i-3,:)=nanmean(DLSconcat_ACQ,1);
+        MeanDMS(k).img(i-3,:)=nanmean(DMSconcat_ACQ,1);
+        if i>10
+            DLSconcat_ACQ_3ds(k).DLS=cat(1,DLSconcat_ACQ_3ds(k).DLS,DLSconcat_ACQ);
+            DMSconcat_ACQ_3ds(k).DMS=cat(1,DMSconcat_ACQ_3ds(k).DMS,DMSconcat_ACQ);
+            
+            DLSsem(k,:)=nanste(DLSconcat_ACQ(:,time),1);
+            DMSsem(k,:)=nanste(DMSconcat_ACQ(:,time),1);
+        end
+    end
+end
+
+%% Proportion of neurons in phasic vs non pHasic classes
+
+for j=1:2
+    for i=4:13
+        DLSselectionP=Ses(i).Coord(:,4)==BrainRegion(1) & (Ses(i).Class(:,1)<4) & Ses(i).Celltype(:,1)==1 ;
+        DMSselectionP=Ses(i).Coord(:,4)==BrainRegion(2) & (Ses(i).Class(:,1)<4) & Ses(i).Celltype(:,1)==1 ;
+
+        if j==1;DSselection=logical(DLSselectionP); else DSselection=logical(DMSselectionP);end
+        
+        DS_Nb_tot_phasic(i-3,j)=sum(DSselection);
+        PercentPhasic(i-3,1+(j-1)*2)=100*sum(DSselection)/sum(Ses(i).Coord(:,4)==BrainRegion(j) & Ses(i).Celltype(:,1)==1);
+        PercentPhasic(i-3,2+(j-1)*2)=100-PercentPhasic(i-3,1+(j-1)*2);
+    end   
+end
+
+subplot(5,3,2) % plot DLS phasic vs non phasic
+bp1=barh(PercentPhasic(:,1:2),1,'stacked');
+axis([0 100 0.5 10.5])
+xlabel('% DLS neurons')
+ylabel('sessions')
+subplot(5,3,3) % plot DMS phasic vs non phasic
+bp2=barh(PercentPhasic(:,3:4),1,'stacked');
+axis([0 100 0.5 10.5])
+xlabel('% DMS neurons')
+ylabel('sessions')
+
+%% heat maps classes of neurons
+Eventnames={'';'LI';'';'LP1';'';'LP2';'';'LP3';'';'LP4';'';'LP5';'';'PE'};
+titlegraphDLS={'DLS - Start','DLS - Stop', 'DLS - middle','DLS - INH','DLS - EXC'};
+titlegraphDMS={'DMS - Start','DMS - Stop', 'DMS - middle','DMS - INH','DMS - EXC'};
+%titlePSTH={'Inhibited neurons', 'Excited neurons'};
+
+BrainRegion=[10 20];
+figure
+nplot=0;
+
+for k=1:max(ClassACQ(:,2)) % loop thru non phasic classes 
+    if k<4 || k==5
+        nplot=nplot+1;
+    elseif k==4
+        nplot=nplot+10;
+    end
+    subplot(6,4,2+(nplot-1))
+    imagesc(time,[1 size(MeanDLS(k).img,1)],MeanDLS(k).img,Clim); %colorbar;axis tight;
+    colormap(jet);
+    %colormap(plasma);
+    set(gca,'XTick',[0:25.5:357]);
+    set(gca,'xticklabel',Eventnames)
+    title(titlegraphDLS(k));
+    ylabel('sessions');
+    hold on
+    plot([51 51],[size(MeanDLS(k).img,1)+0.5 0.5],'k')
+    plot([102 102],[size(MeanDLS(k).img,1)+0.5 0.5],'k')
+    plot([153 153],[size(MeanDLS(k).img,1)+0.5 0.5],'k')
+    plot([204 204],[size(MeanDLS(k).img,1)+0.5 0.5],'k')
+    plot([255 255],[size(MeanDLS(k).img,1)+0.5 0.5],'k')
+    plot([306 306],[size(MeanDLS(k).img,1)+0.5 0.5],'k')
+    hold off
+    
+    
+    subplot(6,4,6+(nplot-1))
+    imagesc(time,[1 size(MeanDMS(k).img,1)],MeanDMS(k).img,Clim); %colorbar;axis tight;
+    colormap(jet);
+    set(gca,'XTick',[0:25.5:357]);
+    set(gca,'xticklabel',Eventnames)
+    title(titlegraphDMS(k));
+    ylabel('sessions');
+    hold on
+    plot([51 51],[size(MeanDMS(k).img,1)+0.5 0.5],'k')
+    plot([102 102],[size(MeanDMS(k).img,1)+0.5 0.5],'k')
+    plot([153 153],[size(MeanDMS(k).img,1)+0.5 0.5],'k')
+    plot([204 204],[size(MeanDMS(k).img,1)+0.5 0.5],'k')
+    plot([255 255],[size(MeanDMS(k).img,1)+0.5 0.5],'k')
+    plot([306 306],[size(MeanDMS(k).img,1)+0.5 0.5],'k')
+    hold off
+    
+    subplot(6,4,10+(nplot-1))
+    
+    DLSpsth=nanmean(DLSconcat_ACQ_3ds(k).DLS,1);
+    DLS_sem=nanste(DLSconcat_ACQ_3ds(k).DLS,1);
+    DLSup=DLSpsth+DLS_sem;
+    DLSdown=DLSpsth-DLS_sem;
+    
+    DMSpsth=nanmean(DMSconcat_ACQ_3ds(k).DMS,1);
+    DMS_sem=nanste(DMSconcat_ACQ_3ds(k).DMS,1);
+    DMSup=DMSpsth+DMS_sem;
+    DMSdown=DMSpsth-DMS_sem;
+    
+    hold on;
+    patch([time,time(end:-1:1)],[DLSup,DLSdown(end:-1:1)],myblue,'EdgeColor','none');alpha(0.5);
+    plot(time,DLSpsth,'Color',myblue,'linewidth',1.5);
+    patch([time,time(end:-1:1)],[DMSup,DMSdown(end:-1:1)],'r','EdgeColor','none');alpha(0.5);
+    plot(time,DMSpsth,'r','linewidth',1.5);
+    
+    yaxis=[round(min(DMSdown)-0.5,0) round(max(DMSup)+0.5,0)];
+    plot([0 357], [0 0],'LineStyle',':','Color','k');
+    plot([51 51],yaxis,'LineStyle',':','Color','k')
+    plot([102 102],yaxis,'LineStyle',':','Color','k')
+    plot([153 153],yaxis,'LineStyle',':','Color','k')
+    plot([204 204],yaxis,'LineStyle',':','Color','k')
+    plot([255 255],yaxis,'LineStyle',':','Color','k')
+    plot([306 306],yaxis,'LineStyle',':','Color','k')
+    axis([Xaxis,yaxis]);
+    set(gca,'XTick',0:25.5:357);
+    set(gca,'xticklabel',Eventnames)
+    ylabel('z-score');
+    hold off
+end
+
+
+%% Proportion of neurons
+for k=1:3
+    for i=4:13   
+        DLSselection=Ses(i).Coord(:,4)==BrainRegion(1) & (Ses(i).Class(:,1)==k) & Ses(i).Celltype(:,1)==1;
+        DMSselection=Ses(i).Coord(:,4)==BrainRegion(2) & (Ses(i).Class(:,1)==k) & Ses(i).Celltype(:,1)==1;
+        
+        DLS_prop_Phasic(i-3,k)=100*sum(DLSselection)/(sum(Ses(i).Coord(:,4)==BrainRegion(1) & Ses(i).Celltype(:,1)==1 & Ses(i).Class(:,1)<4 & Ses(i).Class(:,1)~=0));
+        DMS_prop_Phasic(i-3,k)=100*sum(DMSselection)/(sum(Ses(i).Coord(:,4)==BrainRegion(2) & Ses(i).Celltype(:,1)==1 & Ses(i).Class(:,1)<4 & Ses(i).Class(:,1)~=0));
+        
+        DLSselectionNumber(i-3,k)=sum(DLSselection);
+        DMSselectionNumber(i-3,k)=sum(DMSselection);
+    end
+end
+
+subplot(6,4,1)
+b1=barh(DLS_prop_Phasic,1,'stacked');
+axis([0 100 0.5 10.5])
+xlabel('% DLS neurons')
+ylabel('sessions')
+
+subplot(6,4,5)
+b2=barh(DMS_prop_Phasic,1,'stacked');
+axis([0 100 0.5 10.5])
+xlabel('% DMS neurons')
+ylabel('sessions')
+
+
+for k=1:2
+    for i=4:13   
+        DLSselection=Ses(i).Coord(:,4)==BrainRegion(1) & (Ses(i).Class(:,1)==k+3) & Ses(i).Celltype(:,1)==1;
+        DMSselection=Ses(i).Coord(:,4)==BrainRegion(2) & (Ses(i).Class(:,1)==k+3) & Ses(i).Celltype(:,1)==1;
+        
+        DLS_prop_nonPhasic(i-3,k)=100*sum(DLSselection)/(sum(Ses(i).Coord(:,4)==BrainRegion(1) & Ses(i).Celltype(:,1)==1 & Ses(i).Class(:,1)>3 & Ses(i).Class(:,1)~=0));
+        DMS_prop_nonPhasic(i-3,k)=100*sum(DMSselection)/(sum(Ses(i).Coord(:,4)==BrainRegion(2) & Ses(i).Celltype(:,1)==1 & Ses(i).Class(:,1)>3 & Ses(i).Class(:,1)~=0));
+        
+        DLSselectionNumber(i-3,k)=sum(DLSselection);
+        DMSselectionNumber(i-3,k)=sum(DMSselection);
+    end
+end
+
+subplot(6,4,13)
+b1=barh(DLS_prop_nonPhasic,1,'stacked');
+axis([0 100 0.5 10.5])
+xlabel('% DLS neurons')
+ylabel('sessions')
+
+subplot(6,4,17)
+b2=barh(DMS_prop_nonPhasic,1,'stacked');
+axis([0 100 0.5 10.5])
+xlabel('% DMS neurons')
+ylabel('sessions')
+
+%% %% make table for statistics
+
+selectDT5Session=sessionID_allsession(:,1)>3;
+
+table_ACQ=cat(2,sessionID_allsession(selectDT5Session,1),Coord_allsession(selectDT5Session,1),Celltype(selectDT5Session,2),ClassACQ(selectDT5Session,1),ClassACQ(selectDT5Session,2),DSmeanz_ACQ_allsession(selectDT5Session,:),TRN_allsession(selectDT5Session,1),EXCINH_allsession(:,1), RatID_ACQ(selectDT5Session,1));
+save('table_ACQ.mat','table_ACQ')
+
+
+nbevent=[1:12];
+
+selectionSTAT=table_ACQ(:,2)~=0 & table_ACQ(:,3)==1 & table_ACQ(:,18)~=0;
+sessionCat=categorical(table_ACQ(selectionSTAT,1));
+tableSTAT=array2table(table_ACQ(selectionSTAT,1:20),'VariableNames',{'session','region','celltype','phasicID','class','win1','win2','win3','win4','win5','win6','win7','win8','win9','win10','win11','win12','TRN','EXCINH','habit'});
+tableSTAT.session=categorical(tableSTAT.session);
+tableSTAT.region=categorical(tableSTAT.region);
+tableSTAT.habit=categorical(tableSTAT.habit);
+ 
+ rm = fitrm(tableSTAT,'win1-win12~region*session','WithinDesign',nbevent');
+ ranovatbl_event = ranova(rm);   
+ Btw_ranovatbl_event = anova(rm);
+ SphericityAssumption=mauchly(rm);
+ posthoc = multcompare(rm,'session', 'by', 'region');
+ 
+ 
+ rmh = fitrm(tableSTAT,'win1-win12~region*session*habit','WithinDesign',nbevent');
+ ranovatbl_event_habit = ranova(rmh);   
+ Btw_ranovatbl_event_habit = anova(rmh);
+ SphericityAssumption_habit=mauchly(rmh);
+ posthoc_habit = multcompare(rmh,'region', 'by', 'habit');
+
+ 
+%  for i=1:max(table_ACQ(:,5))
+for i=1:3 %loop thru phasic neurons, problem, not enough neurons in non phasic classes
+    selectionSTAT1=table_ACQ(:,2)~=0 & table_ACQ(:,3)==1 & table_ACQ(:,5)==i; 
+    
+    tableSTAT1=array2table(table_ACQ(selectionSTAT1,1:20),'VariableNames',{'session','region','celltype','phasicID','class','win1','win2','win3','win4','win5','win6','win7','win8','win9','win10','win11','win12','TRN','EXCINH','habit'});
+    tableSTAT1.session=categorical(tableSTAT1.session);
+    tableSTAT1.region=categorical(tableSTAT1.region);
+    tableSTAT1.habit=categorical(tableSTAT1.habit);
+    
+    rm = fitrm(tableSTAT1,'win1-win12~region*session','WithinDesign',nbevent');
+    stat(i).ranovatbl_event = ranova(rm);
+    stat(i).Btw_ranovatbl_event = anova(rm);
+    
+    [~,stat(i).LIevent,~]=anovan(table_ACQ(selectionSTAT1,6),{table_ACQ(selectionSTAT1,2)},'varnames',{'region'}, 'display','off');
+    [~,stat(i).LRevent,~]=anovan(table_ACQ(selectionSTAT1,16),{table_ACQ(selectionSTAT1,2)},'varnames',{'region'}, 'display','off');
+    [~,stat(i).PEevent,~]=anovan(table_ACQ(selectionSTAT1,17),{table_ACQ(selectionSTAT1,2)},'varnames',{'region'}, 'display','off');
+
+    rmh = fitrm(tableSTAT1,'win1-win12~region*session*habit','WithinDesign',nbevent');
+    stat(i).ranovatbl_event_habit = ranova(rmh);
+    stat(i).Btw_ranovatbl_habit = anova(rmh); 
+    stat(i).posthoc_habit = multcompare(rmh,'region', 'by', 'habit');
+end
+
+
+ for i=1:10
+     selectionSTAT2=table_ACQ(:,2)~=0 & table_ACQ(:,3)==1 & table_ACQ(:,1)==i+3 & table_ACQ(:,4)==1; 
+    [chi(i).tbl1,chi(i).chi2,chi(i).p1,chi(i).labels] = crosstab(table_ACQ(selectionSTAT2,2),table_ACQ(selectionSTAT2,5));
+ end
+ 
+ save('Rearlytraining_light.mat','Erefnames','Ses','Tm') 

scripts/H_generateFigure7A_Decoding1.m

@@ -0,0 +1,293 @@
+%% generate figure 7A. Decoding brain region ID without PCA
+
+% code from David Ottenheimer 04/19/2018
+
+%get the data
+clear all
+close all
+
+load('Rearlytraining_light.mat')
+Xaxis1=[-0.25 0.25];
+Ishow=find(Tm>=Xaxis1(1) & Tm<=Xaxis1(2));
+
+
+for i=4:13 %loop thru ACQ session
+
+   DLSselection=Ses(i).Coord(:,4)==10 & Ses(i).Celltype(:,1)==1;
+   DMSselection=Ses(i).Coord(:,4)==20 & Ses(i).Celltype(:,1)==1;
+    
+    DLSconcat_ACQ=[];DMSconcat_ACQ=[];
+    DLSconcat_ACQ = cat(2,DLSconcat_ACQ,Ses(i).Ev(7).Meanz(DLSselection,1));
+    DMSconcat_ACQ = cat(2,DMSconcat_ACQ,Ses(i).Ev(7).Meanz(DMSselection,1));
+    for j=1:5
+        DLSconcat_ACQ = cat(2,DLSconcat_ACQ,Ses(i).Ev(j).MeanzPRE(DLSselection,1));
+        DLSconcat_ACQ = cat(2,DLSconcat_ACQ,Ses(i).Ev(j).Meanz(DLSselection,1));
+        DMSconcat_ACQ = cat(2,DMSconcat_ACQ,Ses(i).Ev(j).MeanzPRE(DMSselection,1));
+        DMSconcat_ACQ = cat(2,DMSconcat_ACQ,Ses(i).Ev(j).Meanz(DMSselection,1));
+    end
+    DLSconcat_ACQ = cat(2,DLSconcat_ACQ,Ses(i).Ev(6).MeanzPRE(DLSselection,1));
+    DMSconcat_ACQ = cat(2,DMSconcat_ACQ,Ses(i).Ev(6).MeanzPRE(DMSselection,1));
+
+    PSTHzDecodeACQ(i-3).DLS=DLSconcat_ACQ(~isnan(mean(DLSconcat_ACQ,2)),:);
+    PSTHzDecodeACQ(i-3).DMS=DMSconcat_ACQ(~isnan(mean(DMSconcat_ACQ,2)),:);
+end
+
+clear Ses
+load('Rextendedtraining_light.mat');
+load('Celltype_extendedTraining.mat');
+
+DLSselection=Coord(:,4)==10 & Celltype(:,1)==1;
+DMSselection=Coord(:,4)==20 & Celltype(:,1)==1;
+
+DLSconcat_OT = Ev(7).Meanz(DLSselection,1);
+DMSconcat_OT = Ev(7).Meanz(DMSselection,1);
+for i=1:5
+    DLSconcat_OT = cat(2,DLSconcat_OT,Ev(i).MeanzPRE(DLSselection,1));
+    DLSconcat_OT = cat(2,DLSconcat_OT,Ev(i).Meanz(DLSselection,1));
+    DMSconcat_OT = cat(2,DMSconcat_OT,Ev(i).MeanzPRE(DMSselection,1));
+    DMSconcat_OT = cat(2,DMSconcat_OT,Ev(i).Meanz(DMSselection,1));
+end
+DLSconcat_OT = cat(2,DLSconcat_OT,Ev(6).MeanzPRE(DLSselection,1));
+DMSconcat_OT = cat(2,DMSconcat_OT,Ev(6).MeanzPRE(DMSselection,1));
+
+MeanzDecodeOT.DLS=DLSconcat_OT;
+MeanzDecodeOT.DMS=DMSconcat_OT;
+
+save('MeanzDecodeACQ','MeanzDecodeOT')
+
+
+
+%% ACQ PCcomponents Comparison sessions. 
+%PCA first
+%first do PCA on an equal number of neurons in each region
+nbsession=0;
+for k=1:10
+    clear rep
+    countRep=0;
+    PSTHzDecode.DLS=PSTHzDecodeACQ(k).DLS;
+    PSTHzDecode.DMS=PSTHzDecodeACQ(k).DMS;
+
+    if length(PSTHzDecode.DLS(:,1))<length(PSTHzDecode.DMS(:,1))
+        PCAneurons=length(PSTHzDecode.DLS(:,1));
+    else
+        PCAneurons=length(PSTHzDecode.DMS(:,1));
+    end
+    
+    for l=1:50 % loop added to account for selection bias in VSel and NSel
+        %pick which neurons to use
+        SetupDMSel=cat(1,ones(PCAneurons,1),zeros(length(PSTHzDecode.DMS(:,1))-PCAneurons,1));
+        VSel=(SetupDMSel(randperm(length(SetupDMSel)))==1);
+        
+        %pick which neurons to use
+        SetupDLSel=cat(1,ones(PCAneurons,1),zeros(length(PSTHzDecode.DLS(:,1))-PCAneurons,1));
+        NSel=(SetupDLSel(randperm(length(SetupDLSel)))==1);
+        
+        rep(l).concat=cat(1,PSTHzDecode.DLS(NSel,:),PSTHzDecode.DMS(VSel,:));
+    end
+    
+    %% Decoding
+    clear Partitions SVMModel prediction actual correct DecodeRegion
+    folds = 10;%PCAneurons; %number of times cross-validated
+    shuffs = 1; %number of times shuffled
+    DiscrimType= 'linear';
+    
+    %setup the event identity matrix for decoding
+    DecodeRegion(1:PCAneurons,1)=1;
+    DecodeRegion((PCAneurons+1):(PCAneurons*2),1)=2;
+    for l=1:50
+        
+            countRep=countRep+1;
+            clear CVacc CVaccSh
+            %Setup spikes matrix for decoding
+            DecodeMeanZ=rep(l).concat;
+            
+            %normal model
+            for r = 1:folds
+                
+                Partitions = cvpartition(DecodeRegion,'KFold',folds);
+                SVMModel = fitcdiscr(DecodeMeanZ(Partitions.training(r),:),DecodeRegion(Partitions.training(r)),'DiscrimType',DiscrimType);
+                prediction = predict(SVMModel,DecodeMeanZ(Partitions.test(r),:));
+                actual = DecodeRegion(Partitions.test(r));
+                correct = prediction - actual;
+                CVacc(r,1) = sum(correct==0) / length(correct);
+                
+            end
+            
+            %shuffled model
+            for q=1:shuffs
+                DecodeRsSh=DecodeRegion(randperm(length(DecodeRegion)));
+                PartitionsSh = cvpartition(DecodeRsSh,'KFold',folds);
+                for s = 1:folds
+                    SVMModelSh = fitcdiscr(DecodeMeanZ(PartitionsSh.training(s),:),DecodeRsSh(PartitionsSh.training(s)),'DiscrimType',DiscrimType);
+                    predictionSh = predict(SVMModelSh,DecodeMeanZ(PartitionsSh.test(s),:));
+                    actualSh = DecodeRsSh(PartitionsSh.test(s));
+                    correctSh = predictionSh - actualSh;
+                    CVaccSh(s,1) = sum(correctSh==0) / length(correctSh);
+                end
+                AccShuff(q,1) = nanmean(CVaccSh);
+            end
+            PCADecodeAccAA{k,1}(countRep,1)=nanmean(CVacc);
+            PCADecodeAccShAA{k,1}(countRep,1)=nanmean(AccShuff);
+            fprintf(['RepSelection #' num2str(countRep) '\n']);
+        end
+    
+end
+
+
+%% plotting
+%colors
+inh=[0 0 0.6];
+exc=[0.8 0 0];
+NumSession=[1 2 3 4 5 6 7 8 9 10];
+
+
+for i=1:10
+    subplot(1,2,1);
+    hold on;
+    errorbar(NumSession(i),nanmean(PCADecodeAccAA{i,1}),nanstd(PCADecodeAccAA{i,1}),'o','color',exc);
+    errorbar(NumSession(i),nanmean(PCADecodeAccShAA{i,1}),nanstd(PCADecodeAccShAA{i,1}),'o','color',inh);
+end
+axis([0 10+1 0.35 0.85]);
+title('Decoding of region across session');
+xlabel('Sessions');
+ylabel('Accuracy');
+legend({'True','Shuffled'},'location','northwest');
+    
+
+
+%% Overtraining as a function of PC and number of neurons included
+
+for k=1:5
+    countRep=0;
+    clear rep
+    
+    PSTHzDecode.DLS=MeanzDecodeOT.DLS;
+    PSTHzDecode.DMS=MeanzDecodeOT.DMS;
+    
+    nbneurons=[30 60 100 200 300];
+    PCAneurons=nbneurons(k);
+    
+    for l=1:50 % loop added to account for selection bias in VSel and NSel
+    SetupDMSel=cat(1,ones(PCAneurons,1),zeros(length(PSTHzDecode.DMS(:,1))-PCAneurons,1));
+    VSel=(SetupDMSel(randperm(length(SetupDMSel)))==1);
+    
+    %pick which neurons to use
+    SetupDLSel=cat(1,ones(PCAneurons,1),zeros(length(PSTHzDecode.DLS(:,1))-PCAneurons,1));
+    NSel=(SetupDLSel(randperm(length(SetupDLSel)))==1);
+    
+    rep(l).concatOT=cat(1,PSTHzDecode.DLS(NSel,:),PSTHzDecode.DMS(VSel,:));
+    end
+    
+    %% Decoding
+    clear Partitions SVMModel prediction actual correct DecodeRegion
+    folds = 10;%PCAneurons; %number of times cross-validated
+    shuffs = 1; %number of times shuffled
+    repetitions=20; %how many times to run the analysis
+    DiscrimType= 'linear';
+    
+    %setup the event identity matrix for decoding
+    DecodeRegion(1:PCAneurons,1)=1;
+    DecodeRegion((PCAneurons+1):(PCAneurons*2),1)=2;
+    for l=1:50
+   
+        countRep=countRep+1;
+
+        clear CVacc CVaccSh
+        %Setup spikes matrix for decoding
+        DecodeMeanZ=rep(l).concatOT;
+        
+        %normal model
+        for r = 1:folds
+            
+            Partitions = cvpartition(DecodeRegion,'KFold',folds);
+            SVMModel = fitcdiscr(DecodeMeanZ(Partitions.training(r),:),DecodeRegion(Partitions.training(r)),'DiscrimType',DiscrimType);
+            prediction = predict(SVMModel,DecodeMeanZ(Partitions.test(r),:));
+            actual = DecodeRegion(Partitions.test(r));
+            correct = prediction - actual;
+            CVacc(r,1) = sum(correct==0) / length(correct);
+            
+        end
+        
+        %shuffled model
+        for q=1:shuffs
+            DecodeRsSh=DecodeRegion(randperm(length(DecodeRegion)));
+            PartitionsSh = cvpartition(DecodeRsSh,'KFold',folds);
+            for s = 1:folds
+                SVMModelSh = fitcdiscr(DecodeMeanZ(PartitionsSh.training(s),:),DecodeRsSh(PartitionsSh.training(s)),'DiscrimType',DiscrimType);
+                predictionSh = predict(SVMModelSh,DecodeMeanZ(PartitionsSh.test(s),:));
+                actualSh = DecodeRsSh(PartitionsSh.test(s));
+                correctSh = predictionSh - actualSh;
+                CVaccSh(s,1) = sum(correctSh==0) / length(correctSh);
+            end
+            AccShuff(q,1) = nanmean(CVaccSh);
+        end
+        PCADecodeAccOT{1,k}(countRep,1)=nanmean(CVacc);
+        PCADecodeAccShOT{1,k}(countRep,1)=nanmean(AccShuff);
+        
+        fprintf(['Rep #' num2str(countRep) '\n']);
+    end
+   
+end
+
+
+%% plotting
+%colors
+inh=[0 0 0.6];
+exc=[0.8 0 0];
+
+Shuffle=[];
+for k=1:5
+    subplot(1,2,2);
+    hold on;
+    errorbar(nbneurons(k),nanmean(PCADecodeAccOT{1,k}),nanstd(PCADecodeAccOT{1,k}),'o','color',exc(1,:));
+    errorbar(nbneurons(k),nanmean(PCADecodeAccShOT{1,k}),nanstd(PCADecodeAccShOT{1,k}),'o','color',inh(1,:));
+end
+
+axis([0.5 300 0.35 0.85]);
+title('Ext training: Decoding DMS vs DLS');
+xlabel('Ext Training - ensemble size');
+ylabel('Accuracy');
+
+%save('Decode_withoutPCA.mat','PCADecodeAccAA','PCADecodeAccOT','PCADecodeAccShAA','PCADecodeAccShOT');
+
+%% stat
+for i=1:10
+   tableAcc(:,i)=cat(1,PCADecodeAccAA{i,1},PCADecodeAccShAA{i,1});
+   Between_factor=cat(1,zeros(50,1),ones(50,1));
+   [p,tbl]=anova1(tableAcc(:,i),Between_factor,'off');
+   p_value(i)=p;
+   stat(i).t=tbl;
+end
+
+for i=1:5
+    tableAccOT(:,i)=cat(1,PCADecodeAccOT{1,i},PCADecodeAccShOT{1,i});
+    [p,tbl]=anova1(tableAccOT(:,i),Between_factor,'off');
+    pOT_value(i)=p;
+    statOT(i).t=tbl;
+end
+
+tableAccACQ_OT=cat(1,PCADecodeAccAA{10,1},PCADecodeAccOT{1,1});
+[p_ACQ_OT,tbl_ACQ_OT]=anova1(tableAccACQ_OT,Between_factor,'off');
+
+%%
+
+for i=1:length(PCADecodeAccAA)
+    allsamples=cat(1,PCADecodeAccAA{i,1},PCADecodeAccShAA{i,1});
+    for p=1:10000
+        shuffsamples=allsamples(randperm(length(allsamples)));
+        shuffdiff(p,1)=abs(mean(shuffsamples(1:length(allsamples)/2))-mean(shuffsamples(length(allsamples)/2+1:end)));
+    end
+    diff=abs(nanmean(PCADecodeAccAA{i,1})-nanmean(PCADecodeAccShAA{i,1}));
+    trainingpval(i,1)=(sum(shuffdiff>diff)+1)/(length(shuffdiff)+1);
+end
+
+%% overtraining
+for i=1:length(PCADecodeAccOT)
+    allsamples=cat(1,PCADecodeAccOT{1,i},PCADecodeAccShOT{1,i});
+    for p=1:10000
+        shuffsamples=allsamples(randperm(length(allsamples)));
+        shuffdiff(p,1)=abs(mean(shuffsamples(1:length(allsamples)/2))-mean(shuffsamples(length(allsamples)/2+1:end)));
+    end
+    diff=abs(nanmean(PCADecodeAccOT{1,i})-nanmean(PCADecodeAccShOT{1,i}));
+    OTpval(1,i)=(sum(shuffdiff>diff)+1)/(length(shuffdiff)+1);
+end

scripts/H_generateFigure7B_Decoding2.m

@@ -0,0 +1,408 @@
+%% Generate figure 7B and figure 7- supplement figure 1
+
+% decoding using PCA; Up left, compare 3 early training sessions across PCs 
+% Up right, compare different ensemble size of extended training across PCs 
+% code from David Ottenheimer 04/19/2018
+
+clear all
+close all
+
+load('Rearlytraining_light.mat')
+Xaxis1=[-0.25 0.25];
+Ishow=find(Tm>=Xaxis1(1) & Tm<=Xaxis1(2));
+
+
+for i=4:13 %loop thru ACQ session
+
+   DLSselection=Ses(i).Coord(:,4)==10 & Ses(i).Celltype(:,1)==1;
+   DMSselection=Ses(i).Coord(:,4)==20 & Ses(i).Celltype(:,1)==1;
+    
+    DLSconcat_ACQ=[];DMSconcat_ACQ=[];
+    DLSconcat_ACQ = cat(2,DLSconcat_ACQ,Ses(i).Ev(7).PSTHz(DLSselection,Ishow));
+    DMSconcat_ACQ = cat(2,DMSconcat_ACQ,Ses(i).Ev(7).PSTHz(DMSselection,Ishow));
+    for j=1:6
+        DLSconcat_ACQ = cat(2,DLSconcat_ACQ,Ses(i).Ev(j).PSTHz(DLSselection,Ishow));
+        DMSconcat_ACQ = cat(2,DMSconcat_ACQ,Ses(i).Ev(j).PSTHz(DMSselection,Ishow));
+    end
+
+    PSTHzDecodeACQ(i-3).DLS=DLSconcat_ACQ(~isnan(mean(DLSconcat_ACQ,2)),:);
+    PSTHzDecodeACQ(i-3).DMS=DMSconcat_ACQ(~isnan(mean(DMSconcat_ACQ,2)),:);
+
+end
+
+clear Ses
+load('Rextendedtraining_light.mat');
+load('Celltype_extendedtraining.mat');
+
+DLSselection=Coord(:,4)==10 & Celltype(:,1)==1;
+DMSselection=Coord(:,4)==20 & Celltype(:,1)==1;
+
+DLSconcat_OT = Ev(7).PSTHz(DLSselection,Ishow);
+DMSconcat_OT = Ev(7).PSTHz(DMSselection,Ishow);
+for i=1:6
+    DLSconcat_OT = cat(2,DLSconcat_OT,Ev(i).PSTHz(DLSselection,Ishow));
+    DMSconcat_OT = cat(2,DMSconcat_OT,Ev(i).PSTHz(DMSselection,Ishow));
+end
+PSTHzDecodeOT.DLS=DLSconcat_OT;
+PSTHzDecodeOT.DMS=DMSconcat_OT;
+
+save('PSTHzDecodeACQ','PSTHzDecodeOT')
+
+
+
+%% ACQ PCcomponents Comparison sessions. 
+%PCA first
+%first do PCA on an equal number of neurons in each region
+nbsession=0;
+for k=1:3:10
+    countRep=0;
+    nbsession=nbsession+1;
+    
+    PSTHzDecode.DLS=PSTHzDecodeACQ(k).DLS;
+    PSTHzDecode.DMS=PSTHzDecodeACQ(k).DMS;
+    
+    if length(PSTHzDecode.DLS(:,1))<length(PSTHzDecode.DMS(:,1))
+        PCAneurons=length(PSTHzDecode.DLS(:,1));
+    else
+        PCAneurons=length(PSTHzDecode.DMS(:,1));
+    end
+    
+    for l=1:50 % loop added to account for selection bias in VSel and NSel
+        %pick which neurons to use
+        SetupDMSel=cat(1,ones(PCAneurons,1),zeros(length(PSTHzDecode.DMS(:,1))-PCAneurons,1));
+        VSel=(SetupDMSel(randperm(length(SetupDMSel)))==1);
+        
+        %pick which neurons to use
+        SetupDLSel=cat(1,ones(PCAneurons,1),zeros(length(PSTHzDecode.DLS(:,1))-PCAneurons,1));
+        NSel=(SetupDLSel(randperm(length(SetupDLSel)))==1);
+        
+        PCAneuronsACQ(nbsession,1)=PCAneurons;
+        
+        concat=cat(1,PSTHzDecode.DLS(NSel,:),PSTHzDecode.DMS(VSel,:))';
+        [coeff,score,~,~,explained] = pca(concat);
+        
+        
+        % Decoding
+        clear Partitions SVMModel prediction actual correct
+        folds = 10;%PCAneurons; %number of times cross-validated
+        shuffs = 1; %number of times shuffled
+        NumCoeffs = [1 2 3 4 5 6 7 8 9 10]; %number of PCs used in decoding
+        repetitions=20; %how many times to run the analysis
+        DiscrimType= 'linear';
+        
+        %setup the event identity matrix for decoding
+        DecodeRegion(1:PCAneurons,1)=1;
+        DecodeRegion((PCAneurons+1):(PCAneurons*2),1)=2;
+        
+        
+        countRep=countRep+1;
+        
+        for i=1:length(NumCoeffs)
+            clear CVacc CVaccSh
+            %Setup spikes matrix for decoding
+            DecodePCs=coeff(:,1:NumCoeffs(i));
+            %ADD IN PCA HERE???
+            
+            %normal model
+            for r = 1:folds
+                
+                Partitions = cvpartition(DecodeRegion,'KFold',folds);
+                SVMModel = fitcdiscr(DecodePCs(Partitions.training(r),:),DecodeRegion(Partitions.training(r)),'DiscrimType',DiscrimType);
+                prediction = predict(SVMModel,DecodePCs(Partitions.test(r),:));
+                actual = DecodeRegion(Partitions.test(r));
+                correct = prediction - actual;
+                CVacc(r,1) = sum(correct==0) / length(correct);
+                
+            end
+            
+            %shuffled model
+            for q=1:shuffs
+                DecodeRsSh=DecodeRegion(randperm(length(DecodeRegion)));
+                PartitionsSh = cvpartition(DecodeRsSh,'KFold',folds);
+                for s = 1:folds
+                    SVMModelSh = fitcdiscr(DecodePCs(PartitionsSh.training(s),:),DecodeRsSh(PartitionsSh.training(s)),'DiscrimType',DiscrimType);
+                    predictionSh = predict(SVMModelSh,DecodePCs(PartitionsSh.test(s),:));
+                    actualSh = DecodeRsSh(PartitionsSh.test(s));
+                    correctSh = predictionSh - actualSh;
+                    CVaccSh(s,1) = sum(correctSh==0) / length(correctSh);
+                end
+                AccShuff(q,1) = nanmean(CVaccSh);
+            end
+            PCADecodeAccACQ{i,nbsession}(countRep,1)=nanmean(CVacc);
+            PCADecodeAccShACQ{i,nbsession}(countRep,1)=nanmean(AccShuff);
+        end
+        fprintf(['RepACQ #' num2str(countRep) '\n']);
+    end  
+end
+
+%plotting
+%colors
+inh=[0 0.3333 1; 0 0.5 1; 0 0.67 1; 0 0.8353 1];
+exc=[1 0.1647 0; 0.9020 0.4510 0; 0.9020 0.6 0; 1 0.8353 0];
+
+%Y = prctile(X,p)
+
+for i=1:length(NumCoeffs)
+    for k=1:4
+        subplot(2,2,1);
+        hold on;
+        errorbar(NumCoeffs(i),nanmean(PCADecodeAccACQ{i,k}),nanstd(PCADecodeAccACQ{i,k}),'o','color',exc(k,:));
+        errorbar(NumCoeffs(i),nanmean(PCADecodeAccShACQ{i,4}),nanstd(PCADecodeAccShACQ{i,4}),'o','color',inh(1,:));
+    end
+end
+
+axis([0 NumCoeffs(end)+1 0.30 0.85]);
+title('Early training: Decoding DMS vs DLS');
+xlabel('Number of PCs');
+ylabel('Accuracy');
+
+%% Overtraining as a function of PC and number of neurons included
+
+for k=1:5
+    countRep=0;
+    
+    PSTHzDecode.DLS=PSTHzDecodeOT.DLS;
+    PSTHzDecode.DMS=PSTHzDecodeOT.DMS;
+    
+    nbneurons=[30 60 100 200 300];
+    PCAneurons=nbneurons(k);
+    
+    for l=1:50 % loop added to account for selection bias in VSel and NSel
+        %pick which neurons to use
+        SetupDMSel=cat(1,ones(PCAneurons,1),zeros(length(PSTHzDecode.DMS(:,1))-PCAneurons,1));
+        VSel=(SetupDMSel(randperm(length(SetupDMSel)))==1);
+        
+        %pick which neurons to use
+        SetupDLSel=cat(1,ones(PCAneurons,1),zeros(length(PSTHzDecode.DLS(:,1))-PCAneurons,1));
+        NSel=(SetupDLSel(randperm(length(SetupDLSel)))==1);
+        
+        concat=cat(1,PSTHzDecode.DLS(NSel,:),PSTHzDecode.DMS(VSel,:))';
+        [coeff,score,~,~,explained] = pca(concat);
+        
+        %% Decoding
+        clear Partitions SVMModel prediction actual correct DecodeRegion
+        folds = 10;%PCAneurons; %number of times cross-validated
+        shuffs = 1; %number of times shuffled
+        NumCoeffs = [1 2 3 4 5 6 7 8 9 10]; %number of PCs used in decoding
+        DiscrimType= 'linear';
+        
+        %setup the event identity matrix for decoding
+        DecodeRegion(1:PCAneurons,1)=1;
+        DecodeRegion((PCAneurons+1):(PCAneurons*2),1)=2;
+        
+        countRep=countRep+1;
+        
+        for i=1:length(NumCoeffs)
+            
+            clear CVacc CVaccSh
+            %Setup spikes matrix for decoding
+            DecodePCs=coeff(:,1:NumCoeffs(i));
+            %ADD IN PCA HERE???
+            
+            %normal model
+            for r = 1:folds
+                
+                Partitions = cvpartition(DecodeRegion,'KFold',folds);
+                SVMModel = fitcdiscr(DecodePCs(Partitions.training(r),:),DecodeRegion(Partitions.training(r)),'DiscrimType',DiscrimType);
+                prediction = predict(SVMModel,DecodePCs(Partitions.test(r),:));
+                actual = DecodeRegion(Partitions.test(r));
+                correct = prediction - actual;
+                CVacc(r,1) = sum(correct==0) / length(correct);
+                
+            end
+            
+            %shuffled model
+            for q=1:shuffs
+                DecodeRsSh=DecodeRegion(randperm(length(DecodeRegion)));
+                PartitionsSh = cvpartition(DecodeRsSh,'KFold',folds);
+                for s = 1:folds
+                    SVMModelSh = fitcdiscr(DecodePCs(PartitionsSh.training(s),:),DecodeRsSh(PartitionsSh.training(s)),'DiscrimType',DiscrimType);
+                    predictionSh = predict(SVMModelSh,DecodePCs(PartitionsSh.test(s),:));
+                    actualSh = DecodeRsSh(PartitionsSh.test(s));
+                    correctSh = predictionSh - actualSh;
+                    CVaccSh(s,1) = sum(correctSh==0) / length(correctSh);
+                end
+                AccShuff(q,1) = nanmean(CVaccSh);
+            end
+            PCADecodeAccOT{i,k}(countRep,1)=nanmean(CVacc);
+            PCADecodeAccShOT{i,k}(countRep,1)=nanmean(AccShuff);
+        end
+        fprintf(['RepOT #' num2str(countRep) '\n']);
+    end
+end
+
+
+%% plotting
+%colors
+inh=[0 0 0.6; 0 0.3333 1; 0 0.5 1; 0 0.67 1; 0 0.8353 1];
+exc=[0.8 0 0; 1 0.1647 0; 0.9020 0.4510 0; 0.9020 0.6 0; 1 0.8353 0];
+
+
+for i=1:length(NumCoeffs)
+    for k=1:5
+        subplot(2,2,2);
+        hold on;
+        errorbar(NumCoeffs(i),nanmean(PCADecodeAccOT{i,k}),nanstd(PCADecodeAccOT{i,k}),'o','color',exc(k,:));
+        errorbar(NumCoeffs(i),nanmean(PCADecodeAccShOT{i,5}),nanstd(PCADecodeAccShOT{i,5}),'o','color',inh(1,:));
+    end
+end
+
+axis([0 NumCoeffs(end)+1 0.30 0.85]);
+title('Ext Training: Decoding DMS vs DLS');
+xlabel('Number of PCs');
+ylabel('Accuracy');
+
+
+
+%% Decoding DMS vs DLS across session, based on weight of 3 first PCs
+%PCA first
+for i=1:10
+    countRep=0;
+    
+    PSTHzDecode.DLS=PSTHzDecodeACQ(i).DLS;
+    PSTHzDecode.DMS=PSTHzDecodeACQ(i).DMS;
+    
+    if length(PSTHzDecode.DLS(:,1))<length(PSTHzDecode.DMS(:,1))
+        PCAneurons=length(PSTHzDecode.DLS(:,1));
+    else
+        PCAneurons=length(PSTHzDecode.DMS(:,1));
+    end
+    
+    for l=1:50
+        SetupDMSel=cat(1,ones(PCAneurons,1),zeros(length(PSTHzDecode.DMS(:,1))-PCAneurons,1));
+        VSel=(SetupDMSel(randperm(length(SetupDMSel)))==1);
+        
+        %pick which neurons to use
+        SetupDLSel=cat(1,ones(PCAneurons,1),zeros(length(PSTHzDecode.DLS(:,1))-PCAneurons,1));
+        NSel=(SetupDLSel(randperm(length(SetupDLSel)))==1);
+        
+        concat=cat(1,PSTHzDecode.DLS(NSel,:),PSTHzDecode.DMS(VSel,:))';
+        [coeff,score,~,~,explained] = pca(concat);
+        
+        %% Decoding
+        clear Partitions SVMModel prediction actual correct DecodeRegion
+        folds = 10;%PCAneurons; %number of times cross-validated
+        shuffs = 1; %number of times shuffled
+        NumCoeffs = 3; %number of PCs used in decoding
+        DiscrimType= 'linear';
+        
+        %setup the event identity matrix for decoding
+        DecodeRegion(1:PCAneurons,1)=1;
+        DecodeRegion((PCAneurons+1):(PCAneurons*2),1)=2;
+        
+        countRep=countRep+1;
+        
+        clear CVacc CVaccSh
+        %Setup spikes matrix for decoding
+        DecodePCs=coeff(:,1:NumCoeffs);
+        %ADD IN PCA HERE???
+        
+        %normal model
+        for r = 1:folds
+            
+            Partitions = cvpartition(DecodeRegion,'KFold',folds);
+            SVMModel = fitcdiscr(DecodePCs(Partitions.training(r),:),DecodeRegion(Partitions.training(r)),'DiscrimType',DiscrimType);
+            prediction = predict(SVMModel,DecodePCs(Partitions.test(r),:));
+            actual = DecodeRegion(Partitions.test(r));
+            correct = prediction - actual;
+            CVacc(r,1) = sum(correct==0) / length(correct);
+            
+        end
+        
+        %shuffled model
+        for q=1:shuffs
+            DecodeRsSh=DecodeRegion(randperm(length(DecodeRegion)));
+            PartitionsSh = cvpartition(DecodeRsSh,'KFold',folds);
+            for s = 1:folds
+                SVMModelSh = fitcdiscr(DecodePCs(PartitionsSh.training(s),:),DecodeRsSh(PartitionsSh.training(s)),'DiscrimType',DiscrimType);
+                predictionSh = predict(SVMModelSh,DecodePCs(PartitionsSh.test(s),:));
+                actualSh = DecodeRsSh(PartitionsSh.test(s));
+                correctSh = predictionSh - actualSh;
+                CVaccSh(s,1) = sum(correctSh==0) / length(correctSh);
+            end
+            AccShuff(q,1) = nanmean(CVaccSh);
+        end
+        PCADecodeAccAA{i,1}(countRep,1)=nanmean(CVacc);
+        PCADecodeAccShAA{i,1}(countRep,1)=nanmean(AccShuff);
+        fprintf(['RepACQses #' num2str(countRep) '\n']);
+    end
+end
+
+
+%% plotting
+%colors
+inh=[0 0 0.6; 0 0.3333 1; 0 0.5 1; 0 0.67 1; 0 0.8353 1];
+exc=[0.8 0 0; 1 0.1647 0; 0.9020 0.4510 0; 0.9020 0.6 0; 1 0.8353 0];
+NumSession=[1 2 3 4 5 6 7 8 9 10];
+
+for i=1:10
+    subplot(2,2,3);
+    hold on;
+    errorbar(NumSession(i),nanmean(PCADecodeAccAA{i,1}),nanstd(PCADecodeAccAA{i,1}),'o','color',exc(1,:));
+    errorbar(NumSession(i),nanmean(PCADecodeAccShAA{i,1}),nanstd(PCADecodeAccShAA{i,1}),'o','color',inh(1,:));
+end
+axis([0 10+1 0.35 0.85]);
+title('Decoding of region across session');
+xlabel('Sessions');
+ylabel('Accuracy');
+
+
+for k=1:5
+subplot(2,2,4);
+hold on;
+errorbar(nbneurons(k),nanmean(PCADecodeAccOT{3,k}),nanstd(PCADecodeAccOT{3,k}),'o','color',exc(1,:));
+errorbar(nbneurons(k),nanmean(PCADecodeAccShOT{3,k}),nanstd(PCADecodeAccShOT{3,k}),'o','color',inh(1,:));
+end
+axis([0 300 0.35 0.85]);
+xlabel('OT');
+ylabel('Accuracy');
+
+save('PCADecode.mat','PCADecodeAccAA','PCADecodeAccACQ','PCADecodeAccOT','PCADecodeAccShAA','PCADecodeAccShACQ','PCADecodeAccShOT');
+
+%% stat
+for i=1:10
+   tableAcc(:,i)=cat(1,PCADecodeAccAA{i,1},PCADecodeAccShAA{i,1});
+   Between_factor=cat(1,zeros(50,1),ones(50,1));
+   [p,tbl]=anova1(tableAcc(:,i),Between_factor,'off');
+   p_value(i)=p;
+   stat(i).t=tbl;
+end
+
+for i=1:5
+    tableAccOT(:,i)=cat(1,PCADecodeAccOT{3,i},PCADecodeAccShOT{3,i});
+    [p,tbl]=anova1(tableAccOT(:,i),Between_factor,'off');
+    pOT_value(i)=p;
+    statOT(i).t=tbl;
+end
+
+tableAccACQ_OT=cat(1,PCADecodeAccAA{10,1},PCADecodeAccOT{3,1});
+[p_ACQ_OT,tbl_ACQ_OT]=anova1(tableAccACQ_OT,Between_factor,'off');
+
+%% Permutation test
+for i=1:length(PCADecodeAccAA)
+    allsamples=cat(1,PCADecodeAccAA{i,1},PCADecodeAccShAA{i,1});
+    for p=1:10000
+        shuffsamples=allsamples(randperm(length(allsamples)));
+        shuffdiff(p,1)=abs(mean(shuffsamples(1:length(allsamples)/2))-mean(shuffsamples(length(allsamples)/2+1:end)));
+    end
+    diff=abs(nanmean(PCADecodeAccAA{i,1})-nanmean(PCADecodeAccShAA{i,1}));
+    trainingpval(i,1)=(sum(shuffdiff>diff)+1)/(length(shuffdiff)+1);
+end
+
+%% overtraining
+figure;
+
+for i=1:size(PCADecodeAccOT,2)
+    allsamples=cat(1,PCADecodeAccOT{3,i},PCADecodeAccShOT{3,i});
+    for p=1:10000
+        shuffsamples=allsamples(randperm(length(allsamples)));
+        shuffdiff(p,1)=abs(mean(shuffsamples(1:length(allsamples)/2))-mean(shuffsamples(length(allsamples)/2+1:end)));
+    end
+    diff=abs(nanmean(PCADecodeAccOT{1,i})-nanmean(PCADecodeAccShOT{1,i}));
+    OTpval(1,i)=(sum(shuffdiff>diff)+1)/(length(shuffdiff)+1);
+    
+    subplot(5,1,i)
+    hist(shuffdiff);
+    hold on;
+    plot([diff diff],[0 4000]);
+end

scripts/I_generateFigure8_SensitivityDeval.m

@@ -0,0 +1,287 @@
+%% generate figure 8
+
+clc;clear;close all;
+load('Rextendedtraining_light.mat');
+load('RatID_extendedtraining.mat'); %% column 1 = rat ID, column 2 = habit ID; 0=goal-directed / 1=habit
+load('Celltype_extendedTraining.mat');
+load('TRN_DT5_extendedtrain.mat')
+load('table_OT.mat')
+
+%% --- --- Plotting color-coded and average PSTHs --- ---
+BrainRegion=[10 20]; %10 for DLS, 20 for DMS
+Xaxis=[0 357];Yaxis=[-2 5];
+Xaxis1=[-0.25 0.25];
+Ishow=find(Tm>=Xaxis1(1) & Tm<=Xaxis1(2));
+
+% ----------- COLORS ---------
+myblue=[0 113/255 188/255];
+mypurple=[200/255 0 255/255];
+mydarkblue=[0 0/255 255/255];
+TickSize=[0.015 0.02];
+
+Eventnames={'';'LI';'';'LP1';'';'LP2';'';'LP3';'';'LP4';'';'LP5';'';'PE'};
+titlegraph={'Extended Training - sensitive', 'Extended Training - insensitive'};
+time=-0.25:0.25:3.25;time2=1:1:357;
+
+%% Extended training - Average psTH
+%analysis restricted on putative MSN. nonTRN included. 
+
+
+%concat with normalized data
+DSconcat_OT = Ev(7).PSTHz(:,Ishow);
+for i=1:6
+    DSconcat_OT = cat(2,DSconcat_OT,Ev(i).PSTHz(:,Ishow));
+end
+
+for i=1:2 %loop through habit/not habit neurons
+DLSselection=Coord(:,4)==10 & Celltype(:,1)==1 & TRN(:,1)~=0 & RatID(:,2)==i-1;
+DMSselection=Coord(:,4)==20 & Celltype(:,1)==1 & TRN(:,1)~=0 & RatID(:,2)==i-1;
+
+LPwin=[52:280];
+binLP=mean(DSconcat_OT(:,LPwin),2);
+SEQwin=[26:331];
+binSEQ=mean(DSconcat_OT(TRN(:,1)~=0,SEQwin),2);
+TRN(TRN(:,1)~=0,2)=sign(binSEQ);
+
+DLSconcat_OT=DSconcat_OT(DLSselection,:);
+DMSconcat_OT=DSconcat_OT(DMSselection,:);
+
+
+subplot(2,3,3+(i-1)*3)
+
+DLSpsth=nanmean(DLSconcat_OT,1);
+DLSsem=nanste(DLSconcat_OT,1);
+DLSup=DLSpsth+DLSsem;
+DLSdown=DLSpsth-DLSsem;
+hold on;
+patch([time2,time2(end:-1:1)],[DLSup,DLSdown(end:-1:1)],myblue,'EdgeColor','none');alpha(0.5);
+g1=plot(time2,DLSpsth,'Color',myblue,'linewidth',1.5);
+
+DMSpsth=nanmean(DMSconcat_OT,1);
+DMSsem=nanste(DMSconcat_OT,1);
+DMSup=DMSpsth+DMSsem;
+DMSdown=DMSpsth-DMSsem;
+patch([time2,time2(end:-1:1)],[DMSup,DMSdown(end:-1:1)],'r','EdgeColor','none');alpha(0.5);
+g2=plot(time2,DMSpsth,'r','linewidth',1.5);
+
+%plot([0 0], [-2 2.5],'LineStyle',':','Color','k');
+plot([0 357], [0 0],'LineStyle',':','Color','k');
+plot([51 51],[-2 5],'LineStyle',':','Color','k')
+plot([102 102],[-2 5],'LineStyle',':','Color','k')
+plot([153 153],[-2 5],'LineStyle',':','Color','k')
+plot([204 204],[-2 5],'LineStyle',':','Color','k')
+plot([255 255],[-2 5],'LineStyle',':','Color','k')
+plot([306 306],[-2 5],'LineStyle',':','Color','k')
+axis([Xaxis,Yaxis]);
+set(gca,'XTick',[0:25.5:357]);
+set(gca,'xticklabel',Eventnames)
+traces=cat(2,g1,g2);
+labels={'DLS','DMS'};
+ylabel('z-score');
+legend(traces,labels)
+title(titlegraph(i))
+end
+
+DSconcat_OT = Ev(7).PSTHz(:,Ishow);
+for i=1:6
+    DSconcat_OT = cat(2,DSconcat_OT,Ev(i).PSTHz(:,Ishow));
+end
+LPwin=[52:280];
+binLP=mean(DSconcat_OT(:,LPwin),2);
+
+%% heatmap Extended training neurons
+figure(2)
+
+titleplot={'DLS sensitive', 'DMS sensitive', 'DLS insensitive', 'DMS insensitive'};
+for i=1:max(table_OT(:,6))
+    nbplottitle=0;
+    for j=1:2 %loop thru habit profile
+        for k=1:2 %loop thru region
+            nbplottitle=nbplottitle+1;
+            selection=Coord(:,4)==k*10 & RatID(:,2)==j-1 & table_OT(:,6)==i;
+            DSconcat_Class_OT=DSconcat_OT(selection,:);
+            
+            %TMPdls=R.Class(DLSselection,1); TMPdms=R.Class(DMSselection,1);
+            TMPds=binLP(selection,1);
+            TMPds(isnan(TMPds))=0;
+            [~,SORTimgDS]=sort(TMPds);[~,SORTimgDS]=sort(TMPds);
+            DSconcat_Class_OT=DSconcat_Class_OT(SORTimgDS,:);
+            
+            subplot(4,5,1+(i-1)+10*(k-1)+5*(j-1))
+            Clim=[-5 5];
+            
+            [nr,nc] = size(DSconcat_Class_OT);
+            imagesc(time,[1 size(DSconcat_Class_OT,1)],DSconcat_Class_OT,Clim); %colorbar;axis tight;
+            colormap(jet);
+            hold on
+            plot([0.25 0.25],[size(DSconcat_Class_OT,1)+0.5 0.5],'k')
+            plot([0.75 0.75],[size(DSconcat_Class_OT,1)+0.5 0.5],'k')
+            plot([1.25 1.25],[size(DSconcat_Class_OT,1)+0.5 0.5],'k')
+            plot([1.75 1.75],[size(DSconcat_Class_OT,1)+0.5 0.5],'k')
+            plot([2.25 2.25],[size(DSconcat_Class_OT,1)+0.5 0.5],'k')
+            plot([2.75 2.75],[size(DSconcat_Class_OT,1)+0.5 0.5],'k')
+            set(gca,'XTick',[0:0.25:3.25]);
+            set(gca,'xticklabel',Eventnames)
+            ylabel('neurons');
+            title(titleplot{nbplottitle})
+            hold off
+        end
+    end
+end
+
+%% EARLY TRAININIG DATASET - average PSTH
+
+clear Ev 
+load ('Rearlytraining_light.mat');
+load ('Celltype_earlyTraining.mat');
+load('TRN_DT5_earlytrain.mat')
+load('RatID_earlytraining.mat');
+load('table_ACQ.mat')
+
+
+for i=1:max(Celltype(:,1))
+    selectionIN=Ses(i).Coord(:,4)~=0;
+    selectsession=Celltype(:,1)==i;
+
+    Ses(i).Celltype(1:length(selectionIN),1)=0;
+    Ses(i).Celltype(selectionIN,1)=Celltype(selectsession,2); 
+    Ses(i).RatID(selectionIN,1)=RatID_ACQ(selectsession,1);  
+end
+
+
+%analysis restricted on putative MSN. nonTRN included
+%concat with normalized data
+nbplot=0;
+figure(1)
+titlegraph={'1st session - sensitive','last session - sensitive', '1st session - insensitive','last session - insensitive'};
+for l=1:2 %loop thru habit profile
+    for j=1:2 %loop thru 1st last session
+        vars={'TMPdls','TMPdms','binSEQ','binLP','DLSselection','DMSselection','DSconcat_ACQ','DLSconcat_ACQ','DMSconcat_ACQ'};
+        clear(vars{:});
+        nbplot=nbplot+1;
+        
+        if j==1
+            k=4;
+        else
+            k=13;
+        end
+        
+        DSconcat_ACQ = Ses(k).Ev(7).PSTHz(:,Ishow);
+        for i=1:6
+            DSconcat_ACQ = cat(2,DSconcat_ACQ,Ses(k).Ev(i).PSTHz(:,Ishow));
+        end
+        
+        DLSselection=Ses(k).Coord(:,4)==10 & Ses(k).Celltype(:,1)==1 & sesTRN(k).TRN(:,1)~=0 & Ses(k).RatID(:,1)==l-1;
+        DMSselection=Ses(k).Coord(:,4)==20 & Ses(k).Celltype(:,1)==1 & sesTRN(k).TRN(:,1)~=0 & Ses(k).RatID(:,1)==l-1;
+        
+        LPwin=[52:280];
+        binLP=mean(DSconcat_ACQ(:,LPwin),2);
+        SEQwin=[26:331];
+        binSEQ=mean(DSconcat_ACQ(:,SEQwin),2);
+        sesTRN(k).TRN(:,2)=sign(binSEQ);
+        
+        DLSconcat_ACQ=DSconcat_ACQ(DLSselection,:);
+        DMSconcat_ACQ=DSconcat_ACQ(DMSselection,:);
+        
+        
+        
+        %% Average OT PSTH
+        traces=[];
+        subplot(2,3,1+(j-1)+3*(l-1))
+        
+        DLSpsth=nanmean(DLSconcat_ACQ,1);
+        DLSsem=nanste(DLSconcat_ACQ,1);
+        DLSup=DLSpsth+DLSsem;
+        DLSdown=DLSpsth-DLSsem;
+        hold on;
+        patch([time2,time2(end:-1:1)],[DLSup,DLSdown(end:-1:1)],myblue,'EdgeColor','none');alpha(0.5);
+        g1=plot(time2,DLSpsth,'Color',myblue,'linewidth',1.5);
+        
+        DMSpsth=nanmean(DMSconcat_ACQ,1);
+        DMSsem=nanste(DMSconcat_ACQ,1);
+        DMSup=DMSpsth+DMSsem;
+        DMSdown=DMSpsth-DMSsem;
+        patch([time2,time2(end:-1:1)],[DMSup,DMSdown(end:-1:1)],'r','EdgeColor','none');alpha(0.5);
+        g2=plot(time2,DMSpsth,'r','linewidth',1.5);
+        
+        %plot([0 0], [-2 2.5],'LineStyle',':','Color','k');
+        plot([0 357], [0 0],'LineStyle',':','Color','k');
+        plot([51 51],[-2 5],'LineStyle',':','Color','k')
+        plot([102 102],[-2 5],'LineStyle',':','Color','k')
+        plot([153 153],[-2 5],'LineStyle',':','Color','k')
+        plot([204 204],[-2 5],'LineStyle',':','Color','k')
+        plot([255 255],[-2 5],'LineStyle',':','Color','k')
+        plot([306 306],[-2 5],'LineStyle',':','Color','k')
+        axis([Xaxis,Yaxis]);
+        set(gca,'XTick',[0:25.5:357]);
+        set(gca,'xticklabel',Eventnames)
+        traces=cat(2,g1,g2);
+        labels={'DLS','DMS'};
+        ylabel('z-score');
+        legend(traces,labels)
+        title(titlegraph(nbplot))
+        
+    end
+end
+
+%% Early training data set - heatmaps
+
+
+%Xaxis=[0 357];Yaxis=[-1 3];
+Xaxis1=[-0.25 0.25];
+Ishow=find(Tm>=Xaxis1(1) & Tm<=Xaxis1(2));
+time=-0.25:0.25:3.25;
+Eventnames={'LI';'';'LP1';'';'LP2';'';'LP3';'';'LP4';'';'LP5';'';'PE';''};
+
+
+for i=1:max(Celltype(:,1))
+    selectionIN=Ses(i).Coord(:,4)~=0;
+    selectsession=Celltype(:,1)==i;
+
+    Ses(i).Celltype(1:length(selectionIN),1)=0;
+    Ses(i).Celltype(selectionIN,1)=Celltype(selectsession,2); 
+    Ses(i).RatID(selectionIN,1)=RatID_ACQ(selectsession,1); 
+end
+
+
+figure
+
+for i=1:max(table_ACQ(:,5))
+    nbtitle=0;
+    for j=1:2 %loop thru habit profile
+        for k=1:2 %loop thru region
+            nbtitle=nbtitle+1;
+            for l=4:13 %loop thru session
+                DSconcat_ACQ=[];
+                selection=Ses(l).Coord(:,4)==k*10 & Ses(l).Celltype(:,1)==1 & Ses(l).RatID(:,1)==j-1 & Ses(l).Class(:,1)==i;
+                DSconcat_ACQ = Ses(l).Ev(7).PSTHz(selection,Ishow);
+                for m=1:6
+                    DSconcat_ACQ = cat(2,DSconcat_ACQ,Ses(l).Ev(m).PSTHz(selection,Ishow));
+                end
+                MeanDS.img(l-3,:)=nanmean(DSconcat_ACQ,1);
+                if j==1
+                    groupGD(k).region(l-3,i)=sum(selection);
+                else
+                    groupHAB(k).region(l-3,i)=sum(selection);
+                end
+            end
+            
+            subplot(4,5,1+(i-1)+10*(k-1)+5*(j-1))
+            Clim=[-5 5];
+            imagesc(time,[1 size(MeanDS.img,1)],MeanDS.img,Clim); %colorbar;axis tight;
+            colormap(jet);
+            hold on
+            plot([0.25 0.25],[size(MeanDS.img,1)+0.5 0.5],'k')
+            plot([0.75 0.75],[size(MeanDS.img,1)+0.5 0.5],'k')
+            plot([1.25 1.25],[size(MeanDS.img,1)+0.5 0.5],'k')
+            plot([1.75 1.75],[size(MeanDS.img,1)+0.5 0.5],'k')
+            plot([2.25 2.25],[size(MeanDS.img,1)+0.5 0.5],'k')
+            plot([2.75 2.75],[size(MeanDS.img,1)+0.5 0.5],'k')
+            set(gca,'XTick',[0:0.25:3.25]);
+            set(gca,'xticklabel',Eventnames)
+            ylabel('neurons');
+            title(titleplot(nbtitle));
+            hold off
+        end
+    end
+end

scripts/J_generateFigure10_correlLPlatency.m

@@ -0,0 +1,128 @@
+%% generate Figure 10 and Figure 10 supplement 1 - Correlation analysis
+
+%firing analysis post lever insertion and correl with response latency
+%Analysis on all neurons without DLS/DMS dissociation
+%histogram for correlated neurons in each regions, in each class
+
+clear all;clc;
+tic
+%global Dura Baseline Tm Tbase BSIZE Tbin
+SAVE_FLAG=1;
+BSIZE=0.01;
+Dura=[-2 4]; Tm=Dura(1):BSIZE:Dura(2);
+Baseline=[-2 0]; Tbase=Baseline(1):BSIZE:Baseline(2); %used to calculate Z-scores
+Tbin=-1:0.005:1; %window used to determine the optimal binsize
+PStat=0.05;
+prewin=[Dura(1) 0];
+postwin=[0 .5];
+postwin2=[Dura(1):0.05:Dura(2)]; %bounds should match Dura
+Slopebounds=[-1.5:0.25:1.5];
+R2Bounds=[0:0.05:0.5];
+PvalBounds=[0:0.01:0.5];
+Struct=10;
+XI=1;
+RunRegress=1;
+SAVEFIG=1;
+%R=[];R.Ninfo={}; 
+NN=0;
+pre=[]; post=[];
+if ~exist('RAW'), load ('RAWextendedtraining.mat'); RAW=RAW; end
+if ~exist('Ev'), load('Rextendedtraining_light.mat'); end
+load('C');
+
+%path='C:\FRED\GALLO\Exps\2009-present GoNoGo\FINAL\PLX Files GoNoGo\GoNoGo NEX files\RESULTgn14.xls';
+%%
+
+if RunRegress
+  % PREALLOCATING
+    for i=1:length(RAW)
+        EvInd=strmatch('LeverInsertion',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
+        Rind=strmatch('First_LP',RAW(i).Einfo(:,2),'exact');  %finds the index of the event
+        
+        if  sum(EvInd)~=0 && sum(Rind)~=0
+            DS(i)=length(RAW(i).Erast{EvInd});
+            Neur(i)=length(RAW(i).Nrast);
+        end
+    end
+    MaxTrials=max(DS); MaxNeur=sum(Neur); MaxWin=length(postwin);
+    C.Lat=NaN(MaxNeur, MaxTrials); % (neurons, trials)
+    C.pre=NaN(MaxNeur,MaxTrials);  % (neurons, trials)
+    C.post=NaN(MaxNeur,MaxTrials); % (neurons, trials)
+    C.postwin=NaN(MaxNeur,MaxTrials,MaxWin); % (neurons, trials, windows)
+    
+    
+    %
+    for i=1:length(RAW) %loops through sessions
+        EvInd=strmatch('LeverInsertion',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
+        Rind=strmatch('First_LP',RAW(i).Einfo(:,2),'exact');  %finds the index of the event
+        
+        if  sum(EvInd)~=0 && sum(Rind)~=0
+            Dcell=MakePSR04(RAW(i).Erast(Rind),RAW(i).Erast{EvInd},[-1 60],{2,'first'});
+            D=cell2mat(Dcell); %inds=find(~isnan(D));
+            D(isnan(D))=60; %trials with no response are set to 10.25 to allow to use them in the correlation
+            for j= 1:size(RAW(i).Nrast,1) %Number of neurons within sessions
+                NN=NN+1;
+                C.Lat(NN,1:length(RAW(i).Erast{EvInd}))=D;
+                [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{2});% makes trial by trail rasters. PSR1 is a cell(neurons, trials)
+                for m=1:size(RAW(i).Erast{Rind},1) %loops through trials
+                    Resplat=RAW(i).Eint{1,8}(m)-RAW(i).Eint{1,7}(m);
+                    C.pre(NN,m)=sum(PSR2{1,m}<prewin(2) & PSR2{1,m}>prewin(1))/(abs(prewin(2)-prewin(1))); %neurons, trials
+                    C.postLAT1(NN,m)=sum(PSR2{1,m}<postwin(2) & PSR2{1,m}>postwin(1))/(abs(postwin(2)-postwin(1)));
+                end %m loop
+                fprintf('Neuron ID # %d\n',NN);
+            end %j loop
+        end %if
+    end %i loop
+    
+    
+    %% Runs the regression analysis
+    for NN=1:MaxNeur
+        RegX=[C.Lat(NN,:)', ones(size(C.Lat,2),1)];
+        [C.SLOPEpre,C.STATSpre]=corr(C.pre(NN,:)',C.Lat(NN,:)','rows','pairwise','type','Spearman');
+        [C.SLOPEpostLAT1(NN,:),C.STATSpostLAT1(NN,:)]= corr(C.postLAT1(NN,:)',C.Lat(NN,:)','rows','pairwise','type','Spearman');
+        [C.SLOPEpostpre(NN,:),C.STATSpostpre(NN,:)]= corr((C.post(NN,:)-C.pre(NN,:))',C.Lat(NN,:)','rows','pairwise','type','Spearman');
+
+        fprintf('CORREL Neuron ID # %d\n',NN);
+    end
+    
+    % Saving data in C structure and excel file
+    save('C.mat', 'C');
+    Cmat=[];
+%    Cmat=cat(2,C.STATSpre, C.SLOPEpre, C.STATSpost, C.SLOPEpost);
+    %xlswrite(path,Cmat,'Correl', 'F3');
+end
+
+%% --------------Plots distribution of correlation variables SINGLE WINDOW -------------
+% POST WINDOW
+
+XI=1; clear A;
+
+for i=1:2
+    figure(1)
+    sel1=Coord(:,4)==i*10;
+    sel2= sel1 & C.STATSpostLAT1(:,1)<PStat;
+    A(1,:)=cat(2, sum(sel1), sum(sel2));     
+    h=subplot(2,2,i);
+    xmin=floor(min(C.SLOPEpostLAT1(sel1,1)));
+    xmax=ceil(max(C.SLOPEpostLAT1(sel1,1)));
+    step=0.05;
+    xcenters=xmin:step:xmax;
+    Nsel1=hist(C.SLOPEpostLAT1(sel1,1),xcenters);
+    Nsel2=hist(C.SLOPEpostLAT1(sel2,1),xcenters);
+    hold on;
+    bar(xcenters, Nsel1, 'w', 'EdgeColor', 'k', 'BarWidth', 1);alpha(0.3);
+    bar(xcenters, Nsel2, 'k', 'EdgeColor', 'k', 'BarWidth', 1);
+    plot([0 0], [0 30],'LineStyle',':','Color','k');
+    %axis([-10 10 0 20]) %note: their is an outlier @-14 for shell
+    pval1=signrank(C.SLOPEpostLAT1(sel1,1));
+    pval2=signrank(C.SLOPEpostLAT1(sel2,1));
+    MEANsel1=mean(C.SLOPEpostLAT1(sel2,1));
+    plot([MEANsel1 MEANsel1], [0 10],'Color','r');
+    %text(-5.8,9.8,sprintf('p=%d', pval1), 'Parent',h);
+    %text(-5.8,8.8,sprintf('p=%d', pval2), 'Parent',h);    
+    subplot(2,2,i+2)
+    [N,xcenters2]=hist(C.STATSpostLAT1(sel2,1),20);
+    bar(xcenters2, N, 'k', 'EdgeColor', 'k', 'BarWidth', 1); alpha(0.3);
+    axis ([0 ceil(max(xcenters2)*10)/10 0 ceil(max(N))+1]);
+end
+save('C.mat', 'C');

scripts/J_generateFigure10_correlPElatency.m

@@ -0,0 +1,123 @@
+%% generate Figure 10 and Figure 10 supplement 1 - Correlation analysis
+
+%firing analysis post lever retraction and correl with PE latency
+%Analysis on all neurons without DLS/DMS dissociation
+%histogram for correlated neurons in each regions, in each class
+
+
+clear all;clc;close all
+tic
+global Dura Baseline Tm Tbase BSIZE Tbin
+SAVE_FLAG=1;
+BSIZE=0.01;
+Dura=[-2 4]; Tm=Dura(1):BSIZE:Dura(2);
+Baseline=[-2 0]; Tbase=Baseline(1):BSIZE:Baseline(2); %used to calculate Z-scores
+Tbin=-1:0.005:1; %window used to determine the optimal binsize
+PStat=0.05;
+prewin=[Dura(1) 0];
+postwin=[0 .5];
+prewinPE=[-0.5 0];
+postwin2=[Dura(1):0.05:Dura(2)]; %bounds should match Dura
+Slopebounds=[-1.5:0.25:1.5];
+R2Bounds=[0:0.05:0.5];
+PvalBounds=[0:0.01:0.5];
+Struct=10;
+XI=1;
+RunRegress=1;
+SAVEFIG=1;
+%R=[];R.Ninfo={}; 
+NN=0;
+pre=[]; post=[];
+if ~exist('RAW'), load ('RAWextendedtraining.mat'); RAW=RAW; end
+if ~exist('Ev'), load('Rextendedtraining_light.mat'); end
+load('C');
+
+
+
+if RunRegress
+  % PREALLOCATING
+    for i=1:length(RAW)
+        EvInd=strmatch('Last_LP',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
+        Rind=strmatch('FirstPE_PostRew',RAW(i).Einfo(:,2),'exact');  %finds the index of the event
+        
+        if  sum(EvInd)~=0 && sum(Rind)~=0
+            DS(i)=length(RAW(i).Erast{EvInd});
+            Neur(i)=length(RAW(i).Nrast);
+        end
+    end
+    MaxTrials=max(DS); MaxNeur=sum(Neur); MaxWin=length(postwin);
+    C.LatPE=NaN(MaxNeur, MaxTrials); % (neurons, trials)
+    C.prePE=NaN(MaxNeur,MaxTrials);  % (neurons, trials)
+    C.postPE=NaN(MaxNeur,MaxTrials); % (neurons, trials)
+    C.postwinPE=NaN(MaxNeur,MaxTrials,MaxWin); % (neurons, trials, windows)
+    
+    
+    %
+    for i=1:length(RAW) %loops through sessions
+        EvInd=strmatch('Last_LP',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
+        Rind=strmatch('FirstPE_PostRew',RAW(i).Einfo(:,2),'exact');  %finds the index of the event
+        
+        if  sum(EvInd)~=0 && sum(Rind)~=0
+            Dcell=MakePSR04(RAW(i).Erast(Rind),RAW(i).Erast{EvInd},[-1 10],{2,'first'});
+            D=cell2mat(Dcell); %inds=find(~isnan(D));
+            %D(isnan(D))=60; %trials with no response are set to 10.25 to allow to use them in the correlation
+            for j= 1:size(RAW(i).Nrast,1) %Number of neurons within sessions
+                NN=NN+1;
+                C.LatPE(NN,1:length(RAW(i).Erast{EvInd}))=D;
+                [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{2});% ref= last lever press / makes trial by trail rasters. 
+                [PSR3,N3]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{Rind},Dura,{2});% ref= post-reward PE / makes trial by trail rasters.
+                for m=1:size(RAW(i).Erast{Rind},1) %loops through trials
+                    PElatency=RAW(i).Erast{Rind}(m)-RAW(i).Erast{EvInd}(m);
+                    C.prePE(NN,m)=sum(PSR2{1,m}<prewin(2) & PSR2{1,m}>prewin(1))/(abs(prewin(2)-prewin(1))); %neurons, trials
+                    C.postPE1(NN,m)=sum(PSR2{1,m}<postwin(2) & PSR2{1,m}>postwin(1))/(abs(postwin(2)-postwin(1))); %When last LP is ref, postwin is 0-0.5s post LP
+                end %m loop
+                fprintf('Neuron ID # %d\n',NN);
+            end %j loop
+        end %if
+    end %i loop
+    
+    
+    %% Runs the regression analysis
+    for NN=1:MaxNeur
+        RegX=[C.LatPE(NN,:)', ones(size(C.LatPE,2),1)];
+        [C.SLOPEprePE,C.STATSprePE]=corr(C.prePE(NN,:)',C.LatPE(NN,:)','rows','pairwise','type','Spearman');
+        [C.SLOPEpostPE1(NN,:),C.STATSpostPE1(NN,:)]= corr(C.postPE1(NN,:)',C.LatPE(NN,:)','rows','pairwise','type','Spearman');
+        [C.SLOPEpostprePE(NN,:),C.STATSpostprePE(NN,:)]= corr((C.postPE(NN,:)-C.prePE(NN,:))',C.LatPE(NN,:)','rows','pairwise','type','Spearman');
+        fprintf('CORREL Neuron ID # %d\n',NN);
+    end
+    
+    save('C.mat', 'C');
+end
+
+%% --------------Plots distribution of correlation variables SINGLE WINDOW -------------
+% POST WINDOW
+figure(1);
+
+for i=1:2
+    sel1=Coord(:,4)==i*10;
+    sel2= sel1 & C.STATSpostPE1(:,1)<PStat;  
+    h=subplot(2,2,i);
+    xmin=floor(min(C.SLOPEpostPE1(sel1,1)));
+    xmax=ceil(max(C.SLOPEpostPE1(sel1,1)));
+    step=0.05;
+    xcenters=xmin:step:xmax;
+    Nsel1=hist(C.SLOPEpostPE1(sel1,1),xcenters);
+    Nsel2=hist(C.SLOPEpostPE1(sel2,1),xcenters);
+    hold on;
+    bar(xcenters, Nsel1, 'w', 'EdgeColor', 'k', 'BarWidth', 1);alpha(0.3);
+    bar(xcenters, Nsel2, 'k', 'EdgeColor', 'k', 'BarWidth', 1);
+    plot([0 0], [0 30],'LineStyle',':','Color','k');
+    %axis([-10 10 0 20]) %note: their is an outlier @-14 for shell
+    pval1=signrank(C.SLOPEpostPE1(sel1,1));
+    pval2=signrank(C.SLOPEpostPE1(sel2,1));
+    MEANsel1=mean(C.SLOPEpostPE1(sel2,1));
+    plot([MEANsel1 MEANsel1], [0 10],'Color','r');
+    %text(-5.8,9.8,sprintf('p=%d', pval1), 'Parent',h);
+    %text(-5.8,8.8,sprintf('p=%d', pval2), 'Parent',h);
+    subplot(2,2,i+2)
+    [N,xcenters2]=hist(C.STATSpostPE1(sel2,1),20);
+    bar(xcenters2, N, 'k', 'EdgeColor', 'k', 'BarWidth', 1); alpha(0.3);
+    axis ([0 ceil(max(xcenters2)*10)/10 0 ceil(max(N))+1]);
+end
+
+save('C.mat', 'C');

scripts/J_generateFigure10_correlRRvsFR.m

@@ -0,0 +1,154 @@
+%% generate Figure 10 and Figure 10 supplement 1 - Correlation analysis - Correlation analysis - RR vs FR
+
+%firing analysis during LP and correl with response rate
+%Analysis on all neurons without DLS/DMS dissociation
+%histogram for correlated neurons in each regions, in each class
+
+clear all;clc;
+tic
+global Dura Baseline Tm Tbase BSIZE Tbin
+SAVE_FLAG=1;
+BSIZE=0.05;
+Dura=[-25 60]; Tm=-1:BSIZE:60;
+Baseline=[-2 0]; Tbase=Baseline(1):BSIZE:Baseline(2); %used to calculate Z-scores
+Tbin=-1:0.005:1; %window used to determine the optimal binsize
+PStat=0.05;
+prewin=[Dura(1) -20];
+postwin=[0 .5];
+EventWin=[-0.5 0.5];
+postwin2=[Dura(1):0.05:Dura(2)]; %bounds should match Dura
+Slopebounds=[-1.5:0.25:1.5];
+R2Bounds=[0:0.05:0.5];
+PvalBounds=[0:0.01:0.5];
+Struct=[10 20];
+XI=1;
+RunRegress=0;
+SAVEFIG=1;
+%R=[];R.Ninfo={}; 
+NN=0;
+pre=[]; post=[];
+if ~exist('RAW'), load ('RAWextendedtraining.mat'); RAW=RAW; end
+if ~exist('Ev'), load('Rextendedtraining_light.mat'); end
+load('C.mat');
+
+%Smoothing
+Smoothing=1; %0 for raw and 1 for smoothing
+SmoothTYPE='lowess';
+SmoothSPAN=5;
+if Smoothing~=1, SmoothTYPE='NoSmoothing';SmoothSPAN=NaN; end
+
+
+%path='C:\FRED\GALLO\Exps\2009-present GoNoGo\FINAL\PLX Files GoNoGo\GoNoGo NEX files\RESULTgn14.xls';
+%%
+
+if RunRegress
+  % PREALLOCATING
+    for i=1:length(RAW)
+        EvInd=strmatch('First_LP',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
+        Rind=strmatch('ReliableLP_DT',RAW(i).Einfo(:,2),'exact');  %finds the index of the event
+        
+        if  sum(EvInd)~=0 && sum(Rind)~=0
+            DS(i)=length(RAW(i).Erast{EvInd});
+            Neur(i)=length(RAW(i).Nrast);
+        end
+    end
+    MaxTrials=max(DS); MaxNeur=sum(Neur); MaxWin=length(postwin);
+%     C.Lat=NaN(MaxNeur, MaxTrials); % (neurons, trials)
+    C.preRR=NaN(MaxNeur,MaxTrials);  % (neurons, trials)
+    C.SEQ=NaN(MaxNeur,MaxTrials); % (neurons, trials)
+    C.postwinRR=NaN(MaxNeur,MaxTrials,MaxWin); % (neurons, trials, windows)
+    
+    
+    %
+    for i=1:length(RAW) %loops through sessions
+        EvInd=strmatch('First_LP',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
+        Rind=strmatch('ReliableLP_DT',RAW(i).Einfo(:,2),'exact');  %finds the index of the event
+        
+        if  sum(EvInd)~=0 && sum(Rind)~=0            
+            Dcell2=MakePSR05(RAW(i).Erast(Rind),RAW(i).Erast{EvInd},[-1 60],{2}); % makes trial by trail 
+                    
+            for j= 1:size(RAW(i).Nrast,1) %Number of neurons within sessions
+                NN=NN+1;
+
+                [PSRNeur,Nneur]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},[-1 60],{1});% collapsed raster
+                [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{2});% makes trial by trail rasters. PSR1 is a cell(neurons, trials)
+                
+                [PTHneur,BW1,~]=MakePTH07(PSRNeur,repmat(N2, size(RAW(i).Nrast{j},1),1),{2,0,BSIZE});
+                PTHneur=smooth(PTHneur,SmoothSPAN,SmoothTYPE)';               
+                C.PSTHraw(NN,1:length(Tm))=PTHneur;
+                Search=find(Tm>=EventWin(1) & Tm<=EventWin(2));
+                [C.MaxVal(NN,1),MaxInd]=max(C.PSTHraw(NN,Search));
+                C.MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
+                
+                for m=1:size(RAW(i).Erast{EvInd},1) %loops through trials
+                    SequenceDur=RAW(i).Eint{1,6}(m)-RAW(i).Eint{1,5}(m);
+                    C.RR(NN,m)=length(Dcell2{1,m})/SequenceDur;
+                    C.preRR(NN,m)=sum(PSR2{1,m}<prewin(2) & PSR2{1,m}>prewin(1))/(abs(prewin(2)-prewin(1))); %neurons, trials
+                    C.SEQ(NN,m)=sum(PSR2{1,m}<SequenceDur & PSR2{1,m}>0)/SequenceDur;
+                end %m loop
+                fprintf('Neuron ID # %d\n',NN);
+            end %j loop
+        end %if
+    end %i loop
+    
+ 
+
+    
+    %% Runs the regression analysis
+    for NN=1:MaxNeur
+        RegX=[C.RR(NN,:)', ones(size(C.RR,2),1)];
+        [C.SLOPEpreRR,C.STATSpreRR]=corr(C.preRR(NN,:)',C.RR(NN,:)','rows','pairwise','type','Spearman');
+        [C.SLOPEseq(NN,:),C.STATSseq(NN,:)]= corr(C.SEQ(NN,:)',C.RR(NN,:)','rows','pairwise','type','Spearman');
+        [C.SLOPEseqpre(NN,:),C.STATSseqpre(NN,:)]= corr((C.SEQ(NN,:)-C.preRR(NN,:))',C.RR(NN,:)','rows','pairwise','type','Spearman');
+        fprintf('CORREL Neuron ID # %d\n',NN);
+    end
+    
+    C.SelectionRR(:,1)=Coord(:,4)==10 & C.STATSseq(:,1)<PStat;
+    C.SelectionRR(:,2)=Coord(:,4)==20 & C.STATSseq(:,1)<PStat;
+    % Saving data in C structure and excel file
+    save('C.mat', 'C');
+
+end
+
+
+%% --------------Plots distribution of correlation variables SINGLE WINDOW -------------
+% POST WINDOW
+figure;
+Struct=[10 20];
+clear A;
+
+for i=1:2 % loop through structure
+    clear sel1;clear sel2;    
+    sel1=Coord(:,4)==i*10; %(R.Class(:,1)==-4); %& R.DSNSStat<PStat & (R.Ev(1).Meanz>R.Ev(2).Meanz);
+    sel2= sel1 & C.STATSseq(:,1)<PStat;
+    
+    
+    A(i,:)=cat(2, sum(Coord(:,4)==i*10), sum(sel1), sum(sel2));
+        
+    h=subplot(2,2,i);
+    xmin=floor(min(C.SLOPEseq(sel1,1)));
+    xmax=ceil(max(C.SLOPEseq(sel1,1)));
+    step=0.05;
+    xcenters=xmin:step:xmax;
+    Nsel1=hist(C.SLOPEseq(sel1,1),xcenters);
+    Nsel2=hist(C.SLOPEseq(sel2,1),xcenters);
+    hold on;
+    bar(xcenters, Nsel1, 'w', 'EdgeColor', 'k', 'BarWidth', 1);alpha(0.3);
+    bar(xcenters, Nsel2, 'k', 'EdgeColor', 'k', 'BarWidth', 1);
+    plot([0 0], [0 30],'LineStyle',':','Color','k');
+    pval1=signrank(C.SLOPEseq(sel1,1));
+    pval2=signrank(C.SLOPEseq(sel2,1));
+    MEANsel1=mean(C.SLOPEseq(sel2,1));
+    plot([MEANsel1 MEANsel1], [0 10],'Color','r');
+
+        
+    subplot(2,2,i+2)
+    [N,xcenters2]=hist(C.STATSseq(sel2,1),20);
+    bar(xcenters2, N, 'k', 'EdgeColor', 'k', 'BarWidth', 1); alpha(0.3);
+    axis ([0 ceil(max(xcenters2)*10)/10 0 ceil(max(N))+1]);
+end
+
+A(:,4)=A(:,3)./A(:,2)*100;
+
+
+save('C.mat', 'C');

scripts/J_generateFigure10_s2_shuffle_LatLP1.m

@@ -0,0 +1,192 @@
+%% generate Figure 10 figure supplement 2
+
+clear all;clc;
+tic
+global Dura Baseline Tm Tbase BSIZE Tbin
+SAVE_FLAG=1;
+BSIZE=0.01;
+Dura=[-2 4]; Tm=Dura(1):BSIZE:Dura(2);
+Baseline=[-2 0]; Tbase=Baseline(1):BSIZE:Baseline(2); %used to calculate Z-scores
+Tbin=-1:0.005:1; %window used to determine the optimal binsize
+PStat=0.05;
+prewin=[Dura(1) 0];
+postwin=[0 .5];
+postwin2=[Dura(1):0.05:Dura(2)]; %bounds should match Dura
+Slopebounds=[-.5:.05:.5];
+R2Bounds=[0:0.05:0.5];
+PvalBounds=[0:0.01:0.5];
+XI=[1 -1];
+RunRegress=1;
+SAVEFIG=1;
+%R=[];R.Ninfo={}; 
+NN=0;
+pre=[]; post=[];
+if ~exist('RAW'), load ('RAWextendedtraining.mat'); RAW=RAW; end
+if ~exist('R'), load('Rextendedtraining.mat'); R=R; end
+if RunRegress==0, load('C'); end 
+load('Celltype_extendedTraining.mat')
+
+%%
+
+   for i=1:length(RAW)
+        EvInd=strmatch('LeverInsertion',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
+        Rind=strmatch('First_LP',RAW(i).Einfo(:,2),'exact');  %finds the index of the event
+        
+        if  sum(EvInd)~=0 && sum(Rind)~=0
+            DS(i)=length(RAW(i).Erast{EvInd});
+            Neur(i)=length(RAW(i).Nrast);
+        end
+   end
+   
+   MaxTrials=max(DS); MaxNeur=sum(Neur); MaxWin=length(postwin);
+
+
+if RunRegress
+  % PREALLOCATING
+
+
+    Cshuffle.Lat=NaN(MaxNeur, MaxTrials); % (neurons, trials)
+    Cshuffle.pre=NaN(MaxNeur,MaxTrials);  % (neurons, trials)
+    Cshuffle.post=NaN(MaxNeur,MaxTrials); % (neurons, trials)
+    Cshuffle.postwin=NaN(MaxNeur,MaxTrials,MaxWin); % (neurons, trials, windows)
+    
+    for i=1:length(RAW) %loops through sessions
+        EvInd=strmatch('LeverInsertion',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
+        Rind=strmatch('First_LP',RAW(i).Einfo(:,2),'exact');  %finds the index of the event
+        
+        if  sum(EvInd)~=0 && sum(Rind)~=0
+            Dcell=MakePSR04(RAW(i).Erast(Rind),RAW(i).Erast{EvInd},[-1 60],{2,'first'});
+            D=cell2mat(Dcell); %inds=find(~isnan(D));
+            D(isnan(D))=60; %trials with no response are set to 60 to allow to use them in the correlation
+            for j= 1:size(RAW(i).Nrast,1) %Number of neurons within sessions
+                NN=NN+1;
+                Cshuffle.Lat(NN,1:length(RAW(i).Erast{EvInd}))=D;
+                [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{2});% makes trial by trail rasters. PSR1 is a cell(neurons, trials)
+                for m=1:size(RAW(i).Erast{EvInd},1) %loops through trials
+                    Cshuffle.pre(NN,m)=sum(PSR2{1,m}<prewin(2) & PSR2{1,m}>prewin(1))/(abs(prewin(2)-prewin(1))); %neurons, trials
+                    Cshuffle.post(NN,m)=sum(PSR2{1,m}<postwin(2) & PSR2{1,m}>postwin(1))/(abs(postwin(2)-postwin(1)));
+                    for k=1:(length(postwin2)-1) %loops through time windows.
+                        Cshuffle.postwin(NN,m,k)=sum(PSR2{1,m}<postwin2(k+1) & PSR2{1,m}>postwin2(k))/(abs(postwin2(k+1)-postwin2(k)));
+                    end %k loop
+                end %m loop
+                fprintf('Neuron ID # %d\n',NN);
+            end %j loop
+        end %if
+    end %i loop
+    
+    
+    
+    %% Runs the regression analysis on XX number of iterations of shuffled data
+    
+XX=1000;
+
+for n=1:XX
+    for NN=1:MaxNeur
+        Cshuffle.LatShuffled(NN,:)=shuffle(Cshuffle.Lat(NN,:));
+        RegX=[Cshuffle.LatShuffled(NN,:)', ones(size(Cshuffle.LatShuffled,2),1)];
+         [Cshuffle.SLOPEpostshuffled(NN,n),Cshuffle.STATSpostshuffled(NN,n)]= corr(Cshuffle.post(NN,:)',Cshuffle.LatShuffled(NN,:)','rows','pairwise','type','Spearman');
+    end
+    fprintf('Shuffled Iteration # %d\n',n);
+    % Saving data in C structure and excel file
+    save('Cshuffle.mat', 'Cshuffle');
+end
+end
+
+%% Calculates number of significant correlations and the mean correlation coefficient for each 
+
+selection1=R.Structure~=0;
+
+for q=1:length(Cshuffle.SLOPEpostshuffled(1,:))
+    Cshuffle.SLOPEmeanshuffled(q)=nanmean(Cshuffle.SLOPEpostshuffled(selection1,q));
+    Cshuffle.SIGshuffled(q)=sum(Cshuffle.STATSpostshuffled(selection1,q)<0.05);
+    Cshuffle.NEGSIGshuffled(q)=sum(Cshuffle.STATSpostshuffled(selection1,q)<0.05 & Cshuffle.SLOPEpostshuffled(selection1,q)<0);
+end
+
+
+
+%% Calculates correlation coefficients and significant values for real data
+
+
+    for NN=1:MaxNeur
+        RegX=[Cshuffle.Lat(NN,:)', ones(size(Cshuffle.Lat,2),1)];
+        [Cshuffle.SLOPEpre,Cshuffle.STATSpre]=corr(Cshuffle.pre(NN,:)',Cshuffle.Lat(NN,:)','rows','pairwise','type','Spearman');
+        [Cshuffle.SLOPEpost(NN,:),Cshuffle.STATSpost(NN,:)]= corr(Cshuffle.post(NN,:)',Cshuffle.Lat(NN,:)','rows','pairwise','type','Spearman');
+        [Cshuffle.SLOPEpostpre(NN,:),Cshuffle.STATSpostpre(NN,:)]= corr((Cshuffle.post(NN,:)-Cshuffle.pre(NN,:))',Cshuffle.Lat(NN,:)','rows','pairwise','type','Spearman'); 
+        fprintf('CORREL Neuron ID # %d\n',NN);
+    end
+    
+    Cshuffle.SelectionLAT(:,1)=R.Structure==10 & R.TRN(:,1)~=0 & Cshuffle.STATSpost(:,1)<PStat;
+    Cshuffle.SelectionLAT(:,2)=R.Structure==20 & R.TRN(:,1)~=0 & Cshuffle.STATSpost(:,1)<PStat;
+    % Saving data in C structure and excel file
+    save('Cshuffle.mat', 'Cshuffle');
+ 
+%% Figure of distributions of means correlation coefficients and number of significant neurons
+selection1=R.Structure==10 & Celltype(:,1)==1;
+selection2=R.Structure==20 & Celltype(:,1)==1;
+figure;
+% Plot of distribution of mean correlation coefficients for all of the
+% shuffled data
+subplot(2,2,3);
+histogram(Cshuffle.SLOPEmeanshuffled,100);
+hold on;
+MEANsel1=nanmean(Cshuffle.SLOPEpost(selection1,1));
+MEANsel2=nanmean(Cshuffle.SLOPEpost(selection2,1));
+plot([MEANsel1 MEANsel1], [0 50],'Color','b');
+plot([MEANsel2 MEANsel2], [0 50],'Color','r');
+xlabel('Mean correlation coefficient') % x-axis label
+ylabel('# of iterations')% y-axis label
+%Plot of distribution of number of significantly correlated neurons across
+%all the shuffled data sets
+subplot(2,2,4);
+histogram((Cshuffle.SIGshuffled*100/sum(Celltype(:,1)==1)),40);
+xlabel('% of significantly correlated unit') % x-axis label
+ylabel('# of iterations')% y-axis label
+hold on;
+
+SigNUM1=100*sum(Cshuffle.STATSpost(selection1,1)<0.05)/(sum(selection1));
+SigNUM2=100*sum(Cshuffle.STATSpost(selection2,1)<0.05)/(sum(selection2));
+plot([SigNUM1 SigNUM1], [0 100],'Color','b');
+plot([SigNUM2 SigNUM2], [0 100],'Color','r');
+propSigSh=Cshuffle.SIGshuffled*100/sum(Celltype(:,1)==1);
+OTpval(1,1)=(sum(propSigSh>SigNUM1)+1)/(length(propSigSh)+1);
+OTpval(2,1)=(sum(propSigSh>SigNUM2)+1)/(length(propSigSh)+1);
+% Plot of correlation coefficients from example shuffled data set
+subplot(2,2,1);
+sel1=R.Structure~=0;
+sel2=sel1 & Cshuffle.STATSpostshuffled(:,1)<PStat;
+xmin=floor(min(Cshuffle.SLOPEpostshuffled(sel1,1)));
+xmax=ceil(max(Cshuffle.SLOPEpostshuffled(sel1,1)));
+step=0.05;
+xcenters=xmin:step:xmax;
+Nsel1=hist(Cshuffle.SLOPEpostshuffled(sel1,1),xcenters);
+Nsel2=hist(Cshuffle.SLOPEpostshuffled(sel2,1),xcenters);
+hold on;
+bar(xcenters, Nsel1, 'w', 'EdgeColor', 'k', 'BarWidth', 1);alpha(0.3);
+bar(xcenters, Nsel2, 'k', 'EdgeColor', 'k', 'BarWidth', 1);
+plot([0 0], [0 30],'LineStyle',':','Color','k');
+MEANsel1=mean(Cshuffle.SLOPEpostshuffled(sel1,1));
+plot([MEANsel1 MEANsel1], [0 10],'Color','r');
+hold on;
+xlabel('% correlation coefficient shuffle data') % x-axis label
+ylabel('# of neurons')% y-axis label
+
+% %Plot of correlation coefficients with significance marked based on
+% %bootstrapping
+subplot(2,2,2);
+sel1=R.Structure~=0;
+sel2= sel1 & Cshuffle.STATSpost(:,1)<PStat;
+xmin=floor(min(Cshuffle.SLOPEpost(sel1,1)));
+xmax=ceil(max(Cshuffle.SLOPEpost(sel1,1)));
+step=0.05;
+xcenters=xmin:step:xmax;
+Nsel1=hist(Cshuffle.SLOPEpost(sel1,1),xcenters);
+Nsel2=hist(Cshuffle.SLOPEpost(sel2,1),xcenters);
+hold on;
+bar(xcenters, Nsel1, 'w', 'EdgeColor', 'k', 'BarWidth', 1);alpha(0.3);
+bar(xcenters, Nsel2, 'k', 'EdgeColor', 'k', 'BarWidth', 1);
+plot([0 0], [0 30],'LineStyle',':','Color','k');
+MEANsel1=mean(Cshuffle.SLOPEpost(sel2,1));
+plot([MEANsel1 MEANsel1], [0 10],'Color','r');
+xlabel('% correlation coefficient real data') % x-axis label
+ylabel('# of neurons')% y-axis label
+

scripts/J_generateFigure10_s3_shuffle_RRFR.m

@@ -0,0 +1,194 @@
+%% generate Figure 10 figure supplement 3
+
+clear all;clc;
+tic
+global Dura Baseline Tm Tbase BSIZE Tbin
+SAVE_FLAG=1;
+BSIZE=0.05;
+Dura=[-25 60]; Tm=-1:BSIZE:60;
+Baseline=[-2 0]; Tbase=Baseline(1):BSIZE:Baseline(2); %used to calculate Z-scores
+Tbin=-1:0.005:1; %window used to determine the optimal binsize
+PStat=0.05;
+prewin=[Dura(1) -20];
+postwin=[0 .5];
+EventWin=[-0.5 0.5];
+postwin2=[Dura(1):0.05:Dura(2)]; %bounds should match Dura
+Slopebounds=[-1.5:0.25:1.5];
+R2Bounds=[0:0.05:0.5];
+PvalBounds=[0:0.01:0.5];
+Struct=[10 20];
+XI=1;
+RunRegress=0;
+SAVEFIG=1;
+%R=[];R.Ninfo={}; 
+NN=0;
+pre=[]; post=[];
+if ~exist('RAW'), load ('RAWextendedtraining.mat'); RAW=RAW; end
+if ~exist('Ev'), load('Rextendedtraining_light.mat'); end
+load('C.mat');
+load('Celltype_extendedTraining.mat')
+% 1 for inclusion of omission in analysis, max time=60s / 0 for exclusion omission
+
+%Smoothing
+Smoothing=1; %0 for raw and 1 for smoothing
+SmoothTYPE='lowess';
+SmoothSPAN=5;
+if Smoothing~=1, SmoothTYPE='NoSmoothing';SmoothSPAN=NaN; end
+
+
+%%
+
+ % PREALLOCATING
+    for i=1:length(RAW)
+        EvInd=strmatch('First_LP',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
+        Rind=strmatch('ReliableLP_DT',RAW(i).Einfo(:,2),'exact');  %finds the index of the event
+        
+        if  sum(EvInd)~=0 && sum(Rind)~=0
+            DS(i)=length(RAW(i).Erast{EvInd});
+            Neur(i)=length(RAW(i).Nrast);
+        end
+    end
+    MaxTrials=max(DS); MaxNeur=sum(Neur); MaxWin=length(postwin);
+    Cshuffle.preRRRR=NaN(MaxNeur,MaxTrials);  % (neurons, trials)
+    Cshuffle.SEQ=NaN(MaxNeur,MaxTrials); % (neurons, trials)
+    Cshuffle.postwinRR=NaN(MaxNeur,MaxTrials,MaxWin); % (neurons, trials, windows)
+    
+    for i=1:length(RAW) %loops through sessions
+        EvInd=strmatch('First_LP',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
+        Rind=strmatch('ReliableLP_DT',RAW(i).Einfo(:,2),'exact');  %finds the index of the event
+        
+        if  sum(EvInd)~=0 && sum(Rind)~=0            
+            Dcell2=MakePSR05(RAW(i).Erast(Rind),RAW(i).Erast{EvInd},[-1 60],{2}); % makes trial by trail 
+                    
+            for j= 1:size(RAW(i).Nrast,1) %Number of neurons within sessions
+                NN=NN+1;
+
+                [PSRNeur,Nneur]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},[-1 60],{1});% collapsed raster
+                [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{2});% makes trial by trail rasters. PSR1 is a cell(neurons, trials)
+                
+                [PTHneur,BW1,~]=MakePTH07(PSRNeur,repmat(N2, size(RAW(i).Nrast{j},1),1),{2,0,BSIZE});
+                PTHneur=smooth(PTHneur,SmoothSPAN,SmoothTYPE)';               
+                Cshuffle.PSTHraw(NN,1:length(Tm))=PTHneur;
+                Search=find(Tm>=EventWin(1) & Tm<=EventWin(2));
+                [Cshuffle.MaxVal(NN,1),MaxInd]=max(Cshuffle.PSTHraw(NN,Search));
+                Cshuffle.MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
+                
+                for m=1:size(RAW(i).Erast{EvInd},1) %loops through trials
+                    SequenceDur=RAW(i).Eint{1,6}(m)-RAW(i).Eint{1,5}(m);
+                    Cshuffle.RR(NN,m)=length(Dcell2{1,m})/SequenceDur;
+                    Cshuffle.preRR(NN,m)=sum(PSR2{1,m}<prewin(2) & PSR2{1,m}>prewin(1))/(abs(prewin(2)-prewin(1))); %neurons, trials
+                    Cshuffle.SEQ(NN,m)=sum(PSR2{1,m}<SequenceDur & PSR2{1,m}>0)/SequenceDur;
+                end %m loop
+                fprintf('Neuron ID # %d\n',NN);
+            end %j loop
+        end %if
+    end %i loop
+    
+    
+    
+    
+    %% Runs the regression analysis on XX number of iterations of shuffled data
+    
+XX=1000;
+
+for n=1:XX
+    for NN=1:MaxNeur
+        Cshuffle.SEQShuffled(NN,:)=shuffle(Cshuffle.SEQ(NN,:));
+        RegX=[Cshuffle.SEQShuffled(NN,:)', ones(size(Cshuffle.SEQShuffled,2),1)];
+        [Cshuffle.SLOPEpostshuffled(NN,n),Cshuffle.STATSpostshuffled(NN,n)]= corr(Cshuffle.RR(NN,:)',Cshuffle.SEQShuffled(NN,:)','rows','pairwise','type','Spearman');
+    end
+    fprintf('Shuffled Iteration # %d\n',n);
+    save('Cshuffle.mat', 'Cshuffle');
+end
+
+
+%% Calculates number of significant correlations and the mean correlation coefficient for each 
+
+selection1=Coord(:,4)~=0;
+
+for q=1:length(Cshuffle.SLOPEpostshuffled(1,:))
+    Cshuffle.SLOPEmeanshuffled(q)=nanmean(Cshuffle.SLOPEpostshuffled(selection1,q));
+    Cshuffle.SIGshuffled(q)=sum(Cshuffle.STATSpostshuffled(selection1,q)<0.05);
+    Cshuffle.NEGSIGshuffled(q)=sum(Cshuffle.STATSpostshuffled(selection1,q)<0.05 & Cshuffle.SLOPEpostshuffled(selection1,q)<0);
+end
+
+
+
+%% Calculates correlation coefficients and significant values for real data
+
+
+    for NN=1:MaxNeur
+        RegX=[Cshuffle.SEQ(NN,:)', ones(size(Cshuffle.SEQ,2),1)];
+        [Cshuffle.SLOPEpre,Cshuffle.STATSpre]=corr(Cshuffle.preRR(NN,:)',Cshuffle.SEQ(NN,:)','rows','pairwise','type','Spearman');
+        [Cshuffle.SLOPEpost(NN,:),Cshuffle.STATSpost(NN,:)]= corr(Cshuffle.RR(NN,:)',Cshuffle.SEQ(NN,:)','rows','pairwise','type','Spearman');
+        [Cshuffle.SLOPEpostpre(NN,:),Cshuffle.STATSpostpre(NN,:)]= corr((Cshuffle.RR(NN,:)-Cshuffle.preRR(NN,:))',Cshuffle.SEQ(NN,:)','rows','pairwise','type','Spearman');
+        fprintf('CORREL Neuron ID # %d\n',NN);
+    end 
+    save('Cshuffle.mat', 'Cshuffle');
+ 
+%% Figure of distributions of means correlation coefficients and number of significant neurons
+selection1=Coord(:,4)==10 & Celltype(:,1)==1;
+selection2=Coord(:,4)==20 & Celltype(:,1)==1;
+figure;
+% Plot of distribution of mean correlation coefficients for all of the
+% shuffled data
+subplot(2,2,3);
+histogram(Cshuffle.SLOPEmeanshuffled,100);
+hold on;
+MEANsel1=nanmean(Cshuffle.SLOPEpost(selection1,1));
+MEANsel2=nanmean(Cshuffle.SLOPEpost(selection2,1));
+plot([MEANsel1 MEANsel1], [0 50],'Color','b');
+plot([MEANsel2 MEANsel2], [0 50],'Color','r');
+xlabel('Mean correlation coefficient') % x-axis label
+ylabel('# of iterations')% y-axis label
+%Plot of distribution of number of significantly correlated neurons across
+%all the shuffled data sets
+subplot(2,2,4);
+histogram((Cshuffle.SIGshuffled*100/849),40);
+xlabel('% of significantly correlated unit') % x-axis label
+ylabel('# of iterations')% y-axis label
+hold on;
+
+SigNUM1=100*sum(Cshuffle.STATSpost(selection1,1)<0.05)/(sum(selection1));
+SigNUM2=100*sum(Cshuffle.STATSpost(selection2,1)<0.05)/(sum(selection2));
+plot([SigNUM1 SigNUM1], [0 100],'Color','b');
+plot([SigNUM2 SigNUM2], [0 100],'Color','r');
+% Plot of correlation coefficients from example shuffled data set
+subplot(2,2,1);
+sel1=Coord(:,4)~=0;
+sel2=sel1 & Cshuffle.STATSpostshuffled(:,1)<PStat;
+xmin=floor(min(Cshuffle.SLOPEpostshuffled(sel1,1)));
+xmax=ceil(max(Cshuffle.SLOPEpostshuffled(sel1,1)));
+step=0.05;
+xcenters=xmin:step:xmax;
+Nsel1=hist(Cshuffle.SLOPEpostshuffled(sel1,1),xcenters);
+Nsel2=hist(Cshuffle.SLOPEpostshuffled(sel2,1),xcenters);
+hold on;
+bar(xcenters, Nsel1, 'w', 'EdgeColor', 'k', 'BarWidth', 1);alpha(0.3);
+bar(xcenters, Nsel2, 'k', 'EdgeColor', 'k', 'BarWidth', 1);
+plot([0 0], [0 30],'LineStyle',':','Color','k');
+MEANsel1=mean(Cshuffle.SLOPEpostshuffled(sel1,1));
+plot([MEANsel1 MEANsel1], [0 10],'Color','r');
+xlabel('% correlation coefficient shuffle data') % x-axis label
+ylabel('# of neurons')% y-axis label
+hold on;
+
+% %Plot of correlation coefficients with significance marked based on
+% %bootstrapping
+subplot(2,2,2);
+sel1=Coord(:,4)~=0;
+sel2= sel1 & Cshuffle.STATSpost(:,1)<PStat;
+xmin=floor(min(Cshuffle.SLOPEpost(sel1,1)));
+xmax=ceil(max(Cshuffle.SLOPEpost(sel1,1)));
+step=0.05;
+xcenters=xmin:step:xmax;
+Nsel1=hist(Cshuffle.SLOPEpost(sel1,1),xcenters);
+Nsel2=hist(Cshuffle.SLOPEpost(sel2,1),xcenters);
+hold on;
+bar(xcenters, Nsel1, 'w', 'EdgeColor', 'k', 'BarWidth', 1);alpha(0.3);
+bar(xcenters, Nsel2, 'k', 'EdgeColor', 'k', 'BarWidth', 1);
+plot([0 0], [0 30],'LineStyle',':','Color','k');
+MEANsel1=mean(Cshuffle.SLOPEpost(sel2,1));
+plot([MEANsel1 MEANsel1], [0 10],'Color','r');
+xlabel('% correlation coefficient real data') % x-axis label
+ylabel('# of neurons')% y-axis label

scripts/J_generateFigure10_s4_shuffle_LatPE.m

@@ -0,0 +1,195 @@
+%% generate Figure 10 figure supplement 4
+
+clear all;clc;
+tic
+global Dura Baseline Tm Tbase BSIZE Tbin
+SAVE_FLAG=1;
+BSIZE=0.01;
+Dura=[-2 4]; Tm=Dura(1):BSIZE:Dura(2);
+Baseline=[-2 0]; Tbase=Baseline(1):BSIZE:Baseline(2); %used to calculate Z-scores
+Tbin=-1:0.005:1; %window used to determine the optimal binsize
+PStat=0.05;
+prewin=[Dura(1) 0];
+postwin=[0 .5];
+postwin2=[Dura(1):0.05:Dura(2)]; %bounds should match Dura
+Slopebounds=[-.5:.05:.5];
+R2Bounds=[0:0.05:0.5];
+PvalBounds=[0:0.01:0.5];
+XI=[1 -1];
+RunRegress=1;
+SAVEFIG=1;
+%R=[];R.Ninfo={}; 
+NN=0;
+pre=[]; post=[];
+if ~exist('RAW'), load ('RAWextendedtraining.mat'); RAW=RAW; end
+if ~exist('R'), load('Rextendedtraining.mat'); R=R; end
+if RunRegress==0, load('C'); end 
+load('Celltype_extendedTraining.mat')
+
+%path='C:\FRED\GALLO\Exps\2009-present GoNoGo\FINAL\PLX Files GoNoGo\GoNoGo NEX files\RESULTgn14.xls';
+%%
+for i=1:length(RAW)
+    EvInd=strmatch('Last_LP',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
+    Rind=strmatch('FirstPE_PostRew',RAW(i).Einfo(:,2),'exact');  %finds the index of the event
+    
+    if  sum(EvInd)~=0 && sum(Rind)~=0
+        DS(i)=length(RAW(i).Erast{EvInd});
+        Neur(i)=length(RAW(i).Nrast);
+    end
+end
+MaxTrials=max(DS); MaxNeur=sum(Neur); MaxWin=length(postwin);
+
+
+if RunRegress
+  % PREALLOCATING
+
+  C.LatPE=NaN(MaxNeur, MaxTrials); % (neurons, trials)
+  C.prePE=NaN(MaxNeur,MaxTrials);  % (neurons, trials)
+  C.postPE=NaN(MaxNeur,MaxTrials); % (neurons, trials)
+  C.postwinPE=NaN(MaxNeur,MaxTrials,MaxWin); % (neurons, trials, windows)
+  
+  Cshuffle.LatPE=NaN(MaxNeur, MaxTrials); % (neurons, trials)
+  Cshuffle.prePE=NaN(MaxNeur,MaxTrials);  % (neurons, trials)
+  Cshuffle.postPE=NaN(MaxNeur,MaxTrials); % (neurons, trials)
+  Cshuffle.postwinPE=NaN(MaxNeur,MaxTrials,MaxWin); % (neurons, trials, windows)
+  
+  for i=1:length(RAW) %loops through sessions
+      EvInd=strmatch('Last_LP',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
+      Rind=strmatch('FirstPE_PostRew',RAW(i).Einfo(:,2),'exact');  %finds the index of the event
+      
+      if  sum(EvInd)~=0 && sum(Rind)~=0
+          Dcell=MakePSR04(RAW(i).Erast(Rind),RAW(i).Erast{EvInd},[-1 10],{2,'first'});
+          D=cell2mat(Dcell); %inds=find(~isnan(D));
+          D(isnan(D))=60; %trials with no response are set to 10.25 to allow to use them in the correlation
+          for j= 1:size(RAW(i).Nrast,1) %Number of neurons within sessions
+              NN=NN+1;
+              Cshuffle.LatPE(NN,1:length(RAW(i).Erast{EvInd}))=D;
+              [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{2});% ref= last lever press / makes trial by trail rasters.
+              for m=1:size(RAW(i).Erast{Rind},1) %loops through trials
+                  PElatency=RAW(i).Erast{Rind}(m)-RAW(i).Erast{EvInd}(m);
+                  Cshuffle.prePE(NN,m)=sum(PSR2{1,m}<prewin(2) & PSR2{1,m}>prewin(1))/(abs(prewin(2)-prewin(1))); %neurons, trials
+                  Cshuffle.postPE1(NN,m)=sum(PSR2{1,m}<postwin(2) & PSR2{1,m}>postwin(1))/(abs(postwin(2)-postwin(1))); %When last LP is ref, postwin is 0-0.5s post LP
+              end %m loop
+              fprintf('Neuron ID # %d\n',NN);
+          end %j loop
+      end %if
+  end %i loop
+    
+    %% Runs the regression analysis on XX number of iterations of shuffled data
+    
+XX=1000;
+
+for n=1:XX
+    for NN=1:MaxNeur
+        Cshuffle.LatShuffled(NN,:)=shuffle(Cshuffle.LatPE(NN,:));
+        RegX=[Cshuffle.LatShuffled(NN,:)', ones(size(Cshuffle.LatShuffled,2),1)];
+        [Cshuffle.SLOPEpostshuffled(NN,n),Cshuffle.STATSpostshuffled(NN,n)]= corr(Cshuffle.postPE1(NN,:)',Cshuffle.LatShuffled(NN,:)','rows','pairwise','type','Spearman');
+    end
+    fprintf('Shuffled Iteration # %d\n',n);
+    % Saving data in C structure and excel file
+    save('Cshuffle.mat', 'Cshuffle');
+end
+end
+
+%% Calculates number of significant correlations and the mean correlation coefficient for each 
+
+selection1=R.Structure~=0;
+
+for q=1:length(Cshuffle.SLOPEpostshuffled(1,:))
+    Cshuffle.SLOPEmeanshuffled(q)=nanmean(Cshuffle.SLOPEpostshuffled(selection1,q));
+    Cshuffle.SIGshuffled(q)=sum(Cshuffle.STATSpostshuffled(selection1,q)<0.05);
+    Cshuffle.NEGSIGshuffled(q)=sum(Cshuffle.STATSpostshuffled(selection1,q)<0.05 & Cshuffle.SLOPEpostshuffled(selection1,q)<0);
+end
+
+
+
+%% Calculates correlation coefficients and significant values for real data
+
+
+    for NN=1:MaxNeur
+        RegX=[Cshuffle.LatPE(NN,:)', ones(size(Cshuffle.LatPE,2),1)];
+        [Cshuffle.SLOPEpre,Cshuffle.STATSpre]=corr(Cshuffle.prePE(NN,:)',Cshuffle.LatPE(NN,:)','rows','pairwise','type','Spearman');
+        [Cshuffle.SLOPEpost(NN,:),Cshuffle.STATSpost(NN,:)]= corr(Cshuffle.postPE1(NN,:)',Cshuffle.LatPE(NN,:)','rows','pairwise','type','Spearman');  
+        [Cshuffle.SLOPEpostpre(NN,:),Cshuffle.STATSpostpre(NN,:)]= corr((Cshuffle.postPE1(NN,:)-Cshuffle.prePE(NN,:))',Cshuffle.LatPE(NN,:)','rows','pairwise','type','Spearman');
+        fprintf('CORREL Neuron ID # %d\n',NN);
+    end
+    
+    Cshuffle.SelectionLAT(:,1)=R.Structure==10 & R.TRN(:,1)~=0 & Cshuffle.STATSpost(:,1)<PStat;
+    Cshuffle.SelectionLAT(:,2)=R.Structure==20 & R.TRN(:,1)~=0 & Cshuffle.STATSpost(:,1)<PStat;
+    % Saving data in C structure and excel file
+    save('Cshuffle.mat', 'Cshuffle');
+  
+%% Figure of distributions of means correlation coefficients and number of significant neurons
+
+selection1=R.Structure==10 & Celltype(:,1)==1;
+selection2=R.Structure==20 & Celltype(:,1)==1;
+figure;
+% Plot of distribution of mean correlation coefficients for all of the
+% shuffled data
+subplot(2,2,3);
+histogram(Cshuffle.SLOPEmeanshuffled,100);
+hold on;
+MEANsel1=nanmean(Cshuffle.SLOPEpost(selection1,1));
+MEANsel2=nanmean(Cshuffle.SLOPEpost(selection2,1));
+
+plot([MEANsel1 MEANsel1], [0 50],'Color','b');
+plot([MEANsel2 MEANsel2], [0 50],'Color','r');
+xlabel('Mean correlation coefficient') % x-axis label
+ylabel('# of iterations')% y-axis label
+%Plot of distribution of number of significantly correlated neurons across
+%all the shuffled data sets
+subplot(2,2,4);
+histogram((Cshuffle.SIGshuffled*100/sum(Celltype(:,1)==1)),40);
+xlabel('% of significantly correlated unit') % x-axis label
+ylabel('# of iterations')% y-axis label
+hold on;
+
+SigNUM1=100*sum(Cshuffle.STATSpost(selection1,1)<0.05)/(sum(selection1));
+SigNUM2=100*sum(Cshuffle.STATSpost(selection2,1)<0.05)/(sum(selection2));
+plot([SigNUM1 SigNUM1], [0 100],'Color','b');
+plot([SigNUM2 SigNUM2], [0 100],'Color','r');
+
+propSigSh=Cshuffle.SIGshuffled*100/sum(Celltype(:,1)==1);
+OTpval(1,1)=(sum(propSigSh>SigNUM1)+1)/(length(propSigSh)+1);
+OTpval(2,1)=(sum(propSigSh>SigNUM2)+1)/(length(propSigSh)+1);
+
+
+% Plot of correlation coefficients from example shuffled data set
+subplot(2,2,1);
+sel1=R.Structure~=0;
+sel2=sel1 & Cshuffle.STATSpostshuffled(:,1)<PStat;
+xmin=floor(min(Cshuffle.SLOPEpostshuffled(sel1,1)));
+xmax=ceil(max(Cshuffle.SLOPEpostshuffled(sel1,1)));
+step=0.05;
+xcenters=xmin:step:xmax;
+Nsel1=hist(Cshuffle.SLOPEpostshuffled(sel1,1),xcenters);
+Nsel2=hist(Cshuffle.SLOPEpostshuffled(sel2,1),xcenters);
+hold on;
+bar(xcenters, Nsel1, 'w', 'EdgeColor', 'k', 'BarWidth', 1);alpha(0.3);
+bar(xcenters, Nsel2, 'k', 'EdgeColor', 'k', 'BarWidth', 1);
+plot([0 0], [0 30],'LineStyle',':','Color','k');
+MEANsel1=mean(Cshuffle.SLOPEpostshuffled(sel1,1));
+plot([MEANsel1 MEANsel1], [0 10],'Color','r');
+xlabel('% correlation coefficient shuffle data') % x-axis label
+ylabel('# of neurons')% y-axis label
+hold on;
+
+% %Plot of correlation coefficients with significance marked based on
+% %bootstrapping
+subplot(2,2,2);
+sel1=R.Structure~=0;
+sel2= sel1 & Cshuffle.STATSpost(:,1)<PStat;
+xmin=floor(min(Cshuffle.SLOPEpost(sel1,1)));
+xmax=ceil(max(Cshuffle.SLOPEpost(sel1,1)));
+step=0.05;
+xcenters=xmin:step:xmax;
+Nsel1=hist(Cshuffle.SLOPEpost(sel1,1),xcenters);
+Nsel2=hist(Cshuffle.SLOPEpost(sel2,1),xcenters);
+hold on;
+bar(xcenters, Nsel1, 'w', 'EdgeColor', 'k', 'BarWidth', 1);alpha(0.3);
+bar(xcenters, Nsel2, 'k', 'EdgeColor', 'k', 'BarWidth', 1);
+plot([0 0], [0 30],'LineStyle',':','Color','k');
+MEANsel1=mean(Cshuffle.SLOPEpost(sel2,1));
+plot([MEANsel1 MEANsel1], [0 10],'Color','r');
+xlabel('% correlation coefficient real data') % x-axis label
+ylabel('# of neurons')% y-axis label