On 05-07-2013, at 17:37, Jannetta Steyn <janne...@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
______________________________________________
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.
