Using ions in multiple mechanisms

[neuromau]

Hello,

I have a few questions about using ions in mechanisms. I was looking through the mechanisms in ModelDB #124513 http://senselab.med.yale.edu/modeldb/Sh ... del=124513 (Santhakumar/Morgan dentate gyrus model). I also emailed the programmers my question, but wanted to get other general viewpoints here regarding defining ions.

1. Almost every mechanism in this model defines its own ion. For example, the K+ channels, rather than each using the 'k' ion, invent their own K+ ions. I think that has the effect of NEURON thinking they are all separate ion types, and not counting them towards the potassium concentration. I can't think of why one should do this, unless the current consisted of multiple ion types and the proportion of each type varied or was unknown. And in that case, shouldn't one just use the NONSPECIFIC current?
Ex: gskch.mod calcium-activated potassium channel (non-voltage-dependent) uses ion: sk
CaGk.mod calcium activated K channel (voltage dependent) uses ion: k
bgka.mod Borg-Graham type generic K-A channel uses ion: k
ichan2.mod (which has K+ fast and slow delayed rectifiers, in addition to other channels, all in one mechanism) uses ions: kf and ks
--> So, for K+, we have sk, k, kf, and ks ions

2. One of the mechanisms defines separate ions for the fast and slow components of its current. Again, I don't know what is to be gained from this. If it's easier to calculate the fast and slow components separately, shouldn't one still add them together at the end and set a single ion's current equal to the sum of the components?
Ex: hyperde3.mod (HCN channel) uses the following ions: hyf, hys, hyhtf, hyhts
(the other thing that I find odd about this mechanism is both the control (hy?) and experimental (hyht?) current types seem as if they would both be active at the same time.)

Thanks for any thoughts/feedback you may have.

~ Marianne

[ted]

Good questions.

1. Almost every mechanism in this model defines its own ion. For example, the K+ channels, rather than each using the 'k' ion, invent their own K+ ions. I think that has the effect of NEURON thinking they are all separate ion types, and not counting them towards the potassium concentration. I can't think of why one should do this, unless the current consisted of multiple ion types and the proportion of each type varied or was unknown. And in that case, shouldn't one just use the NONSPECIFIC current?
Ex: gskch.mod calcium-activated potassium channel (non-voltage-dependent) uses ion: sk
CaGk.mod calcium activated K channel (voltage dependent) uses ion: k
bgka.mod Borg-Graham type generic K-A channel uses ion: k
ichan2.mod (which has K+ fast and slow delayed rectifiers, in addition to other channels, all in one mechanism) uses ions: kf and ks

In general it can be difficult to know a programmer's intent, and in some cases it may even be more difficult to decide what is intentional and what is accidental (especially when the topic is NMODL source code--not least because NMODL's strange syntax tends to invite what might be called "superstitious programming" i.e. odd bits and pieces of code or syntax that serve no evident purpose, but got there by virtue of rote replication of excerpts lifted from prior examples of "working code").

Here's one situation in which it may be useful to define different ionic species: suppose there is a close anatomical association between a particular class of Ca channels and a particular class of Ca-gated K channels, so that when the Ca channels open, there is a very localized accumulation of Ca that affects the K channels. Also suppose that current through other kinds of Ca channels has little effect on the Ca-gated K channels. In this case, it might make sense to define a new ion for the Ca channels that affect the Ca-gated K channels, so it can have its own accumulation mechanism. Of course this opens questions about how Ca is cleared from this special pool (including membrane transport and ion exchanges between this pool of intracelluar Ca and all the other intracellular Ca), but those are different issues. Did the example you cite do something like that?

If not, maybe the programmer's goal was to facilitate identification of which potassium current component was generated by which mechanism. There are better ways to accomplish such a goal--ways that do not interfere with proper management of mass balance. Example:

Code: Select all

