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];
% Fuhs and Touretzky (2006)'s developmental model
% eric zilli - 20111214 - v1.0
% This script demonstrates Fuhs and Touretzky (2006)'s developmental model
% of the formation of symmetric rings of output from a cell onto the other
% cells in the rectangular sheet of cells.
% In the model, sinusoidal gratings called "wave packets" flow across the
% sheet of cells, each beginning with a random orientationg and moving
% from one side to the other, in a direction perpendicular to the
% orientation of its stripes.
% The BCM-like learning rule associates a cell to the cells coactive with
% it, so as the wave packets pass any individual cell, the sinusoidal
% grating gets lightly stamped into its outputs. As gratings of random
% orientations keep passing over it, its synaptic weights contain the
% combination of gratings at every orientation, which eventually produces
% symmetric rings of synaptic output.
% The desired ring pattern should start to appear once 20-30 packets
% have passed (depending on the random orientations they had).
% This code is released into the public domain. Not for use in skynet.

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

% Weight matrix options
useCurrentW = 0; % if W exists in namespace, use that one instead of loading/generating
loadWIfPossible = 0;
saveW = 0; % save W to disk after generating

%% Network/Weight matrix parameters
n = 62;
ncells = n*n; % total number of cells in network

% frequency of wave packets
kappa = 9*pi/31;
wavePerPacket = 3;
% tonic firing rate of cells during development (wave packets push activity
% from 0 to 2)
ftonic = 1;

%% Simulation parameters
dt = 1; % time step, ms
npackets = 1000;
simdur = 1e3; % time per packet simulation, ms
t = 0; % simulation time variable, ms
tind = 0; % time step number
x = 0; % position, cm
y = 0; % position, cm

%% Initial conditions
f = zeros(1,ncells); % initial firing rate

%% Make x a 2-by-ncells vector of the 2D cell positions on the neural sheet
% plot as: figure; imagesc(reshape(cellDists,n,n))
x = linspace(-round(n/2),round(n/2),n);
[X,Y] = meshgrid(x,x);
x = [reshape(X,1,[]); reshape(Y,1,[])];
cellDists = sqrt(x(1,:).^2 + x(2,:).^2); % distance from (0,0)

% convert to polar coordinates:
% get angles of cells around (0,0)
theta = atan2(x(2,:),x(1,:));

%% Optionally load weight matrix
fname = sprintf('data/W_FT2006_dev_n%d.mat',ncells);
if ~(useCurrentW && exist('W','var'))
  if loadWIfPossible && exist(fname,'file')
    fprintf('Attempting to load pre-generated W...\n')
    fprintf('+ Loaded pre-generated W.\n')
    W = zeros(ncells);

%% Make optional figure of sheet of activity
if livePlot
  h = figure('color','w','name','Activity of sheet of cells on brain''s surface');
  set(h,'position',[520 378 1044 420])
%% !! Learn weight matrix (main loop)
fprintf('Development starting. Press ctrl+c to end...\n')
for packet=1:npackets
  t = 0;
  tind = -40;
  % random direction for this packet:
  packetDir = 2*pi*rand;
  % learning rate decreases from 2e3 to 1e5 "units?"
  % but they don't really tell us how. I'm guessing they mean
  % it slowed exponentially toward 1e5 rather than growing exponentially
  % toward it, and since we don't know the growth rate, I'll guess:
  tauw = 2e3 + (1e5-2e3)*(2./(1+exp(-packet/200))-1);
  %% Time within each packet
  while t<simdur
    t = dt*tind;
    tind = tind+1;
    % packet phase of each cell:
    phases = t - kappa*cellDists.*sin(theta - packetDir);

    % firing rate of each cell as a function of packet phase:
    f = ftonic + sin(phases);
    f(phases<0) = ftonic;
    f(phases>2*pi*wavePerPacket) = ftonic;
    % end this run if the packet has left the cells
    if all(phases>0) && all(all(f==ftonic))
    % update weight matrix
    % we use their approximation that mean(f) = ftonic
    W = W + dt/tauw*((f' - ftonic)*f - W);

    if livePlot>0 && (livePlot==1 || mod(tind,livePlot)==1)
        title(sprintf('Activity of sheet of cells; t = %.1f ms, packet %d',t,packet))
        axis square
        axis square
        title({'Weights from center cell onto all other cells'})

%% Optionally save weight matrix
if saveW