A Hodgkin Huxley Model

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
...

(By the way, the semicolons at the ends of lines are unnecessary;
they're harmless, but usually the mark of a recovering MATLAB
user ...)