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: janne...@henning.org
===================================
[[alternative HTML version deleted]]
______________________________________________
R-help@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
______________________________________________
R-help@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.
