[R] A Hodgkin Huxley Model

[Ben Bolker]

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