[R] A Hodgkin Huxley Model
Ben Bolker
bbolker at gmail.com
Sun Feb 10 17:09:55 CET 2013
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 ...)
More information about the R-help
mailing list