Browse Source

Upload files to 'MatlabScripts'

David Ottenheimer 5 years ago
parent
commit
ac93b930da

+ 50 - 0
MatlabScripts/a_NexExtract.m

@@ -0,0 +1,50 @@
+%% create a struct with timestamps of all neurons and events
+
+%before starting:
+%put all nex files in the same folder
+%name of each file starts with "NA1" or corresponding region and rat
+%for water and 3 rewards sessions, "WA1" or "TH1" (all VP)
+
+%need supporting programs "readNexFile" and "myfind"
+
+clear all;clc 
+SAVE_FLAG=1;
+tic
+
+address=['C:\Users\dottenh2\Documents\MATLAB\David\2Rewards Nex Files'];
+
+AF=dir([address,'\\*.nex']);
+Sesnum=0;
+for k=1:length(AF)
+      
+        fname=AF(k).name;
+        [nexfile] = readNexFile([address,'\\',fname]);  fprintf([fname,'\n'])
+        Iind=myfind(nexfile.intervals,'AllFile'); 
+        Sesnum=Sesnum+1;
+        % Get all events timestamps for the selected session
+        NUM_EVENTS=length(nexfile.events);
+        for j=1:NUM_EVENTS
+            evt=nexfile.events{j}.timestamps;
+            RAW(Sesnum).Erast{j,1}=evt(find((evt>nexfile.intervals{Iind}.intStarts) & (evt<nexfile.intervals{Iind}.intEnds)));
+            RAW(Sesnum).Einfo(j,:)={fname,nexfile.events{j}.name};
+        end
+        % start from neuron 1 get ready for the next session
+        
+        % Get the Neuron timestamps and waveforms
+       NUM_NEURONS=length(nexfile.neurons);
+        for j=1:NUM_NEURONS
+            nrn=nexfile.neurons{j}.timestamps;
+            RAW(Sesnum).Nrast{j,1}=nrn(find((nrn>nexfile.intervals{Iind}.intStarts) & (nrn<nexfile.intervals{Iind}.intEnds)));
+            RAW(Sesnum).Ninfo(j,:)={fname,nexfile.neurons{j}.name};
+        end
+        
+        RAW(Sesnum).Type=fname(1:3);
+    
+end
+toc   
+%% SAVING DATA
+if SAVE_FLAG
+    save('RAW.mat','RAW')
+end
+fprintf('\n')
+toc

+ 682 - 0
MatlabScripts/b_EventAnalysis.m

@@ -0,0 +1,682 @@
+%% initial analysis of event-evoked activity
+
+%need to save excel sheets in the specified folder
+%the excel sheets say which events to make PSTHs for and what windows to
+%use for statistical testing
+
+%clear all; clc;
+global Dura Baseline Tm Tbase BSIZE Tbin
+tic
+
+if ~exist('RAW'), load ('RAW.mat'); end
+
+%Main settings
+SAVE_FLAG=1;
+BSIZE=0.01; %Do not change
+Dura=[-5 12]; Tm=Dura(1):BSIZE:Dura(2);
+%Baseline=[-22 0]; Tbase=Baseline(1):BSIZE:Baseline(2); %now defined line 98
+Tbin=-0.5:0.005:0.5; %window used to determine the optimal binsize
+PStat=0.05; %for comparing pre versus post windows, or event A versus event B
+MinNumTrials=5; %how many trials of event there need to be to conduct analysis
+BinSize=0.05;
+
+%Smoothing
+Smoothing=1; %0 for raw and 1 for smoothing
+SmoothTYPE='lowess';
+SmoothSPAN=5;
+if Smoothing~=1, SmoothTYPE='NoSmoothing';SmoothSPAN=NaN; end
+
+%for bins analysis
+BinBase=[-22 -12]; %For normalizing activity
+BinDura=[0 0.6]; %size of bin
+bins = 66; %number of bins
+binint = 0.1; %spacing of bins
+binstart = -2; %start time of first bin relative to event
+
+%%  Standard sucrose and maltodextrin sessions
+
+%start fresh
+R_2R=[];R_2R.Ninfo={};NN=0;Nneurons=0;
+
+% List of events to analyze and analysis windows EXTRACTED from excel file
+path='C:\Users\dottenh2\Documents\MATLAB\David\2Rewards Nex Files\2R_paper.xls';
+[~,Erefnames]=xlsread(path,'Windows','a3:a15'); % cell that contains the event names
+prewin =  xlsread(path,'Windows','b3:c15');
+postwin = xlsread(path,'Windows','d3:e15');
+
+%saves event names for reference later
+R_2R.Erefnames=Erefnames;
+
+%Finds the total number of neurons in 2R and marks them by region/session
+for i=1:length(RAW)
+    if strcmp('NA',RAW(i).Type(1:2)) | strcmp('VP',RAW(i).Type(1:2))
+        R_2R.Ninfo=cat(1,R_2R.Ninfo,RAW(i).Ninfo);
+        Nneurons=Nneurons+size(RAW(i).Nrast,1);
+    end
+end
+for i=1:Nneurons
+    Session=string(R_2R.Ninfo(i,1));
+    Name=char(Session);
+    R_2R.Ninfo(i,3)=cellstr(Name(1:2));
+    R_2R.Ninfo(i,4)=cellstr(Name(1:3));
+end
+
+% preallocating
+R_2R.Param.Tm=Tm;
+R_2R.Param.Tbin=Tbin;
+R_2R.Param.Dura=Dura;
+R_2R.Param.Baseline=Baseline;
+R_2R.Param.PStat=PStat;
+R_2R.Param.MinNumTrials=MinNumTrials;
+R_2R.Param.path=path;
+R_2R.Param.prewin=prewin;
+R_2R.Param.postwin=postwin;
+R_2R.Param.SmoothTYPE=SmoothTYPE;
+R_2R.Param.SmoothSPAN=SmoothSPAN;
+R_2R.Param.BinBase=BinBase;
+R_2R.Param.BinDura=BinDura;
+R_2R.Param.bins = bins;
+R_2R.Param.binint = binint;
+R_2R.Param.binstart = binstart;
+
+for k=1:length(Erefnames)
+    R_2R.Ev(k).PSTHraw(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
+    R_2R.Ev(k).PSTHz(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
+    R_2R.Ev(k).Meanraw(1:Nneurons,1)=NaN;
+    R_2R.Ev(k).Meanz(1:Nneurons,1)=NaN;
+    R_2R.Ev(k).BW(1:Nneurons,1)=NaN;
+    R_2R.Ev(k).signrank(1:Nneurons,1)=NaN;
+    R_2R.Ev(k).RespDir(1:Nneurons,1)=NaN;
+    R_2R.Ev(k).NumberTrials(1:Nneurons,1)=NaN;
+end
+
+% runs the main routine
+
+for i=1:length(RAW) %loops through sessions
+    if strcmp('NA',RAW(i).Type(1:2)) | strcmp('VP',RAW(i).Type(1:2))
+        for j= 1:size(RAW(i).Nrast,1) %Number of neurons per session
+            NN=NN+1; %neuron counter
+            for k=1:length(Erefnames) %loops thorough the events
+                EvInd=strcmp(Erefnames(k),RAW(i).Einfo(:,2)); %find the event id number from RAW
+
+                if sum(EvInd)==0
+                    fprintf('HOWDY, CANT FIND EVENTS FOR ''%s''\n',Erefnames{k});
+                end
+
+                R_2R.Ev(k).NumberTrials(NN,1)=length(RAW(i).Erast{EvInd});
+                Tbase=prewin(k,1):BSIZE:prewin(k,2);
+                if  ~isempty(EvInd) && R_2R.Ev(k).NumberTrials(NN,1)>=MinNumTrials %avoid analyzing sessions where that do not have enough trials
+                    [PSR0,N0]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},prewin(k,:),{1});% makes collapsed rasters for baseline for use in normalizing
+                    [PSR1,N1]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{1});% makes collpased rasters. PSR1 is a cell(neurons)
+                    [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
+                    [PSR3,N0]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},prewin(k,:),{2});% makes trial by trial rasters for baseline for use in response detection
+                    if ~isempty(PSR0{1}) || ~isempty(PSR1{1}) %to avoid errors, added on 12/28 2011
+                        %Optimal bin size
+                        %[PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Nrast{j},1),1),{2,1});%-----DP used here
+                        %[PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(i).Nrast{j},1),1),{1,0,BW1});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
+                        %Fixed bin size
+                        [PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Nrast{j},1),1),{2,0,BinSize});%-----Fixed bin used here
+                        [PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(i).Nrast{j},1),1),{1,0,BinSize});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
+                        PTH1=smooth(PTH1,SmoothSPAN,SmoothTYPE)';
+                        PTH0=smooth(PTH0,SmoothSPAN,SmoothTYPE)';
+
+                        %------------- Fills the R_2R.Ev(k) fields --------------
+                        R_2R.Ev(k).BW(NN,1)=BW1;
+                        R_2R.Ev(k).PSTHraw(NN,1:length(Tm))=PTH1;
+                        R_2R.Ev(k).Meanraw(NN,1)=nanmean(R_2R.Ev(k).PSTHraw(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
+                        if sum(PTH0,2)~=0
+                            R_2R.Ev(k).PSTHz(NN,1:length(Tm))=normalize(PTH1,PTH0,0);
+                            R_2R.Ev(k).Meanz(NN,1)=nanmean(R_2R.Ev(k).PSTHz(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
+                        else
+                            R_2R.Ev(k).PSTHz(NN,1:length(Tm))=NaN(1,length(Tm));
+                            R_2R.Ev(k).Meanz(NN,1)=NaN;
+                        end
+
+                        %------------------ firing (in Hz) per trial in pre/post windows ------------------
+                        %used to make the between events comparisons and Response detection in a single window----
+                        ev(k).pre=NaN(size(RAW(i).Erast{EvInd},1),1);
+                        ev(k).post=NaN(size(RAW(i).Erast{EvInd},1),1);
+                        for m=1:size(RAW(i).Erast{EvInd},1) %loops through trials
+                            ev(k).pre(m)=sum(PSR3{m}<prewin(k,2) & PSR3{m}>prewin(k,1))/(prewin(k,2)-prewin(k,1)); %CHANGED FROM PSR2 to PSR0 here 10/24/17
+                            ev(k).post(m)=sum(PSR2{m}<postwin(k,2) & PSR2{m}>postwin(k,1))/(postwin(k,2)-postwin(k,1));
+                        end
+
+                        %-------------------- signrank on event and direction----
+                        [R_2R.Ev(k).signrank(NN,1),~]=signrank(ev(k).pre, ev(k).post); %Signrank used here because it is a dependant sample test
+                        if R_2R.Ev(k).signrank(NN,1)<PStat
+                            R_2R.Ev(k).RespDir(NN,1)=sign(mean(ev(k).post)-mean(ev(k).pre));
+                        else R_2R.Ev(k).RespDir(NN,1)=0;
+                        end
+
+                    end %if ~isempty(PSR0{1}) || ~isempty(PSR1{1})
+                end %if EvInd=0 OR n(trials) < MinNumTrials fills with NaN
+            end %Events
+
+
+            %----------------------- Check for difference in cue responding after sucrose and maltodextrin trials -----------------------         
+            CueP1=strcmp('CueP1', Erefnames);
+            CueP2=strcmp('CueP2', Erefnames);
+            [R_2R.CueP12Stat(NN,1),~]=ranksum(ev(CueP1).post,ev(CueP2).post); %Ranksum test used becasue it is an independent sample test
+
+
+            %----------------------- Check for difference in PE responding after sucrose and maltodextrin trials -----------------------         
+            PEP1=strcmp('PEP1', Erefnames);
+            PEP2=strcmp('PEP2', Erefnames);
+            [R_2R.PEP12Stat(NN,1),~]=ranksum(ev(PEP1).post,ev(PEP1).post); %Ranksum test used becasue it is an independent sample test
+
+            %----------------------- Check for reward modulation across all the bins surrounding reward delivery-----------------------
+            BHz=0;
+            Trial=0;
+            Reward=0;
+
+            EV1=strmatch('RD1',RAW(i).Einfo(:,2),'exact');
+            EV2=strmatch('RD2',RAW(i).Einfo(:,2),'exact');
+
+            %get mean baseline firing for all reward trials
+            [Bcell1,B1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV1},BinBase,{2});% makes trial by trial rasters for baseline
+            for y= 1:B1n
+                BHz(1,y)=sum(Bcell1{1,y}>BinBase(1))/(BinBase(1,2)-BinBase(1,1));
+                Reward(1,y)=1;
+            end
+            [Bcell2,B2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV2},BinBase,{2});% makes trial by trial rasters for baseline
+            for z= 1:B2n
+                BHz(1,z+B1n)=sum(Bcell2{1,z}>BinBase(1))/(BinBase(1,2)-BinBase(1,1));
+                Reward(1,z+B1n)=2;
+            end
+
+            for k=1:bins %go through each bin for each neuron
+                EVHz=0;
+                %get trial by trial firing rate for sucrose trials
+                [EV1cell,EV1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV1},[(BinDura(1)+(binstart - binint)+k*binint) (BinDura(2)+(binstart - binint)+k*binint)],{2});% makes trial by trial rasters for event
+                for y= 1:EV1n
+                    EVHz(1,y)=sum(EV1cell{1,y}>binstart)/(BinDura(1,2)-BinDura(1,1));
+                end       
+
+                %get trial by trial firing rate for maltodextrin trials
+                [EV2cell,EV2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV2},[(BinDura(1)+(binstart - binint)+k*binint) (BinDura(2)+(binstart - binint)+k*binint)],{2});% makes trial by trial rasters for event
+                for z= 1:EV2n
+                    EVHz(1,z+EV1n)=sum(EV2cell{1,z}>binstart)/(BinDura(1,2)-BinDura(1,1));
+                end       
+
+
+                %run ANOVA on the effects of pre vs post, cued or not,
+                %kind of reward, and the interactions
+                [~,Bintbl,~]=anovan([BHz EVHz],{cat(2,zeros(1,length(Reward)),ones(1,length(Reward))),[Reward Reward]},'varnames',{'Pre vs Post','Reward'},'model','full','display','off');
+                R_2R.BinStatRwrd{k,1}(NN).IntSig=cell2mat(Bintbl(4,7));
+
+                %get mean firing rates
+                R_2R.BinStatRwrd{k,1}(NN).R1Mean=(nanmean(EVHz(Reward==1))-nanmean(BHz(Reward==1)))/std(BHz(Reward==1));
+                R_2R.BinStatRwrd{k,1}(NN).R2Mean=(nanmean(EVHz(Reward==2))-nanmean(BHz(Reward==2)))/std(BHz(Reward==2));
+
+                %signranks for checking if inhibitory or excitatory
+
+                %sucrose trials
+                signranksig=signrank(BHz(Reward==1), EVHz(Reward==1)); %Signrank used here because it is a dependant sample test
+                        if signranksig<PStat
+                            R_2R.BinStatRwrd{k,1}(NN).SucRespDir=sign(nanmean(EVHz(Reward==1))-mean(BHz(Reward==1)));
+                        else
+                            R_2R.BinStatRwrd{k,1}(NN).SucRespDir=0;
+                        end
+
+                %maltodextrin trials
+                signranksig=signrank(BHz(Reward==2), EVHz(Reward==2)); %Signrank used here because it is a dependant sample test
+                        if signranksig<PStat
+                            R_2R.BinStatRwrd{k,1}(NN).MalRespDir=sign(nanmean(EVHz(Reward==2))-mean(BHz(Reward==2)));
+                        else
+                            R_2R.BinStatRwrd{k,1}(NN).MalRespDir=0;
+                        end
+
+            end %bins for bins analysis
+
+
+            fprintf('2R Neuron #%d\n',NN);
+        end %neurons: FOR j= 1:size(RAW(i).Nrast,1)
+    end %if it's the right kind of session
+end %sessions: FOR i=1:length(RAW)
+
+if SAVE_FLAG
+    save('R_2R.mat','R_2R')
+end
+
+toc
+
+
+%%  water and maltodextrin sessions
+
+%start fresh
+R_W=[];R_W.Ninfo={};NN=0;Nneurons=0;
+
+% List of events to analyze and analysis windows EXTRACTED from excel file
+path='C:\Users\dottenh2\Documents\MATLAB\David\2Rewards Nex Files\2R_paper.xls';
+[~,Erefnames]=xlsread(path,'Windows','a3:a15'); % cell that contains the event names
+prewin =  xlsread(path,'Windows','b3:c15');
+postwin = xlsread(path,'Windows','d3:e15');
+
+%saves event names for reference later
+R_W.Erefnames=Erefnames;
+
+%Finds the total number of neurons in 2R and marks them by region/session
+for i=1:length(RAW)
+    if strcmp('WA',RAW(i).Type(1:2))
+        R_W.Ninfo=cat(1,R_W.Ninfo,RAW(i).Ninfo);
+        Nneurons=Nneurons+size(RAW(i).Nrast,1);
+    end
+end
+for i=1:Nneurons
+    Session=string(R_W.Ninfo(i,1));
+    Name=char(Session);
+    R_W.Ninfo(i,3)=cellstr(Name(1:2));
+    R_W.Ninfo(i,4)=cellstr(Name(1:3));
+end
+
+% preallocating
+R_W.Param.Tm=Tm;
+R_W.Param.Tbin=Tbin;
+R_W.Param.Dura=Dura;
+R_W.Param.Baseline=Baseline;
+R_W.Param.PStat=PStat;
+R_W.Param.MinNumTrials=MinNumTrials;
+R_W.Param.path=path;
+R_W.Param.prewin=prewin;
+R_W.Param.postwin=postwin;
+R_W.Param.SmoothTYPE=SmoothTYPE;
+R_W.Param.SmoothSPAN=SmoothSPAN;
+R_W.Param.BinBase=BinBase;
+R_W.Param.BinDura=BinDura;
+R_W.Param.bins = bins;
+R_W.Param.binint = binint;
+R_W.Param.binstart = binstart;
+
+for k=1:length(Erefnames)
+    R_W.Ev(k).PSTHraw(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
+    R_W.Ev(k).PSTHz(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
+    R_W.Ev(k).Meanraw(1:Nneurons,1)=NaN;
+    R_W.Ev(k).Meanz(1:Nneurons,1)=NaN;
+    R_W.Ev(k).BW(1:Nneurons,1)=NaN;
+    R_W.Ev(k).signrank(1:Nneurons,1)=NaN;
+    R_W.Ev(k).RespDir(1:Nneurons,1)=NaN;
+    R_W.Ev(k).NumberTrials(1:Nneurons,1)=NaN;
+end
+
+% runs the main routine
+
+for i=1:length(RAW) %loops through sessions
+    if strcmp('WA',RAW(i).Type(1:2))
+        for j= 1:size(RAW(i).Nrast,1) %Number of neurons per session
+            NN=NN+1; %neuron counter
+            for k=1:length(Erefnames) %loops thorough the events
+                EvInd=strcmp(Erefnames(k),RAW(i).Einfo(:,2)); %find the event id number from RAW
+
+                if sum(EvInd)==0
+                    fprintf('HOWDY, CANT FIND EVENTS FOR ''%s''\n',Erefnames{k});
+                end
+
+                R_W.Ev(k).NumberTrials(NN,1)=length(RAW(i).Erast{EvInd});
+                Tbase=prewin(k,1):BSIZE:prewin(k,2);
+                if  ~isempty(EvInd) && R_W.Ev(k).NumberTrials(NN,1)>=MinNumTrials %avoid analyzing sessions where that do not have enough trials
+                    [PSR0,N0]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},prewin(k,:),{1});% makes collapsed rasters for baseline for use in normalizing
+                    [PSR1,N1]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{1});% makes collpased rasters. PSR1 is a cell(neurons)
+                    [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
+                    [PSR3,N0]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},prewin(k,:),{2});% makes trial by trial rasters for baseline for use in response detection
+                    if ~isempty(PSR0{1}) || ~isempty(PSR1{1}) %to avoid errors, added on 12/28 2011
+                        %optimal bin size
+                        %[PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Nrast{j},1),1),{2,1});%-----DP used here
+                        %[PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(i).Nrast{j},1),1),{1,0,BW1});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
+                        %Fixed bin size
+                        [PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Nrast{j},1),1),{2,0,BinSize});%-----fixed bin size
+                        [PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(i).Nrast{j},1),1),{1,0,BinSize});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
+                        PTH1=smooth(PTH1,SmoothSPAN,SmoothTYPE)';
+                        PTH0=smooth(PTH0,SmoothSPAN,SmoothTYPE)';
+
+                        %------------- Fills the R_W.Ev(k) fields --------------
+                        R_W.Ev(k).BW(NN,1)=BW1;
+                        R_W.Ev(k).PSTHraw(NN,1:length(Tm))=PTH1;
+                        R_W.Ev(k).Meanraw(NN,1)=nanmean(R_W.Ev(k).PSTHraw(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
+                        if sum(PTH0,2)~=0
+                            R_W.Ev(k).PSTHz(NN,1:length(Tm))=normalize(PTH1,PTH0,0);
+                            R_W.Ev(k).Meanz(NN,1)=nanmean(R_W.Ev(k).PSTHz(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
+                        else
+                            R_W.Ev(k).PSTHz(NN,1:length(Tm))=NaN(1,length(Tm));
+                            R_W.Ev(k).Meanz(NN,1)=NaN;
+                        end
+
+                        %------------------ firing (in Hz) per trial in pre/post windows ------------------
+                        %used to make the between events comparisons and Response detection in a single window----
+                        ev(k).pre=NaN(size(RAW(i).Erast{EvInd},1),1);
+                        ev(k).post=NaN(size(RAW(i).Erast{EvInd},1),1);
+                        for m=1:size(RAW(i).Erast{EvInd},1) %loops through trials
+                            ev(k).pre(m)=sum(PSR3{m}<prewin(k,2) & PSR3{m}>prewin(k,1))/(prewin(k,2)-prewin(k,1)); %CHANGED FROM PSR2 to PSR0 here 10/24/17
+                            ev(k).post(m)=sum(PSR2{m}<postwin(k,2) & PSR2{m}>postwin(k,1))/(postwin(k,2)-postwin(k,1));
+                        end
+
+                        %-------------------- signrank on event and direction----
+                        [R_W.Ev(k).signrank(NN,1),~]=signrank(ev(k).pre, ev(k).post); %Signrank used here because it is a dependant sample test
+                        if R_W.Ev(k).signrank(NN,1)<PStat
+                            R_W.Ev(k).RespDir(NN,1)=sign(mean(ev(k).post)-mean(ev(k).pre));
+                        else R_W.Ev(k).RespDir(NN,1)=0;
+                        end
+
+                    end %if ~isempty(PSR0{1}) || ~isempty(PSR1{1})
+                end %if EvInd=0 OR n(trials) < MinNumTrials fills with NaN
+            end %Events
+
+
+            %----------------------- Check for reward modulation across all the bins surrounding reward delivery-----------------------
+            BHz=0;
+            Trial=0;
+            Reward=0;
+
+            EV1=strmatch('RD1',RAW(i).Einfo(:,2),'exact');
+            EV2=strmatch('RD2',RAW(i).Einfo(:,2),'exact');
+
+            %get mean baseline firing for all reward trials
+            [Bcell1,B1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV1},BinBase,{2});% makes trial by trial rasters for baseline
+            for y= 1:B1n
+                BHz(1,y)=sum(Bcell1{1,y}>BinBase(1))/(BinBase(1,2)-BinBase(1,1));
+                Reward(1,y)=1;
+            end
+            [Bcell2,B2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV2},BinBase,{2});% makes trial by trial rasters for baseline
+            for z= 1:B2n
+                BHz(1,z+B1n)=sum(Bcell2{1,z}>BinBase(1))/(BinBase(1,2)-BinBase(1,1));
+                Reward(1,z+B1n)=2;
+            end
+
+            for k=1:bins %go through each bin for each neuron
+                EVHz=0;
+                %get trial by trial firing rate for water trials
+                [EV1cell,EV1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV1},[(BinDura(1)+(binstart - binint)+k*binint) (BinDura(2)+(binstart - binint)+k*binint)],{2});% makes trial by trial rasters for event
+                for y= 1:EV1n
+                    EVHz(1,y)=sum(EV1cell{1,y}>binstart)/(BinDura(1,2)-BinDura(1,1));
+                end       
+
+                %get trial by trial firing rate for maltodextrin trials
+                [EV2cell,EV2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV2},[(BinDura(1)+(binstart - binint)+k*binint) (BinDura(2)+(binstart - binint)+k*binint)],{2});% makes trial by trial rasters for event
+                for z= 1:EV2n
+                    EVHz(1,z+EV1n)=sum(EV2cell{1,z}>binstart)/(BinDura(1,2)-BinDura(1,1));
+                end       
+
+
+                %run ANOVA on the effects of pre vs post, cued or not,
+                %kind of reward, and the interactions
+                [~,Bintbl,~]=anovan([BHz EVHz],{cat(2,zeros(1,length(Reward)),ones(1,length(Reward))),[Reward Reward]},'varnames',{'Pre vs Post','Reward'},'model','full','display','off');
+                R_W.BinStatRwrd{k,1}(NN).IntSig=cell2mat(Bintbl(4,7));
+
+                %get mean firing rates
+                R_W.BinStatRwrd{k,1}(NN).R1Mean=(nanmean(EVHz(Reward==1))-nanmean(BHz(Reward==1)))/std(BHz(Reward==1));
+                R_W.BinStatRwrd{k,1}(NN).R2Mean=(nanmean(EVHz(Reward==2))-nanmean(BHz(Reward==2)))/std(BHz(Reward==2));
+
+                %signranks for checking if inhibitory or excitatory
+
+                %water trials
+                signranksig=signrank(BHz(Reward==1), EVHz(Reward==1)); %Signrank used here because it is a dependant sample test
+                if signranksig<PStat
+                    R_W.BinStatRwrd{k,1}(NN).WatRespDir=sign(nanmean(EVHz(Reward==1))-mean(BHz(Reward==1)));
+                else
+                    R_W.BinStatRwrd{k,1}(NN).WatRespDir=0;
+                end
+
+                %maltodextrin trials
+                signranksig=signrank(BHz(Reward==2), EVHz(Reward==2)); %Signrank used here because it is a dependant sample test
+                if signranksig<PStat
+                    R_W.BinStatRwrd{k,1}(NN).MalRespDir=sign(nanmean(EVHz(Reward==2))-mean(BHz(Reward==2)));
+                else
+                    R_W.BinStatRwrd{k,1}(NN).MalRespDir=0;
+                end
+
+            end %bins for bins analysis
+
+
+            fprintf('Water neuron #%d\n',NN);
+        end %neurons: FOR j= 1:size(RAW(i).Nrast,1)
+    end %if right kind of session
+end %sessions: FOR i=1:length(RAW)
+
+if SAVE_FLAG
+    save('R_W.mat','R_W')
+end
+
+toc
+
+%%  All three rewards sessions
+
+%make this lower because there are only 3 of some current/previous reward combinations
+MinNumTrials=3; %how many trials of event there need to be to conduct analysis
+
+
+%start fresh
+R_3R=[];R_3R.Ninfo={};NN=0;Nneurons=0;
+
+% List of events to analyze and analysis windows EXTRACTED from excel file
+path='C:\Users\dottenh2\Documents\MATLAB\David\2Rewards Nex Files\3R_paper.xls';
+[~,Erefnames]=xlsread(path,'Windows','a3:a23'); % cell that contains the event names
+prewin =  xlsread(path,'Windows','b3:c23');
+postwin = xlsread(path,'Windows','d3:e23');
+
+%saves event names for reference later
+R_3R.Erefnames=Erefnames;
+
+%Finds the total number of neurons in 2R and marks them by region/session
+for i=1:length(RAW)
+    if strcmp('TH',RAW(i).Type(1:2))
+        R_3R.Ninfo=cat(1,R_3R.Ninfo,RAW(i).Ninfo);
+        Nneurons=Nneurons+size(RAW(i).Nrast,1);
+    end
+end
+for i=1:Nneurons
+    Session=string(R_3R.Ninfo(i,1));
+    Name=char(Session);
+    R_3R.Ninfo(i,3)=cellstr(Name(1:2));
+    R_3R.Ninfo(i,4)=cellstr(Name(1:3));
+end
+
+% preallocating
+R_3R.Param.Tm=Tm;
+R_3R.Param.Tbin=Tbin;
+R_3R.Param.Dura=Dura;
+R_3R.Param.Baseline=Baseline;
+R_3R.Param.PStat=PStat;
+R_3R.Param.MinNumTrials=MinNumTrials;
+R_3R.Param.path=path;
+R_3R.Param.prewin=prewin;
+R_3R.Param.postwin=postwin;
+R_3R.Param.SmoothTYPE=SmoothTYPE;
+R_3R.Param.SmoothSPAN=SmoothSPAN;
+R_3R.Param.BinBase=BinBase;
+R_3R.Param.BinDura=BinDura;
+R_3R.Param.bins = bins;
+R_3R.Param.binint = binint;
+R_3R.Param.binstart = binstart;
+
+for k=1:length(Erefnames)
+    R_3R.Ev(k).PSTHraw(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
+    R_3R.Ev(k).PSTHz(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
+    R_3R.Ev(k).Meanraw(1:Nneurons,1)=NaN;
+    R_3R.Ev(k).Meanz(1:Nneurons,1)=NaN;
+    R_3R.Ev(k).BW(1:Nneurons,1)=NaN;
+    R_3R.Ev(k).signrank(1:Nneurons,1)=NaN;
+    R_3R.Ev(k).RespDir(1:Nneurons,1)=NaN;
+    R_3R.Ev(k).NumberTrials(1:Nneurons,1)=NaN;
+end
+
+% runs the main routine
+
+for i=1:length(RAW) %loops through sessions
+    if strcmp('TH',RAW(i).Type(1:2))
+        for j= 1:size(RAW(i).Nrast,1) %Number of neurons per session
+            NN=NN+1; %neuron counter
+            for k=1:length(Erefnames) %loops thorough the events
+                EvInd=strcmp(Erefnames(k),RAW(i).Einfo(:,2)); %find the event id number from RAW
+
+                if sum(EvInd)==0
+                    fprintf('HOWDY, CANT FIND EVENTS FOR ''%s''\n',Erefnames{k});
+                end
+
+                R_3R.Ev(k).NumberTrials(NN,1)=length(RAW(i).Erast{EvInd});
+                Tbase=prewin(k,1):BSIZE:prewin(k,2);
+                if  ~isempty(EvInd) && R_3R.Ev(k).NumberTrials(NN,1)>=MinNumTrials %avoid analyzing sessions where that do not have enough trials
+                    [PSR0,N0]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},prewin(k,:),{1});% makes collapsed rasters for baseline for use in normalizing
+                    [PSR1,N1]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{1});% makes collpased rasters. PSR1 is a cell(neurons)
+                    [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
+                    [PSR3,N0]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},prewin(k,:),{2});% makes trial by trial rasters for baseline for use in response detection
+                    if ~isempty(PSR0{1}) || ~isempty(PSR1{1}) %to avoid errors, added on 12/28 2011
+                        %optimal bin size
+                        %[PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Nrast{j},1),1),{2,1});%-----DP used here
+                        %[PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(i).Nrast{j},1),1),{1,0,BW1});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
+                        %Fixed bin size
+                        [PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Nrast{j},1),1),{2,0,BinSize});%-----Fixed bin used here
+                        [PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(i).Nrast{j},1),1),{1,0,BinSize});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
+                        PTH1=smooth(PTH1,SmoothSPAN,SmoothTYPE)';
+                        PTH0=smooth(PTH0,SmoothSPAN,SmoothTYPE)';
+
+                        %------------- Fills the R_3R.Ev(k) fields --------------
+                        R_3R.Ev(k).BW(NN,1)=BW1;
+                        R_3R.Ev(k).PSTHraw(NN,1:length(Tm))=PTH1;
+                        R_3R.Ev(k).Meanraw(NN,1)=nanmean(R_3R.Ev(k).PSTHraw(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
+                        if sum(PTH0,2)~=0
+                            R_3R.Ev(k).PSTHz(NN,1:length(Tm))=normalize(PTH1,PTH0,0);
+                            R_3R.Ev(k).Meanz(NN,1)=nanmean(R_3R.Ev(k).PSTHz(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
+                        else
+                            R_3R.Ev(k).PSTHz(NN,1:length(Tm))=NaN(1,length(Tm));
+                            R_3R.Ev(k).Meanz(NN,1)=NaN;
+                        end
+
+                        %------------------ firing (in Hz) per trial in pre/post windows ------------------
+                        %used to make the between events comparisons and Response detection in a single window----
+                        ev(k).pre=NaN(size(RAW(i).Erast{EvInd},1),1);
+                        ev(k).post=NaN(size(RAW(i).Erast{EvInd},1),1);
+                        for m=1:size(RAW(i).Erast{EvInd},1) %loops through trials
+                            ev(k).pre(m)=sum(PSR3{m}<prewin(k,2) & PSR3{m}>prewin(k,1))/(prewin(k,2)-prewin(k,1)); %CHANGED FROM PSR2 to PSR0 here 10/24/17
+                            ev(k).post(m)=sum(PSR2{m}<postwin(k,2) & PSR2{m}>postwin(k,1))/(postwin(k,2)-postwin(k,1));
+                        end
+
+                        %-------------------- signrank on event and direction----
+                        [R_3R.Ev(k).signrank(NN,1),~]=signrank(ev(k).pre, ev(k).post); %Signrank used here because it is a dependant sample test
+                        if R_3R.Ev(k).signrank(NN,1)<PStat
+                            R_3R.Ev(k).RespDir(NN,1)=sign(mean(ev(k).post)-mean(ev(k).pre));
+                        else R_3R.Ev(k).RespDir(NN,1)=0;
+                        end
+
+                    end %if ~isempty(PSR0{1}) || ~isempty(PSR1{1})
+                end %if EvInd=0 OR n(trials) < MinNumTrials fills with NaN
+            end %Events
+
+            %----------------------- Check for reward modulation across all the bins surrounding reward delivery-----------------------
+
+            EV1=strmatch('RD1',RAW(i).Einfo(:,2),'exact');
+            EV2=strmatch('RD2',RAW(i).Einfo(:,2),'exact');
+            EV3=strmatch('RD3',RAW(i).Einfo(:,2),'exact');
+
+            Bcell1=0;
+            Bcell2=0;
+            Bcell3=0;
+            BHz=0;
+            Trial=0;
+            Reward=0;
+
+            %get mean baseline firing for all reward trials
+            [Bcell1,B1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV1},BinBase,{2});% makes trial by trial rasters for baseline
+            for y= 1:B1n
+                BHz(1,y)=sum(Bcell1{1,y}>BinBase(1))/(BinBase(1,2)-BinBase(1,1));
+                Trial(1,y)=y;
+                Reward(1,y)=1;
+            end
+            [Bcell2,B2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV2},BinBase,{2});% makes trial by trial rasters for baseline
+            for z= 1:B2n
+                BHz(1,z+B1n)=sum(Bcell2{1,z}>BinBase(1))/(BinBase(1,2)-BinBase(1,1));
+                Trial(1,z+B1n)=z+B1n;
+                Reward(1,z+B1n)=2;
+            end
+            [Bcell3,B3n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV3},BinBase,{2});% makes trial by trial rasters for baseline
+            for x= 1:B3n
+                BHz(1,x+B2n+B1n)=sum(Bcell3{1,x}>BinBase(1))/(BinBase(1,2)-BinBase(1,1));
+                Trial(1,x+B2n+B1n)=x+B2n+B1n;
+                Reward(1,x+B2n+B1n)=3;
+            end
+
+            for k=1:bins %go through each bin for each neuron
+                EVHz=0;
+                EV1cell=0;
+                EV2cell=0;
+                EV3Cell=0;
+
+                %get trial by trial firing rate for water trials
+                [EV1cell,EV1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV1},[(BinDura(1)+(binstart - binint)+k*binint) (BinDura(2)+(binstart - binint)+k*binint)],{2});% makes trial by trial rasters for event
+                for y= 1:EV1n
+                    EVHz(1,y)=sum(EV1cell{1,y}>binstart)/(BinDura(1,2)-BinDura(1,1));
+                end       
+
+                %get trial by trial firing rate for maltodextrin trials
+                [EV2cell,EV2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV2},[(BinDura(1)+(binstart - binint)+k*binint) (BinDura(2)+(binstart - binint)+k*binint)],{2});% makes trial by trial rasters for event
+                for z= 1:EV2n
+                    EVHz(1,z+EV1n)=sum(EV2cell{1,z}>binstart)/(BinDura(1,2)-BinDura(1,1));
+                end
+
+                %get trial by trial firing rate for water trials
+                [EV3cell,EV3n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV3},[(BinDura(1)+(binstart - binint)+k*binint) (BinDura(2)+(binstart - binint)+k*binint)],{2});% makes trial by trial rasters for event
+                for x= 1:EV3n
+                    EVHz(1,x+EV1n+EV2n)=sum(EV3cell{1,x}>binstart)/(BinDura(1,2)-BinDura(1,1));
+                end    
+
+                [~,Bintbl,~]=anovan([BHz EVHz],{[Reward Reward],cat(2,zeros(1,length(Trial)),ones(1,length(Trial)))},'varnames',{'Reward','Pre vs Post'},'model','full','display','off');
+                R_3R.BinStatRwrd{k,1}(NN).IntSig=cell2mat(Bintbl(4,7));
+                
+                %get mean firing rates
+                R_3R.BinStatRwrd{k,1}(NN).R1Mean=(nanmean(EVHz(Reward==1))-nanmean(BHz(Reward==1)))/std(BHz(Reward==1));
+                R_3R.BinStatRwrd{k,1}(NN).R2Mean=(nanmean(EVHz(Reward==2))-nanmean(BHz(Reward==2)))/std(BHz(Reward==2));
+                R_3R.BinStatRwrd{k,1}(NN).R3Mean=(nanmean(EVHz(Reward==3))-nanmean(BHz(Reward==3)))/std(BHz(Reward==3));
+
+                %signranks for checking if inhibitory or excitatory
+
+                %sucrose trials
+                signranksig=signrank(BHz(Reward==1), EVHz(Reward==1)); %Signrank used here because it is a dependant sample test
+                if signranksig<PStat
+                    R_3R.BinStatRwrd{k,1}(NN).SucRespDir=sign(nanmean(EVHz(Reward==1))-mean(BHz(Reward==1)));
+                else
+                    R_3R.BinStatRwrd{k,1}(NN).SucRespDir=0;
+                end
+
+                %maltodextrin trials
+                signranksig=signrank(BHz(Reward==2), EVHz(Reward==2)); %Signrank used here because it is a dependant sample test
+                if signranksig<PStat
+                    R_3R.BinStatRwrd{k,1}(NN).MalRespDir=sign(nanmean(EVHz(Reward==2))-mean(BHz(Reward==2)));
+                else
+                    R_3R.BinStatRwrd{k,1}(NN).MalRespDir=0;
+                end
+
+                %water trials
+                signranksig=signrank(BHz(Reward==3), EVHz(Reward==3)); %Signrank used here because it is a dependant sample test
+                if signranksig<PStat
+                    R_3R.BinStatRwrd{k,1}(NN).WatRespDir=sign(nanmean(EVHz(Reward==3))-mean(BHz(Reward==3)));
+                else
+                    R_3R.BinStatRwrd{k,1}(NN).WatRespDir=0;
+                end
+
+            end %bins for bins analysis
+
+
+            
+            fprintf('3R neuron #%d\n',NN);
+        end %neurons: FOR j= 1:size(RAW(i).Nrast,1)
+    end %if right kind of session
+end %sessions: FOR i=1:length(RAW)
+
+if SAVE_FLAG
+    save('R_3R.mat','R_3R')
+end
+
+toc
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

