Sorry that was supposed to be Berend and not Ben
On 15 February 2013 09:28, Jannetta Steyn <[email protected]> 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?!? > > Kind Regards > Jannetta > > > On 14 February 2013 17:41, Berend Hasselman <[email protected]> wrote: > >> >> Forgot Reply to All. >> >> On 13-02-2013, at 23:30, Jannetta Steyn <[email protected]> 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 >> >> ______________________________________________ >> [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. >> > > > > -- > > =================================== > Web site: http://www.jannetta.com > Email: [email protected] > =================================== > -- =================================== 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.
