h_Fig4RewardHistory.m 22 KB

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