TITLE  virtual Na current
 
UNITS {
	(nS) = (nanosiemens)
	(nA) = (nanoamp)
	(mV) = (millivolt)
}
 
NEURON {
        POINT_PROCESS VirNa
        RANGE  gnavir,ena,del,dur,ivirna,m,h
        RANGE  minfshift,mslopefact,hinfshift,hslopefact
        NONSPECIFIC_CURRENT ivirna
}
 
: INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}
 
PARAMETER {
        gnavir  = 0.0 (nS)
        ena  = 60.0  (mV)
        del = 1000 (ms)
        dur = 1000 (ms)
        hf1 = 0.77 (1)
        minfshift = 0 (mV)
        mslopefact = 1 (mV)
        hinfshift = 0 (mV)
        hslopefact = 1 (mV)
}
 
STATE {
        m h hs
}
 
ASSIGNED {
	v (mV)
	ivirna (nA)
}
 
BREAKPOINT {
	SOLVE states METHOD cnexp
       
	at_time(del)
	at_time(del+dur)
	if (t < del + dur && t > del) {
		ivirna = (0.001)*gnavir*m*m*m*(h+hs)*(v - 60)
	} else {
		ivirna = 0
	}
}
 
UNITSOFF
 
INITIAL {
	:m = boltz(v,-30.0907+minfshift,9.7264*mslopefact)
	m = boltz(v-minfshift,-30.0907,9.7264*mslopefact)
	h = hf1*boltz(v-hinfshift,-54.0289,-10.7665*hslopefact)
	hs = (1-hf1)*boltz(v-hinfshift,-54.0289,-10.7665*hslopefact)
}

DERIVATIVE states {  :Computes state variables m, h, and n 
	LOCAL minf,hinf,hsinf,ninf,pinf,qinf,mtau,htau,hstau,ntau,ptau,qtau
	:minf = boltz(v,-30.0907+minfshift,9.7264*mslopefact)
	minf = boltz(v-minfshift,-30.0907,9.7264*mslopefact)
	hinf = hf1*boltz(v-hinfshift,-54.0289,-10.7665*hslopefact)
	hsinf = (1-hf1)*boltz(v-hinfshift,-54.0289,-10.7665*hslopefact)
	:mtau = 1/(linoid(v,-19.565,-0.50542,0.79992)+exponential(v,3.0212,-7.4630e-3)) + 0.01
	mtau = 1/(linoid(v-minfshift,-19.565,-0.50542,0.79992)+exponential(v-minfshift,3.0212,-7.4630e-3)) + 0.01
	htau = 1/(exponential(v-hinfshift,5.0754e-4,-6.3213e-2)+exponential(v-hinfshift,9.7529,0.13442)) + 0.4
	hstau = 1/(exponential(v-hinfshift,1.3134e-02,8.5636e-04)+exponential(v-hinfshift,13.165,0.14556)) + 0.4
	m' = (minf-m)/mtau
	h' = (hinf-h)/htau
	hs' = (hsinf-hs)/hstau
}
 
FUNCTION boltz(x,y,z) {
	boltz = 1/(1 + exp(-(x - y)/z))
}

FUNCTION linoid(v,a,b,c) {
	linoid = c*(a+b*v)/(exp(a+b*v)-1)
}

FUNCTION exponential(v,d,e) {
	exponential = d*exp(e*v)
}

UNITSON