Non-ohmic Ca current and Ca accumulation

[altaylor]

Folks,

I had a couple of questions about the implementation of a simple
calcium-buffering mechanism in NEURON/NMODL, and would be very grateful
for any help.

I'm implementing a model with a Ca channel and a KCa channel, in which
aspects of the activation/inactivation of both channels depend on the
internal Ca concentration. I've implemented this a certain way, and
wanted to make sure that this was the "standard" way (or one standard
way), and ask a few questions about implementation and performance using
this scheme.

I implemented this using three separate mod files:

ca.mod
kca.mod
caint.mod

Predictably, ca.mod implements the Ca channel, kca.mod implements the
KCa channel, and caint.mod implements the internal Ca accumulation. In
this model, Ca channel activation is voltage-dependent, and Ca channel
inactivation is a function of internal Ca concentration (i.e. h_ca is
not a state variable, it's a function of cai). The KCa channel in the
model has voltage- and calcium-dependent activation, and calcium-
dependent inactivation (m_kca and h_kca are both state variables). The
Ca accumulation model (contained in caint.mod) is very simplified -- the
equation for it is given by:

cai' = (cai_base-cai)/tau_ca - 1/(z_ca*q_e*N_A*d)*ica

where cai_base is the steady-state level of cai in the absence of Ca
current, tau_ca is the time constant of Ca buffering, d is the depth of
the shell into which incoming calcium enters (and instantaneously
mixes), z_ca is the valence of Ca (+2), q_e is the elementary charge,
and N_A is Avogadro's number.

There are other currents in my model, but only the one Ca current, and
none of the other currents are Ca-dependent.

I use the GHK current equation for calculating ica (in ca.mod).

I use the cnexp method for all my SOLVE statements in NMODL, since in
each case the time-derivatives of the state variables are linear
functions of the state variables themselves.

I'm using the reverse Euler integration method, with a largish time step
of 0.1 ms.

So, my questions are:

1) Does everything I've described above sound reasonable? Is this the
"normal" way to implement the mechanisms I've described?

2) Is it appropriate to use the cnexp method for all the solve
statements? For some reason, I have some lingering unease about
this, given that by Ca current is not ohmic. But that's silly, isn't
it?

3) How does neuron deal with non-ohmic currents (like my Ca current)
when integrating the voltage? The descriptions I've read of the
Hines scheme for integrating the membrane potential all assume ohmic
currents. Does NEURON linearize the current using the present values
of cai, cao, and v, and then proceed as normal? Can I speed things
up at all by somehow providing an analytical form for this
linearization?

4) I've tried to use the TABLE statement in all my NMODL files, to make
the simulations run faster. But the steady-state activation of my
KCa current is a function of both voltage and internal Ca, so I don't
think I can use TABLE. Is this right? I know there's a way to do a
2D lookup table for a current in HOC, but if there's a way to do that
inside the .mod file, I think I'd prefer that.

Okay, I think that's about it. My deepest thanks to anyone who is
willing to help with any of these questions!

Best,
Adam Taylor
Brandeis University

[ted]

I moved this from "Getting started" to "Adding new mechanisms and functions to
NEURON", leaving a "ghost topic" in the former.

Good questions, Adam. Everything sounds appropriate. I assume that you have
verified consistency of units with modlunit, and that you have also run tests on each
of your mechanisms to ensure that their dependence on v and cai are what you want.

You could have used the builtin constant FARADAY instead of all those other
constants in
cai' = (cai_base-cai)/tau_ca - 1/(z_ca*q_e*N_A*d)*ica
e.g.

Code: Select all

UNITS {
  FARADAY = (faraday) (coulombs)
  . . . 
}
 . . .
BREAKPOINT { SOLVE state METHOD cnexp }
DERIVATIVE state {
  cai' = (cai_base-cai)/tau_ca - (1e8)*ica/(2*d*FARADAY)
}

where (1e8) is in parentheses so that modlunit will recognize it as a scale factor.

To include instantaneous buffering, just multiply d by a factor that represents what
buffering does to "effective volume of distribution."

