Voltage-stepping schemes for the simulation of spiking neural networks

Abstract

The numerical simulation of spiking neural networks requires particular attention. On the one hand, time-stepping methods are generic but they are prone to numerical errors and need specific treatments to deal with the discontinuities of integrate-and-fire models. On the other hand, event-driven methods are more precise but they are restricted to a limited class of neuron models. We present here a voltage-stepping scheme that combines the advantages of these two approaches and consists of a discretization of the voltage state-space. The numerical simulation is reduced to a local event-driven method that induces an implicit activity-dependent time discretization (time-steps automatically increase when the neuron is slowly varying). We show analytically that such a scheme leads to a high-order algorithm so that it accurately approximates the neuronal dynamics. The voltage-stepping method is generic and can be used to simulate any kind of neuron models. We illustrate it on nonlinear integrate-and-fire models and show that it outperforms time-stepping schemes of Runge-Kutta type in terms of simulation time and accuracy.

Appendix: Order of voltage-stepping schemes

For simplicity, we consider a general neuron model described as follows (notations are defined in Sections 2 and 3)

$$ \frac{dv(t)}{dt}=f\left( v(t)\right) \label{org_sys} $$

(13)

Over V _i , one possible linear interpolation of the function f(v) is achieved by interpolating it at the boundaries v _i and v _i + 1. Other possible linear interpolations will be discussed in the second part of this appendix.

1.1 Appendix A1: Linear interpolation at boundaries (VS2 method)

Let us consider v _Δv the solution of the following dynamical system:

$$ \frac{dv}{dt}=f_{\Delta v}(v), \label{app_sys} $$

(14)

where f _Δv is the piecewise linear function defined as follows:

$$ f_{\Delta v}(v)=\frac{v_{i+1}-v}{\Delta v}f\left( v_{i}\right) +\frac{v-v_{i} }{\Delta v}f\left( v_{i+1}\right) $$

It is straightforward to show that the error of approximation is given by

$$ \left\vert f\left( v\right) -f_{\Delta v}\left( v\right) \right\vert =O\left( \Delta v^{2}\right) \label{boud} $$

(15)

Considering the following theorem:

Theorem 1

Fundamental Inequality (see for instance Hubbard and West (1991)). For a differential equation $\dot{x}=F(x)$ satisfying the Lipschitz condition with K ≠ 0 and if u ₁(t) and u ₂(t) are two continuous, piecewise differentiable functions satisfying $\left\vert \dot{u}_{i}(t)-F(u_{i}(t)) \right\vert \leq \varepsilon _{i}$ for all t at which u ₁(t) and u ₂(t) are differentiable and if $\left\vert u_{1}(0)-u_{2}(0)\right\vert \leq \delta ,$ then

$$ \left\vert u_{1}(t)-u_{2}(t)\right\vert \leq \delta e^{K\left\vert t\right\vert }+\frac{\varepsilon _{1}+\varepsilon _{2}}{K}\left( e^{K\left\vert t\right\vert }-1\right) . $$

Applying Theorem 1 to Eqs. (13), (14) and using Eq. (15), it can be proved that $\vert v-v_{\Delta v} \vert =O \left( \Delta v^{2}\right)$. It follows $\left\vert t_{ex}^f-t_{ap}^f\right\vert=O\left( \Delta v^{2}\right)$ that means that the estimate error on the exact spike time is of order $O\left( \Delta v^{2}\right)$.

Remarks

• At the neural network level, the incoming spikes generated by presynaptic neurons introduced a second order error (since $\left\vert t_{ex}^{f}-t_{ap}^{f}\right\vert =O(\Delta v^{2})$). Noting the fact that f _Δv also introduced a second-order error, the proposed voltage-stepping scheme (VS2) guarantees the same accuracy at the network level as the neuron level, even after considering the effect of propagation of error on spike times.

• For the p-dimensional case, the only difference is to approximate f(v) over $V_{i}\subseteq R^{p}$ by a linear system: $f_{\Delta v}\left( v\right) =A_{i}v+b_{i}$, which is uniquely determined as the linear interpolation vector field of f(v). The same result can be proved using a norm on $\mathcal{R}^p$, $\left\Vert \cdot \right\Vert$, instead of the absolute value $\left\vert \cdot \right\vert$ ( see Girard 2002 for more details).

1.2 Appendix A2: Linear interpolation at gaussian abscissas (VS4 method)

Let us consider (Eq. (13)) and (Eq. (14)) over a voltage interval V _i . Without loss of generality, assume that the voltage interval V _i + 1 is reached. We have

$$ \Delta t_{i}=t_{ex}^{i}-t_{ap}^{i}=\int\nolimits_{v_{i}}^{v_{i+1}}\left( \frac{1}{f\left( v\right) }-\frac{1}{f_{\Delta v}\left( v\right) }\right) dv \label{intGs} $$

(16)

where $t_{ex}^{i}$ represents the exact exit time of V _i , and $t_{ap}^{i}$ represents its approximation. The best choice for f _Δv is those that minimize (Eq. (16)). We have $1/f_{\Delta v}(v) \in C^{\infty }$ (almost everywhere). For 1/f(v) ∈ C ^k, k ≥ 4, and according to Gaussian quadrature rule, the linear interpolation at gaussian abscissas f _Δv satisfies

$$ \int\nolimits_{v_{i}}^{v_{i+1}}\left( \frac{1}{f\left( v\right) }-\frac{1}{ f_{\Delta v}\left( v\right) }\right) dv=O\left( \Delta v^{5}\right) \label{intCond} $$

(17)

The only point is to calculate the gaussian abscissas over V _i . Over $V_{i}=\left[ v_{i},v_{i+1}\right[$, the gaussian abscissas are:

$$\begin{array}{lll} v_{i,1} &=&\frac{\sqrt{3}-1}{2\sqrt{3}}v_{i+1}+\frac{\sqrt{3}+1}{2\sqrt{3}} v_{i} \\ v_{i,2} &=&\frac{\sqrt{3}+1}{2\sqrt{3}}v_{i+1}+\frac{\sqrt{3}-1}{2\sqrt{3}} v_{i} \end{array}$$

based on which, the linear interpolation of f(v) can be described as follows:

$$ f_{\Delta v}(v)=\frac{v_{i,2}-v}{v_{i,2}-v_{i,1}}f\left( v_{i,1}\right) + \frac{v-v_{i,1}}{v_{i,2}-v_{i,1}}f\left( v_{i,2}\right) $$

Over each V _i , the local error (approximation of exit time) is of order $O\left( \Delta v^{5}\right)$. The estimate error on the exact spike time is obtained considering the exit times over the entire voltage-space that gives a global error of $O\left( \Delta v^{4}\right)$. This error estimate agrees exactly with the results of the numerical simulations.

It should be noted that this method requires one-dimensional neural models. The major reason is that Eq. (17) cannot be always fulfilled in high dimensional case. Moreover the methods also failed for neural network simulation. Assume that we can estimate the incoming spikes generated by presynaptic neurons to an accuracy of $O\left( \Delta v^{4}\right)$. Therefore a fourth order error is introduced in Eq. (17) and the error can be no better than $O\left( \Delta v^{4}\right)$ which makes impossible to calculate the exit time with a fifth-order accuracy.