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');


Regards
Jannetta


On 15 February 2013 09:46, Berend Hasselman <bhh@xs4all.nl> wrote:

>
> On 15-02-2013, at 10:28, Jannetta Steyn <jannetta@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
>
>
>


-- 

===================================
Web site: http://www.jannetta.com
Email: jannetta@henning.org
===================================

	[[alternative HTML version deleted]]

