## need help understanding potassium channel code in NMODL

The basics of how to develop, test, and use models.
citron
Posts: 4
Joined: Tue May 13, 2008 2:47 pm

### need help understanding potassium channel code in NMODL

Hello, I am a bit new in NEURON and I am trying to understand code from potassium channel created by Zach Mainen. The full code is here:

Code: Select all

``````COMMENT

kv.mod

Potassium channel, Hodgkin-Huxley style kinetics
Kinetic rates based roughly on Sah et al. and Hamill et al. (1991)

Author: Zach Mainen, Salk Institute, 1995, zach@salk.edu

ENDCOMMENT

INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}

NEURON {
SUFFIX kv
USEION k READ ek WRITE ik
RANGE n, gk, gbar
RANGE ninf, ntau
GLOBAL Ra, Rb
GLOBAL q10, temp, tadj, vmin, vmax
}

UNITS {
(mA) = (milliamp)
(mV) = (millivolt)
(pS) = (picosiemens)
(um) = (micron)
}

PARAMETER {
gbar = 5   	(pS/um2)	: 0.03 mho/cm2
v 		(mV)

tha  = 25	(mV)		: v 1/2 for inf
qa   = 9	(mV)		: inf slope

Ra   = 0.02	(/ms)		: max act rate
Rb   = 0.002	(/ms)		: max deact rate

dt		(ms)
celsius		(degC)
temp = 23	(degC)		: original temp
q10  = 2.3			: temperature sensitivity

vmin = -120	(mV)
vmax = 100	(mV)
}

ASSIGNED {
a		(/ms)
b		(/ms)
ik 		(mA/cm2)
gk		(pS/um2)
ek		(mV)
ninf
ntau (ms)
}

STATE { n }

INITIAL {
trates(v)
n = ninf
}

BREAKPOINT {
SOLVE states
ik = (1e-4) * gk * (v - ek)
}

LOCAL nexp

PROCEDURE states() {   :Computes state variable n
trates(v)      :             at the current v and dt.
n = n + nexp*(ninf-n)
VERBATIM
return 0;
ENDVERBATIM
}

PROCEDURE trates(v) {  :Computes rate and other constants at current v.
:Call once from HOC to initialize inf at resting v.
LOCAL tinc
TABLE ninf, nexp
DEPEND dt, celsius, temp, Ra, Rb, tha, qa

FROM vmin TO vmax WITH 199

rates(v): not consistently executed from here if usetable_hh == 1

nexp = 1 - exp(tinc/ntau)
}

PROCEDURE rates(v) {  :Computes rate and other constants at current v.
:Call once from HOC to initialize inf at resting v.

a = Ra * (v - tha) / (1 - exp(-(v - tha)/qa))
b = -Rb * (v - tha) / (1 - exp((v - tha)/qa))
ntau = 1/(a+b)
ninf = a*ntau
}
``````

Now, what I do not understand are these lines:

Code: Select all

`` n = n + nexp*(ninf-n)``
and

Code: Select all

``nexp = 1 - exp(tinc/ntau)``
I suppose they have something to do with the equation

Code: Select all

``dn/dt = alpha(v) (1-n) - beta(V)n``

and its adjustment to current temperature, but I still do not see why they are so complex. I tried to googlesearch it for quite some time and I am getting nothing. Can anyone explain this to me please ?

ted
Posts: 5770
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:
I suppose they have something to do with the equation
Your guess is correct, they have everything to do with the analytical solution to that
equation, and an old method for taking advantage of that analytical solution in order to
speed up simulations that involve channels with HH-style kinetics.

The trick is to treat the voltage-dependent steady state values and time constants, which
govern the evolution of gating variables, as if they remain constant during a time step.
Then one can update the gating variables from the analytic solutions of the ODEs that
describe their dynamics. You'll see this a lot in "legacy code." It works, but only with fixed
time step integration.

This is an example of deprecated code. "Deprecated" means that it uses an approach
that works, but should not be employed in writing new NMODL code. An up-to-date
implementation of this particular mechanism would have a SOLVE statement that reads
SOLVE states METHOD cnexp
Instead of PROCEDURE states, there would be a DERIVATIVE block

Code: Select all

``````DERIVATIVE states {
trates(v)
n' = (ninf - n)/ntau
}``````
PROCEDURE trates would be rewritten so that it does not use dt. I don't see any
particular reason to separate the code in rates from what's in trates, but the formulas
that calculate a and b should be protected from 0/0 indeterminacy by applying l'Hospital's
rule in the vicinity of v == tha. With these changes, an acceptable implementation of
PROCEDURE trates would look like this:

Code: Select all

``````PROCEDURE trates(v) {
TABLE ninf, taun DEPEND celsius, temp, Ra, Rb, tha, qa
FROM vmin TO vmax WITH 199

a = Ra * qa * efun(-(v - tha)/qa)
b = Rb * qa * efun((v - tha)/qa)

ninf = a/(a+b)
}

FUNCTION efun(z) {
if (fabs(z) < 1e-4) {
efun = 1 - z/2
}else{
efun = z/(exp(z) - 1)
}
}``````
These changes would eliminate the need for dt and nexp. The resulting mod file would
describe a mechanism that works with either fixed time step or adaptive integration.

The INDEPENDENT statement has no effect in NEURON and can be omitted.
The comment in this line
gbar = 5 (pS/um2) : 0.03 mho/cm2
in the PARAMETER block is untrue.

citron
Posts: 4
Joined: Tue May 13, 2008 2:47 pm
Hello Ted,
thanks for the answer, hovewer, it brings up some more questions.

1. As I am new to differential equations as well, I would like to know how exactly did you get this

Code: Select all

``n' = (ninf - n)/ntau``
from the equation

Code: Select all

``dn/dt = alpha(v) (1-n) - beta(V)n``
I would like to undestand the whole derivative process of the equation.

2. in FUNCTION efun, how did you get the expression

Code: Select all

``efun = z/(exp(z) - 1)``
3. and also, why did the equations for a and b change so dramatically compared to the original ones ?

I suppose all my questions are related to the same derivative process, unfortunately, I do not see nor understand individual steps in this derivation.
Your help would be much appreciated. Thanks again.

ted
Posts: 5770
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:
citron wrote:I would like to undestand the whole derivative process of the equation.
Substitute ntau = 1/(a+b) and ninf = a/(a+b) and rearrange the RHS of the first equation.
in FUNCTION efun, how did you get the expression

Code: Select all

``efun = z/(exp(z) - 1)``
The RHS of the original equations for a and b are special instances of this general form.
why did the equations for a and b change so dramatically compared to the original ones ?
Through application of l'Hospital's rule and algebraic rearrangement.

citron
Posts: 4
Joined: Tue May 13, 2008 2:47 pm
Great, I have just derived it myself and I finally understand it.
Thanks a lot Ted.