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');