[R] A Hodgkin Huxley Model
David Winsemius
dwinsemius at comcast.net
Sun Feb 10 19:01:37 CET 2013
On Feb 10, 2013, at 4:22 AM, Jannetta Steyn wrote:
> 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)
Since the function errors out after some warnings, the plot statement
has nothing to work with. Leaving it out allowed one to look at
traceback() which I found rather uninformative. I was rather surprised
not to see u and v being returned from the function. Instead you are
returning dv and du which I would have expected _not_ to be returned
but rather to be calculated at the next step.
( I do get more interesting and error-free results but still get a
warning by following this path. Further reading of the vignette in the
package suggest the returned values are not being used in succeeding
cycle, I think if you are implicitly creating discontinuities, that
you should not be returning gradients.)
Caveat: I am not a mathematician , instead a physician with very
limited training in diffyQ.
Putting in a debugging line lets you follow the "progress" of the
process".
> 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; cat(paste("u= ",u," :: v=", v,"\n"))
+ 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)
u= 0.2 :: v= 8.2
u= 0.2 :: v= 8.2
u= 0.199516713662881 :: v= 8.19951671366288
u= 0.199516648247331 :: v= 8.19951664824733
u= 0.199033296573405 :: v= 8.1990332965734
u= 0.199033231197198 :: v= 8.1990332311972
u= 0.186610529982878 :: v= 8.18661052998288
u= 0.186610696611869 :: v= 8.18661069661187
snipped
There is an initial phase until something (that I don't really
understand) suddenly pushes v to -65
u= -7.50652247221155 :: v= 0.493477527788452
u= -7.50366977501533 :: v= 0.496330224984674
u= -7.50366976442979 :: v= 0.496330235570206
u= -7.49992636486762 :: v= 0.500073635132377
u= -7.49992633949909 :: v= 0.500073660500912
u= -7.49589685117969 :: v= -65
u= -7.49589682789333 :: v= -65
u= -7.49154790150432 :: v= -65
u= -7.49154790376302 :: v= -65
u= -7.48684008116797 :: v= -65
snipped
Then it evolves to a point where some warnings are thrown but
continues for quite a while after that:
u= -4.5758982855538 :: v= -65
u= -4.57755515069294 :: v= -65
u= -4.57755496668311 :: v= -65
DLSODA- Warning..Internal T (=R1) and H (=R2) are
such that in the machine, T + H = T on the next step
(H = step size). Solver will continue anyway.
In above message, R =
[1] 4.511249e+01 3.336971e-15
u= -4.56929471441551 :: v= -65
u= -4.56929462175284 :: v= -65
DLSODA- Warning..Internal T (=R1) and H (=R2) are
such that in the machine, T + H = T on the next step
(H = step size). Solver will continue anyway.
In above message, R =
[1] 4.511249e+01 3.336971e-15
u= -4.56029029422924 :: v= -65
u= -4.56029003859891 :: v= -65
DLSODA- Warning..Internal T (=R1) and H (=R2) are
such that in the machine, T + H = T on the next step
(H = step size). Solver will continue anyway.
In above message, R =
[1] 4.511249e+01 3.336971e-15
u= -4.55039409155327 :: v= -65
u= -4.55039364611721 :: v= -65
u= -4.5524591470321 :: v= -65
u= -4.55245893270376 :: v= -65
DLSODA- Warning..Internal T (=R1) and H (=R2) are
such that in the machine, T + H = T on the next step
(H = step size). Solver will continue anyway.
In above message, R =
[1] 4.511249e+01 2.667571e-15
u= -4.54396212536601 :: v= -65
u= -4.54396204287127 :: v= -65
DLSODA- Warning..Internal T (=R1) and H (=R2) are
such that in the machine, T + H = T on the next step
(H = step size). Solver will continue anyway.
In above message, R =
[1] 4.511249e+01 2.667571e-15
u= -4.53467596517973 :: v= -65
u= -4.5346756233244 :: v= -65
u= -4.53633248846354 :: v= -65
u= -4.53633230445371 :: v= -65
DLSODA- Warning..Internal T (=R1) and H (=R2) are
such that in the machine, T + H = T on the next step
(H = step size). Solver will continue anyway.
In above message, R =
[1] 4.511249e+01 2.209650e-15
u= -4.52807205218611 :: v= -65
u= -4.52807195952345 :: v= -65
DLSODA- Warning..Internal T (=R1) and H (=R2) are
such that in the machine, T + H = T on the next step
(H = step size). Solver will continue anyway.
In above message, R =
[1] 4.511249e+01 2.209650e-15
u= -4.51906763199984 :: v= -65
u= -4.51906737636951 :: v= -65
DLSODA- Warning..Internal T (=R1) and H (=R2) are
such that in the machine, T + H = T on the next step
(H = step size). Solver will continue anyway.
In above message, R =
[1] 4.511249e+01 2.209650e-15
u= -4.50917142932387 :: v= -65
u= -4.50917098388781 :: v= -65
u= -4.51123648480269 :: v= -65
u= -4.51123627047435 :: v= -65
DLSODA- Warning..Internal T (=R1) and H (=R2) are
such that in the machine, T + H = T on the next step
(H = step size). Solver will continue anyway.
In above message, R =
[1] 4.511249e+01 1.766392e-15
u= -4.5027394631366 :: v= -65
u= -4.50273938064186 :: v= -65
DLSODA- Warning..Internal T (=R1) and H (=R2) are
such that in the machine, T + H = T on the next step
(H = step size). Solver will continue anyway.
In above message, R =
[1] 4.511249e+01 1.766392e-15
DLSODA- Above warning has been issued I1 times.
It will not be issued again for this problem.
In above message, I =
[1] 10
u= -4.49345330295032 :: v= -65
u= -4.49345296109499 :: v= -65
u= -4.49510982623412 :: v= -65
u= -4.49510964222429 :: v= -65
snipped
With u rising steadily and v solidly at -65
And a final "terminal phase where u suddenly becomes infinite and
throws an error:
u= 27.4860465093796 :: v= -65
u= 27.4953325870711 :: v= -65
u= 27.4953329289265 :: v= -65
u= 27.4936760637873 :: v= -65
u= 27.4936762477971 :: v= -65
u= 27.5019365000647 :: v= -65
u= 27.5019365927274 :: v= -65
u= 27.510940920251 :: v= -65
u= 27.5109411758814 :: v= -65
u= Inf :: v= -65
Error in if (v >= 30) v <- c else v <- u + d :
missing value where TRUE/FALSE needed
> plot(out)
>
> Can someone perhaps help me out?
>
> Regards
> Jannetta
>
> --
>
> ===================================
> Web site: http://www.jannetta.com
> Email: jannetta at henning.org
> ===================================
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
David Winsemius, MD
Alameda, CA, USA
More information about the R-help
mailing list