Commentary: Accelerating spiking neural network simulations with PymoNNto and PymoNNtorch

Frontiers in Neuroinformatics, vol. 18, pp. 1446620

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) Stationary membrane potential distributions q (∞) (v) for exact integration (simulation: orange, Markov analyis: red) and forward Euler (simulation: light blue, Markov analysis: dark blue). Simulation data are for a single simulation of 1000 s duration; bin width 0.01 mV. The dotted lines in the bottom-left inset mark α = 0.952 for exact integration and α = 1.0 for forward Euler, respectively, indicating the difference in the first update step after reset to v = 0 mV. (B) Relative difference in membrane potential distributions between exact integration and forward Euler. The peak near 1 mV is due to the difference in the first step after reset. (C) Population activtity (1 ms bins) for simulations with NEST (orange) and Brian2 using delta synapses (blue: forward Euler, purple: exponential Euler); data from a single simulation with 10000 neurons.It will converge to the stationary membrane potential distribution q (∞) (v) under reasonable assumptions about p(v ′ |v) (Lasota and Mackey, 1994). The firing probability for an interval h is given by the integral term in Eq. 14, yielding the steady-state firing rater (∞) = h θ -∞ [p(0|v) -p(0|v)] q (∞) (v)dv . (16)PlesserSince no closed-form solution is available for these equations, we discretize the transition probability as a Markov matrix M . This matrix has a single eigenvector with eigenvalue 1 which corresponds to q (∞) (v)up to normalization (Feller, 1970;von Mises, 1964).In the absence of recurrent connections, Fig. 1A shows excellent agreement between simulations and Markov analyis (top) for exact integration and forward Euler, respectively, while Fig. 1B displays clear differences in the membrane potential distributions between the two models. We obtain from the Markov analysis for exact integration r (∞) = 9.93 sp/s vs 9.96 ± 0.04 sp/s from simulations and for forward Euler When recurrent connections are included and simulations performed using the original authors' code (saeedark and gitmv, 2024, git hash c9c59f0), simulation results differ from results obtained with NEST even if using the exponential Euler method in Brian2. Analysis of their code (Brian LIF.py) revealed that this neuron model used exponential post-synaptic currents instead of delta pulses. Once this was corrected, results obtained with Brian2 with exponential integration were consistent with NEST results, as shown in Fig. 1C.Careful analysis of STDP in a network of 1000 neurons revealed errors in both NEST and Brian2. NEST missed about 10% of weight increases, which could be traced to an error in NESTML code generation which is being fixed2 . For Brian2, a small number of synapses experienced up to seven STDP events, while at most three are expected. Details are given in the Supplementary Material.Based on the exact solution (Rotter and Diesmann, 1999) to the model defined by Eqs. 1-2, the analysis above confirms that the discrepancy in firing rates between NEST and the other simulation codes observed by Vieth et al. (2024) is due to the use of the unsuitable forward Euler method in the latter codes: using this method to integrate the membrane potential evolution actually changes the model under study. Indeed, inserting α and β for exact integration into Eqs. 7-8 and solving for τ and C shows that the forward Euler method solved the model for τ ≈ 10.5 ms and C ≈ 1.05 pF. Only NEST, the "odd one out", actually simulated the model as defined mathematically.Interestingly, the investigation of the discrepancy observed by Vieth et al. (2024) also helped us to uncover a subtle bug in NESTML. This shows that careful comparison of results by different simulators provides mutual benefits.The observed difference in firing rates of around 10% most likely will not have more than a 10% effect on the simulation runtimes measured by Vieth et al. (2024) and thus will not affect their overall conclusions.Nonetheless, I find it important to point out that benchmarks are meaningless unless one defines a desired level of precision (Morrison et al., 2007) and confirms that the simulation codes compared produce statistically equivalent results (see, e.g,. Van Albada et al., 2018).Finally, I would like to question the choice of benchmark model as such. The LIF neuron network model defined by Vieth et al. (2024) is peculiar in that it is a purely excitatory all-to-all network with a plasticity rule that only allows synapses to become stronger. The model can therefore only be benchmarked over a limited simulation time before run-away excitation leads to excessive spike rates. Balanced networks such as the two-population model by Brunel (2000) can, in contrast, be simulated for arbitrary periodsCommentary: Accelerating with PymoNNto of time and usually in different dynamic regimes. Furthermore, comparability with other studies might benefit from using the cortical microcircuit model by Potjans and Diesmann (2014) as a benchmark case, since this model has become widely used as a reference model in recent years (Knight and Nowotny, 2018;Van Albada et al., 2018;Golosio et al., 2021;Kauth et al., 2023).

Authors 1

1. Hans Ekkehard Pleßer corresponding Aachen Käte Hamburger Kolleg "Cultures of Research" (Core)

   RWTH Aachen University · Forschungszentrum Jülich · Norwegian University of Life Sciences

   Affiliation as printed

   Department of Data Science, Faculty of Science and Technology, Norwegian University of Life Sciences, Ås, Norway

   Institute for Advanced Simulation (IAS-6), Research Centre Jülich, Jülich, Germany

   Käte Hamburger Kolleg: Cultures of Research (c:o/re), RWTH Aachen University, Aachen, Germany