R help - Feb 2013 - using and event in deSolve

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

Hi Jannetta,

as far as I can see, your implementation was almost ;-) correct, except 
that:

init[2] <- init[1] + d

should be:

init[2] <- init[2] + d



The root function with:

return(init[1]-30)

is correct, because this triggers an event when init[1]-30 crosses the 
zero line, i.e. when init[1] == 30. This can also be seen at the figure 
if we decrease the external step size, see the corrected example below.

Note also, that ";" at the end of line is not required (and should be 
avoided), because R is not C.

You can also send such questions to the special interest group:

r-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models


Thanks for using deSolve

Thomas



-------------- next part --------------
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)
    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)
times <- seq(from = 1, to = 100, by = 0.005)
init  <- c(v = -65, u = 0.2)

root <- function(time, init, parms) {
  return(v = init[1] - 30)
}

event <- function(time, init, parms) {
  with(as.list(parms), {
    init[2] <- init[2] + 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)