Hi Berend
This is the code. Pretty much just changed to what you suggested which is
CaI=1, removing the unnecessary variables and using deSolve:
rm(list = ls())
library(deSolve)
## Hodkin-Huxley model
HH_soma <- function(times, init, parms) {
with(as.list(c(init, parms)),{
# Na only used in Axon
#Naminf <-1/(1+exp(-(v+24.7)/5.29));
#Nataum <- function(v) 1.32 - (1.26/(1+exp(-(v+120)/25)));
#Nahinf <-1/(1+exp((v+489)/5.18));
#Natauh <-(0.67/(1+exp(-(v+62.9)/10))) * (1.5+(1/(1+exp((v+34.9)/36))));
#PD
# mca10
CaTminf <- function(v) 1/(1+exp(-(v+25)/7.2));
# hca10
CaThinf <- function(v) 1/(1+exp(v+36)/7);
# taumca1
CaTtaum <- function(v) 55- (49.5/(1+exp(-v+58)/17));
# tauhca1
CaTtauh <- function(v) 350 - (300/(1+exp(-v+50)/16.9));
#mca20
CaSminf <- function(v) 1/(1+exp(-(v+22)/8.5));
#taumca2
CaStaum <- function(v) 16-(13.1/(1+exp(-(v+25.1)/26.4)));
# mna0
Napminf <- function(v) 1/(1+exp(-(v+26.8)/8.2));
# hna0
Naphinf <- function(v) 1/(1+exp(-(v+48.5)/5.18));
taumna <- function(v) 19.8-(10.7/(1+exp(-(v+26.5)/8.6)));
tauhna <- function(v) 666-(379/(1+exp(-(v+33.6)/11.7)));
# mh0
hminf <- function(v) 1/(1+exp(v+70)/6);
# taumh
htaum <- function(v) 272+(1499/(1+exp(-(v+42.2)/8.73)));
Kminf <- function(v) 1/(1+exp(-(v+14.2)/11.8));
Ktaum <- function(v) 7.2-(6.4/(1+exp(-(v+28.3)/19.2)));
# Reversal potential of intracellular calcium concentration
# Nernst Equation using extracellular concentration of Ca = 13mM
# eca
ECa <- function(Ca2) 12.2396*log(13000/(Ca2));
#ECa <- function(CaI) 12.2396*log(13000/(CaI));
#Sum of all the Ca
# function(v) CaTminf(v) + CaSminf(v);
CaI <- gCaT*CaTm^3*CaTh*(v-ECa(CaI)) + gCaS*CaSm^3*(v-ECa(CaI))
#AB
#dCa2 <- (((-F*Caminf(v))-Caminf(v) + C0)/TauCa)
dCa2 <- (((-F*CaI) - Ca2 + C0)/TauCa)
# mk20
KCaminf <- function(v, Ca2) (Ca2/(Ca2+30))*(1/(1+exp(-(v+51)/8)));
# taumk
KCataum <- function(v) 90.3 - ((75.09/(1+exp(-(v+46)/22.7))));
#AB
Aminf <- function(v) 1/(1+exp(-(v+27)/8.7));
Ahinf <- function(v) 1/(1+exp((v+56.9)/4.9));
Ataum <- function(v) 11.6-(10.4/(1+exp(-(v+32.9)/15.2)));
Atauh <- function(v) 38.6-(29.2*(1+exp(-(v+38.9)/26.5)));
#proc
#mp0
procminf <- function(v) 1/(1+exp((v+56.9)/4));
#taump
proctaum <- function(v) 0.5;
dv <- (-1*(I
+ CaI
+ gNap*Napm^3*Naph*(v-ENap)
+ gh*hm*(v-Eh)
+ gK*Km^4*(v-EK)
+ gKCa * KCam^4*(v-EKCa)
+ gA*Am^4*Ah*(v-EA)
+ gL*(v-EL))
/ C);
dCaTm <- (CaTminf(v) - CaTm)/CaTtaum(v);
dCaTh <- (CaThinf(v) - CaTh)/CaTtauh(v);
dCaSm <- (CaSminf(v) - CaSm)/CaStaum(v);
dNapm <- (Napminf(v) - Napm)/taumna(v);
dNaph <- (Napminf(v) - Naph)/tauhna(v);
dhm <- (hminf(v) - hm)/htaum(v);
dKm <- (Kminf(v) - Km)/Ktaum(v);
dKCam <- (KCaminf(v, Ca2) - KCam)/KCataum(v);
dAm <- (Aminf(v) - Am)/Ataum(v);
dAh <- (Ahinf(v) - Ah)/Atauh(v);
list(c(dv,
dCaTm, dCaTh,
dCaSm,
dNapm, dNaph,
dhm,
dKm,
dKCam,
dAm, dAh))
})
}
parms = c(gCaT=22.5, gCaS=60, gNap=4.38, gh=0.219,
gK=1576.8, gKCa=251.85, gA=39.42, gL=0.105,
ENap=50, Ca2=0.52,
Eh=-20, EK=-80, EL=-55, EKCa=-80, EA=-80,
C=1/12, I=6.5, F=0.418, TauCa=303, C0=0.5, CaI=1);
time = seq(from=0, to=1000, by=0.1);
init = c(v=-65, CaTm=0.52 , CaTh=0.52, CaSm=0.52,
Napm=0.52, Naph=0.52, hm=0.52, Km=0.52,
KCam=0.52, Am=0.52, Ah=0.52);
out<-ode(y=init, times=time, func=HH_soma, parms=parms);
plot(ode)
o<-data.frame(out);
plot(o$time, o$v, type='l');
Regards
Jannetta
On 15 February 2013 09:46, Berend Hasselman <[email protected]> wrote:
>
> On 15-02-2013, at 10:28, Jannetta Steyn <[email protected]> wrote:
>
> > Hi Ben
> >
> > Thank you so much for spending the time to look at my code. I really
> > appreciate it.
> >
> > The unnecessary parameters were artefacts from me working on the model
> and
> > putting things in and taking things out.
> >
> > The model isn't quite producing what I expect yet, but it is definitely
> > starting to look more interesting :-)
> >
> > I am now getting this error:
> >
> >> plot(ode)
> > Error in is.null(times) : 'times' is missing
> >
> > What on earth does it mean because I cant see any 'times' being missing.
> To
> > me it seems that all is there?!?
> >
>
> Difficult to tell as you haven't shown how ode was defined.
>
> I just did this:
>
> out<-ode(y=init, times=times, func=HH_soma, parms=parms)
> plot(out)
>
> and got the plots of all state variables.
>
> Berend
>
>
>
--
===================================
Web site: http://www.jannetta.com
Email: [email protected]
===================================
[[alternative HTML version deleted]]
______________________________________________
[email protected] 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.
