How to use XPP to get the FI curve for the HH equations

1. Run the ODE file, setting I0, the applied current to zero
2. Integrate the equations for a while (go into nUmerics and change
Total to 200 or so)
3. Click on InitCond-Last several times. We want to make sure we are
on an equilibrium point. (KB shortcut, I L I L .....). Do not skimp on
this. Do it ten times if you want!!
4. Steady state is V=-64.999..., m=.0529,h=.596,n=.317

5. Click on File Auto to get the AUTO window
All the rest of the commands are in the AUTO window

6. Click on Parameter and make sure Par1 is I0, the applied current
7. I am assuming here the units are mu A/mm^2 to set up the axes.
8. Click on Axes - Hi-Lo
    Y-axis: V
    Main Parm: I0
    ...
    Xmin: 0
    Ymin: -90
    Xmax: 1
    Ymax: 50

9. Numerics -
    Ntst: 60
    Nmax: 2000
    NPr: 500
    Ds: 0.02
    Dsmin: .000001
    Dsmax: 0.1
    Par Max: 1
    Par Min: 0

all othere as they are

10. Run-steady state

You will see a curve that is thick and then thin. The transition is
labeled 2.

11. Grab: hit the TAB key until you reach point #2 (it says HB at the
bottom of the window) - if you go past it, keep hitting TAB. When you
are on this point (HB), hit enter/return. You have grabbed thr HB
point and now AUTO will let you study solutions that are near this
point. The HB is a Hopf bifurctaion and periodic solutions will
emerge.

12.Click Run. It should say Hopf Pt.  Choose Periodic. A branch of
solutions will emerge wiggling around abit and then turning around and
going to the right. The max and min of the periodic solutions are
being plotted.

13. Click on Axes - Frequency and then accept the defaults
14. Click on Axes-Fit  and then Redraw
15. The FI curve appears - it is firing rate per msec.  The thin line
is stable solutions and the open curcles are unstable.

16. Click on File Postscript to save the plot.

-----------------------------------------
#  hh.ode the Dayan & Abbott dimensions
init v=-65  m=.05  h=0.6  n=.317
par i0=0
par vna=50  vk=-77  vl=-54.4  gna=1.20  gk=.36  gl=0.003  c=.01  phi=1
par ip=0 pon=50 poff=150
is(t)=ip*heav(t-pon)*heav(poff-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'=(I0+is(t) - 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
@ bound=10000,total=200,maxstor=10000
done