Differentiable Spike-Time Solvers

• Differentiable spike-time solvers are methods that compute gradients via implicit differentiation of threshold equations, enabling effective training of spiking neural networks.
• They employ various formulations—timing-based, event-driven, and representation-based—to retain temporal coding while addressing the discontinuity of spike generation.
• These methods enhance SNN performance across neuron models by balancing exact gradient computation with computational efficiency and model generality.

Differentiable spike-time solvers are methods for training spiking neural networks by propagating gradients through the quantities that define spiking dynamics: threshold-crossing times, reset maps, event sequences, or continuous representations of spike trains. In the exact event-driven setting, spike times are treated as implicit functions of parameters and differentiated through threshold equations, causal dependencies, or numerical root-finding; in related sub-differentiable formulations, entire spike trains are mapped to continuous layer representations so that standard backpropagation can proceed without per-time-step surrogate derivatives (Lee et al., 2022, Yamamoto et al., 2022, König et al., 9 Mar 2026, Meng et al., 2022, Holberg et al., 2024). Across these lines of work, the central objective is to retain temporal coding while avoiding the usual incompatibility between hard thresholding and gradient-based optimization.

1. Conceptual scope and historical lines

The modern literature contains three closely related formulations of differentiability in spiking systems. One line derives exact gradients directly from spike times and threshold conditions. Another line differentiates numerical solvers that integrate continuous dynamics and locate discrete events. A third line avoids direct spike-time differentiation by encoding a spike train into a continuous spike representation and proving that the resulting layer map is sub-differentiable (Banerjee, 2014, Lee et al., 2022, Meng et al., 2022, König et al., 9 Mar 2026).

Formulation	Differentiated object	Representative papers
Timing-based exact methods	Implicit spike times and reset paths	(Banerjee, 2014, Lee et al., 2022, Yamamoto et al., 2022)
Representation-based methods	Layerwise spike representation and clamp map	(Meng et al., 2022)
Event-driven autodiff solvers	ODE/SDE trajectories, root-finding, resets	(König et al., 9 Mar 2026, Holberg et al., 2024)

A recurrent misconception in the field has been that exact gradients for SNNs with respect to weights do not exist because spike generation is discontinuous. The forward-propagation framework explicitly challenges this view: by applying the implicit function theorem to discrete spike times, it proves that SNNs have well-defined gradients with respect to their weights, even though they are non-differentiable in time (Lee et al., 2022). Earlier timing-based work had already shown that precise spike-train-to-spike-train learning can be handled in multilayer feedforward networks by differentiating output spike times and propagating those sensitivities through intermediate spikes (Banerjee, 2014).

This suggests that “differentiability” in spiking models should not be identified solely with smooth membrane trajectories. In the literature surveyed here, differentiability is often piecewise, implicit, event-conditioned, or numerical-solver-relative rather than globally smooth in the elementary ANN sense.

2. Deterministic event formulations

A canonical starting point is the leaky integrate-and-fire formulation used in timing-based backpropagation without single-spike restrictions. For neuron $i$,

$\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$

with postsynaptic current

$$I_i(t) = \sum_j w_{ji}\sum_k \kappa(t-t_j^k), \qquad \kappa(t)=\beta_I e^{-t/\tau_I}\theta(t).$$

A spike is emitted at the smallest $t_i^k$ such that $V_i(t_i^k)=V_{th}$, after which $V_i(t_i^{k+})\leftarrow 0$, while $I_i(t)$ does not reset, so past spikes continue to contribute to future current (Yamamoto et al., 2022).

The forward-propagation framework recasts continuous-time LIF neurons in a post-synaptic-kernel form. If $\mathcal F_i=\{t_i^1,t_i^2,\dots\}$ denotes the sorted output spike times of neuron $i$, then

$$y_{ji}(t)=\sum_{f\in\mathcal F_j} h_{ji}(t-f), \qquad h_{ji}(t)=\frac{e^{-\alpha_j t}-e^{-\beta_i t}}{\beta_i-\alpha_j}u(t),$$

and the membrane potential becomes

$\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$0

with each output spike time satisfying the threshold condition $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$1 (Lee et al., 2022).

An earlier multilayer feedforward formulation based on the Spike-Response Model employs a similar event-centered perspective but assigns virtual weights to spikes rather than synapses. Its membrane potential is written as

$\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$2

where $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$3 is the PSP kernel and $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$4 is the afterhyperpolarization kernel. This device permits a perturbation analysis of individual spike times and weights, with online-style updates derived from first-order Taylor expansions (Banerjee, 2014).

