need help understanding potassium channel code in NMODL

[citron]

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)	
	tadj
}
 

STATE { n }

INITIAL { 
	trates(v)
	n = ninf
}

BREAKPOINT {
        SOLVE states
	gk = tadj*gbar*n
	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

        tadj = q10^((celsius - temp)/10)

        tinc = -dt * tadj
        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]

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

  tadj = q10^((celsius - temp)/10)
  a = Ra * qa * efun(-(v - tha)/qa)
  b = Rb * qa * efun((v - tha)/qa)

  ntau = 1/(a+b)/tadj
  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.

Other comments about this mod file:
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]

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]

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]

Great, I have just derived it myself and I finally understand it.
Thanks a lot Ted.