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];
% Burgess 2008's oscillatory interference models
% eric zilli - 20110909 - v1.0
% Burgess 2008 described a wide chunk of theory relating to oscillatory
% interference grid cell models, both reviewing the previous work on
% these models and pointing out many simple variations on the basic model.
% The variations include:
%   * Additive vs. multiplicative interactions of the active
%       oscillators/VCOs
%   * Sinusoidal or "punctate" (spike-like) oscillator output
%   * Combining oscillators with or without leaky integration
%   * Using baseline modulation or not (showing baseline modulation is
%       necessary, at least when the oscillations are summed rather than
%       multiplied, cf. Hasselmo2008.m)
%   * Directional or nondirectional VCOs (whether VCOs output when
%       preferred direction is more than 90 degrees from current heading)
% and we provide options to simulate each of these variations.
% Importantly: the grid cell threshold depends on the exact combination of
% options used, and while I started tryin to pre-program them all, there
% are just too many! It is very easy to find the threshold needed for
% any particular case though:
%   * First, set the simulation duration to 10 seconds or so and run it.
%   * Then do figure; plot(fhist) if it isn't already plotted and select
%       an appropriate threshold value a bit below the highest peaks sin
%       that plot but well above the lowest peaks.
% The manuscript suggests one other variation where the frequency
% of the baseline oscillation is the average of the frequencies
% of the active VCOs over all preferred directions. It is not clear though
% how one oscillator could take the average of the frequencies of a number
% of other oscillators and use that as its own frequency.
% The manuscript includes a precession model using 6 directional VCOs
% (see others in BurgessEtAl2007_precession.m and Hasselmo2008.m). The fix
% isn't so plausible, though. The directional fix used here blocks the
% output from a VCO onto a grid cell if the animal's heading is not within
% 90 degrees of the VCO's preferred direction, but the VCO still integrates
% velocity even when its output is blocked. In theory this could be
% implemented, e.g.
%   A. as a cell that has nonlinear limit cycle subthreshold oscillations
%       that transition into spiking at a frequency equal to the baseline
%       frequency, but it seems unlikely that such an oscillator could
%       have its frequency correctly controlled in the post-spike period,
%       so it may not work.
%   B. the output from this cell could pass through an intermediate cell
%       that is inhibited by the opposite direction head direction input,
%       but adding in extra ad hoc circuiry to produce phase precession
%       is not particularly compelling
% This manuscript also fixes a problem with the Burgess et al. 2007 model
% where the product over the pairs of oscillators was being
% thresholded as a whole, rather than its individual components
% being thresholded as is done here. Though really it is only a problem in
% as much as it lets two negative oscillations multiply together into a
% positive value, which does not seem biologically plausible.
% Finally, this manuscript presents the model's equations in a form
% suitable for arbitrary input trajectories (though Giocomo et al. 2007
% was probably historically the earliest paper to do so).
% Note: This code is not necessarily optimized for speed, but is meant as
% a transparent implementation of the model as described in the manuscript.
% This code is released into the public domain. Not for use in skynet.

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

% are the oscillators combined with a sum or product?
oscillatorInteraction = 'sum'; % 'sum' or 'product'
% are the active VCOs or baseline sinusoidal or spike/delta-shaped?
VCOOutput = 'sine'; % 'sine' or 'spike'
baselineOutput = 'sine'; % 'sine' or 'spike'

% does the grid cell reflect only its immediate inputs or does it integrate
% them over time? integration makes more sense with
% oscillatorInteraction = 'sum' but if it is 'prod', we'll multiply the
% oscillations first and let the result be leakily-integrated. NB that
% though Burgess states that leaky integration is a reasonable assumption
% for entorhinal neurons, it is specifically untrue for stellate cells,
% which are famously resonant and so have a more complex response than
% simple integration.
integrateGridInputs = 0; % if =1, grid cell is a leaky integrator; if =0, grid cell only responds to immediate inputs
baselineModulation = 1; % if =0, no baseline oscillation is used; if =1, a baseline is used
directionalVCOs = 0; % if =0, VCOs output to grid cell for any heading; if =1, VCOs only output if heading is within 90 degrees of their preferred direction

