Numerical Integration

The goal of the present simulation is to illustrate how the user should select the proper time step (i.e., keyword TIMING:, Step-size (h), in the *.smu file) for the numerical integration.

In the current version of SNNAP, the ordinary differential equations (ODEs) are integrated with the forward-Euler rule, which uses a fixed time-step to numerically solve ODEs. The forward-Euler rule is fast and accurate, but it can be stable if the time step is too large. In any given simulation, the integration time step must be small enough to ensure stability of the fastest process (i.e., the smallest time constant).

Four simulations are included in this example. All simulations are of the same neural network and cells. The network is the three cell oscillator that is also in the \examples\Neural_Networks subdirectory. The four simulations differ only in the size of the time step used in the numerical integration. The four simulations are:

osc 1000

This simulation uses a time step of 1 ms (i.e., 1000 us). The results of this simulation are illustrated in osc1000.gif. This is a clear example of a simulation that has become catastrophically unstable.

osc 450

This simulation uses a time step of 450 us. The results of this simulation are illustrated in osc450.gif. This simulation represents a potentially misleading result. At first glance the user may think the results are reasonable; i.e., the pattern of activity is similar to the anticipated result. Closer inspection, however, reveals that the simulation was unstable. First, note that the “after hyperpolarization” of the spike extends well beyond the reversal potential for the potassium current (i.e., -75 mV), which is an non-physiological result. In addition, note that the “kinetics” of the “after hyperpolarization” are instaneous, which also non-physiological. Finally, note that the peak of the “spike” extends beyond the reversal potential for the sodium current (i.e., 76 mV), which is an non-physiological result.

osc 45

This simulation uses a time step of 45 us. The results of this simulation are illustrated in osc45.gif. The results of this simulation appear to be accurate and stable. The waveform of the spikes is physiologically reasonable and evolution of the simulation is smooth (i.e., no instaneous jumps to new values). Note the different number of spikes in the osc_450 and osc_45 simulations. The spike threshold is a highly non linear region and small inaccuracies / instabilities can have a large effect at this point. Thus, as the user changes the size of the time step, the number of spikes may also change.

osc 25

This simulation uses a time step of 25 us. The results of this simulation are illustrated in osc25.gif. The results of this simulation and those of osc_45.smu are identical, which indicates that an acceptable time step is 45 us (i.e., further reduces do not produce additional changes in the outcome of the simulation).

The user should have also noticed that as the size of the time step was reduced the speed of the simulation became slower. This is simply because more and more iterations were required as the time step became smaller. Thus, selecting an appropriate time step is a balancing act between accuracy, stability and speed.

As a general rule the user should select a time step that is about ten times faster than the smallest time constant in the simulation. For example, the activation time constant for the sodium current is generally the fastest process in a HH-type model. In the osc.smu model, this time constant is 450 us. Thus, a time step of 45 us should (and did) yield stable/accurate results.