2) Is it appropriate to use the cnexp method for all the solve
statements? For some reason, I have some lingering unease about
this, given that by Ca current is not ohmic. But that's silly, isn't
it?

It's never silly to be concerned about something if you're not sure.
Quoting from page 225 in The NEURON Book, cnexp "is appropriate when the right
hand side of y'=f(v,y) is linear in y and involves no other states . . . "
So if the DE for the calcium current's activation variable is of the usual form
m' = (minf - m)/mtau
you're OK for that state. But is the DE for h linear in h?

3) How does neuron deal with non-ohmic currents (like my Ca current)
when integrating the voltage? The descriptions I've read of the
Hines scheme for integrating the membrane potential all assume ohmic
currents. Does NEURON linearize the current using the present values
of cai, cao, and v, and then proceed as normal?

That's my understanding of what happens.

Can I speed things
up at all by somehow providing an analytical form for this
linearization?

I don't think so--as I understand it, the linearization is accomplished by numerical
differentiation (perturbing v by a small amount).

4) I've tried to use the TABLE statement in all my NMODL files, to make
the simulations run faster. But the steady-state activation of my
KCa current is a function of both voltage and internal Ca, so I don't
think I can use TABLE. Is this right?

There is a syntax for doing a 2D lookup table, but I have only seen it used with v and
celsius. cai varies over a wide range of magnitudes, and other things may depend on
it in a very nonlinear manner. NMODL's TABLE divides the independent variable
ranges into equally sixed chunks, so I doubt that a reasonable table size (remember,
there will be 100 or 200 values of v for each value of calcium) would span a wide
enough range of cai to be very useful.

[hines]

With regard to question 3, I would encourage you to compare the results and performance of the fixed step method with decreasing dt and secondorder=0 and secondorder=2 and the results of the variable step method with decreasing atol. I am guessing you will be very pleased with the variable step method. Both in performance and if you ever plot any currents.
Note that the SOLVE method pertains only to the fixed step method and in that case the only reasonable option for DERIVATIVE blocks are cnexp and derivimplicit.

With respect to question 4, I would consider representing the kca channel as a kinetic scheme where the ca binding and voltage dependent steps are segregated. Then the voltage dependent step can use a table depending only on voltage. The performance may not be as good but the clarity is attractive. Also the entire thing can be represented with a channelbuilder including the dt dependent table for the fixed step method. Performance is always an experimental question and you may be surprised by how little some things matter.

[altaylor]

Ted, Michael,

Thanks very much for your replies. They were very helpful.

Good questions, Adam. Everything sounds appropriate. I assume that you
have verified consistency of units with modlunit, and that you have also
run tests on each of your mechanisms to ensure that their dependence on
v and cai are what you want.

Thanks! About the units and mechanism tests: Yes, I've done both.

You could have used the builtin constant FARADAY instead of all those other
constants in
cai' = (cai_base-cai)/tau_ca - 1/(z_ca*q_e*N_A*d)*ica

To include instantaneous buffering, just multiply d by a factor that
represents what buffering does to "effective volume of distribution."

Okay, thanks.

2) Is it appropriate to use the cnexp method for all the solve
statements? For some reason, I have some lingering unease about
this, given that by Ca current is not ohmic. But that's silly, isn't
it?

It's never silly to be concerned about something if you're not sure.
Quoting from page 225 in The NEURON Book, cnexp "is appropriate when the
right hand side of y'=f(v,y) is linear in y and involves no other states
. . . " So if the DE for the calcium current's activation variable is of
the usual form
m' = (minf - m)/mtau
you're OK for that state. But is the DE for h linear in h?

For the Ca channel, h is simply a function of cai, with no
time-dependence at all. I.e. h is not a state variable. The code I
used for that is:

Code: Select all

STATE  : also a var-decl block, where one declares state vars local to
       : this particular mechanism
  {
  m
  }

BREAKPOINT
  {
  SOLVE state_change METHOD cnexp
  fh(cai)  : calc h, from cai
  ica = Pbar*m^3*h*ghk(v,cai,cao,z_ca)
  }

DERIVATIVE state_change
  {
  rates(v)
  m' = (minf-m)/taum
  }

INITIAL
  {
  rates(v)
  m = minf
  }

