# [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 =
 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 =
 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 =
 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 =
 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 =
 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 =
 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 =
 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 =
 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 =
 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 =
 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 =
 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