Spike-Aware Propagation Approximation for Conductance-Based LIF Equations

1. Introduction

Large-scale spiking neural network (SNN) simulation connects cellular dynamics with collective circuit behavior. Reduced point-neuron models make such studies computationally tractable, and the leaky integrate-and-fire (LIF) description offers a useful compromise between biophysical detail and numerical efficiency [1,2]. Conductance-based LIF models additionally retain receptor-specific decay and voltage-dependent synaptic driving forces that are absent from purely current-based descriptions [3]. These models support studies ranging from canonical recurrent circuits to data-constrained cortical and brain-scale systems [4,5,6,7,8,9]. As network scale and biological scope increase, the neuronal update kernel is executed ever more frequently. Its numerical accuracy and efficiency therefore place direct limits on feasible simulation duration, parameter exploration, and network size.

The principal difficulty is that a conductance-based LIF network is a hybrid dynamical system, not merely a smooth ordinary differential equation (ODE). Receptor traces decay continuously between synaptic arrivals and jump when spikes arrive. The membrane voltage follows a nonautonomous linear ODE whose coefficients depend on those traces, while threshold crossings trigger resets and refractory clamps. Fine temporal resolution can track these mechanisms but requires repeated updates even when much of the state evolves predictably. Coarser resolution reduces this workload but may displace or miss a threshold crossing and delay its downstream effect. Because event timing changes later resets and synaptic inputs, a local error can propagate through the recurrent network. An effective integration strategy must therefore balance continuous voltage accuracy, discrete event timing, and computational efficiency.

Modern simulators improve throughput through parallel execution, code generation, graphics processing units (GPUs), and communication-aware spike delivery [10,11,12,13,14,15]. The accuracy available at a chosen temporal resolution nevertheless remains controlled by the neuronal integration and spike-detection scheme. Exact propagation is possible for time-invariant linear subthreshold systems [16], while off-grid methods improve spike timing within globally time-driven simulations [17,18]. Event-driven and retrospective tests provide accurate threshold handling for selected point-neuron models [19,20], and analytical conductance-based schemes exploit additional model structure [21,22]. These developments solve complementary parts of the problem, but combining exact receptor-trace updates, economical propagation under time-dependent conductance, and spike localization remains challenging in large recurrent networks.

Here we develop a spike-aware propagation (SAP) approximation for this combined setting. SAP updates visible receptor traces exactly, applies an analytic homogeneous propagator to the membrane equation, and uses Gauss–Legendre quadrature only for the remaining integral. A bisection is then applied to localize a spike when there is a candidate threshold crossing within a scan interval, which separates the roles of ODE updating width and spike-localization tolerance. The analysis establishes an error bound and a conditional convergence result for a prescribed visible-arrival schedule and fixed event itinerary. Large-scale comparisons use sampled-voltage error relative to a fine NEST reference and neuronal ODE-solving time to identify the accuracy-efficiency performance of Euler, SAP, and NEST. The SAP ablations then determine the sensitivity of its main hyperparameters. Together, the analysis and experiments show that SAP can improve the accuracy-efficiency Pareto frontier of conductance-based LIF simulation, highlighting the practical value in large-scale spiking neural network applications.

2. Related Work

2.1. Large-Scale Spiking Network Models

Large-scale SNN models have progressively increased both biological scope and data dependence. Early large recurrent and thalamocortical models established how reduced spiking neurons can reproduce collective regimes that are difficult to study at the single-cell level [2,23]. Subsequent work integrated cell-type-specific connectivity into a full-scale cortical microcircuit [4], reconstructed neocortical microcircuitry at greater cellular detail [5], and linked local microcircuits across multiple macaque cortical areas [6]. Data-constrained mouse visual-cortex models extended this program by jointly using structural and functional measurements [7]. More recent brain-computing studies have pursued human-scale simulation and assimilation on large GPU systems [8,9]. These studies motivate numerical methods that remain economical when the same point-neuron update is executed many billions of times. The present work does not propose a new biological network architecture; it focuses on the integration kernel used to advance conductance-based LIF neurons within such large recurrent simulations.

2.2. Simulation Platforms and Scalability

General-purpose neural simulators balance model expressiveness, reproducibility, and performance in different ways [24]. NEURON and its optimized CoreNEURON engine support detailed neuronal models [25,26], whereas NEST targets large networks of point neurons with shared-memory and distributed parallelism [10]. Brian and Brian 2 use high-level model definitions and code generation to combine flexible specification with compiled execution [11,27]. NESTML extends this code-generation approach to portable descriptions of neuron and synapse models and makes the separation between model equations, event handling, and numerical integration explicit [28].

Hardware-oriented work has reduced the cost of both neuronal updates and spike delivery. GeNN generates accelerator-specific code for SNN simulation [12], while Brian2CUDA supplies a GPU backend for the Brian 2 model language [13]. Procedural connectivity reduces GPU memory pressure by generating connections when required [14], and NEST GPU studies have demonstrated both single-GPU and multi-GPU acceleration [29,30]. On distributed-memory systems, communication and data-structure optimizations have supported petascale and prospective exascale simulation [15,31,32]. Direct NEST–GeNN benchmarks further show that the preferable platform depends on model size, connectivity, hardware, and setup costs [33]. This body of work addresses where and how simulation workloads are executed. SAP is complementary: it changes the per-neuron propagation and threshold-localization procedure while retaining a shared scan schedule suited to parallel execution.

2.3. Time-Driven, Event-Driven, and Exact Integration

Numerical strategies for SNNs are commonly organized by how continuous flows and discrete spike events are scheduled. Globally time-driven methods update neuronal states on a shared grid and are straightforward to parallelize, but a grid-constrained threshold test quantizes spike times. Globally event-driven methods advance from one event to the next and can avoid inactive updates. However, they require reliable prediction of the next event and may introduce model-specific data structures. These are not mutually exclusive categories: hybrid schemes combine global time stepping with exact or off-grid operations inside each step [24].