NEURON {
  SUFFIX sk
  USEION k READ ek WRITE ik
  RANGE g, i : so you can define ASSIGNED i that will be known to hoc as i_sk
 . . .
ASSIGNED {
  i (mA/cm2)
  ik (mA/cm2)
  g (S/cm2)
 . . .
BREAKPOINT {
 . . .
  g = whatever  
  i = g*(e-ek)
  ik = i : because this mechanism has to WRITE ik
}

Note the following:
1. The mechanism WRITEs ik, so it affects charge balance AND k mass balance.
2. The mechanism uses an auxiliary variable i, declared to be RANGE, so that one can discover the current generated by this mechanism (value will be known to hoc as i_sk). This variable does NOT affect charge balance or mass balance.
3. SUFFIX and variable names have been deliberately crafted to be concise yet mnemonic, while avoiding ugly hiccups--i_sk is far preferable to isk_gskch. Much existing NMODL code seems to be afflicted by very poor naming strategies that reduce code legibility.

Code: Select all

2. One of the mechanisms defines separate ions for the fast and slow components of its current. Again, I don't know what is to be gained from this. If it's easier to calculate the fast and slow components separately, shouldn't one still add them together at the end and set a single ion's current equal to the sum of the components?
Ex: hyperde3.mod (HCN channel) uses the following ions: hyf, hys, hyhtf, hyhts
(the other thing that I find odd about this mechanism is both the control (hy?) and experimental (hyht?) current types seem as if they would both be active at the same time.)

That is a bit of a wart, isn't it? Not to mention its unnecessary INDEPENDENT statement and superfluous VERBATIM block. And it uses the long-deprecated "analytical solution" approach to coding HH-style mechanisms, and contains the (unnecessary for this mechanism, and unused) FUNCTION vtrap. I can see that one might want to know the values of the fast and slow components, but this is doable by calculating auxiliary variables that can then be summed to find the total ik.

[neuromau]

In this case, it might make sense to define a new ion for the Ca channels that affect the Ca-gated K channels, so it can have its own accumulation mechanism.

Yes, they also have many different Ca2+ ions, and as there are some Ca2+ dependent K+ channels, your note likely explains the multiple Ca2+ ions. I haven't yet found an explanation of why the various K+ ions are needed though, unless it was, as you say, just re-using old code that worked or wanting to track currents but not taking advantage of the suffixes. Your response was heartening enough that I will now go through each K+-channel mechanism and standardize the ions used!We'll see if that affects the simulation or if the programmer turns out to have a reason for the unique ion types.

I can see that one might want to know the values of the fast and slow components, but this is doable by calculating auxiliary variables that can then be summed to find the total ik.

Just as I was hoping - I will rewrite this one as well.

Not to mention its unnecessary INDEPENDENT statement and superfluous VERBATIM block. And it uses the long-deprecated "analytical solution" approach to coding HH-style mechanisms, and contains the (unnecessary for this mechanism, and unused) FUNCTION vtrap.

I have been cleaning up all of these mechanisms, documenting and renaming them, in hopes of having straightforward mod files in our next model release. My cleaning had been limited to removing STATE_DISCONTINUITYs and those unnecessary VERBATIM blocks, but now that I know these are issues, I'll fix them as well.

Thanks for your help!

[neuromau]

For the sake of completeness: I also realized how it was possible that the mechanism contained currents for both the control and the treated model, without applying both currents to each model:

Ex: hyperde3.mod (HCN channel) uses the following ions: hyf, hys, hyhtf, hyhts
(the other thing that I find odd about this mechanism is both the control (hy?) and experimental (hyht?) current types seem as if they would both be active at the same time.)

Wherever the mechanism is inserted, either the control current conductance or the treated current conductance was set to 0, rendering that type of current 0 as well.

For example, in a treated animal, after inserting the mechanism and setting the treated conductance appropriately, you would cancel out the control one with the following:

Code: Select all

ghyhtsbar = 0

so that the control conductance in the mod file will become 0:

Code: Select all

ghyhts = ghyhtsbar * hyhts* hyhts

and then the calculated control current is also 0.

~Marianne

[ted]

they also have many different Ca2+ ions, and as there are some Ca2+ dependent K+ channels, your note likely explains the multiple Ca2+ ions.

but only if there is an accumulation mechanism for each ion, i.e. something with a
USEION x READ xi (and/or xo), ix WRITE xi (and/or xo
statement in its NEURON block, where x is the name of some ionic species.

I will now go through each K+-channel mechanism and standardize the ions used!We'll see if that affects the simulation or if the programmer turns out to have a reason for the unique ion types.

You might want to check the new mechanism against the old to make sure that both generate numerically identical results when subjected to a voltage clamp scenario (in a single compartment model with SEClamp attached so no other mechanisms are necessary). For Ca-dependent mechanisms it may be necessary to run such tests for two or more values of Ca concentration.

I have been cleaning up all of these mechanisms, documenting and renaming them, in hopes of having straightforward mod files in our next model release. My cleaning had been limited to removing STATE_DISCONTINUITYs and those unnecessary VERBATIM blocks, but now that I know these are issues, I'll fix them as well.

The trick is to do it in a way that doesn't break anything that already works. After all, the virtues of the existing code, warts included, are (1) it does work, and (2) it was useful for addressing a research problem. Will the revised code find its way into ModelDB? Future users would find that to be a great benefit.

[neuromau]

Will the revised code find its way into ModelDB? Future users would find that to be a great benefit.

Absolutely! My time investment, especially in the time taken for documentation, is really only worth it if I share the product with others. I've gotten through most of the code, rewriting, comparing old and new results, and documenting as I go. I put off the mechanisms for last, but I think the end is finally in sight.

You might want to check the new mechanism against the old to make sure that both generate numerically identical results when subjected to a voltage clamp scenario (in a single compartment model with SEClamp attached so no other mechanisms are necessary).
...
And it uses the long-deprecated "analytical solution" approach to coding HH-style mechanisms

When I updated the hyperde3 mechanism to use the DERIVATIVE block (as in http://www.neuron.yale.edu/phpBB/viewto ... 46&start=0), the total current is slightly different than the previous version (when I insert the mechanism into a single compartment with an SEClamp and measure the total current through the mechanism at soma(0.5)). Of the 16,000 time steps in my simulation, 14,000+ of them have a slightly different current. The differences in the current range from -24 to 7 pA, with the biggest differences occurring as the voltage clamp switches to a new level. For context, the currents at these 'timepoints of greatest difference' are -780 nA and -280 nA. Is this expected and acceptable? Or have I done something wrong?

Here is the part of the mechanism that I updated (sorry, the two components are still separated .... I will fix them next):

Code: Select all

BREAKPOINT {
	SOLVE states METHOD cnexp

	ghyf = ghyfbar * hyf*hyf
	ihyf = ghyf * (v-ehyf)
	
	ghys = ghysbar * hys*hys
	ihys = ghys * (v-ehys)
	
	myi = ihyf + ihys : this just sums the two currents for easy readout, till I merge them together
		}
 
UNITSOFF
 
INITIAL {
	trates(v)
	
	hyf = hinf
    hys = hinf
}

DERIVATIVE states {	:computes hyf and hys at current v and dt 
	trates(v)

	hyf' = (hinf-hyf)/hyftau
	hys' = (hinf-hys)/hystau
}
 
LOCAL q10
 
PROCEDURE trates(v) {  :Computes rate and other constants at current v.
                      :Call once from HOC to initialize inf at resting v.
	TABLE hinf, hyftau, hystau	
	DEPEND celsius 
	FROM -120 TO 100 WITH 220
                           
    q10 = 3^((celsius - 6.3)/10)
       
	hinf =  1 / (1 + exp( (v+91)/10 ))

	:"hyf" FAST CONTROL Hype activation system
	hyftau = (14.9 + 14.1 / (1+exp(-(v+95.2)/0.5)))/q10

	:"hys" SLOW CONTROL Hype activation system
	hystau = (80 + 172.7 / (1+exp(-(v+59.3)/-0.83)))/q10

}

UNITSON

[ted]

Comment out the TABLE stuff (all three lines) and see what happens to the discrepancy.

[neuromau]

Without the TABLE lines, the discrepancy increases.

I guess that I can't make analogous changes to what was done in the other post (http://www.neuron.yale.edu/phpBB/viewto ... 46&start=0), probably because the inf and tau expressions in the hyperde3 mechanism are not exactly of the form a/(a+b) and 1/(a+b). So I thought of deriving the equation in the states DERIVATIVE from the equation in the states PROCEDURE; perhaps then I could figure out what I have to do differently? I tried this for the mechanism in the other post and am a bit stuck, working backwards from the target equation and forwards from the starting equation (my math skills are rusty...).

I couldn't get the images I have to post in the preview box, so I emailed them to you separately.

[neuromau]

I tried this for the mechanism in the other post and am a bit stuck, working backwards from the target equation and forwards from the starting equation (my math skills are rusty...)

I think that the -dt*tadj/ntau term would be small enough that I can use the approximation "as h --> 0, exp(h) --> h + 1". In that case, it is easy to double check the translation. Still not sure why there is a discrepancy though.

[ted]

The discrepancy in results is a consequence of using TABLEs in mod files to set up lookup tables that implement piecewise linear approximations for the nonlinear voltage-dependence of time constants and steady-state values of the gating variables. This approximation, which is typically done to speed up simulations, comes at the cost of a slight reduction of accuracy (typically less than 0.1%, which is much smaller than the errors that afflict the experimental data that serve as the empirical basis for computational models). The error it introduces is proportional to the coarseness of the piecewise linear approximation. For example,

Code: Select all

PROCEDURE rates(v) {
  TABLE hinf DEPEND celsius FROM -120 TO 100 WITH 220
  hinf = f(v, celsius)
}

sets up a table with 220 entries for v and the corresponding hinf with values that range from -120 to 100 mV (almost, but not quite, one entry per millivolt--to get 1 entry per mV, the table would need 100-120+1 = 221 entries, i.e. FROM -120 TO 100 WITH 221). To reduce the error of the piecewise linear approximation by a factor of 2, just double the tablesize parameter, i.e. the TABLE statement becomes

Code: Select all

TABLE hinf DEPEND celsius FROM -120 TO 100 WITH 440
: "WITH 441" to get exactly two entries per mV

Returning to the case at hand, when an NMODL-specified mechanism has one or more TABLE statements, it automatically has a variable called usetable that can be employed to activate or deactivate use of the table from hoc. Comparing the currents generated by the new and old versions of the mechanism reveals that the difference is 0 when both use or don't use the TABLE. However, setting either implementation's usetable to 0 while leaving the other one == 1 brings out a difference that emerges slowly over time even when the cell is simply held at -65 mV. Disabling table usage by both implementations makes the difference between their currents vanish.

Now back to more general considerations--

Employing a piecewise linear approximation to the nonlinear voltage-dependence of the time constants and steady-state values of the gating variables introduces a small error. There is a tradeoff with every approximation; here the tradeoff is speed vs. accuracy. The question is always whether it is worth the cost, and the ultimate way to judge is to see if it substantively affects simulation results. My guess for this mechanism is that the error won't have much effect on most simulations. Furthermore, for any mechanism in any empirically-based model, the use of TABLEs won't be the only source of errors). But you can always deactivate use of TABLEs via at most a few hoc statements, so it's easy to check for yourself in any particular situation.

[phendric]

Will the revised code find its way into ModelDB? Future users would find that to be a great benefit.

Absolutely! My time investment, especially in the time taken for documentation, is really only worth it if I share the product with others. I've gotten through most of the code, rewriting, comparing old and new results, and documenting as I go. I put off the mechanisms for last, but I think the end is finally in sight.

Has this since been published to ModelDB? I've not found it if it has been.

[neuromau]

It has not been published yet. I'll post a note in this thread when it gets uploaded.

~ Marianne