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] =================================== [[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.
