TITLE


INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}
UNITS {
	(mA) = (milliamp)
	(mV) = (millivolt)
	(pS) = (picosiemens)
	(um) = (micron)
} 

NEURON {
	SUFFIX snahh
	USEION na READ ena WRITE ina
	RANGE m, h, gna, gbar
        RANGE wnoise_rv,N,Ntot,D,gamma,alpham,betam,alphah,betah
	RANGE minf, hinf, mtau, htau,am,bm,ah,bh
	GLOBAL q10, temp, tadj, vmin, vmax
}

PARAMETER {
	gbar = 1200   	(pS/um2)	: 0.12 mho/cm2
	temp = 6.3	(degC)		: original temp 
	q10  = 3			: temperature sensitivity

	v 		(mV)
	dt		(ms)
	celsius		(degC)
	vmin = -120	(mV)
	vmax = 100	(mV)

        gamma  =  20      (pS)    
        D     =  60      (1/um2)
}

ASSIGNED {
        am		(/ms)
	bm		(/ms)
        ah		(/ms)
	bh		(/ms)

	ina 		(mA/cm2)
	gna		(pS/um2)
	ena		(mV)
	minf 		
        hinf
	mtau (ms)	 
        htau (ms)
	tadj

        N 
        Ntot
        scale_dens (pS/um2) 
        alpham		
	betam		
        alphah		
	betah		
        wnoise_rv
	wnmu 	                    : This is a nasty trick to pass variables into the verbatim code
	wnsig

}
 
STATE { m h 
        m0h0 m0h1 m1h0 m1h1 m2h0 m2h1 m3h0 m3h1
m0h0_m1h0  m1h0_m2h0  m2h0_m3h0  m0h1_m1h1  m1h1_m2h1  m2h1_m3h1  
m3h0_m2h0  m2h0_m1h0  m1h0_m0h0  m3h1_m2h1  m2h1_m1h1  m1h1_m0h1  
m0h0_m0h1 m0h1_m0h0 m1h0_m1h1 m1h1_m1h0 m2h0_m2h1 m2h1_m2h0 m3h0_m3h1 m3h1_m3h0 
}

INITIAL { 
	trates(v)
	m = minf
	h = hinf
        scale_dens = gamma/area
        N   = floor(D*area+0.99) 

        m0h0 = floor((1-m)*(1-m)*(1-m)*(1-h)*N)
        m1h0 = floor(3*m*(1-m)*(1-m)*(1-h)*N)
        m2h0 = floor(3*m*m*(1-m)*(1-h)*N)
        m3h0 = floor(m*m*m*(1-h)*N)

        m0h1 = floor((1-m)*(1-m)*(1-m)*h*N)
        m1h1 = floor(3*m*(1-m)*(1-m)*h*N)
        m2h1 = floor(3*m*m*(1-m)*h*N)
        m3h1 = floor(m*m*m*h*N)
              
  
m0h0_m1h0=0 
m1h0_m2h0=0 
m2h0_m3h0=0 
m0h1_m1h1=0
m1h1_m2h1=0
m2h1_m3h1=0 
m3h0_m2h0=0 
m2h0_m1h0=0
m1h0_m0h0=0 
m3h1_m2h1=0 
m2h1_m1h1=0 
m1h1_m0h1=0

m0h0_m0h1=0 
m0h1_m0h0=0 
m1h0_m1h1=0 
m1h1_m1h0=0 
m2h0_m2h1=0 
m2h1_m2h0=0 
m3h0_m3h1=0 
m3h1_m3h0=0

}

BREAKPOINT {
        SOLVE states
        gna = strap(m3h1)*scale_dens 
	ina = (1e-4) * gna * (v - ena)
} 

LOCAL mexp, hexp 

