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];
% Blair, Gupta, and Zhang 2008's ring attractor oscillatory interference model
% eric zilli - 20110908 - v1.0
% This manuscript deals with both grid and place cells, but this
% implementation focuses only on the grid cell component.
% (Notice this is the 1D version of the model, so the output is just a
%   series of blue bumps. Use BlairEtAl2008_2D.m for a 2D grid as output).
% Blair et al. 2008 overcomes one of the limitations of the Blair et al.
% 2007 model in that the existance of the grid pattern is no longer assumed
% in the inputs. Instead, this model draws inspiration from oscillatory
% interference models(*) and creates the (1D) repetitive spatial pattern from
% a pair of oscillators. (At the very end they sketch out a 2D version
% using one reference oscillator and multiple ring attractors that
% integrate velocity along directions at 60 degree angles, as in other
% 2D interference models, see BlairEtAl2008_2D.m).
% (*) That's sort of putting words in their mouth, but ultimately it does
% seem like ring attractor models originally came out of interference
% models if you go back to O'Keefe 1991, Burgess et al. 1992 and then
% forward to Touretzky et al. 1993 and Redish et al. 1994 who started using
% so-called "sinusoidal arrays" that look like the first steps toward
% ring attractor models, at least in the medial temporal lobe. Perhaps they
% go back much father!
% Like Blair et al. 2007, this model does not actually simulate the
% path integration. Instead, on each time step the ring attractor
% activities are set to the values that would occur if path integration had
% been carried out (calculated by translating the trajectory velocity
% signal to a frequency signal and integrating the frequencies to yield the
% phase of each attractor oscillator as a function of time).
% As a result, this script has little to do and the model amounts to
% literally multiplying together cosine gratings to make the hexagonal
% pattern.
% This code is released into the public domain. Not for use in skynet.

% if =0, uses a line segment as trajectory
% if =1, load trajectory from disk and plot activity only along that
useRealTrajectory = 1;

%% Model parameters
f0 = 7; % reference oscillator freq, Hz
nCellsPerRing = 6;
% phase offsets for each cell in the ring
phi = (1:nCellsPerRing)*2*pi/nCellsPerRing;
% inverse gain factors for ring attractor i's velocity input:
ilambda = [0 1/60];

%% Trajectory
% we're just using 1D trajectories here:
if useRealTrajectory
  % compute activity at each position in trajectory
  % load trajectory from disk:
  load data/HaftingTraj_centimeters_seconds.mat;
  % package the trajectory the way we want it:
  % this is basically the "simulation duration" for this script:
  nTrajectorySamples = length(pos(1,:));
  x = pos(1,1:nTrajectorySamples); % cm
  dT = pos(3,2)-pos(3,1); % s
  clear pos;
  nTrajectorySamples = 10e3;
  x = linspace(0,1000,nTrajectorySamples); % cm
  dT = 0.002; % s

%% Calculate velocity inputs
vel = diff(x)/dT; % cm/s

%% Oscillator frequencies are a baseline plus velocity times gain
% Note F is nTrajectorySamples-by-2 because vel(:) is nTrajectorySamples-by-1
% and lambda is 1-by-2.
F = f0 + vel(:)*ilambda; % Hz

%% Integrate frequencies to get phase offsets for each ring:
alpha = 2*pi*cumsum(F)*dT; % nTrajectorySamples-1--by--2, rad
% Take each ring's base phase offset and make one copy for each cell in each ring
alpha = [repmat(alpha(:,1),1,nCellsPerRing) repmat(alpha(:,2),1,nCellsPerRing)];

%% Calculate ring cell activities at each point in time
% This takes the base phase offset and adds the appropriate value from
% phi to find the phase offset for each individual cell in a ring.
theta = (sin(alpha + repmat(phi,nTrajectorySamples-1,2))+1)/2;

%% Grid output is the product of two ring cell activities
% We'll use the first cell in each ring
gridActivity = theta(:,1).*theta(:,nCellsPerRing+1);

%% Plot the activity and trajectory
figure; plot(x(1:end-1),gridActivity)
title('Activity vs. position')
xlabel('1D Position (cm)')