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 <bhh@xs4all.nl> wrote:

>
> Forgot Reply to All.
>
> On 13-02-2013, at 23:30, Jannetta Steyn <jannetta@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
>
> ______________________________________________
> R-help@r-project.org 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: jannetta@henning.org
===================================

	[[alternative HTML version deleted]]

