GHK current equation

[GTR]

Dear Ted,
I have two problems:
First,
I want to model a high voltage activated calcium channel CaN with ca current of the form:
ica = gmax*q^2*u*ghk(v,cai,cao) where q is the activation and u the inactivation states and ghk is the GHK current equation with the form:

I = (V/Cao)*[(Cai - Caoexp(-zFV/(RT)))/(1 - exp(-zFV/(RT)))]

where Cai ,Cao are the intracellular and extracellular calcium concentrations respectively. So I made the following code based on the camchan.mod of nrn's demos.

Code: Select all

UNITS {
    (mV) = (millivolt)
    (mA) = (milliamp)
}
NEURON {
    SUFFIX CaN
    USEION ca READ eca,cai,cao WRITE ica
    RANGE gmax,ica
}
UNITS {
	:FARADAY = 96520 (coul)
	:R = 8.3134 (joule/degC)
	FARADAY = (faraday) (coulomb)
	R = (k-mole) (joule/degC)
}
PARAMETER {
    gmax  = 1.15e-3 (mho/cm2)
}
ASSIGNED {
    celsius (degC) 
    v (mV)
    eca (mV)
    ica (mA/cm2)
    cai	(mM)
    cao	(mM)
    qinf
    qtau(ms)
    uinf
    utau(ms)
}
STATE {
    q u
}
INITIAL {
	q=qinf
}
BREAKPOINT {
    SOLVE states METHOD cnexp
    ica = gmax*q^2*u*ghk(v,cai,cao)
}
DERIVATIVE states {	: exact when v held constant
	rate(v*1(/mV))
q'=(qinf-q)/qtau
u'=(uinf-u)/utau
}
FUNCTION ghk(v(mV), cai(mM), cao(mM))  (.001 coul/cm3) {
	LOCAL z
	z = 2*FARADAY*v/(R*(celsius+273.15))
	ghk=(v/cao)*((cai-cao*exp(-z))/(1-exp(-z)))
}
UNITSOFF
PROCEDURE rate(v) { LOCAL a, b :rest = -70
	TABLE qinf, qtau,uinf,utau DEPEND celsius FROM -100 TO 100 with 200
	qinf=1/(1+exp((-24.6-v)/11.3))
        qtau=1.25/cosh(-0.03(v+37.1))
        uinf=1/(1+exp((v+60)/12.5))
        utau=98+cosh(0.021(10.1-v))	
}
UNITSON

...and the interpreter finds error(illegal expression).I have no information about cai and cao values but I think the problem is in the FUNCTION. I can't understand the efun(z) function of your camchan.mod quoted here:

Code: Select all

FUNCTION ghk(v(mV), ci(mM), co(mM)) (.001 coul/cm3) {
	LOCAL z, eci, eco
	z = (1e-3)*2*FARADAY*v/(R*(celsius+273.15))
	eco = co*efun(z)
	eci = ci*efun(-z)
	:high cao charge moves inward
	:negative potential charge moves inward
	ghk = (.001)*2*FARADAY*(eci - eco)
}

FUNCTION efun(z) {
	if (fabs(z) < 1e-4) {
		efun = 1 - z/2
	}else{
		efun = z/(exp(z) - 1)
	}
}

so how can I implement the GHK FUNCTION?

and the second question:
what changes in the code if the DERIVATIVE states depend on the intracellular calcium concentration Cai and not on potential V?

rates(cai):
q'(cai)=[qinf(cai)-q(cai)]/qtau(cai)

where [cai] is modelled as:
[cai] ' =(Ical+Ican+IcaT)*c*(cai0-[cai])/tca
c :is the conversion constant
tca:the time constant of the decay
cai0:the basal intracellular calcium level
and Ical,Ican,IcaT three different calcium currents

thanks a lot,
George.

[ted]

Before wrestling with this, I took the liberty of inserting "code" markups so that it is more
readable.

You're on the right track, but the NMODL code does contain many errors.

To start with, please check it with modlunit and correct the units errors. For example,
a mechanism whose channels obey the GHK equation has a density in permeability
units, not mho/cm2

After each correction, check again with modlunit, which will catch many other kinds of
mistakes--although the error messages it produces may seem obscure or even
bizarre. For example,
qtau=1.25/cosh(-0.03(v+37.1))
generates this error message

Code: Select all

Need declaration in UNITS block of the form:
        (v)     (units)
Cannot recognize the units: v at line 78 in file x.mod
        qtau=1.25/cosh(-0.03(v+37.1)<<ERROR>>)

Of course there isn't anything wrong with v, and there definitely does NOT have to be
a declaration (v) (units) in the UNITS block. The mistake here is omission of an
asterisk between the 0.03 and the (v. Easy to forget those things. I myself make at
least 10 such mistakes every morning before breakfast.

The TABLE statement in PROCEDURE rate is also incorrect; just delete it.

I can't understand the efun(z) function of your camchan.mod

so why not replace it with something that doesn't work, right?

It was a good idea to start with camchan.mod and modify it to fit your needs. But you
have to know what things to leave alone. FUNCTION ghk and FUNCTION efun fall
into the category of "leave alone." It ain't broke, so don't fix it. Here's a hint to help you
understand efun: what happens to the numerator and denominator of z/(exp(z)-1) as
z approaches 0, and how do you deal with it? Thank you, Guillaume de l'Hopital rule (or
should we pay gratitude for this to Bernoulli instead?).