nVCOs = 3;

% Threshold to determine when the grid cell has spiked
spikeThreshold = 3.5;

% if =0, just give constant velocity. if =1, load trajectory from disk
useRealTrajectory = 1;
constantVelocity = 1*[.5; 0*0.5]; % m/s

%% Simulation parameters
dt = .02; % time step, s
simdur = 200; % total simulation time, s
tind = 1; % time step number for indexing
t = 0; % simulation time variable, s
x = 0; % position, m
y = 0; % position, m

%% Model parameters
ncells = 1;
% Baseline maintains a fixed frequency
baseFreq = 6; % Hz
% Directional preference of each VCO (this also sets the number of VCOs)
dirPreferences = (0:nVCOs-1)*pi/3;
% Scaling factor relating speed to oscillator frequencies
% NB paper uses 0.05*2pi rad/cm. But we do the conversion to rad later,
% leaving 0.05 Hz/cm = 5 Hz/m which produces very tight field spacing. For
% cosmetic purposes for the trajectory we use here, we'll use beta = 2.
beta = 2; % Hz/(m/s) 
T = 0.010; % s, time constant for integration of oscillations if used
C = 1/0.5; % normalization coefficient. the integral works out to be pi*(2n-3)!!/(2^(n-1)*(n-1)!) by my reckoning
if integrateGridInputs && dt>(0.05*T)
  warning('Your time step dt should be much smaller than the time constant T!')
  dt = 0.05*T;
  fprintf('Using dt = %g, but frankly it should be smaller!\n',dt)

%% History variables
speed = zeros(1,ceil(simdur/dt));
curDir = zeros(1,ceil(simdur/dt));
vhist = zeros(1,ceil(simdur/dt));
fhist = zeros(1,ceil(simdur/dt));
f = 0;

%% Firing field plot variables
nSpatialBins = 60;
minx = -0.90; maxx = 0.90; % m
miny = -0.90; maxy = 0.90; % m
occupancy = zeros(nSpatialBins);
spikes = zeros(nSpatialBins);

spikeTimes = [];
spikeCoords = [];
spikePhases = [];

%% Initial conditions
% Oscillators will start at phase 0:
VCOPhases = zeros(1,length(dirPreferences)); % rad
basePhase = 0; % rad

%% Make optional figure of sheet of activity
if livePlot
  h = figure('color','w','name','Activity of one cell');
  if useRealTrajectory
    set(h,'position',[520 378 1044 420])

