# hhh.ode init v=0 m=0 h=0 n=0 par vna=50 vk=-77 vl=-54.4 gna=120 gk=36 gl=0.3 c=1 i=10 am(v)=.1*(v+40)/(1-exp(-(v+40)/10)) bm(v)=4*exp(-(v+65)/18) ah(v)=.07*exp(-(v+65)/20) bh(v)=1/(1+exp(-(v+35)/10)) an(v)=.01*(v+55)/(1-exp(-(v+55)/10)) bn(v)=.125*exp(-(v+65)/80) v'=(I - gna*h*(v-vna)*m^3-gk*(v-vk)*n^4-gl*(v-vl))/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 doneand the iterative version using the exponential method:
# this is the exponential form for HH equations hhexp.dif init v=0 m=0 h=0 n=0 par vna=50 vk=-77 vl=-54.4 gna=120 gk=36 gl=0.3 c=1 i=10 par delt=.5 am(v)=.1*(v+40)/(1-exp(-(v+40)/10)) bm(v)=4*exp(-(v+65)/18) ah(v)=.07*exp(-(v+65)/20) bh(v)=1/(1+exp(-(v+35)/10)) an(v)=.01*(v+55)/(1-exp(-(v+55)/10)) bn(v)=.125*exp(-(v+65)/80) # x'=-a*x+b ==> x(t+delt)=exp(-a*delt)*(x(t)-b/a)+b/a bv=(I+vna*gna*h*m^3+vk*gk*n^4+gl*vl)/c av=(gna*h*m^3+gk*n^4+gl)/c v'=exp(-av*delt)*(v-bv/av)+bv/av amv=am(v)+bm(v) bmv=am(v)/amv m'=exp(-amv*delt)*(m-bmv)+bmv ahv=ah(v)+bh(v) bhv=ah(v)/ahv h'=exp(-ahv*delt)*(h-bhv)+bhv anv=an(v)+bn(v) bnv=an(v)/anv n'=exp(-anv*delt)*(n-bnv)+bnv doneIn the exponential Euler version, switch the method to Discrete from the numerics menu. Iterate for 500 time steps, which at the step size of 0.5 milliseconds represents (500)(.5)=250. The oscillation occurs at the 85th iterate or about 42 milliseconds. Now repeat this with the true differential equation using a variety of integration methods (Backward Euler, Euler, modified euler, Runge-Kutta, Gear) with different time steps. Compare the results. Try exponential euler by setting the parameter delt=1.
This is a severely abridged version of Sherman's notes. The full version is available as a postscript file by clicking here.