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: [email protected]
===================================
[[alternative HTML version deleted]]
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
