123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590 |
- function Stats=plot_cued(Colors, R, RewHist, bm_RD, DOI, ROI, os, bm_cue, RAW)
- %which neurons are plotted
- IncludedNeurons=ChooseNs(R,DOI,ROI);
- %subselect sessions of interest
- for i=1:length(RAW)
- name=char(RAW(i).Ninfo(1,1));
- sessrat(i,1)=str2num(name(2:3));
- sessday(i,1)=str2num(name(8:9));
- end
- %decide which sessions to include
- inc_sessions = (ismember(sessrat,ROI) & ismember(sessday,DOI));
- RAW=RAW(inc_sessions);
- xaxis=[-2 5];
- %event names
- RD1NP=strcmp('Cue3RD1', R.Erefnames);
- RD1PR=strcmp('Cue1RD', R.Erefnames);
- RD2NP=strcmp('Cue3RD2', R.Erefnames);
- RD2PR=strcmp('Cue2RD', R.Erefnames);
- Cue1=strcmp('Cue1', R.Erefnames);
- Cue2=strcmp('Cue2', R.Erefnames);
- Cue31=strcmp('Cue31', R.Erefnames);
- Cue32=strcmp('Cue32', R.Erefnames);
- %when the model data is only on subselected neurons, make the logicals in terms of all neurons
- NNumber=[1:length(R.Ninfo)]';
- index=NNumber(IncludedNeurons);
- histmasks=zeros(length(R.Ninfo),3);
- histmasks(index,1)=bm_RD.mask_base'|bm_RD.mask_base_cue';
- histmasks(index,2)=bm_RD.mask_curr'|bm_RD.mask_curr_cue';
- histmasks(index,3)=bm_RD.mask_mean'|bm_RD.mask_mean_cue';
- cuemasks=zeros(length(R.Ninfo),2);
- cuemasks(index,1)=bm_RD.mask_base_cue'|bm_RD.mask_curr_cue'|bm_RD.mask_mean_cue';
- cuemasks(index,2)=bm_RD.mask_base'|bm_RD.mask_curr'|bm_RD.mask_mean';
- %analysis windows
- HistWindow = [0.75 1.95];
- CueWindow = [0 0.75];
- figure;
- %% bar graph for history effects
- stacked_data = zeros(2,3);
- for subset=1:3
- stacked_data(1,subset)=sum(IncludedNeurons&histmasks(:,subset)&cuemasks(:,2));
- stacked_data(2,subset)=sum(IncludedNeurons&histmasks(:,subset)&cuemasks(:,1));
- end
- subplot(3,6,1)
- hold on;
- b = bar([1 2],stacked_data./sum(IncludedNeurons),'stacked','edgecolor','none');
- b(3).FaceColor = 'Flat';
- b(3).CData = Colors('mean');
- b(2).FaceColor = 'Flat';
- b(2).CData = Colors('current');
- b(1).FaceColor = 'Flat';
- b(1).CData = Colors('rpe');
- axis([0 7 0 1])
- legend('RPE','Current','Unmod.');
- set(gca,'xtick',[]);
- yticks([0,.2,.4,.6,.8,1])
- ylabel('Fraction of population')
- %chi2 for distribution of cue effects among history effects
- for i = 1:length(histmasks)
- histcategory(i,1)=sum(find(histmasks(i,:)));
- cuecategory(i,1)=sum(find(cuemasks(i,:)));
- end
- [table,Stats.HistCueDistX2,Stats.HistCueDistpval]=crosstab(histcategory(IncludedNeurons),cuecategory(IncludedNeurons));
- %% mean activity
- Sel=IncludedNeurons;
- ymax=5;
- alph=0.3;
- Xaxis=[-2 4];
- Ishow=find(R.Param.Tm>=Xaxis(1) & R.Param.Tm<=Xaxis(2));
- time1=R.Param.Tm(Ishow);
- %plot suc after non-specific cue
- psth1=nanmean(R.Ev(RD1NP).PSTHz(Sel,Ishow),1);
- sem1=nanste(R.Ev(RD1NP).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
- up1=psth1+sem1;
- down1=psth1-sem1;
- %plot suc after sucrose cue
- psth2=nanmean(R.Ev(RD1PR).PSTHz(Sel,Ishow),1);
- sem2=nanste(R.Ev(RD1PR).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
- up2=psth2+sem2;
- down2=psth2-sem2;
- %plot malt after non-specific cue
- psth3=nanmean(R.Ev(RD2NP).PSTHz(Sel,Ishow),1);
- sem3=nanste(R.Ev(RD2NP).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
- up3=psth3+sem3;
- down3=psth3-sem3;
- %plot malt after maltodextrin cue
- psth4=nanmean(R.Ev(RD2PR).PSTHz(Sel,Ishow),1);
- sem4=nanste(R.Ev(RD2PR).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
- up4=psth4+sem4;
- down4=psth4-sem4;
- %make the plot
- subplot(3,3,4);
- hold on;
- title(['All neurons (n=' num2str(sum(Sel)) ')']);
- patch([HistWindow(1) HistWindow(1) HistWindow(2) HistWindow(2)],[-5 12 12 -5],Colors('gray'),'EdgeColor','none');
- a=plot(time1,psth1,'Color',Colors('sucrose'),'linewidth',1);
- b=plot(time1,psth2,'Color',Colors('blockrose'),'linewidth',1);
- c=plot(time1,psth3,'Color',Colors('maltodextrin'),'linewidth',1);
- d=plot(time1,psth4,'Color',Colors('blockrodextrin'),'linewidth',1);
- patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],Colors('sucrose'),'EdgeColor','none');alpha(alph);
- patch([time1,time1(end:-1:1)],[up2,down2(end:-1:1)],Colors('blockrose'),'EdgeColor','none');alpha(alph);
- patch([time1,time1(end:-1:1)],[up3,down3(end:-1:1)],Colors('maltodextrin'),'EdgeColor','none');alpha(alph);
- patch([time1,time1(end:-1:1)],[up4,down4(end:-1:1)],Colors('blockrodextrin'),'EdgeColor','none');alpha(alph);
- plot([-2 5],[0 0],':','color','k','linewidth',0.75);
- plot([0 0],[-3 ymax],'color','k','linewidth',0.5);
- text(0.1,ymax-0.5,'RD','color','k');
- axis([Xaxis(1) Xaxis(2) -0.5 ymax]);
- ylabel('Mean firing (z-score)');
- xlabel('Seconds from RD');
- legend([a b c d],'Non-pr suc','Pred suc','Non-pr mal','Pred mal','location','northeast');
- disp('Mean activity');
- %% reward history plots
- %broken down by previous trial
- ymax=10;
- alph=0.3;
- trialsback=10;
- linecolor{1}=Colors('rpe');
- linecolor{2}=Colors('current');
- linecolor{3}=Colors('mean');
- titles={'RPE neurons (n=';'Current outcome neurons (n=';'Mean neurons (n='};
- for subset=1:3
- Sel=IncludedNeurons&histmasks(:,subset)&cuemasks(:,2);
-
- %calculate model based on all reward selective neurons
- allRDHz=[];
- allpreds=[];
- selindex=find(Sel);
- for i=1:sum(Sel)
- allRDHz=cat(1,allRDHz,RewHist.RDHz{selindex(i),1});
- allpreds=cat(1,allpreds,RewHist.Predictors{selindex(i),1}(:,1:trialsback+1));
- end
- allVPmdl=fitlm(allpreds,allRDHz);
- Stats.trials.subsets{subset,1}=allVPmdl.Coefficients.pValue(2:12);
- %plot coefficients
- subplot(3,6,2);
- hold on;
- errorbar(0:trialsback,allVPmdl.Coefficients.Estimate(2:trialsback+2),allVPmdl.Coefficients.SE(2:trialsback+2),'color',linecolor{subset},'linewidth',1.5);
-
- end
- axis([-1 trialsback+1 -1.5 4]);
- plot([-1 trialsback+1],[0 0],':','color','k','linewidth',0.75);
- xlabel('Trials back');
- ylabel('Coefficient weight');
- title('Trial history regression, subsets');
- legend('RPE','Current','Mean');
- disp('Rew History Responses')
- %% scatterplot of cue values for effects on RD
- subplot(3,6,5);
- sz=25;
- hold on;
- sucvalues=nan(sum(IncludedNeurons),1);
- malvalues=nan(sum(IncludedNeurons),1);
- %mean
- Sel=bm_RD.mask_mean_cue;
- index=find(Sel);
- vsuc = getParameters_ott(os, 'RD', 'mean_cue', 'vsuc');
- vmal = getParameters_ott(os, 'RD', 'mean_cue', 'vmal');
- scatter(-vsuc(Sel),-vmal(Sel),sz,Colors('mean'));
- sucvalues(index)=-vsuc(Sel)';
- malvalues(index)=-vmal(Sel)';
- %current
- Sel=bm_RD.mask_curr_cue;
- index=find(Sel);
- vsuc = getParameters_ott(os, 'RD', 'curr_cue', 'vsuc');
- vmal = getParameters_ott(os, 'RD', 'curr_cue', 'vmal');
- scatter(-vsuc(Sel),-vmal(Sel),sz,Colors('current'));
- sucvalues(index)=-vsuc(Sel)';
- malvalues(index)=-vmal(Sel)';
- %base
- Sel=bm_RD.mask_base_cue;
- index=find(Sel);
- vsuc = getParameters_ott(os, 'RD', 'base_cue', 'vsuc');
- vmal = getParameters_ott(os, 'RD', 'base_cue', 'vmal');
- scatter(-vsuc(Sel),-vmal(Sel),sz,Colors('rpe'));
- sucvalues(index)=-vsuc(Sel)';
- malvalues(index)=-vmal(Sel)';
- axis([-1 1 -1 1]);
- xlabel('Sucrose value');
- ylabel('Maltodextrin value');
- plot([0 0],[-1 1],'color','k');
- plot([-1 1],[0 0],'color','k');
- quadrant=[];
- for i=1:length(sucvalues)
- if sucvalues(i)<0 & malvalues(i)>0
- quadrant(i)=1;
- elseif sucvalues(i)>0 & malvalues(i)>0
- quadrant(i)=2;
- elseif sucvalues(i)<0 & malvalues(i)<0
- quadrant(i)=3;
- elseif sucvalues(i)>0 & malvalues(i)<0
- quadrant(i)=4;
- end
- end
- quadrant=quadrant(find(quadrant));
- text(-0.9,0.8,[num2str(round(100*sum(quadrant==1)/length(quadrant),1)) '%']);
- text(0.6,0.8,[num2str(round(100*sum(quadrant==2)/length(quadrant),1)) '%']);
- text(-0.9,-0.8,[num2str(round(100*sum(quadrant==3)/length(quadrant),1)) '%']);
- text(0.6,-0.8,[num2str(round(100*sum(quadrant==4)/length(quadrant),1)) '%']);
- [h,Stats.RDpval,Stats.RDX2] = chi2gof(quadrant,'Ctrs',[1 2 3 4],'expected',[length(quadrant)/4 length(quadrant)/4 length(quadrant)/4 length(quadrant)/4]);
- if h
- text(0.5,-0.8,'*');
- end
- %% effect of cue plots
- AllSel=zeros(length(IncludedNeurons),1);
- index=find(IncludedNeurons);
- Sel=sucvalues>0 & malvalues<0;
- AllSel(index)=Sel;
- Sel=logical(AllSel);
- ymax=7;
- alph=0.3;
- trialsback=10;
- Xaxis=[-0.5 2.5];
- Ishow=find(R.Param.Tm>=Xaxis(1) & R.Param.Tm<=Xaxis(2));
- time1=R.Param.Tm(Ishow);
- %plot sucrose after non-specific cue
- psth1=nanmean(R.Ev(RD1NP).PSTHz(Sel,Ishow),1);
- sem1=nanste(R.Ev(RD1NP).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
- up1=psth1+sem1;
- down1=psth1-sem1;
- %plot sucrose after sucrose cue
- psth2=nanmean(R.Ev(RD1PR).PSTHz(Sel,Ishow),1);
- sem2=nanste(R.Ev(RD1PR).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
- up2=psth2+sem2;
- down2=psth2-sem2;
- %plot malt after non-specific cue
- psth3=nanmean(R.Ev(RD2NP).PSTHz(Sel,Ishow),1);
- sem3=nanste(R.Ev(RD2NP).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
- up3=psth3+sem3;
- down3=psth3-sem3;
- %plot malt after maltodextrin cue
- psth4=nanmean(R.Ev(RD2PR).PSTHz(Sel,Ishow),1);
- sem4=nanste(R.Ev(RD2PR).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
- up4=psth4+sem4;
- down4=psth4-sem4;
- %make the plot
- subplot(3,6,6);
- hold on;
- title(['Suc>0, Mal<0 (n=' num2str(sum(Sel)) ')']);
- patch([HistWindow(1) HistWindow(1) HistWindow(2) HistWindow(2)],[-5 12 12 -5],Colors('gray'),'EdgeColor','none');
- a=plot(time1,psth1,'Color',Colors('sucrose'),'linewidth',1);
- b=plot(time1,psth2,'Color',Colors('blockrose'),'linewidth',1);
- c=plot(time1,psth3,'Color',Colors('maltodextrin'),'linewidth',1);
- d=plot(time1,psth4,'Color',Colors('blockrodextrin'),'linewidth',1);
- patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],Colors('sucrose'),'EdgeColor','none');alpha(alph);
- patch([time1,time1(end:-1:1)],[up2,down2(end:-1:1)],Colors('blockrose'),'EdgeColor','none');alpha(alph);
- patch([time1,time1(end:-1:1)],[up3,down3(end:-1:1)],Colors('maltodextrin'),'EdgeColor','none');alpha(alph);
- patch([time1,time1(end:-1:1)],[up4,down4(end:-1:1)],Colors('blockrodextrin'),'EdgeColor','none');alpha(alph);
- plot([-2 5],[0 0],':','color','k','linewidth',0.75);
- plot([0 0],[-3 ymax],'color','k','linewidth',0.5);
- text(0.1,ymax-0.5,'RD','color','k');
- axis([Xaxis(1) Xaxis(2) -2 ymax]);
- ylabel('Mean firing (z-score)');
- xlabel('Seconds from RD');
- legend([a b c d],'Non-pr suc','Pred suc','Non-pr mal','Pred mal','location','northeast');
- disp('RD cue impact');
- %% bar graph for cue effects on cue
- %when the model data is only on subselected neurons
- NNumber=[1:length(R.Ninfo)]';
- index=NNumber(IncludedNeurons);
- cuemasksc=zeros(length(R.Ninfo),2);
- cuemasksc(index,1)=bm_cue.mask_base_cue'|bm_cue.mask_mean_cue';
- cuemasksc(index,2)=bm_cue.mask_base'|bm_cue.mask_mean';
- stacked_data = zeros(2,2);
- for subset=1:2
- stacked_data(1,subset)=sum(IncludedNeurons&cuemasksc(:,subset));
- end
- subplot(3,6,10)
- hold on;
- b = bar(stacked_data./sum(stacked_data,2),'stacked','edgecolor','none');
- b(2).FaceColor = 'Flat';
- b(2).CData = [0 0 0];
- b(1).FaceColor = 'Flat';
- b(1).CData = Colors('cueeffect');
- axis([0 7 0 1])
- text(1.5,0.7,'No cue effect','color',[0 0 0]);
- %text(1.5,0.15,'Info','color',[0.3 0 0.8]);
- text(1.5,0.2,'Cue effect','color',Colors('cueeffect'));
- set(gca,'xtick',[]);
- yticks([0,.2,.4,.6,.8,1])
- ylabel('Fraction of population')
- %% scatterplot of cue responses
- subplot(3,6,11);
- sz=25;
- hold on;
- sucvalues=nan(sum(IncludedNeurons),1);
- malvalues=nan(sum(IncludedNeurons),1);
- %mean
- Sel=bm_cue.mask_mean_cue;
- index=find(Sel);
- vsuc = getParameters_ott(os, 'cue', 'mean_cue', 'vsuc');
- vmal = getParameters_ott(os, 'cue', 'mean_cue', 'vmal');
- scatter(vsuc(Sel),vmal(Sel),sz,Colors('cueeffect'));
- sucvalues(index)=vsuc(Sel)';
- malvalues(index)=vmal(Sel)';
- %base
- Sel=bm_cue.mask_base_cue;
- index=find(Sel);
- vsuc = getParameters_ott(os, 'cue', 'base_cue', 'vsuc');
- vmal = getParameters_ott(os, 'cue', 'base_cue', 'vmal');
- scatter(vsuc(Sel),vmal(Sel),sz,Colors('cueeffect'));
- sucvalues(index)=vsuc(Sel)';
- malvalues(index)=vmal(Sel)';
- axis([-1 1 -1 1]);
- xlabel('Sucrose value');
- ylabel('Maltodextrin value');
- plot([0 0],[-1 1],'color','k');
- plot([-1 1],[0 0],'color','k');
- quadrant=[];
- for i=1:length(sucvalues)
- if sucvalues(i)<0 & malvalues(i)>0
- quadrant(i)=1;
- elseif sucvalues(i)>0 & malvalues(i)>0
- quadrant(i)=2;
- elseif sucvalues(i)<0 & malvalues(i)<0
- quadrant(i)=3;
- elseif sucvalues(i)>0 & malvalues(i)<0
- quadrant(i)=4;
- end
- end
- quadrant=quadrant(find(quadrant));
- text(-0.9,0.8,[num2str(round(100*sum(quadrant==1)/length(quadrant),1)) '%']);
- text(0.6,0.8,[num2str(round(100*sum(quadrant==2)/length(quadrant),1)) '%']);
- text(-0.9,-0.8,[num2str(round(100*sum(quadrant==3)/length(quadrant),1)) '%']);
- text(0.6,-0.8,[num2str(round(100*sum(quadrant==4)/length(quadrant),1)) '%']);
- [h,Stats.Cuepval,Stats.CueX2] = chi2gof(quadrant,'Ctrs',[1 2 3 4],'expected',[length(quadrant)/4 length(quadrant)/4 length(quadrant)/4 length(quadrant)/4]);
- if h
- text(0.5,-0.8,'*');
- end
- %% effect of cue plots
- AllSel=zeros(length(IncludedNeurons),1);
- index=find(IncludedNeurons);
- Sel=sucvalues>0 & malvalues<0;
- AllSel(index)=Sel;
- Sel=logical(AllSel);
- ymax=8;
- alph=0.3;
- trialsback=10;
- Xaxis=[-0.5 1.5];
- Ishow=find(R.Param.Tm>=Xaxis(1) & R.Param.Tm<=Xaxis(2));
- time1=R.Param.Tm(Ishow);
- %plot non-specific cue, sucrose trials
- psth1=nanmean(R.Ev(Cue31).PSTHz(Sel,Ishow),1);
- sem1=nanste(R.Ev(Cue31).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
- up1=psth1+sem1;
- down1=psth1-sem1;
- %plot sucrose cue
- psth2=nanmean(R.Ev(Cue1).PSTHz(Sel,Ishow),1);
- sem2=nanste(R.Ev(Cue1).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
- up2=psth2+sem2;
- down2=psth2-sem2;
- %plot non-specific cue, maltodextrin trials
- psth3=nanmean(R.Ev(Cue32).PSTHz(Sel,Ishow),1);
- sem3=nanste(R.Ev(Cue32).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
- up3=psth3+sem3;
- down3=psth3-sem3;
- %plot maltodextrin cue
- psth4=nanmean(R.Ev(Cue2).PSTHz(Sel,Ishow),1);
- sem4=nanste(R.Ev(Cue2).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
- up4=psth4+sem4;
- down4=psth4-sem4;
- %make the plot
- subplot(3,6,12);
- hold on;
- title(['Suc>0, Mal<0 (n=' num2str(sum(Sel)) ')']);
- patch([CueWindow(1) CueWindow(1) CueWindow(2) CueWindow(2)],[-5 12 12 -5],Colors('gray'),'EdgeColor','none');
- a=plot(time1,psth1,'Color',Colors('sucrose'),'linewidth',1);
- b=plot(time1,psth2,'Color',Colors('blockrose'),'linewidth',1);
- c=plot(time1,psth3,'Color',Colors('maltodextrin'),'linewidth',1);
- d=plot(time1,psth4,'Color',Colors('blockrodextrin'),'linewidth',1);
- patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],Colors('sucrose'),'EdgeColor','none');alpha(alph);
- patch([time1,time1(end:-1:1)],[up2,down2(end:-1:1)],Colors('blockrose'),'EdgeColor','none');alpha(alph);
- patch([time1,time1(end:-1:1)],[up3,down3(end:-1:1)],Colors('maltodextrin'),'EdgeColor','none');alpha(alph);
- patch([time1,time1(end:-1:1)],[up4,down4(end:-1:1)],Colors('blockrodextrin'),'EdgeColor','none');alpha(alph);
- plot([-2 5],[0 0],':','color','k','linewidth',0.75);
- plot([0 0],[-3 ymax],'color','r','linewidth',0.5);
- text(Xaxis(1)+0.1,ymax-0.5,'Cue','color','r');
- axis([Xaxis(1) Xaxis(2) -2 ymax]);
- ylabel('Mean firing (z-score)');
- xlabel('Seconds from cue');
- legend([a b c d],'Non-pr suc','Pred suc','Non-pr mal','Pred mal','location','northeast');
- disp('Cue cue impact')
- %% RPE neurons based on different RPE values
- bins=[-3 -1.2;-1.2 -0.7;-0.7 -0.4;-0.4 0;0 0.4;0.4 0.7;0.7 1.2;1.2 3];
- Sel=bm_RD.mask_base';
- mintrials=1;
- global Dura Tm BSIZE
- BSIZE=0.01; %Do not change
- Dura=[-4 10]; Tm=Dura(1):BSIZE:Dura(2); %should match whatever was used to create R for simplicity
- Baseline=[-11 -1]; %for PSTHs, should match whatever is in R
- %smoothing filter for whole PSTH
- PSTHsmoothbins=50; %number of previous bins used to smooth
- halfnormal=makedist('HalfNormal','mu',0,'sigma',20); %std=3.98
- PSTHfilterweights=pdf(halfnormal,0:PSTHsmoothbins);
- %smoothing filter for individual trials
- trialsmoothbins=25; %number of previous bins used to smooth
- halfnormal=makedist('HalfNormal','mu',0,'sigma',10); %std=3.98
- trialfilterweights=pdf(halfnormal,0:trialsmoothbins);
- sessions=unique(R.Ninfo(:,1));
- index=find(Sel);
- RPEs=[];
- for i=1:sum(Sel)
-
- %collect all RPE values for all cells of interest
- RPEs=cat(1,RPEs,os(index(i)).mod_RD.base.RPEs);
-
- neuroninfo = R.Ninfo(IncludedNeurons,:);
- session=neuroninfo(index(i),1); %what's the name of session
- sessnum=ismember(sessions,session); %which session of all sessions
- neuronnum=ismember(R.Ninfo(strcmp(session,R.Ninfo(:,1)),2),neuroninfo(index(i),2)); %which number neuron in the session
- RD=strmatch('RD',RAW(sessnum).Einfo(:,2),'exact'); %get index for reward delivery times
- Cue=strmatch('Cues',RAW(sessnum).Einfo(:,2),'exact'); %get index for reward delivery times
-
- %get mean baseline firing for all reward trials
- [Bcell1,B1n]=MakePSR04(RAW(sessnum).Nrast(neuronnum),RAW(sessnum).Erast{Cue},Baseline,{2});% makes trial by trial rasters for baseline
- basespk=[];
- for y= 1:B1n
- basespk(1,y)=sum(Bcell1{1,y}>Baseline(1));
- end
- FiringRate=basespk/(Baseline(1,2)-Baseline(1,1));
- Bmean=nanmean(FiringRate);
- Bstd=nanstd(FiringRate);
-
- %make PSTHs around reward delivery
- trialtimes=RAW(sessnum).Erast{RD,1}(:,1); %reward delivery times
- for j=1:length(bins)
- RPEs=os(index(i)).mod_RD.base.RPEs;
- RPEsZ = (RPEs - mean(RPEs))/std(RPEs);
- toi=trialtimes(RPEsZ>=bins(j,1) & RPEsZ<=bins(j,2)); %which trials, binned by model-produced RPE
- if length(toi)>=mintrials
- %make PSTHs
- [PSR1,~]=MakePSR04(RAW(sessnum).Nrast(neuronnum),toi,Dura,{2});% makes collpased rasters. PSR1 is a cell(neurons)
- smoothedtrials=[];
- for trial=1:length(PSR1)
- binned=histcounts(PSR1{trial},[Tm Tm(end)+(Tm(end)-Tm(end-1))]);
- for l=1:length(Tm)
- smoothedtrials(trial,l)=sum(binned(1,l-min([l-1 trialsmoothbins]):l).*fliplr(trialfilterweights(1:min([l trialsmoothbins+1]))))/sum(trialfilterweights(1:min([l trialsmoothbins+1])));
- end
- end
- PTH1=mean(smoothedtrials,1)/BSIZE;
- PTH1smooth=[];
- for l=1:length(Tm)
- PTH1smooth(1,l)=sum(PTH1(1,l-min([l-1 PSTHsmoothbins]):l).*fliplr(PSTHfilterweights(1:min([l PSTHsmoothbins+1]))))/sum(PSTHfilterweights(1:min([l PSTHsmoothbins+1])));
- end
- RPEPSTHs{j,1}.PSTHraw(i,:)=PTH1smooth;
- %normalize already smoothed activity
- for l=1:length(PTH1smooth)
- RPEPSTHs{j,1}.PSTHz(i,l)=(PTH1smooth(l)-Bmean)/Bstd;
- end
- else
- RPEPSTHs{j,1}.PSTHraw(i,1:length(Tm))=NaN;
- RPEPSTHs{j,1}.PSTHz(i,1:length(Tm))=NaN;
- end
- end
- end
- %plotting the histograms
- ymax=8;
- Xaxis=[-1.5 4];
- Ishow=find(R.Param.Tm>=Xaxis(1) & R.Param.Tm<=Xaxis(2));
- time1=R.Param.Tm(Ishow);
- alph=0.2;
- psthcolor=[0 0.9 1;...
- 0 0.7 0.75;...
- 0 0.5 0.5;...
- 0 0.3 0.25;...
- 0.25 0.3 0;...
- 0.5 0.5 0;...
- 0.75 0.7 0;...
- 1 0.9 0];
- subplot(3,3,7);
- hold on;
- patch([HistWindow(1) HistWindow(1) HistWindow(2) HistWindow(2)],[-4 ymax ymax -4],Colors('gray'),'edgecolor','none');alpha(alph);
- a=[];
- for i=1:length(RPEPSTHs) %for each bin
-
- %plot activity for each bin
- psth1=nanmean(RPEPSTHs{i,1}.PSTHz(:,Ishow),1);
- sem1=nanste(RPEPSTHs{i,1}.PSTHz(:,Ishow),1); %calculate standard error of the mean
- up1=psth1+sem1;
- down1=psth1-sem1;
- %make the plot
- a{i}=plot(time1,psth1,'Color',psthcolor(i,:),'linewidth',1.5);
- patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],psthcolor(i,:),'EdgeColor','none');alpha(alph);
- labels{i}=[num2str(bins(i,1)) ' < RPE < ' num2str(bins(i,2))];
- end
- plot([-2 5],[0 0],':','color','k','linewidth',0.75);
- plot([0 0],[-5 ymax],'color','k','linewidth',0.5);
- axis([Xaxis(1) Xaxis(2) -1 ymax]);
- ylabel('Mean firing (z-score)');
- xlabel('Seconds from RD');
- colormap(psthcolor);
- h=colorbar;
- set(h,'YTick',[0:0.125:1]);
- set(h,'yticklabel',[-3;-1.2;-.7;-.4;0;0.4;0.7;1.2;3])
- set(get(h,'title'),'string','RPE');
- %legend([a{:}],labels,'location','northeast');
- disp('RPE responses')
|