function [ts,ys] = stepnetwork(yinitial,alpha,tau)
%STEPNETWORK  Solve 6-unit network with step function response.

%    [TS,YS] = STEPNETWORK(YINITIAL,ALPHA,TAU) computes solution YS,
%    as a function of time TS, to the 
%    6-unit network (Edwards et al., 1999):
%
%         dy
%         -- = - y + w*g - tau,  
%         dt  
%
%   where w is a 6x6 matrix containing the 'synaptic efficacies' of the
%   units of the network, and the response function g(y) is a
%   Heaviside step function.
%  
%   Input parameters:
%    ALPHA = scalar 'synaptic connection' parameter; decreasing alpha
%             represents weakening of synaptic efficacies 
%            ( 0 <= alpha <= 1 )
%    TAU = vector of thresholds for the units in the network
%    YINITIAL = initial value of solution y = (y1,y2, ..., y6)
%
%   Reference: R Edwards, A Beuter and L Glass. (1999) Parkinsonian
%   tremor and simplification in network dynamics.  Bulletin of
%   Mathematical Biology, 61(1): 157-177.
%
%   See also SIGMNETWORK.
%

%   Copyright 2003
%   Centre for Nonlinear Dynamics in Physiology and Medicine

dt = 0.01;
lt = 5000;

n = 6;
ts = zeros(1,lt);
ys = zeros(n,lt);

if (nargout ~= 2),
  disp('ERROR: must have 2 output variables in the command syntax')
  disp('       eg. [ts,ys]=stepnetwork(y0,a,t)')
  break
end

if (nargin ~= 3),
  disp('ERROR: must have 3 input arguments in the command syntax')
  disp('       eg. [ts,ys]=stepnetwork(y0,a,t)')
  break
end

y = yinitial(:);
if (length(y) ~= 6),
  disp('ERROR: initial y must be a 6x1 vector')
  break
end

if (length(alpha) ~= 1),
  disp('ERROR: alpha must be a scalar')
  break
end

if ((alpha > 1.0) | (alpha < 0)),
  disp('ERROR: alpha must be between 0 and 1')
  break
end

tau = tau(:);
if (length(tau) ~= 6),
  disp('ERROR: tau must be a 6x1 vector')
  break
end

g = zeros(n,1);
w = zeros(n,n);

g(find(y>=0))=1;
w=-[ 0 1 0 0 0 1; 0 0 0 1 0 1; 0 0 0 1 alpha 0; 1 0 0 0 0 1; 
   1 1 0 0 0 0; 0 0 alpha 0 alpha 0];

tj=0;
while (tj < lt),
  a = w*g-tau;
  b = (a-y)./a;
%   compute first time a variable crosses 0
%   if no crossing, then compute solution to end of time interval
  tm = min(log(b(find(b>1)))); 
  ti = tj+1;
  if (isempty(tm))
    tj = lt;
  else
    tj = tj+ceil(tm/dt);
  end
%   compute solution up to crossing time
  ys(:,ti:tj) = a*ones(1,tj-ti+1) + (y-a)*exp(-(ti:tj)*dt)/exp(-(ti-1)*dt);
%   reset initial value of y and on/off function g
  y = ys(:,tj);
  g = zeros(n,1);
  g(find(y>=0)) = 1;
end

ts = dt:dt:lt*dt;
ys = ys(:,1:lt);

disp('Solution computed')

return





