[R] A Hodgkin Huxley Model

[Ben Bolker]

On 13-02-10 01:22 PM, David Winsemius wrote:
> 
> 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:
>

 [ image probably didn't make it to r-help, but I saw it in direct e-mail ]


  So that should be "if (v>=30) { v <- c; u <- u+d }" which should
certainly be solvable in the framework I suggested above.