Browse Source

Upload files to 'MatlabScripts'

David Ottenheimer 10 months ago
parent
commit
00e1508a89
1 changed files with 866 additions and 0 deletions
  1. 866 0
      MatlabScripts/d_Fig2S2S3S5EventsBins.m

+ 866 - 0
MatlabScripts/d_Fig2S2S3S5EventsBins.m

@@ -0,0 +1,866 @@
+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];
+NAc=[0.5  0.1  0.8];
+VP=[0.3  0.7  0.1];
+
+%% 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(cumsel(NAneurons,bins-1)),'Color', NAc,'linewidth',1.5);
+plot(xaxis,sum(cumsel(VPneurons,2:bins-1),1)/sum(cumsel(VPneurons,bins-1)),'Color', VP,'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 1]);
+legend('NAc','VP','location','southeast');
+xlabel('Seconds from RD');
+ylabel('Cumulative distribution');
+title('Reward-selective onsets');
+
+%statistical test comparing onset of reward-selective responses
+for i=1:length(cumsel)
+    if sum(cumsel(i,:))>0
+        onsetcolumn=min(find(cumsel(i,:)));
+        onset(i,1)=round((onsetcolumn-1)*binint+binstart+BinDura(2)/2,1);
+    else
+        onset(i,1)=NaN;
+    end
+end        
+[~,R_2R.OnsetSig{1,1},R_2R.OnsetSig{1,2}]=anovan(onset,{NAneurons,R_2R.Ninfo(:,4)},'varnames',{'Region','Subject'},'random',[2],'nested',[0 0;1 0],'display','off','model','full');
+if cell2mat(R_2R.OnsetSig{1,1}(2,7))<0.05
+    plot((xaxis(end)-xaxis(1))/2+xaxis(1),0.95,'*','color','k','linewidth',1);
+end
+
+%% 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;
+R_2R.Param.SelectiveBins=[Bin1 Bin2];
+
+
+% 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(RD1).PSTHz(Sel,Ushow),1);
+    sem1=nanste(R_2R.Ev(RD1).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(RD2).PSTHz(Sel,Ushow),1);
+    sem2=nanste(R_2R.Ev(RD2).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