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')
My MatLab expert came to my rescue and both model now produce similar
results :-)
Regards
Jannetta
On 5 July 2013 12:11, Jannetta Steyn <janne...@henning.org> wrote:
>
>
>
> On 5 July 2013 09:44, Berend Hasselman <b...@xs4all.nl> wrote:
>
>>
>> On 05-07-2013, at 09:53, Jannetta Steyn <janne...@henning.org> wrote:
>>
>> >
>> >
>> > >
>> > > I don't quite know how to explain the "doesn't work" in more detail
>> without
>> > > any visual aid.
>> >
>> > You said that R got into an indefinite loop, whatever that maybe.
>> >
>> >
>> > > When I change the solver to ode45 the script never stops running. I
>> have even left it over night at one stage but the next day it was still
>> busy. Is there a way to see what it is doing and to determine why it seems
>> to be in an
>> > > "infinite loop"?
>> >
>> > The script just ran using ode45!! For the first time ever.
>> >
>>
>> Please keep the conversation on R-help.
>> Don't reply to me personally.
>>
>
>
> Apologies for that. It was meant to go to the list.
>
>
>>
>> Which script?
>> With the corrections suggested by me and Peter? (gK_axon_AB=52.5-3 ==>
>> gK_axon_AB=52.5e-3)
>>
>>
> The R script. With the correction of AB=52.5e-3
>
>
>
>> For what it's worth: lsdoda seems quicker.
>>
>
> Yes, lsoda is quicker but I want to use ode45 because that is what they
> used in the paper I am referencing and I am trying to get the same results
> as they did.
>
>
>> Variable v_acon_AB now converges to -60 (the equilibrium state?)
>>
>
> I'm trying to determine now whether this is correct.
>
> Regards
> Jannetta
>
>
> --
>
> ===================================
> Web site: http://www.jannetta.com
> Email: janne...@henning.org
> ===================================
>
--
===================================
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.
