FuhsTouretzky2006_development.m

% 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')
    load(fname);
    fprintf('+ Loaded pre-generated W.\n')
  else
    W = zeros(ncells);
  end
end

%% 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])
  drawnow
end
  
%% !! 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))
      break
    end
    
    % 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)
        figure(h);
        subplot(121);
        image(32*reshape(f,n,n));
        title(sprintf('Activity of sheet of cells; t = %.1f ms, packet %d',t,packet))
        axis square
        set(gca,'ydir','normal')
        subplot(122);
        imagesc(reshape(W(:,ncells/2-round(n/2)),n,n));
        axis square
        set(gca,'ydir','normal')
        title({'Weights from center cell onto all other cells'})
        drawnow
    end
  end
end

%% Optionally save weight matrix
if saveW
  save(fname,'W','-v7.3');
end