function [tn,yn] = sigmnetwork(yinitial,alpha,tau,beta)
%SIGMNETWORK  Solve 6-unit network with sigmoidal response.

%    [TN,YN] = SIGMNETWORK(YINITIAL,ALPHA,TAU,BETA) computes solution
%    YN, as a function of time TN, 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 sigmoid
%   function of the form 
%
%                  1 + tanh(beta*y)
%         g(y) = -------------------
%                        2
%
%   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)
%    BETA = controls slope (gain) of the response function
%
%   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 DERIVS, RUKA4, PLOTRESPONSE, STEPNETWORK.
%
 
%   Copyright 2003
%   Centre for Nonlinear Dynamics in Physiology and Medicine
 
global betav tauv w

h = 0.01;
nt = 5000;

n = 6;
tn = zeros(1,nt);
yn = zeros(n,nt);

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

if (nargin ~= 4),
  disp('ERROR: must have 4 input arguments in the command syntax')
  disp('       eg. [ts,ys]=stepnetwork(y0,a,t,b)')
  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(:);
tauv = tau;
if (length(tauv) ~= 6),
  disp('ERROR: tau must be a 6x1 vector')
  break
end

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

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];

gg = 0.5*(1+tanh(beta*y));

yn(:,1) = y;

disp('Computing solution ...')

for i = 2:nt+1,
  t = (i-1)*h;
  yn(:,i) = ruka4(['derivs'],y,h,t);
  y = yn(:,i);
end

tn = h:h:nt*h;
yn = yn(:,1:nt);

disp('Solution computed')

return

