#  hhh.ode 
init v=-65  m=.05  h=0.6  n=.317  
par vna=50  vk=-77  vl=-54.4  gna=120  gk=36  gl=0.3  c=1  phi=1 i0=0
par ip=0 ton=50 toff=150
i(t)=i0+ip*heav(t-ton)*heav(toff-t)
am(v)=phi*.1*(v+40)/(1-exp(-(v+40)/10))
bm(v)=phi*4*exp(-(v+65)/18)
ah(v)=phi*.07*exp(-(v+65)/20)
bh(v)=phi*1/(1+exp(-(v+35)/10))
an(v)=phi*.01*(v+55)/(1-exp(-(v+55)/10))
bn(v)=phi*.125*exp(-(v+65)/80)
v'=(I(t) - gna*h*(v-vna)*m^3-gk*(v-vk)*n^4-gl*(v-vl)-gsyn*s*(v-vsyn))/c
m'=am(v)*(1-m)-bm(v)*m
h'=ah(v)*(1-h)-bh(v)*h
n'=an(v)*(1-n)-bn(v)*n
s'=sinf(v)*(1-s)-s/tausyn
# track the currents 
sinf(v)=alpha/(1+exp(-v/vshp))
par alpha=2,vshp=5,tausyn=20,gsyn=0,vsyn=0
aux ina=gna*(v-vna)*h*m^3
aux ik=gk*(v-vk)*n^4
aux il=gl*(v-vl)
# track the stimulus
aux stim=i(t)
@ meth=qualrk,total=200,dt=.25,xlo=0,xhi=200,ylo=-75,yhi=30,yp=v,xp=t
@ bound=10000
done