Spatial discretization reduced the cable equation, a partial differential equation with derivatives in space and time, to a set of ordinary differential equations with first order derivatives in time. Selection of an integration method to solve these equations is guided by concerns of stability, accuracy, and efficiency (Hines and Carnevale 1995).
NEURON offers the user a choice of two stable implicit integration methods: backward Euler, and a variant of Crank-Nicholson (C-N). Backward Euler is the default because of its robust numerical stability properties. Backward Euler can be used with extremely large time steps in order to find the steady-state solution for a linear ("passive") system. It produces good qualitative results even with large time steps, and it works even if some or all of the equations are strictly algebraic relations among states.
A method which is more accurate for small time steps is available by setting the global parameter secondorder to 2. NEURON then uses a variant of the C-N method, in which numerical error is proportional to .
Both of these are implicit integration methods, in which all current balance equations must be solved simultaneously. The backward Euler algorithm does not resort to iteration to deal with nonlinearities, since its numerical error is proportional to anyway. The special feature of the C-N variant is its use of a staggered time step algorithm to avoid iteration of nonlinear equations (see 3.3.1 Efficiency). This converts the current balance part of the problem to one that requires only the solution of simultaneous linear equations.
Although the C-N method is formally stable, it is sometimes plagued by spurious large amplitude oscillations (see Fig. 7 in Hines and Carnevale (1995)). This occurs when is too large, as may occur in models that involve fast voltage clamps or that have compartments which are coupled by very small resistances. However, C-N is safe in most situations, and it can be much more efficient than backward Euler for a given accuracy.
These two methods are almost identical in terms of computational cost per time step (see 3.3.1 Efficiency). Since the current balance equations have the structure of a tree (there are no current loops), direct gaussian elimination is optimal for their solution (Hines 1984). This takes exactly thesame number of computer operations as would be required for an unbranched cable with the same number of compartments.
For any particular problem, the best way to determine which is the method of choice is to compare both methods with several values of to see which allows the largest consistent with the desired accuracy. In performing such trials, one must remember that the stability properties of a simulation depend on the entire system that is being modeled. Because of interactions between the "biological" components and any "nonbiological" elements, such as stimulators or voltage-clamps, the time constants of the entire system may be different from those of the biological components alone. A current source (perfect current clamp) does not affect stability because it does not change the time constants. Any other signal source imposes a load on the compartment to which it is attached, changing the time constants and potentially requiring use of a smaller time step to avoid numerical oscillations in the C-N method. The more closely a signal source approximates a voltage source (perfect voltage clamp), the greater this effect will be.
Nonlinear equations generally need to be solved iteratively to maintain second order correctness. However, voltage dependent membrane properties, which are typically formulated in analogy to Hodgkin-Huxley (HH) type channels, allow the cable equation to be cast in a linear form, still second order correct, that can be solved without iterations. A direct solution of the voltage equations at each time step using the linearized membrane current I(V,t) = G (V - E) is sufficient as long as the slope conductance G and the effective reversal potential E are known to second order at time . HH type channels are easy to solve at since the conductance is a function of state variables which can be computed using a separate time step that is offset by with respect to the voltage equation time step. That is, to integrate a state from to we only require a second order correct value for the voltage dependent rates at the midpoint
Figure 3.5 contrasts this approach with the common technique of replacing nonlinear coefficients by their values at the beginning of a time step. For HH equations in a single compartment, the staggered time grid approach converts four simultaneous nonlinear equations at each time step to four independent linear equations that have the same order of accuracy at each time step. Since the voltage dependent rates use the voltage at the midpoint of the integration step, integration of channel states can be done analytically in just a single addition and multiplication operation and two table lookup operations. While this efficient scheme achieves second order accuracy, the tradeoff is that the tables depend on the value of the time step and must be recomputed whenever the time step changes.
|Figure 3.5 A. The equations , are computed using the Crank-Nicholson method. Here and are determined using their values at time t.|
|Figure 3.5 B. Staggered time steps yield decoupled linear equations. is determined using x(t) , after which is determined using .|
Neuronal architecture can also be exploited to increase computational
efficiency. Since neurons generally have a branched tree structure with no
loops, the number of arithmetic operations required to solve the cable equation
by Gaussian elimination is exactly the same as for an unbranched cable with the
same number of compartments. That is, we need only O(N) arithmetic operations
for the equations that describe N compartments connected in the form of a tree,
even though standard Gaussian elimination generally takes O(N3)
operations to solve N equations in N unknowns. The tremendous efficiency
increase results from the fact that, in a tree, one can always find a leaf
compartment i that is connected to only one other compartment j,
1) the equation for compartment i (Eq. 3a) involves only the voltages in compartments i and j,
2) the voltage in leaf compartment i is involved only in the equations for compartments i and j (Eq. 3a and b).
|aii Vi + aij Vj = bi||(3a)|
|aji Vi + ajj Vj + [terms from other compartments] = bj||(3b)|
|Using Eq. 3a to eliminate the Vi term from Eq. 3b, which requires O(1) (instead of N) operations, gives Eq. 4 and leaves N-1 equations in N-1 unknowns.|
|a'jj Vj + [terms from other compartments] = b'j||(4)|
a'jj = ajj
(aij aji / aii)
b'j = bj
(bi aji / aii).
This strategy can be applied until there is only one equation in one unknown.
Assume that we know the solution to these N-1 equations, and in particular that we know Vj. Then we can find Vi from
|Efficient Gaussian elimination requires an ordering that can be found by a simple algorithm: choose the equation with the current minimum number of terms as the equation to use in the elimination step. This minimum degree ordering algorithm is commonly employed in standard sparse matrix solver packages. For example, NEURON's "Matrix" class uses the matrix library written by Stewart and Leyk (1994). This and many other sparse matrix packages are freely available on the Internet at http://www.netlib.org .|
Address questions and inquiries to
Michael Hines or Ted Carnevale
Digital preprint of "The NEURON Simulation Environment" by M.L. Hines and N.T. Carnevale,
Neural Computation, Volume 9, Number 6 (August 15, 1997), pp. 1179-1209.
Copyright © 1997 by the Massachusetts Institute of Technology, all rights reserved.