123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687 |
- function stats=plot_TH(Colors, R, RewHist, bm_RD, bm_cue, RAW, os, plot_sub_figures, bm_RD_BIC)
- %% behavior
- if plot_sub_figures
- figure;
- subplot(3,1,1)
- session_to_use = 3;
- session_end = strcmp('MedEnd',RAW(session_to_use(1)).Einfo(:,2)); %finds the index of the event
- stop = RAW(session_to_use(1)).Erast{session_end,1};
- line([0,3+stop(end)/60],[2.3,2.3],'Color','k')
- %get events for this session
- RD1 = strcmp('RD1',RAW(session_to_use).Einfo(:,2)); %finds the index of the event
- RD2 = strcmp('RD2',RAW(session_to_use).Einfo(:,2)); %finds the index of the event
- RD3 = strcmp('RD3',RAW(session_to_use).Einfo(:,2)); %finds the index of the event
- suc_delivery = RAW(session_to_use).Erast{RD1,1}/60;
- malt_delivery = RAW(session_to_use).Erast{RD2,1}/60;
- wat_delivery = RAW(session_to_use).Erast{RD3,1}/60;
- for trial = 1:length(suc_delivery)
- a=line([suc_delivery(trial),suc_delivery(trial)],[2.2,2.8],'Color',Colors('sucrose'));
- end
- for trial = 1:length(malt_delivery)
- b=line([malt_delivery(trial),malt_delivery(trial)],[2,2.6],'Color',Colors('maltodextrin'));
- end
- for trial = 1:length(wat_delivery)
- c=line([wat_delivery(trial),wat_delivery(trial)],[1.8,2.4],'Color',Colors('water'));
- end
- text(90,3,'Sucrose, unpredictable','Color',Colors('sucrose'));
- text(90,2.25,'Maltodextrin, unpredictable','Color',Colors('maltodextrin'));
- text(90,1.5,'Water, unpredictable','Color',Colors('water'));
- xticks([0 20 40 60 80]);
- axis([0 110 0.9 4]);
- ax1 = gca;
- ax1.YAxis.Visible = 'off';
- xlabel('Session time (minutes)')
- %lick plots
- global Dura Tm BSIZE %Tbin
- BSIZE=0.05; %what's the bin for the PSTH?
- Dura=[-4 15]; Tm=Dura(1):BSIZE:Dura(2);
- %Tbin=-0.5:0.005:0.5; %window used to determine the optimal binsize
- %Smoothing PSTHs
- %smoothing filter
- smoothbins=25; %number of previous bins used to smooth
- halfnormal=makedist('HalfNormal','mu',0,'sigma',8); %std=3.98
- filterweights=pdf(halfnormal,[0:smoothbins]);
- %generate lick PSTHs
- for session = 1:length(RAW)
- % Make logical index vectors to call every event type
- R1 = strcmp('RD1',RAW(session).Einfo(:,2)); %finds the index of the event
- R2 = strcmp('RD2',RAW(session).Einfo(:,2)); %finds the index of the event
- R3 = strcmp('RD3',RAW(session).Einfo(:,2)); %finds the index of the event
- Licks = strcmp('Licks',RAW(session).Einfo(:,2)); %finds the index of the event
- EvList{1}=RAW(session).Erast{R1};
- EvList{2}=RAW(session).Erast{R2};
- EvList{3}=RAW(session).Erast{R3};
- for k=1:3
- [PSR1,N1]=MakePSR04(RAW(session).Erast(Licks),EvList{k},Dura,{1});% makes collpased rasters. PSR1 is a cell(neurons)
- [PTH1,~,~]=MakePTH07(PSR1,repmat(N1,size(RAW(session).Erast(Licks),1),1),{2,0,BSIZE});%-----Fixed bin used here
- PTH1smooth=[];
- for l=1:length(Tm)
- PTH1smooth(1,l)=sum(PTH1(1,l-min([l-1 smoothbins]):l).*fliplr(filterweights(1:min([l smoothbins+1]))))/sum(filterweights(1:min([l smoothbins+1])));
- end
- Lick(k).PSTHraw(session,:)=PTH1smooth;
- end
- end
- %plotting
- Xaxis=[-2 12];
- Ishow=find(Tm>=Xaxis(1) & Tm<=Xaxis(2));
- time1=Tm(Ishow);
- Sel=ones(length(RAW),1)==1;
- psth1=nanmean(Lick(1).PSTHraw(Sel,Ishow),1);
- sem1=nanste(Lick(1).PSTHraw(Sel,Ishow),1); %calculate standard error of the mean
- up1=psth1+sem1;
- down1=psth1-sem1;
- psth2=nanmean(Lick(2).PSTHraw(Sel,Ishow),1);
- sem2=nanste(Lick(2).PSTHraw(Sel,Ishow),1); %calculate standard error of the mean
- up2=psth2+sem2;
- down2=psth2-sem2;
- psth3=nanmean(Lick(3).PSTHraw(Sel,Ishow),1);
- sem3=nanste(Lick(3).PSTHraw(Sel,Ishow),1); %calculate standard error of the mean
- up3=psth3+sem3;
- down3=psth3-sem3;
- %plotting
- subplot(3,3,9);
- hold on;
- plot(time1,psth1,'Color',Colors('sucrose'),'linewidth',1);
- plot(time1,psth2,'Color',Colors('maltodextrin'),'linewidth',1);
- plot(time1,psth3,'Color',Colors('water'),'linewidth',1);
- patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],Colors('sucrose'),'EdgeColor','none');alpha(0.5);
- patch([time1,time1(end:-1:1)],[up2,down2(end:-1:1)],Colors('maltodextrin'),'EdgeColor','none');alpha(0.5);
- patch([time1,time1(end:-1:1)],[up3,down3(end:-1:1)],Colors('water'),'EdgeColor','none');alpha(0.5);
- plot([0 0],[-2 8],':','color','k','linewidth',0.75);
- axis([Xaxis(1) Xaxis(2) 0 8]);
- xlabel('Seconds from reward delivery');
- ylabel('Mean lick rate (licks/s)');
- legend('Suc trials','Mal trials','Wat trials');
- end
- %% main figure
- xaxis=[-2 5];
- time_indicies=find(R.Param.Tm>=xaxis(1) & R.Param.Tm<=xaxis(end)); %returns the indicies where these are both true
- activity_xaxis=R.Param.Tm(time_indicies);
- suc_delivery=strcmp('RD1', R.Erefnames);
- malt_delivery=strcmp('RD2', R.Erefnames);
- wat_delivery=strcmp('RD3', R.Erefnames);
- THcolors{1,1}=[1 0.7 0.4];
- THcolors{2,1}=[.8 0.6 0.3];
- THcolors{3,1}=[.5 0.3 0];
- THcolors{4,1}=[1 0.5 0.8];
- THcolors{5,1}=[.8 0.4 .7];
- THcolors{6,1}=[.5 0.1 .3];
- THcolors{7,1}=[0.2 0.8 0.9];
- THcolors{8,1}=[0.1 0.7 0.8];
- THcolors{9,1}=[0 0.5 0.6];
- HistWindow=RewHist.RDWindow;
- CueWindow=RewHist.CueWindow;
- figure;
- %% raster
- alph=0.3;
- %raster
- window=[0 3];
- sample_neuron_number=228; %237, 228
- subplot(1,4,1);
- hold on;
- sessions=unique(R.Ninfo(:,1));
- session=R.Ninfo(sample_neuron_number,1);
- sample_neuron_session=strcmp(session,sessions);
- neurons=R.Ninfo(strcmp(session,R.Ninfo(:,1)),2);
- sample_neuron_session_number=strcmp(R.Ninfo(sample_neuron_number,2),neurons);
- RD=strcmp('RD', RAW(sample_neuron_session).Einfo(:,2));
- RDtimes=RAW(sample_neuron_session).Erast{RD,1};
- [~,order]=sort(os(sample_neuron_number).mod_RD.base.RPEs);
- RDtimes=RDtimes(order);
- patch([HistWindow(1) HistWindow(1) HistWindow(2) HistWindow(2)],[0 length(RDtimes)*11 length(RDtimes)*11 0],Colors('gray'),'edgecolor','none');alpha(alph);
- PlotRaster(RAW(sample_neuron_session).Nrast{sample_neuron_session_number,1},RDtimes,window);
- ylabel('Trials, sorted by RPE (low to high)');
- yticks([]);
- xticks(0:3);
- xticklabels(0:3);
- % plot([HistWindow(1) HistWindow(1)],[0 length(RDtimes)*11],'color',[0 0.5 0.3],'linewidth',0.85);
- % plot([HistWindow(2) HistWindow(2)],[0 length(RDtimes)*11],'color',[0 0.5 0.3],'linewidth',0.85);
- xlabel('Seconds from reward delivery');
- xlim(window);
- disp('Raster')
- %% activity plots
- ymax=5;
- alph=0.3;
- Sel=R.Rat>0; %all neurons
- %plot sucrose preferring neurons to sucrose
- psth1=nanmean(R.Ev(suc_delivery).PSTHz(Sel,time_indicies),1);
- sem1=nanste(R.Ev(suc_delivery).PSTHz(Sel,time_indicies),1); %calculate standard error of the mean
- up1=psth1+sem1;
- down1=psth1-sem1;
- %plot sucrose preferring neurons to malt
- psth2=nanmean(R.Ev(malt_delivery).PSTHz(Sel,time_indicies),1);
- sem2=nanste(R.Ev(malt_delivery).PSTHz(Sel,time_indicies),1); %calculate standard error of the mean
- up2=psth2+sem2;
- down2=psth2-sem2;
- %plot sucrose preferring neurons to malt
- psth3=nanmean(R.Ev(wat_delivery).PSTHz(Sel,time_indicies),1);
- sem3=nanste(R.Ev(wat_delivery).PSTHz(Sel,time_indicies),1); %calculate standard error of the mean
- up3=psth3+sem3;
- down3=psth3-sem3;
- subplot(3,8,3:5)
- hold on
- patch([HistWindow(1) HistWindow(1) HistWindow(2) HistWindow(2)],[-3 ymax ymax -3],Colors('gray'),'edgecolor','none');alpha(alph);
- a=plot(activity_xaxis,psth1,'Color',Colors('sucrose'),'linewidth',1);
- b=plot(activity_xaxis,psth2,'Color',Colors('maltodextrin'),'linewidth',1);
- c=plot(activity_xaxis,psth3,'Color',Colors('water'),'linewidth',1);
- patch([activity_xaxis,activity_xaxis(end:-1:1)],[up1,down1(end:-1:1)],Colors('sucrose'),'EdgeColor','none');alpha(alph);
- patch([activity_xaxis,activity_xaxis(end:-1:1)],[up2,down2(end:-1:1)],Colors('maltodextrin'),'EdgeColor','none');alpha(alph);
- patch([activity_xaxis,activity_xaxis(end:-1:1)],[up3,down3(end:-1:1)],Colors('water'),'EdgeColor','none');alpha(alph);
- plot([-2 5],[0 0],':','color','k','linewidth',0.75);
- plot([0 0],[-3 ymax],'color','k','linewidth',0.75);
- plot([-0.5 -0.5],[-3 ymax],'color','b','linewidth',0.75);
- axis([-2 4 -3 ymax]);
- title(['Reward response, all neurons (n=' num2str(sum(Sel)) ')']);
- legend([a b c],'Suc trials','Mal trials','Wat trials');
- text(-1.3,ymax-0.5,'PE','color','b');
- xlabel('Seconds from reward delivery');
- disp('Interspersed event responses')
- %% categories and linear models
- trialsback=10;
- %neuron categories
- masks=cat(2,bm_RD.mask_base',bm_RD.mask_curr',bm_RD.mask_mean');
- %masks=cat(2,(bm_RD.mask_base' & bm_RD_BIC.mask_base'==0),bm_RD.mask_curr',bm_RD.mask_mean');
- stacked_data = zeros(2,3);
- for subset=1:3
- stacked_data(1,subset)=sum(masks(:,subset));
- end
- subplot(3,4,4)
- hold on;
- b = bar(stacked_data./sum(stacked_data,2),'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])
- text(1.5,0.8,'Mean','color',Colors('mean'));
- text(1.5,0.5,'Current','color',Colors('current'));
- text(1.5,0.2,'RPE','color',Colors('rpe'));
- set(gca,'xtick',[]);
- yticks([0,.2,.4,.6,.8,1])
- ylabel('Fraction of population')
- linecolor{1}=Colors('rpe');
- linecolor{2}=Colors('current');
- linecolor{3}=Colors('mean');
- for i=1:3
- if i==1 Sel=masks(:,1); end %base neurons
- if i==2 Sel=bm_RD.mask_curr'; end
- if i==3 Sel=bm_RD.mask_mean'; end
- %calculate model based on all reward selective neurons
- allRDHz=[];
- allpreds=[];
- selindex=find(Sel);
- for j=1:sum(Sel)
- allRDHz=cat(1,allRDHz,RewHist.RDHz{selindex(j),1});
- allpreds=cat(1,allpreds,RewHist.Predictors{selindex(j),1}(:,1:trialsback+1));
- end
-
- if i==2 allpreds(allpreds==0.75)=0.82; end %set to 0.75, which is best for base neurons, but 0.82 is best for current only
- allVPmdl=fitlm(allpreds,allRDHz);
- stats.trials.subsets{i,1}=allVPmdl.Coefficients.pValue(2:12);
- %plot coefficients
- subplot(3,4,8);
- hold on;
- errorbar(0:trialsback,allVPmdl.Coefficients.Estimate(2:trialsback+2),allVPmdl.Coefficients.SE(2:trialsback+2),'color',linecolor{i},'linewidth',1.5);
-
- end
- 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');
- axis([-1 trialsback+1 -1.2 7]);
- disp('B&G regressions')
- %% RPE neurons based on different RPE values
- bins=[-3 -1.2;-1.2 -0.7;-0.7 -0.3;-0.3 0;0 0.3;0.3 0.7;0.7 1.2;1.2 3];
- Sel=masks(:,1);
- 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);
-
- session=R.Ninfo(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),R.Ninfo(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('Cue',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=10;
- 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,4,6: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) -4 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;-.3;0;0.3;0.7;1.2;3])
- set(get(h,'title'),'string','RPE');
- %legend([a{:}],labels,'location','northeast');
- if plot_sub_figures
- %% cue responses
- masks=cat(2,bm_cue.mask_base',bm_cue.mask_mean');
- figure;
- %% bar graphs
- stacked_data = zeros(2,3);
- for subset=1:2
- stacked_data(1,subset)=sum(masks(:,subset));
- end
- subplot(2,3,1)
- hold on;
- b = bar(stacked_data./sum(stacked_data,2),'stacked','edgecolor','none');
- b(2).FaceColor = 'Flat';
- b(2).CData = Colors('mean');
- b(1).FaceColor = 'Flat';
- b(1).CData = Colors('rpe');
- axis([0 5 0 1])
- set(gca,'xtick',[]);
- yticks([0,.2,.4,.6,.8,1])
- ylabel('Fraction of population')
- legend('Value','Unmod.','location','northwest');
- %% activity plots for subsets
- ymax=9;
- %ymax=5;
- Xaxis=[-1 2];
- Ishow=find(R.Param.Tm>=Xaxis(1) & R.Param.Tm<=Xaxis(2));
- time1=R.Param.Tm(Ishow);
- suc_cue=strcmp('CueP1', R.Erefnames);
- mal_cue=strcmp('CueP2', R.Erefnames);
- wat_cue=strcmp('CueP3', R.Erefnames);
- titles={'Value neurons (n=';'Unmodulated neurons (n='};
- for i=1:2
- if i==1 Sel=bm_cue.mask_base'; end %base neurons
- if i==2 Sel=bm_cue.mask_mean'; end %mean neurons
- %plot cue after suc
- psth1=nanmean(R.Ev(suc_cue).PSTHz(Sel,Ishow),1);
- sem1=nanste(R.Ev(suc_cue).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
- up1=psth1+sem1;
- down1=psth1-sem1;
- %plot cue after malt
- psth2=nanmean(R.Ev(mal_cue).PSTHz(Sel,Ishow),1);
- sem2=nanste(R.Ev(mal_cue).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
- up2=psth2+sem2;
- down2=psth2-sem2;
- %plot cue after water
- psth3=nanmean(R.Ev(wat_cue).PSTHz(Sel,Ishow),1);
- sem3=nanste(R.Ev(wat_cue).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
- up3=psth3+sem3;
- down3=psth3-sem3;
- %make the plot
- subplot(2,3,1+i);
- hold on;
- title([titles{i} num2str(sum(Sel)) ')'])
- patch([CueWindow(1) CueWindow(1) CueWindow(2) CueWindow(2)],[-5 12 12 -5],Colors('gray'),'EdgeColor','none');
- a=plot(time1,psth2,'Color',Colors('maltodextrin'),'linewidth',1);
- b=plot(time1,psth1,'Color',Colors('sucrose'),'linewidth',1);
- c=plot(time1,psth3,'Color',Colors('water'),'linewidth',1);
- patch([time1,time1(end:-1:1)],[up2,down2(end:-1:1)],Colors('maltodextrin'),'EdgeColor','none');alpha(alph);
- patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],Colors('sucrose'),'EdgeColor','none');alpha(alph);
- patch([time1,time1(end:-1:1)],[up3,down3(end:-1:1)],Colors('water'),'EdgeColor','none');alpha(alph);
- plot([-2 5],[0 0],':','color','k','linewidth',0.75);
- plot([0 0],[-2 ymax],'color','k','linewidth',0.5);
- axis([Xaxis(1) Xaxis(2) -2 ymax]);
- ylabel('Mean firing (z-score)');
- xlabel('Seconds from cue');
- if i==2 legend([b a c],'After suc','After mal','After wat','location','northeast'); end
- end
- disp('Subsets event responses')
- %% B&G for cue subsets
- linecolor{1}=Colors('rpe');
- linecolor{2}=Colors('mean');
- for i=1:2
- if i==1 Sel=bm_cue.mask_base'; end %base neurons
- if i==2 Sel=bm_cue.mask_mean'; end %base neurons
- %calculate model based on all reward selective neurons
- allRDHz=[];
- allpreds=[];
- selindex=find(Sel);
- for j=1:sum(Sel)
- allRDHz=cat(1,allRDHz,RewHist.CueHz{selindex(j),1});
- allpreds=cat(1,allpreds,RewHist.Predictors{selindex(j),1}(:,2:trialsback+1));
- end
-
- %allpreds(allpreds==0.75)=0.5;
- allVPmdl=fitlm(allpreds,allRDHz);
- stats.cue_trials.subsets{i,1}=allVPmdl.Coefficients.pValue(2:11);
- %plot coefficients
- subplot(2,3,4);
- hold on;
- errorbar(1:trialsback,allVPmdl.Coefficients.Estimate(2:trialsback+1),allVPmdl.Coefficients.SE(2:trialsback+1),'color',linecolor{i},'linewidth',1.5);
- end
- plot([-1 trialsback+1],[0 0],':','color','k','linewidth',0.75);
- xlabel('Trials back');
- ylabel('Coefficient weight');
- title('Trial history regression, subsets');
- legend('Value','Unmod.');
- axis([0 trialsback+1 -2.5 4]);
- disp('B&G regressions')
- %% Value neurons based on different values
- bins=[-3 -1.2;-1.2 -0.5;-0.5 0.5;0.5 1.2;1.2 3];
- Sel=bm_cue.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',10); %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',5); %std=3.98
- trialfilterweights=pdf(halfnormal,0:trialsmoothbins);
- sessions=unique(R.Ninfo(:,1));
- index=find(Sel);
- RPEPSTHs={};
- for i=1:sum(Sel)
- session=R.Ninfo(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),R.Ninfo(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('Cue',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 cue onset
- cuetimes=RAW(sessnum).Erast{Cue,1}(:,1); %cue times
- RDtimes=RAW(sessnum).Erast{RD,1}(:,1); %reward delivery times
- trialtimes=[];
- for trial = 1:length(RDtimes)
- trialtimes(trial,1) = max(cuetimes(cuetimes<RDtimes(trial)));
- end
- for j=1:length(bins)
- RPEs=os(index(i)).mod_cue.base.V;
- 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=9;
- Xaxis=[-1 2];
- Ishow=find(R.Param.Tm>=Xaxis(1) & R.Param.Tm<=Xaxis(2));
- time1=R.Param.Tm(Ishow);
- alph=0.2;
- psthcolor=[0 0.1 0;...
- 0 0.3 0.2;...
- 0 0.5 0.3;...
- 0 0.75 0.5;...
- 0 1 0.7];
- subplot(2,2,4);
- hold on;
- patch([CueWindow(1) CueWindow(1) CueWindow(2) CueWindow(2)],[-5 12 12 -5],Colors('gray'),'EdgeColor','none');
- 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],[-3 ymax],'color','k','linewidth',0.5);
- axis([Xaxis(1) Xaxis(2) -2 ymax]);
- ylabel('Mean firing (z-score)');
- xlabel('Seconds from cue');
- colormap(psthcolor);
- h=colorbar;
- set(h,'YTick',[0:0.2:1]);
- set(h,'yticklabel',[-3;-1.2;-.5;0.5;1.2;3])
- set(get(h,'title'),'string','Value');
- %legend([a{:}],labels,'location','northeast');
- end
|