HVCIter.m

function [w xdyn] = HVCIter(p)
% Runs one iteration of the simulation. p is a structure of parameters.
%
% Emily Mackevicius 1/14/2015, heavily copied from Hannah Payne's code
% which builds off Ila Fiete's model, with help from Michale Fee and Tatsuo
% Okubo. 

% redefining params that are used often outside the for loop
nsteps = p.nsteps;
n = p.n;
m = p.m;
w = p.w; 
wmax = p.wmax;
Wmax = wmax*m; 
eta = p.eta;
bdyn=p.input;

% initializing variables
xdyn=zeros(n,nsteps);
oldx = zeros(n,1);
oldy = zeros(n,1);

for t = 1:nsteps
 % Adaptation
 y = oldy + 1/p.tau*(-oldy+oldx);
 Aadapt = p.alpha * y; % adaptation
 
 % Net feedforward input. 
 B = bdyn(:,t); % external input
 AE = w*oldx; % excitatory input 
 AIff = p.beta * sum(oldx); % feed forward inhibition
 Anet = AE - AIff - Aadapt + B; % net feed forward input
 Anet(Anet < 0) = 0; % rectify
 
 % recurrent inhibition
 AIrec = p.gamma * sum(Anet); 
 
 % binary output
 x = (Anet - AIrec) > 0; 
 
 % STDP rule (Fiete et al 2010)
 dw_STDP = eta.*(x*(oldx)'-(oldx)*x');
 
 % Hetersynaptic penalty (Fiete et al 2010)
 dw_hLTDpre = eta*ones(n,1)*max(0, sum(w+dw_STDP,1)-Wmax); % Weights leaving cells (pre)
 dw_hLTDpost = eta*max(0, sum(w+dw_STDP,2)-Wmax)*ones(1,n); % Weights onto cells (post)
 
 % Update weights
 if eta>0
 dwtotal = dw_STDP-p.epsilon*(dw_hLTDpre+dw_hLTDpost);
 w = w + dwtotal;
 w(w > wmax) = wmax; % hard limit on strength of a single synapse
 w(w < 0) = 0; % weights cannot be negative
 w = w.*(~eye(p.n)); % clamp diagonal
 end
 
 oldx = double(x);
 oldy = y;
 xdyn(:,t)=x;
end