Across these deterministic formulations, spike times are not auxiliary outputs; they are the primary state variables through which synaptic interactions, resets, and learning signals are routed.

3. Exact gradient mechanisms

The basic exact-gradient identity in threshold-based SNNs follows from differentiating the firing equation $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$5 at $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$6. For a weight $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$7,

$\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$8

For $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$9, the numerator is $$I_i(t) = \sum_j w_{ji}\sum_k \kappa(t-t_j^k), \qquad \kappa(t)=\beta_I e^{-t/\tau_I}\theta(t).$$0, and the denominator is the time derivative of the full membrane potential, including reset contributions from previous self-spikes. The resulting derivative has a closed form in terms of smooth exponential kernels (Lee et al., 2022).

Banerjee’s multilayer spike-train learning framework begins from an error functional that compares desired and observed spike trains by their putative impact on a virtual postsynaptic neuron. The kernel family is

$$I_i(t) = \sum_j w_{ji}\sum_k \kappa(t-t_j^k), \qquad \kappa(t)=\beta_I e^{-t/\tau_I}\theta(t).$$1

and integrating out the nuisance parameters $$I_i(t) = \sum_j w_{ji}\sum_k \kappa(t-t_j^k), \qquad \kappa(t)=\beta_I e^{-t/\tau_I}\theta(t).$$2 yields a closed-form spike-train discrepancy $$I_i(t) = \sum_j w_{ji}\sum_k \kappa(t-t_j^k), \qquad \kappa(t)=\beta_I e^{-t/\tau_I}\theta(t).$$3 with $$I_i(t) = \sum_j w_{ji}\sum_k \kappa(t-t_j^k), \qquad \kappa(t)=\beta_I e^{-t/\tau_I}\theta(t).$$4 and $$I_i(t) = \sum_j w_{ji}\sum_k \kappa(t-t_j^k), \qquad \kappa(t)=\beta_I e^{-t/\tau_I}\theta(t).$$5 iff the two spike-time vectors coincide. Gradients are then decomposed into $$I_i(t) = \sum_j w_{ji}\sum_k \kappa(t-t_j^k), \qquad \kappa(t)=\beta_I e^{-t/\tau_I}\theta(t).$$6, $$I_i(t) = \sum_j w_{ji}\sum_k \kappa(t-t_j^k), \qquad \kappa(t)=\beta_I e^{-t/\tau_I}\theta(t).$$7, and $$I_i(t) = \sum_j w_{ji}\sum_k \kappa(t-t_j^k), \qquad \kappa(t)=\beta_I e^{-t/\tau_I}\theta(t).$$8, where sums over $$I_i(t) = \sum_j w_{ji}\sum_k \kappa(t-t_j^k), \qquad \kappa(t)=\beta_I e^{-t/\tau_I}\theta(t).$$9 capture indirect through-spike effects (Banerjee, 2014).

Yamamoto et al. derive timing-based backpropagation by changing variables from spike time to

$t_i^k$0

Within each inter-spike interval, the membrane potential admits the closed form

$t_i^k$1

The next spike time $t_i^k$2 is the smallest root of $t_i^k$3, and for $t_i^k$4,

$t_i^k$5

Backpropagation proceeds through downstream spike times and output quantities by the chain rule, with

$t_i^k$6

Weight gradients are then expressed as $t_i^k$7, where $t_i^k$8 is obtained by differentiating the spike-time root formula (Yamamoto et al., 2022).

These derivations share a structural property: the discontinuity of the spike itself is not differentiated directly as a Heaviside surrogate. Instead, the threshold equation defining the event is differentiated, and reset-induced dependencies are handled explicitly.

4. Causality, multiple spikes, and memory structure

A major limitation of many conventional timing-based methods is the single-spike restriction, often associated with time-to-first-spike coding. Yamamoto et al. remove this restriction: each neuron $t_i^k$9 may fire $V_i(t_i^k)=V_{th}$0 times, indexed by $V_i(t_i^k)=V_{th}$1. Every term in the spike-time coefficients $V_i(t_i^k)=V_{th}$2 and $V_i(t_i^k)=V_{th}$3, and in their derivatives, sums over all previous presynaptic spikes and all previous self-spikes $V_i(t_i^k)=V_{th}$4. Gradient backpropagation likewise sums reset paths from each spike to all later spikes of the same neuron, and synaptic paths into many postsynaptic spikes (Yamamoto et al., 2022).

