[R] An extended Hodgkin-Huxley model that doesn't want to work.
Berend Hasselman
bhh at xs4all.nl
Thu Feb 14 18:41:16 CET 2013
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
More information about the R-help
mailing list