PROCEDURE states() {   :Computes state variables m, h, and n 
        trates(v)      :             at the current v and dt.
        m = m + mexp*(minf-m)
        h = h + hexp*(hinf-h)

        alpham = strap(am*dt)
        betam  = strap(bm*dt)

        m0h0_m1h0  = binom(3.0*alpham,m0h0)
        m1h0_m2h0  = binom(2.0*alpham,m1h0)
        m2h0_m3h0  = binom(alpham,m2h0)

        m0h1_m1h1  = binom(3.0*alpham,m0h1)
        m1h1_m2h1  = binom(2.0*alpham,m1h1)
        m2h1_m3h1  = binom(alpham,m2h1)

        m3h0_m2h0  = binom(3.0*betam,m3h0)
        m2h0_m1h0  = binom(2.0*betam,m2h0)
        m1h0_m0h0  = binom(betam,m1h0)

        m3h1_m2h1  = binom(3.0*betam,m3h1)
        m2h1_m1h1  = binom(2.0*betam,m2h1)
        m1h1_m0h1  = binom(betam,m1h1)

	m0h0 = m0h0  - m0h0_m1h0 + m1h0_m0h0
  
	m1h0 = m1h0  - m1h0_m2h0 - m1h0_m0h0  + m2h0_m1h0 + m0h0_m1h0
        m2h0 = m2h0  - m2h0_m3h0 - m2h0_m1h0  + m3h0_m2h0 + m1h0_m2h0
        m3h0 = m3h0  - m3h0_m2h0 + m2h0_m3h0

	m0h1 = m0h1  - m0h1_m1h1 + m1h1_m0h1 
        m1h1 = m1h1  - m1h1_m2h1 - m1h1_m0h1 + m2h1_m1h1 + m0h1_m1h1
        m2h1 = m2h1  - m2h1_m3h1 - m2h1_m1h1 + m3h1_m2h1 + m1h1_m2h1
        m3h1 = m3h1  - m3h1_m2h1 + m2h1_m3h1

        alphah = strap(ah*dt)
        betah  = strap(bh*dt)

        m0h0_m0h1 = binom(alphah,m0h0)
        m0h1_m0h0 = binom(betah,m0h1)
        m1h0_m1h1 = binom(alphah,m1h0)
        m1h1_m1h0 = binom(betah,m1h1)
        m2h0_m2h1 = binom(alphah,m2h0)
        m2h1_m2h0 = binom(betah,m2h1)
        m3h0_m3h1 = binom(alphah,m3h0)
        m3h1_m3h0 = binom(betah,m3h1)

        m0h0 = m0h0 - m0h0_m0h1  + m0h1_m0h0  
        m1h0 = m1h0 - m1h0_m1h1  + m1h1_m1h0 
        m2h0 = m2h0 - m2h0_m2h1  + m2h1_m2h0 
        m3h0 = m3h0 - m3h0_m3h1  + m3h1_m3h0 

        m0h1 = m0h1 - m0h1_m0h0  + m0h0_m0h1 
        m1h1 = m1h1 - m1h1_m1h0  + m1h0_m1h1 
        m2h1 = m2h1 - m2h1_m2h0  + m2h0_m2h1 
        m3h1 = m3h1 - m3h1_m3h0  + m3h0_m3h1 

        VERBATIM
        return 0;
        ENDVERBATIM
}

PROCEDURE trates(v) {  
        LOCAL tinc
        TABLE minf, mexp, hinf, hexp,am,bm,ah,bh
	DEPEND dt, celsius, temp
	
	FROM vmin TO vmax WITH 199

	rates(v): not consistently executed from here if usetable == 1

        tadj = q10^((celsius - temp)/10)
        tinc = -dt * tadj

        mexp = 1 - exp(tinc/mtau)
        hexp = 1 - exp(tinc/htau)
}


PROCEDURE rates(vm) {  
        am = .1 * vtrap(-(vm+40),10)
        bm =   4 * exp(-(vm+65)/18)
	mtau = 1/(am+bm)*tadj
	minf = am*mtau/tadj
        : if you change to the dm//dt format you can say
 	: mdelta = a*(1-m) - b*m + wnoise(1,1)  
		:"h" inactivation 

        ah = .07 * exp(-(vm+65)/20)
        bh =  1 / (exp(-(vm+35)/10) + 1)
	htau = 1/(ah+bh)*tadj
	hinf = ah*htau/tadj
}


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
 	}
}	

FUNCTION strap(x) {  :Traps for negaiv values for number of channels
if (x < 0) {
	strap = 0
	VERBATIM
	fprintf (stderr,"snam95.mod:strap: negative value for state");
	ENDVERBATIM
}else{
	strap = x
        }
}


FUNCTION vtrap(x,y) {  :Traps for 0 in denominator of rate eqns.
        if (fabs(x/y) < 1e-6) {
                vtrap = y*(1 - x/y/2)
        }else{
                vtrap = x/(exp(x/y) - 1)
        }
}


FUNCTION wnoise(mu,sig)
{
VERBATIM
  extern double gasdev();
ENDVERBATIM
wnmu = mu                      : This is a nasty trick to pass variables into the verbatim code
wnsig = sig
  if (wnsig <=0) {
        binom = 0
     } else {
VERBATIM
       wnoise_rv = wnmu+gasdev()*wnsig;   
ENDVERBATIM
  wnoise= wnoise_rv
  }
}
UNITSON