The forward-propagation algorithm organizes these dependencies through a time-sorted event list and a causality graph. For each new spike $V_i(t_i^k)=V_{th}$5, it records the set of input spikes that actually contributed and the last reset time of the neuron. The local linear system

$V_i(t_i^k)=V_{th}$6

is solved by back-substitution for the new row because $V_i(t_i^k)=V_{th}$7 is strictly lower-triangular plus diagonal. This yields forward-time exact gradient accumulation without full matrix inversion (Lee et al., 2022).

These exact timing-based methods incur characteristic path-combinatorial costs. Naïve timing-based backprop explores $V_i(t_i^k)=V_{th}$8 reset paths. In Yamamoto et al., a two-stage loop prunes neurons with no further spikes after $V_i(t_i^k)=V_{th}$9 steps to reduce memory; with two-stage partitioning, memory drops from $V_i(t_i^{k+})\leftarrow 0$0 to roughly $V_i(t_i^{k+})\leftarrow 0$1, and the implementation is GPU-parallelizable over $V_i(t_i^{k+})\leftarrow 0$2 loops (Yamamoto et al., 2022).

The empirical interpretation offered in the same work is that allowing multiple spikes naturally improves the computational capacity of SNNs. A plausible implication is that differentiable spike-time solvers become more expressive when they model not only the first threshold crossing but also the entire reset-mediated event cascade of each neuron.

5. Differentiable numerical solvers for ODE and SDE neuron models

Closed-form spike-time expressions are available only for restricted neuron classes and parameter regimes. Eventax addresses this limitation by combining differentiable numerical ODE solvers with event-based spike handling in JAX. In this framework, each neuron is treated as a hybrid dynamical system whose continuous state $V_i(t_i^{k+})\leftarrow 0$3 evolves according to

$V_i(t_i^{k+})\leftarrow 0$4

while spikes are defined by an event function $V_i(t_i^{k+})\leftarrow 0$5. Diffrax monitors $V_i(t_i^{k+})\leftarrow 0$6 at each step, and if it switches from negative to positive, a built-in root-finder using bisection plus secant locates $V_i(t_i^{k+})\leftarrow 0$7 to any desired tolerance, such as $V_i(t_i^{k+})\leftarrow 0$8. At the exact event time, a model-specific reset mapping $V_i(t_i^{k+})\leftarrow 0$9 is applied, and reverse-mode autodiff propagates through the ODE integration, root-finding, and reset callbacks exactly, without surrogate approximations (König et al., 9 Mar 2026).

Eventax exposes this mechanism through a NeuronModel interface in Equinox with methods such as init_state(n), dynamics(t, y), spike_condition(t, y), input_spike(y, w), and reset_spiked(y, mask). Feed-forward and recurrent networks are assembled through classes such as FFEvNN and EvRNN, while .ttfs(in_spikes) and .state_at_t(in_spikes, T) provide task-level interfaces. The framework supports Leaky Integrate-and-fire, Quadratic Integrate-and-fire, Exponential Integrate-and-fire, Izhikevich, Event-based Gated Recurrent Unit, and a multi-compartment neuron with a model of dendritic spikes in human layer 2/3 cortical Pyramidal neurons (König et al., 9 Mar 2026).

The stochastic extension replaces ODE dynamics by Event SDEs driven by càdlàg rough paths. A neuron state $I_i(t)$0 evolves, between spike times, according to

$I_i(t)$1

with event times

$I_i(t)$2

Under the standing assumptions $I_i(t)$3, $I_i(t)$4, $I_i(t)$5 on $I_i(t)$6, $I_i(t)$7 on $I_i(t)$8, and $I_i(t)$9 on $\mathcal F_i=\{t_i^1,t_i^2,\dots\}$0, one obtains almost sure differentiability of $\mathcal F_i=\{t_i^1,t_i^2,\dots\}$1 with respect to the initial condition. The event-time gradient satisfies

$\mathcal F_i=\{t_i^1,t_i^2,\dots\}$2

and the implementation is provided as part of the diffrax library (Holberg et al., 2024).

Together, these developments move differentiable spike-time solvers from analytically tractable LIF settings toward arbitrary neuron models defined by ODEs and stochastic spiking neural networks with noise affecting both spike timing and network dynamics.

6. Representation-based alternatives and the exact-versus-surrogate debate

