R and MatLab implementations of the same model differs
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