I suppose you already know that this
i = -v*(ghk/cout)*((cin-cout*exp(-v*f/r/temp))/(1-exp(-v*f/r/temp)))
deviates from the GHK equation for current (in particular the ghk/cout bit). The correct formula appears in many places, including
eqn. 5.3 on p 123 of the book "Ion Channels" by Aidley & Stanfield
and
eqn. 5.3.1 on page 123 of the book "Foundations of Cellular Neurophysiology" by Johnston & Wu.
(what a coincidence of equation and page numbers!)
Following the principle of
whenever possible, steal [great] code (Anon.)
I resorted to examining the NMODL examples that accompany NEURON; these are in
nrn/share/examples/nrniv/nmodl
of the expanded nrn*tgz source code, and MSWin users will find them in
c:\nrnxx\examples\nrniv\nmodl
The most relevant files for this case are called fh.mod and standard.inc. I pilfered the essentials from them to create my own kghk.mod:
Code: Select all
: kghk.mod
: potassium leak governed by Goldman-Hodgkin-Katz equation
NEURON {
SUFFIX kghk
USEION k READ ki, ko WRITE ik
RANGE p, i
}
UNITS {
(mV) = (millivolt)
(mA) = (milliamp)
(molar) = (/liter)
(mM) = (millimolar)
FARADAY = (faraday) (coulomb)
R = (k-mole) (joule/degC)
}
PARAMETER {
v (mV)
celsius (degC)
p = 1.2e-3 (cm/s)
}
ASSIGNED {
i (mA/cm2)
ik (mA/cm2)
ko (mM)
ki (mM)
}
INITIAL {
i = p*ghk(v, ki, ko)
ik = i
}
BREAKPOINT {
i = p*ghk(v, ki, ko)
ik = i
}
FUNCTION ghk(v(mV), ci(mM), co(mM)) (0.001 coul/cm3) {
:assume a single charge
LOCAL z, eci, eco
z = (1e-3)*FARADAY*v/(R*(celsius+273.15))
UNITSOFF
eco = co*efun(z)
eci = ci*efun(-z)
UNITSON
ghk = (0.001)*FARADAY*(eci - eco)
}
FUNCTION efun(z) {
if (fabs(z) < 1e-4) {
efun = 1 - z/2
}else{
efun = z/(exp(z) - 1)
}
}
The advantages of this particular code, aside from the fact that it expects one to specify an actual potassium permeability, are (1) it recognizes celsius as the temperature parameter, and (2) it will "play well with other mechanisms." These are important considerations if you ever want to use it in a model that has any other temperature-dependent mechanisms, or that has some other mechanism that WRITEs ki or ko. Otherwise you'll end up with code that has "magic numbers" (to specify temperature and concentrations) buried in various files, which is a code management headache that can lead to absurdities such as multiple mechanisms in the same sections having different and mutually incompatible values for temperature, equilibrium potentials, and concentrations.
"But my model is very simple and lacks any mechanism that WRITEs ki or ko. Can I specify those from hoc?"
Yes, because if a section lacks any mechanism that WRITEs an ionic concentration, that concentration is simply a range variable parameter, not a state variable, and as such it can be specified by a simple assignment statement. Here's an example of such a model:
Code: Select all
load_file("nrngui.hoc")
create soma
access soma
L = 10
diam = 10/PI // area = 100 um2
// so total current in nA = current density in mA/cm2
insert kghk
celsius = 37 // deg C
// note: the following assignment statements
// will suffice for setting intra and extracellular
// k concentrations as long as there is no mechanism
// that WRITEs ki or ko
ki = 145 // mM
ko = 3
And here is a bit of "instrumentation and control code" that will exercise this model, sweeping membrane potential from -100 mV to +100 mV in 2 mV increments and plotting i_kghk at each potential. It also sets up a control panel so you can interactively assign new values to temperature, concentrations, and permeability, and generate new plots to see what happens:
Code: Select all
// argument is membrane potential
// returned value is current at that potential
func ikghk() {
finitialize($1)
return i_kghk
}
objref g
g = new Graph(0)
g.size(-100,100,0,65)
g.view(-100, 0, 200, 65, 0, 105, 300.48, 200.32)
g.addexpr("ikghk(v)", 1, 1, 0.2, 1, 2)
proc doplot() {
v = -100
g.begin()
for (v=-100; v<=100; v=v+2) g.plot(v)
g.flush()
g.exec_menu("View = plot")
}
doplot()
{
xpanel("GHK leak",0)
xvalue("celsius","celsius", 1,"", 0, 1 )
xvalue("ki","soma.ki", 1,"", 0, 0 )
xvalue("ko","soma.ko", 1,"", 0, 0 )
xvalue("p_kghk","soma.p_kghk", 1,"", 0, 0 )
xbutton("New plot","doplot()")
xpanel(327,105)
}