The fact that there are MANY SIGNIFICANT FIGURES in the value of the output DOES NOT necessarily indicate that the value is ACCURATE! There are a variety of possible numerical errors arising from using large compartment sizes, long time steps in integration, or poor numerical integration method.
Because computers were relatively slow in the 1970s, and we wanted to begin simulations of propagating impulses find the fastest way to integrate the HH equations with minimal possible errors, Fidel Ramon and I (Moore & Ramon, 1974) compared the shape of the action potential for several methods of numerical integration over a range of time increments used in the integration. We took the action potential shape to which the simulations converged as we reduced the integration time step to be the "true" or accurate solution. Hoping to find optimal choices, we measured the speed of computation, looking for some balance between a "good " solution taking a long computation time on one hand and a short computer run yielding a "fair" or poorer solution on the other.
We found the
shape of the computed action
potential to be method-invariant but the latency was very dependent on the integration
method,
some reduced the latency and others increased it. Of course, in either case, reductions
in the
integration time step moved the impulse closer to the convergent solution.
Numerical Instabilities. For all of the integration methods tried, the time
step could be
increased to 50uS with little change in impulse shape but for time steps equal to or
greater
than 60uS
numerical
instabilities and
oscillations showed up in several methods and serious deviations appeared
even with the best method (Runga-Kutta) at 80uS. Following a smooth upsweep to the
peak of
the spike, the oscillations were triggered by an inaccurate prediction of the sodium
current
where it drops suddenly, following its first peak value. This frequently overlooked "notch" in the sodium
current
results from rapid changes in its driving voltage, as the membrane potential first
approaches
and then recedes from the sodium reversal potential.This is illustrated in the
figure from our work shown here.
Threshold. In the process. we also searched for the threshold voltage
for a patch of HH
membrane, expecting to confirm the FitzHugh & Antosiewicz
(1959) value of -6.4010523 mV. We
plotted the apparent threshold voltage as a function of the integration step size and
chose, as
"the best value", the intercept on the voltage axis. We were astonished to find a value
of -
6.5077040 mV which differed by 1.5% from their stated value, having a resolution of
one part in
108! However, it was reassuring to see that, while different numerical integration
methods had
different slopes - both positive and negative, they had a common intercept or value for
an
infinitesimal time step. In their simulations FitzHugh & Antosiewicz employed only a
single time
step and single numerical integration method (probably because their computer, an IBM
704,
was still relatively slow and many runs were required to converge on threshold). They
used a
numerical integration method which seemed intuitively very sensible because it
checked on the
size of the membrane voltage increment called for in the next time. If necessary, it
reduced the
size of the time step so that only increments voltage smaller than a predetermined
value were
accepted before advancing to the next integration time point. By the same checking
routine, the
size of the time step was increased when the rate of change voltage with time was
small. The
problem with this method was that, although it was excellent for accurate calculation of
the
SHAPE of the action potential, it was very poor near threshold. Near threshold there are
LONG
LATENCIES as the contest between nearly equal opposing are played out, and their
method
allowed much too large time steps to be taken and lead to the above noted error in
threshold. This problem is illustrated by our figure below where the generation of a spike depended on the integration
step size for the highly regarded Runga-Kutta method.
Furthermore, for a given time step, the generation of the spike sensitive to the
integration method as shown here.
Thus the apparent threshold was sensitive to BOTH the time step and the integration method.
Much attention and care has been paid to these problems in the development of NEURON. In order to minimize such errors as much as possible, Mike Hines made an extensive evaluation of numerical methods. He has incorporated some the most advanced numerical methods available in the design of NEURON to take much of the burden off of the user. Nevertheless the modeler should always be aware of these possible errors and evaluate their magnitude by sensitivity studies to the size of compartmental discretion and integration time step as well as using an alternative integration method.