function [ a, b, c] = traub( t, x, flag, Vinp, gna, gk, ga, gca, gahp)

global	C gl ena ek el eca tca bconst
global    I Ina Ik Ia Ica Iahp

if isempty( flag)	% No flag: return xdot
   % state variables
   i=x(1); m=x(2); h=x(3); n=x(4); p=x(5); q=x(6); s=x(7); r=x(8); u=x(9); ca=x(10);
   % input current
   if size(Vinp,2) == 2	% (t,I) table
      v = interp1(Vinp(:,1),Vinp(:,2),t);
   else			% constant
      v = Vinp(1);
   end
   % activation and inactivation rates
   am = 0.32*(v+46.9)/(1-exp(-(v+46.9)/4));
   bm = 0.28*(v+19.9)/(exp((v+19.9)/5)-1);
   ah = 0.128*exp(-(v+43)/18);
   bh = 4/(1+exp(-(v+20)/5));
   an = 0.016*(v+24.9)/(1-exp(-(v+24.9)/5));
   bn = 0.25*exp(-(v+40)/40);
   ap = 0.02*(v+46.9)/(1-exp(-(v+46.9)/10));
   bp = 0.0175*(v+19.9)/(exp((v+19.9)/10)-1);
   aq = 0.016*exp(-(v+73.0)/18);
   bq = 0.05/(1+exp(-(v+49.9)/5));
   as = 0.912/(exp(-0.72*(v-5))+1);
   bs = 0.0114*(v+8.9)/(exp((v+8.9)/5)-1);
   ar = min(0.005,0.005*exp(-(v+60)/20));
   br = 0.005 - ar;
   uinf =(0.0005*ca)^2;
   tu = 0.0338/(min(0.00001*ca,0.01)+0.001);
   dm = am*(1-m)-bm*m;
   dh = ah*(1-h)-bh*h;
   dn = an*(1-n)-bn*n;
   dp = ap*(1-p)-bp*p;
   dq = aq*(1-q)-bq*q;
   ds = as*(1-s)-bs*s;
   dr = ar*(1-r)-br*r;
   du = (uinf-u)/tu;
   Ina = gna*m^2*h*(v-ena);
   Ik = gk*n*(v-ek);
   Ia = ga*p*q*(v-ek);
   Ica = gca*s^2*r*(v-eca);
   Iahp = gahp*u*(v-ek);
   % xdot
   a = [ gna*(v-ena)*(2*m*dm*h+m^2*dh)+gk*(v-ek)*dn+ga*(v-ek)*(dp*q+p*dq)+gca*(v-eca)*(2*s*ds*r+s^2*dr)+gahp*(v-ek)*du
      dm;
      dh;
      dn;
      dp;
      dq;
      ds;
      dr;
      du;
      -bconst*Ica-ca/tca];
   return;
end
switch( flag)
case 'init'			% Initialization
   % membrane capacitance (uF/cm^2)
   C=1;
   % maximum conductances (uS/cm^2)
   gl = 0.1;		% leak
   % reversal potentials (mV)
   ena = 60;		% sodium
   ek = -75;		% potassium
   eca = 135;      % calcium
   el = -60;		% leak
   % parameters of calcium dynamics
   bconst = 3;
   tca = 60;
   % return default values
   a = [0 100];	% default tspan
   b = [ 0; 0.5; 0.5; 0.5; 0.5; 0.5; 0.5; 0.5; 0.5; 1.0];	% default initial v,m,h,n
   c = [];			% default option made by odeset()
case 'name'			% Variable names
   a = char( 'i', 'm', 'h', 'n','p','q','s','r','u','ca');
otherwise
   error( [ ' invalid option: ', flag]);
end

% end of traub.m