For linear, time-invariant subthreshold dynamics, matrix-exponential propagation gives an exact grid-to-grid update [16]. Morrison et al. combined exact subthreshold integration with continuous spike times in a discrete-time network schedule [17], and Hanuschkin et al. extended precise spike localization to nonlinear point-neuron models by iterative localization [18]. Event-driven algorithms have also been developed for nonlinear integrate-and-fire dynamics [19]. When a trajectory can cross and return below threshold between grid points, retrospective state-space tests can guarantee detection for selected linear subthreshold models [20]. Collectively, these methods show that exact propagation and accurate event timing can be combined with scalable scheduling. Their guarantees nevertheless depend on the subthreshold dynamics and on how incoming events are exposed to the neuron.

2.4. Conductance-Based Analytical Approximations

Conductance-based synapses are more difficult than current-based synapses because the synaptic input multiplies the membrane voltage. Brette developed an exact event-driven treatment for an integrate-and-fire model with exponential synaptic conductances under specific analytic restrictions [21]. Rudolph and Destexhe derived analytical conductance-based LIF models suited to event-driven strategies [22]. Related exact methods for exponential synaptic currents permit multiple time constants but apply to a different, current-based coupling structure [34]. Continued work on reduced conductance models, including simplified NMDA-receptor dynamics, illustrates the ongoing need to balance receptor detail and tractability [35].

SAP takes a different approximation route. It preserves the exponential decay and jumps of each visible receptor trace, so no numerical ODE solver is required for those variables. The resulting membrane equation is treated by variation of constants: the homogeneous factor is evaluated analytically, and only the intractable forcing integral is approximated. This construction draws on the theory of exponential propagation and high-order Gaussian quadrature [36,37,38]. Threshold crossings and visible arrivals are then handled as explicit discontinuities, consistent with the general principle that a smooth solver must be coupled to event treatment for discontinuous ODE problems [39]. The distinction from the earlier conductance-based analytical schemes is therefore not that SAP makes the full hybrid trajectory exact. Rather, it isolates the approximation to a smooth integral and spike localization, which permits a direct decomposition of the resulting errors while retaining a common scan schedule for parallel execution.

3. Methods

3.1. Conductance-Based LIF Equations

For neuron i, let $V_i\left(t\right)$ denote the membrane potential, $C_i$ the membrane capacitance, and $g_{L,i}$ the leak conductance. The leak, threshold, and reset potentials are $V_L$, $V_{\mathrm{th},i}$, and $V_{\mathrm{reset},i}$, respectively, and $T_{\mathrm{ref},i}$ is the absolute refractory period. The active receptor set is $\mathcal{U}=\{\mathrm{AMPA},{\mathrm{GABA}}_A\}$. For $u\in \mathcal{U}$, the receptor trace $J_i^u\left(t\right)$ has decay time ${\tau }_i^u$, conductance scale $g_i^u$, and reversal potential $V^u$. The synaptic weight from neuron j is $w_{ij}^u$, and the kth corresponding arrival occurs at $t_{ij}^{k,u}$. With constant external drive $I_i$, the network dynamics are written in area-normalized units: time is measured in $\mathrm{ms}$, voltage in $\mathrm{mV}$, capacitance density in $\mu \mathrm{F}{\mathrm{cm}}^{-2}$, conductance density in $\mathrm{mS}{\mathrm{cm}}^{-2}$, and current density in $\mu \mathrm{A}{\mathrm{cm}}^{-2}$. The traces $J_i^u$ and multipliers $w_{ij}^u$ are dimensionless.

$$\left\{\begin{array}{cc}\hfill C_i{\displaystyle \frac{\mathrm{d}{V}_i\left(t\right)}{\mathrm{d}t}}& =-\left(g_{L,i}+\sum _{u\in \mathcal{U}}{g}_i^u{J}_i^u\left(t\right)\right)V_i\left(t\right)+g_{L,i}{V}_L+\sum _{u\in \mathcal{U}}{g}_i^u{V}^u{J}_i^u\left(t\right)+I_i,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}{J}_i^u\left(t\right)}{\mathrm{d}t}}& =-{\displaystyle \frac{J_i^u\left(t\right)}{{\tau }_i^u}}+\sum _{j,k}{w}_{ij}^u\delta (t-t_{ij}^{k,u}),u\in \mathcal{U},\hfill \\ \hfill t_i^k& =inf\{t>t_i^{k-1}+T_{\mathrm{ref},i}:V_i\left(t\right)\ge V_{\mathrm{th},i}\},\hfill \\ \hfill V_i\left(t\right)& =V_{\mathrm{reset},i},t\in [t_i^k,t_i^k+T_{\mathrm{ref},i}].\hfill \end{array}\right.$$

(1)

The membrane equation in (1) is equivalently expressed through the synaptic current $I_i^u\left(t\right)=g_i^u(V^u-V_i\left(t\right))J_i^u\left(t\right)$.

Let $[T_0,T_1]$ denote a batch-delivery interval of length $H=T_1-T_0$. The state at $T_0$ comprises the membrane potentials, receptor traces, refractory variables, and all presynaptic arrivals visible under the prescribed schedule. SAP scans the interval through common bins $[a,b]$ of width h, with a possibly shorter final bin. Although all neurons share this scan schedule, propagation for an individual neuron may begin within a bin when its refractory period ends. Under batched spike delivery, spikes emitted in $[T_0,T_1]$ enter the postsynaptic traces at $T_1$. Hence H determines the delivery cadence, whereas h determines the grid used to identify candidate spikes.

3.2. Spike-Aware Propagation Approximation

Consider a predictor interval $[t_0,t]\subseteq [T_0,T_1]$, with $t>t_0$, on which the dynamics are smooth between visible arrivals. For $s\in [t_0,t]$, define

$$p_i\left(s\right)=g_{L,i}+\sum _{u\in \mathcal{U}}{g}_i^u{J}_i^u\left(s\right),q_i\left(s\right)=g_{L,i}{V}_L+\sum _{u\in \mathcal{U}}{g}_i^u{V}^u{J}_i^u\left(s\right)+I_i.$$

(2)

These coefficients reduce the subthreshold voltage dynamics to a scalar nonautonomous linear equation. Once the visible arrival set is fixed, the receptor traces can be updated exactly. Let ${\mathcal{A}}_i^u(t_0,t)=\{(j,k):t_0<t_{ij}^{k,u}\le t\}$ contain the index pairs for arrivals visible to neuron i through receptor u after the initial trace $J_i^u\left(t_0\right)$ has been formed. Then

