Commentary: Accelerating spiking neural network simulations with PymoNNto and PymoNNtorch

Author:
Open Access
Publisher website
Google Scholar
Abstract
In a recent paper introducing the PymoNNto and PymoNNtorch simulators, Vieth et al. (2024) compare their own simulators to several publicly available simulation codes using two different neuronal network models. They report that "the native NEST implementation consistently generates slightly fewer spikes in comparison to the other simulators" for their LIF neuron network model (see their Fig. S2), but do not investigate this issue beyond confirming that the discrepancy is not caused by the use of double precision numerics in NEST.This commentary shows that the observed difference in firing rate is due to an unsuitable numerical integration scheme used in the other simulators while NEST creates the correct model dynamics. Vieth et al. (2024) study a weakly-coupled all-to-all network of K excitatory neurons driven by individual random input. The plasticity mechanism in the network is not relevant to the numerical issues investigated here and thus left out. Generalizing the equations of the original paper slightly, the membrane potential of an individual neuron k is given byv(k) (t) = - v (k) (t) τ + I (k) ext (t) + I (k) syn (t) C (1)with membrane time constant τ and capacitance C. The membrane potential is reset to v (k) = 0 when v (k) = θ and there is no refractory time. The external input current defined byİ(k) ext (t) = (U[I min , I max ) -I) δ(t -jh) (2)is a piecewise constant current with a new amplitude chosen from a uniform distribution at fixed time intervals h. The recurrent synaptic input isI (k) syn (t) = l,i w kl δ(t -t (l,i) -d)(3)with uniformly distributed synaptic weights w kl ∼ U[0, 1/K) mV, the synaptic delay d, and t (l,i) representing the time of the i th spike fired by neuron l.For the parameters used by Vieth et al. (2024) 1 , the mean input current is µ = 1 2 pA with variance σ 2 = 1 12 pA 2 , and the average synaptic weight w = 1 2K mV = 5 • 10 -5 mV. In this regime, spiking activity is mainly driven by external input with an average firing rate between 100 ms and 295 ms of 11.6 sp/s for the full model and of 9.9 sp/s without any connections.We consider first an indvidual neuron driven by external noise only. Dropping unnecessary subscripts and discretizing the dynamics on the fixed time grid t j = jh at which the input current takes on new random values I j , Eq. 1 becomesv j+1 = βv j + αI j .(4)The coefficients α and β depend on the numerical integration method chosen.Eq. 4 will provide the exact solution of Eqs. 1-2 if we use the exact integration method (Rotter and Diesmann, 1999) for whichα = τ C 1 -e -h/τ ≈ 0.952(5)β = e -h/τ ≈ 0.905 .In the specific case of piecewise constant input over the integration step as in the given model, this method is equivalent to the exponential integration by MacGregor (1987, Ch. 14.C.5). We will use this solution as reference.The coefficients for the forward Euler method areα = h C = 1.0 (7) β = 1 - h τ = 0.9 . (8)For first-passage time problems for the process defined by Eq. 4 in combination with a spiking threshold, closed-form solutions exist for noise symmetric to the origin (Larralde, 2004), while only bounds are known for more general cases as studied here (Novikov and Kordzakhia, 2008).1 K = 10 4 , τ = 10 ms, C = 1 pF, h = d = 1 ms, I min = 0 pA, Imax = 1 pA, θ = 6 mVWe thus consider the free membrane potential in the absence of a threshold. Assuming v 0 = 0, we obtain by repeated application of the update equationv j+1 = α j k=0 β k I k ,(9)where we have renumbered the random currents I k , exploiting their independence. Averaging over currents and taking the limit j → ∞, we obtain the mean membrane potential and its standard deviation⟨v⟩ = α 1 -β µ = τ C µ = 5 mV (10)⟨∆v 2 ⟩ = α 1 -β 2 σ 2 =    τ σ C 1-e -h/τ 1+e -h/τ ≈ 0.645 mV (exact integration) τ σ C h 2τ -h ≈ 0.662 mV (forward Euler) .(11)The broader distribution of the free membrane potential in the case of forward Euler integration (8.1% more probability mass above threshold θ) suggests that neurons will fire more frequently. This is supported by simulations below.We confirm this observation quantitatively through numerical Markov analysis. This analysis is in principle continuous in time, but since the noise in the network model switches in intervals h, we consider only time points at which the noise changes. Given v k , the smallest possible v k+1 is obtained for I = I min during the time step and the largest for I = I max . In between the maximum and minimum possible values, any value of v k+1 is attained with equal probability. We can thus write the free transition probability for the membrane potential asp(v ′ |v) =      0 if v ′ < αI min + βv ⇔ v > (v ′ -αI min )/β 1 τ (Imax-I min ) else 0 if v ′ > αI max + βv ⇔ v < (v ′ -αI max )/β (12) = 1 τ (I max -I min ) Θ v - v ′ -αI min β -Θ v - v ′ -αI max β(13)where Θ(x) is the Heaviside step function. To include the effect of the spiking threshold θ, we definep(v ′ |v) =      p(0|v) + ∞ θ p(v ′ |v)dv ′ v ′ = 0 0 v ′ > θ p(v ′ |v) else. (14)Here, the integral in the first clause describes neurons re-inserted at the reset potential after spiking, and the corresponding probability is removed for superthreshold v ′ by the second clause.If q k (v) is the membrane potential distribution at time step k, then the distribution at step k + 1 is given by q k+1 (v ′ ) = ∞ -∞ p(v ′ |v)q k (v)dv . (15)(A)...

This publication has 13 references indexed in Scilit: