i_FigS7CuePERewHist.m 29 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624
  1. %look at impact of previous trials on Cue and PE
  2. %Looking at the average firing rate for a given window in each of 4
  3. %current/previous reward conditions
  4. clear all;
  5. load ('R_2R.mat');
  6. load ('RAW.mat');
  7. %run linear model and stats? 1 is yes, 0 is no
  8. runanalysis=0;
  9. %divide neurons up by region
  10. NAneurons=strcmp(R_2R.Ninfo(:,3),'NA');
  11. VPneurons=strcmp(R_2R.Ninfo(:,3),'VP');
  12. %get parameters
  13. trialsback=6; %how many trials back to look
  14. PEBaseline=R_2R.Param.BinBase+0.5; %For normalizing activity
  15. CueBaseline=[-11 -1];
  16. CueWindow=[0 0.5]; %what window of activity is analyzed
  17. PEWindow=[-0.6 0.7];
  18. BinDura=R_2R.Param.BinDura;
  19. bins=R_2R.Param.bins;
  20. binint=R_2R.Param.binint;
  21. binstart=R_2R.Param.binstart;
  22. %sorting bin -- which bin the neurons' activity is sorted on for heatmap(in seconds)
  23. SortBinTime=1; %seconds
  24. SortBin=(((SortBinTime-BinDura(2)/2)-binstart)/binint); %convert to bin name
  25. %reset
  26. NN=0;
  27. PEEvMeanz=0;
  28. if runanalysis==1
  29. for i=1:length(RAW) %loops through sessions
  30. if strcmp('NA',RAW(i).Type(1:2)) | strcmp('VP',RAW(i).Type(1:2)) %only look at suc v mal sessions
  31. %events
  32. EV3=strmatch('CueP1',RAW(i).Einfo(:,2),'exact');
  33. EV4=strmatch('CueP2',RAW(i).Einfo(:,2),'exact');
  34. EV1=strmatch('PEP1',RAW(i).Einfo(:,2),'exact');
  35. EV2=strmatch('PEP2',RAW(i).Einfo(:,2),'exact');
  36. Cue=strmatch('Cue',RAW(i).Einfo(:,2),'exact');
  37. PE=strmatch('PECue',RAW(i).Einfo(:,2),'exact');
  38. R1=strmatch('Reward1Deliv',RAW(i).Einfo(:,2),'exact');
  39. R2=strmatch('Reward2Deliv',RAW(i).Einfo(:,2),'exact');
  40. %% linear model for impact of previous rewards
  41. %reset
  42. Xcue=[];
  43. Xpe=[];
  44. Y=[];
  45. %set up the matrix with outcome identity for each session
  46. rewards1=cat(2,RAW(i).Erast{R1,1}(:,1),ones(length(RAW(i).Erast{R1,1}(:,1)),1));
  47. rewards2=cat(2,RAW(i).Erast{R2,1}(:,1),zeros(length(RAW(i).Erast{R2,1}(:,1)),1));
  48. rewards=cat(1,rewards1,rewards2);
  49. [rewards(:,1),ind]=sort(rewards(:,1));
  50. rewards(:,2)=rewards(ind,2);
  51. %cue
  52. 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
  53. for k=firstcue+1:length(RAW(i).Erast{Cue,1}(:,1))
  54. time=RAW(i).Erast{Cue,1}(k,1);
  55. entry=find(rewards==max(rewards(rewards(:,1)<time)));
  56. for m=1:trialsback
  57. Xcue(k-trialsback,m)=rewards(entry+1-m,2);
  58. end
  59. end
  60. %PE
  61. for k=trialsback+1:length(RAW(i).Erast{PE,1}(:,1))
  62. time=RAW(i).Erast{PE,1}(k,1);
  63. entry=find(rewards==max(rewards(rewards(:,1)<time)));
  64. for m=1:trialsback
  65. Xpe(k-trialsback,m)=rewards(entry+1-m,2);
  66. end
  67. end
  68. for j= 1:length(RAW(i).Nrast) %Number of neurons within sessions
  69. NN=NN+1;
  70. %cue
  71. rewspk=0;
  72. basespk=0;
  73. %get mean baseline firing for all cues
  74. [Bcell1,B1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{Cue},CueBaseline,{2});% makes trial by trial rasters for baseline
  75. for y= 1:B1n
  76. basespk(1,y)=sum(Bcell1{1,y}>CueBaseline(1));
  77. end
  78. Bhz=basespk/(CueBaseline(1,2)-CueBaseline(1,1));
  79. Bmean=nanmean(Bhz);
  80. Bstd=nanstd(Bhz);
  81. %get trial by trial firing rate for all cue trials
  82. [EVcell,EVn]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{Cue},CueWindow,{2});% makes trial by trial rasters for event
  83. for y= 1:EVn
  84. rewspk(y,1)=sum(EVcell{1,y}>CueWindow(1));
  85. end
  86. Y=((rewspk(trialsback+1:end,1)/(CueWindow(1,2)-CueWindow(1,1)))-Bmean)/Bstd; %normalize the activity to baseline
  87. %true data
  88. linmdl{NN,1}=fitlm(Xcue,Y);
  89. R_2R.RewHist.LinMdlCueWeights(NN,1:trialsback)=linmdl{NN,1}.Coefficients.Estimate(2:trialsback+1)';
  90. R_2R.RewHist.LinMdlCuePVal(NN,1:trialsback)=linmdl{NN,1}.Coefficients.pValue(2:trialsback+1)';
  91. %shuffled
  92. YSh=Y(randperm(length(Y)));
  93. linmdlSh{NN,1}=fitlm(Xcue,YSh);
  94. R_2R.RewHist.LinMdlCueWeightsSh(NN,1:trialsback)=linmdlSh{NN,1}.Coefficients.Estimate(2:trialsback+1)';
  95. R_2R.RewHist.LinMdlCuePValSh(NN,1:trialsback)=linmdlSh{NN,1}.Coefficients.pValue(2:trialsback+1)';
  96. %PE
  97. rewspk=0;
  98. basespk=0;
  99. %get mean baseline firing for all PEs
  100. [Bcell1,B1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{PE},PEBaseline,{2});% makes trial by trial rasters for baseline
  101. for y= 1:B1n
  102. basespk(1,y)=sum(Bcell1{1,y}>PEBaseline(1));
  103. end
  104. Bhz=basespk/(PEBaseline(1,2)-PEBaseline(1,1));
  105. Bmean=nanmean(Bhz);
  106. Bstd=nanstd(Bhz);
  107. %get trial by trial firing rate for all PE trials
  108. [EVcell,EVn]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{PE},PEWindow,{2});% makes trial by trial rasters for event
  109. for y= 1:EVn
  110. rewspk(y,1)=sum(EVcell{1,y}>PEWindow(1));
  111. end
  112. Y=((rewspk(trialsback+1:end,1)/(PEWindow(1,2)-PEWindow(1,1)))-Bmean)/Bstd; %normalize the activity to baseline
  113. %true data
  114. linmdl{NN,1}=fitlm(Xpe,Y);
  115. R_2R.RewHist.LinMdlPEWeights(NN,1:trialsback)=linmdl{NN,1}.Coefficients.Estimate(2:trialsback+1)';
  116. R_2R.RewHist.LinMdlPEPVal(NN,1:trialsback)=linmdl{NN,1}.Coefficients.pValue(2:trialsback+1)';
  117. %shuffled
  118. YSh=Y(randperm(length(Y)));
  119. linmdlSh{NN,1}=fitlm(Xpe,YSh);
  120. R_2R.RewHist.LinMdlPEWeightsSh(NN,1:trialsback)=linmdlSh{NN,1}.Coefficients.Estimate(2:trialsback+1)';
  121. R_2R.RewHist.LinMdlPEPValSh(NN,1:trialsback)=linmdlSh{NN,1}.Coefficients.pValue(2:trialsback+1)';
  122. %% stats comparing effect of current and previous reward on PE
  123. %resetting
  124. Bcell=0;
  125. EV1spk=0;
  126. EV2spk=0;
  127. EV1norm=0;
  128. EV2norm=0;
  129. %get mean baseline firing for all PE trials
  130. [Bcell1,B1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV1},PEBaseline,{2});% makes trial by trial rasters for baseline
  131. for y= 1:B1n
  132. Bcell(1,y)=sum(Bcell1{1,y}>PEBaseline(1));
  133. end
  134. [Bcell2,B2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV2},PEBaseline,{2});% makes trial by trial rasters for baseline
  135. for z= 1:B2n
  136. Bcell(1,z+B1n)=sum(Bcell2{1,z}>PEBaseline(1));
  137. end
  138. Bhz=Bcell/(PEBaseline(1,2)-PEBaseline(1,1));
  139. Bmean=nanmean(Bhz);
  140. Bstd=nanstd(Bhz);
  141. %get trial by trial firing rate for post suc trials
  142. [EV1cell,EV1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV1},CueWindow,{2});% makes trial by trial rasters for event
  143. for y= 1:EV1n
  144. EV1spk(1,y)=sum(EV1cell{1,y}>CueWindow(1));
  145. end
  146. EV1hz=EV1spk/(CueWindow(1,2)-CueWindow(1,1));
  147. for x= 1:EV1n
  148. EV1norm(1,x)=((EV1hz(1,x)-Bmean)/Bstd);
  149. end
  150. %get trial by trial firing rate for post mal trials
  151. [EV2cell,EV2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV2},CueWindow,{2});% makes trial by trial rasters for event
  152. for y= 1:EV2n
  153. EV2spk(1,y)=sum(EV2cell{1,y}>CueWindow(1));
  154. end
  155. EV2hz=EV2spk/(CueWindow(1,2)-CueWindow(1,1));
  156. for x= 1:EV2n
  157. EV2norm(1,x)=((EV2hz(1,x)-Bmean)/Bstd);
  158. end
  159. PEEvMeanz(NN,1)=nanmean(EV1norm);
  160. PEEvMeanz(NN,2)=nanmean(EV2norm);
  161. %% stats comparing effect of current and previous reward on cue
  162. %resetting
  163. Bcell=0;
  164. EV1spk=0;
  165. EV2spk=0;
  166. EV1norm=0;
  167. EV2norm=0;
  168. %get mean baseline firing for all cue trials
  169. [Bcell1,B1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV3},CueBaseline,{2});% makes trial by trial rasters for baseline
  170. for y= 1:B1n
  171. Bcell(1,y)=sum(Bcell1{1,y}>CueBaseline(1));
  172. end
  173. [Bcell2,B2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV4},CueBaseline,{2});% makes trial by trial rasters for baseline
  174. for z= 1:B2n
  175. Bcell(1,z+B1n)=sum(Bcell2{1,z}>CueBaseline(1));
  176. end
  177. Bhz=Bcell/(CueBaseline(1,2)-CueBaseline(1,1));
  178. Bmean=nanmean(Bhz);
  179. Bstd=nanstd(Bhz);
  180. %get trial by trial firing rate for post suc trials
  181. [EV1cell,EV1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV3},CueWindow,{2});% makes trial by trial rasters for event
  182. for y= 1:EV1n
  183. EV1spk(1,y)=sum(EV1cell{1,y}>CueWindow(1));
  184. end
  185. EV1hz=EV1spk/(CueWindow(1,2)-CueWindow(1,1));
  186. for x= 1:EV1n
  187. EV1norm(1,x)=((EV1hz(1,x)-Bmean)/Bstd);
  188. end
  189. %get trial by trial firing rate for post mal trials
  190. [EV2cell,EV2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV4},CueWindow,{2});% makes trial by trial rasters for event
  191. for y= 1:EV2n
  192. EV2spk(1,y)=sum(EV2cell{1,y}>CueWindow(1));
  193. end
  194. EV2hz=EV2spk/(CueWindow(1,2)-CueWindow(1,1));
  195. for x= 1:EV2n
  196. EV2norm(1,x)=((EV2hz(1,x)-Bmean)/Bstd);
  197. end
  198. CueEvMeanz(NN,1)=nanmean(EV1norm);
  199. CueEvMeanz(NN,2)=nanmean(EV2norm);
  200. fprintf('Neuron # %d\n',NN);
  201. end
  202. end
  203. R_2R.RewHist.PrevRewPEMeanz=PEEvMeanz;
  204. R_2R.RewHist.PrevRewCueMeanz=CueEvMeanz;
  205. end
  206. end
  207. %% which neurons to look at for stats and plotting?
  208. % Sel=R_2R.SucN | R_2R.MalN; %only look at reward-selective neurons
  209. Sel=NAneurons | VPneurons; %look at all neurons
  210. %Sel=R_2R.RewHist.LinMdlPVal(:,2)<0.1; %only neurons with significant impact of previous trial
  211. %Sel=R_2R.RewHist.LinMdlPEWeights(:,2)<-1; %only neurons with strong negative impact of previous trial
  212. %% ANOVAs
  213. %setup and run ANOVA for effects of current reward, previous reward, and region on reward firing
  214. PrevRew=cat(2,zeros(length(NAneurons),1),ones(length(NAneurons),1));
  215. Region=cat(2,NAneurons,NAneurons);
  216. Rat=cat(2,R_2R.Ninfo(:,4),R_2R.Ninfo(:,4));
  217. %to look at a specific selection of cells
  218. PEEvMeanz=R_2R.RewHist.PrevRewPEMeanz(Sel,:);
  219. CueEvMeanz=R_2R.RewHist.PrevRewCueMeanz(Sel,:);
  220. PrevRew=PrevRew(Sel,:);
  221. Region=Region(Sel,:);
  222. Rat=Rat(Sel,:);
  223. %cue
  224. %each region individually
  225. [~,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');
  226. [~,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');
  227. %region comparison
  228. [~,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');
  229. %pe
  230. %each region individually
  231. [~,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');
  232. [~,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');
  233. %region comparison
  234. [~,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');
  235. %setup and run ANOVA for effects of shuffle, trial, and region on coefficient
  236. Trial=[];
  237. Region=[];
  238. Rat=[];
  239. for i=1:trialsback
  240. Trial=cat(2,Trial,i*ones(length(NAneurons),1));
  241. Region=cat(2,Region,NAneurons);
  242. Rat=cat(2,Rat,R_2R.Ninfo(:,4));
  243. end
  244. Trial=cat(2,Trial,Trial);
  245. Region=cat(2,Region,Region);
  246. Rat=cat(2,Rat,Rat);
  247. Shuffd=cat(2,zeros(length(NAneurons),trialsback),ones(length(NAneurons),trialsback));
  248. %Cue
  249. Coeffs=cat(2,R_2R.RewHist.LinMdlCueWeights(:,1:trialsback),R_2R.RewHist.LinMdlCueWeightsSh(:,1:trialsback));
  250. [~,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');
  251. [~,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');
  252. %PE
  253. Coeffs=cat(2,R_2R.RewHist.LinMdlPEWeights(:,1:trialsback),R_2R.RewHist.LinMdlPEWeightsSh(:,1:trialsback));
  254. [~,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');
  255. [~,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');
  256. %% plotting
  257. XaxisCue=[-0.5 1];
  258. XaxisPE=[-1 2];
  259. Ishow=find(R_2R.Param.Tm>=XaxisCue(1) & R_2R.Param.Tm<=XaxisCue(2));
  260. Ushow=find(R_2R.Param.Tm>=XaxisPE(1) & R_2R.Param.Tm<=XaxisPE(2));
  261. time1=R_2R.Param.Tm(Ishow);
  262. time2=R_2R.Param.Tm(Ushow);
  263. %color map
  264. [magma,inferno,plasma,viridis]=colormaps;
  265. colormap(plasma);
  266. c=[-100 2000];ClimE=sign(c).*abs(c).^(1/4);%colormap
  267. %colors
  268. sucrose=[.95 0.55 0.15];
  269. maltodextrin=[.9 0.3 .9];
  270. water=[0.00 0.75 0.75];
  271. total=[0.3 0.1 0.8];
  272. inh=[0.1 0.021154 0.6];
  273. exc=[0.9 0.75 0.205816];
  274. NAc=[0.5 0.1 0.8];
  275. VP=[0.3 0.7 0.1];
  276. %extra colors to make a gradient
  277. sucrosem=[.98 0.8 0.35];
  278. sucrosel=[1 1 0.4];
  279. maltodextrinm=[1 0.75 1];
  280. maltodextrinl=[1 0.8 1];
  281. RD1P1=strcmp('CueP1', R_2R.Erefnames);
  282. RD1P2=strcmp('CueP2', R_2R.Erefnames);
  283. RD2P1=strcmp('PEP1', R_2R.Erefnames);
  284. RD2P2=strcmp('PEP2', R_2R.Erefnames);
  285. %% Get mean firing according to previous trial and then plot
  286. %NAc
  287. %plot suc after suc
  288. psth1=nanmean(R_2R.Ev(RD1P1).PSTHz(Sel&NAneurons,Ishow),1);
  289. sem1=nanste(R_2R.Ev(RD1P1).PSTHz(Sel&NAneurons,Ishow),1); %calculate standard error of the mean
  290. up1=psth1+sem1;
  291. down1=psth1-sem1;
  292. %plot suc after malt
  293. psth2=nanmean(R_2R.Ev(RD1P2).PSTHz(Sel&NAneurons,Ishow),1);
  294. sem2=nanste(R_2R.Ev(RD1P2).PSTHz(Sel&NAneurons,Ishow),1); %calculate standard error of the mean
  295. up2=psth2+sem2;
  296. down2=psth2-sem2;
  297. %make the plot
  298. subplot(2,4,1);
  299. hold on;
  300. title(['Cue response, NAc'])
  301. plot(time1,psth1,'Color',sucrose,'linewidth',1);
  302. plot(time1,psth2,'Color',maltodextrin,'linewidth',1);
  303. patch([time1,time1(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
  304. patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
  305. plot([-2 5],[0 0],':','color','k','linewidth',0.75);
  306. plot([0 0],[-2 8],':','color','k','linewidth',0.75);
  307. plot([CueWindow(1) CueWindow(1)],[-2 8],'color','b','linewidth',0.85);
  308. plot([CueWindow(2) CueWindow(2)],[-2 8],'color','b','linewidth',0.85);
  309. axis([XaxisCue(1) XaxisCue(2) min(down2)-0.15 max(up1)+0.2]);
  310. ylabel('Mean firing (z-score)');
  311. xlabel('Seconds from cue');
  312. legend('Post-suc','Post-mal','location','northeast');
  313. if cell2mat(R_2R.RewHist.PrevRewStatsCueNASubj{1,1}(2,7))<0.05
  314. plot(CueWindow(1)+(CueWindow(2)-CueWindow(1))/2,max(up2)+0.1,'*','color','k','markersize',13);
  315. end
  316. subplot(2,4,5);
  317. hold on;
  318. title('PE response, NAc');
  319. %plot malt after suc
  320. psth3=nanmean(R_2R.Ev(RD2P1).PSTHz(Sel&NAneurons,Ushow),1);
  321. sem3=nanste(R_2R.Ev(RD2P1).PSTHz(Sel&NAneurons,Ushow),1); %calculate standard error of the mean
  322. up3=psth3+sem3;
  323. down3=psth3-sem3;
  324. %plot malt after malt
  325. psth4=nanmean(R_2R.Ev(RD2P2).PSTHz(Sel&NAneurons,Ushow),1);
  326. sem4=nanste(R_2R.Ev(RD2P2).PSTHz(Sel&NAneurons,Ushow),1); %calculate standard error of the mean
  327. up4=psth4+sem4;
  328. down4=psth4-sem4;
  329. plot(time2,psth3,'Color',sucrose,'linewidth',1);
  330. plot(time2,psth4,'Color',maltodextrin,'linewidth',1);
  331. patch([time2,time2(end:-1:1)],[up3,down3(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
  332. patch([time2,time2(end:-1:1)],[up4,down4(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
  333. plot([-2 5],[0 0],':','color','k','linewidth',0.75);
  334. plot([0 0],[-2 8],':','color','k','linewidth',0.75);
  335. plot([PEWindow(1) PEWindow(1)],[-2 8],'color','b','linewidth',0.85);
  336. plot([PEWindow(2) PEWindow(2)],[-2 8],'color','b','linewidth',0.85);
  337. axis([XaxisPE(1) XaxisPE(2) min(down3)-0.2 max(up4)+0.2]);
  338. ylabel('Mean firing (z-score)');
  339. xlabel('Seconds from PE');
  340. legend('Post-suc','Post-mal','location','northeast');
  341. if cell2mat(R_2R.RewHist.PrevRewStatsPENASubj{1,1}(2,7))<0.05
  342. plot(PEWindow(1)+(PEWindow(2)-PEWindow(1))/2,max(up3)+0.1,'*','color','k','markersize',13);
  343. end
  344. %VP
  345. %plot suc after suc
  346. psth1=nanmean(R_2R.Ev(RD1P1).PSTHz(Sel&VPneurons,Ishow),1);
  347. sem1=nanste(R_2R.Ev(RD1P1).PSTHz(Sel&VPneurons,Ishow),1); %calculate standard error of the mean
  348. up1=psth1+sem1;
  349. down1=psth1-sem1;
  350. %plot suc after malt
  351. psth2=nanmean(R_2R.Ev(RD1P2).PSTHz(Sel&VPneurons,Ishow),1);
  352. sem2=nanste(R_2R.Ev(RD1P2).PSTHz(Sel&VPneurons,Ishow),1); %calculate standard error of the mean
  353. up2=psth2+sem2;
  354. down2=psth2-sem2;
  355. %make the plot
  356. subplot(2,4,2);
  357. title(['Cue response, VP'])
  358. hold on;
  359. plot(time1,psth1,'Color',sucrose,'linewidth',1);
  360. plot(time1,psth2,'Color',maltodextrin,'linewidth',1);
  361. patch([time1,time1(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
  362. patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
  363. plot([-2 5],[0 0],':','color','k','linewidth',0.75);
  364. plot([0 0],[-2 8],':','color','k','linewidth',0.75);
  365. plot([CueWindow(1) CueWindow(1)],[-2 8],'color','b','linewidth',0.85);
  366. plot([CueWindow(2) CueWindow(2)],[-2 8],'color','b','linewidth',0.85);
  367. axis([XaxisCue(1) XaxisCue(2) min(down2)-0.3 max(up1)+0.3]);
  368. ylabel('Mean firing (z-score)');
  369. xlabel('Seconds from cue');
  370. legend('Post-suc','Post-mal','location','northeast');
  371. if cell2mat(R_2R.RewHist.PrevRewStatsCueVPSubj{1,1}(2,7))<0.05
  372. plot(CueWindow(1)+(CueWindow(2)-CueWindow(1))/2,max(up2)+0.1,'*','color','k','markersize',13);
  373. end
  374. subplot(2,4,6);
  375. title('PE response, VP');
  376. hold on;
  377. %plot malt after suc
  378. psth3=nanmean(R_2R.Ev(RD2P1).PSTHz(Sel&VPneurons,Ushow),1);
  379. sem3=nanste(R_2R.Ev(RD2P1).PSTHz(Sel&VPneurons,Ushow),1); %calculate standard error of the mean
  380. up3=psth3+sem3;
  381. down3=psth3-sem3;
  382. %plot malt after malt
  383. psth4=nanmean(R_2R.Ev(RD2P2).PSTHz(Sel&VPneurons,Ushow),1);
  384. sem4=nanste(R_2R.Ev(RD2P2).PSTHz(Sel&VPneurons,Ushow),1); %calculate standard error of the mean
  385. up4=psth4+sem4;
  386. down4=psth4-sem4;
  387. plot(time2,psth3,'Color',sucrose,'linewidth',1);
  388. plot(time2,psth4,'Color',maltodextrin,'linewidth',1);
  389. patch([time2,time2(end:-1:1)],[up3,down3(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
  390. patch([time2,time2(end:-1:1)],[up4,down4(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
  391. plot([-2 5],[0 0],':','color','k','linewidth',0.75);
  392. plot([0 0],[-2 8],':','color','k','linewidth',0.75);
  393. plot([PEWindow(1) PEWindow(1)],[-2 8],'color','b','linewidth',0.85);
  394. plot([PEWindow(2) PEWindow(2)],[-2 8],'color','b','linewidth',0.85);
  395. axis([XaxisPE(1) XaxisPE(2) min(down4)-0.3 max(up3)+0.3]);
  396. ylabel('Mean firing (z-score)');
  397. xlabel('Seconds from PE');
  398. legend('Post-suc','Post-mal','location','northeast');
  399. if cell2mat(R_2R.RewHist.PrevRewStatsPEVPSubj{1,1}(2,7))<0.05
  400. plot(PEWindow(1)+(PEWindow(2)-PEWindow(1))/2,max(up3)+0.1,'*','color','k','markersize',13);
  401. end
  402. %% plot linear model coefficients
  403. %Plot all neurons
  404. Sel=NAneurons<2;
  405. %coefficients for each trial
  406. subplot(2,4,3);
  407. hold on;
  408. 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);
  409. 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);
  410. 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');
  411. 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');
  412. xlabel('Trials back');
  413. ylabel('Mean coefficient weight');
  414. title('Linear model coefficients');
  415. axis([0 trialsback+1 -0.5 1]);
  416. plot([-1 trialsback+1],[0 0],':','color','k','linewidth',0.75);
  417. legend('NAc','VP','Shuff');
  418. %stats to check if VP and NAc are greater than chance
  419. R_2R.RewHist.LinCoeffCueMultComp=[];
  420. [c,~,~,~]=multcompare(R_2R.RewHist.LinCoeffStatsCueNASubj{1,2},'dimension',[1,2],'display','off');
  421. [d,~,~,~]=multcompare(R_2R.RewHist.LinCoeffStatsCueVPSubj{1,2},'dimension',[1,2],'display','off');
  422. for i=1:trialsback
  423. %NAc vs shuff
  424. Cel=c(:,1)==2*(i-1)+1 & c(:,2)==2*(i-1)+2;
  425. if c(Cel,6)<0.05 R_2R.RewHist.LinCoeffCueMultComp(i,1)=1; else R_2R.RewHist.LinCoeffCueMultComp(i,1)=0; end
  426. R_2R.RewHist.LinCoeffCueMultComp(i,2)=c(Cel,2);
  427. %VP vs shuff
  428. Cel=d(:,1)==2*(i-1)+1 & d(:,2)==2*(i-1)+2;
  429. if d(Cel,6)<0.05 R_2R.RewHist.LinCoeffCueMultComp(i,3)=1; else R_2R.RewHist.LinCoeffCueMultComp(i,3)=0; end
  430. R_2R.RewHist.LinCoeffCueMultComp(i,4)=d(Cel,4);
  431. end
  432. plot([1:trialsback]-0.1,(R_2R.RewHist.LinCoeffCueMultComp(:,1)-0.5)*1.3,'*','color',NAc); %VP vs shuff
  433. plot([1:trialsback]+0.1,(R_2R.RewHist.LinCoeffCueMultComp(:,3)-0.5)*1.3,'*','color',VP); %NAc vs shuff
  434. %number of neurons with significant weights
  435. subplot(2,4,4);
  436. hold on;
  437. plot(1:trialsback,sum(R_2R.RewHist.LinMdlCuePVal(Sel&NAneurons,1:trialsback)<0.05,1)/sum(Sel&NAneurons),'color',NAc);
  438. plot(1:trialsback,sum(R_2R.RewHist.LinMdlCuePVal(Sel&VPneurons,1:trialsback)<0.05,1)/sum(Sel&VPneurons),'color',VP);
  439. plot(1:trialsback,sum(R_2R.RewHist.LinMdlCuePValSh(Sel&NAneurons,1:trialsback)<0.05,1)/sum(Sel&NAneurons),'color',NAc/3);
  440. plot(1:trialsback,sum(R_2R.RewHist.LinMdlCuePValSh(Sel&VPneurons,1:trialsback)<0.05,1)/sum(Sel&VPneurons),'color',VP/3);
  441. axis([0 trialsback+1 0 0.5]);
  442. xlabel('Trials back');
  443. ylabel('Fraction of the population');
  444. title('Outcome-modulated cue resp');
  445. legend('NAc','VP','Shuff');
  446. %Chi squared stat for each trial
  447. for i=1:trialsback
  448. [~,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));
  449. [~,R_2R.RewHist.LinMdlCue2region(i,1),R_2R.RewHist.LinMdlCue2region(i,2)]=crosstab(R_2R.RewHist.LinMdlCuePVal(Sel,i)<0.05,VPneurons);
  450. end
  451. %plot([1:trialsback]-0.1,(R_2R.RewHist.LinMdlCue2all(:,2)<0.05)-0.52,'*','color',NAc);
  452. plot([1:trialsback],(R_2R.RewHist.LinMdlCue2region(:,2)<0.05&R_2R.RewHist.LinMdlCue2all(:,2)<0.05)-0.52,'*','color','k');
  453. %% PE
  454. %coefficients for each trial
  455. subplot(2,4,7);
  456. hold on;
  457. 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);
  458. 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);
  459. 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');
  460. 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');
  461. xlabel('Trials back');
  462. ylabel('Mean coefficient weight');
  463. title('Linear model coefficients');
  464. axis([0 trialsback+1 -0.5 1]);
  465. plot([-1 trialsback+1],[0 0],':','color','k','linewidth',0.75);
  466. legend('NAc','VP','Shuff');
  467. %stats to check if VP and NAc are greater than chance
  468. R_2R.RewHist.LinCoeffPEMultComp=[];
  469. [c,~,~,~]=multcompare(R_2R.RewHist.LinCoeffStatsPENASubj{1,2},'dimension',[1,2],'display','off');
  470. [d,~,~,~]=multcompare(R_2R.RewHist.LinCoeffStatsPEVPSubj{1,2},'dimension',[1,2],'display','off');
  471. for i=1:trialsback
  472. %NAc vs shuff
  473. Cel=c(:,1)==2*(i-1)+1 & c(:,2)==2*(i-1)+2;
  474. if c(Cel,6)<0.05 R_2R.RewHist.LinCoeffPEMultComp(i,1)=1; else R_2R.RewHist.LinCoeffPEMultComp(i,1)=0; end
  475. R_2R.RewHist.LinCoeffPEMultComp(i,2)=c(Cel,2);
  476. %VP vs shuff
  477. Cel=d(:,1)==2*(i-1)+1 & d(:,2)==2*(i-1)+2;
  478. if d(Cel,6)<0.05 R_2R.RewHist.LinCoeffPEMultComp(i,3)=1; else R_2R.RewHist.LinCoeffPEMultComp(i,3)=0; end
  479. R_2R.RewHist.LinCoeffPEMultComp(i,4)=d(Cel,4);
  480. end
  481. plot([1:trialsback]-0.1,(R_2R.RewHist.LinCoeffPEMultComp(:,1)-0.5)*1.3,'*','color',NAc); %VP vs shuff
  482. plot([1:trialsback]+0.1,(R_2R.RewHist.LinCoeffPEMultComp(:,3)-0.5)*1.3,'*','color',VP); %NAc vs shuff
  483. %number of neurons with significant weights
  484. subplot(2,4,8);
  485. hold on;
  486. plot(1:trialsback,sum(R_2R.RewHist.LinMdlPEPVal(Sel&NAneurons,1:trialsback)<0.05,1)/sum(Sel&NAneurons),'color',NAc);
  487. plot(1:trialsback,sum(R_2R.RewHist.LinMdlPEPVal(Sel&VPneurons,1:trialsback)<0.05,1)/sum(Sel&VPneurons),'color',VP);
  488. plot(1:trialsback,sum(R_2R.RewHist.LinMdlPEPValSh(Sel&NAneurons,1:trialsback)<0.05,1)/sum(Sel&NAneurons),'color',NAc/3);
  489. plot(1:trialsback,sum(R_2R.RewHist.LinMdlPEPValSh(Sel&VPneurons,1:trialsback)<0.05,1)/sum(Sel&VPneurons),'color',VP/3);
  490. axis([0 trialsback+1 0 0.5]);
  491. xlabel('Trials back');
  492. ylabel('Fraction of the population');
  493. title('Outcome-modulated PE resp');
  494. legend('NAc','VP','Shuff');
  495. %Chi squared stat for each trial
  496. for i=1:trialsback
  497. [~,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));
  498. [~,R_2R.RewHist.LinMdlPEX2region(i,1),R_2R.RewHist.LinMdlPEX2region(i,2)]=crosstab(R_2R.RewHist.LinMdlPEPVal(Sel,i)<0.05,VPneurons);
  499. end
  500. %plot([0:trialsback]-0.1,(R_2R.RewHist.LinMdlPE2all(:,2)<0.05)-0.52,'*','color',NAc);
  501. plot([1:trialsback],(R_2R.RewHist.LinMdlPEX2region(:,2)<0.05&R_2R.RewHist.LinMdlPEX2all(:,2)<0.05)-0.52,'*','color','k');
  502. %% stats comparing PE coefficient weight in first 2 trials in selective and non-selective neurons in VP
  503. Sel=R_2R.SucN | R_2R.MalN;
  504. NSel=(R_2R.SucN | R_2R.MalN) == 0;
  505. [~,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)),...
  506. {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)),...
  507. 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');
  508. save('R_2R.mat','R_2R');