www.neuron.yale.edu

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

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?

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

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?

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

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).

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

  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.

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.

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.

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

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

Would it still do the numerical differentiation