+ 525 - 0
MatlabScripts/c_Fig1Behavior.m

@@ -0,0 +1,525 @@
+%this program analyzes behavior and puts it in the Lick struct
+%can change what it's saved as here to later compare Int and Blocks
+clear all;
+load ('RAW.mat');
+load ('R_2R.mat');
+
+%which lick intervals to look at?
+startinterval=1;
+endinterval=30;
+
+%time to look at difference in licking activity
+earlylicks=[1 4.5];
+
+%% conduct lick analysis
+global Dura Tm BSIZE Tbin
+tic
+path='C:\Users\dottenh2\Documents\MATLAB\David\2Rewards Nex Files\2RLick_paper.xls';
+if ~exist('RAW'), load ('RAW.mat'); end
+
+%Main settings
+BSIZE=0.01; %Do not change
+Dura=[-22 20]; Tm=Dura(1):BSIZE:Dura(2);
+Tbin=-0.5:0.005:0.5; %window used to determine the optimal binsize
+MinNumTrials=5;
+Lick=[];Lick.Ninfo={};LL=0;NSess=0;
+
+%Smoothing
+Smoothing=1; %0 for raw and 1 for smoothing
+SmoothTYPE='lowess'; %can change this between lowess and rlowess (more robust, ignores outliers more)
+SmoothSPAN=50; %percentage of total data points
+if Smoothing~=1, SmoothTYPE='NoSmoothing';SmoothSPAN=NaN; end
+
+% List of events to analyze and analysis windows EXTRACTED from excel file
+[~,Erefnames]=xlsread(path,'Windows','a3:a10'); % cell that contains the event names
+
+%Finds the total number of suc vs mal sessions
+for i=1:length(RAW)
+    if strcmp('NA',RAW(i).Type(1:2)) | strcmp('VP',RAW(i).Type(1:2))
+        NSess=NSess+1;
+        Lick.Linfo(NSess,1)=RAW(i).Ninfo(1,1);
+    end
+end
+Lick.Erefnames= Erefnames;
+
+
+%preallocating for behavior stats
+Lick.RewDurR1=NaN(NSess,35);
+Lick.RewDurR2=NaN(NSess,35);
+Lick.NumLicksR1=NaN(NSess,35);
+Lick.NumLicksR2=NaN(NSess,35);
+Lick.LicksEarlyR1=NaN(NSess,35);
+Lick.LicksEarlyR2=NaN(NSess,35);
+Lick.LicksR1P1=NaN(NSess,35);
+Lick.LicksR1P2=NaN(NSess,35);
+Lick.LicksR2P1=NaN(NSess,35);
+Lick.LicksR2P2=NaN(NSess,35);
+for i=1:endinterval
+    Lick.ILIR1{i,1}=NaN(NSess,35);
+    Lick.ILIR2{i,1}=NaN(NSess,35);
+end
+
+%preallocating the result matrix
+for k=1:length(Erefnames)
+    Lick.Ev(k).PSTHraw(1:NSess,1:length(Tm))=NaN(NSess,length(Tm));
+    Lick.Ev(k).BW(1:NSess,1)=NaN;
+    Lick.Ev(k).NumberTrials(1:NSess,1)=NaN;
+end
+
+for i=1:length(RAW) %loops through sessions
+    if strcmp('NA',RAW(i).Type(1:2)) | strcmp('VP',RAW(i).Type(1:2))
+        LL=LL+1; %lick session counter
+        for k=1:length(Erefnames) %loops thorough the events
+            EvInd=strcmp(Erefnames(k),RAW(i).Einfo(:,2)); %find the event id number from RAW
+            LickInd=strcmp('Licks',RAW(i).Einfo(:,2)); %find the event id number from RAW
+            if sum(EvInd)==0
+                fprintf('HOWDY, CANT FIND EVENTS FOR ''%s''\n',Erefnames{k});
+            end
+
+            Lick.Ev(k).NumberTrials(LL,1)=length(RAW(i).Erast{EvInd});
+            if  ~isempty(EvInd) && Lick.Ev(k).NumberTrials(LL,1)>MinNumTrials %avoid analyzing sessions where that do not have enough trials
+                [PSR1,N1]=MakePSR04(RAW(i).Erast(LickInd),RAW(i).Erast{EvInd},Dura,{1});% makes collpased rasters. PSR1 is a cell(neurons)
+                if ~isempty(PSR1{1}) %to avoid errors, added on 12/28 2011
+                    %forcing 100ms bin size to keep it consistent across
+                    %sessions (no reason is should be different for licks)
+                    [PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Erast{LickInd},1),1),{2,0,0.1});%these values force bin size to be 100ms
+                    PTH1=smooth(PTH1,SmoothSPAN,SmoothTYPE)';
+
+                    %------------- Fills the R.Ev(k) fields --------------
+                    Lick.Ev(k).BW(LL,1)=BW1;
+                    Lick.Ev(k).PSTHraw(LL,1:length(Tm))=PTH1;
+                end
+            end
+        end
+    
+
+
+
+
+%% this section finds the ILIs, number of licks, and lick duration
+
+
+        RD1=strmatch('RD1',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
+        RD2=strmatch('RD2',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
+        RD1P1=strmatch('RD1P1',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
+        RD1P2=strmatch('RD1P2',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
+        RD2P1=strmatch('RD2P1',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
+        RD2P2=strmatch('RD2P2',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
+        Rind=strmatch('BoutEndRLast',RAW(i).Einfo(:,2),'exact');  %finds the index of the event
+        Lind=strmatch('Licks',RAW(i).Einfo(:,2),'exact');  %finds the index of the event
+
+
+        %trying to make the ISI calculator Peter recommended
+        %unfortunately, there's one trial of maltodextrin with only one
+        %lick
+        %so I need to put this lame code in to avoid an error
+    %         NL2{1,4}=NaN(40,1);
+        %except it doesn't work!
+
+        %find durations of reward consumption for all, R1, and R2 trls
+        RDu1=MakePSR04(RAW(i).Erast(Rind),RAW(i).Erast{RD1},[0 15],{2,'first'});
+        RDu2=MakePSR04(RAW(i).Erast(Rind),RAW(i).Erast{RD2},[0 15],{2,'first'});
+        Lick.RewDurR1(LL,1:length(RDu1))=cell2mat(RDu1);
+        Lick.RewDurR2(LL,1:length(RDu2))=cell2mat(RDu2);
+        %find number of licks for all, R1, and R2 trls
+        NL1=MakePSR04(RAW(i).Erast(Lind),RAW(i).Erast{RD1},[0 15],{2});
+        NL2=MakePSR04(RAW(i).Erast(Lind),RAW(i).Erast{RD2},[0 15],{2});
+        %find number of licks in 1-3 seconds post-reward
+        EL1=MakePSR04(RAW(i).Erast(Lind),RAW(i).Erast{RD1},[earlylicks(1) earlylicks(2)],{2});
+        EL2=MakePSR04(RAW(i).Erast(Lind),RAW(i).Erast{RD2},[earlylicks(1) earlylicks(2)],{2});
+        %break that down by previous reward
+        NL1P1=MakePSR04(RAW(i).Erast(Lind),RAW(i).Erast{RD1P1},[earlylicks(1) earlylicks(2)],{2});
+        NL1P2=MakePSR04(RAW(i).Erast(Lind),RAW(i).Erast{RD1P2},[earlylicks(1) earlylicks(2)],{2});
+        NL2P1=MakePSR04(RAW(i).Erast(Lind),RAW(i).Erast{RD2P1},[earlylicks(1) earlylicks(2)],{2});
+        NL2P2=MakePSR04(RAW(i).Erast(Lind),RAW(i).Erast{RD2P2},[earlylicks(1) earlylicks(2)],{2});
+
+        %make lick matrices
+        for j=1:length(NL1)
+            Lick.NumLicksR1(LL,j)=sum(NL1{1,j}>0);
+            Lick.LicksEarlyR1(LL,j)=sum(EL1{1,j}>0);
+
+            for k=startinterval:endinterval
+                if length(NL1{1,j})>k
+                    Lick.ILIR1{k,1}(LL,j)= NL1{1,j}(k+1,1)-NL1{1,j}(k,1);
+                else
+                    Lick.ILIR1{k,1}(LL,j)= NaN;
+                end
+            end
+        end
+        for j=1:length(NL2)
+            Lick.NumLicksR2(LL,j)=sum(NL2{1,j}>0);
+            Lick.LicksEarlyR2(LL,j)=sum(EL2{1,j}>0);
+
+            for k=startinterval:endinterval
+                if length(NL2{1,j})>k
+                    Lick.ILIR2{k,1}(LL,j)= NL2{1,j}(k+1,1)-NL2{1,j}(k,1);
+                else
+                    Lick.ILIR2{k,1}(LL,j)= NaN;
+                end
+            end
+        end
+
+        %now for the previous reward situation
+        for j=1:length(NL1P1)
+            Lick.LicksR1P1(LL,j)=sum(NL1P1{1,j}>0);
+        end
+        for j=1:length(NL1P2)
+            Lick.LicksR1P2(LL,j)=sum(NL1P2{1,j}>0);
+        end
+        for j=1:length(NL2P1)
+            Lick.LicksR2P1(LL,j)=sum(NL2P1{1,j}>0);
+        end
+        for j=1:length(NL2P2)
+            Lick.LicksR2P2(LL,j)=sum(NL2P2{1,j}>0);
+        end
+
+        %get session and rat #
+        Lick.Sess(LL,1:35)=i;
+        A=char(RAW(i).Ninfo(1,1));
+        Lick.Rat(LL,1:35)=cellstr(A(1:3));
+
+        fprintf('Session # %d\n',LL);
+    end
+end
+
+%% computing stats
+
+%setting up result stat matrix
+R_2R.BehAvs{2,1}='RewDur';
+R_2R.BehAvs{3,1}='NumLicks';
+R_2R.BehAvs{4,1}='NumLicksEarly';
+R_2R.BehAvs{1,2}='R1 Avg';
+R_2R.BehAvs{1,3}='R2 Avg';
+
+%putting all values into one row array and removing NaN, doing 3-way ANOVA
+%with fixed effect of reward and random effects of session nested in rat
+
+%reward duration
+Lick.RDR1=reshape(Lick.RewDurR1,1,numel(Lick.RewDurR1));
+Lick.RDR1Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
+Lick.RDR1Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
+Lick.RDR1Sess(isnan(Lick.RDR1))=[];
+Lick.RDR1Rat(isnan(Lick.RDR1))=[];
+Lick.RDR1(isnan(Lick.RDR1))=[];
+Lick.RDR2=reshape(Lick.RewDurR2,1,numel(Lick.RewDurR2));
+Lick.RDR2Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
+Lick.RDR2Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
+Lick.RDR2Sess(isnan(Lick.RDR2))=[];
+Lick.RDR2Rat(isnan(Lick.RDR2))=[];
+Lick.RDR2(isnan(Lick.RDR2))=[];
+[~,R_2R.RDstats,~]=anovan(cat(2,Lick.RDR1,Lick.RDR2),{(cat(2,zeros(1,length(Lick.RDR1)),ones(1,length(Lick.RDR2)))),(cat(2,Lick.RDR1Sess,Lick.RDR2Sess)),(cat(2,Lick.RDR1Rat,Lick.RDR2Rat))},'varnames',{'reward','session','rat'},'nested',[0 0 0;0 0 1;0 0 0],'random',[2 3],'display','off','model','full'); %nested ANOVA
+
+
+%number of licks
+Lick.NLR1=reshape(Lick.NumLicksR1,1,numel(Lick.NumLicksR1));
+Lick.NLR1Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
+Lick.NLR1Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
+Lick.NLR1Sess(isnan(Lick.NLR1))=[];
+Lick.NLR1Rat(isnan(Lick.NLR1))=[];
+Lick.NLR1(isnan(Lick.NLR1))=[];
+Lick.NLR2=reshape(Lick.NumLicksR2,1,numel(Lick.NumLicksR2));
+Lick.NLR2Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
+Lick.NLR2Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
+Lick.NLR2Sess(isnan(Lick.NLR2))=[];
+Lick.NLR2Rat(isnan(Lick.NLR2))=[];
+Lick.NLR2(isnan(Lick.NLR2))=[];
+[~,R_2R.LicksStats,~]=anovan(cat(2,Lick.NLR1,Lick.NLR2),{(cat(2,zeros(1,length(Lick.NLR1)),ones(1,length(Lick.NLR2)))),(cat(2,Lick.NLR1Sess,Lick.NLR2Sess)),(cat(2,Lick.NLR1Rat,Lick.NLR2Rat))},'varnames',{'reward','session','rat'},'nested',[0 0 0;0 0 1;0 0 0],'random',[2 3],'display','off','model','full'); %nested ANOVA
+
+%number of licks 1-3sec
+Lick.ELR1=reshape(Lick.LicksEarlyR1,1,numel(Lick.LicksEarlyR1));
+Lick.ELR1Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
+Lick.ELR1Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
+Lick.ELR1Sess(isnan(Lick.ELR1))=[];
+Lick.ELR1Rat(isnan(Lick.ELR1))=[];
+Lick.ELR1(isnan(Lick.ELR1))=[];
+Lick.ELR2=reshape(Lick.LicksEarlyR2,1,numel(Lick.LicksEarlyR2));
+Lick.ELR2Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
+Lick.ELR2Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
+Lick.ELR2Sess(isnan(Lick.ELR2))=[];
+Lick.ELR2Rat(isnan(Lick.ELR2))=[];
+Lick.ELR2(isnan(Lick.ELR2))=[];
+[~,R_2R.LicksEarlyStats,~]=anovan(cat(2,Lick.ELR1,Lick.ELR2),{(cat(2,zeros(1,length(Lick.ELR1)),ones(1,length(Lick.ELR2)))),(cat(2,Lick.ELR1Sess,Lick.ELR2Sess)),(cat(2,Lick.ELR1Rat,Lick.ELR2Rat))},'varnames',{'reward','session','rat'},'nested',[0 0 0;0 0 1;0 0 0],'random',[2 3],'display','off','model','full'); %nested ANOVA
+
+%number of licks 1-3sec
+Lick.NLR1P1=reshape(Lick.LicksR1P1,1,numel(Lick.LicksR1P1));
+Lick.NLR1P1Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
+Lick.NLR1P1Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
+Lick.NLR1P1Sess(isnan(Lick.NLR1P1))=[];
+Lick.NLR1P1Rat(isnan(Lick.NLR1P1))=[];
+Lick.NLR1P1(isnan(Lick.NLR1P1))=[];
+Lick.NLR1P2=reshape(Lick.LicksR1P2,1,numel(Lick.LicksR1P2));
+Lick.NLR1P2Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
+Lick.NLR1P2Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
+Lick.NLR1P2Sess(isnan(Lick.NLR1P2))=[];
+Lick.NLR1P2Rat(isnan(Lick.NLR1P2))=[];
+Lick.NLR1P2(isnan(Lick.NLR1P2))=[];
+Lick.NLR2P1=reshape(Lick.LicksR2P1,1,numel(Lick.LicksR2P1));
+Lick.NLR2P1Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
+Lick.NLR2P1Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
+Lick.NLR2P1Sess(isnan(Lick.NLR2P1))=[];
+Lick.NLR2P1Rat(isnan(Lick.NLR2P1))=[];
+Lick.NLR2P1(isnan(Lick.NLR2P1))=[];
+Lick.NLR2P2=reshape(Lick.LicksR2P2,1,numel(Lick.LicksR2P2));
+Lick.NLR2P2Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
+Lick.NLR2P2Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
+Lick.NLR2P2Sess(isnan(Lick.NLR2P2))=[];
+Lick.NLR2P2Rat(isnan(Lick.NLR2P2))=[];
+Lick.NLR2P2(isnan(Lick.NLR2P2))=[];
+[~,R_2R.LicksPrevRewStats,stats]=anovan(cat(2,Lick.NLR1P1,Lick.NLR1P2,Lick.NLR2P1,Lick.NLR2P2),{cat(2,zeros(1,length(cat(2,Lick.NLR1P1,Lick.NLR1P2))),ones(1,length(cat(2,Lick.NLR2P1,Lick.NLR2P2)))),...
+    cat(2,zeros(1,length(Lick.NLR1P1)),ones(1,length(Lick.NLR1P2)),zeros(1,length(Lick.NLR2P1)),ones(1,length(Lick.NLR2P2))),cat(2,Lick.NLR1P1Sess,Lick.NLR1P2Sess,Lick.NLR2P1Sess,Lick.NLR2P2Sess),...
+    cat(2,Lick.NLR1P1Rat,Lick.NLR1P2Rat,Lick.NLR2P1Rat,Lick.NLR2P2Rat)},'varnames',{'current reward','prev reward','session','rat'},'nested',[0 0 0 0;0 0 0 0;0 0 0 1;0 0 0 0],'random',[3 4],'display','off','model','full'); %nested ANOVA
+
+% %without random factor of session
+% [~,R_2R.LicksPrevRewStats,stats]=anovan(cat(2,Lick.NLR1P1,Lick.NLR1P2,Lick.NLR2P1,Lick.NLR2P2),{cat(2,zeros(1,length(cat(2,Lick.NLR1P1,Lick.NLR1P2))),ones(1,length(cat(2,Lick.NLR2P1,Lick.NLR2P2)))),...
+%     cat(2,zeros(1,length(Lick.NLR1P1)),ones(1,length(Lick.NLR1P2)),zeros(1,length(Lick.NLR2P1)),ones(1,length(Lick.NLR2P2))),...
+%     cat(2,Lick.NLR1P1Rat,Lick.NLR1P2Rat,Lick.NLR2P1Rat,Lick.NLR2P2Rat)},'varnames',{'current reward','prev reward','rat'},'random',[3],'display','off','model','full');
+
+
+%ILI
+for m=startinterval:endinterval
+    %first get session averages so I can do an ANOVA across all ILI
+    Lick.ILIR1sessav{m,1}=nanmean(Lick.ILIR1{m,1},2);
+    Lick.ILIR2sessav{m,1}=nanmean(Lick.ILIR2{m,1},2);
+    Lick.ILIsesssess{m,1}=Lick.Sess(:,1);
+    Lick.ILIsessrat{m,1}=Lick.Rat(:,1);
+    
+    
+    %now do it with all trials
+    Lick.ILIR1{m,1}=reshape(Lick.ILIR1{m,1},1,numel(Lick.ILIR1{m,1}));
+    Lick.ILIR1Sess{m,1}=reshape(Lick.Sess,1,numel(Lick.ILIR1{m,1}));
+    Lick.ILIR1Rat{m,1}=reshape(Lick.Rat,1,numel(Lick.ILIR1{m,1}));
+    Lick.ILIR1Sess{m,1}(isnan(Lick.ILIR1{m,1}))=[];
+    Lick.ILIR1Rat{m,1}(isnan(Lick.ILIR1{m,1}))=[];
+    Lick.ILIR1{m,1}(isnan(Lick.ILIR1{m,1}))=[];
+    Lick.ILIR2{m,1}=reshape(Lick.ILIR2{m,1},1,numel(Lick.ILIR2{m,1}));
+    Lick.ILIR2Sess{m,1}=reshape(Lick.Sess,1,numel(Lick.ILIR2{m,1}));
+    Lick.ILIR2Rat{m,1}=reshape(Lick.Rat,1,numel(Lick.ILIR2{m,1}));
+    Lick.ILIR2Sess{m,1}(isnan(Lick.ILIR2{m,1}))=[];
+    Lick.ILIR2Rat{m,1}(isnan(Lick.ILIR2{m,1}))=[];
+    Lick.ILIR2{m,1}(isnan(Lick.ILIR2{m,1}))=[];
+end
+
+%getting effects of interval #, reward, rat, and session on ILI
+Lick.ILIR1all=[];
+Lick.ILIR1all=[];
+Lick.ILIR1Ratall=[];
+Lick.ILIR1Sessall=[];
+Lick.ILIR1Rewall=[];
+Lick.ILIR1Intall=[];
+Lick.ILIR2all=[];
+Lick.ILIR2all=[];
+Lick.ILIR2Ratall=[];
+Lick.ILIR2Sessall=[];
+Lick.ILIR2Rewall=[];
+Lick.ILIR2Intall=[];
+
+%to do session averages for comparing across all intervals
+for m=startinterval:endinterval
+    Lick.ILIR1all=cat(1,Lick.ILIR1all,Lick.ILIR1sessav{m,1});
+    Lick.ILIR1Ratall=cat(1,Lick.ILIR1Ratall,Lick.ILIsessrat{m,1});
+    Lick.ILIR1Sessall=cat(1,Lick.ILIR1Sessall,Lick.ILIsesssess{m,1});
+    Lick.ILIR1Rewall=cat(1,Lick.ILIR1Rewall,zeros(length(Lick.ILIR1sessav{m,1}),1));
+    Lick.ILIR1Intall=cat(1,Lick.ILIR1Intall,(m*ones(length(Lick.ILIR1sessav{m,1}),1)));
+    Lick.ILIR2all=cat(1,Lick.ILIR2all,Lick.ILIR2sessav{m,1});
+    Lick.ILIR2Ratall=cat(1,Lick.ILIR2Ratall,Lick.ILIsessrat{m,1});
+    Lick.ILIR2Sessall=cat(1,Lick.ILIR2Sessall,Lick.ILIsesssess{m,1});
+    Lick.ILIR2Rewall=cat(1,Lick.ILIR2Rewall,ones(length(Lick.ILIR2sessav{m,1}),1));
+    Lick.ILIR2Intall=cat(1,Lick.ILIR2Intall,(m*ones(length(Lick.ILIR2sessav{m,1}),1)));    
+end
+
+%no factor of session
+[~,R_2R.ILItrialstats,~]=anovan(cat(1,Lick.ILIR1all,Lick.ILIR2all),{(cat(1,Lick.ILIR1Rewall,Lick.ILIR2Rewall)),(cat(1,Lick.ILIR1Ratall,Lick.ILIR2Ratall)),(cat(1,Lick.ILIR1Intall,Lick.ILIR2Intall))},'varnames',{'reward','rat','interval'},'random',[2],'display','off','model','full');
+
+%finding values for R_2R.BehAvs
+R_2R.BehAvs{2,2}=mean(Lick.RDR1);
+R_2R.BehAvs{2,3}=mean(Lick.RDR2);
+R_2R.BehAvs{3,2}=mean(Lick.NLR1);
+R_2R.BehAvs{3,3}=mean(Lick.NLR2);
+R_2R.BehAvs{4,2}=mean(Lick.ELR1);
+R_2R.BehAvs{4,3}=mean(Lick.ELR2);
+
+%Finding ILI
+for m=startinterval:endinterval
+    Lick.ILI1(m,1)=nanmean(Lick.ILIR1{m,1}(:));
+    Lick.ILI1(m,2)=nanstd(Lick.ILIR1{m,1}(:))/sqrt(sum(~isnan(Lick.ILIR1{m,1}(:))));
+    Lick.ILI2(m,1)=nanmean(Lick.ILIR2{m,1}(:));
+    Lick.ILI2(m,2)=nanstd(Lick.ILIR2{m,1}(:))/sqrt(sum(~isnan(Lick.ILIR2{m,1}(:))));
+end
+
+%% plotting licks in each reward condition
+Xaxis=[-2 12];
+Ishow=find(Tm>=Xaxis(1) & Tm<=Xaxis(2));
+time1=Tm(Ishow);
+c=[-1000 7500];ClimE=sign(c).*abs(c).^(1/4);%ColorMapExc
+colormap(jet);
+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];
+
+Ev1=strcmp('RD1', Lick.Erefnames);
+Ev2=strcmp('RD2', Lick.Erefnames);
+
+psth1=nanmean(Lick.Ev(Ev1).PSTHraw(:,Ishow),1);
+sem1=nanste(Lick.Ev(Ev1).PSTHraw(:,Ishow),1); %calculate standard error of the mean
+up1=psth1+sem1;
+down1=psth1-sem1;
+
+psthE=nanmean(Lick.Ev(Ev2).PSTHraw(:,Ishow),1);
+semE=nanste(Lick.Ev(Ev2).PSTHraw(:,Ishow),1); %calculate standard error of the mean
+upE=psthE+semE;
+downE=psthE-semE;
+
+%plotting
+subplot(2,2,1);
+hold on;
+plot(time1,psth1,'Color',sucrose,'linewidth',1);
+plot(time1,psthE,'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)],[upE,downE(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
+title('Mean lick rate');
+%plot([-2 10],[0 0],':','color','k');
+plot([0 0],[-2 8],':','color','k','linewidth',0.75);
+plot([earlylicks(1) earlylicks(2)],[7.3 7.3],'color','k','linewidth',0.75);
+axis([-2 12 0 7.5]);
+xlabel('Seconds from reward delivery');
+ylabel('Licks/s');
+legend('Sucrose trials','Maltodextrin trials');
+
+
+%% plotting ILI
+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];
+xaxis=(startinterval:endinterval);
+
+subplot(2,2,3);
+hold on;
+errorbar(xaxis,Lick.ILI1(startinterval:endinterval,1),Lick.ILI1(startinterval:endinterval,2),'*','Color',sucrose,'linewidth',1.5);
+errorbar(xaxis,Lick.ILI2(startinterval:endinterval,1),Lick.ILI2(startinterval:endinterval,2),'*','Color',maltodextrin,'linewidth',1.5);
+xlabel('Lick Interval #');
+ylabel('Duration (s)');
+title(['Interlick intervals'])
+axis([1 30 0.13 0.19]);
+legend('Sucrose trials','Maltodextrin trials','location','northwest');
+
+subplot(4,4,15);
+hold on;
+errorbar(2,nanmean(Lick.ILI1(startinterval:endinterval,1)),nanste(Lick.ILI1(startinterval:endinterval,1),1),'*','Color',sucrose,'linewidth',1.5);
+errorbar(2,nanmean(Lick.ILI2(startinterval:endinterval,1)),nanste(Lick.ILI2(startinterval:endinterval,1),1),'*','Color',maltodextrin,'linewidth',1.5);
+if R_2R.LicksEarlyStats{2,7}<0.05
+    plot(2,nanmean(Lick.ILI2(startinterval:endinterval,1))+0.005,'*','Color','k','MarkerSize',5);
+end
+set(gca,'xtick',[]);
+axis([1.5 2.5 0.14 0.17]);
+
+
+%% plotting sucrose preferences
+
+%Plot sucrose preference
+
+%colors
+NAc=[0.5  0.1  0.8];
+VP=[0.3  0.7  0.1];
+
+%data
+x=[1 2];
+NAcPref{1}=[81.7 81.3];
+NAcPref{2}=[77 73.9];
+NAcPref{3}=[37.4 74.6];
+NAcPref{4}=[71 95.6];
+NAcPref{5}=[89 83.6];
+NAcPref{6}=[94 93];
+VPPref{1}=[59.6 78.8];
+VPPref{2}=[88.8 85.8];
+VPPref{3}=[71.2 89];
+VPPref{4}=[87 76.6];
+VPPref{5}=[85.2 84.1];
+
+subplot(2,4,3);
+hold on;
+%do the first 2 first in order to get the legend right
+plot(x,NAcPref{1},'Marker','o','MarkerFaceColor',NAc,'LineWidth',1,'Color',NAc);
+plot(x,VPPref{1},'Marker','o','MarkerFaceColor',VP,'LineWidth',1,'Color',VP);
+
+for i=2:length(NAcPref)
+    plot(x,NAcPref{i},'Marker','o','MarkerFaceColor',NAc,'LineWidth',1,'Color',NAc);
+end
+
+for i=2:length(VPPref)
+    plot(x,VPPref{i},'Marker','o','MarkerFaceColor',VP,'LineWidth',1,'Color',VP);
+end
+
+axis([0.7 2.2 0 100]);
+title('Two-bottle choice');
+xlabel('Initial       Final');
+ylabel('% sucrose consumption');
+legend('NAc rat','VP rat','location','southeast');
+plot([0 3],[50 50],':','color','k','linewidth',0.75);
+set(gca,'xtick',[])
+
+
+%% licks for each rat
+figure;
+
+Xaxis=[-2 12];
+Ishow=find(Tm>=Xaxis(1) & Tm<=Xaxis(2));
+time1=Tm(Ishow);
+
+%colors
+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];
+inh=[0.1 0.021154 0.6];
+exc=[0.9 0.75 0.205816];
+NAc=[0.5  0.1  0.8];
+VP=[0.3  0.7  0.1];
+
+%split data up into individual rats
+for i=1:length(Lick.Linfo)
+    A=char(Lick.Linfo(i,1));
+    Rats(i,:)=cellstr(A(1:3));
+end
+
+C=unique(Rats);
+
+for i=1:length(C)
+    Sel=strcmp(C(i),Rats);
+    Ev1=strcmp('RD1', Lick.Erefnames);
+    Ev2=strcmp('RD2', Lick.Erefnames);
+
+    psth1=nanmean(Lick.Ev(Ev1).PSTHraw(Sel,Ishow),1);
+    sem1=nanste(Lick.Ev(Ev1).PSTHraw(Sel,Ishow),1); %calculate standard error of the mean
+    up1=psth1+sem1;
+    down1=psth1-sem1;
+
+    psthE=nanmean(Lick.Ev(Ev2).PSTHraw(Sel,Ishow),1);
+    semE=nanste(Lick.Ev(Ev2).PSTHraw(Sel,Ishow),1); %calculate standard error of the mean
+    upE=psthE+semE;
+    downE=psthE-semE;
+
+    %plotting
+    subplot(4,3,i);
+    hold on;
+    plot(time1,psth1,'Color',sucrose,'linewidth',1);
+    plot(time1,psthE,'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)],[upE,downE(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
+    title([cell2mat(C(i)) ' (n=' num2str(sum(Sel)) ')' ]);
+    %plot([-2 10],[0 0],':','color','k');
+    plot([0 0],[-2 8],':','color','k','linewidth',0.75);
+    patch([earlylicks(1) earlylicks(1) earlylicks(2) earlylicks(2)],[5 8 8 5],'k','facecolor','none','edgecolor',[0.7 0.7 0.7]);alpha(0.5);
+    axis([-2 12 0 8]);
+    xlabel('Seconds from reward delivery');
+    ylabel('Licks/s');
+    if i==1 legend('Sucrose trials','Maltodextrin trials'); end
+
+end
+
+save('R_2R.mat','R_2R');

+ 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(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

+ 97 - 0
MatlabScripts/e_SingleUnitDecoding.m

@@ -0,0 +1,97 @@
+%decodes trial identity from spiking across bins for each individual neuron
+clear all;
+load ('R_2R.mat');
+load ('RAW.mat');
+
+folds = 5; %number of times cross-validated:
+shuffs = 1; %number of shuffled models created
+
+%load parameters
+BinDura=R_2R.Param.BinDura;
+bins=R_2R.Param.bins;
+binint=R_2R.Param.binint;
+binstart=R_2R.Param.binstart;
+
+
+
+%setup variables
+UnitDec=[];
+NN = 0;
+
+for i=1:length(RAW) %loops through sessions
+    if strcmp('NA',RAW(i).Type(1:2)) | strcmp('VP',RAW(i).Type(1:2)) %check if it's a suc vs mal session
+        %events being compared
+        RD1=strcmp('RD1', RAW(i).Einfo(:,2));
+        RD2=strcmp('RD2', RAW(i).Einfo(:,2));
+        for j= 1:size(RAW(i).Nrast,1) %Number of neurons per session
+            NN=NN+1; %neuron counter
+            for l=1:bins
+
+                LowHz=zeros(1,1);
+                DecodeSpikes=NaN(1,1);
+                DecodeRs=NaN(1,1);
+
+                [PSR1,N1]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{RD1},[(BinDura(1)+(binstart - binint)+l*binint) (BinDura(2)+(binstart - binint)+l*binint)],{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
+                for m=1:length(PSR1)
+                    DecodeSpikes(m,1)=sum(PSR1{1,m}>(binstart));
+                    DecodeRs(m,1)=1;
+                end
+
+                %get all the spikes from reward 1p2 trials
+                [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{RD2},[(BinDura(1)+(binstart - binint)+l*binint) (BinDura(2)+(binstart - binint)+l*binint)],{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
+                for n=1:length(PSR2)
+                    DecodeSpikes(n+m,1)=sum(PSR2{1,n}>(binstart));
+                    DecodeRs(n+m,1)=2;
+                end
+
+                if sum(DecodeSpikes(1:N1,1))<7 || sum(DecodeSpikes((1+N1):(N1+N2),1))<7 %prevents errors with LDA on conditions with no variance
+                    LowHz(1,1)=1;
+                end
+
+
+                %if one neuron or more has too few spikes, so LDA has error because not
+                %enough variance, this removes those neurons from the analysis
+                if sum(LowHz)==0
+
+                    %creating models
+                    CVacc = NaN(folds,1);
+                    CVaccSh = NaN(folds,1);
+                    %normal model
+                    Partitions = cvpartition(DecodeRs,'KFold',folds);
+                    for r = 1:folds
+
+                        LDAModel = fitcdiscr(DecodeSpikes(Partitions.training(r),:),DecodeRs(Partitions.training(r)));
+                        prediction = predict(LDAModel,DecodeSpikes(Partitions.test(r),:));
+                        actual = DecodeRs(Partitions.test(r));
+                        correct = prediction - actual;
+                        CVacc(r) = sum(correct==0) / length(correct);
+
+                    end
+
+                    %shuffled model
+                    for q=1:shuffs
+                        DecodeRsSh=DecodeRs(randperm(length(DecodeRs)));
+                        PartitionsSh = cvpartition(DecodeRsSh,'KFold',folds);
+                        for s = 1:folds
+                            LDAModelSh = fitcdiscr(DecodeSpikes(PartitionsSh.training(s),:),DecodeRsSh(PartitionsSh.training(s)));
+                            predictionSh = predict(LDAModelSh,DecodeSpikes(PartitionsSh.test(s),:));
+                            actualSh = DecodeRsSh(PartitionsSh.test(s));
+                            correctSh = predictionSh - actualSh;
+                            CVaccSh(s) = sum(correctSh==0) / length(correctSh);
+                        end
+                        AccShuff(q,1) = nanmean(CVaccSh);
+                    end
+                    UnitDec.True(NN,l) = nanmean(CVacc);
+                    UnitDec.Shuff(NN,l) = nanmean(AccShuff);
+                else 
+                    UnitDec.True(NN,l) = NaN;
+                    UnitDec.Shuff(NN,l) = NaN;
+                end
+            end
+
+            fprintf('Neuron ID # %d\n',NN);
+        end %neurons: FOR j= 1:size(RAW(i).Nrast,1)
+    end %checking if the right session type
+end %sessions: FOR i=1:length(RAW)
+
+save('UnitDec_2R.mat','UnitDec');

+ 201 - 0
MatlabScripts/f_PooledEnsDecoding.m

@@ -0,0 +1,201 @@
+%decodes trial identity from spiking across bins for neurons pooled across sessions
+%also does this for only selective and non-selective neurons
+
+%need to run A, B, and D first
+clear all;
+load ('R_2R.mat');
+load ('RAW.mat');
+
+%NAc is 1, VP is 2
+region={'NA';'VP'};
+
+%decoding parameters
+NumNeurons=[10 25 50 100 150]; %matrix of how many neurons to use on each iteration
+repetitions=50; %how many times to run the analysis
+folds = 5; %number of times cross-validated:
+shuffs = 1; %number of shuffled models created
+
+%load parameters
+BinDura=R_2R.Param.BinDura;
+bins=R_2R.Param.bins;
+binint=R_2R.Param.binint;
+binstart=R_2R.Param.binstart;
+
+%Need to keep number of trials used in VP and NAc constant, so have to pick
+%minimum number across all sessions in both regions. Right now with
+%existing sessions it's 20 Ev1s and 20 Ev2s
+Ev1s=20;
+Ev2s=20;
+
+
+
+%find number of trials in each session
+for i=1:length(RAW)
+    
+    %events being compared
+    Ev1=strcmp('RD1', RAW(i).Einfo(:,2));
+    Ev2=strcmp('RD2', RAW(i).Einfo(:,2));
+    
+    Ev1perSess(i)=length(RAW(i).Erast{Ev1,1});
+    Ev2perSess(i)=length(RAW(i).Erast{Ev2,1});
+end
+
+%setup variables
+PoolDec=[];
+NN = 0;
+
+%if you change this to e=1:3, you can also look at nonselective neurons
+for e=1:2 %different selections of neurons
+
+    
+    %pick which set of neurons: all, reward-specific, or non-reward specific
+    if e==1 selection=R_2R.Ninfo; end %all neurons
+    if e==2 selection=R_2R.RSinfo; end %reward-selective neurons
+    if e==3 selection=R_2R.RNSinfo; end %reward-nonselective neurons
+    
+    for f=1:2 %region
+
+        TotalNeurons = 0;
+
+        %total number of neurons in all sessions and events per session
+        for i=1:length(RAW)
+            if strcmp(region(f),RAW(i).Type(1:2))
+                TotalNeurons=TotalNeurons+size(RAW(i).Nrast,1);
+            end
+        end
+
+        for l=1:bins
+            DecodeSpikes=NaN(1,1);
+            AllSpikes=NaN(1,1);
+            NN = 0;
+            for i=1:length(RAW) %loops through sessions
+                if strcmp(region(f),RAW(i).Type(1:2)) && Ev1perSess(i)>=20 && Ev2perSess(i)>=20
+                    
+                    %events being compared
+                    Ev1=strcmp('RD1', RAW(i).Einfo(:,2));
+                    Ev2=strcmp('RD2', RAW(i).Einfo(:,2));
+                    
+                    for j= 1:size(RAW(i).Nrast,1) %Number of neurons per session
+                        if sum(strcmp(RAW(i).Ninfo(j,1),selection(:,1)) & strcmp(RAW(i).Ninfo(j,2),selection(:,2)))==1 %check if this neuron is on the selection list
+                            NN=NN+1; %neuron counter
+
+                            Ev1Spikes=NaN(1,1);
+                            Ev2Spikes=NaN(1,1);
+
+                            %get number of spikes for this neuron in this bin for all
+                            %Ev1 trials
+                            [PSR1,N1]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{Ev1},[(BinDura(1)+(binstart - binint)+l*binint) (BinDura(2)+(binstart - binint)+l*binint)],{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
+                            for m=1:length(PSR1)
+                                Ev1Spikes(m,1)=sum(PSR1{1,m}>(binstart));
+                            end
+
+                            %pick which trials get used for decoding
+                            SetupTrials=cat(1,ones(Ev1s,1),zeros(length(PSR1)-Ev1s,1));
+                            Trials=(SetupTrials(randperm(length(SetupTrials)))==1);
+
+                            %put spikes from those trials in a matrix
+                            AllSpikes(1:Ev1s,NN)=Ev1Spikes(Trials,1);
+
+                            %get all the spikes from reward 1p2 trials
+                            [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{Ev2},[(BinDura(1)+(binstart - binint)+l*binint) (BinDura(2)+(binstart - binint)+l*binint)],{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
+                            for n=1:length(PSR2)
+                                Ev2Spikes(n,1)=sum(PSR2{1,n}>(binstart));
+                            end
+
+                            %pick which trials get used for decoding
+                            SetupTrials2=cat(1,ones(Ev2s,1),zeros(length(PSR2)-Ev2s,1));
+                            Trials2=(SetupTrials2(randperm(length(SetupTrials2)))==1);
+
+                            %put spikes from those trials in a matrix
+                            AllSpikes((Ev1s+1):(Ev1s+Ev2s),NN)=Ev2Spikes(Trials2,1);
+                        end
+                    end %neurons
+                end %checking if enough events in that session
+            end %sessions
+
+            TotalNeurons=NN;    
+            
+            %do the decoding
+            for v = 1:length(NumNeurons)
+                if TotalNeurons>NumNeurons(v)
+                    for u=1:repetitions
+                        LowHz=zeros(1,NumNeurons(v)); %reset the lowHz identifier
+
+                        %pick which neurons to use
+                        SetupSel=cat(1,ones(NumNeurons(v),1),zeros(TotalNeurons-NumNeurons(v),1));
+                        Sel=(SetupSel(randperm(length(SetupSel)))==1);
+
+                        %setup the event identity matrix for decoding
+                        %DecodeRs=zeros(Ev1s+Ev2s,NumNeurons(v));
+                        DecodeRs(1:Ev1s,1)=1;
+                        %DecodeRs(1:Ev1s,end)=1;
+                        %DecodeRs((Ev1s+1):(Ev1s+Ev2s),2:end-1)=1;
+                        DecodeRs((Ev1s+1):(Ev1s+Ev2s),1)=0;
+                    
+                        %setup decode spike matrix
+                        DecodeSpikes=AllSpikes(:,Sel);
+                        
+                        %find neurons with too few spikes
+
+                        for t=1:NumNeurons(v)
+                            if sum(DecodeSpikes(1:Ev1s,t))<7 || sum(DecodeSpikes((1+Ev1s):(Ev1s+Ev2s),t))<7
+                                LowHz(1,t)=1;
+                            end
+                        end
+
+                        %remove neurons with too few spikes
+                        if sum(LowHz)>0
+                            Sel2=LowHz<1;
+                            DecodeSpikes=DecodeSpikes(:,Sel2);
+                        end
+
+                        %set up validation result matrices
+                        CVacc = NaN(folds,1);
+                        CVaccSh = NaN(folds,1);
+
+                        %normal model
+                        Partitions = cvpartition(DecodeRs,'KFold',folds);
+                        for r = 1:folds
+
+                            SVMModel = fitcdiscr(DecodeSpikes(Partitions.training(r),:),DecodeRs(Partitions.training(r)));
+                            prediction = predict(SVMModel,DecodeSpikes(Partitions.test(r),:));
+                            actual = DecodeRs(Partitions.test(r));
+                            correct = prediction - actual;
+                            CVacc(r) = sum(correct==0) / length(correct);
+
+                        end
+
+                        %shuffled model
+                        for q=1:shuffs
+                            DecodeRsSh=DecodeRs(randperm(length(DecodeRs)));
+                            PartitionsSh = cvpartition(DecodeRsSh,'KFold',folds);
+                            for s = 1:folds
+                                SVMModelSh = fitcdiscr(DecodeSpikes(PartitionsSh.training(s),:),DecodeRsSh(PartitionsSh.training(s)));
+                                predictionSh = predict(SVMModelSh,DecodeSpikes(PartitionsSh.test(s),:));
+                                actualSh = DecodeRsSh(PartitionsSh.test(s));
+                                correctSh = predictionSh - actualSh;
+                                CVaccSh(s) = sum(correctSh==0) / length(correctSh);
+                            end
+                            AccShuff(q,1) = nanmean(CVaccSh);
+                        end
+
+                        PoolDec{e,f}.True{v,1}(u,l) = nanmean(CVacc);
+                        PoolDec{e,f}.Shuff{v,1}(u,l) = nanmean(AccShuff);
+
+                        
+                    end %repetitions
+                    fprintf(['Selection #' num2str(e) ', Region #' num2str(f) ', Bin #' num2str(l) ', Condition #' num2str(v) '\n']);
+                else
+                    PoolDec{e,f}.True{v,1} = NaN;
+                    PoolDec{e,f}.Shuff{v,1} = NaN;
+                    fprintf(['Selection #' num2str(e) ', Region #' num2str(f) ', Bin #' num2str(l) ', Condition #' num2str(v) '\n']);
+                end %checking if there are enough neurons
+            end %conditions (number of neurons)
+        end %bins
+    end %region
+end %selections
+
+save('PoolDec_2R.mat','PoolDec');
+
+R_2R.Param.NumNeurons=NumNeurons;
+save('R_2R.mat','R_2R');

+ 421 - 0
MatlabScripts/g_Fig3FigS6Decoding.m

@@ -0,0 +1,421 @@
+%analyze and plot the decoding data
+%NOTE: this program takes a really long time because of large ANOVAs
+
+clear all;
+load ('R_2R.mat');
+load ('RAW.mat');
+load ('PoolDec_2R.mat');
+load ('UnitDec_2R.mat');
+
+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];
+NAc=[0.5  0.1  0.8];
+VP=[0.3  0.7  0.1];
+NAcP{5,1}=[1  0.7  1];
+NAcP{4,1}=[0.85  0.3  0.9];
+NAcP{3,1}=[0.6  0.1  0.8];
+NAcP{2,1}=[0.3  0.05  0.5];
+NAcP{1,1}=[0.2  0  0.3];
+VPP{5,1}=[0.6  1  0.25];
+VPP{4,1}=[0.4  0.9  0.15];
+VPP{3,1}=[0.3  0.7  0.1];
+VPP{2,1}=[0.1  0.4  0.05];
+VPP{1,1}=[0.03  0.2  0];
+NAcShuff=[0.3  0.05  0.5];
+VPShuff=[0.05  0.4  0];
+
+%load parameters
+BinDura=R_2R.Param.BinDura;
+bins=R_2R.Param.bins;
+binint=R_2R.Param.binint;
+binstart=R_2R.Param.binstart;
+NumNeurons=R_2R.Param.NumNeurons;
+repetitions=length(PoolDec{1,1}.True{1,1}(:,1));
+testbins=R_2R.Param.SelectiveBins;
+
+xaxis=linspace(binstart+BinDura(2)/2,binstart+(bins-1)*binint+BinDura(2)/2,bins); %include all bins
+
+%% single unit decoding
+
+%divide neurons up by region
+NAneurons=strcmp(R_2R.Ninfo(:,3),'NA');
+VPneurons=strcmp(R_2R.Ninfo(:,3),'VP');
+
+%if you change this to e=1:3, you can also look at nonselective neurons
+for e=1:2 %different selections of neurons
+    
+    figure;
+    
+    %pick which set of neurons: all, reward-specific, or non-reward specific
+    if e==1 selection=R_2R.SucN<2; end %all neurons
+    if e==2 selection=R_2R.SucN | R_2R.MalN; end %reward-selective neurons
+    if e==3 selection=(R_2R.SucN | R_2R.MalN) == 0; end %reward-nonselective neurons
+
+    %get average accuracy for each bin and run stats comparing region
+    for i = 1:bins
+        AvgAccNAc(1,i)=nanmean(UnitDec.True(NAneurons&selection,i)); %average accuracy NAc
+        SEMAccNAc(1,i)=nanste(UnitDec.True(NAneurons&selection,i),1); %SEM
+        AvgAccNAcShuff(1,i)=nanmean(UnitDec.Shuff(NAneurons&selection,i)); %average accuracy NAcShuff
+        SEMAccNAcShuff(1,i)=nanste(UnitDec.Shuff(NAneurons&selection,i),1); %SEM
+        AvgAccVP(1,i)=nanmean(UnitDec.True(VPneurons&selection,i)); %average accuracy NAcShuff
+        SEMAccVP(1,i)=nanste(UnitDec.True(VPneurons&selection,i),1); %SEM
+        AvgAccVPShuff(1,i)=nanmean(UnitDec.Shuff(VPneurons&selection,i)); %average accuracy NAcShuff
+        SEMAccVPShuff(1,i)=nanste(UnitDec.Shuff(VPneurons&selection,i),1); %SEM
+      
+    end
+
+    %prepare shading
+    upSEMNAc=AvgAccNAc+SEMAccNAc;
+    downSEMNAc=AvgAccNAc-SEMAccNAc;
+    upSEMVP=AvgAccVP+SEMAccVP;
+    downSEMVP=AvgAccVP-SEMAccVP;
+    upSEMNAcShuff=AvgAccNAcShuff+SEMAccNAcShuff;
+    downSEMNAcShuff=AvgAccNAcShuff-SEMAccNAcShuff;
+    upSEMVPShuff=AvgAccVPShuff+SEMAccVPShuff;
+    downSEMVPShuff=AvgAccVPShuff-SEMAccVPShuff;
+
+    %plotting decoder accuracies over time
+    subplot(2,3,1);
+    hold on;
+    plot(xaxis,AvgAccNAc(1:bins),'Color', NAc,'linewidth',1.5); %accumbens
+    plot(xaxis,AvgAccVP(1:bins),'Color', VP,'linewidth',1.5); %vp
+    %shuffled
+    plot(xaxis,AvgAccNAcShuff(1:bins),'Color', 'k','linewidth',1.5);
+    plot(xaxis,AvgAccVPShuff(1:bins),'Color', 'k','linewidth',1.5);
+    %error
+    patch([xaxis,xaxis(end:-1:1)],[upSEMNAc,downSEMNAc(end:-1:1)],NAc,'EdgeColor','none');alpha(0.5);
+    patch([xaxis,xaxis(end:-1:1)],[upSEMNAcShuff,downSEMNAcShuff(end:-1:1)],'k','EdgeColor','none');alpha(0.5);
+    patch([xaxis,xaxis(end:-1:1)],[upSEMVP,downSEMVP(end:-1:1)],VP,'EdgeColor','none');alpha(0.5);
+    patch([xaxis,xaxis(end:-1:1)],[upSEMVPShuff,downSEMVPShuff(end:-1:1)],'k','EdgeColor','none');alpha(0.5);
+    xlabel('Seconds from reward delivery');
+    ylabel('Accuracy');
+    title('Unit decoding accuracy');
+    legend('NAc units','VP units','Shuffled','Location','northwest');
+    axis([xaxis(1) xaxis(end) 0.47 max(AvgAccVP)+0.04]);
+
+    %find and plot bins exceeding confidence interval
+    for i=1:bins
+        %confidence interval for NAc
+        x = UnitDec.Shuff(NAneurons&selection,i);                      % Create Data
+        SEM = nanstd(x)/sqrt(length(x));               % Standard Error
+        ts = tinv([0.005  0.995],length(x)-1);      % T-Score
+        CI = nanmean(x) + ts*SEM; % Confidence Intervals
+        if nanmean(UnitDec.True(NAneurons&selection,i))>CI(2) NAcConf(1,i)=1; else NAcConf(1,i)=0; end
+        
+        
+        %confidence interval for VP
+        x = UnitDec.Shuff(VPneurons&selection,i);                      % Create Data
+        SEM = nanstd(x)/sqrt(length(x));               % Standard Error
+        ts = tinv([0.005  0.995],length(x)-1);      % T-Score
+        CI = nanmean(x) + ts*SEM; % Confidence Intervals
+        if nanmean(UnitDec.True(VPneurons&selection,i))>CI(2) VPConf(1,i)=1; else VPConf(1,i)=0; end
+        
+    end
+    
+    %find consecutive bins
+    R_2R.UnitNAcComp{e,1}=zeros(1,bins);
+    R_2R.UnitVPComp{e,1}=zeros(1,bins);
+    for i=2:bins-1
+        if NAcConf(1,i)==1
+            if NAcConf(1,i-1)==1 | NAcConf(1,i+1)==1
+                R_2R.UnitNAcComp{e,1}(1,i)=1;
+            end
+        end
+        if VPConf(1,i)==1
+            if VPConf(1,i-1)==1 | VPConf(1,i+1)==1
+                R_2R.UnitVPComp{e,1}(1,i)=1;
+            end
+        end
+    end
+    
+    %check to see if the first bin above shuffled data would be different
+    %if using fewer VP neurons to match NAc
+    for j=1:20
+        for i=1:bins
+            matchedneurons=cat(1,ones(sum(NAneurons&selection),1),zeros(sum(VPneurons&selection)-sum(NAneurons&selection),1));
+            matchedneurons=(matchedneurons(randperm(length(matchedneurons)))==1);
+            %confidence interval for VP
+            VPselectionSh=UnitDec.Shuff(VPneurons&selection,i);
+            VPselection=UnitDec.True(VPneurons&selection,i);
+            x = VPselectionSh(matchedneurons,1);                      % Create Data
+            SEM = nanstd(x)/sqrt(length(x));               % Standard Error
+            ts = tinv([0.005  0.995],length(x)-1);      % T-Score
+            CI = nanmean(x) + ts*SEM; % Confidence Intervals
+            if nanmean(VPselection(matchedneurons,1))>CI(2) R_2R.VPaboveshuff{e,1}(j,i)=1; else R_2R.VPaboveshuff{e,1}(j,i)=0; end
+        end
+    end
+        
+        
+    %plot
+    plot(xaxis,R_2R.UnitVPComp{e,1}*0.48,'*','Color',VP);
+    plot(xaxis,R_2R.UnitNAcComp{e,1}*0.485,'*','Color',NAc);
+    
+    %multiple comparisons for NAc vs VP
+    
+    %make a matrix indicating which region each neuron-bin comes from
+    nbinregion=[];
+    binname=[];
+    rat=[];
+    for i=1:bins %total number of bins being tested
+        nbinregion=cat(2,nbinregion,NAneurons);
+        binname=cat(2,binname,i*ones(length(NAneurons),1));
+        rat=cat(2,rat,R_2R.Ninfo(:,4));
+    end
+    
+    testtrue=UnitDec.True(selection,testbins(1)+1:testbins(2)+1); %becaue the bin 3 is actually the 4th entry, etc
+    testshuff=UnitDec.Shuff(selection,testbins(1)+1:testbins(2)+1);
+    testregion=nbinregion(selection,testbins(1)+1:testbins(2)+1);
+    testbinname=binname(selection,testbins(1)+1:testbins(2)+1);
+    testrat=rat(selection,testbins(1)+1:testbins(2)+1);
+    
+    
+    %find effects of real vs shuffled, region, and bins on accuracy
+    [~,R_2R.UnitDecStatsSubj{e,1},R_2R.UnitDecStatsSubj{e,2}]=anovan(cat(1,testtrue(:),testshuff(:)),{cat(1,zeros(length(testtrue(:)),1),ones(length(testshuff(:)),1)),cat(1,testregion(:),testregion(:)),cat(1,testbinname(:),testbinname(:)),cat(1,testrat(:),testrat(:))},'varnames',{'real vs shuffled','region','bins','subject'},'random',[4],'nested',[0 0 0 0;0 0 0 0;0 0 0 0;0 1 0 0],'display','off','model','full');
+
+%     %do post-hoc comparisons
+%     %this is commented out because can't do post-hoc comparison with nested anova (which we're using for including random effect of subject)
+%     %if you do want to look at this in matlab, need to do the anova without subject (as written out here)
+%     
+%     [~,R_2R.UnitDecStats{e,1},R_2R.UnitDecStats{e,2}]=anovan(cat(1,testtrue(:),testshuff(:)),{cat(1,zeros(length(testtrue(:)),1),ones(length(testshuff(:)),1)),cat(1,testregion(:),testregion(:)),cat(1,testbinname(:),testbinname(:))},'varnames',{'real vs shuffled','region','bins'},'display','off','model','full');
+%     [c,~,~,names]=multcompare(R_2R.UnitDecStats{e,2},'Dimension',[1 2 3],'CType','tukey-kramer','display','off');
+% 
+%     %find post-hoc differences
+%     R_2R.UnitNAcVPComp{e,1}=NaN(2,bins);
+%     for i=1:1+testbins(2)-testbins(1)
+%         %NAc vs VP
+%         Sel=c(:,1)==4*(i-1)+1 & c(:,2)==4*(i-1)+3;
+%         if c(Sel,6)<0.05 R_2R.UnitNAcVPComp{e,1}(1,i+testbins(1))=1; else R_2R.UnitNAcVPComp{e,1}(1,i+testbins(1))=0; end
+%         R_2R.UnitNAcVPComp{e,1}(2,i+testbins(1))=c(Sel,6);
+%     end
+%     
+%     %add it to plot
+%     plot(xaxis,R_2R.UnitNAcVPComp{e,1}(1,:)*(max(AvgAccVP)+0.015),'*','Color','k');
+    
+    %plotting CDF at peak bin
+    subplot(2,3,4)
+    hold on;
+    %NAc
+    [~,NAcbin]=max(AvgAccNAc);
+    [cdfNAc,xNAc] = ecdf(UnitDec.True(NAneurons&selection,NAcbin));
+    [cdfNAcSh,xNAcSh] = ecdf(UnitDec.Shuff(NAneurons&selection,NAcbin));
+    plot(xNAc,cdfNAc,'Color',NAc,'linewidth',1.5);
+    %VP
+    [~,VPbin]=max(AvgAccVP);
+    [cdfVP,xVP] = ecdf(UnitDec.True(VPneurons&selection,VPbin));
+    [cdfVPSh,xVPSh] = ecdf(UnitDec.Shuff(VPneurons&selection,VPbin));
+    plot(xVP,cdfVP,'Color',VP,'linewidth',1.5);
+    %shuffled
+    plot(xNAcSh,cdfNAcSh,'Color','k','linewidth',1.5);
+    xlabel('Decoding accuracy')
+    plot(xVPSh,cdfVPSh,'Color','k','linewidth',1.5);
+    axis([0 1 0 1]);
+    plot([0.5 0.5],[0 1],':','color','k','linewidth',0.75);
+    title(['Accuracy at peak bin: ' num2str(((NAcbin-1)*binint)+(binstart+BinDura(2)/2)) 's NAc, ' num2str(((VPbin-1)*binint)+(binstart+BinDura(2)/2)) 's VP']);
+    xlabel('Decoding accuracy');
+    ylabel('Cumulative fraction of population');
+    legend('NAc units','VP units','Shuffled','Location','northwest');
+
+    %stats comparing peak bins
+    [~,R_2R.UnitDecPeakBinSubj{e,1},~]=anovan(cat(1,UnitDec.True(NAneurons&selection,NAcbin),UnitDec.Shuff(NAneurons&selection,NAcbin),UnitDec.True(VPneurons&selection,VPbin),UnitDec.Shuff(VPneurons&selection,VPbin)),...
+        {cat(1,zeros(sum(NAneurons&selection),1),ones(sum(NAneurons&selection),1),zeros(sum(VPneurons&selection),1),ones(sum(VPneurons&selection),1)),...
+        cat(1,zeros(sum(NAneurons&selection),1),zeros(sum(NAneurons&selection),1),ones(sum(VPneurons&selection),1),ones(sum(VPneurons&selection),1)),...
+        cat(1,R_2R.Ninfo(NAneurons&selection,4),R_2R.Ninfo(NAneurons&selection,4),R_2R.Ninfo(VPneurons&selection,4),R_2R.Ninfo(VPneurons&selection,4))},'varnames',{'real vs shuffled','region','subject'},'random',[3],'nested',[0 0 0;0 0 0;0 1 0],'display','off','model','full');
+%% pooled decoding
+
+    %reset matrices for stats analysis
+    A=[];
+    B=[];
+    C=[];
+    D=[];
+    
+     %get average accuracy for each bin
+    for v=1:length(PoolDec{e,1}.True) %each condition (number of neurons used)
+        for i = 1:bins
+            %NAc
+            if length(PoolDec{e,1}.True{v,1})>1 %in case analysis wasn't run because not enough neurons
+                PoolAccNAc{v,1}(1,i)=nanmean(PoolDec{e,1}.True{v,1}(:,i)); %average accuracy NAc
+                PoolSEMNAc{v,1}(1,i)=nanste(PoolDec{e,1}.True{v,1}(:,i),1); %SEM
+                PoolAccNAcShuff{v,1}(1,i)=nanmean(PoolDec{e,1}.Shuff{v,1}(:,i)); %average accuracy NAcShuff
+                PoolSEMNAcShuff{v,1}(1,i)=nanste(PoolDec{e,1}.Shuff{v,1}(:,i),1); %SEM
+            else
+                PoolAccNAc{v,1}(1,i)=NaN;
+                PoolSEMNAc{v,1}(1,i)=NaN;
+                PoolAccNAcShuff{v,1}(1,i)=NaN;
+                PoolSEMNAcShuff{v,1}(1,i)=NaN;
+            end
+            
+            %VP
+            if length(PoolDec{e,2}.True{v,1})>1 %in case analysis wasn't run because not enough neurons
+                PoolAccVP{v,1}(1,i)=nanmean(PoolDec{e,2}.True{v,1}(:,i)); %average accuracy VP
+                PoolSEMVP{v,1}(1,i)=nanste(PoolDec{e,2}.True{v,1}(:,i),1); %SEM
+                PoolAccVPShuff{v,1}(1,i)=nanmean(PoolDec{e,2}.Shuff{v,1}(:,i)); %average accuracy VPshuff
+                PoolSEMVPShuff{v,1}(1,i)=nanste(PoolDec{e,2}.Shuff{v,1}(:,i),1); %SEM
+            else
+                PoolAccVP{v,1}(1,i)=NaN;
+                PoolSEMVP{v,1}(1,i)=NaN;
+                PoolAccVPShuff{v,1}(1,i)=NaN;
+                PoolSEMVPShuff{v,1}(1,i)=NaN;
+            end
+            
+        end
+
+        %find the time of the most accurate bin for each of the repetitions
+        %(Following reward delivery)
+        time0bin=round((((0-BinDura(2)/2)-binstart)/binint));
+        for j = 1:repetitions
+            %NAc
+            if length(PoolDec{e,1}.True{v,1})>1 %in case analysis wasn't run because not enough neurons
+                [~,PeakBinsNAc{v,1}(j,1)]=max(PoolDec{e,1}.True{v,1}(j,time0bin+1:end));
+                PeakBinsNAc{v,1}(j,2)=(PeakBinsNAc{v,1}(j,1)-1)*binint;
+            else
+                PeakBinsNAc{v,1}(j,1)=NaN;
+                PeakBinsNAc{v,1}(j,2)=NaN;
+            end
+            %VP
+            if length(PoolDec{e,2}.True{v,1})>1 %in case analysis wasn't run because not enough neurons
+                [~,PeakBinsVP{v,1}(j,1)]=max(PoolDec{e,2}.True{v,1}(j,time0bin+1:end));
+                PeakBinsVP{v,1}(j,2)=(PeakBinsVP{v,1}(j,1)-1)*binint;
+            else
+                PeakBinsVP{v,1}(j,1)=NaN;
+                PeakBinsVP{v,1}(j,2)=NaN;
+            end
+        end
+
+        %prepare shading
+        PupSEMNAc{v,1}=PoolAccNAc{v,1}+PoolSEMNAc{v,1};
+        PdownSEMNAc{v,1}=PoolAccNAc{v,1}-PoolSEMNAc{v,1};
+        PupSEMVP{v,1}=PoolAccVP{v,1}+PoolSEMVP{v,1};
+        PdownSEMVP{v,1}=PoolAccVP{v,1}-PoolSEMVP{v,1};
+        PupSEMNAcShuff{v,1}=PoolAccNAcShuff{v,1}+PoolSEMNAcShuff{v,1};
+        PdownSEMNAcShuff{v,1}=PoolAccNAcShuff{v,1}-PoolSEMNAcShuff{v,1};
+        PupSEMVPShuff{v,1}=PoolAccVPShuff{v,1}+PoolSEMVPShuff{v,1};
+        PdownSEMVPShuff{v,1}=PoolAccVPShuff{v,1}-PoolSEMVPShuff{v,1};
+
+        %here I plot the main lines for each condition that will have an entry
+        %in the legend (legend goes by order plotted)
+
+        %plotting decoder accuracy over time
+        subplot(2,3,2); %accumbens
+        hold on;
+        plot(xaxis,PoolAccNAc{v,1}(1:bins),'Color', NAcP{v,1},'linewidth',1);
+        %plot(xaxis,PoolAccNAcShuff{v,1}(1:66),'Color', NAcShuff,'linewidth',3);
+        %plot(xaxis,NAcSig{v,1}-0.53,'*','Color', 'k');
+        %patch([xaxis,xaxis(end:-1:1)],[PupSEMNAc{v,1},PdownSEMNAc{v,1}(end:-1:1)],NAc{v,1},'EdgeColor','none');alpha(0.5);
+        %patch([xaxis,xaxis(end:-1:1)],[PupSEMNAcShuff{v,1},PdownSEMNAcShuff{v,1}(end:-1:1)],NAcShuff,'EdgeColor','none');alpha(0.5);
+        xlabel('Seconds post reward delivery');
+        ylabel('Accuracy');
+        title('NAc decoding accuracy');
+        axis([xaxis(1) xaxis(end) 0.4 1]);
+
+        subplot(2,3,3); %VP
+        hold on;
+        plot(xaxis,PoolAccVP{v,1}(1:bins),'Color', VPP{v,1},'linewidth',1);
+        %plot(xaxis,PoolAccVPShuff{v,1}(1:66),'Color', VPShuff,'linewidth',3);
+        %plot(xaxis,VPSig{v,1}-0.53,'*','Color','k');
+        %patch([xaxis,xaxis(end:-1:1)],[PupSEMVP{v,1},PdownSEMVP{v,1}(end:-1:1)],VP{v,1},'EdgeColor','none');alpha(0.5);
+        %patch([xaxis,xaxis(end:-1:1)],[PupSEMVPShuff{v,1},PdownSEMVPShuff{v,1}(end:-1:1)],VPShuff,'EdgeColor','none');alpha(0.5);
+        xlabel('Seconds post reward delivery');
+        ylabel('Accuracy');
+        title('VP decoding accuracy');
+        axis([xaxis(1) xaxis(end) 0.4 1]);
+
+        %comparison at best bin
+        subplot(2,3,5);
+        hold on;
+        [~,NAcbin]=max(PoolAccNAc{v,1});
+        [~,VPbin]=max(PoolAccVP{v,1});
+        
+        if length(PoolDec{e,1}.True{v,1})>1 && length(PoolDec{e,2}.True{v,1})>1
+            %make matrices that will include best bin at all levels
+            A=cat(1,A,PoolDec{e,1}.True{v,1}(:,NAcbin));
+            B=cat(1,B,PoolDec{e,1}.Shuff{v,1}(:,NAcbin));
+            C=cat(1,C,PoolDec{e,2}.True{v,1}(:,VPbin));
+            D=cat(1,D,PoolDec{e,2}.Shuff{v,1}(:,VPbin));
+        end
+
+        errorbar(NumNeurons(v),PoolAccNAc{v,1}(NAcbin),PoolSEMNAc{v,1}(NAcbin),'*','Color', NAcP{v,1},'linewidth',1.5);
+        errorbar(NumNeurons(v),PoolAccVP{v,1}(VPbin),PoolSEMVP{v,1}(VPbin),'*','Color', VPP{v,1},'linewidth',1.5);
+        xlabel('Neurons used');
+        ylabel('Accuracy');
+        title('Accuracy for most predictive bin');
+        legend('NAc','VP','Location','southeast');
+        axis([0 155 0.5 1]); 
+
+
+        %plotting time of most accurate bin
+        subplot(2,3,6);
+        hold on;
+        errorbar(NumNeurons(v),nanmean(PeakBinsNAc{v,1}(:,2)),nanste(PeakBinsNAc{v,1}(:,2),1),'*','Color', NAcP{v,1},'linewidth',1.5);
+        errorbar(NumNeurons(v),nanmean(PeakBinsVP{v,1}(:,2)),nanste(PeakBinsVP{v,1}(:,2),1),'*','Color', VPP{v,1},'linewidth',1.5);
+        xlabel('Neurons used');
+        ylabel('Seconds post reward delivery')
+        title('Time of most accurate bin');
+        %legend('NAc','VP','Location','southeast');
+        axis([0 155 0 2.5]);
+        legend('NAc','VP','location','southeast')
+
+    end %conditions
+
+    %now plot the shading, significance, and the shuffled separately
+    for v=1:length(PoolDec{e,1}.True) %each condition (number of neurons used)
+        %plotting decoder accuracy over time
+        subplot(2,3,2); %accumbens
+        hold on;
+        plot(xaxis,PoolAccNAcShuff{v,1}(1:66),'Color', 'k','linewidth',1);
+        patch([xaxis,xaxis(end:-1:1)],[PupSEMNAc{v,1},PdownSEMNAc{v,1}(end:-1:1)],NAcP{v,1},'EdgeColor','none');alpha(0.5);
+        patch([xaxis,xaxis(end:-1:1)],[PupSEMNAcShuff{v,1},PdownSEMNAcShuff{v,1}(end:-1:1)],'k','EdgeColor','none');alpha(0.5);
+        xlabel('Seconds post reward delivery');
+        title('NAc decoding accuracy');
+        axis([xaxis(1) xaxis(end) 0.4 1]);
+        legend('10 units','25 units','50 units','100 units','150 units','Shuffled','location','northwest')    
+
+        subplot(2,3,3); %VP
+        hold on;
+        plot(xaxis,PoolAccVPShuff{v,1}(1:66),'Color', 'k','linewidth',1);
+        patch([xaxis,xaxis(end:-1:1)],[PupSEMVP{v,1},PdownSEMVP{v,1}(end:-1:1)],VPP{v,1},'EdgeColor','none');alpha(0.5);
+        patch([xaxis,xaxis(end:-1:1)],[PupSEMVPShuff{v,1},PdownSEMVPShuff{v,1}(end:-1:1)],'k','EdgeColor','none');alpha(0.5);
+        xlabel('Seconds post reward delivery');
+        title('VP decoding accuracy');
+        axis([xaxis(1) xaxis(end) 0.4 1]);
+        legend('10 units','25 units','50 units','100 units','150 units','Shuffled','location','northwest')
+        
+    end
+
+    %stats!
+    %ANOVA for effects of number of neurons and region on time of peak accuracy
+    %time
+    PeakTimes=[];
+    for i=1:length(NumNeurons)
+        if length(PoolDec{e,1}.True{i,1})>1 && length(PoolDec{e,2}.True{i,1})>1
+            PeakTimes(1:repetitions,i)=PeakBinsNAc{i,1}(:,2);
+            PeakTimes(1+repetitions:repetitions+repetitions,i)=PeakBinsVP{i,1}(:,2);
+        end
+    end
+    [~,R_2R.PoolDecTimeStats{e,1},R_2R.PoolDecTimeStats{e,2}]=anova2(PeakTimes,repetitions,'off'); %columns is ensemble size, rows is region
+
+    %ANOVA for effect of ensemble size, region, and shuff vs true on accuracy
+    %of peak bin
+    size=[];
+    for i=1:length(NumNeurons)
+        if length(PoolDec{e,1}.True{i,1})>1 && length(PoolDec{e,2}.True{i,1})>1
+            size=cat(1,size,i*ones(repetitions,1));
+        end
+    end
+    
+    %with shuff vs true
+%     size2=cat(1,size,size,size,size);
+%     shuffd=cat(1,ones(length(A),1),zeros(length(B),1),ones(length(C),1),zeros(length(D),1));
+%     region=cat(1,ones(length(A)+length(B),1),zeros(length(C)+length(D),1));
+%     [~,R_2R.AccAnova,~]=anovan(cat(1,A,B,C,D),{shuffd,region,size2},'varnames',{'real vs shuffled','region','ens size'},'display','off','model','full');
+
+    %without shuff vs true
+    size3=cat(1,size,size);
+    region3=cat(1,ones(length(A),1),zeros(length(C),1));
+    [~,R_2R.PoolDecAccStats{e,1},R_2R.PoolDecAccStats{e,2}]=anovan(cat(1,A,C),{region3,size3},'varnames',{'region','ens size'},'display','off','model','full');
+    
+    
+end %selections
+
+save('R_2R.mat','R_2R');

+ 523 - 0
MatlabScripts/h_Fig4RewardHistory.m

@@ -0,0 +1,523 @@
+%Looking at the average firing rate for a given window in each of 4
+%current/previous reward conditions
+
+clear all;
+load ('R_2R.mat');
+load ('RAW.mat');
+
+%run linear model and stats? 1 is yes, 0 is no
+runanalysis=0;
+
+%divide neurons up by region
+NAneurons=strcmp(R_2R.Ninfo(:,3),'NA');
+VPneurons=strcmp(R_2R.Ninfo(:,3),'VP');
+
+%get parameters
+trialsback=6; %how many trials back to look
+Baseline=R_2R.Param.BinBase; %For normalizing activity
+Window=[0.8 1.3]; %what window of activity is analyzed
+BinDura=R_2R.Param.BinDura;
+bins=R_2R.Param.bins;
+binint=R_2R.Param.binint;
+binstart=R_2R.Param.binstart;
+
+%sorting bin -- which bin the neurons' activity is sorted on for heatmap(in seconds)
+SortBinTime=1; %seconds
+SortBin=(((SortBinTime-BinDura(2)/2)-binstart)/binint); %convert to bin name
+
+%reset
+NN=0;
+EvMeanz=0;
+
+if runanalysis==1  
+    for i=1:length(RAW) %loops through sessions
+        if strcmp('NA',RAW(i).Type(1:2)) | strcmp('VP',RAW(i).Type(1:2))
+            %events
+            EV1=strmatch('RD1P2',RAW(i).Einfo(:,2),'exact');
+            EV2=strmatch('RD1P1',RAW(i).Einfo(:,2),'exact');
+            EV3=strmatch('RD2P2',RAW(i).Einfo(:,2),'exact');
+            EV4=strmatch('RD2P1',RAW(i).Einfo(:,2),'exact');
+            RD=strmatch('RD',RAW(i).Einfo(:,2),'exact');
+            R1=strmatch('Reward1Deliv',RAW(i).Einfo(:,2),'exact');
+            R2=strmatch('Reward2Deliv',RAW(i).Einfo(:,2),'exact');
+
+            %% linear model for impact of previous rewards
+            %reset
+            X=[];
+            Y=[];
+
+            %set up the matrix with outcome identity for each session
+            rewards1=cat(2,RAW(i).Erast{R1,1}(:,1),ones(length(RAW(i).Erast{R1,1}(:,1)),1));
+            rewards2=cat(2,RAW(i).Erast{R2,1}(:,1),zeros(length(RAW(i).Erast{R2,1}(:,1)),1));
+            rewards=cat(1,rewards1,rewards2);
+            [rewards(:,1),ind]=sort(rewards(:,1));
+            rewards(:,2)=rewards(ind,2);
+
+            for k=trialsback+1:length(RAW(i).Erast{RD,1}(:,1))
+                time=RAW(i).Erast{RD,1}(k,1);
+                entry=find(rewards(:,1)==time);
+                for m=1:trialsback+1
+                    X(k-trialsback,m)=rewards(entry+1-m,2);
+                end
+            end
+
+            for j= 1:length(RAW(i).Nrast) %Number of neurons within sessions
+
+                NN=NN+1;
+                rewspk=0;
+                basespk=0;
+
+                %get mean baseline firing for all reward trials
+                [Bcell1,B1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{RD},Baseline,{2});% makes trial by trial rasters for baseline
+                for y= 1:B1n
+                    basespk(1,y)=sum(Bcell1{1,y}>Baseline(1));
+                end
+
+                Bhz=basespk/(Baseline(1,2)-Baseline(1,1));
+                Bmean=nanmean(Bhz);
+                Bstd=nanstd(Bhz);
+
+                %get trial by trial firing rate for all reward trials
+                [EVcell,EVn]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{RD},Window,{2});% makes trial by trial rasters for event
+                for y= 1:EVn
+                    rewspk(y,1)=sum(EVcell{1,y}>Window(1));
+                end       
+                Y=((rewspk(trialsback+1:end,1)/(Window(1,2)-Window(1,1)))-Bmean)/Bstd; %normalize the activity to baseline
+
+                %true data
+                linmdl{NN,1}=fitlm(X,Y);
+                R_2R.RewHist.LinMdlWeights(NN,1:trialsback+1)=linmdl{NN,1}.Coefficients.Estimate(2:trialsback+2)';
+                R_2R.RewHist.LinMdlPVal(NN,1:trialsback+1)=linmdl{NN,1}.Coefficients.pValue(2:trialsback+2)';
+
+                %shuffled
+                YSh=Y(randperm(length(Y)));
+                linmdlSh{NN,1}=fitlm(X,YSh);
+                R_2R.RewHist.LinMdlWeightsSh(NN,1:trialsback+1)=linmdlSh{NN,1}.Coefficients.Estimate(2:trialsback+2)';
+                R_2R.RewHist.LinMdlPValSh(NN,1:trialsback+1)=linmdlSh{NN,1}.Coefficients.pValue(2:trialsback+2)';            
+
+                %% stats comparing effect of current and previous reward
+                %resetting
+                Bcell=0;
+                EV1spk=0;
+                EV2spk=0;
+                EV3spk=0;
+                EV4spk=0;
+                EV1norm=0;
+                EV2norm=0;
+                EV3norm=0;
+                EV4norm=0;
+
+                %get mean baseline firing for all reward trials
+                [Bcell1,B1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV1},Baseline,{2});% makes trial by trial rasters for baseline
+                for y= 1:B1n
+                    Bcell(1,y)=sum(Bcell1{1,y}>Baseline(1));
+                end
+                [Bcell2,B2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV2},Baseline,{2});% makes trial by trial rasters for baseline
+                for z= 1:B2n
+                    Bcell(1,z+B1n)=sum(Bcell2{1,z}>Baseline(1));
+                end
+                [Bcell3,B3n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV3},Baseline,{2});% makes trial by trial rasters for baseline
+                for z= 1:B3n
+                    Bcell(1,z+B1n+B2n)=sum(Bcell3{1,z}>Baseline(1));
+                end
+                [Bcell4,B4n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV4},Baseline,{2});% makes trial by trial rasters for baseline
+                for z= 1:B4n
+                    Bcell(1,z+B1n+B2n+B3n)=sum(Bcell4{1,z}>Baseline(1));
+                end
+
+                Bhz=Bcell/(Baseline(1,2)-Baseline(1,1));
+                Bmean=nanmean(Bhz);
+                Bstd=nanstd(Bhz);
+
+                %get trial by trial firing rate for suc/mal trials
+                [EV1cell,EV1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV1},Window,{2});% makes trial by trial rasters for event
+                for y= 1:EV1n
+                    EV1spk(1,y)=sum(EV1cell{1,y}>Window(1));
+                end       
+                EV1hz=EV1spk/(Window(1,2)-Window(1,1));
+                for x= 1:EV1n
+                    EV1norm(1,x)=((EV1hz(1,x)-Bmean)/Bstd);
+                end
+
+                %get trial by trial firing rate for suc/suc trials
+                [EV2cell,EV2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV2},Window,{2});% makes trial by trial rasters for event
+                for y= 1:EV2n
+                    EV2spk(1,y)=sum(EV2cell{1,y}>Window(1));
+                end       
+                EV2hz=EV2spk/(Window(1,2)-Window(1,1));
+                for x= 1:EV2n
+                    EV2norm(1,x)=((EV2hz(1,x)-Bmean)/Bstd);
+                end
+
+                %get trial by trial firing rate for mal/mal trials
+                [EV3cell,EV3n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV3},Window,{2});% makes trial by trial rasters for event
+                for y= 1:EV3n
+                    EV3spk(1,y)=sum(EV3cell{1,y}>Window(1));
+                end       
+                EV3hz=EV3spk/(Window(1,2)-Window(1,1));
+                for x= 1:EV3n
+                    EV3norm(1,x)=((EV3hz(1,x)-Bmean)/Bstd);
+                end
+
+                %get trial by trial firing rate for mal/suc trials
+                [EV4cell,EV4n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV4},Window,{2});% makes trial by trial rasters for event
+                for y= 1:EV4n
+                    EV4spk(1,y)=sum(EV4cell{1,y}>Window(1));
+                end       
+                EV4hz=EV4spk/(Window(1,2)-Window(1,1));
+                for x= 1:EV4n
+                    EV4norm(1,x)=((EV4hz(1,x)-Bmean)/Bstd);
+                end
+
+                EvMeanz(NN,1)=nanmean(EV1norm);
+                EvMeanz(NN,2)=nanmean(EV2norm);
+                EvMeanz(NN,3)=nanmean(EV3norm);
+                EvMeanz(NN,4)=nanmean(EV4norm);
+
+                R_2R.RewHist.PrevRewMeanz=EvMeanz;
+
+                fprintf('Neuron # %d\n',NN);
+            end
+        end
+
+    end
+end
+
+%% which neurons to look at for stats and plotting?
+
+% Sel=R_2R.SucN | R_2R.MalN; %only look at reward-selective neurons
+Sel=NAneurons | VPneurons; %look at all neurons
+%Sel=R_2R.RewHist.LinMdlPVal(:,2)<0.1; %only neurons with significant impact of previous trial
+%Sel=R_2R.RewHist.LinMdlWeights(:,2)<-1; %only neurons with strong negative impact of previous trial
+
+%% ANOVAs
+
+%setup and run ANOVA for effects of current reward, previous reward, and region on reward firing
+CurrRew=cat(2,zeros(length(NAneurons),2),ones(length(NAneurons),2));
+PrevRew=cat(2,zeros(length(NAneurons),1),ones(length(NAneurons),1),zeros(length(NAneurons),1),ones(length(NAneurons),1));
+Region=cat(2,NAneurons,NAneurons,NAneurons,NAneurons);
+Rat=cat(2,R_2R.Ninfo(:,4),R_2R.Ninfo(:,4),R_2R.Ninfo(:,4),R_2R.Ninfo(:,4));
+
+%to look at a specific selection of cells
+EvMeanz=R_2R.RewHist.PrevRewMeanz(Sel,:);
+CurrRew=CurrRew(Sel,:);
+PrevRew=PrevRew(Sel,:);
+Region=Region(Sel,:);
+Rat=Rat(Sel,:);
+
+%each region individually
+[~,R_2R.RewHist.PrevRewStatsVPSubj{1,1},R_2R.RewHist.PrevRewStatsVPSubj{1,2}]=anovan(reshape(EvMeanz(VPneurons,:),[sum(VPneurons)*4 1]),{reshape(CurrRew(VPneurons,:),[sum(VPneurons)*4 1]),reshape(PrevRew(VPneurons,:),[sum(VPneurons)*4 1]),reshape(Rat(VPneurons,:),[sum(VPneurons)*4 1])},'varnames',{'Current Reward','Previous Reward','Rat'},'random',[3],'display','off','model','full');
+[~,R_2R.RewHist.PrevRewStatsNASubj{1,1},R_2R.RewHist.PrevRewStatsNASubj{1,2}]=anovan(reshape(EvMeanz(NAneurons,:),[sum(NAneurons)*4 1]),{reshape(CurrRew(NAneurons,:),[sum(NAneurons)*4 1]),reshape(PrevRew(NAneurons,:),[sum(NAneurons)*4 1]),reshape(Rat(NAneurons,:),[sum(NAneurons)*4 1])},'varnames',{'Current Reward','Previous Reward','Rat'},'random',[3],'display','off','model','full');
+
+%region comparison
+[~,R_2R.RewHist.PrevRewStatsSubj{1,1},R_2R.RewHist.PrevRewStatsSubj{1,2}]=anovan(EvMeanz(:),{CurrRew(:),PrevRew(:),Region(:),Rat(:)},'varnames',{'Current Reward','Previous Reward','Region','Rat'},'random',[4],'nested',[0 0 0 0;0 0 0 0;0 0 0 0;0 0 1 0],'display','off','model','full');
+
+%setup and run ANOVA for effects of shuffle, trial, and region on coefficient
+Trial=[];
+Region=[];
+Rat=[];
+for i=1:trialsback+1
+    Trial=cat(2,Trial,i*ones(length(NAneurons),1));
+    Region=cat(2,Region,NAneurons);
+    Rat=cat(2,Rat,R_2R.Ninfo(:,4));
+end
+Trial=cat(2,Trial,Trial);
+Region=cat(2,Region,Region);
+Rat=cat(2,Rat,Rat);
+Shuffd=cat(2,zeros(length(NAneurons),trialsback+1),ones(length(NAneurons),trialsback+1));
+Coeffs=cat(2,R_2R.RewHist.LinMdlWeights(:,1:trialsback+1),R_2R.RewHist.LinMdlWeightsSh(:,1:trialsback+1));
+
+[~,R_2R.RewHist.LinCoeffStatsVPSubj{1,1},R_2R.RewHist.LinCoeffStatsVPSubj{1,2}]=anovan(reshape(Coeffs(VPneurons,:),[sum(VPneurons)*2*(trialsback+1) 1]),{reshape(Shuffd(VPneurons,:),[sum(VPneurons)*2*(trialsback+1) 1]),reshape(Trial(VPneurons,:),[sum(VPneurons)*2*(trialsback+1) 1]),reshape(Rat(VPneurons,:),[sum(VPneurons)*2*(trialsback+1) 1])},'varnames',{'Shuffd','Trial','Subject'},'random',[3],'display','off','model','full');
+[~,R_2R.RewHist.LinCoeffStatsNASubj{1,1},R_2R.RewHist.LinCoeffStatsNASubj{1,2}]=anovan(reshape(Coeffs(NAneurons,:),[sum(NAneurons)*2*(trialsback+1) 1]),{reshape(Shuffd(NAneurons,:),[sum(NAneurons)*2*(trialsback+1) 1]),reshape(Trial(NAneurons,:),[sum(NAneurons)*2*(trialsback+1) 1]),reshape(Rat(NAneurons,:),[sum(NAneurons)*2*(trialsback+1) 1])},'varnames',{'Shuffd','Trial','Subject'},'random',[3],'display','off','model','full');
+%% plotting
+Xaxis=[-0.5 3];
+Ishow=find(R_2R.Param.Tm>=Xaxis(1) & R_2R.Param.Tm<=Xaxis(2));
+time1=R_2R.Param.Tm(Ishow);
+
+%color map
+[magma,inferno,plasma,viridis]=colormaps;
+colormap(plasma);
+c=[-100 2000];ClimE=sign(c).*abs(c).^(1/4);%colormap
+
+%colors
+sucrose=[.95  0.55  0.15];
+maltodextrin=[.9  0.3  .9];
+water=[0.00  0.75  0.75];
+total=[0.3  0.1  0.8];
+inh=[0.1 0.021154 0.6];
+exc=[0.9 0.75 0.205816];
+NAc=[0.5  0.1  0.8];
+VP=[0.3  0.7  0.1];
+
+%extra colors to make a gradient
+sucrosem=[.98  0.8  0.35];
+sucrosel=[1  1  0.4];
+maltodextrinm=[1  0.75  1];
+maltodextrinl=[1  0.8  1];
+
+RD1P1=strcmp('RD1P1', R_2R.Erefnames);
+RD1P2=strcmp('RD1P2', R_2R.Erefnames);
+RD2P1=strcmp('RD2P1', R_2R.Erefnames);
+RD2P2=strcmp('RD2P2', R_2R.Erefnames);
+
+%% Get mean firing according to previous trial and then plot
+
+%NAc
+
+%plot suc after suc
+psth1=nanmean(R_2R.Ev(RD1P1).PSTHz(Sel&NAneurons,Ishow),1); 
+sem1=nanste(R_2R.Ev(RD1P1).PSTHz(Sel&NAneurons,Ishow),1); %calculate standard error of the mean
+up1=psth1+sem1;
+down1=psth1-sem1;
+
+%plot suc after malt
+psth2=nanmean(R_2R.Ev(RD1P2).PSTHz(Sel&NAneurons,Ishow),1); 
+sem2=nanste(R_2R.Ev(RD1P2).PSTHz(Sel&NAneurons,Ishow),1); %calculate standard error of the mean
+up2=psth2+sem2;
+down2=psth2-sem2;
+
+
+%plot malt after suc
+psth3=nanmean(R_2R.Ev(RD2P1).PSTHz(Sel&NAneurons,Ishow),1); 
+sem3=nanste(R_2R.Ev(RD2P1).PSTHz(Sel&NAneurons,Ishow),1); %calculate standard error of the mean
+up3=psth3+sem3;
+down3=psth3-sem3;
+
+%plot malt after malt
+psth4=nanmean(R_2R.Ev(RD2P2).PSTHz(Sel&NAneurons,Ishow),1); 
+sem4=nanste(R_2R.Ev(RD2P2).PSTHz(Sel&NAneurons,Ishow),1); %calculate standard error of the mean
+up4=psth4+sem4;
+down4=psth4-sem4;
+
+%make the plot
+subplot(2,4,1);
+hold on;
+title(['Reward response, current/prev reward'])
+plot(time1,psth2,'Color',sucrosem,'linewidth',1);
+plot(time1,psth1,'Color',sucrose,'linewidth',1);
+plot(time1,psth4,'Color',maltodextrinm,'linewidth',1);
+plot(time1,psth3,'Color',maltodextrin,'linewidth',1);
+patch([time1,time1(end:-1:1)],[up2,down2(end:-1:1)],sucrosem,'EdgeColor','none');alpha(0.5);
+patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
+patch([time1,time1(end:-1:1)],[up3,down3(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
+patch([time1,time1(end:-1:1)],[up4,down4(end:-1:1)],maltodextrinm,'EdgeColor','none');alpha(0.5);
+plot([-2 5],[0 0],':','color','k','linewidth',0.75);
+plot([0 0],[-2 8],':','color','k','linewidth',0.75);
+plot([Window(1) Window(1)],[-2 8],'color','b','linewidth',0.85);
+plot([Window(2) Window(2)],[-2 8],'color','b','linewidth',0.85);
+axis([Xaxis(1) Xaxis(2) min(down3)-0.15 max(up2)+0.2]);
+ylabel('Mean firing (z-score)');
+xlabel('Seconds from RD');
+legend('Suc/mal','Suc/suc','Mal/mal','Mal/suc','location','northeast');
+
+if cell2mat(R_2R.RewHist.PrevRewStatsNASubj{1,1}(3,7))<0.05
+    plot(Window(1)+(Window(2)-Window(1))/2,max(up2)+0.1,'*','color','k','markersize',13);
+end
+
+%VP
+
+%plot suc after suc
+psth1=nanmean(R_2R.Ev(RD1P1).PSTHz(Sel&VPneurons,Ishow),1); 
+sem1=nanste(R_2R.Ev(RD1P1).PSTHz(Sel&VPneurons,Ishow),1); %calculate standard error of the mean
+up1=psth1+sem1;
+down1=psth1-sem1;
+
+%plot suc after malt
+psth2=nanmean(R_2R.Ev(RD1P2).PSTHz(Sel&VPneurons,Ishow),1); 
+sem2=nanste(R_2R.Ev(RD1P2).PSTHz(Sel&VPneurons,Ishow),1); %calculate standard error of the mean
+up2=psth2+sem2;
+down2=psth2-sem2;
+
+%plot malt after suc
+psth3=nanmean(R_2R.Ev(RD2P1).PSTHz(Sel&VPneurons,Ishow),1); 
+sem3=nanste(R_2R.Ev(RD2P1).PSTHz(Sel&VPneurons,Ishow),1); %calculate standard error of the mean
+up3=psth3+sem3;
+down3=psth3-sem3;
+
+%plot malt after malt
+psth4=nanmean(R_2R.Ev(RD2P2).PSTHz(Sel&VPneurons,Ishow),1); 
+sem4=nanste(R_2R.Ev(RD2P2).PSTHz(Sel&VPneurons,Ishow),1); %calculate standard error of the mean
+up4=psth4+sem4;
+down4=psth4-sem4;
+
+%make the plot
+subplot(2,4,5);
+hold on;
+title(['Reward response, current/prev reward'])
+plot(time1,psth2,'Color',sucrosem,'linewidth',1);
+plot(time1,psth1,'Color',sucrose,'linewidth',1);
+plot(time1,psth4,'Color',maltodextrinm,'linewidth',1);
+plot(time1,psth3,'Color',maltodextrin,'linewidth',1);
+patch([time1,time1(end:-1:1)],[up2,down2(end:-1:1)],sucrosem,'EdgeColor','none');alpha(0.5);
+patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
+patch([time1,time1(end:-1:1)],[up3,down3(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
+patch([time1,time1(end:-1:1)],[up4,down4(end:-1:1)],maltodextrinm,'EdgeColor','none');alpha(0.5);
+plot([-2 5],[0 0],':','color','k','linewidth',0.75);
+plot([0 0],[-2 8],':','color','k','linewidth',0.75);
+plot([Window(1) Window(1)],[-2 8],'color','b','linewidth',0.85);
+plot([Window(2) Window(2)],[-2 8],'color','b','linewidth',0.85);
+axis([Xaxis(1) Xaxis(2) min(down3)-0.3 max(up2)+0.3]);
+ylabel('Mean firing (z-score)');
+xlabel('Seconds from RD');
+legend('Suc/mal','Suc/suc','Mal/mal','Mal/suc','location','northeast');
+
+if cell2mat(R_2R.RewHist.PrevRewStatsVPSubj{1,1}(3,7))<0.05
+    plot(Window(1)+(Window(2)-Window(1))/2,max(up2)+0.1,'*','color','k','markersize',13);
+end
+
+%% Plotting heatmaps of each condition
+
+%NAc
+
+%sucrose responses following sucrose
+subplot(2,24,[10 11]); %heatmap of water preferring neurons' response to water
+Neurons=R_2R.Ev(RD1P1).PSTHz(Sel&NAneurons,Ishow); %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(time1,[1,sum(Sel&NAneurons,1)],Neurons,ClimE);
+title(['Suc after suc']);
+xlabel('Seconds from RD');
+hold on;
+plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
+
+%following maltodextrin
+subplot(2,24,[7 8]); %heatmap of sucrose preferring neurons' response to maltodextrin
+Neurons=R_2R.Ev(RD1P2).PSTHz(Sel&NAneurons,Ishow); %get the firing rates of neurons of interest
+Neurons=Neurons(SORTimg,:); %sort the neurons same as before
+imagesc(time1,[1,sum(Sel&NAneurons,1)],Neurons,ClimE);
+title(['Suc after mal']);
+ylabel('All neurons plotted individually');
+hold on;
+plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
+
+%maltodextrin responses following sucrose
+subplot(2,24,[16 17]); %heatmap of water preferring neurons' response to water
+Neurons=R_2R.Ev(RD2P1).PSTHz(Sel&NAneurons,Ishow); %get the firing rates of neurons of interest
+Neurons=Neurons(SORTimg,:); %sort the neurons as before
+imagesc(time1,[1,sum(Sel&NAneurons,1)],Neurons,ClimE);
+title(['Mal after suc']);
+hold on;
+plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
+
+%following maltodextrin
+subplot(2,24,[13 14]); %heatmap of sucrose preferring neurons' response to maltodextrin
+Neurons=R_2R.Ev(RD2P2).PSTHz(Sel&NAneurons,Ishow); %get the firing rates of neurons of interest
+Neurons=Neurons(SORTimg,:); %sort the neurons same as before
+imagesc(time1,[1,sum(Sel&NAneurons,1)],Neurons,ClimE);
+title(['Mal after mal']);
+hold on;
+plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
+
+%VP
+
+%sucrose responses following sucrose
+subplot(2,24,[34 35]); %heatmap of water preferring neurons' response to water
+Neurons=R_2R.Ev(RD1P1).PSTHz(Sel&VPneurons,Ishow); %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(time1,[1,sum(Sel&VPneurons,1)],Neurons,ClimE);
+title(['Suc after suc']);
+xlabel('Seconds from RD');
+hold on;
+plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
+
+%following maltodextrin
+subplot(2,24,[31 32]); %heatmap of sucrose preferring neurons' response to maltodextrin
+Neurons=R_2R.Ev(RD1P2).PSTHz(Sel&VPneurons,Ishow); %get the firing rates of neurons of interest
+Neurons=Neurons(SORTimg,:); %sort the neurons same as before
+imagesc(time1,[1,sum(Sel&VPneurons,1)],Neurons,ClimE);
+title(['Suc after mal']);
+ylabel('All neurons plotted individually');
+hold on;
+plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
+
+%maltodextrin responses following sucrose
+subplot(2,24,[40 41]); %heatmap of water preferring neurons' response to water
+Neurons=R_2R.Ev(RD2P1).PSTHz(Sel&VPneurons,Ishow); %get the firing rates of neurons of interest
+Neurons=Neurons(SORTimg,:); %sort the neurons as before
+imagesc(time1,[1,sum(Sel&VPneurons,1)],Neurons,ClimE);
+title(['Mal after suc']);
+hold on;
+plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
+
+%following maltodextrin
+subplot(2,24,[37 38]); %heatmap of sucrose preferring neurons' response to maltodextrin
+Neurons=R_2R.Ev(RD2P2).PSTHz(Sel&VPneurons,Ishow); %get the firing rates of neurons of interest
+Neurons=Neurons(SORTimg,:); %sort the neurons same as before
+imagesc(time1,[1,sum(Sel&VPneurons,1)],Neurons,ClimE);
+title(['Mal after mal']);
+hold on;
+plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
+
+%% plot linear model coefficients
+
+%Plot all neurons
+Sel=NAneurons<2;
+
+%coefficients for each trial
+subplot(2,4,4);
+hold on;
+errorbar(0:trialsback,nanmean(R_2R.RewHist.LinMdlWeights(Sel&NAneurons,1:trialsback+1),1),nanste(R_2R.RewHist.LinMdlWeights(Sel&NAneurons,1:trialsback+1),1),'color',NAc,'linewidth',1);
+errorbar(0:trialsback,nanmean(R_2R.RewHist.LinMdlWeights(Sel&VPneurons,1:trialsback+1),1),nanste(R_2R.RewHist.LinMdlWeights(Sel&VPneurons,1:trialsback+1),1),'color',VP,'linewidth',1);
+errorbar(0:trialsback,nanmean(R_2R.RewHist.LinMdlWeightsSh(Sel&NAneurons,1:trialsback+1),1),nanste(R_2R.RewHist.LinMdlWeights(Sel&NAneurons,1:trialsback+1),1),'color','k','linewidth',1);
+errorbar(0:trialsback,nanmean(R_2R.RewHist.LinMdlWeightsSh(Sel&VPneurons,1:trialsback+1),1),nanste(R_2R.RewHist.LinMdlWeights(Sel&VPneurons,1:trialsback+1),1),'color','k','linewidth',1);
+xlabel('Trials back');
+ylabel('Mean coefficient weight');
+title('Linear model coefficients');
+legend('NAc','VP','Shuff');
+
+axis([-1 trialsback+1 -1 2.8]);
+plot([-1 trialsback+1],[0 0],':','color','k','linewidth',0.75);
+
+%stats to check if VP and NAc are greater than chance
+R_2R.RewHist.LinCoeffMultComp=[];
+[c,~,~,~]=multcompare(R_2R.RewHist.LinCoeffStatsNASubj{1,2},'dimension',[1,2],'display','off');
+[d,~,~,~]=multcompare(R_2R.RewHist.LinCoeffStatsVPSubj{1,2},'dimension',[1,2],'display','off');
+for i=1:trialsback+1
+
+    %NAc vs shuff
+    Cel=c(:,1)==2*(i-1)+1 & c(:,2)==2*(i-1)+2;
+    if c(Cel,6)<0.05 R_2R.RewHist.LinCoeffMultComp(i,1)=1; else R_2R.RewHist.LinCoeffMultComp(i,1)=0; end
+    R_2R.RewHist.LinCoeffMultComp(i,2)=c(Cel,2);
+
+    %VP vs shuff
+    Cel=d(:,1)==2*(i-1)+1 & d(:,2)==2*(i-1)+2;
+    if d(Cel,6)<0.05 R_2R.RewHist.LinCoeffMultComp(i,3)=1; else R_2R.RewHist.LinCoeffMultComp(i,3)=0; end
+    R_2R.RewHist.LinCoeffMultComp(i,4)=d(Cel,4);
+end
+subplot(2,4,4);
+hold on;
+plot([0:trialsback]-0.1,(R_2R.RewHist.LinCoeffMultComp(:,1)-0.5)*5,'*','color',NAc); %VP vs shuff
+plot([0:trialsback]+0.1,(R_2R.RewHist.LinCoeffMultComp(:,3)-0.5)*5,'*','color',VP); %NAc vs shuff
+
+
+%number of neurons with significant weights
+subplot(2,4,8);
+hold on;
+plot(0:trialsback,sum(R_2R.RewHist.LinMdlPVal(Sel&NAneurons,1:trialsback+1)<0.05,1)/sum(Sel&NAneurons),'color',NAc,'linewidth',1);
+plot(0:trialsback,sum(R_2R.RewHist.LinMdlPVal(Sel&VPneurons,1:trialsback+1)<0.05,1)/sum(Sel&VPneurons),'color',VP,'linewidth',1);
+plot(0:trialsback,sum(R_2R.RewHist.LinMdlPValSh(Sel&NAneurons,1:trialsback+1)<0.05,1)/sum(Sel&NAneurons),'color',NAc/3,'linewidth',1);
+plot(0:trialsback,sum(R_2R.RewHist.LinMdlPValSh(Sel&VPneurons,1:trialsback+1)<0.05,1)/sum(Sel&VPneurons),'color',VP/3,'linewidth',1);
+axis([-1 trialsback+1 0 0.5]);
+xlabel('Trials back');
+ylabel('Fraction of the population');
+title('Outcome-modulated neurons');
+
+%Chi squared stat for each trial
+R_2R.RewHist.LinMdlX2all=[];
+R_2R.RewHist.LinMdlX2region=[];
+
+for i=1:trialsback+1
+    [~,R_2R.RewHist.LinMdlX2all(i,1),R_2R.RewHist.LinMdlX2all(i,2)]=crosstab(cat(1,R_2R.RewHist.LinMdlPVal(Sel,i)<0.05,R_2R.RewHist.LinMdlPValSh(Sel,i)<0.05),cat(1,VPneurons,VPneurons+2));
+    [~,R_2R.RewHist.LinMdlX2region(i,1),R_2R.RewHist.LinMdlX2region(i,2)]=crosstab(R_2R.RewHist.LinMdlPVal(Sel,i)<0.05,VPneurons);
+end
+%plot([0:trialsback]-0.1,(R_2R.LinMdlX2all(:,2)<0.05)-0.52,'*','color',NAc);
+plot([0:trialsback],(R_2R.RewHist.LinMdlX2region(:,2)<0.05&R_2R.RewHist.LinMdlX2all(:,2)<0.05)-0.52,'*','color','k');
+
+save('R_2R.mat','R_2R');

+ 624 - 0
MatlabScripts/i_FigS7CuePERewHist.m

@@ -0,0 +1,624 @@
+%look at impact of previous trials on Cue and PE
+
+%Looking at the average firing rate for a given window in each of 4
+%current/previous reward conditions
+
+clear all;
+load ('R_2R.mat');
+load ('RAW.mat');
+
+%run linear model and stats? 1 is yes, 0 is no
+runanalysis=0;
+
+%divide neurons up by region
+NAneurons=strcmp(R_2R.Ninfo(:,3),'NA');
+VPneurons=strcmp(R_2R.Ninfo(:,3),'VP');
+
+%get parameters
+trialsback=6; %how many trials back to look
+PEBaseline=R_2R.Param.BinBase+0.5; %For normalizing activity
+CueBaseline=[-11 -1];
+CueWindow=[0 0.5]; %what window of activity is analyzed
+PEWindow=[-0.6 0.7];
+BinDura=R_2R.Param.BinDura;
+bins=R_2R.Param.bins;
+binint=R_2R.Param.binint;
+binstart=R_2R.Param.binstart;
+
+
+%sorting bin -- which bin the neurons' activity is sorted on for heatmap(in seconds)
+SortBinTime=1; %seconds
+SortBin=(((SortBinTime-BinDura(2)/2)-binstart)/binint); %convert to bin name
+
+%reset
+NN=0;
+PEEvMeanz=0;
+
+
+if runanalysis==1  
+    for i=1:length(RAW) %loops through sessions
+        if strcmp('NA',RAW(i).Type(1:2)) | strcmp('VP',RAW(i).Type(1:2)) %only look at suc v mal sessions
+            %events
+            EV3=strmatch('CueP1',RAW(i).Einfo(:,2),'exact');
+            EV4=strmatch('CueP2',RAW(i).Einfo(:,2),'exact');
+            EV1=strmatch('PEP1',RAW(i).Einfo(:,2),'exact');
+            EV2=strmatch('PEP2',RAW(i).Einfo(:,2),'exact');
+            Cue=strmatch('Cue',RAW(i).Einfo(:,2),'exact');
+            PE=strmatch('PECue',RAW(i).Einfo(:,2),'exact');
+            R1=strmatch('Reward1Deliv',RAW(i).Einfo(:,2),'exact');
+            R2=strmatch('Reward2Deliv',RAW(i).Einfo(:,2),'exact');
+
+            %% linear model for impact of previous rewards
+            %reset
+            Xcue=[];
+            Xpe=[];
+            Y=[];
+
+            %set up the matrix with outcome identity for each session
+            rewards1=cat(2,RAW(i).Erast{R1,1}(:,1),ones(length(RAW(i).Erast{R1,1}(:,1)),1));
+            rewards2=cat(2,RAW(i).Erast{R2,1}(:,1),zeros(length(RAW(i).Erast{R2,1}(:,1)),1));
+            rewards=cat(1,rewards1,rewards2);
+            [rewards(:,1),ind]=sort(rewards(:,1));
+            rewards(:,2)=rewards(ind,2);
+
+            %cue
+            firstcue=find(RAW(i).Erast{Cue,1}==min(RAW(i).Erast{Cue,1}(RAW(i).Erast{Cue,1}(:,1)>rewards(trialsback),1))); %find first cue with at least 5 rewards prior
+            for k=firstcue+1:length(RAW(i).Erast{Cue,1}(:,1))
+                time=RAW(i).Erast{Cue,1}(k,1);
+                entry=find(rewards==max(rewards(rewards(:,1)<time)));
+                for m=1:trialsback
+                    Xcue(k-trialsback,m)=rewards(entry+1-m,2);
+                end
+            end
+
+            %PE
+            for k=trialsback+1:length(RAW(i).Erast{PE,1}(:,1))
+                time=RAW(i).Erast{PE,1}(k,1);
+                entry=find(rewards==max(rewards(rewards(:,1)<time)));
+                for m=1:trialsback
+                    Xpe(k-trialsback,m)=rewards(entry+1-m,2);
+                end
+            end    
+
+            for j= 1:length(RAW(i).Nrast) %Number of neurons within sessions
+
+                NN=NN+1;
+
+                %cue
+                rewspk=0;
+                basespk=0;
+
+                %get mean baseline firing for all cues
+                [Bcell1,B1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{Cue},CueBaseline,{2});% makes trial by trial rasters for baseline
+                for y= 1:B1n
+                    basespk(1,y)=sum(Bcell1{1,y}>CueBaseline(1));
+                end
+
+                Bhz=basespk/(CueBaseline(1,2)-CueBaseline(1,1));
+                Bmean=nanmean(Bhz);
+                Bstd=nanstd(Bhz);
+
+                %get trial by trial firing rate for all cue trials
+                [EVcell,EVn]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{Cue},CueWindow,{2});% makes trial by trial rasters for event
+                for y= 1:EVn
+                    rewspk(y,1)=sum(EVcell{1,y}>CueWindow(1));
+                end       
+                Y=((rewspk(trialsback+1:end,1)/(CueWindow(1,2)-CueWindow(1,1)))-Bmean)/Bstd; %normalize the activity to baseline
+
+                %true data
+                linmdl{NN,1}=fitlm(Xcue,Y);
+                R_2R.RewHist.LinMdlCueWeights(NN,1:trialsback)=linmdl{NN,1}.Coefficients.Estimate(2:trialsback+1)';
+                R_2R.RewHist.LinMdlCuePVal(NN,1:trialsback)=linmdl{NN,1}.Coefficients.pValue(2:trialsback+1)';
+
+                %shuffled
+                YSh=Y(randperm(length(Y)));
+                linmdlSh{NN,1}=fitlm(Xcue,YSh);
+                R_2R.RewHist.LinMdlCueWeightsSh(NN,1:trialsback)=linmdlSh{NN,1}.Coefficients.Estimate(2:trialsback+1)';
+                R_2R.RewHist.LinMdlCuePValSh(NN,1:trialsback)=linmdlSh{NN,1}.Coefficients.pValue(2:trialsback+1)';
+
+                %PE
+                rewspk=0;
+                basespk=0;
+
+                %get mean baseline firing for all PEs
+                [Bcell1,B1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{PE},PEBaseline,{2});% makes trial by trial rasters for baseline
+                for y= 1:B1n
+                    basespk(1,y)=sum(Bcell1{1,y}>PEBaseline(1));
+                end
+
+                Bhz=basespk/(PEBaseline(1,2)-PEBaseline(1,1));
+                Bmean=nanmean(Bhz);
+                Bstd=nanstd(Bhz);
+
+                %get trial by trial firing rate for all PE trials
+                [EVcell,EVn]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{PE},PEWindow,{2});% makes trial by trial rasters for event
+                for y= 1:EVn
+                    rewspk(y,1)=sum(EVcell{1,y}>PEWindow(1));
+                end       
+                Y=((rewspk(trialsback+1:end,1)/(PEWindow(1,2)-PEWindow(1,1)))-Bmean)/Bstd; %normalize the activity to baseline
+
+                %true data
+                linmdl{NN,1}=fitlm(Xpe,Y);
+                R_2R.RewHist.LinMdlPEWeights(NN,1:trialsback)=linmdl{NN,1}.Coefficients.Estimate(2:trialsback+1)';
+                R_2R.RewHist.LinMdlPEPVal(NN,1:trialsback)=linmdl{NN,1}.Coefficients.pValue(2:trialsback+1)';
+
+                %shuffled
+                YSh=Y(randperm(length(Y)));
+                linmdlSh{NN,1}=fitlm(Xpe,YSh);
+                R_2R.RewHist.LinMdlPEWeightsSh(NN,1:trialsback)=linmdlSh{NN,1}.Coefficients.Estimate(2:trialsback+1)';
+                R_2R.RewHist.LinMdlPEPValSh(NN,1:trialsback)=linmdlSh{NN,1}.Coefficients.pValue(2:trialsback+1)';  
+
+                %% stats comparing effect of current and previous reward on PE
+                %resetting
+                Bcell=0;
+                EV1spk=0;
+                EV2spk=0;
+                EV1norm=0;
+                EV2norm=0;
+
+                %get mean baseline firing for all PE trials
+                [Bcell1,B1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV1},PEBaseline,{2});% makes trial by trial rasters for baseline
+                for y= 1:B1n
+                    Bcell(1,y)=sum(Bcell1{1,y}>PEBaseline(1));
+                end
+                [Bcell2,B2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV2},PEBaseline,{2});% makes trial by trial rasters for baseline
+                for z= 1:B2n
+                    Bcell(1,z+B1n)=sum(Bcell2{1,z}>PEBaseline(1));
+                end
+
+                Bhz=Bcell/(PEBaseline(1,2)-PEBaseline(1,1));
+                Bmean=nanmean(Bhz);
+                Bstd=nanstd(Bhz);
+
+                %get trial by trial firing rate for post suc trials
+                [EV1cell,EV1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV1},CueWindow,{2});% makes trial by trial rasters for event
+                for y= 1:EV1n
+                    EV1spk(1,y)=sum(EV1cell{1,y}>CueWindow(1));
+                end       
+                EV1hz=EV1spk/(CueWindow(1,2)-CueWindow(1,1));
+                for x= 1:EV1n
+                    EV1norm(1,x)=((EV1hz(1,x)-Bmean)/Bstd);
+                end
+
+                %get trial by trial firing rate for post mal trials
+                [EV2cell,EV2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV2},CueWindow,{2});% makes trial by trial rasters for event
+                for y= 1:EV2n
+                    EV2spk(1,y)=sum(EV2cell{1,y}>CueWindow(1));
+                end       
+                EV2hz=EV2spk/(CueWindow(1,2)-CueWindow(1,1));
+                for x= 1:EV2n
+                    EV2norm(1,x)=((EV2hz(1,x)-Bmean)/Bstd);
+                end
+
+
+                PEEvMeanz(NN,1)=nanmean(EV1norm);
+                PEEvMeanz(NN,2)=nanmean(EV2norm);
+
+                %% stats comparing effect of current and previous reward on cue
+                %resetting
+                Bcell=0;
+                EV1spk=0;
+                EV2spk=0;
+                EV1norm=0;
+                EV2norm=0;
+
+                %get mean baseline firing for all cue trials
+                [Bcell1,B1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV3},CueBaseline,{2});% makes trial by trial rasters for baseline
+                for y= 1:B1n
+                    Bcell(1,y)=sum(Bcell1{1,y}>CueBaseline(1));
+                end
+                [Bcell2,B2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV4},CueBaseline,{2});% makes trial by trial rasters for baseline
+                for z= 1:B2n
+                    Bcell(1,z+B1n)=sum(Bcell2{1,z}>CueBaseline(1));
+                end
+
+                Bhz=Bcell/(CueBaseline(1,2)-CueBaseline(1,1));
+                Bmean=nanmean(Bhz);
+                Bstd=nanstd(Bhz);
+
+                %get trial by trial firing rate for post suc trials
+                [EV1cell,EV1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV3},CueWindow,{2});% makes trial by trial rasters for event
+                for y= 1:EV1n
+                    EV1spk(1,y)=sum(EV1cell{1,y}>CueWindow(1));
+                end       
+                EV1hz=EV1spk/(CueWindow(1,2)-CueWindow(1,1));
+                for x= 1:EV1n
+                    EV1norm(1,x)=((EV1hz(1,x)-Bmean)/Bstd);
+                end
+
+                %get trial by trial firing rate for post mal trials
+                [EV2cell,EV2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV4},CueWindow,{2});% makes trial by trial rasters for event
+                for y= 1:EV2n
+                    EV2spk(1,y)=sum(EV2cell{1,y}>CueWindow(1));
+                end       
+                EV2hz=EV2spk/(CueWindow(1,2)-CueWindow(1,1));
+                for x= 1:EV2n
+                    EV2norm(1,x)=((EV2hz(1,x)-Bmean)/Bstd);
+                end
+
+
+                CueEvMeanz(NN,1)=nanmean(EV1norm);
+                CueEvMeanz(NN,2)=nanmean(EV2norm);
+
+                fprintf('Neuron # %d\n',NN);
+            end
+        end
+        R_2R.RewHist.PrevRewPEMeanz=PEEvMeanz;
+        R_2R.RewHist.PrevRewCueMeanz=CueEvMeanz;
+
+    end
+end
+
+%% which neurons to look at for stats and plotting?
+
+% Sel=R_2R.SucN | R_2R.MalN; %only look at reward-selective neurons
+Sel=NAneurons | VPneurons; %look at all neurons
+%Sel=R_2R.RewHist.LinMdlPVal(:,2)<0.1; %only neurons with significant impact of previous trial
+%Sel=R_2R.RewHist.LinMdlPEWeights(:,2)<-1; %only neurons with strong negative impact of previous trial
+
+%% ANOVAs
+
+%setup and run ANOVA for effects of current reward, previous reward, and region on reward firing
+PrevRew=cat(2,zeros(length(NAneurons),1),ones(length(NAneurons),1));
+Region=cat(2,NAneurons,NAneurons);
+Rat=cat(2,R_2R.Ninfo(:,4),R_2R.Ninfo(:,4));
+
+%to look at a specific selection of cells
+PEEvMeanz=R_2R.RewHist.PrevRewPEMeanz(Sel,:);
+CueEvMeanz=R_2R.RewHist.PrevRewCueMeanz(Sel,:);
+PrevRew=PrevRew(Sel,:);
+Region=Region(Sel,:);
+Rat=Rat(Sel,:);
+
+%cue
+%each region individually
+[~,R_2R.RewHist.PrevRewStatsCueVPSubj{1,1},R_2R.RewHist.PrevRewStatsCueVPSubj{1,2}]=anovan(reshape(CueEvMeanz(VPneurons,:),[sum(VPneurons)*2 1]),{reshape(PrevRew(VPneurons,:),[sum(VPneurons)*2 1]),reshape(Rat(VPneurons,:),[sum(VPneurons)*2 1])},'varnames',{'Previous Reward','Rat'},'random',[2],'display','off','model','full');
+[~,R_2R.RewHist.PrevRewStatsCueNASubj{1,1},R_2R.RewHist.PrevRewStatsCueNASubj{1,2}]=anovan(reshape(CueEvMeanz(NAneurons,:),[sum(NAneurons)*2 1]),{reshape(PrevRew(NAneurons,:),[sum(NAneurons)*2 1]),reshape(Rat(NAneurons,:),[sum(NAneurons)*2 1])},'varnames',{'Previous Reward','Rat'},'random',[2],'display','off','model','full');
+
+%region comparison
+[~,R_2R.RewHist.PrevRewStatsCueSubj{1,1},R_2R.RewHist.PrevRewStatsCueSubj{1,2}]=anovan(CueEvMeanz(:),{PrevRew(:),Region(:),Rat(:)},'varnames',{'Previous Reward','Region','Rat'},'random',[3],'nested',[0 0 0;0 0 0;0 1 0],'display','off','model','full');
+
+%pe
+%each region individually
+[~,R_2R.RewHist.PrevRewStatsPEVPSubj{1,1},R_2R.RewHist.PrevRewStatsPEVPSubj{1,2}]=anovan(reshape(PEEvMeanz(VPneurons,:),[sum(VPneurons)*2 1]),{reshape(PrevRew(VPneurons,:),[sum(VPneurons)*2 1]),reshape(Rat(VPneurons,:),[sum(VPneurons)*2 1])},'varnames',{'Previous Reward','Rat'},'random',[2],'display','off','model','full');
+[~,R_2R.RewHist.PrevRewStatsPENASubj{1,1},R_2R.RewHist.PrevRewStatsPENASubj{1,2}]=anovan(reshape(PEEvMeanz(NAneurons,:),[sum(NAneurons)*2 1]),{reshape(PrevRew(NAneurons,:),[sum(NAneurons)*2 1]),reshape(Rat(NAneurons,:),[sum(NAneurons)*2 1])},'varnames',{'Previous Reward','Rat'},'random',[2],'display','off','model','full');
+
+%region comparison
+[~,R_2R.RewHist.PrevRewStatsPESubj{1,1},R_2R.RewHist.PrevRewStatsPESubj{1,2}]=anovan(PEEvMeanz(:),{PrevRew(:),Region(:),Rat(:)},'varnames',{'Previous Reward','Region','Rat'},'random',[3],'nested',[0 0 0;0 0 0;0 1 0],'display','off','model','full');
+
+
+%setup and run ANOVA for effects of shuffle, trial, and region on coefficient
+Trial=[];
+Region=[];
+Rat=[];
+for i=1:trialsback
+    Trial=cat(2,Trial,i*ones(length(NAneurons),1));
+    Region=cat(2,Region,NAneurons);
+    Rat=cat(2,Rat,R_2R.Ninfo(:,4));
+end
+Trial=cat(2,Trial,Trial);
+Region=cat(2,Region,Region);
+Rat=cat(2,Rat,Rat);
+Shuffd=cat(2,zeros(length(NAneurons),trialsback),ones(length(NAneurons),trialsback));
+
+%Cue
+Coeffs=cat(2,R_2R.RewHist.LinMdlCueWeights(:,1:trialsback),R_2R.RewHist.LinMdlCueWeightsSh(:,1:trialsback));
+[~,R_2R.RewHist.LinCoeffStatsCueVPSubj{1,1},R_2R.RewHist.LinCoeffStatsCueVPSubj{1,2}]=anovan(reshape(Coeffs(VPneurons,:),[sum(VPneurons)*2*(trialsback) 1]),{reshape(Shuffd(VPneurons,:),[sum(VPneurons)*2*(trialsback) 1]),reshape(Trial(VPneurons,:),[sum(VPneurons)*2*(trialsback) 1]),reshape(Rat(VPneurons,:),[sum(VPneurons)*2*(trialsback) 1])},'varnames',{'Shuffd','Trial','Subject'},'random',[3],'display','off','model','full');
+[~,R_2R.RewHist.LinCoeffStatsCueNASubj{1,1},R_2R.RewHist.LinCoeffStatsCueNASubj{1,2}]=anovan(reshape(Coeffs(NAneurons,:),[sum(NAneurons)*2*(trialsback) 1]),{reshape(Shuffd(NAneurons,:),[sum(NAneurons)*2*(trialsback) 1]),reshape(Trial(NAneurons,:),[sum(NAneurons)*2*(trialsback) 1]),reshape(Rat(NAneurons,:),[sum(NAneurons)*2*(trialsback) 1])},'varnames',{'Shuffd','Trial','Subject'},'random',[3],'display','off','model','full');
+
+
+
+%PE
+Coeffs=cat(2,R_2R.RewHist.LinMdlPEWeights(:,1:trialsback),R_2R.RewHist.LinMdlPEWeightsSh(:,1:trialsback));
+[~,R_2R.RewHist.LinCoeffStatsPEVPSubj{1,1},R_2R.RewHist.LinCoeffStatsPEVPSubj{1,2}]=anovan(reshape(Coeffs(VPneurons,:),[sum(VPneurons)*2*(trialsback) 1]),{reshape(Shuffd(VPneurons,:),[sum(VPneurons)*2*(trialsback) 1]),reshape(Trial(VPneurons,:),[sum(VPneurons)*2*(trialsback) 1]),reshape(Rat(VPneurons,:),[sum(VPneurons)*2*(trialsback) 1])},'varnames',{'Shuffd','Trial','Subject'},'random',[3],'display','off','model','full');
+[~,R_2R.RewHist.LinCoeffStatsPENASubj{1,1},R_2R.RewHist.LinCoeffStatsPENASubj{1,2}]=anovan(reshape(Coeffs(NAneurons,:),[sum(NAneurons)*2*(trialsback) 1]),{reshape(Shuffd(NAneurons,:),[sum(NAneurons)*2*(trialsback) 1]),reshape(Trial(NAneurons,:),[sum(NAneurons)*2*(trialsback) 1]),reshape(Rat(NAneurons,:),[sum(NAneurons)*2*(trialsback) 1])},'varnames',{'Shuffd','Trial','Subject'},'random',[3],'display','off','model','full');
+
+
+%% plotting
+XaxisCue=[-0.5 1];
+XaxisPE=[-1 2];
+Ishow=find(R_2R.Param.Tm>=XaxisCue(1) & R_2R.Param.Tm<=XaxisCue(2));
+Ushow=find(R_2R.Param.Tm>=XaxisPE(1) & R_2R.Param.Tm<=XaxisPE(2));
+time1=R_2R.Param.Tm(Ishow);
+time2=R_2R.Param.Tm(Ushow);
+
+%color map
+[magma,inferno,plasma,viridis]=colormaps;
+colormap(plasma);
+c=[-100 2000];ClimE=sign(c).*abs(c).^(1/4);%colormap
+
+%colors
+sucrose=[.95  0.55  0.15];
+maltodextrin=[.9  0.3  .9];
+water=[0.00  0.75  0.75];
+total=[0.3  0.1  0.8];
+inh=[0.1 0.021154 0.6];
+exc=[0.9 0.75 0.205816];
+NAc=[0.5  0.1  0.8];
+VP=[0.3  0.7  0.1];
+
+%extra colors to make a gradient
+sucrosem=[.98  0.8  0.35];
+sucrosel=[1  1  0.4];
+maltodextrinm=[1  0.75  1];
+maltodextrinl=[1  0.8  1];
+
+RD1P1=strcmp('CueP1', R_2R.Erefnames);
+RD1P2=strcmp('CueP2', R_2R.Erefnames);
+RD2P1=strcmp('PEP1', R_2R.Erefnames);
+RD2P2=strcmp('PEP2', R_2R.Erefnames);
+
+%% Get mean firing according to previous trial and then plot
+
+%NAc
+
+%plot suc after suc
+psth1=nanmean(R_2R.Ev(RD1P1).PSTHz(Sel&NAneurons,Ishow),1); 
+sem1=nanste(R_2R.Ev(RD1P1).PSTHz(Sel&NAneurons,Ishow),1); %calculate standard error of the mean
+up1=psth1+sem1;
+down1=psth1-sem1;
+
+%plot suc after malt
+psth2=nanmean(R_2R.Ev(RD1P2).PSTHz(Sel&NAneurons,Ishow),1); 
+sem2=nanste(R_2R.Ev(RD1P2).PSTHz(Sel&NAneurons,Ishow),1); %calculate standard error of the mean
+up2=psth2+sem2;
+down2=psth2-sem2;
+
+
+
+
+%make the plot
+subplot(2,4,1);
+hold on;
+title(['Cue response, NAc'])
+plot(time1,psth1,'Color',sucrose,'linewidth',1);
+plot(time1,psth2,'Color',maltodextrin,'linewidth',1);
+
+patch([time1,time1(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
+patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
+
+plot([-2 5],[0 0],':','color','k','linewidth',0.75);
+plot([0 0],[-2 8],':','color','k','linewidth',0.75);
+plot([CueWindow(1) CueWindow(1)],[-2 8],'color','b','linewidth',0.85);
+plot([CueWindow(2) CueWindow(2)],[-2 8],'color','b','linewidth',0.85);
+axis([XaxisCue(1) XaxisCue(2) min(down2)-0.15 max(up1)+0.2]);
+ylabel('Mean firing (z-score)');
+xlabel('Seconds from cue');
+legend('Post-suc','Post-mal','location','northeast');
+
+if cell2mat(R_2R.RewHist.PrevRewStatsCueNASubj{1,1}(2,7))<0.05
+    plot(CueWindow(1)+(CueWindow(2)-CueWindow(1))/2,max(up2)+0.1,'*','color','k','markersize',13);
+end
+
+
+subplot(2,4,5);
+hold on;
+title('PE response, NAc');
+%plot malt after suc
+psth3=nanmean(R_2R.Ev(RD2P1).PSTHz(Sel&NAneurons,Ushow),1); 
+sem3=nanste(R_2R.Ev(RD2P1).PSTHz(Sel&NAneurons,Ushow),1); %calculate standard error of the mean
+up3=psth3+sem3;
+down3=psth3-sem3;
+
+%plot malt after malt
+psth4=nanmean(R_2R.Ev(RD2P2).PSTHz(Sel&NAneurons,Ushow),1); 
+sem4=nanste(R_2R.Ev(RD2P2).PSTHz(Sel&NAneurons,Ushow),1); %calculate standard error of the mean
+up4=psth4+sem4;
+down4=psth4-sem4;
+plot(time2,psth3,'Color',sucrose,'linewidth',1);
+plot(time2,psth4,'Color',maltodextrin,'linewidth',1);
+patch([time2,time2(end:-1:1)],[up3,down3(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
+patch([time2,time2(end:-1:1)],[up4,down4(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
+
+plot([-2 5],[0 0],':','color','k','linewidth',0.75);
+plot([0 0],[-2 8],':','color','k','linewidth',0.75);
+plot([PEWindow(1) PEWindow(1)],[-2 8],'color','b','linewidth',0.85);
+plot([PEWindow(2) PEWindow(2)],[-2 8],'color','b','linewidth',0.85);
+axis([XaxisPE(1) XaxisPE(2) min(down3)-0.2 max(up4)+0.2]);
+ylabel('Mean firing (z-score)');
+xlabel('Seconds from PE');
+legend('Post-suc','Post-mal','location','northeast');
+
+if cell2mat(R_2R.RewHist.PrevRewStatsPENASubj{1,1}(2,7))<0.05
+    plot(PEWindow(1)+(PEWindow(2)-PEWindow(1))/2,max(up3)+0.1,'*','color','k','markersize',13);
+end
+
+
+%VP
+
+%plot suc after suc
+psth1=nanmean(R_2R.Ev(RD1P1).PSTHz(Sel&VPneurons,Ishow),1); 
+sem1=nanste(R_2R.Ev(RD1P1).PSTHz(Sel&VPneurons,Ishow),1); %calculate standard error of the mean
+up1=psth1+sem1;
+down1=psth1-sem1;
+
+%plot suc after malt
+psth2=nanmean(R_2R.Ev(RD1P2).PSTHz(Sel&VPneurons,Ishow),1); 
+sem2=nanste(R_2R.Ev(RD1P2).PSTHz(Sel&VPneurons,Ishow),1); %calculate standard error of the mean
+up2=psth2+sem2;
+down2=psth2-sem2;
+
+%make the plot
+subplot(2,4,2);
+title(['Cue response, VP'])
+hold on;
+plot(time1,psth1,'Color',sucrose,'linewidth',1);
+plot(time1,psth2,'Color',maltodextrin,'linewidth',1);
+
+patch([time1,time1(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
+patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
+
+plot([-2 5],[0 0],':','color','k','linewidth',0.75);
+plot([0 0],[-2 8],':','color','k','linewidth',0.75);
+plot([CueWindow(1) CueWindow(1)],[-2 8],'color','b','linewidth',0.85);
+plot([CueWindow(2) CueWindow(2)],[-2 8],'color','b','linewidth',0.85);
+axis([XaxisCue(1) XaxisCue(2) min(down2)-0.3 max(up1)+0.3]);
+ylabel('Mean firing (z-score)');
+xlabel('Seconds from cue');
+legend('Post-suc','Post-mal','location','northeast');
+
+if cell2mat(R_2R.RewHist.PrevRewStatsCueVPSubj{1,1}(2,7))<0.05
+    plot(CueWindow(1)+(CueWindow(2)-CueWindow(1))/2,max(up2)+0.1,'*','color','k','markersize',13);
+end
+
+subplot(2,4,6);
+title('PE response, VP');
+hold on;
+%plot malt after suc
+psth3=nanmean(R_2R.Ev(RD2P1).PSTHz(Sel&VPneurons,Ushow),1); 
+sem3=nanste(R_2R.Ev(RD2P1).PSTHz(Sel&VPneurons,Ushow),1); %calculate standard error of the mean
+up3=psth3+sem3;
+down3=psth3-sem3;
+
+%plot malt after malt
+psth4=nanmean(R_2R.Ev(RD2P2).PSTHz(Sel&VPneurons,Ushow),1); 
+sem4=nanste(R_2R.Ev(RD2P2).PSTHz(Sel&VPneurons,Ushow),1); %calculate standard error of the mean
+up4=psth4+sem4;
+down4=psth4-sem4;
+
+plot(time2,psth3,'Color',sucrose,'linewidth',1);
+plot(time2,psth4,'Color',maltodextrin,'linewidth',1);
+patch([time2,time2(end:-1:1)],[up3,down3(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
+patch([time2,time2(end:-1:1)],[up4,down4(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
+
+plot([-2 5],[0 0],':','color','k','linewidth',0.75);
+plot([0 0],[-2 8],':','color','k','linewidth',0.75);
+plot([PEWindow(1) PEWindow(1)],[-2 8],'color','b','linewidth',0.85);
+plot([PEWindow(2) PEWindow(2)],[-2 8],'color','b','linewidth',0.85);
+axis([XaxisPE(1) XaxisPE(2) min(down4)-0.3 max(up3)+0.3]);
+ylabel('Mean firing (z-score)');
+xlabel('Seconds from PE');
+legend('Post-suc','Post-mal','location','northeast');
+
+if cell2mat(R_2R.RewHist.PrevRewStatsPEVPSubj{1,1}(2,7))<0.05
+    plot(PEWindow(1)+(PEWindow(2)-PEWindow(1))/2,max(up3)+0.1,'*','color','k','markersize',13);
+end
+
+
+%% plot linear model coefficients
+
+%Plot all neurons
+Sel=NAneurons<2;
+
+%coefficients for each trial
+subplot(2,4,3);
+hold on;
+errorbar(1:trialsback,nanmean(R_2R.RewHist.LinMdlCueWeights(Sel&NAneurons,1:trialsback),1),nanste(R_2R.RewHist.LinMdlCueWeights(Sel&NAneurons,1:trialsback),1),'color',NAc);
+errorbar(1:trialsback,nanmean(R_2R.RewHist.LinMdlCueWeights(Sel&VPneurons,1:trialsback),1),nanste(R_2R.RewHist.LinMdlCueWeights(Sel&VPneurons,1:trialsback),1),'color',VP);
+errorbar(1:trialsback,nanmean(R_2R.RewHist.LinMdlCueWeightsSh(Sel&NAneurons,1:trialsback),1),nanste(R_2R.RewHist.LinMdlCueWeights(Sel&NAneurons,1:trialsback),1),'color','k');
+errorbar(1:trialsback,nanmean(R_2R.RewHist.LinMdlCueWeightsSh(Sel&VPneurons,1:trialsback),1),nanste(R_2R.RewHist.LinMdlCueWeights(Sel&VPneurons,1:trialsback),1),'color','k');
+xlabel('Trials back');
+ylabel('Mean coefficient weight');
+title('Linear model coefficients');
+axis([0 trialsback+1 -0.5 1]);
+plot([-1 trialsback+1],[0 0],':','color','k','linewidth',0.75);
+legend('NAc','VP','Shuff');
+
+%stats to check if VP and NAc are greater than chance
+R_2R.RewHist.LinCoeffCueMultComp=[];
+[c,~,~,~]=multcompare(R_2R.RewHist.LinCoeffStatsCueNASubj{1,2},'dimension',[1,2],'display','off');
+[d,~,~,~]=multcompare(R_2R.RewHist.LinCoeffStatsCueVPSubj{1,2},'dimension',[1,2],'display','off');
+for i=1:trialsback
+
+    %NAc vs shuff
+    Cel=c(:,1)==2*(i-1)+1 & c(:,2)==2*(i-1)+2;
+    if c(Cel,6)<0.05 R_2R.RewHist.LinCoeffCueMultComp(i,1)=1; else R_2R.RewHist.LinCoeffCueMultComp(i,1)=0; end
+    R_2R.RewHist.LinCoeffCueMultComp(i,2)=c(Cel,2);
+
+    %VP vs shuff
+    Cel=d(:,1)==2*(i-1)+1 & d(:,2)==2*(i-1)+2;
+    if d(Cel,6)<0.05 R_2R.RewHist.LinCoeffCueMultComp(i,3)=1; else R_2R.RewHist.LinCoeffCueMultComp(i,3)=0; end
+    R_2R.RewHist.LinCoeffCueMultComp(i,4)=d(Cel,4);
+end
+plot([1:trialsback]-0.1,(R_2R.RewHist.LinCoeffCueMultComp(:,1)-0.5)*1.3,'*','color',NAc); %VP vs shuff
+plot([1:trialsback]+0.1,(R_2R.RewHist.LinCoeffCueMultComp(:,3)-0.5)*1.3,'*','color',VP); %NAc vs shuff
+
+
+%number of neurons with significant weights
+subplot(2,4,4);
+hold on;
+plot(1:trialsback,sum(R_2R.RewHist.LinMdlCuePVal(Sel&NAneurons,1:trialsback)<0.05,1)/sum(Sel&NAneurons),'color',NAc);
+plot(1:trialsback,sum(R_2R.RewHist.LinMdlCuePVal(Sel&VPneurons,1:trialsback)<0.05,1)/sum(Sel&VPneurons),'color',VP);
+plot(1:trialsback,sum(R_2R.RewHist.LinMdlCuePValSh(Sel&NAneurons,1:trialsback)<0.05,1)/sum(Sel&NAneurons),'color',NAc/3);
+plot(1:trialsback,sum(R_2R.RewHist.LinMdlCuePValSh(Sel&VPneurons,1:trialsback)<0.05,1)/sum(Sel&VPneurons),'color',VP/3);
+axis([0 trialsback+1 0 0.5]);
+xlabel('Trials back');
+ylabel('Fraction of the population');
+title('Outcome-modulated cue resp');
+legend('NAc','VP','Shuff');
+
+%Chi squared stat for each trial
+for i=1:trialsback
+    [~,R_2R.RewHist.LinMdlCue2all(i,1),R_2R.RewHist.LinMdlCue2all(i,2)]=crosstab(cat(1,R_2R.RewHist.LinMdlCuePVal(Sel,i)<0.05,R_2R.RewHist.LinMdlCuePValSh(Sel,i)<0.05),cat(1,VPneurons,VPneurons+2));
+    [~,R_2R.RewHist.LinMdlCue2region(i,1),R_2R.RewHist.LinMdlCue2region(i,2)]=crosstab(R_2R.RewHist.LinMdlCuePVal(Sel,i)<0.05,VPneurons);
+end
+%plot([1:trialsback]-0.1,(R_2R.RewHist.LinMdlCue2all(:,2)<0.05)-0.52,'*','color',NAc);
+plot([1:trialsback],(R_2R.RewHist.LinMdlCue2region(:,2)<0.05&R_2R.RewHist.LinMdlCue2all(:,2)<0.05)-0.52,'*','color','k');
+
+%% PE
+
+%coefficients for each trial
+subplot(2,4,7);
+hold on;
+errorbar(1:trialsback,nanmean(R_2R.RewHist.LinMdlPEWeights(Sel&NAneurons,1:trialsback),1),nanste(R_2R.RewHist.LinMdlPEWeights(Sel&NAneurons,1:trialsback),1),'color',NAc);
+errorbar(1:trialsback,nanmean(R_2R.RewHist.LinMdlPEWeights(Sel&VPneurons,1:trialsback),1),nanste(R_2R.RewHist.LinMdlPEWeights(Sel&VPneurons,1:trialsback),1),'color',VP);
+errorbar(1:trialsback,nanmean(R_2R.RewHist.LinMdlPEWeightsSh(Sel&NAneurons,1:trialsback),1),nanste(R_2R.RewHist.LinMdlPEWeights(Sel&NAneurons,1:trialsback),1),'color','k');
+errorbar(1:trialsback,nanmean(R_2R.RewHist.LinMdlPEWeightsSh(Sel&VPneurons,1:trialsback),1),nanste(R_2R.RewHist.LinMdlPEWeights(Sel&VPneurons,1:trialsback),1),'color','k');
+xlabel('Trials back');
+ylabel('Mean coefficient weight');
+title('Linear model coefficients');
+axis([0 trialsback+1 -0.5 1]);
+plot([-1 trialsback+1],[0 0],':','color','k','linewidth',0.75);
+legend('NAc','VP','Shuff');
+
+
+%stats to check if VP and NAc are greater than chance
+R_2R.RewHist.LinCoeffPEMultComp=[];
+[c,~,~,~]=multcompare(R_2R.RewHist.LinCoeffStatsPENASubj{1,2},'dimension',[1,2],'display','off');
+[d,~,~,~]=multcompare(R_2R.RewHist.LinCoeffStatsPEVPSubj{1,2},'dimension',[1,2],'display','off');
+for i=1:trialsback
+
+    %NAc vs shuff
+    Cel=c(:,1)==2*(i-1)+1 & c(:,2)==2*(i-1)+2;
+    if c(Cel,6)<0.05 R_2R.RewHist.LinCoeffPEMultComp(i,1)=1; else R_2R.RewHist.LinCoeffPEMultComp(i,1)=0; end
+    R_2R.RewHist.LinCoeffPEMultComp(i,2)=c(Cel,2);
+
+    %VP vs shuff
+    Cel=d(:,1)==2*(i-1)+1 & d(:,2)==2*(i-1)+2;
+    if d(Cel,6)<0.05 R_2R.RewHist.LinCoeffPEMultComp(i,3)=1; else R_2R.RewHist.LinCoeffPEMultComp(i,3)=0; end
+    R_2R.RewHist.LinCoeffPEMultComp(i,4)=d(Cel,4);
+end
+plot([1:trialsback]-0.1,(R_2R.RewHist.LinCoeffPEMultComp(:,1)-0.5)*1.3,'*','color',NAc); %VP vs shuff
+plot([1:trialsback]+0.1,(R_2R.RewHist.LinCoeffPEMultComp(:,3)-0.5)*1.3,'*','color',VP); %NAc vs shuff
+
+
+%number of neurons with significant weights
+subplot(2,4,8);
+hold on;
+plot(1:trialsback,sum(R_2R.RewHist.LinMdlPEPVal(Sel&NAneurons,1:trialsback)<0.05,1)/sum(Sel&NAneurons),'color',NAc);
+plot(1:trialsback,sum(R_2R.RewHist.LinMdlPEPVal(Sel&VPneurons,1:trialsback)<0.05,1)/sum(Sel&VPneurons),'color',VP);
+plot(1:trialsback,sum(R_2R.RewHist.LinMdlPEPValSh(Sel&NAneurons,1:trialsback)<0.05,1)/sum(Sel&NAneurons),'color',NAc/3);
+plot(1:trialsback,sum(R_2R.RewHist.LinMdlPEPValSh(Sel&VPneurons,1:trialsback)<0.05,1)/sum(Sel&VPneurons),'color',VP/3);
+axis([0 trialsback+1 0 0.5]);
+xlabel('Trials back');
+ylabel('Fraction of the population');
+title('Outcome-modulated PE resp');
+legend('NAc','VP','Shuff');
+
+
+%Chi squared stat for each trial
+for i=1:trialsback
+    [~,R_2R.RewHist.LinMdlPEX2all(i,1),R_2R.RewHist.LinMdlPEX2all(i,2)]=crosstab(cat(1,R_2R.RewHist.LinMdlPEPVal(Sel,i)<0.05,R_2R.RewHist.LinMdlPEPValSh(Sel,i)<0.05),cat(1,VPneurons,VPneurons+2));
+    [~,R_2R.RewHist.LinMdlPEX2region(i,1),R_2R.RewHist.LinMdlPEX2region(i,2)]=crosstab(R_2R.RewHist.LinMdlPEPVal(Sel,i)<0.05,VPneurons);
+end
+%plot([0:trialsback]-0.1,(R_2R.RewHist.LinMdlPE2all(:,2)<0.05)-0.52,'*','color',NAc);
+plot([1:trialsback],(R_2R.RewHist.LinMdlPEX2region(:,2)<0.05&R_2R.RewHist.LinMdlPEX2all(:,2)<0.05)-0.52,'*','color','k');
+
+%% stats comparing PE coefficient weight in first 2 trials in selective and non-selective neurons in VP
+Sel=R_2R.SucN | R_2R.MalN;
+NSel=(R_2R.SucN | R_2R.MalN) == 0;
+[~,R_2R.RewHist.SelectiveHistory{1,1},R_2R.RewHist.SelectiveHistory{1,2}]=anovan(cat(1,R_2R.RewHist.LinMdlPEWeights(Sel&VPneurons,1),R_2R.RewHist.LinMdlPEWeights(NSel&VPneurons,1),R_2R.RewHist.LinMdlPEWeights(Sel&VPneurons,2),R_2R.RewHist.LinMdlPEWeights(NSel&VPneurons,2)),...
+{cat(1,ones(sum(Sel&VPneurons),1),zeros(sum(NSel&VPneurons),1),ones(sum(Sel&VPneurons),1),zeros(sum(NSel&VPneurons),1)),cat(1,ones(sum(VPneurons),1),zeros(sum(VPneurons),1)),...
+cat(1,R_2R.Ninfo(Sel&VPneurons,4),R_2R.Ninfo(NSel&VPneurons,4),R_2R.Ninfo(Sel&VPneurons,4),R_2R.Ninfo(NSel&VPneurons,4))},'varnames',{'Selective','Trial','Subject'},'random',3,'model','full','display','off');
+
+save('R_2R.mat','R_2R');

+ 454 - 0
MatlabScripts/j_Fig5Water.m

@@ -0,0 +1,454 @@
+clear all;
+load ('R_W.mat');
+load ('RAW.mat');
+
+%get parameters
+BinBase=R_W.Param.BinBase;
+BinDura=R_W.Param.BinDura;
+bins=R_W.Param.bins;
+binint=R_W.Param.binint;
+binstart=R_W.Param.binstart;
+PStatBins=0.01; %using more stringent cutoff to reduce pre-delivery noise
+
+%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; %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_W.Ninfo) %runs through neurons
+        if R_W.BinStatRwrd{i,1}(k).IntSig < PStatBins %if neuron is significant for this bin
+            if R_W.BinStatRwrd{i-1,1}(k).IntSig < PStatBins || R_W.BinStatRwrd{i+1,1}(k).IntSig < PStatBins %either previous or subsequent bin must be significant
+                if R_W.BinStatRwrd{i,1}(k).R1Mean > R_W.BinStatRwrd{i,1}(k).R2Mean %if firing is greater on sucrose trials than maltodextrin
+                    R_W.BinRewPref{i,1}(k,1)=1;
+                else
+                    R_W.BinRewPref{i,1}(k,1)=2; %otherwise maltodextrin must be greater than sucrose
+                end
+            else
+                R_W.BinRewPref{i,1}(k,1)=0; %if not significant in 2 consecutive bins
+            end
+        else
+            R_W.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_W.BinRewPref{i,1}(:,1)==1); %sucrose pref
+    NN2perBin(i,1)=sum(R_W.BinRewPref{i,1}(:,1)==2); %malto pref
+    NNperBin(i,1)=sum(R_W.BinRewPref{i,1}(:,1)>0); %either
+    %normalize to number of neurons in population
+    NN1norm=NN1perBin./sum(length(R_W.Ninfo));
+    NN2norm=NN2perBin./sum(length(R_W.Ninfo));
+    NNnorm=NNperBin./sum(length(R_W.Ninfo));
+   
+
+end
+
+%Cumulative reward selectivity
+cumsel=zeros(length(R_W.Ninfo),bins); %set up result matrix
+cumsel1=zeros(length(R_W.Ninfo),bins); %set up result matrix for sucrose
+cumsel2=zeros(length(R_W.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_W.Ninfo) %go through each neuron
+    for j=2:bins-1 %look in each bin we did analysis on
+        if R_W.BinRewPref{j,1}(i,1)>0 || cumsel(i,j-1)==1
+            cumsel(i,j) = 1;
+        end
+        if R_W.BinRewPref{j,1}(i,1)==1 || cumsel1(i,j-1)==1
+            cumsel1(i,j) = 1;
+        end
+        if R_W.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
+subplot(3,2,2);
+hold on;
+plot(xaxis,NNnorm(2:bins-1),'Color', total,'linewidth',1.5);
+plot(xaxis,NN1norm(2:bins-1),'Color',water,'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.7]);
+legend('Total','Wat > mal','Mal > wat','Location','northeast')
+ylabel('Fraction of population');
+title('Reward-selective neurons');
+xlabel('Seconds from RD');
+
+
+
+%% plotting maltodextrin-selective neurons
+
+%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_W.Erefnames);
+ev2=strcmp('RD2', R_W.Erefnames);
+
+%setting up parameters
+Xaxis=[-2 5];
+inttime=find(R_W.Param.Tm>=Xaxis(1) & R_W.Param.Tm<=Xaxis(2));
+paramtime=R_W.Param.Tm(inttime);
+
+%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
+    Pref1(:,i)=R_W.BinRewPref{Bin1+i}==2; %get neurons that have greater firing for water in any of the bins bounded above
+    Resp11(:,i)=Pref1(:,i)&cat(1,R_W.BinStatRwrd{Bin1+i,1}.WatRespDir)==-1; %get neurons with inhibition to water
+    Resp12(:,i)=Pref1(:,i)&cat(1,R_W.BinStatRwrd{Bin1+i,1}.MalRespDir)==1;%get neurons with excitation to maltodextrin
+end
+
+Mel=sum(Pref1,2)>0;  %all neurons selective for mal in any bin
+Mel2=sum(Resp11,2)>0; %all selective neurons water inhibited in any bin
+Mel3=sum(Resp12,2)>0; %all selective neurons malto excited in any bin
+
+subplot(9,40,[(15:23)+120 (15:23)+160]); %heatmap of maltodextrin preferring neurons' response to water
+Neurons=R_W.Ev(ev1).PSTHz(Mel,inttime); %get the firing rates of neurons of interest
+MalResp=cat(1,R_W.BinStatRwrd{SortBin+1,1}.R2Mean); %water responses
+TMP=MalResp(Mel); %find the magnitude of the excitations for water 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,1)],Neurons,ClimE);
+title(['Water trials']);
+ylabel('Individual neurons');
+hold on;
+plot([0 0],[0 sum(Mel)],':','color','k','linewidth',0.75);
+
+subplot(9,40,[(27:35)+120 (27:35)+160]); %heatmap of maltodextrin preferring neurons' response to maltodextrin
+Neurons=R_W.Ev(ev2).PSTHz(Mel,inttime); %get the firing rates of neurons of interest
+Neurons=Neurons(SORTimg,:); %sort the neurons same as before
+imagesc(paramtime,[1,sum(Mel,1)],Neurons,ClimE);
+title(['Maltodextrin trials']);
+xlabel('Seconds from RD');
+hold on;
+plot([0 0],[0 sum(Mel)],':','color','k','linewidth',0.75);
+
+%plot water preferring neurons to water
+psth1=nanmean(R_W.Ev(ev1).PSTHz(Mel,inttime),1); 
+sem1=nanste(R_W.Ev(ev1).PSTHz(Mel,inttime),1); %calculate standard error of the mean
+up1=psth1+sem1;
+down1=psth1-sem1;
+
+%plot water preferring neurons to malt
+psth2=nanmean(R_W.Ev(ev2).PSTHz(Mel,inttime),1); 
+sem2=nanste(R_W.Ev(ev2).PSTHz(Mel,inttime),1); %calculate standard error of the mean
+up2=psth2+sem2;
+down2=psth2-sem2;
+
+%plotting
+subplot(9,3,[10 13]);
+title(['Mal>wat (n=' num2str(sum(Mel)) ' of ' num2str(length(Mel)) ')'])
+hold on;
+plot(paramtime,psth1,'Color',water,'linewidth',1);
+plot(paramtime,psth2,'Color',maltodextrin,'linewidth',1);
+patch([paramtime,paramtime(end:-1:1)],[up1,down1(end:-1:1)],water,'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 10],':','color','k','linewidth',0.75);
+axis([-2 5 -1.5 3.5]);
+ylabel('Mean firing (z-score)');
+legend('Wat trials','Mal trials');
+xlabel('Seconds from RD');
+
+
+%display inhibitions and excitations
+both = sum(Mel2&Mel3);
+excite = sum(Mel3)-both;
+inhib = sum(Mel2)-both;
+
+subplot(9,40,[160 200]);
+b = bar([inhib both excite; 0 0 0],'stacked');
+b(1,1).FaceColor=water;
+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',[])
+
+
+%% emergence of maltodextrin and water responses
+
+%select which time point we're examining
+Window=[0.8 1.8]; %For looking at emergence of excitation
+NN=0; %neuron counter
+Sess=0; %session counter
+
+%plotting
+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];
+
+%get firing rates on each trial
+for i=1:length(RAW) %loops through sessions
+    if strcmp('WA',RAW(i).Type(1:2))
+        Sess=Sess+1; %session counter
+        SN{Sess}=0; %selective neuron counter
+        ev1=strmatch('RD1',RAW(i).Einfo(:,2),'exact');
+        ev2=strmatch('RD2',RAW(i).Einfo(:,2),'exact');
+        for j= 1:size(RAW(i).Nrast,1) %Number of neurons within sessions
+            NN=NN+1;
+            if Mel(NN)==1
+                
+                SN{Sess}=SN{Sess}+1;
+                Bcell1=0;
+                Bcell2=0;
+                Bcell=0;
+                ev1spk=0;
+                Ev2spk=0;
+
+                %get mean baseline firing for all reward trials
+                [Bcell1,B1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{ev1},BinBase,{2});% makes trial by trial rasters for baseline
+                for y= 1:B1n
+                    Bcell(1,y)=sum(Bcell1{1,y}>BinBase(1));
+                end
+                [Bcell2,B2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{ev2},BinBase,{2});% makes trial by trial rasters for baseline
+                for z= 1:B2n
+                    Bcell(1,z+B1n)=sum(Bcell2{1,z}>BinBase(1));
+                end
+                Bhz=Bcell/(BinBase(1,2)-BinBase(1,1));
+                Bmean=nanmean(Bhz);
+                Bstd=nanstd(Bhz);
+
+                %get trial by trial firing rate for water trials
+                [ev1cell,ev1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{ev1},Window,{2});% makes trial by trial rasters for event
+                for y= 1:ev1n
+                    ev1spk(SN{Sess},y)=sum(ev1cell{1,y}>Window(1));
+                end       
+                ev1hz=ev1spk(SN{Sess},:)/(Window(1,2)-Window(1,1));
+                for x= 1:ev1n
+                    ev1norm{Sess}(SN{Sess},x)=((ev1hz(1,x)-Bmean)/Bstd);
+                end
+
+                %get trial by trial firing rate for maltodextrin trials
+                [ev2cell,ev2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{ev2},Window,{2});% makes trial by trial rasters for event
+                for y= 1:ev2n
+                    ev2spk(1,y)=sum(ev2cell{1,y}>Window(1));
+                end       
+                ev2hz=ev2spk/(Window(1,2)-Window(1,1));
+                for x= 1:ev2n
+                    ev2norm{Sess}(SN{Sess},x)=((ev2hz(1,x)-Bmean)/Bstd);
+                end
+            end
+        end %neurons in each session
+        
+        trialmeanz1{Sess}=NaN(2,1);
+        trialmeanz2{Sess}=NaN(2,1);
+        ind=0;
+        order=0;
+
+        [~,ind]=sort(cat(1,RAW(i).Erast{ev1},RAW(i).Erast{ev2}));
+        for j=1:length(ind)
+            order(j,1)=find(ind==j);
+        end
+        xaxis1{Sess}=order(1:length(ev1norm{Sess}(1,:)));
+        xaxis2{Sess}=order(1+length(ev1norm{Sess}(1,:)):end);
+
+        for k=1:length(ev1norm{Sess}(1,:))
+            trialmeanz1{Sess}(1,k)=nanmean(ev1norm{Sess}(:,k));
+            trialmeanz1{Sess}(2,k)=nanste(ev1norm{Sess}(:,k),1);
+        end
+
+        for k=1:length(ev2norm{Sess}(1,:))
+            trialmeanz2{Sess}(1,k)=nanmean(ev2norm{Sess}(:,k));
+            trialmeanz2{Sess}(2,k)=nanste(ev2norm{Sess}(:,k),1);
+        end
+        
+        %stats checking for interaction between # of reward presentations
+        %and reward on firing rate
+        numtrials=min([length(xaxis1{Sess}) length(xaxis2{Sess})]); %number of trials for stats (fewest number of trials for the 2 outcomes)
+        trials1=repmat(1:numtrials,length(ev1norm{Sess}(:,1)),1);
+        trials2=repmat(1:numtrials,length(ev2norm{Sess}(:,1)),1);
+        reward1=zeros(size(trials1));
+        reward2=ones(size(trials2));
+        data1=ev1norm{Sess}(:,1:numtrials);
+        data2=ev2norm{Sess}(:,1:numtrials);
+
+        [~,R_W.EmergenceStats{Sess,1},R_W.EmergenceStats{Sess,2}]=anovan(cat(1,data1(:),data2(:)),{cat(1,trials1(:),trials2(:)),cat(1,reward1(:),reward2(:))},'varnames',{'trials','reward'},'display','off','model','full');
+
+     end %checking if a water session
+     
+     
+    
+end %session
+
+%plotting
+for i=1:Sess
+    subplot(3,2,4+i);
+    hold on;
+    errorbar(xaxis1{i},trialmeanz1{i}(1,:),trialmeanz1{i}(2,:),'*','Color',water,'linewidth',1.5);
+    errorbar(xaxis2{i},trialmeanz2{i}(1,:),trialmeanz2{i}(2,:),'*','Color',maltodextrin,'linewidth',1.5);
+    %plot(xaxis1,(cell2mat(ev1stat{i}(1,:))*13-16.9),'*','Color',water); %(12.5-trialmeanz242(1,1:trials1)-trialmeanz242(2,1:trials1))
+    %plot(xaxis2,(cell2mat(ev2stat{i}(1,:))*13-16.7),'*','Color',maltodextrin); %(12.5-trialmeanz242(1,1:trials1)-trialmeanz242(2,1:trials1))
+    %plot(xaxis1,(cell2mat(trialstat42(1,1:trials1))*13-17.1),'*','Color','k'); %(12.5-trialmeanz242(1,1:trials1)-trialmeanz242(2,1:trials1))
+    plot([0 60],[0 0],':','color','k','linewidth',0.75);
+    ylabel(['Mean firing ' num2str(Window(1,1)) '-' num2str(Window(1,2)) 's (z-score)']);
+    title(['Rat ' num2str(i) ' (n=' num2str(SN{i}) ')']);
+    axis([1 max([xaxis1{i};xaxis2{i}]) min(trialmeanz1{i}(1,:))*1.3 max(trialmeanz2{i}(1,:))*1.3]);
+    xlabel('Trial #');
+    legend('Wat','Mal','location','northwest');
+
+end
+
+
+%% conduct lick analysis
+global Dura Tm BSIZE Tbin
+path='C:\Users\dottenh2\Documents\MATLAB\David\2Rewards Nex Files\2RLick_paper.xls';
+
+%Main settings
+BSIZE=0.01; %Do not change
+Dura=[-22 20]; Tm=Dura(1):BSIZE:Dura(2);
+Tbin=-0.5:0.005:0.5; %window used to determine the optimal binsize
+MinNumTrials=5;
+Lick=[];Lick.Ninfo={};LL=0;Nlick=0;
+
+%Smoothing
+Smoothing=1; %0 for raw and 1 for smoothing
+SmoothTYPE='lowess'; %can change this between lowess and rlowess (more robust, ignores outliers more)
+SmoothSPAN=50; %percentage of total data points
+if Smoothing~=1, SmoothTYPE='NoSmoothing';SmoothSPAN=NaN; end
+
+% List of events to analyze and analysis windows EXTRACTED from excel file
+[~,Erefnames]=xlsread(path,'Windows','a3:a10'); % cell that contains the event names
+
+%Finds the total number of sessions
+for i=1:length(RAW)
+    if strcmp('WA',RAW(i).Type(1:2))
+        Nlick=Nlick+1;
+        Lick.Linfo(i,1)=RAW(i).Ninfo(1,1);
+    end
+end
+Lick.Erefnames= Erefnames;
+
+%preallocating the result matrix
+for k=1:length(Erefnames)
+    Lick.Ev(k).PSTHraw(1:Nlick,1:length(Tm))=NaN(Nlick,length(Tm));
+    Lick.Ev(k).BW(1:Nlick,1)=NaN;
+    Lick.Ev(k).NumberTrials(1:Nlick,1)=NaN;
+end
+
+for i=1:length(RAW) %loops through sessions
+    if strcmp('WA',RAW(i).Type(1:2))
+        LL=LL+1; %lick session counter
+        for k=1:length(Erefnames) %loops thorough the events
+            EvInd=strcmp(Erefnames(k),RAW(i).Einfo(:,2)); %find the event id number from RAW
+            LickInd=strcmp('Licks',RAW(i).Einfo(:,2)); %find the event id number from RAW
+            if sum(EvInd)==0
+                fprintf('HOWDY, CANT FIND EVENTS FOR ''%s''\n',Erefnames{k});
+            end
+
+            Lick.Ev(k).NumberTrials(LL,1)=length(RAW(i).Erast{EvInd});
+            if  ~isempty(EvInd) && Lick.Ev(k).NumberTrials(LL,1)>MinNumTrials %avoid analyzing sessions where that do not have enough trials
+                [PSR1,N1]=MakePSR04(RAW(i).Erast(LickInd),RAW(i).Erast{EvInd},Dura,{1});% makes collpased rasters. PSR1 is a cell(neurons)
+                if ~isempty(PSR1{1}) %to avoid errors, added on 12/28 2011
+                    %forcing 100ms bin size to keep it consistent across
+                    %sessions (no reason is should be different for licks)
+                    [PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Erast{LickInd},1),1),{2,0,0.1});%these values force bin size to be 100ms
+                    PTH1=smooth(PTH1,SmoothSPAN,SmoothTYPE)';
+
+                    %------------- Fills the R.Ev(k) fields --------------
+                    Lick.Ev(k).BW(LL,1)=BW1;
+                    Lick.Ev(k).PSTHraw(LL,1:length(Tm))=PTH1;
+                end
+            end
+        end
+    end
+end
+
+Xaxis=[-2 12];
+Ishow=find(Tm>=Xaxis(1) & Tm<=Xaxis(2));
+time1=Tm(Ishow);
+c=[-1000 7500];ClimE=sign(c).*abs(c).^(1/4);%ColorMapExc
+colormap(jet);
+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];
+
+Ev1=strcmp('RD1', Lick.Erefnames);
+Ev2=strcmp('RD2', Lick.Erefnames);
+
+psth1=nanmean(Lick.Ev(Ev1).PSTHraw(:,Ishow),1);
+sem1=nanste(Lick.Ev(Ev1).PSTHraw(:,Ishow),1); %calculate standard error of the mean
+up1=psth1+sem1;
+down1=psth1-sem1;
+
+psthE=nanmean(Lick.Ev(Ev2).PSTHraw(:,Ishow),1);
+semE=nanste(Lick.Ev(Ev2).PSTHraw(:,Ishow),1); %calculate standard error of the mean
+upE=psthE+semE;
+downE=psthE-semE;
+
+%plotting
+subplot(3,2,1);
+hold on;
+plot(time1,psth1,'Color',water,'linewidth',1);
+plot(time1,psthE,'Color',maltodextrin,'linewidth',1);
+patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],water,'EdgeColor','none');alpha(0.5);
+patch([time1,time1(end:-1:1)],[upE,downE(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
+title('Mean lick rate');
+%plot([-2 10],[0 0],':','color','k');
+plot([0 0],[-2 8],':','color','k','linewidth',0.75);
+axis([-2 12 0 7.5]);
+xlabel('Seconds from reward delivery');
+ylabel('Licks/s');
+legend('Water trials','Maltodextrin trials');
+%color map
+[magma,inferno,plasma,viridis]=colormaps;
+colormap(plasma);
+c=[-100 2000];ClimE=sign(c).*abs(c).^(1/4);%colormap
+
+%% creating and saving new variables
+R_W.MalN=Mel;
+
+% proportion summary and chi squared test for reward selective neurons in
+% sucrose and water sessions in VP
+
+%find all reward-selective neurons
+for i = 1:(Bin2-Bin1+1)
+    %the added +1 is necessary because bin 20 is the 21st entry in the matrix
+    Selective(:,i)=R_W.BinRewPref{Bin1+i}>0; %get neurons that have greater firing for water in any of the bins bounded above
+end
+WSel=sum(Pref1,2)>0;  %all neurons selective for mal in any bin
+
+%get reward selective neurons from these two rats
+load ('R_2R.mat');
+selection=(R_2R.SucN | R_2R.MalN); %selective neurons in suc vs mal sessions
+SSel=selection(strcmp('VP2',R_2R.Ninfo(:,4)) | strcmp('VP5',R_2R.Ninfo(:,4)),1);
+
+%get proportion and stats
+R_W.SucvsWatX2(1,1)=sum(SSel)/length(SSel); %selective neurons in maltodextrin sessions
+R_W.SucvsWatX2(1,2)=sum(WSel)/length(WSel); %how many selective neurons in water sessions
+[~,R_W.SucvsWatX2(2,1),R_W.SucvsWatX2(2,2)]=crosstab(cat(1,SSel,WSel),cat(1,ones(length(SSel),1),zeros(length(WSel),1))); %find chi squared value for this distribution
+
+save('R_W.mat','R_W');

+ 339 - 0
MatlabScripts/k_Fig6ThreeRewards.m

@@ -0,0 +1,339 @@
+%plotting 3rewards data
+
+clear all;
+load ('R_3R.mat');
+load ('RAW.mat');
+
+%get parameters
+BinBase=R_3R.Param.BinBase;
+BinDura=R_3R.Param.BinDura;
+bins=R_3R.Param.bins;
+binint=R_3R.Param.binint;
+binstart=R_3R.Param.binstart;
+PStatBins=0.01; %using more stringent cutoff to reduce pre-delivery noise
+
+%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];
+
+%% conduct lick analysis
+global Dura Tm BSIZE Tbin
+path='C:\Users\dottenh2\Documents\MATLAB\David\2Rewards Nex Files\3RLick_paper.xls';
+
+%Main settings
+BSIZE=0.01; %Do not change
+Dura=[-22 20]; Tm=Dura(1):BSIZE:Dura(2);
+Tbin=-0.5:0.005:0.5; %window used to determine the optimal binsize
+MinNumTrials=5;
+Lick=[];Lick.Ninfo={};LL=0;Nlick=0;
+
+%Smoothing
+Smoothing=1; %0 for raw and 1 for smoothing
+SmoothTYPE='lowess'; %can change this between lowess and rlowess (more robust, ignores outliers more)
+SmoothSPAN=50; %percentage of total data points
+if Smoothing~=1, SmoothTYPE='NoSmoothing';SmoothSPAN=NaN; end
+
+% List of events to analyze and analysis windows EXTRACTED from excel file
+[~,Erefnames]=xlsread(path,'Windows','a3:a8'); % cell that contains the event names
+
+%Finds the total number of sessions
+for i=1:length(RAW)
+    if strcmp('TH',RAW(i).Type(1:2))
+        Nlick=Nlick+1;
+        Lick.Linfo(i,1)=RAW(i).Ninfo(1,1);
+    end
+end
+Lick.Erefnames= Erefnames;
+
+%preallocating the result matrix
+for k=1:length(Erefnames)
+    Lick.Ev(k).PSTHraw(1:Nlick,1:length(Tm))=NaN(Nlick,length(Tm));
+    Lick.Ev(k).BW(1:Nlick,1)=NaN;
+    Lick.Ev(k).NumberTrials(1:Nlick,1)=NaN;
+end
+
+for i=1:length(RAW) %loops through sessions
+    if strcmp('TH',RAW(i).Type(1:2))
+        LL=LL+1; %lick session counter
+        for k=1:length(Erefnames) %loops thorough the events
+            EvInd=strcmp(Erefnames(k),RAW(i).Einfo(:,2)); %find the event id number from RAW
+            LickInd=strcmp('Licks',RAW(i).Einfo(:,2)); %find the event id number from RAW
+            if sum(EvInd)==0
+                fprintf('HOWDY, CANT FIND EVENTS FOR ''%s''\n',Erefnames{k});
+            end
+
+            Lick.Ev(k).NumberTrials(LL,1)=length(RAW(i).Erast{EvInd});
+            if  ~isempty(EvInd) && Lick.Ev(k).NumberTrials(LL,1)>MinNumTrials %avoid analyzing sessions where that do not have enough trials
+                [PSR1,N1]=MakePSR04(RAW(i).Erast(LickInd),RAW(i).Erast{EvInd},Dura,{1});% makes collpased rasters. PSR1 is a cell(neurons)
+                if ~isempty(PSR1{1}) %to avoid errors, added on 12/28 2011
+                    %forcing 100ms bin size to keep it consistent across
+                    %sessions (no reason is should be different for licks)
+                    [PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Erast{LickInd},1),1),{2,0,0.1});%these values force bin size to be 100ms
+                    PTH1=smooth(PTH1,SmoothSPAN,SmoothTYPE)';
+
+                    %------------- Fills the R.Ev(k) fields --------------
+                    Lick.Ev(k).BW(LL,1)=BW1;
+                    Lick.Ev(k).PSTHraw(LL,1:length(Tm))=PTH1;
+                end
+            end
+        end
+    end
+end
+
+Xaxis=[-2 12];
+Ishow=find(Tm>=Xaxis(1) & Tm<=Xaxis(2));
+time1=Tm(Ishow);
+c=[-1000 7500];ClimE=sign(c).*abs(c).^(1/4);%ColorMapExc
+colormap(jet);
+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];
+
+Ev1=strcmp('RD1', Lick.Erefnames);
+Ev2=strcmp('RD2', Lick.Erefnames);
+Ev3=strcmp('RD3', Lick.Erefnames);
+
+psth1=nanmean(Lick.Ev(Ev1).PSTHraw(:,Ishow),1);
+sem1=nanste(Lick.Ev(Ev1).PSTHraw(:,Ishow),1); %calculate standard error of the mean
+up1=psth1+sem1;
+down1=psth1-sem1;
+
+psthE=nanmean(Lick.Ev(Ev2).PSTHraw(:,Ishow),1);
+semE=nanste(Lick.Ev(Ev2).PSTHraw(:,Ishow),1); %calculate standard error of the mean
+upE=psthE+semE;
+downE=psthE-semE;
+
+psth3=nanmean(Lick.Ev(Ev3).PSTHraw(:,Ishow),1);
+sem3=nanste(Lick.Ev(Ev3).PSTHraw(:,Ishow),1); %calculate standard error of the mean
+up3=psth3+sem3;
+down3=psth3-sem3;
+
+%plotting
+subplot(2,3,1);
+hold on;
+plot(time1,psth1,'Color',sucrose,'linewidth',1);
+plot(time1,psthE,'Color',maltodextrin,'linewidth',1);
+plot(time1,psth3,'Color',water,'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)],[upE,downE(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
+patch([time1,time1(end:-1:1)],[up3,down3(end:-1:1)],water,'EdgeColor','none');alpha(0.5);
+title('Mean lick rate');
+%plot([-2 10],[0 0],':','color','k');
+plot([0 0],[-2 8],':','color','k','linewidth',0.75);
+axis([-2 12 0 7.5]);
+xlabel('Seconds from reward delivery');
+ylabel('Licks/s');
+legend('Sucrose trials','Maltodextrin trials','Water trials');
+
+
+%% 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_3R.Ninfo) %runs through neurons
+        if R_3R.BinStatRwrd{i,1}(k).IntSig < PStatBins %if neuron is significant for this bin
+            if R_3R.BinStatRwrd{i-1,1}(k).IntSig < PStatBins || R_3R.BinStatRwrd{i+1,1}(k).IntSig < PStatBins %either previous or subsequent bin must be significant
+                if R_3R.BinStatRwrd{i,1}(k).R1Mean >= R_3R.BinStatRwrd{i,1}(k).R2Mean && R_3R.BinStatRwrd{i,1}(k).R1Mean >= R_3R.BinStatRwrd{i,1}(k).R3Mean %if firing is greater on sucrose trials than maltodextrin and water. if there's a tie, make it sucrose (this is highly unlikely)
+                    R_3R.BinRewPref{i,1}(k,1)=1; %sucrose preferring
+                elseif R_3R.BinStatRwrd{i,1}(k).R2Mean > R_3R.BinStatRwrd{i,1}(k).R1Mean && R_3R.BinStatRwrd{i,1}(k).R2Mean >= R_3R.BinStatRwrd{i,1}(k).R3Mean %if mal is greater than suc and water
+                    R_3R.BinRewPref{i,1}(k,1)=2; %maltodextrin preferring
+                else
+                    R_3R.BinRewPref{i,1}(k,1)=3; %water preferring
+                end
+            else
+                R_3R.BinRewPref{i,1}(k,1)=0; %if not significant in 2 consecutive bins
+            end
+        else
+            R_3R.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_3R.BinRewPref{i,1}(:,1)==1); %sucrose pref
+    NN2perBin(i,1)=sum(R_3R.BinRewPref{i,1}(:,1)==2); %malto pref
+    NN3perBin(i,1)=sum(R_3R.BinRewPref{i,1}(:,1)==3); %malto pref
+    NNperBin(i,1)=sum(R_3R.BinRewPref{i,1}(:,1)>0); %any
+    %normalize to number of neurons in population
+    NN1norm=NN1perBin./length(R_3R.Ninfo);
+    NN2norm=NN2perBin./length(R_3R.Ninfo);
+    NN3norm=NN3perBin./length(R_3R.Ninfo);
+    NNnorm=NNperBin./length(R_3R.Ninfo);
+    
+end
+
+
+%plotting number of significantly modulated neurons across time
+%NAc
+subplot(2,3,2);
+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(xaxis,NN3norm(2:bins-1),'Color',water,'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.85]);
+legend('Total','Suc > mal & wat','Mal > suc & wat','Water > suc & mal','Location','northeast')
+ylabel('Fraction of population');
+xlabel('Seconds from RD');
+title('Reward-selective neurons');
+
+
+%% plotting sucrose-selective neurons
+
+%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
+RD1=strcmp('RD1', R_3R.Erefnames);
+RD2=strcmp('RD2', R_3R.Erefnames);
+RD3=strcmp('RD3', R_3R.Erefnames);
+
+%setting up parameters
+Xaxis=[-2 5];
+inttime=find(R_3R.Param.Tm>=Xaxis(1) & R_3R.Param.Tm<=Xaxis(2));
+paramtime=R_3R.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_3R.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_3R.BinStatRwrd{Bin1+i,1}.SucRespDir)==1; %get neurons with excitation to sucrose
+    Resp12(:,i)=Pref1(:,i)&cat(1,R_3R.BinStatRwrd{Bin1+i,1}.MalRespDir)==1;%get neurons with inhibition to maltodextrin
+    Resp13(:,i)=Pref1(:,i)&cat(1,R_3R.BinStatRwrd{Bin1+i,1}.WatRespDir)==-1;%get neurons with inhibition to maltodextrin
+end
+
+Sel=sum(Pref1,2)>0;  %all neurons selective in any bin
+Sel1=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
+Sel6=sum(Resp13,2)>0; %all selective neurons malto inhibited in any bin
+
+subplot(2,4,5); %heatmap of suc preferring neurons' response to sucrose
+Neurons=R_3R.Ev(RD1).PSTHz(Sel,inttime); %get the firing rates of neurons of interest
+SucResp=cat(1,R_3R.BinStatRwrd{SortBin+1,1}.R1Mean); %sucrose responses
+TMP=SucResp(Sel); %find the magnitude of the excitations for water for this binTMP(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,1)],Neurons,ClimE);
+title(['Sucrose trials']);
+xlabel('Seconds from RD');
+ylabel('Individual neurons');
+hold on;
+plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
+
+
+subplot(2,4,6); %heatmap of suc preferring neurons' response to maltodextrin
+Neurons=R_3R.Ev(RD2).PSTHz(Sel,inttime); %get the firing rates of neurons of interest
+Neurons=Neurons(SORTimg,:); %sort the neurons same as before
+imagesc(paramtime,[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);
+
+subplot(2,4,7); %heatmap of suc preferring neurons' response to water
+Neurons=R_3R.Ev(RD3).PSTHz(Sel,inttime); %get the firing rates of neurons of interest
+Neurons=Neurons(SORTimg,:); %sort the neurons same as before
+imagesc(paramtime,[1,sum(Sel,1)],Neurons,ClimE);
+title(['Water trials']);
+xlabel('Seconds from RD');
+hold on;
+plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
+
+%plot suc preferring neurons to suc
+psthE=nanmean(R_3R.Ev(RD1).PSTHz(Sel,inttime),1); 
+semE=nanste(R_3R.Ev(RD1).PSTHz(Sel,inttime),1); %calculate standard error of the mean
+upE=psthE+semE;
+downE=psthE-semE;
+
+%plot suc preferring neurons to malt
+psth2=nanmean(R_3R.Ev(RD2).PSTHz(Sel,inttime),1); 
+sem2=nanste(R_3R.Ev(RD2).PSTHz(Sel,inttime),1); %calculate standard error of the mean
+up2=psth2+sem2;
+down2=psth2-sem2;
+
+%plot suc preferring neurons to water
+psth3=nanmean(R_3R.Ev(RD3).PSTHz(Sel,inttime),1); 
+sem3=nanste(R_3R.Ev(RD3).PSTHz(Sel,inttime),1); %calculate standard error of the mean
+up3=psth3+sem3;
+down3=psth3-sem3;
+
+%plotting
+subplot(2,3,3);
+hold on;
+plot(paramtime,psthE,'Color',sucrose,'linewidth',1);
+plot(paramtime,psth2,'Color',maltodextrin,'linewidth',1);
+plot(paramtime,psth3,'Color',water,'linewidth',1);
+patch([paramtime,paramtime(end:-1:1)],[upE,downE(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);
+patch([paramtime,paramtime(end:-1:1)],[up3,down3(end:-1:1)],water,'EdgeColor','none');alpha(0.5);
+legend('Suc trials','Mal trials','Wat trials');
+plot([-2 5],[0 0],':','color','k','linewidth',0.75);
+plot([0 0],[-2 8],':','color','k','linewidth',0.75);
+axis([-2 5 -2 4.7]);
+ylabel('Mean firing (z-score)');
+title(['Suc>mal&wat (n=' num2str(sum(Sel)) ' of ' num2str(length(Sel)) ')'])
+xlabel('Seconds from RD');
+
+G = sum(Sel1&Sel3&Sel6);
+F = sum(Sel3&Sel6);
+E = sum(Sel1&Sel6);
+D = sum(Sel1&Sel3);
+A = sum(Sel1);
+B = sum(Sel3);
+C = sum(Sel6);
+
+%vertical venn diagram
+x = [0 0 1 1];
+y1 = [C-E C-E+A C-E+A C-E];
+y2 = [C-F C-F+B C-F+B C-F];
+y3 = [0 C C 0];
+
+subplot(2,35,64);
+hold on;
+s = patch(x,y1,sucrose);
+m = patch(x,y2,maltodextrin);
+w = patch(x,y3,water);
+alpha(s,0.7);
+alpha(w,0.5);
+alpha(m,0.5);
+set(gca,'xtick',[]);
+ylabel('Distribution of neurons');
+axis([0 1 0 C-E+A]);
+
+
+%% stats on reward firing averaged together
+%for simplicity, just looking at average activity in sort bin
+%because I already have that data collected
+
+SucResp=cat(1,R_3R.BinStatRwrd{SortBin+1,1}.R1Mean); %suc responses
+MalResp=cat(1,R_3R.BinStatRwrd{SortBin+1,1}.R2Mean); %mal responses
+WatResp=cat(1,R_3R.BinStatRwrd{SortBin+1,1}.R3Mean); %wat responses
+
+[~,R_3R.RewRespStat{1,1},R_3R.RewRespStat{1,2}]=anovan(cat(1,SucResp(Sel),MalResp(Sel),WatResp(Sel)),{cat(1,zeros(sum(Sel),1),ones(sum(Sel),1),2*ones(sum(Sel),1)),cat(1,R_3R.Ninfo(Sel,4),R_3R.Ninfo(Sel,4),R_3R.Ninfo(Sel,4))},'varnames',{'Reward','Subject'},'random',2,'Display','off','model','full');
+R_3R.RewRespStat{1,3}=multcompare(R_3R.RewRespStat{1,2},'display','off');
+
+save('R_3R.mat','R_3R');