BlairEtAl2007.m

% Blair, Welday, and Zhang 2007's Moire oscillatory interference model
% eric zilli - 20110907 - v1.0
%
% This model produce the entorhinal grid pattern by assuming
% that elsewhere in the brain, grid patterns already exist but at
% much higher spatial frequencies (much smaller field spacings). The
% model then suggest that the interference pattern produced by two such
% "theta grids" can produce a larger scale grid. The theta grids in this
% implementation have a fixed relationship to position, so the activity
% needn't be updated in a manner to produce path integration: the grid cell
% activity at each location along a trajectory can all be calculated
% directly in one step.
%
% Because of this, this script does not perform path integration but rather
% plots either the total spatial pattern or the firing that
% would occur along the specified trajectory (which is just a subset of the
% overall pattern). The output is a grid cell "firing rate" that is the
% thresholded sum of the activity of the input theta grids. The manuscript
% uses two passes of a convolution filter to smooth/blur the fields, but
% they do not specify which filter was used, so we use a simple box-car
% filter.
%
% Since the model already assumes the grids exist, perhaps the primary
% contribution of the manuscript is the demonstration of how, if one
% assumes one of the theta grids is the generator of LFP theta, theta
% phase precession may or may not occur depending on how the interfering
% minigrids are chosen (i.e. the scale vs. rotation rules).
%
% 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.

% Note: to generate the BlairEtAl2007 component of my Figure 1, set:
% useRealTrajectory = 0; useLengthScalingRule = 1; alpha = 0.1;
% nSpatialBins = 1000; minx = -3; maxx = 18; miny = -3; maxy = 18; % cm


% if =0, plot whole 2D pattern.
% if =1, load trajectory from disk and plot activity only along that
useRealTrajectory = 0;

% if =0, uses rotation scaling rule (fixed phase)
% if =1, uses length scaling rule (precessing)
useLengthScalingRule = 1;

%% Model parameters
mu_M = 4; % threshold for cell to fire
c = [0; 0]; % center of one grid field; ignored right now
theta = 0; % orientation of hexagonal pattern
% for length scaling rule:
alpha = 0.05; % difference in scale of theta grids, dorsoventral values: 0.1429 >= lalpha >= 0.0667
% for rotation scaling rule:
phii = 1; % one-half difference in orientation of two theta grids, degrees
lambda = 5; % baseline spacing of theta grids, cm
omegaDirs = theta + [0 pi/3 2*pi/3]; % orientation of component gratings
omega = lambda*[cos(omegaDirs); sin(omegaDirs)]';

%% Gain function g parameters
a = 0.3;
b = -3/2;
g = @(x)(exp(a*(x-b))-1); % gain function

%% Activation function for single grid
% r is a 2-by-n vector of n coordinates
G = @(r)(g(sum(cos(omega*r))));

%% Microgrid rotation matrix for product scaling rule
R = @(phi)([cos(phi) sin(phi); -sin(phi) cos(phi)]);

%% Firing field plot variables
nSpatialBins = 1000;
minx = -100; maxx = 100; % cm
miny = -100; maxy = 100; % cm

%% Compute spatial firing
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:
  r = [pos(1,:); pos(2,:)];
  clear pos;
else
  % compute activity at every position in environment
  nSpatialBins = 1000;
  minx = 0; maxx = 100; % cm
  miny = 0; maxy = 100; % cm
  % This is the "environment" to use for my Figure 1:
%   nSpatialBins = 1000;
%   minx = -3; maxx = 18; % cm
%   miny = -3; maxy = 18; % cm
  xs = linspace(minx,maxx,nSpatialBins);
  ys = linspace(miny,maxy,nSpatialBins);
  [X,Y] = meshgrid(xs,ys);
  xs = reshape(X,1,[]);
  ys = reshape(Y,1,[]);
  r = [xs; ys];
  clear X Y;
end

if useLengthScalingRule
  M = G(r) + G(r*(1+alpha)) - mu_M;
else
  M = G(R(-phii)*r) + G(R(phii)*r) - mu_M;
end

% Thresholding:
M = M.*(M>0);

if useRealTrajectory
  figure; plot(r(1,:),r(2,:));
  hold on;
  plot(r(1,M>0), r(2,M>0),'r.')
  title('Spikes (red) and trajectory (blue)')
else
  M = reshape(M,nSpatialBins,[]);
  figure; imagesc(M); title('Raw firing fields');
  figure; imagesc(conv2(M,ones(6)/36)); title('Boxcar-smoothed firing fields');
  
  % No thresholding so fields smoothly fade out:
%   mM = 64*ones(size(M,1),size(M,2),3);
%   scalefact = 64/G([0; 0]);
%   mM(:,:,1) = mM(:,:,1)-scalefact*reshape(G(r*(1+alpha)),nSpatialBins,[]);
%   mM(:,:,2) = mM(:,:,2)-scalefact*reshape(G(r),nSpatialBins,[]);
%   mM(:,:,3) = mM(:,:,3)-scalefact*reshape(G(r),nSpatialBins,[]);
%   figure; image(mM/64); title('Microgrids (red and green)')
  
  % Threshold colors in plot to make fields clearer:
  mM = 64*ones(size(M,1),size(M,2),3);
  scalefact = 64;
  mM(:,:,1) = mM(:,:,1)-scalefact*(reshape(G(r*(1+alpha))>1,nSpatialBins,[]));
  mM(:,:,2) = mM(:,:,2)-scalefact*(reshape(G(r)>1,nSpatialBins,[]));
  mM(:,:,3) = mM(:,:,3)-scalefact*(reshape(G(r)>1,nSpatialBins,[]));
  figure; image(mM/64); title('Microgrids (red and green)')

end