modlunit Conversion Factor

NMODL and the Channel Builder.
Post Reply
AliAli
Posts: 1
Joined: Sun Oct 09, 2022 3:13 pm

modlunit Conversion Factor

Post by AliAli »

Hello,

I have been trying to implement a delayed rectifying potassium GHK current mechanism with following equations (Neiman et al. 2011):

Image

(image address: https://ibb.co/1R6xSSn )

I get the following error after checking it with modlunit:

The previous expression needs the conversion factor (6.02214e+26)
at line 45 in file drk.mod
ik = p_drk * (v * Faraday^2 / (R * T)) * ((ki - ko * exp(-Faraday * v / (R * T)<<ERROR>>))/(1 - exp(-Faraday * v / (R * T)))) * m^2

I don't quite understand where the conversion factor comes from (does it have to do with the way I have defined Faraday constant in my code?), and I have tried putting it in various places in the equation to solve the issue, but it hasn't worked. I'd appreciate any help.

Here's the code I have written for the mechanism:

Code: Select all

TITLE Delayed Rectifier Current

UNITS {
        (M)     = (mole)
        (molar) = (1/liter)
	(mM)	= (millimolar)
	(S)  	= (siemens)
	(mA) 	= (milliamp)
	(mV) 	= (millivolt)
        (L)     = (liter)
        (kel)   = (kelvin)
}

NEURON {
    SUFFIX DRK
    USEION k READ ki, ko WRITE ik
    RANGE ki, ko, ik
}

PARAMETER {
    p_drk   = 2.4e-14 (L/s)
    ki      = 112 (mM)
    ko      = 2 (mM)    
    Faraday = 96520 (coulombs/M)         :(faraday) (kilocoulombs)

    R      = 8.313424   (joule/degC)
    T      = 295.15     (kel)
    z      = 2
}

ASSIGNED { 
    v (mV)
    ik (mA/cm2)
    m_inf 
    alpha_drk (1/s)
    beta_drk (1/s)
}

STATE {
    m
}

BREAKPOINT {
    SOLVE states METHOD cnexp
    ik  = p_drk * (v * Faraday^2 / (R * T)) * ((ki - ko * exp(-Faraday * v / (R * T)))/(1 - exp(-Faraday * v / (R * T)))) * m^2
}

INITIAL {
    settables(v)
    m = m_inf
}

DERIVATIVE states {  
    settables(v)      
    m' = (m_inf - m) * (alpha_drk + beta_drk)
}

UNITSOFF

PROCEDURE settables(v) {

    m_inf     = 1.0/((1 + exp((v + 48.3)/4.19))^(1/2))
    alpha_drk = 1.0/(3.2 * exp(-v/20.9) + 3)
    beta_drk  = 1.0/(1467 * exp(v/5.96) + 9)

}

UNITSON
Thank you.
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: modlunit Conversion Factor

Post by ted »

This is a tough one--raises a lot of questions for me, even before we get to your question.

Neiman et al. 2011 doesn't ring a bell. PubMed search wasn't helpful. Maybe something in that article would explain why it was necessary to express K current in terms of the GHK current equation, and why P is in units of L/s instead of cm/s (this in particular is important).

The message modlunit gave you
previous expression needs the conversion factor (6.02214e+26)
contains a number that is 1e3 times larger than Avogadro's constant. That's unlikely to be a coincidence. Why Avogadro's constant?

Aiming for a sanity check, I inserted some UNITSOFF . . . UNITSON pairs into the mod file you sent so that modlunit would skip over things I didn't want it to complain about, then inserted these statements into the INITIAL block
printf("numerical value of Faraday is %g\n", Faraday)
printf("numerical value of Faraday/(R * T) is %g\n", Faraday/(R*T))
and got these reassuring results by calling h.finitialize():
numerical value of Faraday is 96520
numerical value of Faraday/(R * T) is 39.3364

But calling h.finitialize() also produced these error messages
exp(2556.83) out of range, returning exp(700)
exp(2556.83) out of range, returning exp(700)
exp(2556.87) out of range, returning exp(700)
exp(2556.87) out of range, returning exp(700)

Of course these are caused by the
exp(-Faraday * v / (R * T))
calls in the rhs of
ik = p_drk etc.
in the BREAKPOINT block. How do I know? Because (1) 2556.87 is the product of 65 (initial value of v) and 39.3364 (numerical value of Faraday/(R*T)), and (2) the error messages go away when I replace exp's argument with a number like 2. Which makes sense.

But here's a little bit of weirdness: if I just declare an ASSIGNED variable called bla, then stick this assignment statement in the BREAKPOINT block, modlunit says:

Code: Select all

The previous primary expression with units: 1.66054-27 
is missing a conversion factor and should read:
  (1.66054e-27)*()
 at line 53 in file xdrk.mod
bla = Faraday * v / (R * T)<<ERROR>>
Really? You may (not) recall (I certainly didn't) that 1.66054e-27 kg is numerically very close to the "atomic mass unit" defined as 1/12 of the mass of carbon 12. Great. Why does it crop up here?

So is there something strange about using the string Faraday in NMODL source code? At least for now, I think I'd use some other name for this constant.

Aside from the above, your NMODL code also has some actually understandable problems that should be fixed. But if the calculation of ik blows up when v lies in the normal range of values, the whole thing needs to be rethought. Maybe there's a grouping of variables and parameters that avoids exponential overflows.
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: modlunit Conversion Factor

Post by ted »

Rather than writing something from scratch, it is generally quicker and easier to proceed by modifying existing NMODL code that works properly, does something similar to what one wants, and follows best practices.

Among the files included with NEURON's source code (downloadable as a zip file from https://github.com/neuronsimulator/nrn) is cachan.mod, which is located in nrn/share/examples/nrniv/nmodl. This file provides an NMODL definition of a voltage-gated calcium channel that has an IV relationship described by the GHK current equation. I have examined this file and verified that its implementation of the GHK current equation is mathematically equivalent to the GHK current equation presented on p. 28 of Johnston and Wu's "Foundations of Cellular Neurophysiology" (Eq. 2.7.17).

You'll want to change a few things--here's a partial list:
the declarations in the NEURON, ASSIGNED, and STATE blocks
the names of the permeability parameter and intra- and extracellular ionic concentrations
the names of the DERIVATIVE block and the FUNCTIONs oca_ss and oca_tau
--but don't change anything in FUNCTION efun, and in FUNCTION ghk the only thing you should do is change every occurrence of 2*FARADAY to FARADAY (because the original file described a calcium channel, but you want a potassium channel and a k ion has +1 charge).

For now I'd suggest consolidating the two UNITS blocks into a single block, but don't change the declarations of FARADAY and R. Your aim at this point should be to produce something that passes modlunit and works properly.

There are some other things that I would change that have more to do with code readability than whether the code works properly (like never write a numeric constant in the range -1..1 with a leading decimal point, e.g. .1 or -.3.; instead put a 0 in front of the decimal point).

Finally, a couple of important notes about ionic concentrations. For any ion x, unless a segment has a mechanism that that WRITEs xi or xo, those concentrations will be treated as constants, with values that can be changed by hoc or python assignment statements. Furthermore, statements in the PARAMETER block that try to specify xi or xo will have no effect. When mod files are compiled, if such an assignment statement is encountered, a warning message will appear that says the assignment statement is ignored. nrnivmodl can print messages so fast that it's easy to overlook those warnings.
Post Reply