Differentiation on Spike Representation (DSR) is not a spike-time solver in the narrow event-equation sense, but it is directly relevant because it replaces non-differentiable spike trains by continuous spike representations and shows that the induced layer mapping is sub-differentiable. For a spike train $\mathcal F_i=\{t_i^1,t_i^2,\dots\}$3, DSR defines a layer representation $\mathcal F_i=\{t_i^1,t_i^2,\dots\}$4. In the IF model this is a simple average rate, while in the LIF model it is an exponentially weighted rate with $\mathcal F_i=\{t_i^1,t_i^2,\dots\}$5. Under mild stationarity and bounded-membrane-voltage assumptions, as $\mathcal F_i=\{t_i^1,t_i^2,\dots\}$6,

$\mathcal F_i=\{t_i^1,t_i^2,\dots\}$7

where $\mathcal F_i=\{t_i^1,t_i^2,\dots\}$8 for IF and $\mathcal F_i=\{t_i^1,t_i^2,\dots\}$9 for LIF, and the clamp is piecewise-linear and everywhere sub-differentiable (Meng et al., 2022).

This framework explicitly contrasts with both back-propagation through time and per-step surrogate gradients. DSR trains weights and per-layer thresholds through gradients of the clamp mapping, analyzes finite-$i$0 representation error into quantization error $i$1 and deviation error $i$2, and establishes the bound

$i$3

To reduce quantization error, it trains the threshold $i$4 in each layer and introduces a firing threshold offset $i$5, replacing $i$6 with $i$7. In the IF constant-input case, $i$8 makes the firing-rate quantization error minimal (Meng et al., 2022).

The relationship between exact event-time methods and surrogate gradients is treated explicitly in the forward-propagation work. There, surrogate-gradient methods are described as replacing the hard threshold $i$9 with a smooth proxy $$y_{ji}(t)=\sum_{f\in\mathcal F_j} h_{ji}(t-f), \qquad h_{ji}(t)=\frac{e^{-\alpha_j t}-e^{-\beta_i t}}{\beta_i-\alpha_j}u(t),$$0 during backprop, yielding an approximate $$y_{ji}(t)=\sum_{f\in\mathcal F_j} h_{ji}(t-f), \qquad h_{ji}(t)=\frac{e^{-\alpha_j t}-e^{-\beta_i t}}{\beta_i-\alpha_j}u(t),$$1, whereas exact methods obtain the needed sensitivity from the actual post-synaptic and refractory kernels. The same source notes that in single-spike or low-rate regimes, many surrogate shapes coincide roughly with the actual implicit-function denominator and numerator, which helps explain why some surrogates perform well (Lee et al., 2022).

Accordingly, the controversy is not simply exactness versus usefulness. Exact event-time solvers offer analytical or solver-relative correctness, but their scope may depend on model structure, regularity assumptions, or event-detection overhead. Representation-based and surrogate approaches relax different parts of the problem, often improving applicability or latency.

7. Empirical behavior, operating regimes, and limitations

Timing-based exact methods have demonstrated competitive empirical performance. On MNIST with a $$y_{ji}(t)=\sum_{f\in\mathcal F_j} h_{ji}(t-f), \qquad h_{ji}(t)=\frac{e^{-\alpha_j t}-e^{-\beta_i t}}{\beta_i-\alpha_j}u(t),$$2 feedforward architecture, the multi-spike timing-based SNN of Yamamoto et al. reaches $$y_{ji}(t)=\sum_{f\in\mathcal F_j} h_{ji}(t-f), \qquad h_{ji}(t)=\frac{e^{-\alpha_j t}-e^{-\beta_i t}}{\beta_i-\alpha_j}u(t),$$3, outperforming its single-spike TTFS variant at approximately $$y_{ji}(t)=\sum_{f\in\mathcal F_j} h_{ji}(t-f), \qquad h_{ji}(t)=\frac{e^{-\alpha_j t}-e^{-\beta_i t}}{\beta_i-\alpha_j}u(t),$$4 and matching the best rate-coded SNNs. In practice, $$y_{ji}(t)=\sum_{f\in\mathcal F_j} h_{ji}(t-f), \qquad h_{ji}(t)=\frac{e^{-\alpha_j t}-e^{-\beta_i t}}{\beta_i-\alpha_j}u(t),$$5 yields approximately $$y_{ji}(t)=\sum_{f\in\mathcal F_j} h_{ji}(t-f), \qquad h_{ji}(t)=\frac{e^{-\alpha_j t}-e^{-\beta_i t}}{\beta_i-\alpha_j}u(t),$$6 spikes per neuron on MNIST; very small $$y_{ji}(t)=\sum_{f\in\mathcal F_j} h_{ji}(t-f), \qquad h_{ji}(t)=\frac{e^{-\alpha_j t}-e^{-\beta_i t}}{\beta_i-\alpha_j}u(t),$$7 leads to explosion of spikes, very large $$y_{ji}(t)=\sum_{f\in\mathcal F_j} h_{ji}(t-f), \qquad h_{ji}(t)=\frac{e^{-\alpha_j t}-e^{-\beta_i t}}{\beta_i-\alpha_j}u(t),$$8 yields nearly single-spike behavior, and empirically $$y_{ji}(t)=\sum_{f\in\mathcal F_j} h_{ji}(t-f), \qquad h_{ji}(t)=\frac{e^{-\alpha_j t}-e^{-\beta_i t}}{\beta_i-\alpha_j}u(t),$$9 achieves best test accuracy. The same study reports that allowing multiple spikes produces an optimal $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$00 at which performance peaks, which was not seen in conventional single-spike TTFS coding (Yamamoto et al., 2022).

