TITLE na3
: Na current 
: modified from Jeff Magee. M.Migliore may97
: added sh to account for higher threshold M.Migliore, Apr.2002
: MM: in parallel to updated na.mod from Mainen & Sejnowski,1996
: trap0() was replaced with efun(), "which is a better approximation
: in the vicinity of a singularity"
: MM, mtaufac added to explore the effect of faster activation

NEURON {
    SUFFIX na3
    USEION na READ ena WRITE ina
    RANGE  gbar, ar, sh, shx, mtaufac, htaufac, qa, ina3, thegna, m, h, s, m3h,qinf, thinf
    GLOBAL minf, hinf, mtau, htau, sinf, taus
}

PARAMETER {
    sh   = 0    (mV)        : leftward shift of minf and mtau
    shx  = 0    (mV)        : rightward shift of hinf and htau
    gbar = 0.010    (mho/cm2)
    mtaufac = 1     (1)     : factor that is multiplied with the expression for mtau    
    htaufac = 1      (1)    : factor that is multiplied with the expression for htau 
                                
    tha  =  -30 (mV)        : v 1/2 for act 
    qa   = 7.2  (mV)        : act slope (4.5)       
    Ra   = 0.4  (/ms)       : open (v)      
    Rb   = 0.124    (/ms)       : close (v)     

    thi1  = -45 (mV)        : v 1/2 for inact   
    thi2  = -45     (mV)        : v 1/2 for inact   
    qd   = 1.5  (mV)            : inact tau slope
    qg   = 1.5      (mV)
    mmin=0.02   
    hmin=0.5            
    q10=2
    Rg   = 0.01     (/ms)       : inact recov (v)   
    Rd   = .03  (/ms)       : inact (v) 
    qq   = 10        (mV)
    tq   = -55      (mV)

    thinf  = -62    (mV)        : inact inf V_1/2 - originally -50 (changed to -62)
    qinf  =  6.9    (mV)        : inact inf slope - originally 4 (changed to 6.9)

        vhalfs=-60  (mV)        : slow inact.
        a0s=0.0003  (ms)        : a0s=b0s
        zetas=12    (1)
        gms=0.2     (1)
        smax=10     (ms)
        vvh=-58     (mV) 
        vvs=2       (mV)
        ar=1        (1)     : 1=no inact., 0=max inact.
    ena     (mV)            : must be explicitly def. in hoc
    celsius
    v       (mV)
}


UNITS {
    (mA) = (milliamp)
    (mV) = (millivolt)
    (pS) = (picosiemens)
    (um) = (micron)
} 

ASSIGNED {
    ina         (mA/cm2)
    ina3        (mA/cm2)
    thegna      (mho/cm2)
    minf        hinf 
    m3h     
    mtau (ms)   htau (ms)   
    sinf (ms)   taus (ms)
}
 

STATE { m h s}

BREAKPOINT {
        SOLVE states METHOD cnexp
        thegna = gbar*m*m*m*h*s
        m3h = m*m*m*h*s
    ina = thegna * (v - ena)
    ina3 = thegna * (v - ena)
    ina=ina3
} 

INITIAL {
    trates(v,ar,sh, shx, mtaufac, htaufac, qa, qinf)
    m=minf  
    h=hinf
    s=sinf
}


FUNCTION alpv(v(mV)) {
         alpv = 1/(1+exp((v-vvh-sh)/vvs))
}
        
FUNCTION alps(v(mV)) {  
  alps = exp(1.e-3*zetas*(v-vhalfs-sh)*9.648e4/(8.315*(273.16+celsius)))
}

FUNCTION bets(v(mV)) {
  bets = exp(1.e-3*zetas*gms*(v-vhalfs-sh)*9.648e4/(8.315*(273.16+celsius)))
}

LOCAL mexp, hexp, sexp

DERIVATIVE states {   
        trates(v,ar,sh, shx, mtaufac, htaufac, qa, qinf)      
        m' = (minf-m)/mtau
        h' = (hinf-h)/htau
        s' = (sinf - s)/taus
}

: efun() is a better approx than trap0 in vicinity of singularity--

PROCEDURE trates(vm,a2,sh2, sh3, taufac, taufac2, qa, qinf) {  
        LOCAL  a, b, c, qt
        qt=q10^((celsius-24)/10)
:   a = trap0(vm,tha+sh2,Ra,qa)
    a = Ra * qa * efun((tha+sh2 - vm)/qa)
:   b = trap0(-vm,-tha-sh2,Rb,qa)
    b = Rb * qa * efun((vm - tha-sh2)/qa)

:    mtau = taufac*0.001
:   mtau = 1/(a+b)/qt
    mtau = taufac/(a+b)/qt
        if (mtau<mmin) {mtau=taufac*mmin}
:        if (mtau<mmin) {mtau=mmin} 
    minf = a/(a+b)

:   a = trap0(vm,thi1+sh2,Rd,qd)
    a = Rd * qd * efun((thi1+sh3 - vm)/qd)
:   b = trap0(-vm,-thi2-sh2,Rg,qg)
    b = Rg * qg * efun((vm - thi2-sh3)/qg)
:   htau =  1/(a+b)/qt
    htau =  taufac2/(a+b)/qt
:        if (htau<hmin) {htau=hmin}
        if (htau<hmin) {htau=taufac2*hmin}
    hinf = 1/(1+exp((vm-thinf-sh3)/qinf))
    c=alpv(vm)
        sinf = c+a2*(1-c)
        taus = bets(vm)/(a0s*(1+alps(vm)))
        if (taus<smax) {taus=smax}
}


COMMENT
FUNCTION trap0(v,th,a,q) {
    if (fabs(v-th) > 1e-6) {
            trap0 = a * (v - th) / (1 - exp(-(v - th)/q))
    } else {
            trap0 = a * q
    }
}   
ENDCOMMENT

FUNCTION efun(z) {
    if (fabs(z) < 1e-6) {
        efun = 1 - z/2
    }else{
        efun = z/(exp(z) - 1)
    }
}