TITLE n-calcium channel
: n-type calcium channel
: Tuomo 2023:
: -Changed GHK formalism to a simpler multiplication by membrane potential minus reversal potential. Included eca in the READ block.
: -Parametrized the half-activation and half-inactivation voltages

UNITS {
	(mA) = (milliamp)
	(mV) = (millivolt)

	FARADAY = 96520 (coul)
	R = 8.3134 (joule/degC)
	KTOMV = .0853 (mV/degC)
}

PARAMETER {
	v (mV)
	celsius 		(degC)
	gcanbar=.0003 (mho/cm2)
	ki=.001 (mM)
	cai=50.e-6 (mM)
	cao = 2  (mM)
	q10=5
	mmin = 0.2
	hmin = 3
	a0mWil =0.03 
	a0mRel =0.023 : Based on the observation that time constant is 30.67 percent larger in reluctant (PMID 25425625). a0m is inversely proportional to the time constant, and thus a0mRel is smaller than a0mWil
	zetam = 2
	vhalfm = -14
	gmm=0.1
	offmaWil = 19.88
	offmbWil = 0.0
	offha = 0.0
	offhb = 39.0
        offmaRel=29.88 : Based on the observation that voltage-dependence is shifted by +10mV by Gbg activation (PMC 2217198, PMID 25425625)
        offmbRel=10.0  : Based on the observation that voltage-dependence is shifted by +10mV by Gbg activation (PMC 2217198, PMID 25425625)
	
}


NEURON {
	SUFFIX CaNtypeWillingReluctant
	USEION ca READ eca, cai WRITE ica
        USEION VGCC READ VGCCi CHARGE 0
        USEION VGCCGibg READ VGCCGibgi CHARGE 0
        RANGE gcanbar, ica, gcan, offmaWil, offmbWil, offha, offhb, offmaRel, offmbRel, a0mWil, a0mRel
        GLOBAL hinf,minfWil,minfRel,taumWil,taumRel,tauh
}

STATE {
	mWil mRel h 
}

ASSIGNED {
	ica (mA/cm2)
	eca (mV)
        gcan  (mho/cm2) 
        minfWil
        minfRel
        hinf
        taumWil
        taumRel
        tauh
        VGCCi       (mM)
        VGCCGibgi   (mM)
}

INITIAL {
        rates(v)
        mRel = minfRel
        mWil = minfWil
        h = hinf
}

BREAKPOINT {
	SOLVE states METHOD cnexp
	gcan = gcanbar*(VGCCi*mWil*mWil+VGCCGibgi*mRel*mRel)/(VGCCi+VGCCGibgi)*h*h2(cai)
	ica = gcan*(v-eca)

}

UNITSOFF
FUNCTION h2(cai(mM)) {
	h2 = ki/(ki+cai)
}


FUNCTION alph(v(mV)) {
	alph = 1.6e-4*exp((offha-v)/48.4)
}

FUNCTION beth(v(mV)) {
	beth = 1/(exp((offhb-v)/10.)+1.)
}

FUNCTION alpmWil(v(mV)) {
	alpmWil = 0.1967*(offmaWil-v)/(exp((offmaWil-v)/10.0)-1.0)
}

FUNCTION betmWil(v(mV)) {
	betmWil = 0.046*exp((offmbWil-v)/20.73)
}

FUNCTION alpmRel(v(mV)) {
	alpmRel = 0.1967*(offmaRel-v)/(exp((offmaRel-v)/10.0)-1.0)
}

FUNCTION betmRel(v(mV)) {
	betmRel = 0.046*exp((offmbRel-v)/20.73)
}

FUNCTION alpmt(v(mV)) {
  alpmt = exp(0.0378*zetam*(v-vhalfm)) 
}

FUNCTION betmt(v(mV)) {
  betmt = exp(0.0378*zetam*gmm*(v-vhalfm)) 
}

UNITSON

DERIVATIVE states {     : exact when v held constant; integrates over dt step
        rates(v)
        mWil' = (minfWil - mWil)/taumWil
        mRel' = (minfRel - mRel)/taumRel
        h' = (hinf - h)/tauh
}

PROCEDURE rates(v (mV)) { :callable from hoc
        LOCAL a, b, qt
        qt=q10^((celsius-25)/10)
        a = alpmWil(v)
        b = 1/(a + betmWil(v))
        minfWil = a*b

	a = alpmRel(v)
        b = 1/(a + betmRel(v))
        minfRel = a*b

	taumWil = betmt(v)/(qt*a0mWil*(1+alpmt(v)))
	if (taumWil<mmin/qt) {taumWil=mmin/qt}
	taumRel = betmt(v)/(qt*a0mRel*(1+alpmt(v)))
	if (taumRel<mmin/qt) {taumRel=mmin/qt}
        a = alph(v)
        b = 1/(a + beth(v))
        hinf = a*b
:	tauh=b/qt
	tauh= 80
	if (tauh<hmin) {tauh=hmin}
}