Grid cell spatial firing models (Zilli 2012)

 Download zip file 
Help downloading and running models
This package contains MATLAB implementations of most models (published from 2005 to 2011) of the hexagonal firing field arrangement of grid cells.
1 . Zilli EA (2012) Models of grid cell spatial firing published 2005-2011. Front Neural Circuits 6:16 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism:
Cell Type(s): Entorhinal cortex stellate cell; Abstract integrate-and-fire leaky neuron;
Gap Junctions:
Simulation Environment: MATLAB;
Model Concept(s): Oscillations; Attractor Neural Network; Spatial Navigation; Grid cell;
Implementer(s): Zilli, Eric [zilli at];
% Mhatre, Gorchetchnikov, Grossberg 2010 model
% eric zilli - 20111102 - v1.0
% Mhatre et al. 2010 gave a self-organizing map model of grid cells.
% The input to the grid cells is provided by unbiased ring attractors.
% The firing activity of a cell in a ring is a set of parallel stripes over
% the environment. Each ring has a different orientation, but the same
% spatial scale. Each cell in a ring has a systematically offset spatial
% phase.
% These cells project onto a set of grid cells that form the
% self-organizing map. The grid cells compete for activity so that only one
% can have a high activity at any time. Plasticity is dependent on
% activity, so the weights from the active stripe cell inputs to that
% highly active grid cell will be increased in strength.
% Mhatre et al. gave a geometric argument that a pair of stripe cells will
% be co-active most often if the angle of the two is a multiple of 60
% degrees (well, they say +/-60 degrees +/- 180 degrees which, well, what
% does that mean? I don't know. Could be a typo, I guess.) It seems to be
% true that the network results in strong synapses onto a grid cell from
% inputs close to a 60 degree angle.
% After the simulation the rate maps, spatial autocorrelations, and spatial
% power spectral densities are plotted. I included the latter because I
% think might be a better way to measure the global gridness, spacing,
% and orientation.
% I am not able to reproduce their results using spacings other than 20
% degrees between the orientations of each ring (even setting the
% stripe spacing to 30 cm and b0 to 1.125, per their Figure 8). Possibly
% some other value was also changed to allow the model to learn with more
% closely spaced orientations.
% I was unable to successfully implement the model based on the original
% manuscript, so this script was written in part based on source code from
% the authors themselves (and Dr. Praveen Pilly who is carrying on
% development of this model). In retrospect, part of the problem was that
% the model is dependent on the specific trajectory used (or some aspect of
% it, e.g. the mean velocity) but they did not mention this in the paper so
% I spent a lot of time using the wrong trajectory. Another issue was that
% they failed to give units for their parameters.
% This code is released into the public domain. Not for use in skynet.

%% Simulation variables
% they ran the same 600 s trajectory 5 times
simdur = 600; % s
trajdt = 0.020; % sampling rate for output; s (20 ms is actual sampling rate)
dt = 0.002; % s
ntrials = 5; 10;
tind=1; % index for current time step

% save history of all variables?
% faster and less memory if not
saveHistory = 0;

% if saving history, optionally downsample it
% store value if mod(tind,historySkip)==0
historySkip = 40;
histind = 1;

% if 1, don't clear W
continueSimulation = 0;

%% Trajectory
% This model is trajectory-dependent so we must use the same trajectory
% they did, which came from Sargolini et al. 2006:
% Careful, there are some NaNs in there.
load data/11207-21060501+02_t6c1

% Upsample trajectory if needed:
pos_x = interp1(t(~isnan(x1)),x1(~isnan(x1)),0:dt:t(end));
pos_y = interp1(t(~isnan(y1)),y1(~isnan(y1)),0:dt:t(end));

%% Model variables
% number of ECII cells, which one hopes will become grid cells:
nECII = 5;

% stripe angles:
stripePhis = deg2rad(-80:20:80); % as in paper, 9 unique orientations
% stripePhis = deg2rad(-80:10:80); % they say this works, but it doesn't for me
nstripeOrientations = length(stripePhis);
nstripePhases = 4;
nstripes = nstripeOrientations*nstripePhases;

% Set the orientation of each stripe input
angles = stripePhis;
stripePhis = [];
for phase=1:nstripeOrientations
  stripePhis = [stripePhis angles(phase)*ones(1,nstripes/nstripeOrientations)];

% They didn't bother to give us the units, but they are:
A=3; % decay rate, 1/s
B=1; % excitatory reversal potential, voltage units
D=1.5; % negative of inhibitory reversal potential, voltage units
alpha=17.5; % self-excitatory feedback gain coefficient, unitless
p=1.5; % inhibitory connection strength
lambda_g = 10; % not in paper, but multiplies a few equations -- some sort of rate parameters? 1/s
lambda_w=0.025; % learning rate, wrong value in paper, 1/s
eta=0.4; % habituative transmitter rate, 1/s
beta=0.2; % scales influence of input and feedback on rate of habituation, unitless

spacing = 20; % cm
b0 = 0.0884*spacing; % cm; value from Praveen's script

if ~continueSimulation
  % weight matrix from stripes to ECII
  W = 0.1*rand(nECII,nstripes);
%   figure; imagesc(W) % plot weights before learning

% ECII inhibitory matrix
P = p*(ones(nECII)-eye(nECII));

%% Firing field plot variables
nSpatialBins = 40;
minx = -50; maxx = 50; % cm
miny = -50; maxy = 50; % cm
occupancy = zeros(nSpatialBins);
spikes = zeros(nSpatialBins,nSpatialBins,nECII);

%% Precompute stripe activities
% Project the trajectory velocity onto the preferred direction vectors
H = [cos(deg2rad(angles')) sin(deg2rad(angles'))];
projTraj = H*[pos_x; pos_y];
% Make a copy of the directional displacement for each spatial phase
projTraj = repmat(projTraj, [1 1 nstripePhases]);
% Rearrange dimensions
projTraj = permute(projTraj, [1 3 2]);
% Spatial phases:
mus = repmat(spacing*((1:nstripePhases)-1)/nstripePhases, [nstripeOrientations 1 length(pos_x)]);
% Now evaluate the stripe cells at each point along the trajectory:
strip_fun = @(mu,sigma,lambda,pos)(exp(-(min(mod(pos-mu,lambda),lambda-mod(pos-mu,lambda)).^2)/(2*sigma^2)));
g = strip_fun(mus, b0, spacing, projTraj);

%% Main simulation loop
for trial=1:ntrials
  t=0; % current time in simulation, s
  tind=1; % index for current time step
  if saveHistory
    x = zeros(nstripes,round(simdur/dt)+1);
    v = zeros(nECII,round(simdur/dt)+1);
    z = zeros(nECII,round(simdur/dt)+1);
    pos = zeros(2,round(simdur/dt)+1);
    x = zeros(nstripes,1);
    v = zeros(nECII,1);
    z = zeros(nECII,1);
    pos = zeros(2,1);
  spikeTimes = [];
  spikeCoords = [];
  fprintf('\nStarting trial %d/%d\n',trial,ntrials)
  while t<simdur
    t = t+dt;
    tind = tind+1;
    if mod(t,round(simdur)/10)<=dt
      disp(sprintf('t = %d, %g elapsed',t,toc));
    pos = [pos_x(tind) pos_y(tind)];
    oldx = x(:);
    oldv = v(:);
    oldz = z(:);
    % update stripe activity for current location:
    x = reshape(g(:,:,tind),[],1);

    % voltage
    oldvsquared = oldv.^2;
    v = v + lambda_g*dt*(-A*oldv + (B-oldv).*((W*oldx) + alpha*oldv.^2).*oldz - (D+oldv).*(P*(oldv'.^2)'));

    % "habituative transmitter gate"
    z = z + lambda_g*dt*(eta*((1-oldz) - beta*oldz.*((W*oldx) + alpha*oldv.^2).^2));

    % learning rule: "competitive instar"
    W = W + dt*lambda_w*repmat(oldv.^2,1,nstripes).*((oldx*(2-sum(W,2)'))' - W.*repmat(sum(oldx)-oldx',nECII,1));
    % that rule is a little odd: each synapse knows the true value of all
    % the inputs to the cell
    % Save spiking for ratemap
    xindex = floor((pos(1)-minx)/(maxx-minx)*nSpatialBins)+1;
    yindex = floor((pos(2)-miny)/(maxy-miny)*nSpatialBins)+1;
    occupancy(yindex,xindex) = occupancy(yindex,xindex) + dt;
    spikes(yindex,xindex,:) = squeeze(spikes(yindex,xindex,:)) + oldvsquared;
    if saveHistory && mod(tind,historySkip)==0
      histind = histind+1;
      poshist(:,histind) = pos;
      xhist(:,histind) = x;
      vhist(:,histind) = v;
      zhist(:,histind) = z;

%% Plot weight matrix at end. 60 degree differences should be noticeable.
title('Stripes to ECII weight matrix');
ylabel('ECII cell');
xlabel('Stripe input (grouped by phase)');

%% Plot where voltage of each cell exceeded a threshold
if saveHistory
  figure('name','trajectory and ECII cell firing');
  for i=1:5
    axis equal

%% Plot rate map, spatial autocorrelation, and spatial power spectral density
% matrix to normalize 2d autocorrelation to yield the unbiased estimator
xcorr2NormMatrix = [1:length(occupancy) (length(occupancy)-1):-1:1];
xcorr2NormMatrix = xcorr2NormMatrix'*xcorr2NormMatrix;
for i=1:5
  figure('position',[520 378 560 225],'name',['Cell ' num2str(i)]);
  gaussWindow = fspecial('gaussian',[5 5],1);
  smoothedRateMap = conv2(spikes(:,:,i)./(occupancy+eps),gaussWindow,'same');
  title('rate map');
  axis equal;
  xlim([0 40]); ylim([0 40])
  unbiasedSpatialAutoCorr = xcorr2(smoothedRateMap)./xcorr2NormMatrix;
  axis equal;
  xlim([0 79]); ylim([0 79])
  title('spatial autocorrelogram');
  spatialPSD = abs(fftshift(fft2(unbiasedSpatialAutoCorr)));
  title({'spatial power spectral','density amplitude'});
  axis equal;
  xlim([0 79]); ylim([0 79])