Synthesis of spatial tuning functions from theta cell spike trains (Welday et al., 2011)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:129067
A single compartment model reproduces the firing rate maps of place, grid, and boundary cells by receiving inhibitory inputs from theta cells. The theta cell spike trains are modulated by the rat's movement velocity in such a way that phase interference among their burst pattern creates spatial envelope function which simulate the firing rate maps.
Reference:
1 . Welday AC, Shlifer IG, Bloom ML, Zhang K, Blair HT (2011) Cosine directional tuning of theta cell burst frequencies: evidence for spatial coding by oscillatory interference. J Neurosci 31:16157-76 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s): Hippocampus CA1 pyramidal GLU cell; Hippocampus CA3 pyramidal GLU cell; Entorhinal cortex stellate cell;
Channel(s): I Na,p;
Gap Junctions:
Receptor(s): GabaA; AMPA;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEURON; MATLAB;
Model Concept(s): Synchronization; Envelope synthesis; Grid cell; Place cell/field;
Implementer(s): Blair, Hugh T.;
Search NeuronDB for information about:  Hippocampus CA1 pyramidal GLU cell; Hippocampus CA3 pyramidal GLU cell; GabaA; AMPA; I Na,p; Gaba; Glutamate;
%%%CREATES TWELVE SIMULATED THETA CELL SPIKE TRAINS FROM COSINE VCOs
%%%WHICH WILL INTERFERE WITH ONE ANOTHER TO FORM THE BORDER CELL IN FIG. 7F

load 'trackingdata';
pixpercm = 5; %pixels per cm in the tracking data

dX=diff(rsX); 
dY=diff(rsY); 
dX=[dX(1); dX]';
dY=[dY(1); dY]';
rsS=sqrt(dX.^2 + dY.^2)/(pixpercm*(rsTS(2)-rsTS(1))); %compute speed for modulating base frequency

%%%%% upsample tracking rate from 30 Hz to 500 Hz
tvec=interp1(rsTS, rsTS, [rsTS(1):.002:(rsTS(length(rsTS)))]); %rat path tracking timestamp vector 
pvec_x=interp1(rsTS, rsX, tvec); clear rsX; %rat path tracking X coordinate 
pvec_y=interp1(rsTS, rsY, tvec); clear rsY; %rat path tracking Y coordinate 
ratspeed=interp1(rsTS, rsS, tvec); clear rsS; %rat path tracking instantaneous speeds
clear rsTS; %delete 30 Hz timestamp vector and speed vector


base_freq = 7; %VCO base frequency in Hz
speedslope = .025; %slope for linear dependence of base_freq upon running speed, in units of Hz/cm/s

%preferred vector lengths
bbb=1.8; %base preferred movement vector length
rhovec = [bbb bbb^2 bbb^3 bbb^4 bbb bbb^2 bbb^3 bbb^4 bbb bbb^2 bbb^3 bbb^4];
%preferred vector orientations in radians
dirvec = [0 0 0 0  .2618*2 .2618*2 .2618*2 .2618*2 -.2618*2 -.2618*2 -.2618*2 -.2618*2];
%VCO starting phases in radians    
startphases = [-0.0000   -0.0000   -0.0000   -0.0000   -0.2500   -0.5000   -1.0000   -2.0000    -0.2500   -0.5000   -1.0000   -2.0000];
%shvec is a polar vector [r theta] specifying distance and direction by which to  the envelope function ([0 0] for no shift)
shvec = [40 0];
%adjust startphases to shift envelope as specified by shvec
startphases=startphases-[rhovec*shvec(1)*.0125].*cos([dirvec-deg2rad(shvec(2))]); 

%center the tracking data at the origin of the plane
center_x=(min(pvec_x)+max(pvec_x))/2;
center_y=(min(pvec_y)+max(pvec_y))/2;

%convert coordinate values from pixels to cm
pvec_x=(pvec_x-center_x)/pixpercm;
pvec_y=(pvec_y-center_y)/pixpercm;

%%%make sure position vectors are column vectors
if (size(pvec_x,1)==1)
    pvec_x=pvec_x';
end
if (size(pvec_y,1)==1)
    pvec_y=pvec_y';
end

dt=tvec(2)-tvec(1); %integration time step (assumed constant)

%%compute the speed-dependent base frequency (lfpfreq) and phase (lfpphase) at each time step
lfpfreq=ones(1,length(pvec_x)-1).*((base_freq+ratspeed(1:(length(ratspeed)-1))*speedslope)*dt); 
lfpphase=cumsum([0 2*pi*lfpfreq]);

cullperc=.25; %theta cell spiking rate is determined by the percentage of "culled" spikes

for i=1:length(rhovec) %loop through VCOs
     pvec=[pvec_x pvec_y]*[[cos(dirvec(i)) sin(dirvec(i))];[-sin(dirvec(i)) cos(dirvec(i))]]; %pvec is rat's position along the VCO's preferred vector
     pvec=pvec(:,1); %drop y coordinates (x is position along the preferred vector)
     fvvec=diff(pvec)*dt; %speed along preferred vector
     foff=rhovec(i)*fvvec; %frequency of each VCO is offset from base_freq by a velocity-dependent value
     thfreq=lfpfreq'+foff;

    VCO(i).thphase=cumsum([startphases(i) [2*pi*(lfpfreq'+foff)]']);
    VCO(i).spiketrain=(((cos(VCO(i).thphase)+1).*(rand(size(VCO(i).thphase,1),size(VCO(i).thphase,2))))>.5); 
    ggg=find(VCO(i).spiketrain>0); ggh=rand(1,length(ggg)); ggg=ggg(find(ggh<cullperc)); 
    VCO(i).spiketrain(ggg)=0;
    VCO(i).spiketrain=tvec(find(VCO(i).spiketrain>0))-tvec(1);
    VCO(i).spiketrain=VCO(i).spiketrain(find(VCO(i).spiketrain>0));
    ggh=rand(1,length(VCO(i).spiketrain)); ggg=find(ggh<.5);
    VCO(i).spiketrain=VCO(i).spiketrain(ggg); %knock down to 50 Hz firing rate for theta cell
    if i==1
        mempot=cos(VCO(i).thphase);
    else
        mempot=mempot+cos(VCO(i).thphase);
    end
end


%save theta cell spike train time stamps (in units of milliseconds for NEURON)
 if (1)
     for i=1:length(rhovec)
         spkts0_r1=VCO(i).spiketrain'*1000; 
         save(['boundarytheta_' num2str(i) '.dat'], 'spkts0_r1', '-ascii');
     end
 end

%plot graph approximating what NEURON simulation output should look like
threshlevel = 0.75;
spos=[];
 tout=10;
 thresh=max(mempot)*threshlevel;
 for t=1:length(tvec)
     if (mempot(t)>thresh)
         if (tout==10)
             spos=[spos; [pvec_x(t) pvec_y(t)]];
             tout=1;
         else
             tout=tout+1;
         end
     else
         tout=10;
     end
 end

figure(1);
hold off;
plot(pvec_x,pvec_y,'-k');
hold on;
scatter(spos(:,1),spos(:,2),'.r')
axis tight; axis square;


Loading data, please wait...