Converting mechanisms to differential equations

[sl]

Many mechanisms in ModelDB and elsewhere do not utilize the DERIVATIVE block for one reason or another, instead opting to use the 'dt' variable directly. I got the feeling that this practice is not optimal, so I tried to convert any and all mechanisms like that to use the derivative block. However, I have come upon forms that do not have apparent equivalent differential forms. For example, an rough excerpt from one mechanism I am using:

Code: Select all

PROCEDURE states() {	: exact when v held constant
	evaluate_fct(v)
	n = n + n_exp * (n_inf - n)
}

UNITSOFF
PROCEDURE evaluate_fct(v(mV)) { LOCAL a,b,v2

	v2 = v - vtraub

	a = 0.032 * (15-v2) / ( exp((15-v2)/5) - 1)
	b = 0.5 * exp((10-v2)/40)
	tau_n = 1 / (a + b) / tadj
	n_inf = a / (a + b)

	n_exp = 1 - exp(-dt/tau_n)
}
UNITSON

By looking at the specification, it is apparent that n', which is what we want to get to put into the derivative block, has this form:

Code: Select all

n' = (1 - exp(-dt/tau_n)) *  (n_inf - n) / dt

This trick of using the n_exp variable is seen in many mechanisms, and I see no way of solving for dt in the above form.

So, my question(s):

Is there a way to convert such a mechanism to a differential equation form?

Where does this approach come from (perhaps a paper/book)? This is just so I can look at the original problem specification and see if I can re derive the equations myself.

Thanks.

[sl]

Yikes, my apologies. I forgot the definition of the differential.

All one has to do is take the limit as dt goes to 0. For the above equation this will get you:

Code: Select all

n' = (n_inf - n) / tau_n

Which is the standard exponential relationship. I have verified that both implementations result in identical outputs in NEURON.

[ted]

Thanks for catching and following up on your own error.


What you are writing about is an old trick for reducing the computational burden required
to do fixed time step integration of the 1st order differential equations (DEs) for the gating
variables of HH-style voltage-gated currents. The strategy was to assume that the time
constants and steady state values for gating variables remained constant during a time
step. This allowed the gating variables to be advanced to the next time step from the
analytical solution to the DEs, which is faster than numerical integration.

There is a lot of old NMODL code floating around that uses this trick, which works fine
for fixed time step integration. But it doesn't work for adaptive integration (variable order,
variable time step integration), which can be useful for dealing with complex models that
run too slowly with fixed step integration, or models that become unstable unless
unreasonably small values of dt are specified, or when you have model in which precise
timing of threshold crossings or other events is important (because you don't want
events to be forced onto dt boundaries).

When NMODL encounters a mod file that contains a DERIVATIVE block, it generates an
output file that contains multiple compiled versions of the same mechanism--one for each
of the various numerical integrators that NEURON can use. And one of these is indeed
the old "fixed time step trick for quickly simluating HH mechanisms." NEURON of course
is smart enough to know which one to use with which integrator. All of this goes on behind
the scenes and requires no user intervention--it just works.

[ted]

A separate item that I forgot to mention:

It often happens that the formula for a rate constant is of a form that may reduce to 0/0.
This is the case with the formula for your "a" rate constant, which is of the form
z / (exp(z) - 1)
and becomes 0/0 when v == 15.

When 0/0 happens under MSWin your simulation may stop because of a floating point
exception; the floating point libraries used in UNIX/Linux/OS X deal with this much more
gracefully. However, even if you use NEURON under UNIX/Linux/OS X, never forget that
most NEURON users are also MSWin users. And it is almost always very easy to fix
rate constant formulas to avoid 0/0, because the troublesome ones can almost always
be rewritten in the form z/(exp(z) - 1), which is amenable to Lhospital's rule as
expressed by this function

Code: Select all

FUNCTION efun(z) {
  if (fabs(z) < 1e-4) {
    efun = 1 - z/2
  } else {
    efun = z/(exp(z) - 1)
  }
}

Then you can fix this particular example by replacing

Code: Select all

   a = 0.032 * (15-v2) / ( exp((15-v2)/5) - 1)

with

Code: Select all

   a = 0.032 * 5 * efun((15-v2)/5)