Delete 'MatlabScripts/d_Fig2S2S3S5EventsBins.m'

MatlabScripts/d_Fig2S2S3S5EventsBins.m

@@ -1,848 +0,0 @@
-clear all;
-load ('R_2R.mat');
-
-%get parameters
-BinBase=R_2R.Param.BinBase;
-BinDura=R_2R.Param.BinDura;
-bins=R_2R.Param.bins;
-binint=R_2R.Param.binint;
-binstart=R_2R.Param.binstart;
-PStatBins=0.01; %using more stringent cutoff to reduce pre-delivery noise
-
-%divide neurons up by region
-NAneurons=strcmp(R_2R.Ninfo(:,3),'NA');
-VPneurons=strcmp(R_2R.Ninfo(:,3),'VP');
-
-%which bins bound the area examined for reward-selectivity? (in seconds)
-Time1=0.4; %seconds
-Time2=3; %seconds
-Bin1=(((Time1-BinDura(2)/2)-binstart)/binint); %convert to bin name
-Bin2=(((Time2-BinDura(2)/2)-binstart)/binint); %convert to bin name
-
-%sorting bin -- which bin the neurons' activity is sorted on for heatmap(in seconds)
-SortBinTime=1.1; %seconds
-SortBin=round((((SortBinTime-BinDura(2)/2)-binstart)/binint)); %convert to bin name
-
-sucrose=[1  0.6  0.1];
-maltodextrin=[.9  0.3  .9];
-water=[0.00  0.75  0.75];
-total=[0.3  0.1  0.8];
-exc=[0 113/255 188/255];
-inh=[240/255 0 50/255];
-
-%% bins analysis
-
-%first and last bin aren't included because you can't compare to the previous/subsequent bin
-%this axis plots the bins on their centers
-xaxis=linspace(binstart+binint+BinDura(2)/2,binstart+(bins-2)*binint+BinDura(2)/2,bins-2);
-
-for i=2:(bins-1) %no including first or last bin because can't compare to previous/subsequent bin
-    
-    %finds out whether firing is stronger (high excitation or lower inhibition) for 1 or 2
-    for k=1:length(R_2R.Ninfo) %runs through neurons
-        if R_2R.BinStatRwrd{i,1}(k).IntSig < PStatBins %if neuron is significant for this bin
-            if R_2R.BinStatRwrd{i-1,1}(k).IntSig < PStatBins || R_2R.BinStatRwrd{i+1,1}(k).IntSig < PStatBins %either previous or subsequent bin must be significant
-                if R_2R.BinStatRwrd{i,1}(k).R1Mean > R_2R.BinStatRwrd{i,1}(k).R2Mean %if firing is greater on sucrose trials than maltodextrin
-                    R_2R.BinRewPref{i,1}(k,1)=1;
-                else
-                    R_2R.BinRewPref{i,1}(k,1)=2; %otherwise maltodextrin must be greater than sucrose
-                end
-            else
-                R_2R.BinRewPref{i,1}(k,1)=0; %if not significant in 2 consecutive bins
-            end
-        else
-            R_2R.BinRewPref{i,1}(k,1)=0; %if no sig reward modulation
-        end
-        
-    end
-    %find how many NAc neurons have significant reward modulation in each bin
-    NN1perBin(i,1)=sum(R_2R.BinRewPref{i,1}(NAneurons,1)==1); %sucrose pref
-    NN2perBin(i,1)=sum(R_2R.BinRewPref{i,1}(NAneurons,1)==2); %malto pref
-    NNperBin(i,1)=sum(R_2R.BinRewPref{i,1}(NAneurons,1)>0); %either
-    %normalize to number of neurons in population
-    NN1norm=NN1perBin./sum(NAneurons);
-    NN2norm=NN2perBin./sum(NAneurons);
-    NNnorm=NNperBin./sum(NAneurons);
-    
-    %find how many VP neurons have significant reward modulation in each bin
-    NV1perBin(i,1)=sum(R_2R.BinRewPref{i,1}(VPneurons,1)==1); %sucrose pref
-    NV2perBin(i,1)=sum(R_2R.BinRewPref{i,1}(VPneurons,1)==2); %malto pref
-    NVperBin(i,1)=sum(R_2R.BinRewPref{i,1}(VPneurons,1)>0); %either
-    %normalize to number of neurons in population
-    NV1norm=NV1perBin./sum(VPneurons);
-    NV2norm=NV2perBin./sum(VPneurons);
-    NVnorm=NVperBin./sum(VPneurons);
-   
-
-end
-
-%Cumulative reward selectivity
-cumsel=zeros(length(R_2R.Ninfo),bins); %set up result matrix
-cumsel1=zeros(length(R_2R.Ninfo),bins); %set up result matrix for sucrose
-cumsel2=zeros(length(R_2R.Ninfo),bins); %set up result matrix for maltodextrin
-%note, this has to be corrected to +1 because of the way the bin matrix is
-%layed out (bin 0 is the first column, not the 0th column)
-for i=1:length(R_2R.Ninfo) %go through each neuron
-    for j=2:bins-1 %look in each bin we did analysis on
-        if R_2R.BinRewPref{j,1}(i,1)>0 || cumsel(i,j-1)==1
-            cumsel(i,j) = 1;
-        end
-        if R_2R.BinRewPref{j,1}(i,1)==1 || cumsel1(i,j-1)==1
-            cumsel1(i,j) = 1;
-        end
-        if R_2R.BinRewPref{j,1}(i,1)==2 || cumsel2(i,j-1)==1
-            cumsel2(i,j) = 1;
-        end        
-    end
-end
-
-%plotting number of significantly modulated neurons across time
-%NAc
-subplot(2,3,1);
-hold on;
-plot(xaxis,NNnorm(2:bins-1),'Color', total,'linewidth',1.5);
-plot(xaxis,NN1norm(2:bins-1),'Color',sucrose,'linewidth',1.5);
-plot(xaxis,NN2norm(2:bins-1),'Color',maltodextrin,'linewidth',1.5);
-plot([Time1 Time1],[-1 1],':','color','k','linewidth',0.75);
-plot([Time2 Time2],[-1 1],':','color','k','linewidth',0.75);
-axis([xaxis(1) xaxis(end) 0 0.4]);
-legend('Total','Suc > mal','Mal > suc','Location','northeast')
-ylabel('Fraction of population');
-title('NAc reward-selective neurons');
-
-%VP
-subplot(2,3,2);
-hold on;
-plot(xaxis,NVnorm(2:bins-1),'Color', total,'linewidth',1.5);
-plot(xaxis,NV1norm(2:bins-1),'Color',sucrose,'linewidth',1.5);
-plot(xaxis,NV2norm(2:bins-1),'Color',maltodextrin,'linewidth',1.5);
-plot([Time1 Time1],[-1 1],':','color','k','linewidth',0.75);
-plot([Time2 Time2],[-1 1],':','color','k','linewidth',0.75);
-axis([xaxis(1) xaxis(end) 0 0.4]);
-legend('Total','Suc > mal','Mal > suc','Location','northeast')
-ylabel('Fraction of population');
-title('VP reward-selective neurons');
-
-%plotting cumulative selectivity in each region
-%NAc
-subplot(2,3,4);
-hold on;
-plot(xaxis,sum(cumsel(NAneurons,2:bins-1),1)/sum(NAneurons),'Color', total,'linewidth',1.5);
-plot(xaxis,sum(cumsel1(NAneurons,2:bins-1),1)/sum(NAneurons),'Color', sucrose,'linewidth',1.5);
-plot(xaxis,sum(cumsel2(NAneurons,2:bins-1),1)/sum(NAneurons),'Color', maltodextrin,'linewidth',1.5);
-plot([Time1 Time1],[0 1],':','color','k','linewidth',0.75);
-plot([Time2 Time2],[0 1],':','color','k','linewidth',0.75);
-axis([xaxis(1) xaxis(end) 0 0.75]);
-xlabel('Seconds from reward delivery');
-ylabel('Cumulative fraction');
-
-%VP
-subplot(2,3,5);
-hold on;
-plot(xaxis,sum(cumsel(VPneurons,2:bins-1),1)/sum(VPneurons),'Color', total,'linewidth',1.5);
-plot(xaxis,sum(cumsel1(VPneurons,2:bins-1),1)/sum(VPneurons),'Color', sucrose,'linewidth',1.5);
-plot(xaxis,sum(cumsel2(VPneurons,2:bins-1),1)/sum(VPneurons),'Color', maltodextrin,'linewidth',1.5);
-plot([Time1 Time1],[0 1],':','color','k','linewidth',0.75);
-plot([Time2 Time2],[0 1],':','color','k','linewidth',0.75);
-axis([xaxis(1) xaxis(end) 0 0.75]);
-xlabel('Seconds from RD');
-ylabel('Cumulative fraction');
-
-%plot difference between the two regions for each measure
-%reward selective neurons per bin
-subplot(2,3,3);
-hold on;
-plot(xaxis,NNnorm(2:bins-1)-NVnorm(2:bins-1),'Color', total,'linewidth',1.5);
-plot(xaxis,NN1norm(2:bins-1)-NV1norm(2:bins-1),'Color',sucrose,'linewidth',1.5);
-plot(xaxis,NN2norm(2:bins-1)-NV2norm(2:bins-1),'Color',maltodextrin,'linewidth',1.5);
-axis([xaxis(1) xaxis(end) -0.3 0.3]);
-plot([xaxis(1) xaxis(end)],[0 0],'color','k','linewidth',0.75);
-plot([Time1 Time1],[-1 1],':','color','k','linewidth',0.75);
-plot([Time2 Time2],[-1 1],':','color','k','linewidth',0.75);
-legend('Total','Suc > mal','Mal > suc','Location','northeast')
-title('Difference (NAc - VP)');
-ylabel('Fraction of population');
-
-%cumulative selectivity
-subplot(2,3,6);
-hold on;
-plot(xaxis,sum(cumsel(NAneurons,2:bins-1),1)/sum(NAneurons)-sum(cumsel(VPneurons,2:bins-1),1)/sum(VPneurons),'Color', total,'linewidth',1.5);
-plot(xaxis,sum(cumsel1(NAneurons,2:bins-1),1)/sum(NAneurons)-sum(cumsel1(VPneurons,2:bins-1),1)/sum(VPneurons),'Color', sucrose,'linewidth',1.5);
-plot(xaxis,sum(cumsel2(NAneurons,2:bins-1),1)/sum(NAneurons)-sum(cumsel2(VPneurons,2:bins-1),1)/sum(VPneurons),'Color', maltodextrin,'linewidth',1.5);
-plot([Time1 Time1],[-1 1],':','color','k','linewidth',0.75);
-plot([Time2 Time2],[-1 1],':','color','k','linewidth',0.75);
-axis([xaxis(1) xaxis(end) -Time1 Time1]);
-plot([xaxis(1) xaxis(end)],[0 0],'color','k','linewidth',0.75);
-xlabel('Seconds from RD');
-ylabel('Cumulative fraction');
-
-%% individual rat data
-figure;
-
-%split data up into individual rats
-for i=1:length(R_2R.Ninfo)
-    A=char(R_2R.Ninfo(i,4));
-    Rats(i,:)=cellstr(A);
-end
-
-C=unique(Rats);
-
-for i=1:length(C)
-    Sel=strcmp(C(i),Rats);
-    for j=2:bins-1
-        NRat1(j,1)=sum(R_2R.BinRewPref{j,1}(Sel,1)==1)/sum(Sel);
-        NRat2(j,1)=sum(R_2R.BinRewPref{j,1}(Sel,1)==2)/sum(Sel);
-        NRat(j,1)=sum(R_2R.BinRewPref{j,1}(Sel,1)>0)/sum(Sel);
-    end
-    
-    %graphing
-    subplot(4,3,i);
-    hold on;
-    plot(xaxis,NRat(2:bins-1),'Color', total,'linewidth',1.5);
-    plot(xaxis,NRat1(2:65),'Color',sucrose,'linewidth',1.5);
-    plot(xaxis,NRat2(2:65),'Color',maltodextrin,'linewidth',1.5);
-    xlabel('Seconds from reward delivery');
-    title([cell2mat(C(i)) ' (n=' num2str(sum(Sel)) ')'])
-    axis([xaxis(1) xaxis(end) 0 0.5]);
-    legend('Total neurons','Suc > mal Hz','Mal > suc Hz','Location','northeast')
-    ylabel('Fraction of total population');
-
-end
-
-%% plotting sucrose-selective neurons
-
-figure;
-
-%color map
-[magma,inferno,plasma,viridis]=colormaps;
-colormap(plasma);
-c=[-100 2000];ClimE=sign(c).*abs(c).^(1/4);%colormap
-
-%events we're looking at
-ev1=strcmp('RD1', R_2R.Erefnames);
-ev2=strcmp('RD2', R_2R.Erefnames);
-
-%setting up parameters
-Xaxis=[-2 5];
-inttime=find(R_2R.Param.Tm>=Xaxis(1) & R_2R.Param.Tm<=Xaxis(2));
-paramtime=R_2R.Param.Tm(inttime);
-
-%find all neurons with greater firing for sucrose
-for i = 1:(Bin2-Bin1+1)
-    %the added +1 is necessary because bin 20 is the 21st entry in the matrix
-    Pref1(:,i)=R_2R.BinRewPref{Bin1+i}==1; %get neurons that have greater firing for sucrose in any of the bins bounded above
-    Resp11(:,i)=Pref1(:,i)&cat(1,R_2R.BinStatRwrd{Bin1+i,1}.SucRespDir)==1; %get neurons with excitation to sucrose
-    Resp12(:,i)=Pref1(:,i)&cat(1,R_2R.BinStatRwrd{Bin1+i,1}.MalRespDir)==-1;%get neurons with inhibition to maltodextrin
-end
-
-Sel=sum(Pref1,2)>0;  %all neurons selective for sucrose in any bin
-Sel2=sum(Resp11,2)>0; %all selective neurons sucrose excited in any bin
-Sel3=sum(Resp12,2)>0; %all selective neurons malto inhibited in any bin
-
-%-----------------NAc--------------------------------------------
-subplot(4,40,[15:23]); %heatmap of sucrose preferring neurons' response to sucrose
-Neurons=R_2R.Ev(ev1).PSTHz(Sel&NAneurons,inttime); %get the firing rates of neurons of interest
-SucResp=cat(1,R_2R.BinStatRwrd{SortBin+1,1}.R1Mean); %sucrose responses
-TMP=SucResp(Sel&NAneurons); %find the magnitude of the excitations for sucrose for this bin
-TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
-[~,SORTimg]=sort(TMP);
-Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
-imagesc(paramtime,[1,sum(Sel&NAneurons,1)],Neurons,ClimE);
-title(['Sucrose trials']);
-ylabel('Individual neurons');
-hold on;
-plot([0 0],[0 sum(Sel&NAneurons)],':','color','k','linewidth',0.75);
-
-subplot(4,40,[27:35]); %heatmap of sucrose preferring neurons' response to maltodextrin
-Neurons=R_2R.Ev(ev2).PSTHz(Sel&NAneurons,inttime); %get the firing rates of neurons of interest
-Neurons=Neurons(SORTimg,:); %sort the neurons same as before
-imagesc(paramtime,[1,sum(Sel&NAneurons,1)],Neurons,ClimE);
-title(['Maltodextrin trials']);
-hold on;
-plot([0 0],[0 sum(Sel&NAneurons)],':','color','k','linewidth',0.75);
-
-%plot sucrose preferring neurons to sucrose
-psth1=nanmean(R_2R.Ev(ev1).PSTHz(Sel&NAneurons,inttime),1); 
-sem1=nanste(R_2R.Ev(ev1).PSTHz(Sel&NAneurons,inttime),1); %calculate standard error of the mean
-up1=psth1+sem1;
-down1=psth1-sem1;
-
-%plot sucrose preferring neurons to malt
-psth2=nanmean(R_2R.Ev(ev2).PSTHz(Sel&NAneurons,inttime),1); 
-sem2=nanste(R_2R.Ev(ev2).PSTHz(Sel&NAneurons,inttime),1); %calculate standard error of the mean
-up2=psth2+sem2;
-down2=psth2-sem2;
-
-%plotting
-subplot(4,3,1);
-title(['Suc>mal (n=' num2str(sum(Sel&NAneurons)) ' of ' num2str(sum(NAneurons)) ')'])
-hold on;
-plot(paramtime,psth1,'Color',sucrose,'linewidth',1);
-plot(paramtime,psth2,'Color',maltodextrin,'linewidth',1);
-patch([paramtime,paramtime(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
-patch([paramtime,paramtime(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
-plot([-2 5],[0 0],':','color','k','linewidth',0.75);
-plot([0 0],[-1 10],':','color','k','linewidth',0.75);
-axis([-2 5 -1 4]);
-ylabel('Mean firing (z-score)');
-legend('Suc trials','Mal trials');
-
-%display inhibitions and excitations
-both = sum(Sel2&Sel3&NAneurons);
-excite = sum(Sel2&NAneurons)-both;
-inhib = sum(Sel3&NAneurons)-both;
-
-subplot(4,40,40);
-b = bar([inhib both excite; 0 0 0],'stacked');
-b(1,1).FaceColor=maltodextrin;
-b(1,2).FaceColor=total;
-b(1,3).FaceColor=sucrose;
-axis([1 1.2 0 both+excite+inhib]);
-ylabel('Inh / Both / Excited');
-set(gca,'xtick',[])
-
-
-%----------------VP----------------
-subplot(4,40,[80+15:80+23]); %heatmap of sucrose preferring neurons' response to sucrose
-Neurons=R_2R.Ev(ev1).PSTHz(Sel&VPneurons,inttime); %get the firing rates of neurons of interest
-SucResp=cat(1,R_2R.BinStatRwrd{SortBin+1,1}.R1Mean); %sucrose responses
-TMP=SucResp(Sel&VPneurons); %find the magnitude of the excitations for sucrose for this bin
-TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
-[~,SORTimg]=sort(TMP);
-Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
-imagesc(paramtime,[1,sum(Sel&VPneurons,1)],Neurons,ClimE);
-title(['Sucrose trials']);
-ylabel('Individual neurons');
-hold on;
-plot([0 0],[0 sum(Sel&VPneurons)],':','color','k','linewidth',0.75);
-
-subplot(4,40,[80+27:80+35]); %heatmap of sucrose preferring neurons' response to maltodextrin
-Neurons=R_2R.Ev(ev2).PSTHz(Sel&VPneurons,inttime); %get the firing rates of neurons of interest
-Neurons=Neurons(SORTimg,:); %sort the neurons same as before
-imagesc(paramtime,[1,sum(Sel&VPneurons,1)],Neurons,ClimE);
-title(['Maltodextrin trials']);
-hold on;
-plot([0 0],[0 sum(Sel&VPneurons)],':','color','k','linewidth',0.75);
-
-%plot sucrose preferring neurons to sucrose
-psth1=nanmean(R_2R.Ev(ev1).PSTHz(Sel&VPneurons,inttime),1); 
-sem1=nanste(R_2R.Ev(ev1).PSTHz(Sel&VPneurons,inttime),1); %calculate standard error of the mean
-up1=psth1+sem1;
-down1=psth1-sem1;
-
-%plot sucrose preferring neurons to malt
-psth2=nanmean(R_2R.Ev(ev2).PSTHz(Sel&VPneurons,inttime),1); 
-sem2=nanste(R_2R.Ev(ev2).PSTHz(Sel&VPneurons,inttime),1); %calculate standard error of the mean
-up2=psth2+sem2;
-down2=psth2-sem2;
-
-%plotting
-subplot(4,3,7);
-title(['Suc>mal (n=' num2str(sum(Sel&VPneurons)) ' of ' num2str(sum(VPneurons)) ')'])
-hold on;
-plot(paramtime,psth1,'Color',sucrose,'linewidth',1);
-plot(paramtime,psth2,'Color',maltodextrin,'linewidth',1);
-patch([paramtime,paramtime(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
-patch([paramtime,paramtime(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
-plot([-2 5],[0 0],':','color','k','linewidth',0.75);
-plot([0 0],[-1 10],':','color','k','linewidth',0.75);
-axis([-2 5 -0.6 4]);
-ylabel('Mean firing (z-score)');
-legend('Suc trials','Mal trials');
-
-%display inhibitions and excitations
-both = sum(Sel2&Sel3&VPneurons);
-excite = sum(Sel2&VPneurons)-both;
-inhib = sum(Sel3&VPneurons)-both;
-
-subplot(4,40,120);
-b = bar([inhib both excite; 0 0 0],'stacked');
-b(1,1).FaceColor=maltodextrin;
-b(1,2).FaceColor=total;
-b(1,3).FaceColor=sucrose;
-axis([1 1.2 0 both+excite+inhib]);
-ylabel('Inh / Both / Excited');
-set(gca,'xtick',[])
-
-%% plotting maltodextrin-selective neurons
-
-%find all neurons with greater firing for maltodextrin
-for i = 1:(Bin2-Bin1+1)
-    %the added +1 is necessary because bin 20 is the 21st entry in the matrix
-    Pref2(:,i)=R_2R.BinRewPref{Bin1+i}==2; %get neurons that have greater firing for maltodextrin in any of the bins bounded above
-    Resp21(:,i)=Pref2(:,i)&cat(1,R_2R.BinStatRwrd{Bin1+i,1}.SucRespDir)==-1; %get neurons with inhibition to sucrose
-    Resp22(:,i)=Pref2(:,i)&cat(1,R_2R.BinStatRwrd{Bin1+i,1}.MalRespDir)==1;%get neurons with excitation to maltodextrin
-end
-
-Mel=sum(Pref2,2)>0;  %all neurons selective for malotdextrin in any bin
-Mel2=sum(Resp21,2)>0; %all selective neurons sucrose inhibited in any bin
-Mel3=sum(Resp22,2)>0; %all selective neurons malto excited in any bin
-
-%-----------------NAc--------------------------------------------
-subplot(4,40,[40+15:40+23]); %heatmap of maltodextrin preferring neurons' response to sucrose
-Neurons=R_2R.Ev(ev1).PSTHz(Mel&NAneurons,inttime); %get the firing rates of neurons of interest
-MalResp=cat(1,R_2R.BinStatRwrd{SortBin+1,1}.R2Mean); %sucrose responses
-TMP=MalResp(Mel&NAneurons); %find the magnitude of the excitations for sucrose for this bin
-TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
-[~,SORTimg]=sort(TMP);
-Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
-imagesc(paramtime,[1,sum(Mel&NAneurons,1)],Neurons,ClimE);
-title(['Sucrose trials']);
-xlabel('Seconds from RD');
-ylabel('Individual neurons');
-hold on;
-plot([0 0],[0 sum(Mel&NAneurons)],':','color','k','linewidth',0.75);
-
-subplot(4,40,[40+27:40+35]); %heatmap of maltodextrin preferring neurons' response to maltodextrin
-Neurons=R_2R.Ev(ev2).PSTHz(Mel&NAneurons,inttime); %get the firing rates of neurons of interest
-Neurons=Neurons(SORTimg,:); %sort the neurons same as before
-imagesc(paramtime,[1,sum(Mel&NAneurons,1)],Neurons,ClimE);
-title(['Maltodextrin trials']);
-xlabel('Seconds from RD');
-hold on;
-plot([0 0],[0 sum(Mel&NAneurons)],':','color','k','linewidth',0.75);
-
-%plot sucrose preferring neurons to sucrose
-psth1=nanmean(R_2R.Ev(ev1).PSTHz(Mel&NAneurons,inttime),1); 
-sem1=nanste(R_2R.Ev(ev1).PSTHz(Mel&NAneurons,inttime),1); %calculate standard error of the mean
-up1=psth1+sem1;
-down1=psth1-sem1;
-
-%plot sucrose preferring neurons to malt
-psth2=nanmean(R_2R.Ev(ev2).PSTHz(Mel&NAneurons,inttime),1); 
-sem2=nanste(R_2R.Ev(ev2).PSTHz(Mel&NAneurons,inttime),1); %calculate standard error of the mean
-up2=psth2+sem2;
-down2=psth2-sem2;
-
-%plotting
-subplot(4,3,4);
-title(['Mal>suc (n=' num2str(sum(Mel&NAneurons)) ' of ' num2str(sum(NAneurons)) ')'])
-hold on;
-plot(paramtime,psth1,'Color',sucrose,'linewidth',1);
-plot(paramtime,psth2,'Color',maltodextrin,'linewidth',1);
-patch([paramtime,paramtime(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
-patch([paramtime,paramtime(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
-plot([-2 5],[0 0],':','color','k','linewidth',0.75);
-plot([0 0],[-2 7],':','color','k','linewidth',0.75);
-axis([-2 5 -1.2 3]);
-ylabel('Mean firing (z-score)');
-xlabel('Seconds from RD');
-legend('Suc trials','Mal trials','Location','northeast');
-
-%display inhibitions and excitations
-both = sum(Mel2&Mel3&NAneurons);
-excite = sum(Mel3&NAneurons)-both;
-inhib = sum(Mel2&NAneurons)-both;
-
-subplot(4,40,80);
-b = bar([inhib both excite; 0 0 0],'stacked');
-b(1,1).FaceColor=sucrose;
-b(1,2).FaceColor=total;
-b(1,3).FaceColor=maltodextrin;
-axis([1 1.2 0 both+excite+inhib]);
-ylabel('Inh / Both / Excited');
-set(gca,'xtick',[]);
-
-
-%----------------VP----------------
-subplot(4,40,[120+15:120+23]); %heatmap of maltodextrin preferring neurons' response to sucrose
-Neurons=R_2R.Ev(ev1).PSTHz(Mel&VPneurons,inttime); %get the firing rates of neurons of interest
-MalResp=cat(1,R_2R.BinStatRwrd{SortBin+1,1}.R2Mean); %sucrose responses
-TMP=MalResp(Mel&VPneurons); %find the magnitude of the excitations for sucrose for this bin
-TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
-[~,SORTimg]=sort(TMP);
-Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
-imagesc(paramtime,[1,sum(Mel&VPneurons,1)],Neurons,ClimE);
-title(['Sucrose trials']);
-xlabel('Seconds from RD');
-ylabel('Individual neurons');
-hold on;
-plot([0 0],[0 sum(Mel&VPneurons)],':','color','k','linewidth',0.75);
-
-subplot(4,40,[120+27:120+35]); %heatmap of maltodextrin preferring neurons' response to maltodextrin
-Neurons=R_2R.Ev(ev2).PSTHz(Mel&VPneurons,inttime); %get the firing rates of neurons of interest
-Neurons=Neurons(SORTimg,:); %sort the neurons same as before
-imagesc(paramtime,[1,sum(Mel&VPneurons,1)],Neurons,ClimE);
-title(['Maltodextrin trials']);
-xlabel('Seconds from RD');
-hold on;
-plot([0 0],[0 sum(Mel&VPneurons)],':','color','k','linewidth',0.75);
-
-%plot maltodextrin preferring neurons to sucrose
-psth1=nanmean(R_2R.Ev(ev1).PSTHz(Mel&VPneurons,inttime),1); 
-sem1=nanste(R_2R.Ev(ev1).PSTHz(Mel&VPneurons,inttime),1); %calculate standard error of the mean
-up1=psth1+sem1;
-down1=psth1-sem1;
-
-%plot malt preferring neurons to malt
-psth2=nanmean(R_2R.Ev(ev2).PSTHz(Mel&VPneurons,inttime),1); 
-sem2=nanste(R_2R.Ev(ev2).PSTHz(Mel&VPneurons,inttime),1); %calculate standard error of the mean
-up2=psth2+sem2;
-down2=psth2-sem2;
-
-%plotting
-subplot(4,3,10);
-title(['Mal>suc (n=' num2str(sum(Mel&VPneurons)) ' of ' num2str(sum(VPneurons)) ')'])
-hold on;
-plot(paramtime,psth1,'Color',sucrose,'linewidth',1);
-plot(paramtime,psth2,'Color',maltodextrin,'linewidth',1);
-patch([paramtime,paramtime(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
-patch([paramtime,paramtime(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
-plot([-2 5],[0 0],':','color','k','linewidth',0.75);
-plot([0 0],[-2 7],':','color','k','linewidth',0.75);
-axis([-2 5 -1.5 2.5]);
-ylabel('Mean firing (z-score)');
-xlabel('Seconds from RD');
-legend('Suc trials','Mal trials','Location','northeast');
-
-%display inhibitions and excitations
-both = sum(Mel2&Mel3&VPneurons);
-excite = sum(Mel3&VPneurons)-both;
-inhib = sum(Mel2&VPneurons)-both;
-
-subplot(4,40,160);
-b = bar([inhib both excite; 0 0 0],'stacked');
-b(1,1).FaceColor=sucrose;
-b(1,2).FaceColor=total;
-b(1,3).FaceColor=maltodextrin;
-axis([1 1.2 0 both+excite+inhib]);
-ylabel('Inh / Both / Excited');
-set(gca,'xtick',[]);
-
-%% creating and saving new variables
-R_2R.RSinfo=R_2R.Ninfo(Sel|Mel,1:2); %reward-selective neurons
-R_2R.RNSinfo=R_2R.Ninfo(Sel==0 & Mel==0,1:2); %reward-nonselective neurons
-R_2R.SucN=Sel;
-R_2R.MalN=Mel;
-
-% proportion summary and chi squared test for reward selective neurons
-R_2R.RegSelX2(1,1)=sum((Sel | Mel) & NAneurons)/sum(NAneurons); %how many selective NAc neurons
-R_2R.RegSelX2(1,2)=sum((Sel | Mel) & VPneurons)/sum(VPneurons); %how many selective VP neurons
-[~,R_2R.RegSelX2(2,1),R_2R.RegSelX2(2,2)]=crosstab(Sel | Mel,NAneurons);
-
-% proportion summary and chi squared test for max neurons in any bin
-R_2R.RegMaxSelX2(1,1)=max(NNnorm); %how many selective NAc neurons
-R_2R.RegMaxSelX2(1,2)=max(NVnorm); %how many selective VP neurons
-[~,R_2R.RegMaxSelX2(2,1),R_2R.RegMaxSelX2(2,2)]=crosstab(cat(1,ones(max(NNperBin),1),zeros(sum(NAneurons)-max(NNperBin),1),ones(max(NVperBin),1),zeros(sum(VPneurons)-max(NVperBin),1)),NAneurons);
-
-
-save('R_2R.mat','R_2R');
-
-%divide neurons up by region
-NAneurons=strcmp(R_2R.Ninfo(:,3),'NA');
-VPneurons=strcmp(R_2R.Ninfo(:,3),'VP');
-
-Xaxis=[-1 2];
-Ishow=find(R_2R.Param.Tm>=Xaxis(1) & R_2R.Param.Tm<=Xaxis(2));
-time1=R_2R.Param.Tm(Ishow);
-Xaxis2=[-2 2];
-Ushow=find(R_2R.Param.Tm>=Xaxis2(1) & R_2R.Param.Tm<=Xaxis2(2));
-time2=R_2R.Param.Tm(Ushow);
-
-inh=[0.1 0.021154 0.6];
-exc=[0.9 0.75 0.205816];
-
-Cue=strcmp('Cue', R_2R.Erefnames);
-PE=strcmp('PECue', R_2R.Erefnames);
-RD=strcmp('RD', R_2R.Erefnames); %should this be LickR or RD?
-%% Plotting heatmaps of all neurons to each event
-figure;
-
-%color map
-[magma,inferno,plasma,viridis]=colormaps;
-colormap(plasma);
-c=[-100 2000];ClimE=sign(c).*abs(c).^(1/4);%colormap
-
-for i=1:2
-    
-    if i==1 Reg=NAneurons; end %1 is NAc
-    if i==2 Reg=VPneurons; end %2 is VP
-    
-    Sel=Reg&R_2R.Ev(Cue).RespDir<2; %gets all neurons
-
-    %sort cue heatmap by magnitude of response
-    Neurons=R_2R.Ev(Cue).PSTHz(Sel,Ishow); %get the firing rates of neurons of interest
-    TMP=R_2R.Ev(Cue).Meanz(Sel); %find the magnitude of the inhibitions for this event
-    TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
-    [~,SORTimg]=sort(TMP);
-    Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
-
-    %heatmap
-    subplot(4,6,[1 7]+(i-1)*12);
-    imagesc(time1,[1,sum(Sel,1)],Neurons,ClimE); title('Cue responses');
-    xlabel('Seconds post cue');
-    ylabel('Individual neurons sorted by response strength');
-    hold on;
-    plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
-
-    %sort PE heatmap by magnitude of response
-    Neurons=R_2R.Ev(PE).PSTHz(Sel,Ushow); %get the firing rates of neurons of interest
-    TMP=R_2R.Ev(PE).Meanz(Sel); %find the magnitude of the inhibitions for this event
-    TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
-    [~,SORTimg]=sort(TMP);
-    Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
-
-    %heatmap
-    subplot(4,6,[3 9]+(i-1)*12);
-    imagesc(time2,[1,sum(Sel,1)],Neurons,ClimE); title('PE responses');
-    xlabel('Seconds post PE');
-    hold on;
-    plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
-
-    %sort RD heatmap by magnitude of response
-    Neurons=R_2R.Ev(RD).PSTHz(Sel,Ushow); %get the firing rates of neurons of interest
-    TMP=R_2R.Ev(RD).Meanz(Sel); %find the magnitude of the inhibitions for this event
-    TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
-    [~,SORTimg]=sort(TMP);
-    Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
-
-    %heatmap
-    subplot(4,6,[5 11]+(i-1)*12);
-    imagesc(time2,[1,sum(Sel,1)],Neurons,ClimE); title('Reward responses');
-    xlabel('Seconds post RD');
-    hold on;
-    plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
-
-    %% Plotting neurons that respond to the cue
-
-    %inhibitions
-    Sel=Reg&R_2R.Ev(Cue).RespDir<0; %Find neurons that have an inhibitory response to this event
-
-    %average firing rate
-    subplot(4,6,2+(i-1)*12);
-    psthI=nanmean(R_2R.Ev(Cue).PSTHz(Sel,Ishow),1); %find the average firing rate of the neurons at the time of the event
-    semI=nanste(R_2R.Ev(Cue).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
-    upI=psthI+semI;
-    downI=psthI-semI;
-    hold on;
-    patch([time1,time1(end:-1:1)],[upI,downI(end:-1:1)],inh,'EdgeColor','none');alpha(0.5);
-    plot(time1,psthI,'Color',inh,'linewidth',1); title('Mean firing (z-score)'); %create plot of avg firing rate
-    plot([-1 2],[0 0],':','color','k','linewidth',0.75);
-    plot([0 0],[-3 9],':','color','k','linewidth',0.75);
-    axis([-1 2 min(downI)-0.2 max(upI)+0.2]);
-    title(['Cue inhibitions (' num2str(round(100*sum(Sel)/sum(Reg))) '%)'])
-    ylabel('Z-score');
-
-    %excitation
-    Sel=Reg&R_2R.Ev(Cue).RespDir>0; %Find neurons that have an excitatory response to this event
-
-    %average firing rate
-    subplot(4,6,8+(i-1)*12);
-    psthE=nanmean(R_2R.Ev(Cue).PSTHz(Sel,Ishow),1);
-    semE=nanste(R_2R.Ev(Cue).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
-    upE=psthE+semE;
-    downE=psthE-semE;
-    hold on;
-    patch([time1,time1(end:-1:1)],[upE,downE(end:-1:1)],exc,'EdgeColor','none');alpha(0.5);
-    plot(time1,psthE,'Color',exc,'linewidth',1);
-    plot([-1 2],[0 0],':','color','k','linewidth',0.75);
-    plot([0 0],[-3 9],':','color','k','linewidth',0.75);
-    axis([-1 2 min(downE)-0.2 max(upE)+0.3]);
-    title(['Cue excitations (' num2str(round(100*sum(Sel)/sum(Reg))) '%)'])
-    ylabel('Z-score');
-
-    %% Plotting neurons that respond to PE
-
-    %inhibitions
-    Sel=Reg&R_2R.Ev(PE).RespDir<0; %Find neurons that have an inhibitory response to this event
-
-    %average firing rate
-    subplot(4,6,4+(i-1)*12);
-    psthI=nanmean(R_2R.Ev(PE).PSTHz(Sel,Ushow),1); %find the average firing rate of the neurons at the time of the event
-    semI=nanste(R_2R.Ev(PE).PSTHz(Sel,Ushow),1); %calculate standard error of the mean
-    upI=psthI+semI;
-    downI=psthI-semI;
-    hold on;
-    patch([time2,time2(end:-1:1)],[upI,downI(end:-1:1)],inh,'EdgeColor','none');alpha(0.5);
-    plot(time2,psthI,'Color',inh,'linewidth',1); title('Mean firing (z-score)'); %create plot of avg firing rate
-    plot([-2 2],[0 0],':','color','k','linewidth',0.75);
-    plot([0 0],[-3 9],':','color','k','linewidth',0.75);
-    axis([-2 2 min(downI)-0.2 max(upI)+0.2]);
-    title(['PE inhibitions (' num2str(round(100*sum(Sel)/sum(Reg))) '%)'])
-    ylabel('Z-score');
-
-    %excitation
-    Sel=Reg&R_2R.Ev(PE).RespDir>0; %Find neurons that have an excitatory response to this event
-
-    %average firing rate
-    subplot(4,6,10+(i-1)*12);
-    psthE=nanmean(R_2R.Ev(PE).PSTHz(Sel,Ushow),1);
-    semE=nanste(R_2R.Ev(PE).PSTHz(Sel,Ushow),1); %calculate standard error of the mean
-    upE=psthE+semE;
-    downE=psthE-semE;
-    hold on;
-    patch([time2,time2(end:-1:1)],[upE,downE(end:-1:1)],exc,'EdgeColor','none');alpha(0.5);
-    plot(time2,psthE,'Color',exc,'linewidth',1);
-    plot([-2 2],[0 0],':','color','k','linewidth',0.75);
-    plot([0 0],[-3 9],':','color','k','linewidth',0.75);
-    axis([-2 2 min(downE)-0.2 max(upE)+0.3]);
-    title(['PE excitations (' num2str(round(100*sum(Sel)/sum(Reg))) '%)'])
-    ylabel('Z-score');
-
-    %% Plotting neurons that respond to reward
-
-    %inhibitions
-    Sel=Reg&R_2R.Ev(RD).RespDir<0; %Find neurons that have an inhibitory response to this event
-
-    %average firing rate
-    subplot(4,6,6+(i-1)*12);
-    psthI=nanmean(R_2R.Ev(RD).PSTHz(Sel,Ushow),1); %find the average firing rate of the neurons at the time of the event
-    semI=nanste(R_2R.Ev(RD).PSTHz(Sel,Ushow),1); %calculate standard error of the mean
-    upI=psthI+semI;
-    downI=psthI-semI;
-    hold on;
-    patch([time2,time2(end:-1:1)],[upI,downI(end:-1:1)],inh,'EdgeColor','none');alpha(0.5);
-    plot(time2,psthI,'Color',inh,'linewidth',1); title('Mean firing (z-score)'); %create plot of avg firing rate
-    plot([-1 2],[0 0],':','color','k','linewidth',0.75);
-    plot([0 0],[-3 9],':','color','k','linewidth',0.75);
-    axis([-1 2 min(downI)-0.2 max(upI)+0.2]);
-    title(['RD inhibitions (' num2str(round(100*sum(Sel)/sum(Reg))) '%)'])
-    ylabel('Z-score');
-
-    %excitation
-    Sel=Reg&R_2R.Ev(RD).RespDir>0; %Find neurons that have an excitatory response to this event
-
-    %average firing rate
-    subplot(4,6,12+(i-1)*12);
-    psthE=nanmean(R_2R.Ev(RD).PSTHz(Sel,Ushow),1);
-    semE=nanste(R_2R.Ev(RD).PSTHz(Sel,Ushow),1); %calculate standard error of the mean
-    upE=psthE+semE;
-    downE=psthE-semE;
-    hold on;
-    patch([time2,time2(end:-1:1)],[upE,downE(end:-1:1)],exc,'EdgeColor','none');alpha(0.5);
-    plot(time2,psthE,'Color',exc,'linewidth',1);
-    plot([-1 2],[0 0],':','color','k','linewidth',0.75);
-    plot([0 0],[-3 9],':','color','k','linewidth',0.75);
-    axis([-1 2 min(downE)-0.2 max(upE)+0.3]);
-    title(['RD excitations (' num2str(round(100*sum(Sel)/sum(Reg))) '%)'])
-    ylabel('Z-score');
-end
-
-%% plotting firing rate on sucrose and maltodextrin trials
-figure;
-
-%color map
-[magma,inferno,plasma,viridis]=colormaps;
-colormap(plasma);
-c=[-100 2000];ClimE=sign(c).*abs(c).^(1/4);%colormap
-
-Xaxis=[-1 4];
-Ishow=find(R_2R.Param.Tm>=Xaxis(1) & R_2R.Param.Tm<=Xaxis(2));
-time1=R_2R.Param.Tm(Ishow);
-Xaxis2=[-1 4];
-Ushow=find(R_2R.Param.Tm>=Xaxis2(1) & R_2R.Param.Tm<=Xaxis2(2));
-time2=R_2R.Param.Tm(Ushow);
-
-sucrose=[1  0.6  0.1];
-maltodextrin=[.9  0.3  .9];
-water=[0.00  0.75  0.75];
-total=[0.3  0.1  0.8];
-exc=[0 113/255 188/255];
-inh=[240/255 0 50/255];
-
-RD=strcmp('RD', R_2R.Erefnames);
-RD1=strcmp('RD1', R_2R.Erefnames);
-RD2=strcmp('RD2', R_2R.Erefnames);
-
-for i=1:2
-    
-    if i==1 Reg=NAneurons; end %1 is NAc
-    if i==2 Reg=VPneurons; end %2 is VP
-
-    %% Plotting heatmaps of sucrose and maltodextrin response
-    Sel=Reg&R_2R.Ev(RD).RespDir<2; %gets all neurons
-
-    %sort RD1 heatmap by magnitude of response
-    Neurons=R_2R.Ev(RD1).PSTHz(Sel,Ushow); %get the firing rates of neurons of interest
-    TMP=R_2R.Ev(RD).Meanz(Sel); %find the magnitude of the inhibitions for this event
-    TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
-    [~,SORTimg]=sort(TMP);
-    Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
-
-    %heatmap
-    subplot(2,40,[1:9]+(i-1)*40);
-    imagesc(time2,[1,sum(Sel,1)],Neurons,ClimE); title('Sucrose trials');
-    ylabel('Individual neurons sorted by RD response');
-    xlabel('Seconds from RD');
-    hold on;
-    plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
-
-    %sort RD2 heatmap the same way
-    Neurons=R_2R.Ev(RD2).PSTHz(Sel,Ushow); %get the firing rates of neurons of interest
-    [~,SORTimg]=sort(TMP);
-    Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
-    xlabel('Seconds from RD');
-
-
-    %heatmap
-    subplot(2,40,[13:21]+(i-1)*40);
-    imagesc(time2,[1,sum(Sel,1)],Neurons,ClimE); title('Maltodextrin trials');
-    xlabel('Seconds from RD');
-    hold on;
-    plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
-
-    %% Plotting average firing rate of reward excitations
-    Sel=Reg&(R_2R.Ev(RD).RespDir>0); %| R_2R.Ev(LickR1).RespDir>0 | R_2R.Ev(LickR2).RespDir>0); %Find neurons that have an inhibitory response to this event
-
-    %average firing rate to sucrose
-    psth1=nanmean(R_2R.Ev(RD1).PSTHz(Sel,Ishow),1);
-    sem1=nanste(R_2R.Ev(RD1).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
-    up1=psth1+sem1;
-    down1=psth1-sem1;
-
-    %average firing rate to maltodextrin
-    psth2=nanmean(R_2R.Ev(RD2).PSTHz(Sel,Ishow),1);
-    sem2=nanste(R_2R.Ev(RD2).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
-    up2=psth2+sem2;
-    down2=psth2-sem2;
-
-    %plotting
-    subplot(4,5,[9 10]+(i-1)*10);
-    hold on;
-    plot(time1,psth1,'Color',sucrose,'linewidth',1);
-    plot(time1,psth2,'Color',maltodextrin,'linewidth',1);
-    patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
-    patch([time1,time1(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
-    plot([-1 4],[0 0],':','color','k','linewidth',0.75);
-    plot([0 0],[-3 9],':','color','k','linewidth',0.75);
-    axis([-1 4 min(down1)-0.2 max(up1)+0.3]);
-    title(['Reward-excited (n=' num2str(sum(Sel)) ')'])
-    ylabel('Mean firing (z-score)');
-    xlabel('Seconds from RD');
-    legend('Suc trials','Mal trials');
-
-    % %onset
-    % subplot(2,4,4);
-    % cumul(R_2R.Ev(38).Onsets(Sel,1),[R_2R.Param.RespWin(1,1):0.04:R_2R.Param.RespWin(38,2)],2,maltodextrin);
-
-
-    %% Plotting average firing rate of reward inhibitions
-    Sel=Reg&(R_2R.Ev(RD).RespDir<0); %| R_2R.Ev(LickR1).RespDir<0 | R_2R.Ev(LickR2).RespDir<0); %Find neurons that have an inhibitory response to this event
-
-    %avg firing rate to sucrose
-    psth1=nanmean(R_2R.Ev(12).PSTHz(Sel,Ushow),1);
-    sem1=nanste(R_2R.Ev(12).PSTHz(Sel,Ushow),1); %calculate standard error of the mean
-    up1=psth1+sem1;
-    down1=psth1-sem1;
-
-    %avg firing rate to maltodextrin
-    psth2=nanmean(R_2R.Ev(13).PSTHz(Sel,Ushow),1);
-    sem2=nanste(R_2R.Ev(13).PSTHz(Sel,Ushow),1); %calculate standard error of the mean
-    up2=psth2+sem2;
-    down2=psth2-sem2;
-
-    %plotting
-    subplot(4,5,[4 5]+(i-1)*10);
-    hold on;
-    plot(time2,psth1,'Color',sucrose,'linewidth',1);
-    plot(time2,psth2,'Color',maltodextrin,'linewidth',1);
-    patch([time2,time2(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
-    patch([time2,time2(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
-    plot([-1 4],[0 0],':','color','k','linewidth',0.75);
-    plot([0 0],[-3 9],':','color','k','linewidth',0.75);
-    axis([-1 4 min(down1)-0.2 max(up1)+0.2]);
-    title(['Reward-inhibited (n=' num2str(sum(Sel)) ')'])
-    ylabel('Mean firing (z-score)');
-    xlabel('Seconds from RD');
-    legend('Suc trials','Mal trials','location','southeast');
-end