Hi All Having been pointed the use of events and roots in deSolve, I was able to implement the Izchikevich model of spiking neurons. However, I'm not too sure of defining the event. The deSolve documentation says:
An event is triggered when the ball hits the ground (height = 0) Then velocity (y2) is reversed and reduced by 10 percent. The root function, y[1] = 0, triggers the event: > root <- function(t, y, parms) y[1] Firstly I couldn't see where y[1] became 0, but I implemented Izchikevich as follows: 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=2, I=10); times=seq(from=1, to=1000, by=0.1); init=c(v=-65, u=0.2); root <- function(time, init, parms) { return(init[1]-30) } event <- function(time, init, parms) { with(as.list(c(init, parms)), { init[2] <- init[1] + d init[1] <- c return(init) }) } out<-ode(y=init, times=times, func=Izhikevich, parms=parms, events=list(func=event, root=TRUE), rootfun=root) plot(out) The reasoning behind my implementation was that if y[1] is 0 without being set then init[1] will be 0 without being set. I need for the event to trigger when init[1] is >= 30. Setting the root to init[1]+30 did not produce the desired result, making it init[1]-30, produced something better but I'm not sure that it is right yet. The actual equations are: v' = 0.04v^2 + 5v + 140 - u + I u- = a(bv-u) If v=30mV then v-, u-u+d Any help would be appreciated. Kind 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.