%ConteScriptHasselmo2008EJN.m - This script generates the figures presented in a 2008 %European Journal of Neuroscience paper. %CITATION: Hasselmo, M.E. (2008) Temporally structured replay of neural activity in a %model of entorhinal cortex, hippocampus and postsubiculum. Eur. J. Neurosci. 28:1301-1315. %Script written by %Michael E. Hasselmo, Boston University, copyright 2007, 2008, 2009, 2010. %This script does a single run as in the earlier figures of that paper (but can be reset %to run multiple runs of 10 as in later figures. Note that the connectivity from grid %cells to place cells is reset randomly, so the pattern is different on each run. %Activity during encoding (waking) is plotted in graphics on left. %Activity during Retrieval (REM) is plotted on the right. %Each row has different information: %1. Location and speed versus time are plotted in the top plots. %Location looks like a sine wave because it is one dimension of location in the circle. %Speed is a flat line with downward deflections for slowing at reward. %2. Place cell spiking activity vs time is in the middle plots (Encoding and Retrieval), %3. head direction cell activity vs. time are in the bottom plots, %4. head direction vectors vs. space are in the smaller square plots next to %head direction vs. time. %5. The plot on the upper right shows the perceived x,y position during encoding (gray line) %and during retrieval (black line). This shows where the rat "thinks" it is located %based on the inverse of the grid cell oscillation phase. %Can increase to 10 examples by setting NumPtestsTot=10 instead of 1. %Note that there is a "while" function that makes it continue running tests until %the internal retrieval (right) matches the number of cycles of running on the track during %encoding. %The pattern of grid cell activity vs space can be seen by setting plotGrids=1. %Note that there are 75 grid cells plotted during creation, encoding and retrieval, %so this slows the script down. %The pattern of place cell activity vs spcace can be seen by setting plotPlaceCells=1. %This plots an array of numPlotPlace rows and numPlotPlace columns for cells during encoding and %retrieval (not during creation). %Setting plotAllPlace=1 plots the sum of all 400 place cells at once. %Note that the plots of place cell spiking versus time plot the number of %cells determined by numPlaceToPlot=10. close all clear all tic %TIME fullT=1200; %START. FULLT. kplot=60; k=0; oldt=1; %Plotting kplace=1; kp=kplace; %Interval of place cell formation. Set kp to kplace so that first place will be same? %HEAD DIRECTION CELL PARAMETERS dt=0.02; h=[0; pi/3; 2*pi/3; pi; 4*pi/3; 5*pi/3]; HHst=[cos(h) sin(h)]; hd=zeros(length(h),fullT); %Head direction cell activity at each time. HDnum=length(h); %RUN REPEATED TESTS WITH DIFFERENT PARAMETERS numPsizesTot=1; %Number of different sizes of numP or other variable. Get memory errors when numP>36. numPtestsTot=1; % 3%Number of individual tests at a single variable value. numCircRun=zeros(numPsizesTot,numPtestsTot+1); numPlaceToPlot=10; %This sets number of place cells that will be plotted in place cell plots. %FIGURE FOR PLOTTING TRAJECTORIES figSize=numPsizesTot*numPtestsTot; figure('name','x y','Numbertitle','off','Position',[910 700 160 100+70*figSize]); ax14=gca; axes(ax14); rc=[1 0 0; 1 1 0; 1 0 1; 1 1 1; 0 1 1; 0 1 0.5; 0 1 0; 1 1 0; 1 0.5 0; 1 0 0]; colorSpikes=0; figure('name','x,y,speed vs time','Numbertitle','off','Position',[50 700 300 100+70]); ax111=gca; figure('name','x,y,speed vs time retrieval','Numbertitle','off','Position',[600 700 300 100+70]); ax112=gca; figure('name','Place Enc v time','Numbertitle','off','Position',[50 400 300 100+70*figSize]); ax16b=gca; figure('name','Place Ret v time','Numbertitle','off','Position',[600 400 300 100+70*figSize]); ax21b=gca; figure('name','HD encoding vs time','Numbertitle','off','Position',[50 100 300 100+70]); ax18b=gca; figure('name','HD retrieval vs time','Numbertitle','off','Position',[600 100 300 100+70*figSize]); ax19b=gca; figure('name','HD vectors Enc','Numbertitle','off','Position',[360 100 160 100+70]); ax18=gca; figure('name','HD vectors Ret','Numbertitle','off','Position',[910 100 160 100+70*figSize]); ax19=gca; strengHD=1.0; for numPsizes=1:numPsizesTot; numP=20+numPsizes*2; %Started with 4+numPsizes*2 to give 6,8,10,12,14,16.... %Get good function with numP=18+numPsizes*2 (e.g. starting point of numP=20!! for numPtests=1:numPtestsTot; while numCircRun(numPsizes,numPtests+1)0).*(HHst*xv(:,1)); end; cirangle=atan2(xv(2,t),xv(1,t)); dxv(1,t+1)=(1.0+0.0*rand)*sin(cirangle-1/96)/dt; %/dt; %1/4 larger inward vs. 1/ 64 larger outward. Use 1/32 to get larger inward movement? dxv(2,t+1)=-(1.0+0.0*rand)*cos(cirangle-1/96)/dt; %/dt; %For 20*sin -1/5 is best. for 10*sin xv(:,t+1) = xv(:,t) + dt*dxv(:,t+1) ; speed=sqrt(dxv(1,t+1)^2+dxv(2,t+1)^2); sumspeed=sumspeed+speed; end sumspeed*dt %Print total distance covered. fullT*dt %Print total time run. avespeed=sumspeed/fullT %Print average speed. pos_x=xv(1,:); pos_y=xv(2,:); elseif datainput==5; %Circle with stopping points. dxv=[0; -0.35]; reward=zeros(1,fullT); speed=zeros(1,fullT); dt=0.02; sumspeed=0; xmax=60; xmin=-60; ymax=60; ymin=-60; %Environment size for circular track. radius=47.5; stopangle=0; ccstop=1; for t=1:fullT if t==1 xv=[radius; 0]; phi(:,1)=(HHst*xv(:,1)>0).*(HHst*xv(:,1)); end; cirangle=atan2(xv(2,t),xv(1,t)); if t>1 cirangleOLD=atan2(xv(2,t-1),xv(1,t-1)); dangle=cirangleOLD-cirangle; if abs(dangle)2 && xv(1,t)>20 ccstop=ccstop+1; dxv(1,t+1)=(0.5+0.0*rand)*sin(cirangle-1/96)/dt; %This causes Slower movement dxv(2,t+1)=-(0.5+0.0*rand)*cos(cirangle-1/96)/dt; xv(:,t+1) = xv(:,t) + dt*dxv(:,t+1) ; %This updates with slow movement [ccstop t+1 xv(1,t) xv(2,t)] if ccstop==30; reward(1,t)=1; stopangle=0; ccstop=1; end else dxv(1,t+1)=(1.0+0.0*rand)*sin(cirangle-1/96)/dt; %/dt; %1/4 larger inward vs. 1/ 64 larger outward. Use 1/32 to get larger inward movement? dxv(2,t+1)=-(1.0+0.0*rand)*cos(cirangle-1/96)/dt; %/dt; %For 20*sin -1/5 is best. for 10*sin xv(:,t+1) = xv(:,t) + dt*dxv(:,t+1) ; end speed(1,t+1)=sqrt(dxv(1,t+1)^2+dxv(2,t+1)^2); sumspeed=sumspeed+speed(1,t+1); end sumspeed*dt %Print total distance covered. fullT*dt %Print total time run. avespeed=sumspeed/fullT %Print average speed. pos_x=xv(1,:); pos_y=xv(2,:); end axes(ax111); tt=[1:fullT+1]; plot(tt,pos_x,tt,speed); xlim([1 fullT]); ylim([-55 55]); %SOME GRAPHICS PARAMETERS plotH=1; %Sets whether head direction vectors are plotted at end. gy=0.7; %Sets trajectory color to gray. figsize=200; plotGrids=0; %Plot grid cell firing patterns (slow). plotPlaceCells=0; %Plot a subset of place cells (slow). numPlotPlace=5; %Number of rows and columns of place cells to plot when plotPlaceCells=1; plotAllPlace=1; %Plot the sum of all place cell firing. plotEnc=1; plotRet=0; %******INITIAL CELL CREATION************************ %CREATE GRID CELLS FOR INITIAL PLACE CELL CREATION phi=zeros(HDnum,fullT); %Vector for phase of dendrites (in radians). phis=zeros(HDnum/2,fullT); %Same soma phase in all parts of vector, but easier to compute with this. T=zeros(1,fullT); T(1,1)=1; iphi=[0; 0; 0]; %Vector of initial phases of dendrites (in radians) %iphi=[pi -0.5*pi -0.5*pi]; B=0.00385; %Scaling factor for effect of head direction on membrane oscillation frequency. %strengHD=0; %Strength of head direction input during retrieval. numF=3; %Number of different frequencies (usually 3 in most simulations) baseF=0; %Baseline frequency for range. (usually 0 in most simulations using cos(phi-phis) multFF=1; %Addition for going through frequencies. numPhi=5; %Number of different starting phases in each dimension (default=5 for most simulations) %numP=20; %Number of rows (and columns) of place cells. thresh=0.3; %Threshold for grid cells. pOff=0.2; %Grid cell phase offset. stdlim=2.7; %2.7 best for circle %Determines size of accepted place fields. xspike=zeros(2,numF,numPhi^2,fullT); %Location of grid cellspiking tspike=zeros(numF*numPhi^2,fullT); %Spiking of grid cells during creation. tspike2=zeros(1,numF*numPhi^2,fullT); %Spiking of grid cells during encoding. timep(1,numF*numPhi^2).spike=[]; allPlotPlaceX=zeros(1,1); %fullT allPlotPlaceY=zeros(1,1); sumPlace=zeros(1,fullT); if plotGrids==1 figure('name','diff grids','Numbertitle','off','Position',[250 50 numPhi*150 numPhi*150]); ax15=gca; end %CREATION OF GRID CELLS AND PLACE CELLS. useDXV=1; %Cycle through frequencies and phases for creating grid cells. for ff=1:numF; %GRID CELL INITIAL CREATION. f=baseF+multFF*ff; for pp=1:numPhi; for p2=1:numPhi; ix=[xv(1,1)+(pp+ff+pOff)*(sqrt(3)/(f*B*(numPhi))); xv(2,1)+(p2+ff+pOff)*(sqrt(3)/(f*B*(numPhi)))]; %f*G=H=sqrt(3)/2*B. G=sqrt(3)/2*f*B hx=((HHst*ix)>0).*(HHst*(ix)); LH=length(h)/2; MH=length(h); if useDXV==1; phis(:,1)=2*pi*f*dt; phi(:,1)=2*pi*f*(dt+B*hx); end; for t=2:fullT; if useDXV==1; phis(:,t)=phis(:,t-1)+2*pi*f*dt; phi(:,t)=phi(:,t-1)+2*pi*f*dt*(1+B*(HHst*dxv(:,t)>0).*(HHst*dxv(:,t))); elseif useDXV==0; phis(:,t)=2*pi*f*t; phi(:,t)=2*pi*f*(t+B*HHst*(xv(:,t)-ix)); end; %Creation of phase representation. end LH=length(h)/2; MH=length(h); g=prod(cos((phi(1:LH,:)-phi(LH+1:MH,:)))); spiket=g>thresh; %Spiking of grid cell. tspike((ff-1)*numPhi^2+(pp-1)*numPhi+p2,1:fullT)=spiket; %Store spiking of all grid cells. if plotGrids==1 %OPTIONAL GRID CELL FIGURE axes(ax15); axes(ax15); plotSpikeX=xv(1,1:fullT).*spiket; plotSpikeY=xv(2,1:fullT).*spiket; subplot(numPhi,numPhi,(p2-1)*numPhi+pp,'replace'); ax15=gca; plot(plotSpikeX,plotSpikeY,'r.'); xlim([xmin-1 xmax+1]) ylim([ymin-1 ymax+1]) set(gca,'XTickLabel',{' '}) set(gca,'YTickLabel',{' '}) drawnow end %End if plotGrids==1 toc end %End for p2=1:numPhi end %End for pp=1:numPhi end %End for f=1:numF %INITIAL PLACE CELLS CREATION - place cells. %numP=6; Number set with numP above. placeSpike=zeros(numP^2,fullT); WHP=zeros(HDnum,numP^2); %Creates matrix from place cells to HD cells. WPG=zeros(numP^2,numF*(numPhi^2)); %WPG matrix will store associations between grid cells and place cells. numGtoP=3; nP=1; nnnn=1; while nP<(numP^2+1) nnnn=nnnn+1; if nnnn>50000; break; end; gp=1; pPhi=ones(gp,1); for ppp=1:numF*numPhi^2; %SELECT THREE GRID CELLS AT RANDOM WITHOUT OVERLAP. pRand=rand>0.9; %Needs to be 0.9 or larger. 0.8 made it VERY slow - too many. if pRand==1; pPhi(gp,1)=ppp; gp=gp+1; end %Sequentially activates elements. if gp==numGtoP+1; break; end; %ppp=1:numF*numPhi^2; end; end pspike=prod(tspike(pPhi,1:fullT)); %Multiplies together random set of grid cells. plotPlaceX=xv(1,1:fullT).*pspike; plotPlaceY=xv(2,1:fullT).*pspike; if length(find(pspike))>3 && gp==numGtoP+1; %Only take place cells with at least three spikes from three grid inputs gp. stX=std(plotPlaceX(1,find(pspike))); stY=std(plotPlaceY(1,find(pspike))); if stX0.1 && stY>0.1; %Selects only place cells with three grid inputs and >3 spikes placeSpike(nP,:)=pspike; %Keep track of times of spikes of selected place cells. sumPlace=sumPlace+pspike; WPG(nP,pPhi)=1; %For a successful place cell, adds connectivity to WPG for input from grids. nP=nP+1; end %End std loop end %End find(pspike) loop end %End nP loop. nnnn toc if plotAllPlace==1; %Always plot the coverage of all place cells figure('name','all place cells','Numbertitle','off','Position',[1000 400 200 200]); ax17=gca; set(gcf,'DefaultAxesColorOrder',[gy gy gy]) axes(ax17); plot(xv(1,1:fullT).*(sumPlace>0.5),xv(2,1:fullT).*(sumPlace>0.5),'r.'); %Plots the sum of all place cell firing to see coverage. drawnow end %********ENCODING OF TRAJECTORY WITH PLACE, HD AND GRID************************ eeeT=fullT; %100; %fullT; tspike2=zeros(numF*numPhi^2,eeeT); placeSpike2=zeros(numP^2,eeeT); %Grid cell phase arrays ephis=zeros(length(h)/2,numF*numPhi^2); %,eeeT); ephi=zeros(length(h),numF*numPhi^2); %,eeeT); ephiOLD=zeros(length(h),numF*numPhi^2); ephisOLD=zeros(length(h)/2,numF*numPhi^2); useDXV=1; %Cycle through frequencies and phases for creating grid cells. for ff=1:numF; %GRID CELLS DURING ENCODING (USING SAME AS DURING CREATION. f=baseF+multFF*ff; for pp=1:numPhi; for p2=1:numPhi; ix=[xv(1,1)+(pp+ff+pOff)*(sqrt(3)/(f*B*(numPhi))); xv(2,1)+(p2+ff+pOff)*(sqrt(3)/(f*B*(numPhi)))]; %f*G=H=sqrt(3)/2*B. G=sqrt(3)/2*f*B hx=((HHst*ix)>0).*(HHst*(ix)); LH=length(h)/2; MH=length(h); if useDXV==1; ephis(:,1)=2*pi*f*dt; ephi(:,1)=2*pi*f*(dt+B*hx); end; for t=2:eeeT; if useDXV==1; ephis(:,t)=ephis(:,t-1)+2*pi*f*dt; ephi(:,t)=ephi(:,t-1)+2*pi*f*dt*(1+B*(HHst*dxv(:,t)>0).*(HHst*dxv(:,t))); elseif useDXV==0; ephis(:,t)=2*pi*f*t; ephi(:,t)=2*pi*f*(t+B*HHst*(xv(:,t)-ix)); end; %Creation of phase representation. end LH=length(h)/2; MH=length(h); g=prod(cos((ephi(1:LH,:)-ephi(LH+1:MH,:)))); spiket=g>thresh; %Spiking of grid cell. tspike2((ff-1)*numPhi^2+(pp-1)*numPhi+p2,1:fullT)=spiket; %Store spiking of all grid cells. end %End for p2=1:numPhi end %End for pp=1:numPhi end %End for f=1:numF for t=1:eeeT; %PLACE CELLS DURING ENCODING placeSpike2(:,t)=((WPG*tspike2(:,t))./numGtoP)>0.9; %Spiking of place cells during encoding. speed(1,t)=sqrt(dxv(1,t)^2+dxv(2,t)^2); %HD CELLS DURING ENCODING hd(:,t)=(HHst*dxv(:,t)>0).*(HHst*dxv(:,t)); %Add the head direction cell activity during encoding. hd(:,t); plothd(:,t)=hd(1:LH,t)-hd(LH+1:MH,t); norm=sum(placeSpike2(:,t)); if norm~=0; for xx=1:numP^2; if placeSpike2(xx,t)==1; WHP(:,xx)=(WHP(:,xx)+hd(:,t).*placeSpike2(xx,t))/2; %Create connections from place to head direction. %WHP(:,xx)=(WHP(:,xx)+hd(:,t).*placeSpike2(xx,t)./norm)/2; %Create connections from place to head direction. end %This procedure takes elements with presynaptic input and averages the existing and new weight. end end end %End for t=1:eeeT if plotGrids==1 figure('name','grid during encoding','Numbertitle','off','Position',[250 50 numPhi*150 numPhi*150]); ax15b=gca; axes(ax15b); %Grid cell figure %Plot Grid cells during encoding. for ff=1:numF; f=baseF+multFF*ff; for pp=1:numPhi; for p2=1:numPhi; plotSpikeX=xv(1,1:eeeT).*tspike2((ff-1)*numPhi^2+(pp-1)*numPhi+p2,1:eeeT); plotSpikeY=xv(2,1:eeeT).*tspike2((ff-1)*numPhi^2+(pp-1)*numPhi+p2,1:eeeT); axes(ax15b); %subplot(numF*numPhi,numPhi,(f-1)*numPhi^2+(pp-1)*numPhi+p2); ax15b=gca; subplot(numPhi,numPhi,(pp-1)*numPhi+p2,'replace'); ax15b=gca; plot(plotSpikeX,plotSpikeY,'r.'); xlim([xmin-1 xmax+1]) ylim([ymin-1 ymax+1]) set(gca,'XTickLabel',{' '}) set(gca,'YTickLabel',{' '}) drawnow end %End for p2=1:numPhi end %End for pp=1:numPhi end %End for f=1:numF end %End if plotGrids==1 if plotPlaceCells==1 figure('name','place cells Encoding','Numbertitle','off','Position',[800 50 numPlotPlace*70 numPlotPlace*70]); ax16bbb=gca; %Plot Place Cells during encoding. for pp=1:numPlotPlace^2; plotPlaceX=xv(1,1:eeeT).*placeSpike2(pp,:); plotPlaceY=xv(2,1:eeeT).*placeSpike2(pp,:); axes(ax16bbb); subplot(numPlotPlace,numPlotPlace,pp,'replace'); ax16bbb=gca; plot(plotPlaceX,plotPlaceY,'r.'); xlim([xmin-1 xmax+1]) ylim([ymin-1 ymax+1]) set(gca,'XTickLabel',{' '}) set(gca,'YTickLabel',{' '}) end %End for pp end %End if plotPlaceCells==1 %PLOT PLACE CELLS VS TIME DURING ENCODING spike times with colors. axes(ax16b); subplot(numPsizesTot*numPtestsTot,1,(numPsizes-1)*numPtestsTot+numPtests,'replace'); ax16b=gca; numSING=numPlaceToPlot; for sing=1:numSING; fps=find(placeSpike2(sing,:)); for mm=1:length(fps); if colorSpikes==1; plot([fps(mm) fps(mm)],[sing sing-0.5],'Color',[rc(sing,:)],'LineWidth',1); hold on; elseif colorSpikes==0; plot([fps(mm) fps(mm)],[sing sing-0.5],'Color','k','LineWidth',1); hold on; end end end if colorSpikes==1; set(gca,'Color','k'); elseif colorSpikes==0 set(gca,'Color','w'); end xlim([1 eeeT]); ylim([0 numSING+0.5]); drawnow if plotEnc==1; %PLOT HD cells during encoding axes(ax18); for m=1:5:eeeT-1; %Plot the head direction cell vectors. %If want to plot those during encoding start m=1: denom=cos(h(1,1))*sin(h(2,1))-sin(h(1,1))*cos(h(2,1)); hdx=((plothd(1,m+1))*sin(h(2,1))-(plothd(2,m+1))*sin(h(1,1))); %./(denom*wd*B); hdy=((plothd(1,m+1))*-cos(h(2,1))+(plothd(2,m+1))*cos(h(1,1))); %./(denom*wd*B); AA=[xv(1,m) xv(1,m)+hdx]; BB=[xv(2,m) xv(2,m)+hdy]; CC=[xv(1,m) xv(1,m)+dxv(1,m+1)*dt]; DD=[xv(2,m) xv(2,m)+dxv(2,m+1)*dt]; if sum(AA)~=0 && sum(BB)~=0; plot(AA,BB,'k-',CC,DD,'r-','LineWidth',1); %AA,BB,'k-', hold on end end %PLOT HD CELLS VERSUS TIME DURING ENCODING axes(ax18b) spikehd=zeros(length(h),length(hd)); for ht=1:length(hd); maxhd=mean(max(hd)); spikehd(:,ht)=rand(length(h),1)>(maxhd*ones(length(h),1)-hd(:,ht))/maxhd; end numSING=length(h); for sing=1:numSING; fps=find(spikehd(sing,:)); for mm=1:length(fps); %plot([fps(mm) fps(mm)],[sing sing-0.5],'Color',[rc(sing,:)],'LineWidth',1); hold on; if colorSpikes==1; plot([fps(mm) fps(mm)],[sing sing-0.5],'Color',[rc(sing,:)],'LineWidth',1); hold on; elseif colorSpikes==0; plot([fps(mm) fps(mm)],[sing sing-0.5],'Color','k','LineWidth',1); hold on; end end end if colorSpikes==1; set(gca,'Color','k'); elseif colorSpikes==0 set(gca,'Color','w'); end xlim([1 eeeT]); ylim([0 numSING+0.5]); drawnow end %End plotEnc==1 toc %********RETRIEVAL OF TRAJECTORY WITH PLACE, HD AND GRID************************ %Add RETRIEVAL period. rrrT=eeeT; %eeeT; hd3=zeros(length(h),rrrT); hd3(:,1:5)=hd(:,5:9); %Set first step of retrieval hd to same as encoding hd. tspike3=zeros(numF*numPhi^2,rrrT); tspike3(:,1:5)=tspike2(:,5:9); %Set initial grid cell spiking to be same as during encoding %Grid cell spiking initial conditions don't matter. Can set to tspike2(:,21:22) and still function. But does initial phase set by ix matter? placeSpike3=zeros(numP^2,rrrT); placeSpike3(:,1:5)=placeSpike2(:,5:9); %set initial place cell to be same as during encoding. rphis=zeros(length(h)/2,numF*numPhi^2); rphi=zeros(length(h),numF*numPhi^2); plotrphi=zeros(length(h)/2,rrrT); plotiG=zeros(length(h),1); rphiOLD=zeros(length(h),numF*numPhi^2); rphisOLD=zeros(length(h)/2,numF*numPhi^2); if useDXV==1; %Set initial conditions for grid cells. for ff=1:numF; f=baseF+multFF*ff; for pp=1:numPhi; for p2=1:numPhi; cG=(ff-1)*numPhi^2+(pp-1)*numPhi+p2; %Number of current grid cell. ix=[xv(1,1)+(pp+ff+pOff)*ceil(sqrt(3)/(f*B*(numPhi))); xv(2,1)+(p2+ff+pOff)*ceil(sqrt(3)/(f*B*(numPhi)))]; %f*G=H=sqrt(3)/2*B. G=sqrt(3)/2*f*B iG=((HHst*ix)>0).*(HHst*ix); for nG=1:length(h); rphis(nG,cG)=2*pi*f*dt; rphi(nG,cG)=2*pi*f*(dt+B*iG(nG)); end if cG==1 plotiG=(HHst*(ix-xv(:,1))>0).*(HHst*(ix-xv(:,1))); end end %End for p2=1:numPhi end %End for pp=1:numPhi end %End for f=1:numF end %End if useDXV. %RETRIEVAL ACTIVITY OF GRID CELLS, PLACE CELLS AND HD CELLS for t=2:rrrT; %GRID CELLS for ff=1:numF; f=baseF+multFF*ff; for pp=1:numPhi; for p2=1:numPhi; cG=(ff-1)*numPhi^2+(pp-1)*numPhi+p2; %Number of current grid cell. ix=[xv(1,1)+(pp+ff+pOff)*ceil(sqrt(3)/(f*B*(numPhi))); xv(2,1)+(p2+ff+pOff)*ceil(sqrt(3)/(f*B*(numPhi)))]; %f*G=H=sqrt(3)/2*B. G=sqrt(3)/2*f*B if useDXV==1; rphiOLD=rphi; rphisOLD=rphis; rphis(:,cG)=rphisOLD(:,cG)+2*pi*f*dt; rphi(:,cG)=rphiOLD(:,cG)+2*pi*f*dt*(1+B*strengHD*hd3(:,t)); %Use head direction vector hd3 retrieved from WHP. elseif useDXV==0; rphis(:,cG)=2*pi*f*t; rphi(:,cG)=2*pi*f*(t+B*hd3(:,t)); end LH=length(h)/2; MH=length(h); tspike3(cG,t)=prod(cos(rphi(1:LH,cG)-rphi(LH+1:MH,cG)))>thresh; if cG==1; plotrphi(:,t)=rphi(1:LH,cG)-rphi(LH+1:MH,cG); plothd3(:,t)=hd3(1:LH,t)-hd3(LH+1:MH,t); end end %End for p2=1:numPhi end %End for pp=1:numPhi end %End for f=1:numF %PLACE CELLS DURING RETRIEVAL placeSpike3(:,t)=((WPG*tspike3(:,t))./numGtoP)>0.9; %Create place cell activity during encoding. %HD CELLS if max(placeSpike3(:,t))==0 hd3(:,t+1)=hd3(:,t); else norm=sum(placeSpike3(:,t)); hd3(:,t+1)=WHP*placeSpike3(:,t)/norm; end end %End for t=1:rrrT %PLOT INVERSE TRANSFORM TO SEE WHERE RAT "THINKS" IT IS LOCATED. axes(ax14); %Go back to x y plot and plot inverse. %subplot(numPsizesTot, numPtestsTot, (numPsizes-1)*numPtestsTot+numPtests); ax14=gca; %This does rows of tests. subplot(numPsizesTot*numPtestsTot,1,(numPsizes-1)*numPtestsTot+numPtests,'replace'); ax14=gca; %This does one column of all. %PLOT ACTUAL TRAJECTORY IN GRAY set(gcf,'DefaultAxesColorOrder',[gy gy gy]) plot(pos_x(minT:minT+fullT), pos_y(minT:minT+fullT)); hold on %PLOT RETRIEVED TRAJECTORY IN BLACK set(gcf,'DefaultAxesColorOrder',[0 0 0]) wd=(baseF+multFF*1)*2*pi; %Calculate inverse transform of grid cell phase (plotrphi(1,:) and plotrphi(2,:) denom=cos(h(1,1))*sin(h(2,1))-sin(h(1,1))*cos(h(2,1)); xph=((plotrphi(1,:)-B*wd*plotiG(1,1))*sin(h(2,1))-(plotrphi(2,:)-B*wd*plotiG(2,1))*sin(h(1,1)))./(denom*wd*B); yph=((plotrphi(1,:)-B*wd*plotiG(1,1))*-cos(h(2,1))+(plotrphi(2,:)-B*wd*plotiG(2,1))*cos(h(1,1)))./(denom*wd*B); plot(xph(1,2:rrrT),yph(1,2:rrrT),'k','LineWidth',2); %xlim([xmin+30 xmax+30]); ylim([ymin-60 ymax-60]); xlim([xmin-1 xmax+1]) ylim([ymin-1 ymax+1]) drawnow toc %PLOT GRID CELLS DURING RETRIEVAL. if plotGrids==1 figure('name','grid Retrieval','Numbertitle','off','Position',[250 50 numPhi*150 numF*150]); ax20=gca; for ff=1:numF; f=baseF+multFF*ff; for pp=1:numPhi; for p2=1:numPhi; plotSpikeX=xph(1,1:rrrT).*tspike3((ff-1)*numPhi^2+(pp-1)*numPhi+p2,1:rrrT); plotSpikeY=yph(1,1:rrrT).*tspike3((ff-1)*numPhi^2+(pp-1)*numPhi+p2,1:rrrT); axes(ax20); subplot(numF,numPhi,(ff-1)*numPhi+pp,'replace'); ax20=gca; plot(plotSpikeX,plotSpikeY,'r.'); xlim([xmin-1 xmax+1]) ylim([ymin-1 ymax+1]) drawnow end %End for p2=1:numPhi end %End for pp=1:numPhi end %End for f=1:numF end %End if plotGrids==1 %Plot Place Cells during retrieval. if plotPlaceCells==1 figure('name','place cells Retrieval','Numbertitle','off','Position',[800 50 numPlotPlace*70 numPlotPlace*70]); ax21=gca; for pp=1:numPlotPlace^2; plotPlaceX=xph(1,1:rrrT).*placeSpike3(pp,:); plotPlaceY=yph(1,1:rrrT).*placeSpike3(pp,:); axes(ax21); subplot(numPlotPlace,numPlotPlace,pp,'replace'); ax21=gca; plot(plotPlaceX,plotPlaceY,'r.'); xlim([xmin-1 xmax+1]) ylim([ymin-1 ymax+1]) set(gca,'XTickLabel',{' '}) set(gca,'YTickLabel',{' '}) end %End for pp end %End if plotPlaceCells==1 %PLOT PLACE VS TIME DURING RETRIEVAL. axes(ax21b); subplot(numPsizesTot*numPtestsTot,1,(numPsizes-1)*numPtestsTot+numPtests,'replace'); ax21b=gca; numSING=numPlaceToPlot; for sing=1:numSING; %for mm=1:length(timep(sing).spike); fps=find(placeSpike3(sing,:)); for mm=1:length(fps); if colorSpikes==1; plot([fps(mm) fps(mm)],[sing sing-0.5],'Color',[rc(sing,:)],'LineWidth',1); hold on; elseif colorSpikes==0; plot([fps(mm) fps(mm)],[sing sing-0.5],'Color','k','LineWidth',1); hold on; end %plot([fps(mm) fps(mm)],[sing sing-0.5],'Color',[rc(sing,:)],'LineWidth',1); hold on; %plot([fps(mm) fps(mm)],[sing sing-0.5],'Color','k','LineWidth',1); hold on; end end if colorSpikes==1; set(gca,'Color','k'); elseif colorSpikes==0 set(gca,'Color','w'); end xlim([1 rrrT]); ylim([0 numSING+0.5]); drawnow %PLOT HD CELLS VS SPACE during retrieval %figure('name','HD cells Retrieval','Numbertitle','off','Position',[750 400 300 300]); ax19=gca; axes(ax19); subplot(numPsizesTot*numPtestsTot,1,(numPsizes-1)*numPtestsTot+numPtests,'replace'); ax19=gca; %This does one column of all. plot(xph(1,2:rrrT),yph(1,2:rrrT),'r-','LineWidth',1); %xlim([xmin+30 xmax+30]); ylim([ymin-60 ymax-60]); hold on for m=2:5:rrrT-1; %Plot the head direction cell vectors. %If want to plot those during encoding start m=1: denom=cos(h(1,1))*sin(h(2,1))-sin(h(1,1))*cos(h(2,1)); hdx=((plothd3(1,m+1))*sin(h(2,1))-(plothd3(2,m+1))*sin(h(1,1))); %./(denom*wd*B); hdy=((plothd3(1,m+1))*-cos(h(2,1))+(plothd3(2,m+1))*cos(h(1,1))); %./(denom*wd*B); %hdx=hd3(1,m+1); %hdy=-cot(h(2,1))*hd3(1,m+1) + hd3(2,m+1)/sin(h(2,1)); %Translates back out of head direction coordinates. AA=[xph(1,m) xph(1,m)+hdx]; BB=[yph(1,m) yph(1,m)+hdy]; if sum(AA)~=0 && sum(BB)~=0; plot(AA,BB,'k-','LineWidth',1); %AA,BB,'k-', xlim([xmin-10 xmax+10]) ylim([ymin-10 ymax+10]) hold on end end drawnow; %PLOT HD CELLS VS TIME DURING RETRIEVAL axes(ax19b); %Plot HD retrieval versus time. subplot(numPsizesTot*numPtestsTot,1,(numPsizes-1)*numPtestsTot+numPtests,'replace'); ax19b=gca; spikehd3=zeros(length(h),length(hd3)); for ht=1:5:length(hd3); maxhd3=mean(max(hd3)); %Create a spiking version of the HD activity in hd3. spikehd3(:,ht)=rand(length(h),1)>(maxhd3*ones(length(h),1)-hd3(:,ht))/maxhd3; %Generate a spike when rand is greater than 1-hd3. end numSING=length(h); for sing=1:numSING; fps=find(spikehd3(sing,:)); for mm=1:length(fps); %plot([fps(mm) fps(mm)],[sing sing-0.5],'Color','k','LineWidth',1); hold on; if colorSpikes==1; plot([fps(mm) fps(mm)],[sing sing-0.5],'Color',[rc(sing,:)],'LineWidth',1); hold on; elseif colorSpikes==0; plot([fps(mm) fps(mm)],[sing sing-0.5],'Color','k','LineWidth',1); hold on; end end end if colorSpikes==1; set(gca,'Color','k'); elseif colorSpikes==0 set(gca,'Color','w'); end xlim([1 rrrT]); ylim([0 numSING+0.5]); drawnow %PLOT X POS AND RUNNING SPEED FROM INVERSE TRANSFORM OF PHASE DURING %RETRIEVAL rspeed=zeros(1,rrrT); for rt=3:rrrT; rspeed(1,rt)=sqrt((xph(rt)-xph(rt-1))^2 + (yph(rt)-yph(rt-1))^2)/dt; end axes(ax112); tt=[1:rrrT]; plot(tt,xph(1,1:rrrT),tt,rspeed); xlim([1 rrrT]); ylim([-55 55]); halfrun=0 numCircRun(numPsizes,1)=numP; for xyt=1:rrrT-1 if yph(xyt)<-5 && yph(xyt+1)>-5 && xph(xyt)<-10 && halfrun==0; numCircRun(numPsizes,numPtests+1)=numCircRun(numPsizes,numPtests+1)+0.5; halfrun=1; end if yph(xyt)>5 && yph(xyt+1)<5 && xph(xyt)>10 && halfrun==1; numCircRun(numPsizes,numPtests+1)=numCircRun(numPsizes,numPtests+1)+0.5; halfrun=0; end end strengHD numCircRun end %End while numCircRun<2.0 end %End numPsizes end %End numPtests fid=fopen('numCircRunAA','a'); fprintf(fid, '%g\t',numCircRun); fprintf(fid, '\n'); fclose(fid); %END OF SCRIPT