TITLE Markov sodium channel of NaV 1.6 

COMMENT

from Balbi et al 2017, modified by Carol Upchurch and Christopher Knowlton

ENDCOMMENT

NEURON {
       SUFFIX NaMark
       USEION na WRITE ina
       RANGE gnabar, gna,  ina, o1, shift, scale,hshift,dist, oinit,DISABLE,c1c2,c2c1,o1c2,c2o1,hinit, BALBI, aiscorr, FAST, allcorr
       :GLOBAL zhalf
}

UNITS {
	(mV) = (millivolt)
	(mA) = (milliamp)
	(nA) = (nanoamp)
	(pA) = (picoamp)
	(S)  = (siemens)
	(uS) = (microsiemens)
	(nS) = (nanosiemens)
	(pS) = (picosiemens)
	(um) = (micron)
	(molar) = (1/liter)
	(mM) = (millimolar)		
}

CONSTANT {
}

PARAMETER {
	v (mV)
	celsius (degC)
	ena = 50 (mV)
	gnabar = 550.0e-6
	shift = 0
	dist = 100
	scale = 1
	hshift = 0
	oinit=-1
	DISABLE = 1
	BALBI = 0
	aiscorr=0
	allcorr=0
	FAST=1e-6
}

ASSIGNED {
	ina (mA/cm2)
    gna (pS/um2) 

    i1c1
    c1i1

    o1i1
    i1o1 : 0
 
    o1c1
    c1o1
    
    i2o1
    o1i2
    
    i2c1
    c1i2
    
    i1i2
    i2i1
    
    c2c1
    c1c2
}

STATE {
	c1   FROM 0 TO 1
	i1   FROM 0 TO 1
	i2   FROM 0 TO 1
	o1   FROM 0 TO 1
	c2   FROM 0 TO 1
}

INITIAL {
	rates(v)
	if(oinit < 0) { 
		SOLVE kin
		STEADYSTATE sparse
	} else {
	:c1 = 0.5-2e-6
	c1 = 1-oinit
	i1 = 0
	i2 = 0
	o1 = oinit
	c2=0
	}
}

BREAKPOINT {
	SOLVE kin METHOD sparse
    gna = gnabar*o1  
	ina = gna * (v - ena)
}

KINETIC kin {
	rates(v)
	~ o1 <-> c2 (o1c1,c1o1) : opening
	:~ o1 <-> c1 (o1c1,c1o1) : opening
	
	~ c1 <-> c2 (c1c2,c2c1) : 10 1
	~ o1 <-> i1 (o1i1,i1o1) : inactivating (fast)
	
	~ c2 <-> i2 (c1i2,0) : closed state inactivation
	:~ c2 <-> i1 (0.03*c1i1,0) : closed state inactivation
	
	~ o1 <-> i2 (o1i2,i2o1) : inactivating (slow)
	~ i1 <-> c1 (i1c1,c1i1)	: de-inactivating
	
	~ c1 <-> i2 (c1i2,i2c1)
	
	:~ i1 <-> i2 (i1i2,i2i1)

	CONSERVE  o1 + i1 + i2 +c1 +c2 = 1
}


PROCEDURE rates( v (mV) ) {
	LOCAL vm, inact, slowfrac, faster
	vm = v+shift
	
	slowfrac=dist
	faster= FAST

	c1o1 = 5*scale*DoubSig(vm,0,0,0,14.0, -2+allcorr,-3.8-aiscorr) : 3.5 for single, 3.0 for realistic morphology or reverse that idk
	o1c1 = 5*scale*DoubSig(vm,3.0,-48.0,9.0,0,0,0) : 18


	inact = DISABLE*0.5*scale*DoubSig(vm,0.2,-42.0+hshift,6.0,20.0, 10.0+hshift,-12.0) : 10 for better morphological model
	
	o1i1 = inact*(1-slowfrac)
	i1o1 = (1-slowfrac)*1e-8
	
	c1c2 = DoubSig(vm,0.3,-20,10.0,0,0,0)
	c2c1 = DoubSig(vm,0.3,-20,-10.0,0,0,0)

	o1i2 = slowfrac*inact
	i2o1 = slowfrac*1e-8
	
	
	i1c1 = scale*DoubSig(vm,faster,-65.0,10.0,0,0,0)
	c1i1 = DISABLE*scale*DoubSig(vm,0,0,0,faster,-65.0,-11.0)
	
	i2c1 = 20.0*DoubSig(vm,1.8e-3,-60.0,5.0,0,0,0)
	:c1i2 = slowfrac*DoubSig(vm,0.1,-25.0,-5.0,0,0,0) : actually c2i2
	c1i2 = 20.0*dist*DoubSig(vm,0,0,0,2.0e-2,-25.0,-5.0) :1000x : or here :above -50, it accumulates

	
	:i1i2 = 0:i1c1*c1i2
	:i2i1 = 0:i2c1*c1i1
	
	i1i2 = 20.0*dist*DoubSig(vm,0,0,0,2.0e-2,-25.0,-5.0) :1000x : or here :above -50, it accumulates
	i2i1 = 2.0*DoubSig(vm,1.8e-3,-50.0,5.0,0,0,0) : 100x : add parameter here 50/5
	
	:c1i2 = o1i2
	:o1i2=0
	
}

FUNCTION boltz(x,y,z) {
               LOCAL arg
                arg= -(x-y)/z
                if (arg > 50) {boltz = 0}
				else {if (arg < -50) {boltz = 1}
                else {boltz = 1.0/(1.0 + exp(arg))}}
}

FUNCTION DoubSig(v,bh,vh,kh,bd,vd,kd){
	LOCAL boltz1, boltz2
	if(bh > 0){
		boltz1 = bh*boltz(v,vh,-kh)
	} else {
	    boltz1 = 0
	}
	if(bd > 0){
		boltz2 = bd*boltz(v,vd,-kd)
	} else {
	    boltz2 = 0
	}
	DoubSig = boltz1 + boltz2
	
}