PROCEDURE fh(cai(mM))
  {
  TABLE h FROM 0 TO 0.100 WITH 10001 
  h=1/(1+(cai/12.57e-3(mM)))
  }

PROCEDURE rates(v(mV)) 
  {  
  TABLE minf, taum FROM -150 TO 150 WITH 301
  minf=1/(1+exp(-(v+15.2(mV))/15.6(mV)))
  taum=1.8(ms)+(4.1(ms)-1.8(ms))/(1+exp(-(v+40.2(mV))/20.7(mV)))
  }

From what you've said, I gather that this is okay, since h does not need
to be integated at all.

3) How does neuron deal with non-ohmic currents (like my Ca current)
when integrating the voltage? The descriptions I've read of the
Hines scheme for integrating the membrane potential all assume ohmic
currents. Does NEURON linearize the current using the present values
of cai, cao, and v, and then proceed as normal?

That's my understanding of what happens.

Can I speed things up at all by somehow providing an analytical form for
this linearization?

I don't think so--as I understand it, the linearization is accomplished
by numerical differentiation (perturbing v by a small amount).

Okay... Sorry to pester, but how does NEURON/NMODL establish whether the
current is linear or nonlinear? Does it do some kind of syntactic
analysis on, for example, the line

Code: Select all

  ica = Pbar*m^3*h*ghk(v,cai,cao,z_ca)

from above? Would it still do the numerical differentiation if that
line looked like

Code: Select all

  ica = Pbar*m^3*h*(v-e_ca)

?

With regard to question 3, I would encourage you to compare the results
and performance of the fixed step method with decreasing dt and
secondorder=0 and secondorder=2 and the results of the variable step
method with decreasing atol. I am guessing you will be very pleased with
the variable step method. Both in performance and if you ever plot any
currents. Note that the SOLVE method pertains only to the fixed step
method and in that case the only reasonable option for DERIVATIVE blocks
are cnexp and derivimplicit.

Okay, I will experiment with different integration methods. I haven't
done that thus far.

4) I've tried to use the TABLE statement in all my NMODL files, to make
the simulations run faster. But the steady-state activation of my
KCa current is a function of both voltage and internal Ca, so I don't
think I can use TABLE. Is this right?

There is a syntax for doing a 2D lookup table, but I have only seen it
used with v and celsius. cai varies over a wide range of magnitudes, and
other things may depend on it in a very nonlinear manner. NMODL's TABLE
divides the independent variable ranges into equally sixed chunks, so I
doubt that a reasonable table size (remember, there will be 100 or 200
values of v for each value of calcium) would span a wide enough range of
cai to be very useful.

Hmmm. Okay. But just to humor me, could you point me to an example of
the syntax for a 2D lookup table? And is there an hard limit
to the size of an NMODL lookup table?

With respect to question 4, I would consider representing the kca
channel as a kinetic scheme where the ca binding and voltage dependent
steps are segregated. Then the voltage-dependent step can use a table
depending only on voltage. The performance may not be as good but the
clarity is attractive. Also the entire thing can be represented with a
channelbuilder including the dt dependent table for the fixed step
method. Performance is always an experimental question and you may be
surprised by how little some things matter.

Thank you for the suggestion. I see what you're saying about the
clarity of using a kinetic scheme for a voltage-and-calcium-dependent
process. I think I may not be able to do that for this model, but
possibly I could for a future model. A lot of work went into developing
even just the model of the KCa kinetics such that they fit our
experimental data, and developing a kinetic scheme to fit that data
would be a non-trivial undertaking, I think.

Thank you again for all your help!

Adam

[ted]

About the units and mechanism tests: Yes, I've done both.

Excellent. I wish everybody was as careful.

For the Ca channel, h is simply a function of cai, with no
time-dependence at all. I.e. h is not a state variable. The code I
used for that is

perfectly fine.

Would it still do the numerical differentiation

Yes. See the detailed discussion of fadvance() on pages 164-171 of The NEURON Book,
especially pages 167-169.

[altaylor]

Ted, Michael,

I will look at those sections of the NEURON book.

Again, thank you both for all your help!

Adam