R help - Feb 2013 - An extended Hodgkin-Huxley model that doesn't want to work.

Hi All

I have been struggling with this model for some time now and I just can't
get it to work correctly. The messages I get when running the code is:

DLSODA-  Warning..Internal T (=R1) and H (=R2) are
      such that in the machine, T + H = T on the next step
     (H = step size). Solver will continue anyway.
In above message, R [1] 0 0
DINTDY-  T (=R1) illegal
In above message, R [1] 0.1
      T not in interval TCUR - HU (= R1) to TCUR (=R2)
In above message, R [1] 0 0
DINTDY-  T (=R1) illegal
In above message, R [1] 0.2
      T not in interval TCUR - HU (= R1) to TCUR (=R2)
In above message, R [1] 0 0
DLSODA-  Trouble in DINTDY.  ITASK = I1, TOUT = R1
In above message, I [1] 1
In above message, R [1] 0.2
Error in lsoda(y, times, func, parms, ...) :
  illegal input detected before taking any integration steps - see written
message



I'll first paste the formulae and then I'll paste my code. If anyone can
spot something wrong with my implementation it would really make my day.

(1)
dV/dt = (I_ext - I_int-I_coup)/C
I_ext = injected current
I_int = Sum of all ion currents
I_coup = coupling current (but we're not using it here )

(2)
I_i = g_i * m_i^pi * h_i^pi(V-E)
i identifies the ion, thus I_K would be Potassium current.

(3)
dm/dt = (m_inf*V - m)/tau_m

(4)
dh/dt = (h_inf*V-h)/tau_h

(5)
The Nernst equation is used to calculate reversal potential for Ca:
Eca = 12.2396 * log(13000/Ca2+)

(6)
d[Ca_2+]/dt = (F*I_Ca - [Ca2+] + C0)/Tau_Ca


tau_m, tau_h, m_inf and h_inf are all calculated according to formulae
provided in a paper. In my code these are calculated for the different
channels into the following variables:

CaTminf, CaThinf, CaTtaum, CaTtauh, CaSminf, CaStaum, Napminf, Naphinf,
taumna, tauhna, hminf, htaum, Kminf and Ktaum

The E (reversal potential) values for all the channels are given, except
for CaT and CaS which uses Eca as calculated in (5).

Current for Ca is calculated by summing the CaT and CaS currents, hence
CaI = gCaT*CaTm^3*CaTh*(v-Eca(v)) + gCaS*CaSm^3(v-ECa(v)


Here is the code:

library(simecol)
## Hodkin-Huxley model
HH_soma <-  function(time, 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,
           dCa2,
           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=10, F=0.418, TauCa=303, C0=0.5, CaI=0);
times = seq(from=0, to=400, 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, ECa=-80);


out<-ode(y=init, times=times, func=HH_soma, parms=parms);
o<-data.frame(out);
plot(o$time, o$v, type='l');

Please ask if any further information is required.

Many thanks
Jannetta



-- 

==================================Web site: http://www.jannetta.com
Email: jannetta@henning.org
==================================
	[[alternative HTML version deleted]]

Forgot Reply to All.

On 13-02-2013, at 23:30, Jannetta Steyn <jannetta at henning.org> wrote:
> Hi All
> 
> I have been struggling with this model for some time now and I just
can't
> get it to work correctly. The messages I get when running the code is:
> 
> DLSODA-  Warning..Internal T (=R1) and H (=R2) are
>     such that in the machine, T + H = T on the next step
>    (H = step size). Solver will continue anyway.
> In above message, R > [1] 0 0
> DINTDY-  T (=R1) illegal
> In above message, R > [1] 0.1
>     T not in interval TCUR - HU (= R1) to TCUR (=R2)
> In above message, R > [1] 0 0
> DINTDY-  T (=R1) illegal
> In above message, R > [1] 0.2
>     T not in interval TCUR - HU (= R1) to TCUR (=R2)
> In above message, R > [1] 0 0
> DLSODA-  Trouble in DINTDY.  ITASK = I1, TOUT = R1
> In above message, I > [1] 1
> In above message, R > [1] 0.2
> Error in lsoda(y, times, func, parms, ...) :
> illegal input detected before taking any integration steps - see written
> message
> 
> 
> 
> I'll first paste the formulae and then I'll paste my code. If
anyone can
> spot something wrong with my implementation it would really make my day.
> 
> (1)
> dV/dt = (I_ext - I_int-I_coup)/C
> I_ext = injected current
> I_int = Sum of all ion currents
> I_coup = coupling current (but we're not using it here )
> 
> (2)
> I_i = g_i * m_i^pi * h_i^pi(V-E)
> i identifies the ion, thus I_K would be Potassium current.
> 
> (3)
> dm/dt = (m_inf*V - m)/tau_m
> 
> (4)
> dh/dt = (h_inf*V-h)/tau_h
> 
> (5)
> The Nernst equation is used to calculate reversal potential for Ca:
> Eca = 12.2396 * log(13000/Ca2+)
> 
> (6)
> d[Ca_2+]/dt = (F*I_Ca - [Ca2+] + C0)/Tau_Ca
> 
> 
> tau_m, tau_h, m_inf and h_inf are all calculated according to formulae
> provided in a paper. In my code these are calculated for the different
> channels into the following variables:
> 
> CaTminf, CaThinf, CaTtaum, CaTtauh, CaSminf, CaStaum, Napminf, Naphinf,
> taumna, tauhna, hminf, htaum, Kminf and Ktaum
> 
> The E (reversal potential) values for all the channels are given, except
> for CaT and CaS which uses Eca as calculated in (5).
> 
> Current for Ca is calculated by summing the CaT and CaS currents, hence
> CaI = gCaT*CaTm^3*CaTh*(v-Eca(v)) + gCaS*CaSm^3(v-ECa(v)
> 
> 
> Here is the code:
> 
> library(simecol)
> ## Hodkin-Huxley model
> HH_soma <-  function(time, 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,
>          dCa2,
>          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=10, F=0.418, TauCa=303, C0=0.5, CaI=0);
> times = seq(from=0, to=400, 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, ECa=-80);
> 
> 
> out<-ode(y=init, times=times, func=HH_soma, parms=parms);
> o<-data.frame(out);
> plot(o$time, o$v, type='l');
> 
1. why are you suing simecol when only deSolve is necessary?

2. there are some errors in your model. The return list from function HH_soma
does not correspond with the state variables.
 You return dCa2 but there is no state variable Ca2. ECa is in the initial state
variable init but there is no dECa derivative in the return list of function
HH_SOMA.  Finally you have set parameter CaI to 0 in the parms vector.
But that parameter is only used in calling function ECa with argument CaI, where
you are now dividing by 0 and taking the logarithm of the result (which is Inf).

So change the return list of function HH_soma to

  list(c(dv,
         dCaTm, dCaTh,
         dCaSm,
         dNapm, dNaph,
         dhm,
         dKm,
         dKCam,
         # dCa2,    <==== Not needed not in init
         dAm, dAh))  

Change parms into (change CaI to 1; if that is a sensible value if for you to
decide).

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=10, F=0.418, TauCa=303, C0=0.5, CaI=1)

and change the initial state to (leaving out ECa)

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)#, ECa=-80);


Finally  you don't need to end R expressions with ; (only needed to join
several expressions on a single line)

I uses library(deSolve) and got an interesting plot.

Berend