Eventax reports that accuracy on Yin–Yang and MNIST matches or exceeds prior continuous-time event-based methods and outperforms discrete-time surrogate-gradient methods when temporal precision matters. LIF, QIF, EIF, and Izhikevich models all train to greater than $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$01 on Yin–Yang; QIF with integral loss reaches greater than $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$02; and latency-encoded MNIST reaches $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$03. Throughput scales linearly with batch size via jax.vmap, while smaller ODE step sizes increase accuracy but cost more CPU/GPU time. Using Diffrax adaptive solvers and recursive-checkpoint adjoint lowers memory overhead compared to naïve backprop-through-time, at the cost of solver-call overhead (König et al., 9 Mar 2026).

DSR emphasizes the latency–accuracy regime. Because the representation error shrinks as $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$04, DSR-trained SNNs can reach ANN-comparable accuracy, for example within $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$05, with only $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$06 time steps on CIFAR-10, CIFAR-100, and ImageNet; typical measured firing sparsity is at most $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$07 across layers; and empirically DSR models lose less than $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$08 accuracy when going from $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$09 down to $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$10. Its analysis also states the principal limitations: the infinite-time mapping is exact only as $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$11, quantization error may impair learning for very small $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$12 such as $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$13, and one must choose $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$14 to justify the LIF averaging (Meng et al., 2022).

The earlier spike-train-to-spike-train learning results are more diagnostic than benchmark-oriented. In a single-layer setting with one neuron and 10 synapses, homogeneous 10 Hz Poisson inputs led to approximately $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$15 of random witness-learning pairs improving MAPE after 50 k updates, while inhomogeneous $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$16 Hz inputs yielded approximately $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$17 improvement after 10 k updates and complete convergence on 50 test neurons. In a two-layer $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$18 network with 30 synapses, 32 out of 50 witness-learning pairs fully converged and 18 diverged; the reported causes were non-informative spike-rate extremes, with too high a rate leading to non-uniqueness and too low a rate producing vanishing gradients. The same experiments identified an excitatory–inhibitory asymmetry: when some hidden neurons were inhibitory, their synapses converged at least $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$19 slower because inhibitory PSPs rise slowly at the threshold crossing, making $\tau_m \frac{dV_i(t)}{dt} = -V_i(t) + I_i(t),$20 smaller (Banerjee, 2014).

For stochastic spiking systems, the rough-path framework states a different class of limitations. Exact pathwise gradients are guaranteed only under the stated regularity, transversality, threshold-orthogonality, and reset-diffusion commutation assumptions. Within those conditions, the framework provides exact pathwise gradients of solution trajectories and event times, a general-purpose loss function based on signature kernels indexed on càdlàg rough paths, and an end-to-end autodifferentiable solver for Event SDEs. It is described as, to the authors’ knowledge, the first framework enabling gradient-based training of stochastic spiking neural networks with noise affecting both spike timing and network dynamics (Holberg et al., 2024).

Taken together, these results show that differentiable spike-time solvers now span exact analytical timing gradients, differentiable event-aware ODE solvers, and Event SDE formulations, while adjacent spike-representation methods provide a sub-differentiable alternative for low-latency regimes. A plausible implication is that future work will continue to negotiate the same three-way trade-off visible across the literature: event-level exactness, neuron-model generality, and computational efficiency.