A Hodgkin Huxley Model

On Feb 10, 2013, at 8:09 AM, Ben Bolker wrote:

Jannetta Steyn <jannetta <at> henning.org> writes:

Hi All

It has been suggested to me that the folks in this list might be able to
help me out of my misery.

As part of my learning curve to building a rather detailed model of a
small
neurone complex I am implementing some existing models in R.

For the moment I want to implement the Izhikevich model as described
in his
2003 paper. The equations are as follows:

v' = 0.04v^2 + 5v + 140 - u - I
u' = a(bv - u)
if v=30 then v<-c else u<-u+d

I don't know how to implement the if condition in my R code. I'm using
deSolve. Here is my code so far.

library(deSolve);
Izhikevich <- function(time, init, parms) {
 with(as.list(c(init, parms)),{
   dv <- (0.04*v^2)+(5*v)+140-u+I;
   du <- a*(b*v-u);
   if (v>=30) v<-c else v<-u+d;
   list( c(dv, du))
 })}
parms=c( a=0.02, b=0.2, c=-65, d=8,I=10);
times=seq(from=1, to=1000, by=1);
init=c(v=-65, u=0.2);

out<-ode(y=init, times=times, func=Izhikevich, parms=parms)
plot(out)

Can someone perhaps help me out?

Regards
Jannetta

I'm not sure what "if v=30 then v<-c else u<-u+d" means. The first clause
is OK: I can imagine you want to implement a resetting event where
when v crosses 30, it gets reset to c: to do this part, see ?events in
the deSolve package (you are looking for "roots" in particular), and
the examples in ?lsodar.  The "else" clause doesn't make much sense to
me, as v will essentially *always* (except at particular, non-generic
time points) be not-equal to 30 ... it would probably make sense to me
if I went back to read the original paper, but that's too much work
...

I was interested enough to already have done that work. The else clause
is incorrect: