d_Fig2S2S3S5EventsBins.m 36 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866
  1. clear all;
  2. load ('R_2R.mat');
  3. %get parameters
  4. BinBase=R_2R.Param.BinBase;
  5. BinDura=R_2R.Param.BinDura;
  6. bins=R_2R.Param.bins;
  7. binint=R_2R.Param.binint;
  8. binstart=R_2R.Param.binstart;
  9. PStatBins=0.01; %using more stringent cutoff to reduce pre-delivery noise
  10. %divide neurons up by region
  11. NAneurons=strcmp(R_2R.Ninfo(:,3),'NA');
  12. VPneurons=strcmp(R_2R.Ninfo(:,3),'VP');
  13. %which bins bound the area examined for reward-selectivity? (in seconds)
  14. Time1=0.4; %seconds
  15. Time2=3; %seconds
  16. Bin1=(((Time1-BinDura(2)/2)-binstart)/binint); %convert to bin name
  17. Bin2=(((Time2-BinDura(2)/2)-binstart)/binint); %convert to bin name
  18. %sorting bin -- which bin the neurons' activity is sorted on for heatmap(in seconds)
  19. SortBinTime=1.1; %seconds
  20. SortBin=round((((SortBinTime-BinDura(2)/2)-binstart)/binint)); %convert to bin name
  21. sucrose=[1 0.6 0.1];
  22. maltodextrin=[.9 0.3 .9];
  23. water=[0.00 0.75 0.75];
  24. total=[0.3 0.1 0.8];
  25. exc=[0 113/255 188/255];
  26. inh=[240/255 0 50/255];
  27. NAc=[0.5 0.1 0.8];
  28. VP=[0.3 0.7 0.1];
  29. %% bins analysis
  30. %first and last bin aren't included because you can't compare to the previous/subsequent bin
  31. %this axis plots the bins on their centers
  32. xaxis=linspace(binstart+binint+BinDura(2)/2,binstart+(bins-2)*binint+BinDura(2)/2,bins-2);
  33. for i=2:(bins-1) %no including first or last bin because can't compare to previous/subsequent bin
  34. %finds out whether firing is stronger (high excitation or lower inhibition) for 1 or 2
  35. for k=1:length(R_2R.Ninfo) %runs through neurons
  36. if R_2R.BinStatRwrd{i,1}(k).IntSig < PStatBins %if neuron is significant for this bin
  37. if R_2R.BinStatRwrd{i-1,1}(k).IntSig < PStatBins || R_2R.BinStatRwrd{i+1,1}(k).IntSig < PStatBins %either previous or subsequent bin must be significant
  38. if R_2R.BinStatRwrd{i,1}(k).R1Mean > R_2R.BinStatRwrd{i,1}(k).R2Mean %if firing is greater on sucrose trials than maltodextrin
  39. R_2R.BinRewPref{i,1}(k,1)=1;
  40. else
  41. R_2R.BinRewPref{i,1}(k,1)=2; %otherwise maltodextrin must be greater than sucrose
  42. end
  43. else
  44. R_2R.BinRewPref{i,1}(k,1)=0; %if not significant in 2 consecutive bins
  45. end
  46. else
  47. R_2R.BinRewPref{i,1}(k,1)=0; %if no sig reward modulation
  48. end
  49. end
  50. %find how many NAc neurons have significant reward modulation in each bin
  51. NN1perBin(i,1)=sum(R_2R.BinRewPref{i,1}(NAneurons,1)==1); %sucrose pref
  52. NN2perBin(i,1)=sum(R_2R.BinRewPref{i,1}(NAneurons,1)==2); %malto pref
  53. NNperBin(i,1)=sum(R_2R.BinRewPref{i,1}(NAneurons,1)>0); %either
  54. %normalize to number of neurons in population
  55. NN1norm=NN1perBin./sum(NAneurons);
  56. NN2norm=NN2perBin./sum(NAneurons);
  57. NNnorm=NNperBin./sum(NAneurons);
  58. %find how many VP neurons have significant reward modulation in each bin
  59. NV1perBin(i,1)=sum(R_2R.BinRewPref{i,1}(VPneurons,1)==1); %sucrose pref
  60. NV2perBin(i,1)=sum(R_2R.BinRewPref{i,1}(VPneurons,1)==2); %malto pref
  61. NVperBin(i,1)=sum(R_2R.BinRewPref{i,1}(VPneurons,1)>0); %either
  62. %normalize to number of neurons in population
  63. NV1norm=NV1perBin./sum(VPneurons);
  64. NV2norm=NV2perBin./sum(VPneurons);
  65. NVnorm=NVperBin./sum(VPneurons);
  66. end
  67. %Cumulative reward selectivity
  68. cumsel=zeros(length(R_2R.Ninfo),bins); %set up result matrix
  69. cumsel1=zeros(length(R_2R.Ninfo),bins); %set up result matrix for sucrose
  70. cumsel2=zeros(length(R_2R.Ninfo),bins); %set up result matrix for maltodextrin
  71. %note, this has to be corrected to +1 because of the way the bin matrix is
  72. %layed out (bin 0 is the first column, not the 0th column)
  73. for i=1:length(R_2R.Ninfo) %go through each neuron
  74. for j=2:bins-1 %look in each bin we did analysis on
  75. if R_2R.BinRewPref{j,1}(i,1)>0 || cumsel(i,j-1)==1
  76. cumsel(i,j) = 1;
  77. end
  78. if R_2R.BinRewPref{j,1}(i,1)==1 || cumsel1(i,j-1)==1
  79. cumsel1(i,j) = 1;
  80. end
  81. if R_2R.BinRewPref{j,1}(i,1)==2 || cumsel2(i,j-1)==1
  82. cumsel2(i,j) = 1;
  83. end
  84. end
  85. end
  86. %plotting number of significantly modulated neurons across time
  87. %NAc
  88. subplot(2,3,1);
  89. hold on;
  90. plot(xaxis,NNnorm(2:bins-1),'Color', total,'linewidth',1.5);
  91. plot(xaxis,NN1norm(2:bins-1),'Color',sucrose,'linewidth',1.5);
  92. plot(xaxis,NN2norm(2:bins-1),'Color',maltodextrin,'linewidth',1.5);
  93. plot([Time1 Time1],[-1 1],':','color','k','linewidth',0.75);
  94. plot([Time2 Time2],[-1 1],':','color','k','linewidth',0.75);
  95. axis([xaxis(1) xaxis(end) 0 0.4]);
  96. legend('Total','Suc > mal','Mal > suc','Location','northeast')
  97. ylabel('Fraction of population');
  98. title('NAc reward-selective neurons');
  99. %VP
  100. subplot(2,3,2);
  101. hold on;
  102. plot(xaxis,NVnorm(2:bins-1),'Color', total,'linewidth',1.5);
  103. plot(xaxis,NV1norm(2:bins-1),'Color',sucrose,'linewidth',1.5);
  104. plot(xaxis,NV2norm(2:bins-1),'Color',maltodextrin,'linewidth',1.5);
  105. plot([Time1 Time1],[-1 1],':','color','k','linewidth',0.75);
  106. plot([Time2 Time2],[-1 1],':','color','k','linewidth',0.75);
  107. axis([xaxis(1) xaxis(end) 0 0.4]);
  108. legend('Total','Suc > mal','Mal > suc','Location','northeast')
  109. ylabel('Fraction of population');
  110. title('VP reward-selective neurons');
  111. %plotting cumulative selectivity in each region
  112. %NAc
  113. subplot(2,3,4);
  114. hold on;
  115. plot(xaxis,sum(cumsel(NAneurons,2:bins-1),1)/sum(NAneurons),'Color', total,'linewidth',1.5);
  116. plot(xaxis,sum(cumsel1(NAneurons,2:bins-1),1)/sum(NAneurons),'Color', sucrose,'linewidth',1.5);
  117. plot(xaxis,sum(cumsel2(NAneurons,2:bins-1),1)/sum(NAneurons),'Color', maltodextrin,'linewidth',1.5);
  118. plot([Time1 Time1],[0 1],':','color','k','linewidth',0.75);
  119. plot([Time2 Time2],[0 1],':','color','k','linewidth',0.75);
  120. axis([xaxis(1) xaxis(end) 0 0.75]);
  121. xlabel('Seconds from reward delivery');
  122. ylabel('Cumulative fraction');
  123. %VP
  124. subplot(2,3,5);
  125. hold on;
  126. plot(xaxis,sum(cumsel(VPneurons,2:bins-1),1)/sum(VPneurons),'Color', total,'linewidth',1.5);
  127. plot(xaxis,sum(cumsel1(VPneurons,2:bins-1),1)/sum(VPneurons),'Color', sucrose,'linewidth',1.5);
  128. plot(xaxis,sum(cumsel2(VPneurons,2:bins-1),1)/sum(VPneurons),'Color', maltodextrin,'linewidth',1.5);
  129. plot([Time1 Time1],[0 1],':','color','k','linewidth',0.75);
  130. plot([Time2 Time2],[0 1],':','color','k','linewidth',0.75);
  131. axis([xaxis(1) xaxis(end) 0 0.75]);
  132. xlabel('Seconds from RD');
  133. ylabel('Cumulative fraction');
  134. %plot difference between the two regions for each measure
  135. %reward selective neurons per bin
  136. subplot(2,3,3);
  137. hold on;
  138. plot(xaxis,NNnorm(2:bins-1)-NVnorm(2:bins-1),'Color', total,'linewidth',1.5);
  139. plot(xaxis,NN1norm(2:bins-1)-NV1norm(2:bins-1),'Color',sucrose,'linewidth',1.5);
  140. plot(xaxis,NN2norm(2:bins-1)-NV2norm(2:bins-1),'Color',maltodextrin,'linewidth',1.5);
  141. axis([xaxis(1) xaxis(end) -0.3 0.3]);
  142. plot([xaxis(1) xaxis(end)],[0 0],'color','k','linewidth',0.75);
  143. plot([Time1 Time1],[-1 1],':','color','k','linewidth',0.75);
  144. plot([Time2 Time2],[-1 1],':','color','k','linewidth',0.75);
  145. legend('Total','Suc > mal','Mal > suc','Location','northeast')
  146. title('Difference (NAc - VP)');
  147. ylabel('Fraction of population');
  148. %cumulative selectivity
  149. subplot(2,3,6);
  150. hold on;
  151. plot(xaxis,sum(cumsel(NAneurons,2:bins-1),1)/sum(cumsel(NAneurons,bins-1)),'Color', NAc,'linewidth',1.5);
  152. plot(xaxis,sum(cumsel(VPneurons,2:bins-1),1)/sum(cumsel(VPneurons,bins-1)),'Color', VP,'linewidth',1.5);
  153. plot([Time1 Time1],[-1 1],':','color','k','linewidth',0.75);
  154. plot([Time2 Time2],[-1 1],':','color','k','linewidth',0.75);
  155. axis([xaxis(1) xaxis(end) 0 1]);
  156. legend('NAc','VP','location','southeast');
  157. xlabel('Seconds from RD');
  158. ylabel('Cumulative distribution');
  159. title('Reward-selective onsets');
  160. %statistical test comparing onset of reward-selective responses
  161. for i=1:length(cumsel)
  162. if sum(cumsel(i,:))>0
  163. onsetcolumn=min(find(cumsel(i,:)));
  164. onset(i,1)=round((onsetcolumn-1)*binint+binstart+BinDura(2)/2,1);
  165. else
  166. onset(i,1)=NaN;
  167. end
  168. end
  169. [~,R_2R.OnsetSig{1,1},R_2R.OnsetSig{1,2}]=anovan(onset,{NAneurons,R_2R.Ninfo(:,4)},'varnames',{'Region','Subject'},'random',[2],'nested',[0 0;1 0],'display','off','model','full');
  170. if cell2mat(R_2R.OnsetSig{1,1}(2,7))<0.05
  171. plot((xaxis(end)-xaxis(1))/2+xaxis(1),0.95,'*','color','k','linewidth',1);
  172. end
  173. %% individual rat data
  174. figure;
  175. %split data up into individual rats
  176. for i=1:length(R_2R.Ninfo)
  177. A=char(R_2R.Ninfo(i,4));
  178. Rats(i,:)=cellstr(A);
  179. end
  180. C=unique(Rats);
  181. for i=1:length(C)
  182. Sel=strcmp(C(i),Rats);
  183. for j=2:bins-1
  184. NRat1(j,1)=sum(R_2R.BinRewPref{j,1}(Sel,1)==1)/sum(Sel);
  185. NRat2(j,1)=sum(R_2R.BinRewPref{j,1}(Sel,1)==2)/sum(Sel);
  186. NRat(j,1)=sum(R_2R.BinRewPref{j,1}(Sel,1)>0)/sum(Sel);
  187. end
  188. %graphing
  189. subplot(4,3,i);
  190. hold on;
  191. plot(xaxis,NRat(2:bins-1),'Color', total,'linewidth',1.5);
  192. plot(xaxis,NRat1(2:65),'Color',sucrose,'linewidth',1.5);
  193. plot(xaxis,NRat2(2:65),'Color',maltodextrin,'linewidth',1.5);
  194. xlabel('Seconds from reward delivery');
  195. title([cell2mat(C(i)) ' (n=' num2str(sum(Sel)) ')'])
  196. axis([xaxis(1) xaxis(end) 0 0.5]);
  197. legend('Total neurons','Suc > mal Hz','Mal > suc Hz','Location','northeast')
  198. ylabel('Fraction of total population');
  199. end
  200. %% plotting sucrose-selective neurons
  201. figure;
  202. %color map
  203. [magma,inferno,plasma,viridis]=colormaps;
  204. colormap(plasma);
  205. c=[-100 2000];ClimE=sign(c).*abs(c).^(1/4);%colormap
  206. %events we're looking at
  207. ev1=strcmp('RD1', R_2R.Erefnames);
  208. ev2=strcmp('RD2', R_2R.Erefnames);
  209. %setting up parameters
  210. Xaxis=[-2 5];
  211. inttime=find(R_2R.Param.Tm>=Xaxis(1) & R_2R.Param.Tm<=Xaxis(2));
  212. paramtime=R_2R.Param.Tm(inttime);
  213. %find all neurons with greater firing for sucrose
  214. for i = 1:(Bin2-Bin1+1)
  215. %the added +1 is necessary because bin 20 is the 21st entry in the matrix
  216. Pref1(:,i)=R_2R.BinRewPref{Bin1+i}==1; %get neurons that have greater firing for sucrose in any of the bins bounded above
  217. Resp11(:,i)=Pref1(:,i)&cat(1,R_2R.BinStatRwrd{Bin1+i,1}.SucRespDir)==1; %get neurons with excitation to sucrose
  218. Resp12(:,i)=Pref1(:,i)&cat(1,R_2R.BinStatRwrd{Bin1+i,1}.MalRespDir)==-1;%get neurons with inhibition to maltodextrin
  219. end
  220. Sel=sum(Pref1,2)>0; %all neurons selective for sucrose in any bin
  221. Sel2=sum(Resp11,2)>0; %all selective neurons sucrose excited in any bin
  222. Sel3=sum(Resp12,2)>0; %all selective neurons malto inhibited in any bin
  223. %-----------------NAc--------------------------------------------
  224. subplot(4,40,[15:23]); %heatmap of sucrose preferring neurons' response to sucrose
  225. Neurons=R_2R.Ev(ev1).PSTHz(Sel&NAneurons,inttime); %get the firing rates of neurons of interest
  226. SucResp=cat(1,R_2R.BinStatRwrd{SortBin+1,1}.R1Mean); %sucrose responses
  227. TMP=SucResp(Sel&NAneurons); %find the magnitude of the excitations for sucrose for this bin
  228. TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
  229. [~,SORTimg]=sort(TMP);
  230. Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
  231. imagesc(paramtime,[1,sum(Sel&NAneurons,1)],Neurons,ClimE);
  232. title(['Sucrose trials']);
  233. ylabel('Individual neurons');
  234. hold on;
  235. plot([0 0],[0 sum(Sel&NAneurons)],':','color','k','linewidth',0.75);
  236. subplot(4,40,[27:35]); %heatmap of sucrose preferring neurons' response to maltodextrin
  237. Neurons=R_2R.Ev(ev2).PSTHz(Sel&NAneurons,inttime); %get the firing rates of neurons of interest
  238. Neurons=Neurons(SORTimg,:); %sort the neurons same as before
  239. imagesc(paramtime,[1,sum(Sel&NAneurons,1)],Neurons,ClimE);
  240. title(['Maltodextrin trials']);
  241. hold on;
  242. plot([0 0],[0 sum(Sel&NAneurons)],':','color','k','linewidth',0.75);
  243. %plot sucrose preferring neurons to sucrose
  244. psth1=nanmean(R_2R.Ev(ev1).PSTHz(Sel&NAneurons,inttime),1);
  245. sem1=nanste(R_2R.Ev(ev1).PSTHz(Sel&NAneurons,inttime),1); %calculate standard error of the mean
  246. up1=psth1+sem1;
  247. down1=psth1-sem1;
  248. %plot sucrose preferring neurons to malt
  249. psth2=nanmean(R_2R.Ev(ev2).PSTHz(Sel&NAneurons,inttime),1);
  250. sem2=nanste(R_2R.Ev(ev2).PSTHz(Sel&NAneurons,inttime),1); %calculate standard error of the mean
  251. up2=psth2+sem2;
  252. down2=psth2-sem2;
  253. %plotting
  254. subplot(4,3,1);
  255. title(['Suc>mal (n=' num2str(sum(Sel&NAneurons)) ' of ' num2str(sum(NAneurons)) ')'])
  256. hold on;
  257. plot(paramtime,psth1,'Color',sucrose,'linewidth',1);
  258. plot(paramtime,psth2,'Color',maltodextrin,'linewidth',1);
  259. patch([paramtime,paramtime(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
  260. patch([paramtime,paramtime(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
  261. plot([-2 5],[0 0],':','color','k','linewidth',0.75);
  262. plot([0 0],[-1 10],':','color','k','linewidth',0.75);
  263. axis([-2 5 -1 4]);
  264. ylabel('Mean firing (z-score)');
  265. legend('Suc trials','Mal trials');
  266. %display inhibitions and excitations
  267. both = sum(Sel2&Sel3&NAneurons);
  268. excite = sum(Sel2&NAneurons)-both;
  269. inhib = sum(Sel3&NAneurons)-both;
  270. subplot(4,40,40);
  271. b = bar([inhib both excite; 0 0 0],'stacked');
  272. b(1,1).FaceColor=maltodextrin;
  273. b(1,2).FaceColor=total;
  274. b(1,3).FaceColor=sucrose;
  275. axis([1 1.2 0 both+excite+inhib]);
  276. ylabel('Inh / Both / Excited');
  277. set(gca,'xtick',[])
  278. %----------------VP----------------
  279. subplot(4,40,[80+15:80+23]); %heatmap of sucrose preferring neurons' response to sucrose
  280. Neurons=R_2R.Ev(ev1).PSTHz(Sel&VPneurons,inttime); %get the firing rates of neurons of interest
  281. SucResp=cat(1,R_2R.BinStatRwrd{SortBin+1,1}.R1Mean); %sucrose responses
  282. TMP=SucResp(Sel&VPneurons); %find the magnitude of the excitations for sucrose for this bin
  283. TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
  284. [~,SORTimg]=sort(TMP);
  285. Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
  286. imagesc(paramtime,[1,sum(Sel&VPneurons,1)],Neurons,ClimE);
  287. title(['Sucrose trials']);
  288. ylabel('Individual neurons');
  289. hold on;
  290. plot([0 0],[0 sum(Sel&VPneurons)],':','color','k','linewidth',0.75);
  291. subplot(4,40,[80+27:80+35]); %heatmap of sucrose preferring neurons' response to maltodextrin
  292. Neurons=R_2R.Ev(ev2).PSTHz(Sel&VPneurons,inttime); %get the firing rates of neurons of interest
  293. Neurons=Neurons(SORTimg,:); %sort the neurons same as before
  294. imagesc(paramtime,[1,sum(Sel&VPneurons,1)],Neurons,ClimE);
  295. title(['Maltodextrin trials']);
  296. hold on;
  297. plot([0 0],[0 sum(Sel&VPneurons)],':','color','k','linewidth',0.75);
  298. %plot sucrose preferring neurons to sucrose
  299. psth1=nanmean(R_2R.Ev(ev1).PSTHz(Sel&VPneurons,inttime),1);
  300. sem1=nanste(R_2R.Ev(ev1).PSTHz(Sel&VPneurons,inttime),1); %calculate standard error of the mean
  301. up1=psth1+sem1;
  302. down1=psth1-sem1;
  303. %plot sucrose preferring neurons to malt
  304. psth2=nanmean(R_2R.Ev(ev2).PSTHz(Sel&VPneurons,inttime),1);
  305. sem2=nanste(R_2R.Ev(ev2).PSTHz(Sel&VPneurons,inttime),1); %calculate standard error of the mean
  306. up2=psth2+sem2;
  307. down2=psth2-sem2;
  308. %plotting
  309. subplot(4,3,7);
  310. title(['Suc>mal (n=' num2str(sum(Sel&VPneurons)) ' of ' num2str(sum(VPneurons)) ')'])
  311. hold on;
  312. plot(paramtime,psth1,'Color',sucrose,'linewidth',1);
  313. plot(paramtime,psth2,'Color',maltodextrin,'linewidth',1);
  314. patch([paramtime,paramtime(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
  315. patch([paramtime,paramtime(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
  316. plot([-2 5],[0 0],':','color','k','linewidth',0.75);
  317. plot([0 0],[-1 10],':','color','k','linewidth',0.75);
  318. axis([-2 5 -0.6 4]);
  319. ylabel('Mean firing (z-score)');
  320. legend('Suc trials','Mal trials');
  321. %display inhibitions and excitations
  322. both = sum(Sel2&Sel3&VPneurons);
  323. excite = sum(Sel2&VPneurons)-both;
  324. inhib = sum(Sel3&VPneurons)-both;
  325. subplot(4,40,120);
  326. b = bar([inhib both excite; 0 0 0],'stacked');
  327. b(1,1).FaceColor=maltodextrin;
  328. b(1,2).FaceColor=total;
  329. b(1,3).FaceColor=sucrose;
  330. axis([1 1.2 0 both+excite+inhib]);
  331. ylabel('Inh / Both / Excited');
  332. set(gca,'xtick',[])
  333. %% plotting maltodextrin-selective neurons
  334. %find all neurons with greater firing for maltodextrin
  335. for i = 1:(Bin2-Bin1+1)
  336. %the added +1 is necessary because bin 20 is the 21st entry in the matrix
  337. Pref2(:,i)=R_2R.BinRewPref{Bin1+i}==2; %get neurons that have greater firing for maltodextrin in any of the bins bounded above
  338. Resp21(:,i)=Pref2(:,i)&cat(1,R_2R.BinStatRwrd{Bin1+i,1}.SucRespDir)==-1; %get neurons with inhibition to sucrose
  339. Resp22(:,i)=Pref2(:,i)&cat(1,R_2R.BinStatRwrd{Bin1+i,1}.MalRespDir)==1;%get neurons with excitation to maltodextrin
  340. end
  341. Mel=sum(Pref2,2)>0; %all neurons selective for malotdextrin in any bin
  342. Mel2=sum(Resp21,2)>0; %all selective neurons sucrose inhibited in any bin
  343. Mel3=sum(Resp22,2)>0; %all selective neurons malto excited in any bin
  344. %-----------------NAc--------------------------------------------
  345. subplot(4,40,[40+15:40+23]); %heatmap of maltodextrin preferring neurons' response to sucrose
  346. Neurons=R_2R.Ev(ev1).PSTHz(Mel&NAneurons,inttime); %get the firing rates of neurons of interest
  347. MalResp=cat(1,R_2R.BinStatRwrd{SortBin+1,1}.R2Mean); %sucrose responses
  348. TMP=MalResp(Mel&NAneurons); %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(paramtime,[1,sum(Mel&NAneurons,1)],Neurons,ClimE);
  353. title(['Sucrose trials']);
  354. xlabel('Seconds from RD');
  355. ylabel('Individual neurons');
  356. hold on;
  357. plot([0 0],[0 sum(Mel&NAneurons)],':','color','k','linewidth',0.75);
  358. subplot(4,40,[40+27:40+35]); %heatmap of maltodextrin preferring neurons' response to maltodextrin
  359. Neurons=R_2R.Ev(ev2).PSTHz(Mel&NAneurons,inttime); %get the firing rates of neurons of interest
  360. Neurons=Neurons(SORTimg,:); %sort the neurons same as before
  361. imagesc(paramtime,[1,sum(Mel&NAneurons,1)],Neurons,ClimE);
  362. title(['Maltodextrin trials']);
  363. xlabel('Seconds from RD');
  364. hold on;
  365. plot([0 0],[0 sum(Mel&NAneurons)],':','color','k','linewidth',0.75);
  366. %plot sucrose preferring neurons to sucrose
  367. psth1=nanmean(R_2R.Ev(ev1).PSTHz(Mel&NAneurons,inttime),1);
  368. sem1=nanste(R_2R.Ev(ev1).PSTHz(Mel&NAneurons,inttime),1); %calculate standard error of the mean
  369. up1=psth1+sem1;
  370. down1=psth1-sem1;
  371. %plot sucrose preferring neurons to malt
  372. psth2=nanmean(R_2R.Ev(ev2).PSTHz(Mel&NAneurons,inttime),1);
  373. sem2=nanste(R_2R.Ev(ev2).PSTHz(Mel&NAneurons,inttime),1); %calculate standard error of the mean
  374. up2=psth2+sem2;
  375. down2=psth2-sem2;
  376. %plotting
  377. subplot(4,3,4);
  378. title(['Mal>suc (n=' num2str(sum(Mel&NAneurons)) ' of ' num2str(sum(NAneurons)) ')'])
  379. hold on;
  380. plot(paramtime,psth1,'Color',sucrose,'linewidth',1);
  381. plot(paramtime,psth2,'Color',maltodextrin,'linewidth',1);
  382. patch([paramtime,paramtime(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
  383. patch([paramtime,paramtime(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
  384. plot([-2 5],[0 0],':','color','k','linewidth',0.75);
  385. plot([0 0],[-2 7],':','color','k','linewidth',0.75);
  386. axis([-2 5 -1.2 3]);
  387. ylabel('Mean firing (z-score)');
  388. xlabel('Seconds from RD');
  389. legend('Suc trials','Mal trials','Location','northeast');
  390. %display inhibitions and excitations
  391. both = sum(Mel2&Mel3&NAneurons);
  392. excite = sum(Mel3&NAneurons)-both;
  393. inhib = sum(Mel2&NAneurons)-both;
  394. subplot(4,40,80);
  395. b = bar([inhib both excite; 0 0 0],'stacked');
  396. b(1,1).FaceColor=sucrose;
  397. b(1,2).FaceColor=total;
  398. b(1,3).FaceColor=maltodextrin;
  399. axis([1 1.2 0 both+excite+inhib]);
  400. ylabel('Inh / Both / Excited');
  401. set(gca,'xtick',[]);
  402. %----------------VP----------------
  403. subplot(4,40,[120+15:120+23]); %heatmap of maltodextrin preferring neurons' response to sucrose
  404. Neurons=R_2R.Ev(ev1).PSTHz(Mel&VPneurons,inttime); %get the firing rates of neurons of interest
  405. MalResp=cat(1,R_2R.BinStatRwrd{SortBin+1,1}.R2Mean); %sucrose responses
  406. TMP=MalResp(Mel&VPneurons); %find the magnitude of the excitations for sucrose for this bin
  407. TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
  408. [~,SORTimg]=sort(TMP);
  409. Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
  410. imagesc(paramtime,[1,sum(Mel&VPneurons,1)],Neurons,ClimE);
  411. title(['Sucrose trials']);
  412. xlabel('Seconds from RD');
  413. ylabel('Individual neurons');
  414. hold on;
  415. plot([0 0],[0 sum(Mel&VPneurons)],':','color','k','linewidth',0.75);
  416. subplot(4,40,[120+27:120+35]); %heatmap of maltodextrin preferring neurons' response to maltodextrin
  417. Neurons=R_2R.Ev(ev2).PSTHz(Mel&VPneurons,inttime); %get the firing rates of neurons of interest
  418. Neurons=Neurons(SORTimg,:); %sort the neurons same as before
  419. imagesc(paramtime,[1,sum(Mel&VPneurons,1)],Neurons,ClimE);
  420. title(['Maltodextrin trials']);
  421. xlabel('Seconds from RD');
  422. hold on;
  423. plot([0 0],[0 sum(Mel&VPneurons)],':','color','k','linewidth',0.75);
  424. %plot maltodextrin preferring neurons to sucrose
  425. psth1=nanmean(R_2R.Ev(ev1).PSTHz(Mel&VPneurons,inttime),1);
  426. sem1=nanste(R_2R.Ev(ev1).PSTHz(Mel&VPneurons,inttime),1); %calculate standard error of the mean
  427. up1=psth1+sem1;
  428. down1=psth1-sem1;
  429. %plot malt preferring neurons to malt
  430. psth2=nanmean(R_2R.Ev(ev2).PSTHz(Mel&VPneurons,inttime),1);
  431. sem2=nanste(R_2R.Ev(ev2).PSTHz(Mel&VPneurons,inttime),1); %calculate standard error of the mean
  432. up2=psth2+sem2;
  433. down2=psth2-sem2;
  434. %plotting
  435. subplot(4,3,10);
  436. title(['Mal>suc (n=' num2str(sum(Mel&VPneurons)) ' of ' num2str(sum(VPneurons)) ')'])
  437. hold on;
  438. plot(paramtime,psth1,'Color',sucrose,'linewidth',1);
  439. plot(paramtime,psth2,'Color',maltodextrin,'linewidth',1);
  440. patch([paramtime,paramtime(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
  441. patch([paramtime,paramtime(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
  442. plot([-2 5],[0 0],':','color','k','linewidth',0.75);
  443. plot([0 0],[-2 7],':','color','k','linewidth',0.75);
  444. axis([-2 5 -1.5 2.5]);
  445. ylabel('Mean firing (z-score)');
  446. xlabel('Seconds from RD');
  447. legend('Suc trials','Mal trials','Location','northeast');
  448. %display inhibitions and excitations
  449. both = sum(Mel2&Mel3&VPneurons);
  450. excite = sum(Mel3&VPneurons)-both;
  451. inhib = sum(Mel2&VPneurons)-both;
  452. subplot(4,40,160);
  453. b = bar([inhib both excite; 0 0 0],'stacked');
  454. b(1,1).FaceColor=sucrose;
  455. b(1,2).FaceColor=total;
  456. b(1,3).FaceColor=maltodextrin;
  457. axis([1 1.2 0 both+excite+inhib]);
  458. ylabel('Inh / Both / Excited');
  459. set(gca,'xtick',[]);
  460. %% creating and saving new variables
  461. R_2R.RSinfo=R_2R.Ninfo(Sel|Mel,1:2); %reward-selective neurons
  462. R_2R.RNSinfo=R_2R.Ninfo(Sel==0 & Mel==0,1:2); %reward-nonselective neurons
  463. R_2R.SucN=Sel;
  464. R_2R.MalN=Mel;
  465. R_2R.Param.SelectiveBins=[Bin1 Bin2];
  466. % proportion summary and chi squared test for reward selective neurons
  467. R_2R.RegSelX2(1,1)=sum((Sel | Mel) & NAneurons)/sum(NAneurons); %how many selective NAc neurons
  468. R_2R.RegSelX2(1,2)=sum((Sel | Mel) & VPneurons)/sum(VPneurons); %how many selective VP neurons
  469. [~,R_2R.RegSelX2(2,1),R_2R.RegSelX2(2,2)]=crosstab(Sel | Mel,NAneurons);
  470. % proportion summary and chi squared test for max neurons in any bin
  471. R_2R.RegMaxSelX2(1,1)=max(NNnorm); %how many selective NAc neurons
  472. R_2R.RegMaxSelX2(1,2)=max(NVnorm); %how many selective VP neurons
  473. [~,R_2R.RegMaxSelX2(2,1),R_2R.RegMaxSelX2(2,2)]=crosstab(cat(1,ones(max(NNperBin),1),zeros(sum(NAneurons)-max(NNperBin),1),ones(max(NVperBin),1),zeros(sum(VPneurons)-max(NVperBin),1)),NAneurons);
  474. save('R_2R.mat','R_2R');
  475. %divide neurons up by region
  476. NAneurons=strcmp(R_2R.Ninfo(:,3),'NA');
  477. VPneurons=strcmp(R_2R.Ninfo(:,3),'VP');
  478. Xaxis=[-1 2];
  479. Ishow=find(R_2R.Param.Tm>=Xaxis(1) & R_2R.Param.Tm<=Xaxis(2));
  480. time1=R_2R.Param.Tm(Ishow);
  481. Xaxis2=[-2 2];
  482. Ushow=find(R_2R.Param.Tm>=Xaxis2(1) & R_2R.Param.Tm<=Xaxis2(2));
  483. time2=R_2R.Param.Tm(Ushow);
  484. inh=[0.1 0.021154 0.6];
  485. exc=[0.9 0.75 0.205816];
  486. Cue=strcmp('Cue', R_2R.Erefnames);
  487. PE=strcmp('PECue', R_2R.Erefnames);
  488. RD=strcmp('RD', R_2R.Erefnames); %should this be LickR or RD?
  489. %% Plotting heatmaps of all neurons to each event
  490. figure;
  491. %color map
  492. [magma,inferno,plasma,viridis]=colormaps;
  493. colormap(plasma);
  494. c=[-100 2000];ClimE=sign(c).*abs(c).^(1/4);%colormap
  495. for i=1:2
  496. if i==1 Reg=NAneurons; end %1 is NAc
  497. if i==2 Reg=VPneurons; end %2 is VP
  498. Sel=Reg&R_2R.Ev(Cue).RespDir<2; %gets all neurons
  499. %sort cue heatmap by magnitude of response
  500. Neurons=R_2R.Ev(Cue).PSTHz(Sel,Ishow); %get the firing rates of neurons of interest
  501. TMP=R_2R.Ev(Cue).Meanz(Sel); %find the magnitude of the inhibitions for this event
  502. TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
  503. [~,SORTimg]=sort(TMP);
  504. Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
  505. %heatmap
  506. subplot(4,6,[1 7]+(i-1)*12);
  507. imagesc(time1,[1,sum(Sel,1)],Neurons,ClimE); title('Cue responses');
  508. xlabel('Seconds post cue');
  509. ylabel('Individual neurons sorted by response strength');
  510. hold on;
  511. plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
  512. %sort PE heatmap by magnitude of response
  513. Neurons=R_2R.Ev(PE).PSTHz(Sel,Ushow); %get the firing rates of neurons of interest
  514. TMP=R_2R.Ev(PE).Meanz(Sel); %find the magnitude of the inhibitions for this event
  515. TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
  516. [~,SORTimg]=sort(TMP);
  517. Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
  518. %heatmap
  519. subplot(4,6,[3 9]+(i-1)*12);
  520. imagesc(time2,[1,sum(Sel,1)],Neurons,ClimE); title('PE responses');
  521. xlabel('Seconds post PE');
  522. hold on;
  523. plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
  524. %sort RD heatmap by magnitude of response
  525. Neurons=R_2R.Ev(RD).PSTHz(Sel,Ushow); %get the firing rates of neurons of interest
  526. TMP=R_2R.Ev(RD).Meanz(Sel); %find the magnitude of the inhibitions for this event
  527. TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
  528. [~,SORTimg]=sort(TMP);
  529. Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
  530. %heatmap
  531. subplot(4,6,[5 11]+(i-1)*12);
  532. imagesc(time2,[1,sum(Sel,1)],Neurons,ClimE); title('Reward responses');
  533. xlabel('Seconds post RD');
  534. hold on;
  535. plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
  536. %% Plotting neurons that respond to the cue
  537. %inhibitions
  538. Sel=Reg&R_2R.Ev(Cue).RespDir<0; %Find neurons that have an inhibitory response to this event
  539. %average firing rate
  540. subplot(4,6,2+(i-1)*12);
  541. psthI=nanmean(R_2R.Ev(Cue).PSTHz(Sel,Ishow),1); %find the average firing rate of the neurons at the time of the event
  542. semI=nanste(R_2R.Ev(Cue).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
  543. upI=psthI+semI;
  544. downI=psthI-semI;
  545. hold on;
  546. patch([time1,time1(end:-1:1)],[upI,downI(end:-1:1)],inh,'EdgeColor','none');alpha(0.5);
  547. plot(time1,psthI,'Color',inh,'linewidth',1); title('Mean firing (z-score)'); %create plot of avg firing rate
  548. plot([-1 2],[0 0],':','color','k','linewidth',0.75);
  549. plot([0 0],[-3 9],':','color','k','linewidth',0.75);
  550. axis([-1 2 min(downI)-0.2 max(upI)+0.2]);
  551. title(['Cue inhibitions (' num2str(round(100*sum(Sel)/sum(Reg))) '%)'])
  552. ylabel('Z-score');
  553. %excitation
  554. Sel=Reg&R_2R.Ev(Cue).RespDir>0; %Find neurons that have an excitatory response to this event
  555. %average firing rate
  556. subplot(4,6,8+(i-1)*12);
  557. psthE=nanmean(R_2R.Ev(Cue).PSTHz(Sel,Ishow),1);
  558. semE=nanste(R_2R.Ev(Cue).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
  559. upE=psthE+semE;
  560. downE=psthE-semE;
  561. hold on;
  562. patch([time1,time1(end:-1:1)],[upE,downE(end:-1:1)],exc,'EdgeColor','none');alpha(0.5);
  563. plot(time1,psthE,'Color',exc,'linewidth',1);
  564. plot([-1 2],[0 0],':','color','k','linewidth',0.75);
  565. plot([0 0],[-3 9],':','color','k','linewidth',0.75);
  566. axis([-1 2 min(downE)-0.2 max(upE)+0.3]);
  567. title(['Cue excitations (' num2str(round(100*sum(Sel)/sum(Reg))) '%)'])
  568. ylabel('Z-score');
  569. %% Plotting neurons that respond to PE
  570. %inhibitions
  571. Sel=Reg&R_2R.Ev(PE).RespDir<0; %Find neurons that have an inhibitory response to this event
  572. %average firing rate
  573. subplot(4,6,4+(i-1)*12);
  574. psthI=nanmean(R_2R.Ev(PE).PSTHz(Sel,Ushow),1); %find the average firing rate of the neurons at the time of the event
  575. semI=nanste(R_2R.Ev(PE).PSTHz(Sel,Ushow),1); %calculate standard error of the mean
  576. upI=psthI+semI;
  577. downI=psthI-semI;
  578. hold on;
  579. patch([time2,time2(end:-1:1)],[upI,downI(end:-1:1)],inh,'EdgeColor','none');alpha(0.5);
  580. plot(time2,psthI,'Color',inh,'linewidth',1); title('Mean firing (z-score)'); %create plot of avg firing rate
  581. plot([-2 2],[0 0],':','color','k','linewidth',0.75);
  582. plot([0 0],[-3 9],':','color','k','linewidth',0.75);
  583. axis([-2 2 min(downI)-0.2 max(upI)+0.2]);
  584. title(['PE inhibitions (' num2str(round(100*sum(Sel)/sum(Reg))) '%)'])
  585. ylabel('Z-score');
  586. %excitation
  587. Sel=Reg&R_2R.Ev(PE).RespDir>0; %Find neurons that have an excitatory response to this event
  588. %average firing rate
  589. subplot(4,6,10+(i-1)*12);
  590. psthE=nanmean(R_2R.Ev(PE).PSTHz(Sel,Ushow),1);
  591. semE=nanste(R_2R.Ev(PE).PSTHz(Sel,Ushow),1); %calculate standard error of the mean
  592. upE=psthE+semE;
  593. downE=psthE-semE;
  594. hold on;
  595. patch([time2,time2(end:-1:1)],[upE,downE(end:-1:1)],exc,'EdgeColor','none');alpha(0.5);
  596. plot(time2,psthE,'Color',exc,'linewidth',1);
  597. plot([-2 2],[0 0],':','color','k','linewidth',0.75);
  598. plot([0 0],[-3 9],':','color','k','linewidth',0.75);
  599. axis([-2 2 min(downE)-0.2 max(upE)+0.3]);
  600. title(['PE excitations (' num2str(round(100*sum(Sel)/sum(Reg))) '%)'])
  601. ylabel('Z-score');
  602. %% Plotting neurons that respond to reward
  603. %inhibitions
  604. Sel=Reg&R_2R.Ev(RD).RespDir<0; %Find neurons that have an inhibitory response to this event
  605. %average firing rate
  606. subplot(4,6,6+(i-1)*12);
  607. psthI=nanmean(R_2R.Ev(RD).PSTHz(Sel,Ushow),1); %find the average firing rate of the neurons at the time of the event
  608. semI=nanste(R_2R.Ev(RD).PSTHz(Sel,Ushow),1); %calculate standard error of the mean
  609. upI=psthI+semI;
  610. downI=psthI-semI;
  611. hold on;
  612. patch([time2,time2(end:-1:1)],[upI,downI(end:-1:1)],inh,'EdgeColor','none');alpha(0.5);
  613. plot(time2,psthI,'Color',inh,'linewidth',1); title('Mean firing (z-score)'); %create plot of avg firing rate
  614. plot([-1 2],[0 0],':','color','k','linewidth',0.75);
  615. plot([0 0],[-3 9],':','color','k','linewidth',0.75);
  616. axis([-1 2 min(downI)-0.2 max(upI)+0.2]);
  617. title(['RD inhibitions (' num2str(round(100*sum(Sel)/sum(Reg))) '%)'])
  618. ylabel('Z-score');
  619. %excitation
  620. Sel=Reg&R_2R.Ev(RD).RespDir>0; %Find neurons that have an excitatory response to this event
  621. %average firing rate
  622. subplot(4,6,12+(i-1)*12);
  623. psthE=nanmean(R_2R.Ev(RD).PSTHz(Sel,Ushow),1);
  624. semE=nanste(R_2R.Ev(RD).PSTHz(Sel,Ushow),1); %calculate standard error of the mean
  625. upE=psthE+semE;
  626. downE=psthE-semE;
  627. hold on;
  628. patch([time2,time2(end:-1:1)],[upE,downE(end:-1:1)],exc,'EdgeColor','none');alpha(0.5);
  629. plot(time2,psthE,'Color',exc,'linewidth',1);
  630. plot([-1 2],[0 0],':','color','k','linewidth',0.75);
  631. plot([0 0],[-3 9],':','color','k','linewidth',0.75);
  632. axis([-1 2 min(downE)-0.2 max(upE)+0.3]);
  633. title(['RD excitations (' num2str(round(100*sum(Sel)/sum(Reg))) '%)'])
  634. ylabel('Z-score');
  635. end
  636. %% plotting firing rate on sucrose and maltodextrin trials
  637. figure;
  638. %color map
  639. [magma,inferno,plasma,viridis]=colormaps;
  640. colormap(plasma);
  641. c=[-100 2000];ClimE=sign(c).*abs(c).^(1/4);%colormap
  642. Xaxis=[-1 4];
  643. Ishow=find(R_2R.Param.Tm>=Xaxis(1) & R_2R.Param.Tm<=Xaxis(2));
  644. time1=R_2R.Param.Tm(Ishow);
  645. Xaxis2=[-1 4];
  646. Ushow=find(R_2R.Param.Tm>=Xaxis2(1) & R_2R.Param.Tm<=Xaxis2(2));
  647. time2=R_2R.Param.Tm(Ushow);
  648. sucrose=[1 0.6 0.1];
  649. maltodextrin=[.9 0.3 .9];
  650. water=[0.00 0.75 0.75];
  651. total=[0.3 0.1 0.8];
  652. exc=[0 113/255 188/255];
  653. inh=[240/255 0 50/255];
  654. RD=strcmp('RD', R_2R.Erefnames);
  655. RD1=strcmp('RD1', R_2R.Erefnames);
  656. RD2=strcmp('RD2', R_2R.Erefnames);
  657. for i=1:2
  658. if i==1 Reg=NAneurons; end %1 is NAc
  659. if i==2 Reg=VPneurons; end %2 is VP
  660. %% Plotting heatmaps of sucrose and maltodextrin response
  661. Sel=Reg&R_2R.Ev(RD).RespDir<2; %gets all neurons
  662. %sort RD1 heatmap by magnitude of response
  663. Neurons=R_2R.Ev(RD1).PSTHz(Sel,Ushow); %get the firing rates of neurons of interest
  664. TMP=R_2R.Ev(RD).Meanz(Sel); %find the magnitude of the inhibitions for this event
  665. TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
  666. [~,SORTimg]=sort(TMP);
  667. Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
  668. %heatmap
  669. subplot(2,40,[1:9]+(i-1)*40);
  670. imagesc(time2,[1,sum(Sel,1)],Neurons,ClimE); title('Sucrose trials');
  671. ylabel('Individual neurons sorted by RD response');
  672. xlabel('Seconds from RD');
  673. hold on;
  674. plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
  675. %sort RD2 heatmap the same way
  676. Neurons=R_2R.Ev(RD2).PSTHz(Sel,Ushow); %get the firing rates of neurons of interest
  677. [~,SORTimg]=sort(TMP);
  678. Neurons=Neurons(SORTimg,:); %sort the neurons by magnitude
  679. xlabel('Seconds from RD');
  680. %heatmap
  681. subplot(2,40,[13:21]+(i-1)*40);
  682. imagesc(time2,[1,sum(Sel,1)],Neurons,ClimE); title('Maltodextrin trials');
  683. xlabel('Seconds from RD');
  684. hold on;
  685. plot([0 0],[0 sum(Sel)],':','color','k','linewidth',0.75);
  686. %% Plotting average firing rate of reward excitations
  687. Sel=Reg&(R_2R.Ev(RD).RespDir>0); %| R_2R.Ev(LickR1).RespDir>0 | R_2R.Ev(LickR2).RespDir>0); %Find neurons that have an inhibitory response to this event
  688. %average firing rate to sucrose
  689. psth1=nanmean(R_2R.Ev(RD1).PSTHz(Sel,Ishow),1);
  690. sem1=nanste(R_2R.Ev(RD1).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
  691. up1=psth1+sem1;
  692. down1=psth1-sem1;
  693. %average firing rate to maltodextrin
  694. psth2=nanmean(R_2R.Ev(RD2).PSTHz(Sel,Ishow),1);
  695. sem2=nanste(R_2R.Ev(RD2).PSTHz(Sel,Ishow),1); %calculate standard error of the mean
  696. up2=psth2+sem2;
  697. down2=psth2-sem2;
  698. %plotting
  699. subplot(4,5,[9 10]+(i-1)*10);
  700. hold on;
  701. plot(time1,psth1,'Color',sucrose,'linewidth',1);
  702. plot(time1,psth2,'Color',maltodextrin,'linewidth',1);
  703. patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
  704. patch([time1,time1(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
  705. plot([-1 4],[0 0],':','color','k','linewidth',0.75);
  706. plot([0 0],[-3 9],':','color','k','linewidth',0.75);
  707. axis([-1 4 min(down1)-0.2 max(up1)+0.3]);
  708. title(['Reward-excited (n=' num2str(sum(Sel)) ')'])
  709. ylabel('Mean firing (z-score)');
  710. xlabel('Seconds from RD');
  711. legend('Suc trials','Mal trials');
  712. % %onset
  713. % subplot(2,4,4);
  714. % cumul(R_2R.Ev(38).Onsets(Sel,1),[R_2R.Param.RespWin(1,1):0.04:R_2R.Param.RespWin(38,2)],2,maltodextrin);
  715. %% Plotting average firing rate of reward inhibitions
  716. Sel=Reg&(R_2R.Ev(RD).RespDir<0); %| R_2R.Ev(LickR1).RespDir<0 | R_2R.Ev(LickR2).RespDir<0); %Find neurons that have an inhibitory response to this event
  717. %avg firing rate to sucrose
  718. psth1=nanmean(R_2R.Ev(RD1).PSTHz(Sel,Ushow),1);
  719. sem1=nanste(R_2R.Ev(RD1).PSTHz(Sel,Ushow),1); %calculate standard error of the mean
  720. up1=psth1+sem1;
  721. down1=psth1-sem1;
  722. %avg firing rate to maltodextrin
  723. psth2=nanmean(R_2R.Ev(RD2).PSTHz(Sel,Ushow),1);
  724. sem2=nanste(R_2R.Ev(RD2).PSTHz(Sel,Ushow),1); %calculate standard error of the mean
  725. up2=psth2+sem2;
  726. down2=psth2-sem2;
  727. %plotting
  728. subplot(4,5,[4 5]+(i-1)*10);
  729. hold on;
  730. plot(time2,psth1,'Color',sucrose,'linewidth',1);
  731. plot(time2,psth2,'Color',maltodextrin,'linewidth',1);
  732. patch([time2,time2(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
  733. patch([time2,time2(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
  734. plot([-1 4],[0 0],':','color','k','linewidth',0.75);
  735. plot([0 0],[-3 9],':','color','k','linewidth',0.75);
  736. axis([-1 4 min(down1)-0.2 max(up1)+0.2]);
  737. title(['Reward-inhibited (n=' num2str(sum(Sel)) ')'])
  738. ylabel('Mean firing (z-score)');
  739. xlabel('Seconds from RD');
  740. legend('Suc trials','Mal trials','location','southeast');
  741. end