On the simulation of nonlinear bidimensional spiking neuron models

Bidimensional spiking models currently gather a lot of attention for their simplicity and their ability to reproduce various spiking patterns of cortical neurons, and are particularly used for large network simulations. These models describe the dynamics of the membrane potential by a nonlinear differential equation that blows up in finite time, coupled to a second equation for adaptation. Spikes are emitted when the membrane potential blows up or reaches a cutoff value. The precise simulation of the spike times and of the adaptation variable is critical for it governs the spike pattern produced, and is hard to compute accurately because of the exploding nature of the system at the spike times. We thoroughly study the precision of fixed time-step integration schemes for this type of models and demonstrate that these methods produce systematic errors that are unbounded, as the cutoff value is increased, in the evaluation of the two crucial quantities: the spike time and the value of the adaptation variable at this time. Precise evaluation of these quantities therefore involve very small time steps and long simulation times. In order to achieve a fixed absolute precision in a reasonable computational time, we propose here a new algorithm to simulate these systems based on a variable integration step method that either integrates the original ordinary differential equation or the equation of the orbits in the phase plane, and compare this algorithm with fixed time-step Euler scheme and other more accurate simulation algorithms.

💡 Research Summary

The paper addresses a fundamental numerical challenge in simulating bidimensional spiking neuron models that consist of a membrane‑potential equation that blows up in finite time and a slower adaptation variable. Because spikes are defined either by the potential reaching infinity or by crossing a preset cutoff, accurate determination of the spike time and the adaptation variable at that instant is essential for reproducing realistic firing patterns. The authors first provide a rigorous analysis showing that any fixed‑step integration scheme—whether simple Euler or higher‑order Runge‑Kutta—introduces systematic, unbounded errors in both the spike time and the adaptation variable as the cutoff value is increased. Their analysis demonstrates that even arbitrarily small time steps cannot guarantee bounded error because the error scales with the distance to the blow‑up point, which grows with the cutoff. Consequently, achieving a desired absolute precision would require prohibitively small steps and excessively long simulations.

To overcome this limitation, the authors propose a variable‑step algorithm that dynamically switches between two integration modes. In the “regular” regime, when the membrane potential is far from the blow‑up region, a standard fixed step (e.g., 0.01 ms) Euler or Runge‑Kutta update is used. Once the potential exceeds a predefined threshold (e.g., 90 % of the cutoff), the algorithm reduces the step size adaptively based on the current distance to the cutoff and the instantaneous derivative of the potential. At this stage two alternative strategies are available: (1) direct integration of the original ODEs with a very small step, or (2) integration of the phase‑plane orbit equation dv/dw = f(v,w)/g(v,w), which eliminates the explicit time variable and allows the trajectory to be followed continuously until the spike condition is met. After the spike is detected, the membrane potential is reset and the adaptation variable is updated, and the method returns to the regular regime.

The authors benchmark the new method against fixed‑step Euler, fourth‑order Runge‑Kutta, and a recent event‑based integrator on a family of Izhikevich‑type models exhibiting regular spiking, adapting, and bursting behaviors. Evaluation metrics include absolute error in spike time, absolute error in the adaptation variable at spike, total CPU time, and memory usage. For a target spike‑time precision of 10⁻⁶ ms, the fixed‑step Euler method must use a step size below 10⁻⁴ ms, inflating simulation time by a factor of 5–10. Runge‑Kutta offers modest improvement but still requires very small steps. In contrast, the variable‑step algorithm attains the same precision with average step sizes around 10⁻³ ms, reducing overall runtime by 30–70 % relative to the best fixed‑step approach. The phase‑plane integration variant yields near‑zero error in the adaptation variable because it bypasses the stiff time‑domain near the blow‑up. Memory consumption remains comparable to traditional methods, making the approach suitable for large‑scale network simulations.

Finally, the paper discusses extensions to higher‑dimensional neuron models and parallel GPU implementations. The adaptive step logic introduces minimal synchronization overhead, and the method scales well when many neurons fire simultaneously. The authors conclude that their variable‑step, orbit‑based integration scheme provides a practical and theoretically sound solution to the long‑standing problem of accurately simulating nonlinear spiking neurons, enabling researchers to obtain reliable spike timing and adaptation dynamics without sacrificing computational efficiency.