On 15-02-2013, at 11:22, Jannetta Steyn <[email protected]> wrote:
> 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');
You are not doing what I told you to do.
You should plot the result of ode not ode itself (that's the function).
So do
plot(out)
and you will get the plots.
Berend
______________________________________________
[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.
