GaussierEtAl2007.m

% 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));
end

%% 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));
         interp1(pos(3,:),pos(2,:),0:dt:pos(3,stopAt));
         interp1(pos(3,:),pos(3,:),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
else
  'oops not done yet'
end

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

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;
  end

  %% 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)
      continue
    end
    
    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);
  end
  
  % 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];
    end
    spikeTimes(spikeind(1)) = t;
    spikeCoords(spikeind(1),:) = [x(tind) y(tind)];
    spikeind(1) = spikeind(1)+1;
  end
  if O(watchCell+1,tind)>spikeThreshold
    if mod(spikeind(2),1000)==1
      spikeCoords2(spikeind(2)+1000,:) = [0 0];
    end
    spikeCoords2(spikeind(2),:) = [x(tind) y(tind)];
    spikeind(2) = spikeind(2)+1;
  end
  if O(watchCell+2,tind)>spikeThreshold
    if mod(spikeind(3),1000)==1
      spikeCoords3(spikeind(3)+1000,:) = [0 0];
    end
    spikeCoords3(spikeind(3),:) = [x(tind) y(tind)];
    spikeind(3) = spikeind(3)+1;
  end
  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);
  end
  
  
  if livePlot>0 && (livePlot==1 || mod(tind,livePlot)==1)
      figure(figh);
      subplot(131);
      imagesc(O)
      title('Activity');
      xlabel('Time (s)')
      ylabel('Grid cell #')
      axis square
      set(gca,'ydir','normal')
      title(sprintf('t = %.1f s',t))
      
      subplot(132);
      imagesc(Wmo);
      axis square
      set(gca,'ydir','normal')
      title({'Weights from modulo pop.','to grid cells O'});
      xlabel('Modulo cell #')
      ylabel('Grid cell #')

      subplot(133);
      plot(pos(1,1:tind),pos(2,1:tind))
      hold on;
      if ~isempty(spikeCoords2)
        plot(spikeCoords2(1:spikeind(2)-1,1),spikeCoords2(1:spikeind(2)-1,2),'g.')
      end
      if ~isempty(spikeCoords3)
        plot(spikeCoords3(1:spikeind(3)-1,1),spikeCoords3(1:spikeind(3)-1,2),'m.')
      end
      if ~isempty(spikeCoords)
        plot(spikeCoords(1:spikeind(1)-1,1),spikeCoords(1:spikeind(1)-1,2),'r.')
      end
      axis square
      title({'Trajectory (blue) and','spikes of 3 grid cells'})
      drawnow
      
  end  

end