% 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
% 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];
% we're just using 1D trajectories here:
% compute activity at each position in trajectory
% load trajectory from disk:
% 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
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
title('Activity vs. position')
xlabel('1D Position (cm)')