function [ a, b, c] = hh( t, x, flag, I)
% HH: dynamics of Hodgkin-Huxley model
% August 1999 by K. Doya

global	C gna gk gl ena ek el

if isempty( flag)	% No flag: return xdot
	% state variables
	v=x(1); m=x(2); h=x(3); n=x(4);
	% input current
	if size(I,2) == 2	% (t,I) table
		It = interp1(I(:,1),I(:,2),t);
	else			% constant
		It = I(1);
	end
	% activation and inactivation rates
	am = 0.1*(v+40)/(1-exp(-(v+40)/10));
	bm = 4*exp(-(v+65)/18);
	ah = 0.07*exp(-(v+65)/20);
	bh = 1/(1+exp(-(v+35)/10));
	an = 0.01*(v+55)/(1-exp(-(v+55)/10));
	bn = 0.125*exp(-(v+65)/80);
	% xdot
	a = [ (It - gna*m^3*h*(v-ena) - gk*n^4*(v-ek) - gl*(v-el))/C;
		am*(1-m)-bm*m;
		ah*(1-h)-bh*h;
		an*(1-n)-bn*n];
	return;
end

switch( flag)
case 'init'			% Initialization
	% membrane capacitance (uF/cm^2)
	C=1;
	% maximum conductances (uS/cm^2)
	gna = 120;		% sodium
	gk = 36;		% potassium
	gl = 0.3;		% leak
	% reversal potentials (mV)
	ena = 50;		% sodium
	ek = -77;		% potassium
	el = -54.4;		% leak
	% return default values
	a = [0 100];	% default tspan
	b = [ -65; 0.052; 0.596; 0.317];	% default initial v,m,h,n
	c = [];			% default option made by odeset()
case 'name'			% Variable names
	a = char( 'v', 'm', 'h', 'n');
otherwise
	error( [ ' invalid option: ', flag]);
end

% end of hh.m