Aah. Thanks. I remember someone mentioned a named vector and I meant to ask
what that was, but I forgot. I have now fixed it.
On 5 July 2013 17:18, Berend Hasselman <b...@xs4all.nl> wrote:
>
> 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
>
>
>
>
>
>
>
--
===================================
Web site: http://www.jannetta.com
Email: janne...@henning.org
===================================
[[alternative HTML version deleted]]
______________________________________________
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.
