function deriv = hhfield(t,state)

% Evaluate the Hodgkin-Huxley vector field.


% parameters:
C = 1;     % uF
E_Na = 50; % mV
E_K = -77;
E_leak = -54.387;

g_Na = 120.;
g_K = 36.;
g_leak = 0.3;

I = 10.;


% Grab state variables.
v = state(1);
m = state(2);
h = state(3);
n = state(4);


% The equations themselves.  If you're going to add an ion
% channel, make sure to update the deriv line below.

dv = ( -g_Na * m^3 * h * ( v - E_Na ) ...
       -g_K * n^4 * ( v - E_K ) ...
       -g_leak * ( v - E_leak ) ...
       + I) / C;

dm = a_m(v) * ( 1-m ) - b_m(v) * m;
dh = a_h(v) * ( 1-h ) - b_h(v) * h;
dn = a_n(v) * ( 1-n ) - b_n(v) * n;


% Make sure the order here matches the "Grab state
% variables" part above.
deriv = [dv ; dm ; dh ; dn];



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Note the shift by 65 mV is necessary because the version of
%the HH equations in Gerstner-Kistler follows the Hodgkin
%and Huxley, i.e. the resting membrane voltage was set to 0.

function y=a_m(v)
v = v + 65;
y=(2.5 - 0.1*v) / ( exp(2.5 - 0.1*v) -1);

function y=b_m(v)
v = v + 65;
y=4*exp(-v/18);

function y=a_h(v)
v = v + 65;
y=0.07*exp(-v/20);

function y=b_h(v)
v = v + 65;
y=1/(exp(3-0.1*v)+1);

function y=a_n(v)
v = v + 65;
y=( 0.1 - 0.01*v ) / ( exp(1-0.1*v) - 1.);

function y=b_n(v)
v = v + 65;
y = 0.125 * exp( -v/80.0);
