TITLE snavb.mod

COMMENT
Sodium channel, kinetics adapted from Vanderberg & Benzanilla (1991) 
Biophys. J. 
ENDCOMMENT


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

NEURON {
	SUFFIX snavb
	USEION na READ ena WRITE ina
	RANGE m, h, gna, gbar
        RANGE binom_rv,N,Ntot,D
        RANGE a,b,c,d,f,g,i,j,y,z,vexp
        RANGE gamma
	GLOBAL q10, temp, tadj, vmin, vmax,RTF
}

PARAMETER {
	gbar = 1200   	(pS/um2)	: 0.12 mho/cm2
	temp = 5	(degC)		: original temp 
	q10  = 2.5			: temperature sensitivity
        RTF  = 24                       : RT/F at 5 degC
	v 		(mV)
	dt		(ms)
	celsius		(degC)
	vmin = -120	(mV)
	vmax = 100	(mV)

        gamma  = 35      (pS)    
        D     =  60      (1/um2)        : ? will have to check if to scale
}

ASSIGNED {
        a		
	b		
        c		
	d		
        f		
	g		
        i		
	j		
        y		
	z		
        vexp            (mV)
	ina 		(mA/cm2)
	gna		(pS/um2)
	ena		(mV)

	tadj

        N 
        Ntot
        scale_dens (pS/um2) 
        binp 
        binn   
        binom_rv
}
 
STATE { 	   
                       I4 I5  I
	      C1 C2 C3 C4 C5  O 

		       c4i4   oi
		       i4c4   io
		        i4i5 i5i
         	        i5i4 ii5
	c1c2 c2c3  c3c4 c4c5 c5o
	c2c1 c3c2  c4c3 c5c4 oc5
}

INITIAL { 
        N   = floor(D*area+0.99) 
        scale_dens = gamma/area
        a=0 b=0 c=0 d=0 f=0 g=0 i=0 j=0 y=0 z=0 vexp=0
        C1 = N
	C2 = 0 C3 = 0  C4 = 0 C5 =0   
        O = 0  
        I4 = 0  I5 = 0  I = 0
   
	c1c2=0 c2c3=0  c3c4=0 c4c5=0 c5o=0
	c2c1=0 c3c2=0  c4c3=0 c5c4=0 oc5=0
        c4i4=0 oi =0   i4c4=0 io =0
        i4i5=0 i5i=0   i5i4=0 ii5=0
}

BREAKPOINT {
        SOLVE states
        ina = (1e-4)*strap(O)*vexp*scale_dens 
} 

PROCEDURE states() {
: Note there is a better way to do it using a scheme like in the paper
: with a single binom call for each state and than divide the poplulation 
: according to the ratio of rtes  

	rates(v)         : Compute exchenge rates 
	c1c2 = binom(y,C1)
   	c2c3 = binom(y,C2) 
	c3c4 = binom(y,C3) 
	c4c5 = binom(a,C4)
	c5o  = binom(c,C5)
	c2c1 = binom(z,C2)
	c3c2 = binom(z,C3)
	c4c3 = binom(z,C4)
	c5c4 = binom(b,C5)
	oc5  = binom(d,O)
	c4i4 = binom(g,C4)
	oi   = binom(f,O)
	i4c4 = binom(j,I4)  
	io   = binom(i,I)
	i4i5 = binom(a,I4) 
	i5i  = binom(c,I5)
        i5i4 = binom(b,I5)
	ii5  = binom(d,I)

        O  = O  - oc5  - oi   + c5o  + io
        I4 = I4 - i4i5 - i4c4 + i5i4 + c4i4
        I5 = I5 - i5i4 - i5i  + i4i5 + ii5
        I  = I  - ii5  - io   + oi   + i5i 
        C2 = C2 - c2c1 - c2c3 + c3c2 + c1c2  
        C3 = C3 - c3c2 - c3c4 + c2c3 + c4c3   
        C4 = C4 - c4c3 - c4i4 - c4c5 + c3c4 + i4c4 + c5c4
        C5 = C5 - c5c4 - c5o  + c4c5 + oc5 
	C1 = N - (C2+C3+C4+C5+O+I4+I5+I)

        VERBATIM
        return 0;
        ENDVERBATIM
}

LOCAL tinc
PROCEDURE rates(v) {
        TABLE  vexp,a,b,c,d,f,g,i,j,y,z
	DEPEND dt, celsius, temp
	FROM vmin TO vmax WITH 219

        tadj = q10^((celsius - temp)/10)
        tinc = dt * tadj*1e-3	
	tabrates(v): not consistently executed from here if usetable == 1
        VERBATIM
        return 0;
        ENDVERBATIM
}


PROCEDURE tabrates(vm) {  
        y = expr(16609,1.5,0.22,RTF,vm)*tinc
        z = expr(971,-1.5,0.78,RTF,vm)*tinc
        a = expr(5750,0.42,0.99,RTF,vm)*tinc
        b = expr(4325,-0.42,0.01,RTF,vm)*tinc
        c = expr(15669,1.91,0.75,RTF,vm)*tinc
        d = expr(1361,-1.91,0.25,RTF,vm)*tinc
        f = expr(432,0.91,0.001,RTF,vm)*tinc
        i = expr(4,-0.91,0.999,RTF,vm)*tinc
        g = expr(770,0.91,0.001,RTF,vm)*tinc
        j = (g*i)/f
        vexp = vm*(exp((vm-ena)/RTF)-1)/(exp(vm/RTF)-1)
        VERBATIM
        return 0;
        ENDVERBATIM
}


FUNCTION expr(A,B,C,rtf,vm) {
     expr = A*exp(B*C*vm/rtf)
}


FUNCTION strap(x) {  :Traps for negaive 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 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
  }
}
UNITSON