# Hodgkin-Huxley Model
#
dV/dt=(-gna*(m^3)*h*(V-Ena)-gk*(n^4)*(V-Ek)-gl*(V-El)+Io*f(t))/Cm
#
dm/dt=alpha_m(V)*(1-m)-beta_m(V)*m
dh/dt=alpha_h(V)*(1-h)-beta_h(V)*h
dn/dt=alpha_n(V)*(1-n)-beta_n(V)*n
#
# activation and inactivation functions
alpha_m(V)=-am*(V+40.0)/(exp(-(V+40)/10.0)-1)
beta_m(V)=20.0*exp(-(V+65.0)/18.0)
alpha_h(V)=0.07*exp(-(V+67.0)/100.0)
beta_h(V)=1.0/(exp(-(V+37.0)/10.0)+1)
alpha_n(V)=-0.01*(V+55.0)/(exp(-(V+55.0)/10.0)-1)
beta_n(V)=bn*exp(-(V+65.0)/80.0)
#
f(t)=heav(t-t_on)*heav(t_off-t)
#
# parameter values
par gna=120 gk=36 gl=0.3
par Ena=50 Ek=-77.5 El=-54.3
par Cm=1
par am=0.5  bn=0.125
par Io=0 t_on=5 t_off=100
#
# initial conditions
V(0)=-65.2152 
m(0)=0.0516
h(0)=0.5504  
n(0)=0.3144
#
aux Iapp=Io*f(t)
#
@ total=100,method=stiff,bound=500,dt=.1,dtmin=1e-12
@ ylo=-80,yhi=50,xlo=0,xhi=100,maxstor=10000
@ bounds=10000
done