[R] An extended Hodgkin-Huxley model that doesn't want to work.

[Berend Hasselman]

On 15-02-2013, at 11:22, Jannetta Steyn <jannetta at henning.org> 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