$$J_i^u\left(t\right)=J_i^u\left(t_0\right)e^{-\eta /{\tau }_i^u}+\sum _{(j,k)\in {\mathcal{A}}_i^u(t_0,t)}{w}_{ij}^u{e}^{-(t-t_{ij}^{k,u})/{\tau }_i^u},\eta =t-t_0.$$

(3)

Equation (3) is exact for the prescribed visible arrivals. In the parallel implementation, spikes emitted during a batch are accumulated and become visible to postsynaptic traces at $T_1$. Propagation inside the batch therefore uses the visible-arrival set fixed at its left boundary.

Applying variation of constants to the membrane equation gives

$$V_i\left(t\right)=F_i(t,t_0)V_i\left(t_0\right)+{\int }_{t_0}^t{\displaystyle \frac{q_i\left(s\right)}{C_i}}{F}_i(t,s)ds,$$

(4)

where $F_i(t,s)=exp[-C_i^{-1}{\int }_s^t{p}_i\left(r\right)dr]$. On an arrival-free segment, $J_i^u\left(r\right)=J_i^u\left(s\right)e^{-(r-s)/{\tau }_i^u}$, and hence

$$F_i(t,s)=exp\left[-{\displaystyle \frac{g_{L,i}}{C_i}}(t-s)-{\displaystyle \frac{1}{C_i}}\sum _{u\in \mathcal{U}}{g}_i^u{\tau }_i^u{J}_i^u\left(s\right)\left(1-e^{-(t-s)/{\tau }_i^u}\right)\right].$$

Visible arrivals partition the predictor interval into arrival-free segments, on each of which this propagator is exact. SAP therefore evaluates the homogeneous contribution analytically. Numerical approximation is confined to the inhomogeneous integral and spike localization.

SAP approximates the remaining integral by an M-node Gauss–Legendre rule on $[t_0,t]$:

$${\widehat{V}}_i(t;t_0)=F_i(t,t_0)V_i\left(t_0\right)+{\displaystyle \frac{\eta }{2}}{\displaystyle \sum _{m=1}^M}{\omega }_m{\displaystyle \frac{q_i\left(s_m\right)}{C_i}}{F}_i(t,s_m),s_m={\displaystyle \frac{t_0+t}{2}}+{\displaystyle \frac{\eta }{2}}{\xi }_m.$$

(5)

Equation (5) approximates only the smooth forcing integral; receptor decay and homogeneous attenuation remain exact. The analysis assumes $0<h<{min}_i{T}_{\mathrm{ref},i}$ and the endpoint-screening condition of Theorem 1. On each active scan bin, the no-reset predictor must either remain subthreshold or possess one unique threshold zero followed by a suprathreshold right endpoint. The endpoint test then detects that crossing, while $h<T_{\mathrm{ref},i}$ prevents a second spike after the reset within the same bin. Within a bin $[a,b]$, neuron i remains clamped until $t_i^{\mathrm{eff}}=max\{a,t_i^{\mathrm{last}}+T_{\mathrm{ref},i}\}$. If $t_i^{\mathrm{eff}}<b$, Equation (5) propagates the neuron over $[t_i^{\mathrm{eff}},b]$ from traces decayed to $t_i^{\mathrm{eff}}$. Otherwise, the neuron remains clamped. A suprathreshold endpoint prediction marks a candidate spike, and the same predictor localizes its time within the bin. The localized spike time determines the next refractory release, and the bin-end voltage is reset. Because $h<T_{\mathrm{ref},i}$, the neuron cannot spike again in the unscanned remainder of that bin. Propagation resumes from the neuron-specific refractory release time in a later bin. Emitted spikes are accumulated and inserted into postsynaptic traces at $T_1$.

Algorithm 1 summarizes the resulting SAP workflow. For each scan bin, it constructs refractory-aware initial states, advances the neurons with Equation (5), localizes candidate spikes, and accumulates emitted spikes for delivery at the batch boundary.

Algorithm 1: SAP: spike-aware propagation approximation

Input:  States $V,J,t^{\mathrm{last}}$; neuronal parameters; final time T Batch length H; scan width $0<h<{min}_i{T}_{\mathrm{ref},i}$; quadrature order M; root tolerance ${\epsilon }_t$; iteration cap B Endpoint-screen completeness on every active bin, as specified in Theorem 1 Output:  Updated states $V,J,t^{\mathrm{last}}$, with emitted spikes delivered at batch boundaries   1: for each batch-delivery interval $[T_0,T_1]$ do   2:   Store $J\left(T_0\right)$ and initialize the outgoing counts $S\leftarrow 0$   3:   for each common scan bin $[a,b]\subseteq [T_0,T_1]$ of nominal width h do   3:   Construct refractory-aware initial states   4:     for all neurons i do   5:      $t_i^{\mathrm{eff}}\leftarrow max\{a,t_i^{\mathrm{last}}+T_{\mathrm{ref},i}\}$ and ${\Delta }_i\leftarrow max\{0,b-t_i^{\mathrm{eff}}\}$   6:      $V_i^0\leftarrow V_{\mathrm{reset},i}$ if $t_i^{\mathrm{eff}}>a$; otherwise, $V_i^0\leftarrow V_i$   7:      Decay $J_i^u\left(T_0\right)$ to $J_i^{u,0}=J_i^u\left(t_i^{\mathrm{eff}}\right)$ for every receptor type u   8:      $V_i^{\mathrm{pred}}\leftarrow {\widehat{V}}_i(b;t_i^{\mathrm{eff}})$ if ${\Delta }_i>0$; otherwise, $V_i^{\mathrm{pred}}\leftarrow V_{\mathrm{reset},i}$   9:     end for   Detect and localize within-bin spikes 10:     Set $\mathcal{C}\leftarrow \{i:{\Delta }_i>0,V_i^{\mathrm{pred}}\ge V_{\mathrm{th},i}\}$ 11:     for all $i\in \mathcal{C}$ do 12:      With the same predictor, localize ${\widehat{t}}_i\in [t_i^{\mathrm{eff}},b]$ to tolerance ${\epsilon }_t$ or iteration cap B 13:      Set $t_i^{\mathrm{last}}\leftarrow {\widehat{t}}_i$, $S_i\leftarrow S_i+1$, and $V_i^{\mathrm{pred}}\leftarrow V_{\mathrm{reset},i}$ 14:     end for 15:     Accept $V_i\leftarrow V_i^{\mathrm{pred}}$ for every neuron 16:   end for 17:   Decay visible traces to $T_1$ and deliver the accumulated spike batch S 18: end for

