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];
% Gaussier, Banquet, Sargolini, Giovannangeli, Save, Poucet 2007's grid cell model
% eric zilli - 20111110 - v1.0
% In this model positional information is first encoded (and path
% integrated) in the firing rate of individual neurons D with each with its
% own preferred direction theta. These cells are leak-less integrators of
% their directional velocity inputs. A cell with preferred direction theta
% starts at firing rate 0 and increases its firing rate in proportion to
% the distance the animal travels along theta (and decreases when moving
% opposite theta where activities seem to be allowed to go negative). Two
% of these cells are used here, with preferred directions at 0 and 60
% degrees.
% Each of these cells' activity is discretized into second populations of
% cells E (one E population for each D cell). This is carried out by
% defining a lower and upper limit of a cell's activity (the size of the
% largest area that can be path integrated), say 0 and 10, and a number of
% cells in the discretization population E, say 20. Then the active E cell
% is based directly on the firing rate of the input D cell. If the input
% activity were between 0 and 0.5 then E cell 1 could be active. With input
% between 0.5 and 1, cell 1, etc. up to inputs rates 9.5 to 10 activating E
% cell 20. This reads out the single-valued position from the input cell D
% into a population-level encoding of position E. One one cell in each E
% is active at any time, so the spatial firing pattern of an E cell is
% a single stripe parallel to the preferred direction of its input D cell
% at a distance from the animal's starting point according to the firing
% rate range it encodes from its input D cell.
% The stripe-like population-level encoding of position in E is then
% "folded" down into a smaller set of cells M in a manner like the modulo
% operator. When folding down cells to a set of four cells, input cell 1,
% 5, 9, 13, etc. all synaptically project to modulo cell 1. Input cells 2,
% 6, 10, 14, etc. to modulo cell 2, and so forth for the inputs of modulo
% cells 3 and 4. This compresses the stripe-like population code E into a
% modular code M.
% These cells M would fire like band or stripe cells used in other models,
% with spatial fields in the form of a series of thin parallel stripes
% across the environment, where the spacing depends on the relative sizes
% of E and M
% The activity of the grid cells is then given by assigning each grid
% cell a preferred cell in each modular population and multiplying the
% inputs so the grid cell is only active when both of its inputs are.
% The manuscript gives a method for self-organizing the specific spatial
% phases learned by each grid cell in appendix B.2. The model uses a very
% simple one-shot learning mechanism that adds a new grid cell every time
% the current combination of modulo population activities has not been
% previously learned. 
% Note that in the manuscript they seemed to discretize the velocity inputs
% themselves, which seemed to produce a large enough path integration error
% that they had to add a place-driven resetting mechanism, but I'm not
% going to discretize the input, so I won't use the resetting (to keep this
% simpler).
% One quirk of the paper: they seem to focus on entorhinal layer II pyramidal
% cells rather than the more populous stellate cells (perhaps just out of
% habit of referring to principal cells in cortex as pyramidals).
% This code is released into the public domain. Not for use in skynet.

clear all

% if >0, plots the sheet of activity during the simulation on every livePlot'th step
livePlot = 20;

% if =0, uses a line segment as trajectory
% if =1, load trajectory from disk and plot activity only along that
useRealTrajectory = 1;
constantVelocity = 1*[.5; 0*0.5]; % m/s

%% Simulation variables
simdur = 200; % s
dt = .02; % s
t = 0; % s
tind = 1;
spikeind = [1 1 1];

%% Model parameters

preferredDirections = [0 2*pi/3]';
nDirs = length(preferredDirections);
H = [cos(preferredDirections) sin(preferredDirections)];

% Minimum and maximum "activity" of D for discretizing into E:
minD = -2;
maxD = 2;
D = zeros(nDirs,1+ceil(simdur/dt));
% Number of cells used to discretize each D(i)
nE = 32;
E = zeros(nE*nDirs,1+ceil(simdur/dt));
% Number of modulo cells (spacing is nM/nO)
nM = 4;
M = zeros(nM*nDirs,1+ceil(simdur/dt));
% Number of grid cells (one for each pair of inputs from the mod populations)
% TODO should this be larger for nDirs>2?
nO = nM*nM;
O = zeros(nO,1+ceil(simdur/dt));

% tuning factor that controls which grid cells can be active
T = zeros(nO,1);
% vigilance parameter for recruiting new cells
vig = 0.9;

% O threshold for saving a spike for the spikes+trajectory plot
spikeThreshold = 0.8;

%% Synaptic weight matrices
% Weights from modulo populations to grid cells start random here
Wmo = .01*rand(nO,nDirs*nM);
% Modulo weights from discretized activity cells to modulo populations
% gotta be a better way of doing this!
Wem = [];
for i=1:nDirs
  Wem = blkdiag(Wem, repmat(eye(nM),1,nE/nM));

