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)
	ratioGibg = 0.0 : How large fraction is Gibg-bound. In the original model, this is calculated as VGCCGibg/(VGCC+VGCCGibg)	
}


NEURON {
	SUFFIX CaNtypeWillingReluctantFixedRatio
	USEION ca READ eca, cai WRITE ica
        RANGE gcanbar, ica, gcan, offmaWil, offmbWil, offha, offhb, offmaRel, offmbRel, a0mWil, a0mRel, ratioGibg
        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
}

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

BREAKPOINT {
	SOLVE states METHOD cnexp
	gcan = gcanbar*((1-ratioGibg)*mWil*mWil+ratioGibg*mRel*mRel)*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}
}