4. Error Analysis

4.1. One-Batch Voltage-Error and Conditional Convergence Analysis

To isolate the propagation and spike-localization errors, fix one neuron i and one batch interval $[T_0,T_1]$. Let $V_i$ denote the exact voltage trajectory and let ${\widehat{V}}_i$ denote the SAP trajectory. The two trajectories start from the same voltage, receptor traces, and refractory status at $T_0$, and they use the same prescribed visible-arrival data during the batch. Values at resets are taken right-continuously. A receptor trace jump at $T_1$ does not instantaneously change the voltage, so the endpoint error is evaluated at $T_1^{-}$. This construction keeps the theorem focused on the observable membrane voltage while retaining the event times needed to control reset and refractory effects.

Let ${\mathcal{P}}_i$ be the collection of active predictor intervals obtained by cutting the scan bins at visible arrivals and refractory releases. Thus every $K=[a_K,b_K]\in {\mathcal{P}}_i$ is smooth, $\left|K\right|\le h$, and ${\sum }_{K\in {\mathcal{P}}_i}\left|K\right|\le H$. For $a_K<x\le b_K$, define

$${\mathcal{G}}_{i,K}^{\left(2M\right)}=\underset{a_K<x\le b_K}{sup}{∥{\partial }_s^{2M}\left\{{\displaystyle \frac{q_i\left(s\right)}{C_i}}{F}_i(x,s)\right\}∥}_{L^{\infty }(a_K,x)}$$

(6)

and

$$d_{i,K}={\beta }_M{\left|K\right|}^{2M+1}{\mathcal{G}}_{i,K}^{\left(2M\right)},{\beta }_M={\displaystyle \frac{{(M!)}^4}{(2M+1){[\left(2M\right)!]}^3}}.$$

(7)

Theorem 1

(One-batch voltage error for spike-aware propagation). Fix $0<h<T_{\mathrm{ref},i}$ and an integer $B\ge 1$. Assume the following conditions on $[T_0,T_1]$.

1. 

The visible-arrival sequence is finite, and the functions in Equation (6) have bounded derivatives through order $2M$ on every $K\in {\mathcal{P}}_i$.

2. 

The exact and SAP trajectories remain in a bounded neighborhood with the same finite event itinerary. In particular, they have the same number and ordering of threshold crossings, resets, and refractory releases, and no such event coincides with a scan or batch boundary.

3. 

The endpoint screen is complete on this itinerary. On every active scan bin, the SAP no-reset predictor either remains below threshold throughout the bin, or it has exactly one threshold zero and is suprathreshold at the right endpoint. Every exact crossing has one corresponding predictor zero in a common initial bracket ${\mathcal{B}}_{i,\ell }\subseteq K_{i,\ell }$ of width ${\kappa }_{i,\ell }\le \left|K_{i,\ell }\right|\le h$.

4. 

Let ${\varphi }_{i,\ell }$ be the smooth unreset continuation of the exact pre-spike voltage on ${\mathcal{B}}_{i,\ell }$. Its threshold crossing ${\tau }_{i,\ell }$ is transversal, and

$${\displaystyle \frac{\mathrm{d}{\varphi }_{i,\ell }}{\mathrm{d}t}}\left(t\right)\ge {\gamma }_{i,\ell }>0,t\in {\mathcal{B}}_{i,\ell }.$$

The corresponding SAP predictor zero also lies in ${\mathcal{B}}_{i,\ell }$.

Let ${\mathcal{E}}_i$ index the crossings in the batch, and define the bisection remainder

$${\epsilon }_{i,\ell }^{\mathrm{bis}}=max\left\{{\epsilon }_t,{\kappa }_{i,\ell }{2}^{-B}\right\}.$$

(8)

Then there exist constants ${\Lambda }_i^{\mathrm{q}}\ge 0$ and ${\Lambda }_i^{\mathrm{t}}\ge 0$ such that

$$\left|V_i\left(T_1^{-}\right)-{\widehat{V}}_i\left(T_1^{-}\right)\right|\le {\Lambda }_i^{\mathrm{q}}\sum _{K\in {\mathcal{P}}_i}{d}_{i,K}+{\Lambda }_i^{\mathrm{t}}\sum _{\ell \in {\mathcal{E}}_i}{\epsilon }_{i,\ell }^{\mathrm{bis}}.$$

(9)

The constants are uniform on the stated event-itinerary neighborhood. They depend on the model parameters, H, the transversality bounds, and local flow sensitivities, but not on h, ${\epsilon }_t$, or B. Here ${\Lambda }_i^{\mathrm{q}}$ is dimensionless, whereas ${\Lambda }_i^{\mathrm{t}}$ has units of voltage per unit time.

If ${\mathcal{G}}_{i,K}^{\left(2M\right)}\le {\overline{\mathcal{G}}}_i^{\left(2M\right)}$ uniformly and ${max}_{\ell }{\epsilon }_{i,\ell }^{\mathrm{bis}}\le c_i{h}^{2M}$ along an itinerary-preserving refinement sequence, then

