Skip to content

A Hodgkin Huxley Model

8 messages · David Winsemius, Ben Bolker, Jannetta Steyn

Ben Bolker
# 
Jannetta Steyn <jannetta <at> henning.org> writes:
I'm not sure what "if v=30 then v<-c else u<-u+d" means. The first clause
is OK: I can imagine you want to implement a resetting event where
when v crosses 30, it gets reset to c: to do this part, see ?events in
the deSolve package (you are looking for "roots" in particular), and
the examples in ?lsodar.  The "else" clause doesn't make much sense to
me, as v will essentially *always* (except at particular, non-generic
time points) be not-equal to 30 ... it would probably make sense to me
if I went back to read the original paper, but that's too much work
...

(By the way, the semicolons at the ends of lines are unnecessary;
they're harmless, but usually the mark of a recovering MATLAB
user ...)
David Winsemius
#
On Feb 10, 2013, at 4:22 AM, Jannetta Steyn wrote:

            

      
      
      
    

        
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

            

      
      
      
    

        

      
      
    

        
David Winsemius, MD
Alameda, CA, USA

  
  


  


      

    

    

      
      

        


  

        

      David Winsemius
    

    

      
      #
    

  


  

      
On Feb 10, 2013, at 8:09 AM, Ben Bolker wrote:

        

            

      
      
      
    

        
I was interested enough to already have done that work. The else  
clause is incorrect:

-------------- next part --------------
A non-text attachment was scrubbed...
Name: Izhikevich.png
Type: image/png
Size: 47689 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20130210/eb7ca7cc/attachment.png>
-------------- next part --------------

            

      
      
      
    

        

      
      
    

        
David Winsemius, MD
Alameda, CA, USA

  
  


  


      

    

    

      
      

        


  

        

      Ben Bolker
    

    

      
      #
    

  


  

      
On 13-02-10 01:22 PM, David Winsemius wrote:

            

      
      
      
    

        
[ image probably didn't make it to r-help, but I saw it in direct e-mail ]


  So that should be "if (v>=30) { v <- c; u <- u+d }" which should
certainly be solvable in the framework I suggested above.

  
  


  


      

    

    

      
      

        
      

    

    

      
      

        


  

        

      David Winsemius
    

    

      
      #
    

  


  

      
On Feb 10, 2013, at 3:35 PM, Jannetta Steyn wrote:

        

            

      
      
      
    

        
The output certainly looks more sensible. You commented out the if- 
reset strategy and are using events. Can you explain how the event  
code is triggered? I no longer see a definition of a v>=30 trigger.  
The vignette says "Events occur when the values of state variables are  
instantaneously changed. They can be specified as a data.frame, or in  
a function. Events can also be triggered by a root function." Is the  
code detecting the sudden rise to infinite values and using that as  
the trigger for a reset?