: Reference: Colbert and Pan 2002

NEURON	{
	SUFFIX NaTa
	USEION na READ ena WRITE ina
	RANGE gbar, g, ina
}

UNITS	{
	(S) = (siemens)
	(mV) = (millivolt)
	(mA) = (milliamp)
}

PARAMETER	{
	gbar = 0.00001 (S/cm2)

	malphaF = 0.182
	mbetaF = 0.124
	mvhalf = -48 (mV)
	mk = 6 (mV)

	halphaF = 0.015
	hbetaF = 0.015
	hvhalf = -69 (mV)
	hk = 6 (mV)
}

ASSIGNED	{
	v	(mV)
	ena	(mV)
	ina	(mA/cm2)
	g	(S/cm2)
	celsius (degC)
	mInf
	mTau
	mAlpha
	mBeta
	hInf
	hTau
	hAlpha
	hBeta
}

STATE	{
	m
	h
}

BREAKPOINT	{
	SOLVE states METHOD cnexp
	g = gbar*m*m*m*h
	ina = g*(v-ena)
}

DERIVATIVE states	{
	rates()
	m' = (mInf-m)/mTau
	h' = (hInf-h)/hTau
}

INITIAL{
	rates()
	m = mInf
	h = hInf
}

PROCEDURE rates(){
  LOCAL qt
  qt = 2.3^((celsius-23)/10)

	UNITSOFF
		mAlpha = malphaF * vtrap(-(v - mvhalf), mk)
		mBeta = mbetaF * vtrap((v - mvhalf), mk)

		mInf = mAlpha/(mAlpha + mBeta)
		mTau = (1/(mAlpha + mBeta))/qt

		hAlpha = halphaF * vtrap(v - hvhalf, hk) : ng - adjusted this to match actual Colbert & Pan values for soma model
		hBeta = hbetaF * vtrap(-(v - hvhalf), hk) : ng - adjusted this to match actual Colbert & Pan values for soma model

		hInf = hAlpha/(hAlpha + hBeta)
		hTau = (1/(hAlpha + hBeta))/qt
	UNITSON
}

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