%% Firing field plot variables
watchCell = 4; % which cell's spatial activity will be plotted?
nSpatialBins = 60;
minx = -0.90; maxx = 0.90; % m
miny = -0.90; maxy = 0.90; % m
occupancy = zeros(nSpatialBins);
spikes = zeros(nSpatialBins);
spikeTimes = [];
spikeCoords = [];
spikeCoords2 = [];
spikeCoords3 = [];

%% Trajectory
if useRealTrajectory
  load data/HaftingTraj_centimeters_seconds.mat;
  % interpolate down to simulation time step
  stopAt = find(pos(3,:)>simdur,1,'first');
  pos = [interp1(pos(3,:),pos(1,:),0:dt:pos(3,stopAt));
  pos(1:2,:) = pos(1:2,:)/100; % cm to m
  vels = [diff(pos(1,:)); diff(pos(2,:))]/dt; % m/s
  x = pos(1,:); % m
  y = pos(2,:); % m
  'oops not done yet'

if livePlot
  figh = figure('position',[520 463 929 335]);

while t<simdur
  tind = tind+1;
  t = tind*dt;

  %% Integrate velocity inputs along preferredDirections matrix H
%   D = cumsum(H*vels)*dt; % do all at once
  D(:,tind) = D(:,tind-1) + H*vels(:,tind)*dt;

  %% Discretize activity of cells D(i) into populations E1 and E2
  for i=1:nDirs
    E(ceil(interp1([minD maxD],[1 nE],D(i,tind)))+(i-1)*nE,tind) = 1;

  %% Project E onto the modulo populations
  M(:,tind) = Wem*E(:,tind);
  % Make a copy of the modulo inputs for each cell, for comparing to
  % its weights.
  Ii = repmat(M(:,tind)',nO,1);
  % Activation rule (equation B.4)
  thresh = @(x)(x.*double(x>0));
  O(:,tind) = T.*thresh(1-sum(abs(Wmo-Ii),2)/(nDirs*nM));
  %% Apply local competition:
  Oprime = O(:,tind).*(O(:,tind)>vig);
  % TODO local competition, they don't give a value for d_max... but it
  % works without local competition anyhow.

  if max(O(:,tind))<vig
    % "The learning and recruitment of a new neuron is triggered only if
    % the activity of the winner neuron is inferior to a vigilance
    % parameter vig."
    newCell = find(T==0,1,'first');
    if isempty(newCell)
    T(newCell) = 1;
    Oprime(newCell,tind) = 1;
    % This is equation B.7. It is not clear whether only the newly
    % recruited neuron learns or whether all cells do. It sounds like
    % only the winner.
    % Winner only:
    Wmo(newCell,:) = Wmo(newCell,:) + (Ii(newCell,:)-Wmo(newCell,:)).*Oprime(newCell,tind);
%     % All cells:
%     Wmo = Wmo + (Ii-Wmo).*O(:,tind);
  % Save firing field information
  if O(watchCell,tind)>spikeThreshold
    if mod(spikeind(1),1000)==1
      spikeTimes(spikeind(1)+1000) = 0;
      spikeCoords(spikeind(1)+1000,:) = [0 0];
    spikeTimes(spikeind(1)) = t;
    spikeCoords(spikeind(1),:) = [x(tind) y(tind)];
    spikeind(1) = spikeind(1)+1;
  if O(watchCell+1,tind)>spikeThreshold
    if mod(spikeind(2),1000)==1
      spikeCoords2(spikeind(2)+1000,:) = [0 0];
    spikeCoords2(spikeind(2),:) = [x(tind) y(tind)];
    spikeind(2) = spikeind(2)+1;
  if O(watchCell+2,tind)>spikeThreshold
    if mod(spikeind(3),1000)==1
      spikeCoords3(spikeind(3)+1000,:) = [0 0];
    spikeCoords3(spikeind(3),:) = [x(tind) y(tind)];
    spikeind(3) = spikeind(3)+1;
  if useRealTrajectory
    xindex = floor((x(tind)-minx)/(maxx-minx)*nSpatialBins)+1;
    yindex = floor((y(tind)-miny)/(maxy-miny)*nSpatialBins)+1;
    occupancy(yindex,xindex) = occupancy(yindex,xindex) + dt;
    spikes(yindex,xindex) = spikes(yindex,xindex) + O(watchCell,tind);
  if livePlot>0 && (livePlot==1 || mod(tind,livePlot)==1)
      xlabel('Time (s)')
      ylabel('Grid cell #')
      axis square
      title(sprintf('t = %.1f s',t))
      axis square
      title({'Weights from modulo pop.','to grid cells O'});
      xlabel('Modulo cell #')
      ylabel('Grid cell #')

      hold on;
      if ~isempty(spikeCoords2)
      if ~isempty(spikeCoords3)
      if ~isempty(spikeCoords)
      axis square
      title({'Trajectory (blue) and','spikes of 3 grid cells'})