An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20130213/60ab15cd/attachment.pl>
An extended Hodgkin-Huxley model that doesn't want to work.
7 messages · Berend Hasselman, Jannetta Steyn
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
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20130215/7f0d5da9/attachment.pl>
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20130215/0591c71c/attachment.pl>
On 15-02-2013, at 10:28, Jannetta Steyn <jannetta at henning.org> wrote:
Hi Ben Thank you so much for spending the time to look at my code. I really appreciate it. The unnecessary parameters were artefacts from me working on the model and putting things in and taking things out. The model isn't quite producing what I expect yet, but it is definitely starting to look more interesting :-) I am now getting this error:
plot(ode)
Error in is.null(times) : 'times' is missing What on earth does it mean because I cant see any 'times' being missing. To me it seems that all is there?!?
Difficult to tell as you haven't shown how ode was defined. I just did this: out<-ode(y=init, times=times, func=HH_soma, parms=parms) plot(out) and got the plots of all state variables. Berend
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20130215/0e339918/attachment.pl>
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