$$\left|V_i\left(T_1^{-}\right)-{\widehat{V}}_i\left(T_1^{-}\right)\right|\le \left[{\Lambda }_i^{\mathrm{q}}{\beta }_{M}H{\overline{\mathcal{G}}}_i^{\left(2M\right)}+{\Lambda }_i^{\mathrm{t}}{c}_i\#{\mathcal{E}}_i\right]h^{2M}.$$

(10)

Thus the one-batch voltage error is $O\left(h^{2M}\right)$ under the stated refinement conditions.

Proof. 

We first estimate the defect of one smooth predictor. Fix $K=[a_K,b_K]\in {\mathcal{P}}_i$ and $x\in (a_K,b_K]$. For an initial voltage v, the exact unreset flow is

$${\Phi }_{i,K,x}\left(v\right)=F_i(x,a_K)v+{\int }_{a_K}^x{\displaystyle \frac{q_i\left(s\right)}{C_i}}{F}_i(x,s)\mathrm{d}s.$$

(11)

The SAP map replaces only the integral by the M-node Gauss–Legendre rule,

$${\widehat{\Phi }}_{i,K,x}\left(v\right)=F_i(x,a_K)v+{\displaystyle \frac{x-a_K}{2}}{\displaystyle \sum _{m=1}^M}{\omega }_m{\displaystyle \frac{q_i\left(s_{m,x}\right)}{C_i}}{F}_i(x,s_{m,x}),$$

(12)

where

$$s_{m,x}={\displaystyle \frac{a_K+x}{2}}+{\displaystyle \frac{x-a_K}{2}}{\xi }_m.$$

The Gauss–Legendre remainder formula and $x-a_K\le \left|K\right|$ give

$$\left|{\Phi }_{i,K,x}\left(v\right)-{\widehat{\Phi }}_{i,K,x}\left(v\right)\right|\le d_{i,K}.$$

(13)

This estimate is uniform in x, so it applies both to accepted endpoint updates and to predictor evaluations made during bisection.

If the exact and numerical segment maps start from different voltages v and $\widehat{v}$, then

$$\left|{\Phi }_{i,K,x}\left(v\right)-{\widehat{\Phi }}_{i,K,x}\left(\widehat{v}\right)\right|\le |F_i(x,a_K)\left|\right|v-\widehat{v}|+d_{i,K}.$$

(14)

On the bounded neighborhood in the theorem, $|F_i(x,a_K)|\le exp(L_i\left|K\right|)$ for some finite $L_i$. Products of these smooth propagation factors over the batch are therefore bounded by $exp\left(L_{i}H\right)$, independently of the number of scan intervals.

We next examine a threshold crossing. Let ${\tau }_{i,\ell }$ be the zero of the exact unreset continuation ${\varphi }_{i,\ell }$, let ${\overline{\tau }}_{i,\ell }$ be the corresponding zero of the SAP predictor before bisection truncation, and let ${\widehat{\tau }}_{i,\ell }$ be the returned bisection time. Denote the uniform discrepancy between the exact unreset continuation and the numerical predictor on their common bracket by $r_{i,\ell }$. Because ${\varphi }_{i,\ell }\left({\tau }_{i,\ell }\right)=V_{\mathrm{th},i}$ and the SAP predictor equals $V_{\mathrm{th},i}$ at ${\overline{\tau }}_{i,\ell }$, transversality and the mean-value theorem yield

$$\begin{array}{ccc}\hfill {\gamma }_{i,\ell }|{\overline{\tau }}_{i,\ell }-{\tau }_{i,\ell }|& \le \hfill & |{\varphi }_{i,\ell }\left({\overline{\tau }}_{i,\ell }\right)-{\varphi }_{i,\ell }\left({\tau }_{i,\ell }\right)|\hfill \\ \hfill & =\hfill & |{\varphi }_{i,\ell }\left({\overline{\tau }}_{i,\ell }\right)-V_{\mathrm{th},i}|\le r_{i,\ell }.\hfill \end{array}$$

(15)

This argument uses the smooth unreset continuation, not the physical voltage after its reset.

The retained bisection bracket contains ${\overline{\tau }}_{i,\ell }$. Its width is at most ${\epsilon }_t$ when the requested tolerance is reached and at most ${\kappa }_{i,\ell }{2}^{-B}$ when the iteration cap is reached. Hence

$$|{\widehat{\tau }}_{i,\ell }-{\overline{\tau }}_{i,\ell }|\le {\epsilon }_{i,\ell }^{\mathrm{bis}},$$

(16)

and Equations (15) and (16) give

$$|{\widehat{\tau }}_{i,\ell }-{\tau }_{i,\ell }|\le {\displaystyle \frac{r_{i,\ell }}{{\gamma }_{i,\ell }}}+{\epsilon }_{i,\ell }^{\mathrm{bis}}.$$

(17)

It remains to propagate these defects through the batch with constants that do not grow with the number of scan intervals. Order the physical events in the common itinerary and regard the intervening members of ${\mathcal{P}}_i$ as smooth stages. At a regular endpoint of a smooth stage, let $e_r$ be the voltage discrepancy and let ${\theta }_r$ be the magnitude of the crossing- or refractory-release-time discrepancy that is currently carried by the two trajectories; set ${\theta }_r=0$ when no such offset is active. Bounded derivatives of the flow on the itinerary-preserving neighborhood give uniform constants ${\overline{L}}_i,c_i,D_i\ge 0$ such that, for a smooth stage $K_r$ of length ${\Delta }_r=\left|K_r\right|$,

$$\begin{array}{c}{e}_{r+1}\le exp\left({\overline{L}}_i{\Delta }_r\right)e_r+c_i{\Delta }_r{\theta }_r+D_i{d}_{i,K_r},\end{array}$$

(18)

$$\begin{array}{c}{\theta }_{r+1}={\theta }_r.\end{array}$$

(19)

Indeed, the first coefficient is the usual smooth-flow stability factor, so $A_r\le exp({\overline{L}}_i|K_r\left|\right)$. The dependence on a displaced event endpoint enters through the vector field integrated over $K_r$; hence its coefficient satisfies $B_r\le c_i|K_r|=O(|K_r\left|\right)$. The final term follows from Equation (14). These quantitative scalings are essential because the number of smooth stages may increase as h decreases.

Choose a fixed conversion factor ${\nu }_i>0$ with units of voltage per unit time, define $z_r=e_r+{\nu }_i{\theta }_r$, and choose ${\tilde{L}}_i\ge max\{{\overline{L}}_i,c_i/{\nu }_i\}$. Since $e^x\ge 1+x$ for $x\ge 0$, Equations (18) and (19) imply

$$z_{r+1}\le exp\left({\tilde{L}}_i{\Delta }_r\right)z_r+D_i{d}_{i,K_r}.$$

(20)

Consequently, discrete Gronwall, or direct telescoping of Equation (20), over any consecutive smooth block $r=p,\dots ,q-1$ gives

$$z_q\le exp\left({\tilde{L}}_i{\displaystyle \sum _{r=p}^{q-1}}{\Delta }_r\right)\left(z_p+D_i{\displaystyle \sum _{r=p}^{q-1}}{d}_{i,K_r}\right).$$

(21)

The sum of all smooth-stage lengths in the batch is at most H. Thus their combined amplification is at most $exp\left({\tilde{L}}_{i}H\right)$, independently of h and of the number of scan intervals.

We now separate the physical event maps from these smooth stages. The fixed event itinerary contains a finite number $J_i$ of visible-arrival jumps, crossing/reset pairs, and refractory releases; $J_i$ is independent of the scan partition. For a crossing event, uniform smooth-stage estimates on its common bracket and Equation (13) bound the predictor discrepancy by

$$r_{i,\ell }\le c_{i,\ell }^{\mathrm{pred}}{z}^{-}+D_{i,\ell }^{\mathrm{pred}}{d}_{i,K_{i,\ell }}.$$

Equation (17) then controls the new timing error. Both trajectories receive the same reset voltage, and the timing error is carried to the corresponding refractory release. At that release, boundedness of the vector field gives Lipschitz dependence of the restarted voltage on the release time. Visible-arrival jumps occur at the same prescribed times and are also locally Lipschitz. It follows that the jth physical event map satisfies

$$z_j^{+}\le Q_{i,j}{z}_j^{-}+R_{i,j}{d}_{i,K\left(j\right)}+{\mathbf{1}}_{\left\{j\mathrm{is}\mathrm{a}\mathrm{crossing}\right\}}{S}_{i,j}{\epsilon }_{i,\ell \left(j\right)}^{\mathrm{bis}},$$

(22)

where the defect term is taken as zero for an event with no associated predictor interval. The constants $Q_{i,j}$ and $R_{i,j}$ are dimensionless, whereas $S_{i,j}$ has units of voltage per unit time. They are uniform on the stated neighborhood: in particular, the factors ${\gamma }_{i,\ell }^{-1}$ from Equation (17) are bounded there.

Finally, telescope the alternating smooth-block estimates (21) and the $J_i$ event-map estimates (22). With

$$P_i=exp\left({\tilde{L}}_{i}H\right){\displaystyle \prod _{j=1}^{J_i}}max\{1,Q_{i,j}\},R_i^{*}=\underset{1\le j\le J_i}{max}{R}_{i,j},S_i^{*}=\underset{1\le j\le J_i}{max}{S}_{i,j},$$

where an empty maximum is zero, the common initial data $z_0=e_0+{\nu }_i{\theta }_0=0$ yield

$$e_{\mathrm{final}}\le z_{\mathrm{final}}\le P_i(D_i+R_i^{*})\sum _{K\in {\mathcal{P}}_i}{d}_{i,K}+P_i{S}_i^{*}\sum _{\ell \in {\mathcal{E}}_i}{\epsilon }_{i,\ell }^{\mathrm{bis}}.$$

(23)

The fixed finite event count, the total-length bound, and the uniform local coefficients show explicitly that neither multiplier depends on h, ${\epsilon }_t$, or B. Taking ${\Lambda }_i^{\mathrm{q}}=P_i(D_i+R_i^{*})$ and ${\Lambda }_i^{\mathrm{t}}=P_i{S}_i^{*}$ proves Equation (9).

Finally,

$$\sum _{K\in {\mathcal{P}}_i}{\left|K\right|}^{2M+1}\le h^{2M}\sum _{K\in {\mathcal{P}}_i}\left|K\right|\le H{h}^{2M}.$$

Combining this inequality with the uniform derivative bound and the assumed bisection scaling proves Equation (10).    □

Equation (9) is a conditional one-batch voltage estimate. It applies to any batch whose initial data and event itinerary satisfy the stated conditions. The theorem provides an error and convergence analysis of the SAP voltage update; the experiments below separately evaluate the tested accuracy-efficiency behavior of the recurrent implementation.

5. Results and Discussion

5.1. Million-Neuron Network Experiments

We evaluated Euler, SAP, and NEST [10] across four constant-drive regimes. The input current densities $I\in \{0.85,1,2,3\}\mu \mathrm{A}{\mathrm{cm}}^{-2}$ produced same-input NEST reference firing rates of $9.16$, $14.69$, $57.97$, and $88.75\mathrm{Hz}$, respectively. This range allowed the solver comparison to include both weakly and strongly driven network activity. All main comparisons retained the same three-solver suite; the subsequent parameter sweeps varied SAP only. The complete configuration-level records are provided in Supplementary File S1.

5.1.1. Experimental Design and Evaluation Metrics

The topology, neuronal and synaptic parameters, external drive, initialization, and output-sampling grid were fixed across solver configurations. The experimental network contained $1,000,000$ neurons, comprising $800,000$ excitatory and $200,000$ inhibitory neurons. Each neuron had fixed average in-degree $d=100$. Each experiment covered $T=1000\mathrm{ms}$ with a batch-delivery interval of $H=1\mathrm{ms}$ under one of four constant drives, $I\in \{0.85,1,2,3\}\mu \mathrm{A}{\mathrm{cm}}^{-2}$.

The same fixed topology artifact was reused by every solver. Each target neuron received 100 distinct presynaptic neurons sampled without replacement, with self-connections excluded; the graph therefore contained neither autapses nor repeated source–target pairs. The first $800,000$ source indices were excitatory and entered the AMPA channel, while the remaining $\mathrm{200,000}$ were inhibitory and entered the ${\mathrm{GABA}}_A$ channel. Each active connection used a fixed dimensionless multiplier $w_{ij}\sim \mathcal{U}[0,1)$ drawn during topology construction.

The area-normalized membrane parameters were $C_i=0.75\mu \mathrm{F}{\mathrm{cm}}^{-2}$, $g_{L,i}=1/30000\mathrm{S}{\mathrm{cm}}^{-2}$, and $V_L=-75\mathrm{mV}$. All neurons started at $V_i\left(0\right)=-57.5\mathrm{mV}$ with zero receptor traces and were not initially refractory. The threshold and reset potentials were $-50$ and $-65\mathrm{mV}$, respectively, and $T_{\mathrm{ref}}=5\mathrm{ms}$. AMPA and ${\mathrm{GABA}}_A$ had reversal potentials 0 and $-70\mathrm{mV}$ and decay times 2 and $10\mathrm{ms}$, respectively. Their base conductance densities were $g^{\mathrm{AMPA}}=1/55\mathrm{mS}{\mathrm{cm}}^{-2}$ and $g^{{\mathrm{GABA}}_A}=0.1\mathrm{mS}{\mathrm{cm}}^{-2}$, followed by the common scale $0.1$. Thus, a connection produced conductance-density increments $0.1g^{\mathrm{AMPA}}{w}_{ij}\mathrm{mS}{\mathrm{cm}}^{-2}$ and $0.1g^{{\mathrm{GABA}}_A}{w}_{ij}\mathrm{mS}{\mathrm{cm}}^{-2}$ in the two channels.

The NEST parameter conversion used an effective membrane area of ${10}^{-3}{\mathrm{cm}}^2$: $C_m=750\mathrm{pF}$, $g_L=100/3\mathrm{nS}$, and $I_e=1000I\mathrm{pA}$. Excitatory and inhibitory NEST connection weights were $+100g^{\mathrm{AMPA}}{w}_{ij}$ and $-100g^{{\mathrm{GABA}}_A}{w}_{ij}\mathrm{nS}$, respectively. The custom solvers stored voltage, receptor traces, quadrature data, and event times in double precision, while the sparse topology and weight arrays used single precision. Sampled voltage arrays were recorded in single precision.

The numerical settings varied only by solver and by the planned SAP ablations. Euler and SAP were implemented as graphics processing unit (GPU) solvers and were evaluated at $h\in \{1,0.1,0.01,0.001\}\mathrm{ms}$. The main SAP configuration used $M=2$, ${\epsilon }_t=0.01\mathrm{ms}$, and at most 15 bisection iterations. The ablations varied $M\in \{1,2,3,4\}$ and ${\epsilon }_t\in \{0.1,0.01,0.001\}\mathrm{ms}$. Every tested scan width satisfied $h<T_{\mathrm{ref}}=5\mathrm{ms}$. NEST 3.10.0-rc.2 used the iaf_cond_exp model with a $1\mathrm{ms}$ delay. For each drive, the NEST simulation at $h=0.001\mathrm{ms}$ served as the numerical reference.

All configurations were evaluated on the same Dell Precision 7960 workstation, equipped with an Intel Xeon w5-3525 processor, 125 GiB of system memory, and three NVIDIA GeForce RTX 4090 GPUs with 24 GiB of memory each. Each Euler or SAP configuration used one RTX 4090 and PyTorch 2.11.0 with CUDA 12.6. NEST 3.10.0-rc.2 ran on the same workstation with eight CPU threads.

The interpretation of h depends on the solver. For SAP, h is the common scan width, although propagation may begin later within a bin when a refractory period ends. For NEST, h is the simulation resolution that bounds each interval advanced by the adaptive Runge–Kutta–Fehlberg $4\left(5\right)$ integrator in the GNU Scientific Library.

We quantified voltage accuracy by the sampled-neuron root-mean-square error (RMSE). We sampled $n_s=1000$ neurons at intervals of $\Delta s=1\mathrm{ms}$. Let $\mathcal{C}$ denote the resulting sampling grid and ${\mathcal{P}}_s$ the sampled-neuron set. For numerical voltage $V_{j,m}^{\mathrm{alg}}$ and reference voltage $V_{j,m}^{★}$, the metric is

$${\mathrm{RMSE}}_{V,\mathrm{sample}}={\left[{\displaystyle \frac{1}{\left|\mathcal{C}\right|}}\sum _{m\in \mathcal{C}}{\displaystyle \frac{1}{n_s}}\sum _{j\in {\mathcal{P}}_s}{\left(V_{j,m}^{\mathrm{alg}}-V_{j,m}^{★}\right)}^2\right]}^{1/2}.$$

(24)

We quantified computational cost by neuronal ordinary differential equation (ODE)-solving time. For every solver, this metric retains only work used to advance the neuronal differential equations. The SAP value consists of two parts: the pure ODE-solving time in the main propagation stage and the pure ODE-solving time of the predictor evaluations performed during bisection for spike localization. Data movement, input/output, recording, post-processing, synaptic delivery, and NEST voltage polling are excluded. The metric therefore compares neuronal ODE work rather than end-to-end wall-clock performance.

5.1.2. Sampled-Neuron Voltage Dynamics

Before comparing aggregate error and cost, we examined the membrane-voltage trajectory of one sampled neuron. Figure 1 compares the high-activity $I=3\mu \mathrm{A}{\mathrm{cm}}^{-2}$ regime across four Euler/SAP step sizes over the first $100\mathrm{ms}$. The same $1\mathrm{ms}$ recording grid and the fixed NEST $h=0.001\mathrm{ms}$ reference were used in every row.

At $h=1\mathrm{ms}$, the fixed trace shows several early reset events for which SAP and NEST are recorded at nearby sampling points, whereas Euler records the corresponding events later. The remaining rows show how the recorded Euler and SAP trajectories change as their step sizes are reduced. The repeated threshold approaches, abrupt resets, and $5\mathrm{ms}$ plateaus provide a qualitative illustration of the hybrid LIF dynamics for one prespecified neuron. The aggregate comparison below uses all 1000 sampled neurons over the complete simulation.

5.1.3. Accuracy-Efficiency Pareto Frontiers Across Activity Regimes

The central empirical comparison concerns sampled-neuron voltage RMSE and pure neuronal ODE-solving time over the tested configurations. Figure 2 plots every evaluated operating point under these two metrics. Their lower-left envelope is the tested Pareto frontier: moving toward it reduces error, cost, or both, and a point directly dominates a comparator when it improves both quantities.

Within the fixed-seed experiments, the clearest joint reductions occurred in the two high-activity regimes shown in panels (c) and (d). At $h=1\mathrm{ms}$, SAP had both lower RMSE and lower pure ODE-solving time than the corresponding NEST $h=1\mathrm{ms}$ configurations. The RMSE reductions were $0.85\mathrm{mV}$ and $1.67\mathrm{mV}$, while the corresponding NEST-to-SAP ODE-time ratios were $3.46$ and $4.75$. These points are direct lower-left improvements over those same-resolution NEST configurations under the reported metrics.

The same tested ordering appeared in selected cross-resolution comparisons with Euler. SAP at $h=1\mathrm{ms}$ had lower RMSE than Euler at $h=0.01\mathrm{ms}$ in both high-activity regimes, despite using a 100-fold larger scan width. The RMSE reductions were $0.23\mathrm{mV}$ and $0.64\mathrm{mV}$, and SAP required less than half the pure ODE-solving time. When SAP at $h=0.1\mathrm{ms}$ was compared with Euler at $h=0.001\mathrm{ms}$, it reduced RMSE by $0.21\mathrm{mV}$ and $0.56\mathrm{mV}$ while reducing pure ODE-solving time by factors exceeding four.

Selected approximately cost-matched points showed the same ordering in the two high-activity runs. Relative to Euler at $h=0.001\mathrm{ms}$, SAP at $h=0.01\mathrm{ms}$ reduced RMSE from $6.88$ to $6.64\mathrm{mV}$ at $I=2\mu \mathrm{A}{\mathrm{cm}}^{-2}$ and from $6.51$ to $5.90\mathrm{mV}$ at $I=3\mu \mathrm{A}{\mathrm{cm}}^{-2}$. Thus, the observed frontier ordering was not restricted to points with markedly different pure ODE-solving costs.

The lower-activity regimes bound this result. SAP at $h=1\mathrm{ms}$ was used approximately $1/2.20$ of the pure ODE-solving time of Euler at $h=0.01\mathrm{ms}$, but their RMSE values differed by only $0.019\mathrm{mV}$ in either direction. The ordering against NEST also changed as h decreased. At $h=0.1\mathrm{ms}$, SAP had lower RMSE only in the highest-activity regime; at $h=0.01\mathrm{ms}$, NEST had lower RMSE in all four regimes while SAP remained less expensive. Figure 2 therefore supports an activity- and configuration-dependent frontier improvement rather than uniform solver dominance.

5.1.4. SAP Hyperparameter Ablations

We next examined which SAP hyperparameters produced useful operating points using RMSE and pure ODE-solving time as the respective accuracy and efficiency metrics. The sweep separated the effects of quadrature order, spike-localization tolerance, and scan width. Table 1 and Table 2 report the $h=1\mathrm{ms}$ slices because they isolate the two local parameters at the coarsest tested scan width. The supplementary records contain all four scan widths. Figure 3 complements the absolute values in the tables by showing the accuracy change and ODE-time multiplier relative to the coarsest setting in each sweep.

Increasing M from 1 to 2 reduced RMSE in all four regimes. This change increased pure ODE-solving time by approximately $32\%$, yielding an accuracy-efficiency tradeoff rather than free improvement. Beyond $M=2$, the errors changed by at most $0.022\mathrm{mV}$ and were not consistently smaller. In particular, $M=4$ required approximately 49–$50\%$ more pure ODE-solving time than $M=2$. Its RMSE was equal or larger in three regimes, while at $I=0.85\mu \mathrm{A}{\mathrm{cm}}^{-2}$ it improved RMSE by only $0.016\mathrm{mV}$. Thus, $M=2$ remained the practical compromise: it captured most of the order-dependent gain without moving to the high-cost, low-return end of the SAP frontier, as summarized in Figure 3a,d.

Tightening ${\epsilon }_t$ from $0.1$ to $0.01\mathrm{ms}$ reduced RMSE by $0.232\mathrm{mV}$ at $I=2\mu \mathrm{A}{\mathrm{cm}}^{-2}$ and by $0.742\mathrm{mV}$ at $I=3\mu \mathrm{A}{\mathrm{cm}}^{-2}$. The corresponding pure ODE-solving times increased by approximately $59\%$. The same tolerance change altered RMSE by only $-0.014$ and $+0.002\mathrm{mV}$ in the two lower-drive cases. Further tightening to $0.001\mathrm{ms}$ increased cost by another 37–$38\%$ without a consistent RMSE reduction. Hence, the tolerance sweep places ${\epsilon }_t=0.01\mathrm{ms}$ near the knee of the tested SAP frontier for the two high-activity regimes; the tighter setting buys little additional accuracy (Figure 3b,e).

Across the tested scan widths, the accuracy-efficiency ordering depended on the activity regime. For example, at $I=3\mu \mathrm{A}{\mathrm{cm}}^{-2}$ with $M=2$ and ${\epsilon }_t=0.01\mathrm{ms}$, RMSE was $5.917$, $5.950$, $5.895$, and $5.967\mathrm{mV}$ as h decreased from 1 to $0.001\mathrm{ms}$. The associated pure ODE-solving time increased from $33.72$ to $7365\mathrm{s}$. At $I=0.85\mu \mathrm{A}{\mathrm{cm}}^{-2}$, moving across the same settings changed RMSE by only $0.025\mathrm{mV}$ while increasing pure ODE-solving time by more than 200-fold. Thus, the lower-left operating points in this sweep were concentrated at intermediate or coarse settings (Figure 3c,f).

5.1.5. Discussion

Within the tested sweep, selected coarse SAP configurations produced favorable sampled-voltage accuracy-efficiency points, particularly in the two higher-activity regimes. The quadrature-order and localization-tolerance ablations showed that intermediate settings captured the useful operating points without incurring the largest ODE-solving costs. These observations motivate the use of $M=2$ and ${\epsilon }_t=0.01\mathrm{ms}$ in the main SAP comparisons.

The NEST run at $h=0.001\mathrm{ms}$ serves consistently as the numerical reference. All reported rankings concern sampled-neuron voltage RMSE and pure neuronal ODE-solving time for the fixed-topology million-neuron network, four constant-current densities, and explicitly tested configurations. Within this scope, intermediate SAP settings occupied favorable points of the tested accuracy-efficiency frontier.

6. Conclusions

We introduced SAP for conductance-based LIF equations by combining analytic homogeneous propagation, Gauss–Legendre quadrature, and spike localization. For a prescribed visible-arrival schedule and fixed bracketed transversal spiking event itinerary, the accompanying analysis bounds the one-batch voltage error and gives conditional convergence under consistent refinement. In the million-neuron network experiments, selected SAP configurations in the two higher-activity regimes yielded lower sampled-voltage RMSE at comparable pure ODE-solving time, or lower pure ODE-solving time at comparable RMSE, than the tested Euler and NEST configurations. The ablations identified $M=2$ and ${\epsilon }_t\approx 0.01\mathrm{ms}$ as practical intermediate settings for the higher-activity cases. The results show that SAP can improve the accuracy-efficiency Pareto frontier of conductance-based LIF simulation, highlighting its practical value in large-scale spiking neural network applications.