COMMENT

skhh.mod

K+ channel, Hodgkin-Huxley style kinetics. Stochastic version 

ENDCOMMENT

INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}

NEURON {
	SUFFIX skhh
	USEION k READ ek WRITE ik
	RANGE n, gk, gbar
        RANGE binom_rv,N,D,gamma,alpha,beta
	RANGE ninf, ntau,a,b
	RANGE ninf,ntau
	GLOBAL q10, temp, tadj, vmin, vmax
}

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

PARAMETER {
	gbar = 360   	(pS/um2)	: 0.036 mho/cm2
	ek  = -77.5	(mV)		 	
	temp = 6.3	(degC)		: original temp 
	q10  = 3			: temprature temperature sensitivity
                                        :  
	v 		(mV)
	dt		(ms)
	celsius		(degC)
	vmin = -120	(mV)
	vmax = 100	(mV)
        D        = 18   (1/um2)
        gamma    = 20   (pS)
}

ASSIGNED {
	ik 		(mA/cm2)
	gk		(pS/um2)
	ninf 	   
        ntau
	tadj

	a		(/ms)
	b		(/ms)
        N 
        scale_dens (pS/um2) 
        alpha		
	beta		
        binp 
        binn   
        binom_rv
}
 

STATE { n
        N0 N1 N2 N3 N4      : N states populations
        n0_n1 n1_n0 : number of channels moving from one state to the other 
        n1_n2 n2_n1 
        n2_n3 n3_n2 
        n3_n4 n4_n3 
}

INITIAL { 
	trates(v)
	n = ninf
        N   = floor(D*area+0.99) 
        N4  = floor((n*n*n*n)*N) 
        N3  = floor(4*(n*n*n)*(1-n)*N) 
        N2  = floor(6*(n*n)*((1-n)*(1-n))*N) 
        N1  = floor(4*n*(1-n)*(1-n)*(1-n)*N) 
        N0   = N - (N1+N2+N3+N4)

        n0_n1=0        n1_n0=0
        n1_n2=0        n2_n1=0 
        n2_n3=0        n3_n2=0 
        n3_n4=0        n4_n3=0 

        scale_dens = gamma/area
}

BREAKPOINT {
        SOLVE states
        gk = strap(N4)*scale_dens
	ik = (1e-4) * gk * (v - ek)
} 

LOCAL nexp

PROCEDURE states() {   :Computes state variables n 
        trates(v)      :             at the current v and dt.
        n = n + nexp*(ninf-n)

        alpha = strap(a*dt)
        beta  = strap(b*dt)

        n0_n1 = binom(4.0*alpha,N0)
        n1_n0 = binom(beta,N1)

        n1_n2 = binom(3.0*alpha,N1)
        n2_n1 = binom(2.0*beta,N2)

        n2_n3 = binom(2.0*alpha,N2)
        n3_n2 = binom(3.0*beta,N3)

        n3_n4 = binom(alpha,N3)
        n4_n3 = binom(4.0*beta,N4)

        N0   = N0 - n0_n1 + n1_n0
        N1   = N1 - n1_n0 - n1_n2 + n0_n1 + n2_n1
        N2   = N2 - n2_n1 - n2_n3 + n1_n2 + n3_n2 
        N3   = N3 - n3_n2 - n3_n4 + n2_n3 + n4_n3
        N4   = N4 - n4_n3 + n3_n4


        VERBATIM
        return 0;
        ENDVERBATIM
}

PROCEDURE trates(v) {  
                      
        LOCAL tinc
        TABLE ninf, nexp,a,b
	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

        nexp = 1 - exp(tinc/ntau)
}

PROCEDURE rates(vm) {  
        a = .01*vtrap(-(vm+55),10) 
        b = .125*exp(-(vm+65)/80)
        ntau = 1/(a+b) * tadj
	ninf = a*ntau/tadj
}


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 strap(x) {  :Traps for negaiv values in denominator of rate eqns.
if (x < 0) {
	strap = 0
	VERBATIM
	fprintf (stderr,"skm95.mod:strap: negative state");
	ENDVERBATIM

}else{
	strap = x
        }
}

FUNCTION binom(ppr,nnr)
{
VERBATIM
  extern double bnldev();
ENDVERBATIM
binp = ppr                      : This is a nasty trick to pass variables into the verbatim code
binn = nnr
  if (ppr<=0 || nnr <=0) {
        binom = 0
     } else {
VERBATIM
       binom_rv = bnldev(binp,(int)binn);
ENDVERBATIM
  binom = binom_rv
  }
}