%% Possibly load trajectory from disk
if useRealTrajectory
  load data/HaftingTraj_centimeters_seconds.mat;
  % interpolate down to simulation time step
  pos = [interp1(pos(3,:),pos(1,:),0:dt:pos(3,end));
  pos(1:2,:) = pos(1:2,:)/100; % cm to m
  vels = [diff(pos(1,:)); diff(pos(2,:))]/dt; % m/s
  x = pos(1,1); % m
  y = pos(2,1); % m

%% !! Main simulation loop
fprintf('Simulation starting. Press ctrl+c to end...\n')
while t<simdur
  tind = tind+1;
  t = dt*tind;
  % Velocity input
  if ~useRealTrajectory
    v = constantVelocity; % m/s
    v = vels(:,tind); % m/s
  curDir(tind) = atan2(v(2),v(1)); % rad
  speed(tind) = sqrt(v(1)^2+v(2)^2); % m/s
  x(tind) = x(tind-1)+v(1)*dt; % m
  y(tind) = y(tind-1)+v(2)*dt; % m
  % VCO frequencies are pushed up or down from the baseline frequency
  % depending on the speed and head direction, with a scaling factor beta
  % that sets the spacing between the spatial grid fields.
  VCOFreqs = baseFreq + beta*speed(tind)*cos(curDir(tind)-dirPreferences); % Hz
  % Advance oscillator phases
  % Radial frequency (2pi times frequency in Hz) is the time derivative of phase.
  VCOPhases = VCOPhases + dt*2*pi*VCOFreqs; % rad
  basePhase = basePhase + dt*2*pi*baseFreq; % rad
  % Compute active oscillation values from current phase
  if strcmp(VCOOutput,'sine')
    VCOs = cos(VCOPhases);
  elseif strcmp(VCOOutput,'spike')
    VCOs = (1/2+cos(VCOPhases)/2).^50;
  if directionalVCOs
    % "Directional" VCOs do not output unless they
    % are within 90 degress of the VCO's preferred direction
    VCOs = VCOs.*(cos(curDir(tind)-dirPreferences)>0);
  % Compute baseline oscillation from current phase
  if strcmp(baselineOutput,'sine')
    baseline = baselineModulation*cos(basePhase);
  elseif strcmp(baselineOutput,'spike')
    baseline = baselineModulation*(1/2+cos(VCOPhases)/2).^50;
  % Sum each dendritic oscillation separately with the baseline oscillation
  VCOsPlusBaseline = VCOs + baseline;
  % Threshold individual oscillator pairs
  VCOsPlusBaseline = VCOsPlusBaseline.*(VCOsPlusBaseline>0);
  % Final activity is the sum or product of the thresholded oscillations.
  if strcmp(oscillatorInteraction,'sum')
    act = sum(VCOsPlusBaseline);
  elseif strcmp(oscillatorInteraction,'product')
    act = prod(VCOsPlusBaseline);
  if integrateGridInputs
    % final output f is a leaky integration of the combined oscillations
    f = f + dt*(-f/T + act*C);
    % final output f is simply the current value of the combined oscillations
    f = act;
  % Save for later
  fhist(tind) = f;
  % Save firing field information
  if f>spikeThreshold
    spikeTimes = [spikeTimes; t];
    spikeCoords = [spikeCoords; x(tind) y(tind)];
    spikePhases = [spikePhases; basePhase];
  if useRealTrajectory
    xindex = round((x(tind)-minx)/(maxx-minx)*nSpatialBins)+1;
    yindex = round((y(tind)-miny)/(maxy-miny)*nSpatialBins)+1;
    occupancy(yindex,xindex) = occupancy(yindex,xindex) + dt;
    spikes(yindex,xindex) = spikes(yindex,xindex) + double(f>spikeThreshold);
  if livePlot>0 && (livePlot==1 || mod(tind,livePlot)==1)
    if ~useRealTrajectory
      xlabel('Time (s)')
      axis square
      title(sprintf('t = %.1f s',t))
      hold on;
      if ~isempty(spikeCoords)
        cmap = jet;
        cmap = [cmap((end/2+1):end,:); cmap(1:end/2,:)];
        phaseInds = mod(spikePhases,2*pi)*(length(cmap)-1)/2/pi;
        pointColors = cmap(ceil(phaseInds)+1,:);
        scatter3(spikeCoords(:,1), ...
                 spikeCoords(:,2), ...
                 zeros(size(spikeCoords(:,1))), ...
                 30*ones(size(spikeCoords(:,1))), ...
                 pointColors, ...
      axis square
      title({'Trajectory (blue) and',...
             'spikes (colored by theta phase',...
             'blues before baseline peak, reds after)'})
      hold on;
      plot([0 tind-1]*dt,[spikeThreshold spikeThreshold],'r')
      title('Activity (blue) and threshold (red)');
      xlabel('Time (s)')
      axis square
      axis square
      title({'Rate map',sprintf('t = %.1f s',t)})
      hold on;
      if ~isempty(spikeCoords)
        cmap = jet;
        cmap = [cmap((end/2+1):end,:); cmap(1:end/2,:)];
        phaseInds = mod(spikePhases,2*pi)*(length(cmap)-1)/2/pi;
        pointColors = cmap(ceil(phaseInds)+1,:);
        scatter3(spikeCoords(:,1), ...
                 spikeCoords(:,2), ...
                 zeros(size(spikeCoords(:,1))), ...
                 30*ones(size(spikeCoords(:,1))), ...
                 pointColors, ...
      axis square
      title({'Trajectory (blue) and',...
             'spikes (colored by theta phase',...
             'blues before baseline peak, reds after)'})