[R] R and MatLab implementations of the same model differs

[Berend Hasselman]

On 05-07-2013, at 17:37, Jannetta Steyn <jannetta at henning.org> wrote:

> Thank you very much to all of you for your help. This model now works as it
> should (I believe). This is the final code:
> 
> rm(list=ls())
> 
> library(deSolve)
> 
> ST <-  function(time, init, parms) {
>  with(as.list(c(init, parms)),{
> 
>    #functions to calculate activation m and inactivation h of the currents
>    #Axon
>    mNax <- function(v) 1/(1+exp(-(v+24.7)/5.29))
>    taumNa <- function(v) 1.32-(1.26/(1+exp(-(v+120)/25)))
>    hNax <- function(v) 1/(1+exp((v+48.9)/5.18))
>    tauhNa <- function(v)
> 0.67/(1+exp(-(v+62.9)/10))*(1.5+(1/(1+exp((v+34.9)/3.6))))
>    mKx <- function(v) 1/(1+exp(-(v+14.2)/11.8))
>    taumK <- function(v) 7.2-(6.4/(1+exp(-(v+28.3)/19.2)))
> 
> 
>    # Currents as product of maximal conducatance(g), activation(m) and
> inactivation(h)
>    # Driving force (v-E) where E is the reversal potential of the
> particular ion
> 
>    # AB axon
>    iNa_axon_AB <-
> gNa_axon_AB*mNa_axon_AB^3*hNa_axon_AB*(v_axon_AB-ENa_axon_AB)
>    iK_axon_AB <- gK_axon_AB*mK_axon_AB^4*(v_axon_AB-EK_axon_AB)
>    iLeak_axon_AB <- gLeak_axon_AB*(v_axon_AB-ELeak_axon_AB)
> 
>    # AB Axon equations
>    dv_axon_AB <- (0 -iNa_axon_AB -iK_axon_AB -iLeak_axon_AB)/C_axon_AB
>    dmNa_axon_AB <- (mNax(v_axon_AB)-mNa_axon_AB)/taumNa(v_axon_AB)
>    dhNa_axon_AB <- (hNax(v_axon_AB)-hNa_axon_AB)/tauhNa(v_axon_AB)
>    dmK_axon_AB <- (mKx(v_axon_AB)-mK_axon_AB)/taumK(v_axon_AB)
> 
>    list(c(dv_axon_AB,dmNa_axon_AB, dhNa_axon_AB, dmK_axon_AB
>    ))
> 
>  })}
> ## Set initial state
> init = c(v_axon_AB=-55,mNa_axon_AB=0,hNa_axon_AB=1,mK_axon_AB=0)
> ## Set parameters
> # AB
> 
> gNa_axon_AB=300e-3
> gK_axon_AB=52.5e-3
> gLeak_axon_AB=0.0018e-3
> 
> # AB Axon Reversal potentials
> ENa_axon_AB=50
> EK_axon_AB=-80
> ELeak_axon_AB=-60
> 
> C_axon_AB=1.5e-3
> 
> 
> I=0
> 
> parms =
> c(ENa_axon_AB,EK_axon_AB,ELeak_axon_AB,gNa_axon_AB,gK_axon_AB,gLeak_axon_AB,C_axon_AB)
> ## Set integrations times
> times = seq(from=0, to=5000, by = 0.05);
> 
> out<-ode(y=init, times=times, func=ST, parms=parms, method="ode45")
> 
> o<-data.frame(out)
> 
> plot(o$v_axon_AB~o$time,type='l')
> 

You must make the parms vector a named vector because the with(.. in your function ST is not using the global values for those entities available in the global environment.

You can check this by inserting this line after the line with times =  and before the line out <- ode(…

rm(list=c("ENa_axon_AB","EK_axon_AB","ELeak_axon_AB","gNa_axon_AB","gK_axon_AB","gLeak_axon_AB","C_axon_AB"))

You will get an error.

So create the parms vector like this

parms =
c(ENa_axon_AB=ENa_axon_AB,EK_axon_AB=EK_axon_AB,ELeak_axon_AB=ELeak_axon_AB,gNa_axon_AB=gNa_axon_AB,
    gK_axon_AB=gK_axon_AB,gLeak_axon_AB=gLeak_axon_AB,C_axon_AB=C_axon_AB)

or better

parms <- c(ENa_axon_AB=50,EK_axon_AB=-80,ELeak_axon_AB=-60,gNa_axon_AB=300e-3,
                  gK_axon_AB=52.5e-3,gLeak_axon_AB=0.0018e-3,C_axon_AB=1.5e-3)

and remove the expressions with the separate assignments to the variables.


And use <- !!

Berend